在许多应用中，时间序列集合是按层级组织的。例如，存在地理层级、产品或类别，它们定义了不同类型的聚合。在这种情况下，预测者通常需要为所有分解和聚合序列提供预测。一个自然的需求是这些预测应是“一致的”，也就是说，底层序列的预测值精确地加总等于聚合序列的预测值。

在本 notebook 中，我们将展示如何使用 HierarchicalForecast 在时间层级之间生成一致的预测。我们将使用经典的澳大利亚国内旅游 (Tourism) 数据集，该数据集包含澳大利亚各州游客数量的月度时间序列。

我们将首先加载 Tourism 数据，并使用 StatsForecast 中的 AutoETS 模型生成基础预测。然后，我们根据时间层级，使用 HierarchicalForecast 中的几种协调算法协调预测结果。

您可以使用 Google Colab，通过 CPU 或 GPU 运行这些实验。

!pip install hierarchicalforecast statsforecast

​
1. 加载和处理数据

在本例中，我们将使用来自 《预测：原理与实践》 一书的 Tourism 数据集。

该数据集只包含最低层级的时间序列，因此我们需要为所有层级创建时间序列。

import numpy as np
import pandas as pd

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

​
2. 时间协调

首先，我们向数据添加一个 unique_id。

Y_df["unique_id"] = Y_df["Country"] + "/" + Y_df["State"] + "/" + Y_df["Region"] + "/" + Y_df["Purpose"]

​
2a. 划分训练/测试集

我们使用最后两年（8 个季度）作为测试集。因此，我们的预测范围 horizon=8。

horizon = 8

Y_test_df = Y_df.groupby("unique_id", as_index=False).tail(horizon)
Y_train_df = Y_df.drop(Y_test_df.index)

​
2a. 根据时间层级聚合数据集

我们首先定义时间聚合规范。该规范是一个字典，其中键是聚合的名称，值是在该聚合中应聚合的底层时间步长的数量。例如，year 包含 12 个月，因此我们定义一个键值对 "yearly":12。我们可以对我们感兴趣的其他聚合做类似的事情。

在本例中，我们选择 year、semiannual 和 quarter 的时间聚合。底层时间步长是季度频率。

spec_temporal = {"year": 4, "semiannual": 2, "quarter": 1}

接下来，我们使用 aggregate_temporal 函数计算时间聚合的训练集和测试集。请注意，由于测试集包含训练集中不包含的时间层级，因此训练集和测试集的聚合矩阵 S 不同。

from hierarchicalforecast.utils import aggregate_temporal

Y_train_df, S_train_df, tags_train = aggregate_temporal(df=Y_train_df, spec=spec_temporal)
Y_test_df, S_test_df, tags_test = aggregate_temporal(df=Y_test_df,  spec=spec_temporal)

tags_train

{'year': array(['year-1', 'year-2', 'year-3', 'year-4', 'year-5', 'year-6',
        'year-7', 'year-8', 'year-9', 'year-10', 'year-11', 'year-12',
        'year-13', 'year-14', 'year-15', 'year-16', 'year-17', 'year-18'],
       dtype=object),
 'semiannual': array(['semiannual-1', 'semiannual-2', 'semiannual-3', 'semiannual-4',
        'semiannual-5', 'semiannual-6', 'semiannual-7', 'semiannual-8',
        'semiannual-9', 'semiannual-10', 'semiannual-11', 'semiannual-12',
        'semiannual-13', 'semiannual-14', 'semiannual-15', 'semiannual-16',
        'semiannual-17', 'semiannual-18', 'semiannual-19', 'semiannual-20',
        'semiannual-21', 'semiannual-22', 'semiannual-23', 'semiannual-24',
        'semiannual-25', 'semiannual-26', 'semiannual-27', 'semiannual-28',
        'semiannual-29', 'semiannual-30', 'semiannual-31', 'semiannual-32',
        'semiannual-33', 'semiannual-34', 'semiannual-35', 'semiannual-36'],
       dtype=object),
 'quarter': array(['quarter-1', 'quarter-2', 'quarter-3', 'quarter-4', 'quarter-5',
        'quarter-6', 'quarter-7', 'quarter-8', 'quarter-9', 'quarter-10',
        'quarter-11', 'quarter-12', 'quarter-13', 'quarter-14',
        'quarter-15', 'quarter-16', 'quarter-17', 'quarter-18',
        'quarter-19', 'quarter-20', 'quarter-21', 'quarter-22',
        'quarter-23', 'quarter-24', 'quarter-25', 'quarter-26',
        'quarter-27', 'quarter-28', 'quarter-29', 'quarter-30',
        'quarter-31', 'quarter-32', 'quarter-33', 'quarter-34',
        'quarter-35', 'quarter-36', 'quarter-37', 'quarter-38',
        'quarter-39', 'quarter-40', 'quarter-41', 'quarter-42',
        'quarter-43', 'quarter-44', 'quarter-45', 'quarter-46',
        'quarter-47', 'quarter-48', 'quarter-49', 'quarter-50',
        'quarter-51', 'quarter-52', 'quarter-53', 'quarter-54',
        'quarter-55', 'quarter-56', 'quarter-57', 'quarter-58',
        'quarter-59', 'quarter-60', 'quarter-61', 'quarter-62',
        'quarter-63', 'quarter-64', 'quarter-65', 'quarter-66',
        'quarter-67', 'quarter-68', 'quarter-69', 'quarter-70',
        'quarter-71', 'quarter-72'], dtype=object)}

我们的聚合矩阵将最低的时间粒度（季度）聚合到年份。

S_train_df.iloc[:5, :5]
temporal_idquarter-1quarter-2quarter-3quarter-4
0year-11.01.01.01.0
1year-20.00.00.00.0
2year-30.00.00.00.0
3year-40.00.00.00.0
4year-50.00.00.00.0
S_test_df.iloc[:5, :5]
temporal_idquarter-1quarter-2quarter-3quarter-4
0year-11.01.01.01.0
1year-20.00.00.00.0
2semiannual-11.01.00.00.0
3semiannual-20.00.01.01.0
4semiannual-30.00.00.00.0

如果您的测试集不可用（这在进行预测时通常是这种情况），则需要创建一个 future dataframe，其中包含正确的底层 unique_ids 和时间戳，以便它们可以进行时间聚合。我们可以使用 make_future_dataframe 辅助函数来完成此操作。

from hierarchicalforecast.utils import make_future_dataframe

Y_test_df_new = make_future_dataframe(Y_train_df, freq="QS", h=horizon)

然后可以在 aggregate_temporal 中使用 Y_test_df_new 来构建时间聚合的结构。

Y_test_df_new, S_test_df_new, tags_test_new = aggregate_temporal(df=Y_test_df_new,  spec=spec_temporal)

并且我们可以验证我们拥有相同的时间聚合测试集，只是 Y_test_df_new 不包含地面真值 y。

S_test_df_new
temporal_idquarter-1quarter-2quarter-3quarter-4quarter-5quarter-6quarter-7quarter-8
0year-11.01.01.01.00.00.00.00.0
1year-20.00.00.00.01.01.01.01.0
2semiannual-11.01.00.00.00.00.00.00.0
3semiannual-20.00.01.01.00.00.00.00.0
4semiannual-30.00.00.00.01.01.00.00.0
5semiannual-40.00.00.00.00.00.01.01.0
6quarter-11.00.00.00.00.00.00.00.0
7quarter-20.01.00.00.00.00.00.00.0
8quarter-30.00.01.00.00.00.00.00.0
9quarter-40.00.00.01.00.00.00.00.0
10quarter-50.00.00.00.01.00.00.00.0
11quarter-60.00.00.00.00.01.00.00.0
12quarter-70.00.00.00.00.00.01.00.0
13quarter-80.00.00.00.00.00.00.01.0
Y_test_df
temporal_idunique_iddsy
0year-1Australia/ACT/Canberra/Business2016-10-01754.139245
1year-2Australia/ACT/Canberra/Business2017-10-01809.950839
2year-1Australia/ACT/Canberra/Holiday2016-10-01735.365896
3year-2Australia/ACT/Canberra/Holiday2017-10-01834.717900
4year-1Australia/ACT/Canberra/Other2016-10-01175.239916
……………
4251quarter-4Australia/Western Australia/Experience Perth/V…2016-10-01439.699451
4252quarter-5Australia/Western Australia/Experience Perth/V…2017-01-01356.867038
4253quarter-6Australia/Western Australia/Experience Perth/V…2017-04-01302.296119
4254quarter-7Australia/Western Australia/Experience Perth/V…2017-07-01373.442070
4255quarter-8Australia/Western Australia/Experience Perth/V…2017-10-01455.316702
Y_test_df_new
temporal_idunique_idds
0year-1Australia/ACT/Canberra/Business2016-10-01
1year-2Australia/ACT/Canberra/Business2017-10-01
2year-1Australia/ACT/Canberra/Holiday2016-10-01
3year-2Australia/ACT/Canberra/Holiday2017-10-01
4year-1Australia/ACT/Canberra/Other2016-10-01
…………
4251quarter-4Australia/Western Australia/Experience Perth/V…2016-10-01
4252quarter-5Australia/Western Australia/Experience Perth/V…2017-01-01
4253quarter-6Australia/Western Australia/Experience Perth/V…2017-04-01
4254quarter-7Australia/Western Australia/Experience Perth/V…2017-07-01
4255quarter-8Australia/Western Australia/Experience Perth/V…2017-10-01

​
3b. 计算基础预测

现在，我们需要计算每个时间聚合的基础预测。以下单元格使用 AutoETS 模型计算 Y_train_df 中每个时间聚合的 基础预测。请注意，Y_hat_df 包含预测结果，但它们不一致。

另请注意，每个时间聚合的频率和范围都不同。在本例中，最低层级的频率是季度，范围为 8（构成 2 年）。因此，year 聚合的频率是年度，范围为 2。

当然，您也可以为时间聚合中的每个层级选择不同的模型 - 您可以尽情发挥创意！

from statsforecast.models import AutoETS
from statsforecast.core import StatsForecast

Y_hat_dfs = []
id_cols = ["unique_id", "temporal_id", "ds", "y"]
# We will train a model for each temporal level
for level, temporal_ids_train in tags_train.items():
    # Filter the data for the level
    Y_level_train = Y_train_df.query("temporal_id in @temporal_ids_train")
    temporal_ids_test = tags_test[level]
    Y_level_test = Y_test_df.query("temporal_id in @temporal_ids_test")
    # For each temporal level we have a different frequency and forecast horizon
    freq_level = pd.infer_freq(Y_level_train["ds"].unique())
    horizon_level = Y_level_test["ds"].nunique()
    # Train a model and create forecasts
    fcst = StatsForecast(models=[AutoETS(model='ZZZ')], freq=freq_level, n_jobs=-1)
    Y_hat_df_level = fcst.forecast(df=Y_level_train[["ds", "unique_id", "y"]], h=horizon_level, level=[80, 90])
    # Add the test set to the forecast
    Y_hat_df_level = Y_hat_df_level.merge(Y_level_test, on=["ds", "unique_id"], how="left")
    # Put cols in the right order (for readability)
    Y_hat_cols = id_cols + [col for col in Y_hat_df_level.columns if col not in id_cols]
    Y_hat_df_level = Y_hat_df_level[Y_hat_cols]
    # Append the forecast to the list
    Y_hat_dfs.append(Y_hat_df_level)

Y_hat_df = pd.concat(Y_hat_dfs, ignore_index=True)

​
3c. 协调预测结果

我们可以使用 HierarchicalReconciliation 类来协调预测结果。在本例中，我们使用 BottomUp 和 MinTrace。请注意，我们必须在 reconcile 函数中设置 temporal=True。

请注意，时间协调目前不支持样本内协调方法，例如 MinTrace(method='mint_shrink')。

from hierarchicalforecast.methods import BottomUp, MinTrace
from hierarchicalforecast.core import HierarchicalReconciliation

reconcilers = [
    BottomUp(),
    MinTrace(method="ols"),
]
hrec = HierarchicalReconciliation(reconcilers=reconcilers)
Y_rec_df = hrec.reconcile(Y_hat_df=Y_hat_df, 
                          S=S_test_df, 
                          tags=tags_test, 
                          temporal=True, 
                          level=[80, 90])

​
4. 评估

HierarchicalForecast 包包含 evaluate 函数，用于评估不同的层级。

我们评估 所有时间聚合 的时间聚合预测结果。

from hierarchicalforecast.evaluation import evaluate
from utilsforecast.losses import mae, scaled_crps

evaluation = evaluate(df = Y_rec_df.drop(columns = 'unique_id'),
                      tags = tags_test,
                      metrics = [mae, scaled_crps],
                      level = [80, 90],
                      id_col='temporal_id')

evaluation.columns = ['level', 'metric', 'Base', 'BottomUp', 'MinTrace(ols)']
numeric_cols = evaluation.select_dtypes(include="number").columns
evaluation[numeric_cols] = evaluation[numeric_cols].map('{:.3}'.format).astype(np.float64)

evaluation
层级指标基础BottomUpMinTrace(ols)
0年mae47.000050.800046.7000
1年scaled_crps0.05620.06200.0666
2半年mae29.500030.500029.1000
3半年scaled_crps0.06430.06810.0727
4季度mae19.400019.400018.7000
5季度scaled_crps0.08760.08760.0864
6总计mae26.200027.100025.7000
7总计scaled_crps0.07650.07840.0797

MinTrace(ols) 是整体最佳的点方法，在 year 和 semiannual 聚合预测以及 quarter 底层聚合预测上得分最低的 mae。然而，Base 方法在概率指标 crps 方面整体更好，其得分最低，这表明在此示例中，使用 Base 方法预测的不确定性水平更好。

​
附录：绘制 S 矩阵

from hierarchicalforecast.utils import HierarchicalPlot

我们绘制了测试集的求和矩阵。它相当直观：测试集中有两年，每年包含 4 个季度。* S 矩阵的第一行显示聚合 2016 如何通过对 2016 年的 4 个季度求和得到。* S 矩阵的第二行显示聚合 2017 如何通过对 2017 年的 4 个季度求和得到。* 接下来的 4 行显示如何获得半年聚合。* 最后几行是每个季度的单位矩阵，表示底部时间层级（每个季度）。

hplot = HierarchicalPlot(S=S_test_df, tags=tags_test, S_id_col="temporal_id")
hplot.plot_summing_matrix()