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

​
1. 层级时间序列

在许多应用中，一组时间序列是按层级组织的。例如，地理层级、产品或类别可以定义不同类型的聚合。

在这种情况下，通常要求预测者为所有分解和聚合系列提供预测。一个很自然的愿望是这些预测应该是“一致的”，即底层序列的预测值精确相加等于聚合序列的预测值。

上图显示了一个简单的层级结构，其中我们有四个底层序列、两个中间层序列和代表总聚合的顶层。其层级聚合或一致性约束如下：

yTotal,τ=yβ1,τ+yβ2,τ+yβ3,τ+yβ4,τy[a],τ=[yTotal,τ,  yβ1,τ+yβ2,τ,  yβ3,τ+yβ4,τ]⊺y[b],τ=[yβ1,τ,  yβ2,τ,  yβ3,τ,  yβ4,τ]⊺y_{\mathrm{Total},\tau} = y_{\beta_{1},\tau}+y_{\beta_{2},\tau}+y_{\beta_{3},\tau}+y_{\beta_{4},\tau} \qquad \qquad \qquad \qquad \qquad \\ \mathbf{y}_{[a],\tau}=\left[y_{\mathrm{Total},\tau},\; y_{\beta_{1},\tau}+y_{\beta_{2},\tau},\;y_{\beta_{3},\tau}+y_{\beta_{4},\tau}\right]^{\intercal} \qquad \mathbf{y}_{[b],\tau}=\left[ y_{\beta_{1},\tau},\; y_{\beta_{2},\tau},\; y_{\beta_{3},\tau},\; y_{\beta_{4},\tau} \right]^{\intercal}]⊺

幸运的是，这些约束可以用以下矩阵紧凑地表示

S[a,b][b]=[A[a][b]I[b][b]]=[1111110000111000010000100001] \mathbf{S}_{[a,b][b]} = \begin{bmatrix} \mathbf{A}_{\mathrm{[a][b]}} \\ \\ \\ \mathbf{I}_{\mathrm{[b][b]}} \\ \\ \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}

其中 A[a,b][b]\mathbf{A}_{[a,b][b]} 将底层序列聚合到上层，而 I[b][b]\mathbf{I}_{\mathrm{[b][b]}} 是单位矩阵。层级序列的表示形式如下：

y[a,b],τ=S[a,b][b]y[b],τ \mathbf{y}_{[a,b],\tau} = \mathbf{S}_{[a,b][b]} \mathbf{y}_{[b],\tau}

为了可视化一个示例，在图 2 中，可以将层级时间序列结构中的层级视为代表不同的地理聚合。例如，在图 2 中，顶层是一个国家内序列的总聚合，中间层是其州，底层是其区域。

​
2. 层级预测

为了实现“一致性”，大多数针对层级预测挑战的统计解决方案都采用了两阶段的预测修正过程。

  1. 首先，我们获得一组基础预测 y^[a,b],τ\mathbf{\hat{y}}_{[a,b],\tau}

  2. 然后，我们将它们修正为一致的预测 y~[a,b],τ\mathbf{\tilde{y}}_{[a,b],\tau}。

大多数层级预测修正方法可以通过以下变换来表示

y~[a,b],τ=S[a,b][b]P[b][a,b]y^[a,b],τ\tilde{\mathbf{y}}_{[a,b],\tau} = \mathbf{S}_{[a,b][b]} \mathbf{P}_{[b][a,b]} \hat{\mathbf{y}}_{[a,b],\tau}

HierarchicalForecast 库提供了一系列 Python 预测修正方法、数据集、评估和可视化工具。其可用的预测修正方法包括 BottomUp、TopDown、MiddleOut、MinTrace、ERM。其概率一致性方法包括 Normality、Bootstrap、PERMBU。

​
3. 最小示例

!pip install hierarchicalforecast statsforecast datasetsforecast

​
数据处理

import numpy as np
import pandas as pd

我们将创建一个合成数据集，以说明图 1 所示的层级时间序列结构。

我们将创建一个包含四个底层序列的两层结构，其中序列的聚合关系显而易见。

# Create Figure 1. synthetic bottom data
ds = pd.date_range(start='2000-01-01', end='2000-08-01', freq='MS')
y_base = np.arange(1,9)
r1 = y_base * (10**1)
r2 = y_base * (10**1)
r3 = y_base * (10**2)
r4 = y_base * (10**2)

ys = np.concatenate([r1, r2, r3, r4])
ds = np.tile(ds, 4)
unique_ids = ['r1'] * 8 + ['r2'] * 8 + ['r3'] * 8 + ['r4'] * 8
top_level = 'Australia'
middle_level = ['State1'] * 16 + ['State2'] * 16
bottom_level = unique_ids

bottom_df = dict(ds=ds,
                 top_level=top_level, 
                 middle_level=middle_level, 
                 bottom_level=bottom_level,
                 y=ys)
bottom_df = pd.DataFrame(bottom_df)
bottom_df.groupby('bottom_level').head(2)
ds顶层中间层底层y
02000-01-01澳大利亚州 1r110
12000-02-01澳大利亚州 1r120
82000-01-01澳大利亚州 1r210
92000-02-01澳大利亚州 1r220
162000-01-01澳大利亚州 2r3100
172000-02-01澳大利亚州 2r3200
242000-01-01澳大利亚州 2r4100
252000-02-01澳大利亚州 2r4200

之前介绍的层级序列 y[a,b]τ\mathbf{y}_{[a,b]\tau} 被捕获在 Y_hier_df 数据框中。

聚合约束矩阵 S[a][b]\mathbf{S}_{[a][b]} 被捕获在 S_df 数据框中。

最后，tags 包含 Y_hier_df 中的一个列表，用于构成每个层级，例如 tags['top_level'] 包含 Australia 的聚合序列索引。

from hierarchicalforecast.utils import aggregate

# Create hierarchical structure and constraints
hierarchy_levels = [['top_level'],
                    ['top_level', 'middle_level'],
                    ['top_level', 'middle_level', 'bottom_level']]
Y_hier_df, S_df, tags = aggregate(df=bottom_df, spec=hierarchy_levels)
print('S_df.shape', S_df.shape)
print('Y_hier_df.shape', Y_hier_df.shape)
print("tags['top_level']", tags['top_level'])

S_df.shape (7, 5)
Y_hier_df.shape (56, 3)
tags['top_level'] ['Australia']

Y_hier_df.groupby('unique_id').head(2)
唯一IDdsy
0澳大利亚2000-01-01220
1澳大利亚2000-02-01440
8澳大利亚/州12000-01-0120
9澳大利亚/州12000-02-0140
16澳大利亚/州22000-01-01200
17澳大利亚/州22000-02-01400
24澳大利亚/州1/r12000-01-0110
25澳大利亚/州1/r12000-02-0120
32澳大利亚/州1/r22000-01-0110
33澳大利亚/州1/r22000-02-0120
40澳大利亚/州2/r32000-01-01100
41澳大利亚/州2/r32000-02-01200
48澳大利亚/州2/r42000-01-01100
49澳大利亚/州2/r42000-02-01200
S_df
唯一ID澳大利亚/州1/r1澳大利亚/州1/r2澳大利亚/州2/r3澳大利亚/州2/r4
0澳大利亚1.01.01.01.0
1澳大利亚/州11.01.00.00.0
2澳大利亚/州20.00.01.01.0
3澳大利亚/州1/r11.00.00.00.0
4澳大利亚/州1/r20.01.00.00.0
5澳大利亚/州2/r30.00.01.00.0
6澳大利亚/州2/r40.00.00.01.0

​
基础预测

接下来，我们使用 naive 模型计算每个时间序列的基础预测。请注意，Y_hat_df 包含预测结果，但它们并不一致。

from statsforecast.models import Naive
from statsforecast.core import StatsForecast

# Split train/test sets
Y_test_df  = Y_hier_df.groupby('unique_id', as_index=False).tail(4)
Y_train_df = Y_hier_df.drop(Y_test_df.index)

# Compute base Naive predictions
# Careful identifying correct data freq, this data monthly 'M'
fcst = StatsForecast(models=[Naive()],
                     freq='MS', n_jobs=-1)
Y_hat_df = fcst.forecast(df=Y_train_df, h=4, fitted=True)
Y_fitted_df = fcst.forecast_fitted_values()

​
预测修正

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

# You can select a reconciler from our collection
reconcilers = [BottomUp()] # MinTrace(method='mint_shrink')
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)
Y_rec_df.groupby('unique_id').head(2)
唯一IDds朴素法朴素法/自下而上
0澳大利亚2000-05-01880.0880.0
1澳大利亚2000-06-01880.0880.0
4澳大利亚/州12000-05-0180.080.0
5澳大利亚/州12000-06-0180.080.0
8澳大利亚/州22000-05-01800.0800.0
9澳大利亚/州22000-06-01800.0800.0
12澳大利亚/州1/r12000-05-0140.040.0
13澳大利亚/州1/r12000-06-0140.040.0
16澳大利亚/州1/r22000-05-0140.040.0
17澳大利亚/州1/r22000-06-0140.040.0
20澳大利亚/州2/r32000-05-01400.0400.0
21澳大利亚/州2/r32000-06-01400.0400.0
24澳大利亚/州2/r42000-05-01400.0400.0
25澳大利亚/州2/r42000-06-01400.0400.0

​
参考文献