地理聚合（旅游）

关于澳大利亚旅游数据的地理层级预测

在许多应用中，一组时间序列是以层次结构组织的。示例包括定义不同类型聚合的地理级别、产品或类别。在这种情况下，预测人员通常需要为所有分散和聚合系列提供预测。一个自然的愿望是这些预测能够 “一致”，即底层系列的总和恰好等于聚合系列的预测。

在本笔记本中，我们展示了如何使用 HierarchicalForecast 在地理层级之间生成一致的预测。我们将使用经典的澳大利亚国内旅游（Tourism）数据集，该数据集包含每个澳大利亚州每月访客数量的时间序列。

我们将首先加载 Tourism 数据，并使用来自 StatsForecast 的 ETS 模型生成基本预测，然后使用来自 HierarchicalForecast 的几个调和算法对预测进行协调。最后，我们显示其性能与 Forecasting: Principles and Practice 中使用 R 包 fable 报告的结果相当。

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

%%capture
!pip install hierarchicalforecast
!pip install -U statsforecast numba

1. 加载和处理数据

在这个例子中，我们将使用旅游数据集，该数据集来自预测：原理与实践书籍。

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

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.to_datetime(Y_df['ds'])
Y_df.head()

	Country	Region	State	Purpose	ds	y
0	Australia	Adelaide	South Australia	Business	1998-01-01	135.077690
1	Australia	Adelaide	South Australia	Business	1998-04-01	109.987316
2	Australia	Adelaide	South Australia	Business	1998-07-01	166.034687
3	Australia	Adelaide	South Australia	Business	1998-10-01	127.160464
4	Australia	Adelaide	South Australia	Business	1999-01-01	137.448533

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

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

使用HierarchicalForecast中的aggregate函数，我们可以获得完整的时间序列集。

from hierarchicalforecast.utils import aggregate

%%capture
Y_df, S_df, tags = aggregate(Y_df, spec)
Y_df = Y_df.reset_index()

Y_df.head()

	unique_id	ds	y
0	Australia	1998-01-01	23182.197269
1	Australia	1998-04-01	20323.380067
2	Australia	1998-07-01	19826.640511
3	Australia	1998-10-01	20830.129891
4	Australia	1999-01-01	22087.353380

S_df.iloc[:5, :5]

	Australia/ACT/Canberra/Business	Australia/ACT/Canberra/Holiday	Australia/ACT/Canberra/Other	Australia/ACT/Canberra/Visiting	Australia/New South Wales/Blue Mountains/Business
Australia	1.0	1.0	1.0	1.0	1.0
Australia/ACT	1.0	1.0	1.0	1.0	0.0
Australia/New South Wales	0.0	0.0	0.0	0.0	1.0
Australia/Northern Territory	0.0	0.0	0.0	0.0	0.0
Australia/Queensland	0.0	0.0	0.0	0.0	0.0

tags['Country/Purpose']

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

分割训练/测试集

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

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

Y_test_df = Y_test_df.set_index('unique_id')
Y_train_df = Y_train_df.set_index('unique_id')

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

2. 计算基础预测

以下单元计算 Y_df 中每个时间序列的 基础预测，使用的是 ETS 模型。请注意，Y_hat_df 包含预测结果，但它们并不连贯。

%%capture
from statsforecast.models import ETS
from statsforecast.core import StatsForecast

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

3. 调和预测

以下单元使用 HierarchicalReconciliation 类使先前的预测保持一致。由于层级结构不是严格的，我们无法使用 TopDown 或 MiddleOut 等方法。在这个例子中，我们使用 BottomUp 和 MinTrace。

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

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)

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

Y_rec_df.head()

	ds	ETS	ETS/BottomUp	ETS/MinTrace_method-mint_shrink	ETS/MinTrace_method-ols
unique_id
Australia	2016-01-01	25990.068359	24379.679688	25438.888351	25894.418893
Australia	2016-04-01	24458.490234	22902.664062	23925.188541	24357.230480
Australia	2016-07-01	23974.056641	22412.984375	23440.310338	23865.929521
Australia	2016-10-01	24563.455078	23127.638672	24101.001833	24470.783968
Australia	2017-01-01	25990.068359	24516.175781	25556.667616	25901.382401

4. 评估

HierarchicalForecast 包含 HierarchicalEvaluation 类，用于评估不同的层次结构，并能够计算与基准模型相比的标准化指标。

from hierarchicalforecast.evaluation import HierarchicalEvaluation

def rmse(y, y_hat):
    return np.mean(np.sqrt(np.mean((y-y_hat)**2, axis=1)))

def mase(y, y_hat, y_insample, seasonality=4):
    errors = np.mean(np.abs(y - y_hat), axis=1)
    scale = np.mean(np.abs(y_insample[:, seasonality:] - y_insample[:, :-seasonality]), axis=1)
    return np.mean(errors / scale)

eval_tags = {}
eval_tags['Total'] = tags['Country']
eval_tags['Purpose'] = tags['Country/Purpose']
eval_tags['State'] = tags['Country/State']
eval_tags['Regions'] = tags['Country/State/Region']
eval_tags['Bottom'] = tags['Country/State/Region/Purpose']
eval_tags['All'] = np.concatenate(list(tags.values()))

evaluator = HierarchicalEvaluation(evaluators=[rmse, mase])
evaluation = evaluator.evaluate(
        Y_hat_df=Y_rec_df, Y_test_df=Y_test_df,
        tags=eval_tags, Y_df=Y_train_df
)
evaluation = evaluation.drop('Overall')
evaluation.columns = ['Base', 'BottomUp', 'MinTrace(mint_shrink)', 'MinTrace(ols)']
evaluation = evaluation.applymap('{:.2f}'.format)

/var/folders/rp/97y9_3ns23v01hdn0rp9ndw40000gp/T/ipykernel_46857/2768439279.py:22: PerformanceWarning: dropping on a non-lexsorted multi-index without a level parameter may impact performance.
  evaluation = evaluation.drop('Overall')

RMSE

下表显示了对每种调和方法在各级别上使用RMSE测量的性能。

evaluation.query('metric == "rmse"')

		Base	BottomUp	MinTrace(mint_shrink)	MinTrace(ols)
level	metric
Total	rmse	1743.29	3028.93	2102.47	1818.94
Purpose	rmse	534.75	791.28	574.84	515.53
State	rmse	308.15	413.43	315.89	287.34
Regions	rmse	51.65	55.14	46.48	46.29
Bottom	rmse	19.37	19.37	17.78	18.19
All	rmse	45.19	54.95	44.59	42.71

MASE

下表显示了使用MASE测量的各个层级的每种协调方法的性能。

evaluation.query('metric == "mase"')

		Base	BottomUp	MinTrace(mint_shrink)	MinTrace(ols)
level	metric
Total	mase	1.59	3.16	2.05	1.67
Purpose	mase	1.32	2.28	1.48	1.25
State	mase	1.39	1.90	1.39	1.25
Regions	mase	1.12	1.19	1.01	0.99
Bottom	mase	0.98	0.98	0.94	1.01
All	mase	1.03	1.08	0.98	1.02

比较 fable

请注意，我们可以恢复 Forecasting: Principles and Practice 中报告的结果。原始结果是使用 R 包 fable 计算得出的。

参考文献

Hyndman, R.J., & Athanasopoulos, G. (2021). “预测：原理与实践，第三版：第11章：分层与分组系列的预测。” OTexts: 墨尔本，澳大利亚。OTexts.com/fpp3 访问于2022年7月。
Rob Hyndman, Alan Lee, Earo Wang, Shanika Wickramasuriya, 和维护者Earo Wang (2021). “hts: 分层与分组时间序列”. URL https://CRAN.R-project.org/package=hts. R包版本0.3.1.
Mitchell O’Hara-Wild, Rob Hyndman, Earo Wang, Gabriel Caceres, Tim-Gunnar Hensel, 和 Timothy Hyndman (2021). “fable: 用于整洁时间序列的预测模型”. URL https://CRAN.R-project.org/package=fable. R包版本6.0.2.

If you find the code useful, please ⭐ us on Github