本 notebook 提供了创建分层预测管道的分步指南。

在此管道中，我们将使用 NeuralForecast 和 HINT 类来创建拟合、预测和协调预测。

我们将使用 TourismL 数据集，该数据集汇总了澳大利亚全国游客大型调查数据。

目录
1. 安装软件包
2. 加载分层数据集
3. 拟合和预测 HINT
4. 基准方法
5. 预测评估

您可以使用 Google Colab 在 GPU 上运行这些实验。

​
1. 安装软件包

!pip install datasetsforecast hierarchicalforecast neuralforecast statsforecast

​
2. 加载分层数据集

这份详细的澳大利亚旅游数据集来自澳大利亚旅游研究机构管理的全国游客调查，包含 1998 年至 2016 年的 555 个月度序列，按地理和旅行目的组织。自然地理层次结构包括七个州，进一步细分为 27 个区域和 76 个地区。旅行目的类别包括度假、探亲访友 (VFR)、商务及其他。MinT (Wickramasuriya et al., 2019) 等分层预测研究过去曾使用过该数据集。该数据集可在 MinT 协调网页上访问，尽管也有其他来源可用。

地理划分每个划分的序列数每个目的的序列数总计
澳大利亚145
州72835
区域27108135
地区76304380
总计111444555
import pandas as pd

from datasetsforecast.hierarchical import HierarchicalData
from hierarchicalforecast.utils import aggregate, HierarchicalPlot
from neuralforecast.utils import augment_calendar_df
from utilsforecast.plotting import plot_series

# Load hierarchical dataset
Y_df, S_df, tags = HierarchicalData.load('./data', 'TourismLarge')
Y_df['ds'] = pd.to_datetime(Y_df['ds'])
Y_df, _ = augment_calendar_df(df=Y_df, freq='M')
S_df = S_df.reset_index(names="unique_id")

从数学上讲，分层多元时间序列可以用向量 y[a,b],t\mathbf{y}_{[a,b],t} 表示，该向量由以下聚合约束定义

y[a,b],t=S[a,b][b]y[b],t⇔[y[a],ty[b],t]=[A[a][b]I[b][b]]y[b],t \mathbf{y}_{[a,b],t} = \mathbf{S}_{[a,b][b]} \mathbf{y}_{[b],t} \quad \Leftrightarrow \quad \begin{bmatrix}\mathbf{y}_{[a],t} \\ %\hline \mathbf{y}_{[b],t}\end{bmatrix} = \begin{bmatrix} \mathbf{A}_{[a][b]}\\ %\hline \mathbf{I}_{[b][b]} \end{bmatrix} \mathbf{y}_{[b],t}

其中 y[a],t\mathbf{y}_{[a],t} 是聚合序列，y[b],t\mathbf{y}_{[b],t} 是底层序列，S[a,b][b]\mathbf{S}_{[a,b][b]} 是分层聚合约束。

# Here we plot the hierarchical constraints matrix
hplot = HierarchicalPlot(S=S_df, tags=tags)
hplot.plot_summing_matrix()


plot_series(forecasts_df=Y_df[["unique_id", "ds", "y"]], ids=['TotalAll'])

​
3. 拟合和预测 HINT

分层预测网络 (HINT) 将三个组件组合成一个易于使用的模型
1. 最先进的神经预测模型。
2. 一种高效灵活的多元概率分布。
3. 内置协调功能。

import logging

import numpy as np

from neuralforecast import NeuralForecast
from neuralforecast.models import NHITS, HINT
from neuralforecast.losses.pytorch import GMM, sCRPS

# Train test splits
horizon = 12
Y_test_df  = Y_df.groupby('unique_id', observed=True).tail(horizon)
Y_train_df = Y_df.drop(Y_test_df.index)

# Horizon and quantiles
level = np.arange(0, 100, 2)
qs = [[50-lv/2, 50+lv/2] if lv!=0 else [50] for lv in level]
quantiles = np.sort(np.concatenate(qs)/100)

# HINT := BaseNetwork + Distribution + Reconciliation
nhits = NHITS(h=horizon,
              input_size=24,
              loss=GMM(n_components=10, quantiles=quantiles),
              hist_exog_list=['month'],
              max_steps=2000,
              early_stop_patience_steps=10,
              val_check_steps=50,
              scaler_type='robust',
              learning_rate=1e-3,
              valid_loss=sCRPS(quantiles=quantiles))

model = HINT(h=horizon, S=S_df.drop(columns='unique_id').values,
             model=nhits,  reconciliation='BottomUp')

INFO:lightning_fabric.utilities.seed:Seed set to 1

logging.getLogger('pytorch_lightning').setLevel(logging.ERROR)

Y_df['y'] = Y_df['y'] * (Y_df['y'] > 0)
nf = NeuralForecast(models=[model], freq='MS')
nf.fit(df=Y_train_df, val_size=12)
Y_hat_df = nf.predict()

Y_hat_df = Y_hat_df.rename(columns=lambda x: x.replace('.0', ''))

plot_series(
    Y_df,
    Y_hat_df.drop(columns='NHITS-median'),
    ids=['TotalAll'],
    level=[90],
    max_insample_length=12*5,
)

​
4. 基准方法

我们将其与 AutoARIMA 进行比较，AutoARIMA 是 StatsForecast 软件包中一种成熟的传统预测方法，我们使用 HierarchicalForecast 对其预测结果进行协调。

from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA
from hierarchicalforecast.methods import BottomUp, MinTrace
from hierarchicalforecast.core import HierarchicalReconciliation

我们定义模型，并创建预测。

sf = StatsForecast(models=[AutoARIMA()], 
                     freq='MS', n_jobs=-1)
Y_hat_df_arima = sf.forecast(df=Y_train_df, 
                             h=12, 
                             fitted=True, 
                             X_df=Y_test_df.drop(columns="y"), 
                             level = np.arange(2, 100, 2))
Y_fitted_df_arima = sf.forecast_fitted_values()

接下来，我们使用 BottomUp 和 MinTrace 协调技术对预测进行协调

reconcilers = [
    BottomUp(),
    MinTrace(method='mint_shrink'),
]
hrec = HierarchicalReconciliation(reconcilers=reconcilers)
Y_rec_df = hrec.reconcile(Y_hat_df=Y_hat_df_arima, 
                          Y_df=Y_fitted_df_arima, 
                          S=S_df, 
                          tags=tags, 
                          level = np.arange(2, 100, 2), 
                          intervals_method="bootstrap")

​
5. 预测评估

为了评估一致的概率预测，我们使用缩放的连续秩概率分数 (sCRPS)，定义如下

CRPS(F^[a,b],τ,y[a,b],τ)=2Na+Nb∑i∫01QL(F^i,τ,yi,τ)qdq \mathrm{CRPS}(\hat{F}_{[a,b],\tau},\mathbf{y}_{[a,b],\tau}) = \frac{2}{N_{a}+N_{b}} \sum_{i} \int^{1}_{0} \mathrm{QL}(\hat{F}_{i,\tau}, y_{i,\tau})_{q} dq sCRPS(F^[a,b ],τ,y[a,b ],τ)=CRPS(F^[a,b ],τ,y[a,b ],τ)∑i∣yi,τ∣ \mathrm{sCRPS}(\hat{F}_{[a,b\,],\tau},\mathbf{y}_{[a,b\,],\tau}) = \frac{\mathrm{CRPS}(\hat{F}_{[a,b\,],\tau},\mathbf{y}_{[a,b\,],\tau})}{\sum_{i} | y_{i,\tau} |}

正如您所见，HINT 模型（使用 NHITS 作为基础模型）在最少调优的情况下也能有效地达到最先进的精度。

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

df_metrics = Y_hat_df.merge(Y_test_df.drop(columns="month"), on=['unique_id', 'ds'])
df_metrics = df_metrics.merge(Y_rec_df, on=['unique_id', 'ds'])

metrics = evaluate(df = df_metrics,
                    tags = tags,
                    metrics = [scaled_crps],
                    models= ["NHITS", "AutoARIMA"],
                    level = np.arange(2, 100, 2),
                    train_df = Y_train_df.drop(columns="month"),
                    )

metrics
层级指标NHITSAutoARIMA
0国家scaled_crps0.0444310.131136
1国家/州scaled_crps0.0634110.147516
2国家/州/区域scaled_crps0.1060600.174071
3国家/州/区域/地区scaled_crps0.1519880.205654
4国家/目的scaled_crps0.0758210.133664
5国家/州/目的scaled_crps0.1146740.181850
6国家/州/区域/目的scaled_crps0.1804910.244324
7国家/州/区域/地区/目的scaled_crps0.2454660.310656
8总计scaled_crps0.1227930.191109

​
参考文献