多重季节性

提示

对于此任务，StatsForecast 的 MSTL 比 Prophet 和 NeuralProphet 准确率高 68%，速度快 600%。(在此重现实验)

多重季节性数据是指具有不止一种明显季节性的时间序列。多重季节性通常出现在低频采样的数据中。例如，每小时电力数据显示出每日和每周季节性。这意味着一天中的特定小时（如下午 6:00 对比凌晨 3:00）或特定日期（如周日 对比 周五）存在明显的电力消费模式。

传统统计模型无法建模多种季节性长度。在此示例中，我们将展示如何使用 LOESS（局部加权散点图平滑）进行多重季节性趋势分解 (MSTL) 来高效地建模这两种季节性。

在此示例中，我们将使用来自宾夕法尼亚州、新泽西州和马里兰州 (PJM) 的每小时电力负荷数据。原始数据可在此处找到。（点击此处了解 PJM 信息）

首先，我们将加载数据，然后使用 StatsForecast.fit 和 StatsForecast.predict 方法预测接下来的 24 小时。然后，我们将把时间序列的不同元素分解为趋势及其多重季节性。最后，您将使用 StatsForecast.forecast 进行生产级预测。

大纲

1. 安装库
2. 加载和探索数据
3. 拟合多重季节性模型
4. 将序列分解为趋势和季节性
5. 预测接下来的 24 小时
6. 可选：在生产环境中预测

提示

您可以使用 Colab 交互式运行此 Notebook

​

安装库

我们假设您已经安装了 StatsForecast。请查看本指南了解如何安装 StatsForecast 的说明。

使用 pip install statsforecast 安装必要的包 “

​

加载数据

StatsForecast 的输入始终是长格式的数据框，包含三列：unique_id、ds 和 y

• unique_id（字符串、整数或类别）代表序列的标识符。

• ds（日期戳或整数）列应为时间索引的整数或日期戳，理想格式为日期的 YYYY-MM-DD 或时间戳的 YYYY-MM-DD HH:MM:SS。

• y（数值）代表我们希望预测的测量值。我们将重命名

您将使用 pandas 读取数据并更改必要的名称。此步骤大约需要 2 秒。

import pandas as pd

df = pd.read_csv('https://raw.githubusercontent.com/panambY/Hourly_Energy_Consumption/master/data/PJM_Load_hourly.csv')
df.columns = ['ds', 'y']
df.insert(0, 'unique_id', 'PJM_Load_hourly')
df = df.sort_values(['unique_id', 'ds']).reset_index(drop=True)
df.tail()

StatsForecast 可以处理未排序的数据，但是为了绘图方便，最好对数据框进行排序。

使用 StatsForecast 类中的 plot 方法绘制序列。此方法从数据集中打印最多 8 个随机序列，对于基本 EDA 非常有用。在本例中，由于只有一个 unique_id，它将只打印一个序列。

注意

StatsForecast.plot 方法默认使用 matplotlib 作为引擎。您可以通过设置 engine="plotly" 来切换到 plotly。

from statsforecast import StatsForecast

StatsForecast.plot(df)

时间序列表现出季节性模式。此外，该时间序列包含 32,896 个观测值，因此有必要使用计算效率非常高的方法。

​

拟合 MSTL 模型

由 Kasun Bandara、Rob J Hyndman 和 Christoph Bergmeir 最初开发的 MSTL（使用 LOESS 的多重季节趋势分解）模型使用局部多项式回归 (LOESS) 将时间序列分解为多重季节性。然后，它使用非季节模型预测趋势，并使用 SeasonalNaive 模型预测每个季节性。您可以选择用于预测 MSTL 模型趋势分量的非季节模型。在此示例中，我们将使用 AutoARIMA。

导入您需要的模型。

from statsforecast.models import MSTL, AutoARIMA

首先，我们必须定义模型参数。如前所述，电力负荷每 24 小时（每小时）和每 24 * 7 小时（每天）呈现季节性。因此，我们将使用 [24, 24 * 7] 作为季节长度。趋势分量将使用 AutoARIMA 模型进行预测。（您也可以尝试使用：AutoTheta、AutoCES 和 AutoETS）

# Create a list of models and instantiation parameters

models = [MSTL(
    season_length=[24, 24 * 7], # seasonalities of the time series 
    trend_forecaster=AutoARIMA() # model used to forecast trend
)]

我们通过实例化一个新的 StatsForecast 对象来拟合模型，该对象具有以下必需参数

• models：模型列表。从模型中选择您想要的模型并导入它们。

• freq：表示数据频率的字符串。（参见 panda 的可用频率。）

任何设置都会传递到构造函数中。然后调用其 fit 方法并传入历史数据框。

sf = StatsForecast(
    models=models, # model used to fit each time series 
    freq='h', # frequency of the data
)

提示

StatsForecast 还支持此可选参数。

• n_jobs：n_jobs：int，并行处理中使用的作业数，使用 -1 表示所有核心。（默认值：1）

• fallback_model：如果模型失败时使用的模型。（默认值：无）

使用 fit 方法将每个模型拟合到每个时间序列。在本例中，我们只将一个模型拟合到一个序列。请查看本指南了解如何将多个模型拟合到多个序列。

注意

StatsForecast 通过 Numba 的 JIT 编译实现了惊人的速度。第一次调用 statsforecast 类时，fit 方法大约需要 10 秒。第二次（一旦 Numba 编译了您的设置）应该会少于 5 秒。

sf = sf.fit(df=df)

​

分解序列

模型拟合后，可以使用 StatsForecast 的 fitted_ 属性访问分解结果。此属性存储每个时间序列拟合模型的所有相关信息。

在本例中，我们为单个时间序列拟合单个模型，因此通过访问 fitted_ 位置 [0, 0]，我们将找到模型的相关信息。MSTL 类会生成一个 model_ 属性，其中包含序列的分解方式。

sf.fitted_[0, 0].model_

我们将使用 matplotlib 可视化序列的不同组成部分。

import matplotlib.pyplot as plt

sf.fitted_[0, 0].model_.tail(24 * 28).plot(subplots=True, grid=True)
plt.tight_layout()
plt.show()

我们观察到明显的上升趋势（橙线）以及每天（24H）和每周（168H）重复的季节性。

​

预测接下来的 24 小时

带置信区间的概率预测

使用 predict 方法生成预测。

predict 方法接受两个参数：预测未来 h（预测范围）和 level。

• h (int)：表示未来 h 步的预测。在此例中，预测未来 12 个月。

• level（浮点数列表）：此可选参数用于概率预测。设置预测区间的 level（或置信度百分位数）。例如，level=[90] 表示模型预期实际值有 90% 的时间落在此区间内。

此处的预测对象是一个新的数据框，包含一列模型名称和 y hat 值，以及不确定性区间的列。

此步骤应少于 1 秒。

forecasts = sf.predict(h=24, level=[90])
forecasts.head()

您可以通过调用 StatsForecast.plot 方法并传入预测数据框来绘制预测结果。

sf.plot(df, forecasts, max_insample_length=24 * 7)

​

在生产环境中预测

如果您想在具有多个序列或模型的生产环境中提高速度，我们建议使用 StatsForecast.forecast 方法代替 .fit 和 .predict。

主要区别在于 .forecast 不存储拟合值，并且在分布式环境中具有高度可伸缩性。

forecast 方法接受两个参数：预测未来 h（预测范围）和 level。

• h (int)：表示未来 h 步的预测。在此例中，预测未来 12 个月。

• level（浮点数列表）：此可选参数用于概率预测。设置预测区间的 level（或置信度百分位数）。例如，level=[90] 表示模型预期实际值有 90% 的时间落在此区间内。

此处的预测对象是一个新的数据框，包含一列模型名称和 y hat 值，以及不确定性区间的列。根据您的计算机，此步骤大约需要 1 分钟。（如果您想加快速度到几秒钟，请移除 AutoModels，例如 ARIMA 和 Theta）

注意

StatsForecast 通过 Numba 的 JIT 编译实现了惊人的速度。第一次调用 statsforecast 类时，fit 方法大约需要 10 秒。第二次（一旦 Numba 编译了您的设置）应该会少于 5 秒。

forecasts_df = sf.forecast(df=df, h=24, level=[90])
forecasts_df.head()

​

参考资料

• Bandara, Kasun & Hyndman, Rob & Bergmeir, Christoph. (2021). “MSTL: 一种用于具有多重季节性时间序列的季节趋势分解算法”.

​

下一步

• 学习如何使用交叉验证评估模型的稳健性。