​

目录

• 引言
• ETS 模型
• ETS 估计
• 模型选择
• 加载库和数据
• 使用 plot 方法探索数据
• 将数据分为训练集和测试集
• 使用 StatsForecast 实现 AutoETS
• 交叉验证
• 模型评估
• 参考文献

​

引言

商业中经常需要对大量单变量时间序列进行自动预测。通常有超过一千条产品线需要至少每月进行预测。即使只需要少量预测，也可能没有人接受过足够的时间序列模型使用培训来生成它们。在这种情况下，自动预测算法是一个必不可少的工具。自动预测算法必须确定合适的时间序列模型，估计参数并计算预测。它们必须对异常时间序列模式具有鲁棒性，并且无需用户干预即可应用于大量序列。最流行的自动预测算法基于指数平滑或 ARIMA 模型。

​

指数平滑

虽然指数平滑方法自 20 世纪 50 年代就已出现，但直到最近才开发出包含模型选择过程的建模框架。Ord, Koehler 和 Snyder (1997)、Hyndman, Koehler, Snyder 和 Grose (2002) 以及 Hyndman, Koehler, Ord 和 Snyder (2005b) 的研究表明，所有指数平滑方法（包括非线性方法）都是创新状态空间模型的最优预测。

指数平滑方法最初由 Pegels’ (1969) 分类法进行分类。后来由 Gardner (1985) 扩展，由 Hyndman 等人 (2002) 修改，并由 Taylor (2003) 再次扩展，总共给出了下表中所示的十五种方法。

其中一些方法以其他名称更广为人知。例如，单元格 (N,N) 描述了 简单指数平滑 (或 SES) 方法，单元格 (A,N) 描述了 Holt 线性方法，单元格 (Ad,N) 描述了 阻尼趋势方法。加法 Holt-Winters 方法由单元格 (A,A) 给出，乘法 Holt-Winters 方法由单元格 (A,M) 给出。其他单元格对应于不太常用但类似的方法。

​

所有方法的点预测

我们用 $y_1,y_2,...,y_n$ 表示观测到的时间序列。基于截止到时间 $t$ 的所有数据的 $y_{t+h}$ 预测表示为 $\hat y_{t+h|t}$。为了说明该方法，我们给出方法 (A,A)（即 Holt-Winters 加法方法）的点预测和更新方程

其中 $m$ 是季节性周期长度（例如，一年中的月数或季度数），$\ell_{t}$ 表示序列的水平，$b_t$ 表示增长，$s_t$ 是季节性分量，$\hat y_{t+h|t}$ 是超前 $h$ 期的预测，并且 $h_{m}^{+} = [(h − 1) mod \ m] + 1$。要使用方法 (1)，我们需要初始状态 $\ell_{0}$、$b_0$ 和 $s_{1−m}, . . . , s_0$ 的值，以及平滑参数 $\alpha, \beta^{*}$ 和 $\gamma$ 的值。所有这些都将从观测数据中估计得出。

方程 (1c) 与 Makridakis 等人 (1998) 或 Bowerman, O’Connell 和 Koehler (2005) 中的常用 Holt-Winters 方程略有不同。这些作者将 (1c) 替换为 $s_{t} = \gamma^* (y_{t}-\ell_{t})+ (1-\gamma^*)s_{t-m}.$。

如果使用 (1a) 代入 $\ell_{t}$，我们得到

$s_{t} = \gamma^*(1-\alpha) (y_{t}-\ell_{t-1}-b_{t-1})+ [1-\gamma^*(1-\alpha)]s_{t-m},$

因此，通过将 (1c) 中的 $\gamma$ 替换为 $\gamma^{*} (1-\alpha)$，我们使用此方法获得相同的预测。Ord 等人 (1997) 提出了 (1c) 中的修改，以简化状态空间表达。它等同于 Archibald (1990) 对 Holt-Winters 方法的变体。

​

创新状态空间模型

对于 其他 ETS 模型 中的每种指数平滑方法，Hyndman 等人 (2008b) 描述了两种可能的创新状态空间模型，一种对应于加法误差模型，另一种对应于乘法误差模型。如果使用相同的参数值，这两种模型会给出相同的点预测，但预测区间不同。因此，在此分类中描述了 30 种潜在模型。

从历史上看，误差分量的性质经常被忽略，因为加法误差和乘法误差之间的区别对点预测没有影响。

我们谨慎区分指数平滑方法与底层状态空间模型。指数平滑方法是一种仅用于生成点预测的算法。底层的随机状态空间模型提供相同的点预测，但也提供了计算预测区间和其他属性的框架。

为了区分加法误差模型和乘法误差模型（以及区分模型与方法），我们在上面表格中的分类基础上添加了第三个字母。我们将每个状态空间模型标记为 ETS(⋅,.,.)，分别代表（误差 Error、趋势 Trend、季节性 Seasonal）。这个标签也可以理解为 ExponenTial Smoothing（指数平滑）。使用上面表格中相同的符号，每个分量（或状态）的可能性如下：Error ={ A,M }，趋势 ={N,A,Ad}，季节性 ={ N,A,M }。

指定模型后，我们可以研究该序列未来值的概率分布，并例如找到已知历史信息下未来观测值的条件均值。我们将其表示为 $\mu_{t+h|t} = E(y_{t+h | xt})$，其中 xt 包含未观测分量，例如 $\ell_t$、$b_t$ 和 $s_t$。对于 $h = 1$，我们使用 $\mu_t ≡ \mu_{t+1|t}$ 作为简写符号。对于许多模型，这些条件均值将与表 其他 ETS 模型 中给出的点预测值相同，因此 $\mu_{t+h|t} = \hat y_{t+h|t}$。然而，对于其他模型（具有乘法趋势或乘法季节性的模型），当 $h ≥ 2$ 时，条件均值和点预测值会略有不同。

我们将使用 Gardner 和 McKenzie (1985) 的阻尼趋势方法来阐述这些思想。

每个模型都由一个描述观测数据的观测方程和一些描述未观测分量或状态（水平项、趋势项、季节项）如何随时间变化的状态方程组成。因此，这些模型被称为状态空间模型。

对于每种方法，存在两种模型：一种带有加法误差，另一种带有乘法误差。如果使用相同的平滑参数值，这些模型产生的点预测值是相同的。但是，它们会生成不同的预测区间。

为了区分带有加法误差的模型和带有乘法误差的模型（以及区分模型与方法），我们在上表中所示的分类中添加了第三个字母。我们将每个状态空间模型标记为 ETS(⋅,.,.)，分别表示 (误差、趋势、季节)。这个标签也可以理解为指数平滑（ExponenTial Smoothing）。使用与上表相同的符号，每个分量（或状态）的可能性为：误差 = {A,M}，趋势 = {N,A,Ad}，季节 = {N,A,M}。

​

ETS(A,N,N)：具有加法误差的简单指数平滑

回顾简单指数平滑的分量形式

如果我们重新排列水平的平滑方程，我们得到“误差修正”形式，

其中 $e_{t}=y_{t}-\ell_{t-1}=y_{t}-\hat{y}_{t|t-1}$ 是时间 $t$ 的残差。

训练数据误差导致在 $t=1,\dots,T$ 的整个平滑过程中对估计的水平进行调整。例如，如果时间 $t$ 的误差为负，则 $y_t < \hat{y}_{t|t-1}$，因此时间 $t-1$ 的水平被高估了。新的水平 $\ell_{t}$ 则是向下调整的先前水平 $\ell_{t-1}$。$\alpha$ 越接近一，水平的估计就越“粗糙”（发生大的调整）。$\alpha$ 越小，水平就越“平滑”（发生小的调整）。

我们也可以写成 $y_t = \ell_{t-1} + e_t$，这样每个观测值都可以用先前的水平加上一个误差来表示。要将其转化为创新状态空间模型，我们只需要指定 $e_t$ 的概率分布。对于加法误差模型，我们假设残差（一步训练误差）$e_t$ 服从均值为 0、方差为 $\sigma^2$ 的正态分布白噪声。其简写符号为 $e_t = \varepsilon_t\sim\text{NID}(0,\sigma^2)$；NID 代表“正态独立分布”（normally and independently distributed）。

那么模型的方程可以写成

我们将 (2) 称为测量（或观测）方程，将 (3) 称为状态（或转移）方程。这两个方程，连同误差的统计分布，构成一个完全指定的统计模型。具体来说，这些构成了简单指数平滑背后的创新状态空间模型。

“创新”（innovations）一词源于所有方程都使用相同的随机误差过程 $\varepsilon_t$。出于同样的原因，这种表述也称为“单一误差源”模型。存在我们在此不介绍的其他多重误差源表述。

观测方程显示了观测值和未观测状态之间的关系。在这种情况下，观测值 $y_t$ 是水平项 $\ell_{t-1}$、$y_t$ 的可预测部分和误差项 $\varepsilon_t$（$y_t$ 的不可预测部分）的线性函数。对于其他创新状态空间模型，这种关系可能是非线性的。

状态方程显示了状态随时间的演变。平滑参数 $\alpha$ 的影响与之前讨论的方法相同。例如，$\alpha$ 控制相继水平项的变化量：$\alpha$ 值高允许水平项快速变化；$\alpha$ 值低导致平滑变化。如果 $\alpha=0$，则序列的水平项不随时间变化；如果 $\alpha=1$，模型简化为随机游走模型，$y_t=y_{t-1}+\varepsilon_t$。

​

ETS(M,N,N)：带有乘法误差的简单指数平滑

类似地，我们可以通过将一步超前训练误差写成相对误差来指定带有乘法误差的模型

$\varepsilon_t = \frac{y_t-\hat{y}_{t|t-1}}{\hat{y}_{t|t-1}}$

其中 $\varepsilon_t \sim \text{NID}(0,\sigma^2)$。将 $\hat{y}_{t|t-1}=\ell_{t-1}$ 代入可得 $y_t = \ell_{t-1}+\ell_{t-1}\varepsilon_t$ 和 $e_t = y_t - \hat{y}_{t|t-1} = \ell_{t-1}\varepsilon_t$。

然后我们可以将状态空间模型的乘法形式写成

​

ETS(A,A,N)：带有加法误差的 Holt 线性方法

对于此模型，我们假设一步超前训练误差由下式给出

$\varepsilon_t=y_t-\ell_{t-1}-b_{t-1} \sim \text{NID}(0,\sigma^2)$

将其代入 Holt 线性方法的误差校正方程，我们得到

其中为简单起见，我们设置了 $\beta=\alpha \beta^*$。

​

ETS(M,A,N)：带有乘法误差的 Holt 线性方法

将一步超前训练误差指定为相对误差，使得

$\varepsilon_t=\frac{y_t-(\ell_{t-1}+b_{t-1})}{(\ell_{t-1}+b_{t-1})}$

并遵循与上述类似的方法，指定了带有乘法误差的 Holt 线性方法基础上的创新状态空间模型，如下所示

其中再次设置 $\beta=\alpha \beta^*$ 且 $\varepsilon_t \sim \text{NID}(0,\sigma^2)$。

​

估计 ETS 模型

除了最小化平方误差和来估计参数外，另一种方法是最大化“似然”。似然是指数据由指定模型产生的概率。因此，大的似然与好的模型相关联。对于加法误差模型，最大化似然（假设误差符合正态分布）得到的结果与最小化平方误差和相同。然而，对于乘法误差模型，将得到不同的结果。在本节中，我们将通过最大化似然来估计平滑参数 $\alpha, \beta, \gamma$ 和 $\phi$ 以及初始状态 $\ell_0, b_0, s_0,s_{-1},\dots,s_{-m+1}$（包括残差方差）。

平滑参数可能的取值受到限制。传统上，参数被限制在 0 到 1 之间，以便这些方程可以被解释为加权平均。也就是说，$0< \alpha,\beta^*,\gamma^*,\phi<1$。对于状态空间模型，我们设置了 $\beta=\alpha\beta^*$ 和 $\gamma=(1-\alpha)\gamma^*$。因此，传统限制转化为 $0< \alpha <1, 0 < \beta < \alpha$ 和 $0< \gamma < 1-\alpha$。实践中，阻尼参数 $\phi$ 通常会受到更严格的限制，以防止模型估计出现数值困难。

另一种看待参数的方式是通过考虑状态空间模型的数学性质。对参数进行限制是为了防止遥远的过去观测值持续影响当前的预测值。这导致了一些参数的可容许约束，这些约束通常（但并非总是）比传统约束区域宽松（Hyndman 等，2008，第 149-161 页）。例如，对于 ETS(A,N,N) 模型，传统参数区域是 $0< \alpha <1$，但可容许区域是 $0< \alpha <2$。对于 ETS(A,A,N) 模型，传统参数区域是 $0<\alpha<1$ 和 $0<\beta<\alpha$，但可容许区域是 $0<\alpha<2$ 和 $0<\beta<4-2\alpha$。

​

模型选择

ETS 统计框架的一个巨大优势是可以使用信息准则进行模型选择。这里可以使用 AIC、AIC_c 和 BIC 来确定哪种 ETS 模型最适合给定的时间序列。

对于 ETS 模型，赤池信息准则（AIC）定义为

$\text{AIC} = -2\log(L) + 2k,$

其中 $L$ 是模型的似然，$k$ 是已估计参数和初始状态的总数（包括残差方差）。

修正小样本偏差后的 AIC（AIC_c）定义为

$AIC_c = AIC + \frac{2k(k+1)}{T-k-1}$

贝叶斯信息准则（BIC）是

$\text{BIC} = \text{AIC} + k[\log(T)-2]$

(误差、趋势、季节) 的三种组合可能导致数值困难。具体来说，可能导致不稳定性的模型是 ETS(A,N,M)、ETS(A,A,M) 和 ETS(A,Ad,M)，这是由于状态方程中存在除以可能接近零的值。在选择模型时，我们通常不考虑这些特定的组合。

带有乘法误差的模型在数据严格为正时非常有用，但当数据包含零或负值时，数值稳定性较差。因此，如果时间序列不严格为正，则不考虑乘法误差模型。在这种情况下，只应用六种完全加法模型。

​

加载库和数据

提示

需要 Statsforecast。要安装，请参阅 说明。

接下来，我们导入绘图库并配置绘图风格。

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf

plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#212946',
    'axes.facecolor': '#212946',
    'savefig.facecolor':'#212946',
    'axes.grid': True,
    'axes.grid.which': 'both',
    'axes.spines.left': False,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.spines.bottom': False,
    'grid.color': '#2A3459',
    'grid.linewidth': '1',
    'text.color': '0.9',
    'axes.labelcolor': '0.9',
    'xtick.color': '0.9',
    'ytick.color': '0.9',
    'font.size': 12 }
plt.rcParams.update(dark_style)

from pylab import rcParams
rcParams['figure.figsize'] = (18,7)

​

读取数据

df = pd.read_csv("https://raw.githubusercontent.com/Naren8520/Serie-de-tiempo-con-Machine-Learning/main/Data/Esperanza_vida.csv", usecols=[1,2])
df.head()

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

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

• ds（时间戳）列应为 Pandas 期望的格式，日期理想为 YYYY-MM-DD，时间戳理想为 YYYY-MM-DD HH:MM:SS。

• y（数值）表示我们希望预测的测量值。

df["unique_id"]="1"
df.columns=["ds", "y", "unique_id"]
df.head()

print(df.dtypes)

ds            object
y            float64
unique_id     object
dtype: object

我们需要将 ds 从 object 类型转换为 datetime。

df["ds"] = pd.to_datetime(df["ds"])

​

使用 plot 方法探索数据

使用 StatsForecast 类的 plot 方法绘制一些时间序列。此方法会打印数据集中的随机时间序列，对于基本探索性数据分析 (EDA) 非常有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

自相关图

fig, axs = plt.subplots(nrows=1, ncols=2)

plot_acf(df["y"],  lags=20, ax=axs[0],color="fuchsia")
axs[0].set_title("Autocorrelation");

# Plot
plot_pacf(df["y"],  lags=20, ax=axs[1],color="lime")
axs[1].set_title('Partial Autocorrelation')

plt.show();

​

时间序列分解

如何以及为何分解时间序列？

在时间序列分析中，为了预测新值，了解历史数据非常重要。更正式地说，了解值随时间变化的模式非常重要。许多原因可能导致我们的预测值方向错误。基本上，时间序列由四个分量组成。这些分量的变化导致时间序列模式的变化。这些分量是

• 水平项 (Level)：这是随时间变化的平均主要值。
• 趋势项 (Trend)：趋势项是导致时间序列出现递增或递减模式的值。
• 季节项 (Seasonality)：这是一种短期内出现在时间序列中的周期性事件，导致时间序列短期出现递增或递减模式。
• 残差/噪声 (Residual/Noise)：这些是时间序列中的随机变化。

这些分量随时间组合形成时间序列。大多数时间序列由水平项和噪声/残差组成，趋势项或季节项是可选值。

如果季节项和趋势项是时间序列的一部分，那么它们将对预测值产生影响。因为预测时间序列的模式可能与之前的时间序列不同。

时间序列中分量的组合可以是两种类型：* 加法 (Additive) * 乘法 (Multiplicative)

​

加法时间序列

如果时间序列的分量相加构成时间序列，则称为加法时间序列。通过可视化，我们可以说，如果时间序列的递增或递减模式在整个序列中相似，则该时间序列是加法的。任何加法时间序列的数学函数都可以表示为：$y(t) = level + Trend + seasonality + noise$

​

乘法时间序列

如果时间序列的分量相乘，则称为乘法时间序列。对于可视化，如果时间序列随时间呈指数增长或下降，则该时间序列可以被视为乘法时间序列。乘法时间序列的数学函数可以表示为。

$y(t) = Level * Trend * seasonality * Noise$

from statsmodels.tsa.seasonal import seasonal_decompose 
a = seasonal_decompose(df["y"], model = "add", period=1)
a.plot();

将时间序列分解为其组成部分有助于我们识别正在分析的时间序列的行为。此外，它还有助于我们了解可以应用的模型类型。以我们的预期寿命数据集为例，我们可以观察到我们的时间序列全年呈现上升趋势，另一方面，也可以观察到该时间序列没有季节性。

通过查看上一张图并了解每个组成部分，我们可以知道可以应用哪种模型： * 存在趋势 * 没有季节性

​

将数据分为训练集和测试集

让我们将数据分成数据集

1. 用于训练模型的数据。
2. 用于测试模型的数据。

对于测试数据，我们将使用最后 6 年的数据来测试和评估模型的性能。

train = df[df.ds<='2013-01-01'] 
test = df[df.ds>'2013-01-01']

train.shape, test.shape

((54, 3), (6, 3))

sns.lineplot(train,x="ds", y="y", label="Train")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.show()

​

使用 StatsForecast 实现 AutoETS

from statsforecast import StatsForecast
from statsforecast.models import AutoETS

​

实例化模型

sf = StatsForecast(models=[AutoETS(model="AZN")], freq='YS')

​

拟合模型

sf.fit(df=train)

StatsForecast(models=[AutoETS])

​

模型预测

y_hat = sf.predict(h=6)
y_hat

sf.plot(train, y_hat)

让我们为预测添加置信区间。

y_hat = sf.predict(h=6, level=[80,90,95])
y_hat

sf.plot(train, y_hat, level=[95])

​

预测方法

内存高效的指数平滑预测。

此方法避免了对象存储带来的内存负担。它类似于不存储信息的 fit_predict。它假定您事先知道预测范围。

y_hat = sf.forecast(df=train, h=6, fitted=True)
y_hat

​

样本内预测

访问拟合的指数平滑样本内预测。

sf.forecast_fitted_values()

​

模型评估

现在我们将使用预测结果评估模型，我们将使用不同类型的指标 MAE、MAPE、MASE、RMSE、SMAPE 来评估准确性。

from functools import partial

import utilsforecast.losses as ufl
from utilsforecast.evaluation import evaluate

evaluate(
    y_hat.merge(test),
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=1), ufl.rmse, ufl.smape],
    train_df=train,
)

​

参考资料

1. Nixtla 自动预测
2. Rob J. Hyndman 和 George Athanasopoulos (2018)。“预测：原理与实践 (第 3 版)”.