在许多情况下，只有层次结构最低级别的时间序列（底部时间序列）可用。HierarchicalForecast 提供了创建所有层次结构时间序列的工具，还允许您计算所有层次结构的预测区间。在本笔记本中，我们将了解如何实现这一点。

!pip install hierarchicalforecast statsforecast

import pandas as pd

# compute base forecast no coherent
from statsforecast.models import AutoETS
from statsforecast.core import StatsForecast

#obtain hierarchical reconciliation methods and evaluation
from hierarchicalforecast.methods import BottomUp, MinTrace
from hierarchicalforecast.utils import aggregate, HierarchicalPlot
from hierarchicalforecast.core import HierarchicalReconciliation

​
聚合底部时间序列

在本示例中，我们将使用来自《预测：原理与实践》一书的 旅游 数据集。该数据集只包含最低级别的时间序列，因此我们需要为所有层次结构创建时间序列。

Y_df = pd.read_csv('https://raw.githubusercontent.com/Nixtla/transfer-learning-time-series/main/datasets/tourism.csv')
Y_df = Y_df.rename({'Trips': 'y', 'Quarter': 'ds'}, axis=1)
Y_df.insert(0, 'Country', 'Australia')
Y_df = Y_df[['Country', 'Region', 'State', 'Purpose', 'ds', 'y']]
Y_df['ds'] = Y_df['ds'].str.replace(r'(\d+) (Q\d)', r'\1\2', regex=True)
Y_df['ds'] = pd.PeriodIndex(Y_df["ds"], freq='Q').to_timestamp()
Y_df.head()
国家区域州目的dsy
0澳大利亚阿德莱德南澳大利亚商业1998-01-01135.077690
1澳大利亚阿德莱德南澳大利亚商业1998-04-01109.987316
2澳大利亚阿德莱德南澳大利亚商业1998-07-01166.034687
3澳大利亚阿德莱德南澳大利亚商业1998-10-01127.160464
4澳大利亚阿德莱德南澳大利亚商业1999-01-01137.448533

数据集可以按以下非严格层次结构进行分组。

spec = [
    ['Country'],
    ['Country', 'State'], 
    ['Country', 'Purpose'], 
    ['Country', 'State', 'Region'], 
    ['Country', 'State', 'Purpose'], 
    ['Country', 'State', 'Region', 'Purpose']
]

使用 HierarchicalForecast 中的 aggregate 函数，我们可以生成：1. Y_df：层次结构序列 y[a,b]τ\mathbf{y}_{[a,b]\tau} 2. S_df：聚合约束数据框 S[a,b]S_{[a,b]} 3. tags：一个列表，包含构成每个聚合级别的“unique_ids”。

Y_df, S_df, tags = aggregate(df=Y_df, spec=spec)

Y_df.head()
unique_iddsy
0澳大利亚1998-01-0123182.197269
1澳大利亚1998-04-0120323.380067
2澳大利亚1998-07-0119826.640511
3澳大利亚1998-10-0120830.129891
4澳大利亚1999-01-0122087.353380
S_df.iloc[:5, :5]
unique_idAustralia/ACT/Canberra/BusinessAustralia/ACT/Canberra/HolidayAustralia/ACT/Canberra/OtherAustralia/ACT/Canberra/Visiting
0澳大利亚1.01.01.01.0
1Australia/ACT1.01.01.01.0
2Australia/New South Wales0.00.00.00.0
3Australia/Northern Territory0.00.00.00.0
4Australia/Queensland0.00.00.00.0
tags['Country/Purpose']

array(['Australia/Business', 'Australia/Holiday', 'Australia/Other',
       'Australia/Visiting'], dtype=object)

我们可以使用 HierarchicalPlot 类可视化 S_df 数据框和 Y_df，如下所示。

hplot = HierarchicalPlot(S=S_df, tags=tags)

hplot.plot_summing_matrix()


hplot.plot_hierarchically_linked_series(
    bottom_series='Australia/ACT/Canberra/Holiday',
    Y_df=Y_df
)

​
拆分训练/测试集

我们使用最后两年（8个季度）作为测试集。

Y_test_df = Y_df.groupby('unique_id', as_index=False).tail(8)
Y_train_df = Y_df.drop(Y_test_df.index)

Y_train_df.groupby('unique_id').size()

unique_id
Australia                                                72
Australia/ACT                                            72
Australia/ACT/Business                                   72
Australia/ACT/Canberra                                   72
Australia/ACT/Canberra/Business                          72
                                                         ..
Australia/Western Australia/Experience Perth/Other       72
Australia/Western Australia/Experience Perth/Visiting    72
Australia/Western Australia/Holiday                      72
Australia/Western Australia/Other                        72
Australia/Western Australia/Visiting                     72
Length: 425, dtype: int64

​
计算基础预测

以下单元计算 Y_df 中每个时间序列使用 AutoETS 模型得到的基础预测。请注意，Y_hat_df 包含预测值，但它们不一致。由于我们使用 bootstrapping 计算预测区间，因此只需要模型的拟合值。

fcst = StatsForecast(models=[AutoETS(season_length=4, model='ZAA')],
                     freq='QS', n_jobs=-1)
Y_hat_df = fcst.forecast(df=Y_train_df, h=8, fitted=True)
Y_fitted_df = fcst.forecast_fitted_values()

​
协调基础预测

以下单元使用 HierarchicalReconciliation 类使之前的预测一致。由于层次结构不是严格的，我们不能使用 TopDown 或 MiddleOut 等方法。在本示例中，我们使用 BottomUp 和 MinTrace。如果要计算预测区间，必须按如下所示使用 level 参数并设置 intervals_method='bootstrap'。

reconcilers = [
    BottomUp(),
    MinTrace(method='mint_shrink'),
    MinTrace(method='ols')
]
hrec = HierarchicalReconciliation(reconcilers=reconcilers)
Y_rec_df = hrec.reconcile(Y_hat_df=Y_hat_df, Y_df=Y_fitted_df, S=S_df, 
                          tags=tags, level=[80, 90], 
                          intervals_method='bootstrap')

数据框 Y_rec_df 包含协调后的预测。

Y_rec_df.head()
unique_iddsAutoETSAutoETS/BottomUpAutoETS/BottomUp-lo-90AutoETS/BottomUp-lo-80AutoETS/BottomUp-hi-80AutoETS/BottomUp-hi-90AutoETS/MinTrace_method-mint_shrinkAutoETS/MinTrace_method-mint_shrink-lo-90AutoETS/MinTrace_method-mint_shrink-lo-80AutoETS/MinTrace_method-mint_shrink-hi-80AutoETS/MinTrace_method-mint_shrink-hi-90AutoETS/MinTrace_method-olsAutoETS/MinTrace_method-ols-lo-90AutoETS/MinTrace_method-ols-lo-80AutoETS/MinTrace_method-ols-hi-80AutoETS/MinTrace_method-ols-hi-90
0澳大利亚2016-01-0126080.87848824487.15250323242.75731123332.59296825379.82948625424.13913725521.55170624407.44271224698.93147926357.02435426466.74068226034.13209124914.19903825100.47050227102.74606527176.467048
1澳大利亚2016-04-0124587.01211523068.31429221823.91910021910.61505723945.98294924278.68324324106.52247923185.40363423283.90225125098.33234225473.23994924567.45791323483.98381423640.62712625709.79287025809.220444
2澳大利亚2016-07-0124147.30774422686.98393321293.52944921526.52561023697.85993124150.87978923717.61050122603.50150722802.77130824802.97326025228.79562924150.11124623030.17819323154.97243625359.91799325404.792198
3澳大利亚2016-10-0124794.04077923428.03763722034.58315322273.82695724241.84044024438.91363524472.93911523361.28551223584.82587125338.71399525469.42662324831.54072123725.92746323836.40191125900.15469525977.249268
4澳大利亚2017-01-0126283.99865424939.63761623695.21755423903.39571325815.63868225973.16460726029.32272424948.33979525144.17903026900.06846127119.07316026348.22975825254.68223425487.51809827410.89415827477.330557

​
绘制预测图

然后我们可以使用以下函数绘制概率预测图。

plot_df = Y_df.merge(Y_rec_df, on=['unique_id', 'ds'], how="outer")

​
绘制单个时间序列图

hplot.plot_series(
    series='Australia',
    Y_df=plot_df, 
    models=['y', 'AutoETS', 'AutoETS/MinTrace_method-ols', 'AutoETS/MinTrace_method-mint_shrink'],
    level=[80]
)


# Since we are plotting a bottom time series
# the probabilistic and mean forecasts
# differ due to bootstrapping
hplot.plot_series(
    series='Australia/Western Australia/Experience Perth/Visiting',
    Y_df=plot_df, 
    models=['y', 'AutoETS', 'AutoETS/BottomUp'],
    level=[80]
)

​
绘制层次关联的时间序列图

hplot.plot_hierarchically_linked_series(
    bottom_series='Australia/Western Australia/Experience Perth/Visiting',
    Y_df=plot_df, 
    models=['y', 'AutoETS', 'AutoETS/MinTrace_method-ols', 'AutoETS/BottomUp'],
    level=[80]
)


# ACT only has Canberra
hplot.plot_hierarchically_linked_series(
    bottom_series='Australia/ACT/Canberra/Other',
    Y_df=plot_df, 
    models=['y', 'AutoETS/MinTrace_method-mint_shrink'],
    level=[80, 90]
)

​
参考文献