​
目录

​
引言

指数平滑法一直是组织中用于支持各种决策的最流行的预测方法之一，应用于库存管理、排程、收入管理等领域。尽管其相对简单和透明使其在研究和实践中非常有吸引力，但识别潜在趋势仍然具有挑战性，并对最终精度产生重大影响。这导致了各种趋势模型的修改，引入了模型选择问题。为了解决这个问题，我们基于复变函数理论提出了复杂指数平滑法（CES）。基本的CES方法只涉及两个参数，不需要模型选择过程。尽管进行了这些简化，CES证明与现有方法相比具有竞争力，甚至更优越。我们展示了CES相对于传统指数平滑模型具有几个优势：它可以对平稳和非平稳过程进行建模和预测，并且CES可以捕捉到传统指数平滑分类中定义的水平和趋势情况。CES在几个预测竞赛数据集上进行了评估，显示出比现有基准更好的性能。我们得出结论，CES具有时间序列建模的理想特性，并为研究开辟了新的有前景的途径。

​
复杂指数平滑法

​
方法与模型

利用时间序列的复值表示，我们类比传统的指数平滑方法提出了CES。考虑简单指数平滑法

y^t=αyt−1+(1−α)y^t−1 \begin{equation} {\hat{y}}_t=\alpha {y}_{t-1}+\left(1-\alpha \right){\hat{y}}_{t-1} \tag{1} \end{equation} (1)​

其中 α\alpha 是平滑参数，而 y^t{\hat{y}}_t 是序列的估计值。如果我们将 y^t−1{\hat{y}}_{t-1} 替换为带有一个索引的公式 (1)，则同样的方法可以表示为先前实际观测值的加权平均值（Brown, 1956）。

y^t=α∑j=1t−1(1−α)j−1yt−j \begin{equation} {\hat{y}}_t=\alpha \sum \limits_{j=1}^{t-1}{\left(1-\alpha \right)}^{j-1}{y}_{t-j} \tag{2} \end{equation}

这种表示的目的是展示权重 α(1−α)j−1\alpha {\left(1-\alpha \right)}^{j-1} 在我们的样本中如何随时间分布。如果平滑参数 α∈(0,1)\alpha \in \left(0,1\right) 那么权重随时间的增加呈指数衰减。如果它位于所谓的“允许界限”（Brenner et al., 1968）内，即 α∈(0,2)\alpha \in \left(0,2\right) 则权重以振荡方式衰减。传统界限和允许界限在实践和学术文献中都得到了有效应用（后者的应用例如见 Gardner & Diaz-Saiz, 2008; Snyder et al., 2017）。然而，在现实生活中，权重的分布可能更复杂，呈现谐波而非指数衰减，这意味着一些过去的观测值可能比最近的观测值更重要。为了实现这种权重分布，我们基于公式(2)，通过将公式(2)中的实变量替换为复变量，引入复杂的动态相互作用。首先，我们将 yt−j{y}_{t-j} 替换为复变量 yt−j+iet−j{y}_{t-j}+{ie}_{t-j}，其中 et{e}_t 是模型的误差项，而 ii 是虚数单位（满足方程 i2=−1{i}^2=-1）。这背后的想法是让实际值和误差在过去对每个观测值的影响作用于最终预测。其次，我们将 α\alpha 替换为复变量 α0+iα1{\alpha}_0+i{\alpha}_1，并将 1 替换为 1+i1+i，以引入谐波衰减权重。根据复杂平滑参数的值，权重分布将随时间呈现各种轨迹，包括指数、振荡和谐波。最后，两个复数相乘的结果将是另一个复数，因此我们将 y^t−j{\hat{y}}_{t-j} 替换为 y^t−j+ie^t−j{\hat{y}}_{t-j}+i{\hat{e}}_{t-j}，其中 e^t−j{\hat{e}}_{t-j} 是误差项的代理。由此获得的CES可以写成

y^t+ie^t=(α0+iα1)∑j=1t−1(1+i−(α0+iα1))j−1(yt−j+iet−j) \begin{equation} {\hat{y}}_t+i{\hat{e}}_t=\left({\alpha}_0+i{\alpha}_1\right)\sum \limits_{j=1}^{t-1}{\left(1+i-\left({\alpha}_0+i{\alpha}_1\right)\right)}^{j-1}\left({y}_{t-j}+{ie}_{t-j}\right) \tag{3} \end{equation}

在公式 (3) 中通过替换，可以得到一个更短的形式

y^t−1+ie^t−1=(α0+iα1)∑j=2t−1(1+i−(α0+iα1))j−1(yt−j+iet−j){\displaystyle \begin{array}{cc}& {\hat{y}}_{t-1}+i{\hat{e}}_{t-1}\\ {}& \kern1em =\left({\alpha}_0+i{\alpha}_1\right)\sum \limits_{j=2}^{t-1}{\left(1+i-\left({\alpha}_0+i{\alpha}_1\right)\right)}^{j-1}\left({y}_{t-j}+{ie}_{t-j}\right)\end{array}}

得到

y^t+ie^t=(α0+iα1)(yt−1+iet−1)+(1−α0+i−iα1)(y^t−1+ie^t−1). \begin{equation} {\displaystyle \begin{array}{cc}{\hat{y}}_t+i{\hat{e}}_t& =\left({\alpha}_0+i{\alpha}_1\right)\left({y}_{t-1}+{ie}_{t-1}\right)\\ {}& \kern1em +\left(1-{\alpha}_0+i-i{\alpha}_1\right)\left({\hat{y}}_{t-1}+i{\hat{e}}_{t-1}\right).\end{array}} \tag 4 \end{equation}

注意，对于时间序列分析和预测目的来说，e^t{\hat{e}}_t 并不重要，但它被用作一个容器，包含有关该方法先前误差的信息。在公式(4)中使用复变量而不是实变量，可以获取实际值和预测误差的指数加权值。通过改变 α0+iα1{\alpha}_0+i{\alpha}_1 的值，我们可以调节实际值和预测误差的哪些比例应该用于未来的预测。

将复值函数表示为两个实值函数的系统，得到

y^t=(α0yt−1+(1−α0)y^t−1)−(α1et−1+(1−α1)e^t−1)e^t=(α1yt−1+(1−α1)y^t−1)+(α0et−1+(1−α0)e^t−1). \begin{equation} {\displaystyle \begin{array}{ll}& {\hat{y}}_t=\left({\alpha}_0{y}_{t-1}+\left(1-{\alpha}_0\right){\hat{y}}_{t-1}\right)-\left({\alpha}_1{e}_{t-1}+\left(1-{\alpha}_1\right){\hat{e}}_{t-1}\right)\\ {}& {\hat{e}}_t=\left({\alpha}_1{y}_{t-1}+\left(1-{\alpha}_1\right){\hat{y}}_{t-1}\right)+\left({\alpha}_0{e}_{t-1}+\left(1-{\alpha}_0\right){\hat{e}}_{t-1}\right).\end{array}} \tag 5 \end{equation}

CES 引入了实部和虚部之间的交互。方程 (6) 通过彼此的先前值连接，导致随时间推移的交互，由复平滑参数值定义。

但该方法本身具有限制性，不易产生预测区间和推导似然函数。理解 CES 背后的统计模型也很重要。该模型可以用以下状态空间形式书写：

yt=lt−1+ϵtlt=lt−1−(1−α1)ct−1+(α0−α1)ϵtct=lt−1+(1−α0)ct−1+(α0+α1)ϵt, \begin{equation} {\displaystyle \begin{array}{ll}& {y}_t={l}_{t-1}+{\epsilon}_t\\ {}& {l}_t={l}_{t-1}-\left(1-{\alpha}_1\right){c}_{t-1}+\left({\alpha}_0-{\alpha}_1\right){\epsilon}_t\\ {}& {c}_t={l}_{t-1}+\left(1-{\alpha}_0\right){c}_{t-1}+\left({\alpha}_0+{\alpha}_1\right){\epsilon}_t,\end{array}} \tag 6 \end{equation}

其中 ϵt{\epsilon}_t 是白噪声误差项，lt{l}_t 是水平分量，ct{c}_t 是观测时的非线性趋势分量。请注意，时间序列中的依赖关系具有交互结构，时间序列中不存在显式趋势分量，因为该模型不需要像 ETS 那样将序列人为地分解为水平和趋势。尽管我们将 ct{c}_t 分量称为“非线性趋势”，但它与传统趋势分量不同，因为它包含先前 ct−1{c}_{t-1} 和水平 lt−1{l}_{t-1} 的信息。此外，请注意，我们在 (6) 中使用 ϵt{\epsilon}_t 而不是 et{e}_t，这意味着 CES 仅在没有误设定误差时才将 (6) 作为其潜在统计模型。在估计此模型的情况下，ϵt{\epsilon}_t 将被 et{e}_t 替代，这将引导我们回到原始公式 (4)。

这个想法允许以更简短通用的方式重写 (6)，类似于一般的单误差源 (SSOE) 状态空间框架：

yt=w′vt−1+ϵtvt=Fvt−1+gϵt, \begin{equation} {\displaystyle \begin{array}{ll}& {y}_t={\mathbf{w}}^{\prime }{\mathbf{v}}_{t-1}+{\epsilon}_t\\ {}& {\mathbf{v}}_t={\mathbf{Fv}}_{t-1}+\mathbf{g}{\epsilon}_t,\end{array}} \tag 7 \end{equation}

其中 vt=(ltct){\mathbf{v}}_t=\left(\begin{array}{c}{l}_t\\ {}{c}_t\end{array}\right) 是状态向量

F=(1−(1−α1)11−α0)\mathbf{F}=\left(\begin{array}{cc}1& -\left(1-{\alpha}_1\right)\\ {}1& 1-{\alpha}_0\end{array}\right) 是转移矩阵

g=(α0−α1α0+α1)\mathbf{g}=\left(\begin{array}{c}{\alpha}_0-{\alpha}_1\\ {}{\alpha}_0+{\alpha}_1\end{array}\right) 是持久性向量，且

w=(10)\mathbf{w}=\left(\begin{array}{c}1\\ {}0\end{array}\right) 是测量向量。

状态空间形式 (7) 允许以类似于 ETS 的方式扩展 CES，以包含季节性或外生变量的附加状态。模型 (7) 与传统 ETS 的主要区别在于，(7) 中的转移矩阵包含平滑参数，而这并非 ETS 模型的标准特征。此外，持久性向量包含复平滑参数的相互作用，而不是平滑参数本身。

(6) 中的误差项是加性的，因此 CES 的似然函数是简单的，与加性指数平滑模型中的类似 (Hyndman 等人，2008 年，第 68 页)：

L(g,v0,σ2∣Y)=(1σ2π)Texp(−12∑t=1T(ϵtσ)2), \begin{equation} \mathrm{\mathcal{L}}\left(\mathbf{g},{\mathbf{v}}_0,{\sigma}^2\mid \mathbf{Y}\right)={\left(\frac{1}{\sigma \sqrt{2\pi }}\right)}^T\exp \left(-\frac{1}{2}\sum \limits_{t=1}^T{\left(\frac{\epsilon_t}{\sigma}\right)}^2\right), \tag 8 \end{equation}

其中 v0{\mathbf{v}}_0 是初始状态向量，σ2{\sigma}^2 是误差项的方差，Y\mathbf{Y} 是所有样本内观测值的向量。

​
CES 的平稳性和稳定性条件

为了理解 CES 的特性，我们需要研究其平稳性和稳定性条件。前者对于状态空间形式 (8) 中的一般指数平滑成立，当所有特征值都位于单位圆内时 (Hyndman 等人，2008 年，第 38 页)。与总是非平稳的 ETS 模型不同，CES 可以是平稳的或非平稳的，这取决于复平滑参数值。计算 CES 的特征值得到以下根：

λ=2−α0±α02+4α1−42. \begin{equation} \lambda =\frac{2-{\alpha}_0\pm \sqrt{\alpha_0^2+4{\alpha}_1-4}}{2}. \tag 9 \end{equation}

如果两个根的绝对值都小于 1，则估计的 CES 是平稳的。

当 α1>1{\alpha}_1>1 时，其中一个特征值将始终大于一。在这种情况下，两个特征值都是实数，CES 会产生非平稳轨迹。当 α1=1{\alpha}_1=1 时，CES 等同于 ETS(A,N,N)。最后，当模型变得平稳时，

{α1<5−2α0α1<1α1>1−α0 \begin{equation} \left\{\begin{array}{l}{\alpha}_1<5-2{\alpha}_0\\ {}{\alpha}_1<1\\ {}{\alpha}_1>1-{\alpha}_0\end{array}\right. \tag {10} \end{equation} (10)​

请注意，我们并未用条件 (10) 来限制 CES，我们仅展示该模型将如何根据复平滑参数的值来表现。CES 的这一特性意味着它能够对平稳或非平稳过程进行建模，而无需在它们之间切换。对于每个单独的时间序列，CES 的属性取决于平滑参数的值。

从 (7) 得出的另一个重要属性是 CES 的稳定性条件。使用 ϵt=yt−lt−1\epsilon_t={y}_t-{l}_{t-1} 得到以下结果

yt=lt−1+ϵt(ltct)=(1−α0+α1−(1−α1)1−α0−α11−α0)(lt−1ct−1)+(α0−α1α1+α0)yt. \begin{equation} {\displaystyle \begin{array}{ll}{y}_t& ={l}_{t-1}+{\epsilon}_t\\ {}\left(\begin{array}{c}{l}_t\\ {}{c}_t\end{array}\right)& =\left(\begin{array}{cc}1-{\alpha}_0+{\alpha}_1& -\left(1-{\alpha}_1\right)\\ {}1-{\alpha}_0-{\alpha}_1& 1-{\alpha}_0\end{array}\right)\left(\begin{array}{c}{l}_{t-1}\\ {}{c}_{t-1}\end{array}\right)\\ {}& \kern1em +\left(\begin{array}{c}{\alpha}_0-{\alpha}_1\\ {}{\alpha}_1+{\alpha}_0\end{array}\right){y}_t.\end{array}} \tag {11} \end{equation} (11)​

矩阵 D=(1−α0+α1−(1−α1)1−α0−α11−α0)\mathbf{D}=\left(\begin{array}{cc}1-{\alpha}_0+{\alpha}_1& -\left(1-{\alpha}_1\right)\\ {}1-{\alpha}_0-{\alpha}_1& 1-{\alpha}_0\end{array}\right) 称为折现矩阵，可以写成一般形式

D=F−gw′. \begin{equation} \mathbf{D}=\mathbf{F}-\mathbf{g}{\mathbf{w}}^{\prime }. \tag {12} \end{equation}

如果 (12) 的所有特征值都位于单位圆内，则称该模型稳定。这比模型的平稳性条件更重要，因为它确保了复权重随时间衰减，并且旧观测值的权重小于新观测值，这是传统 ETS 模型的主要特征之一。特征值由以下公式给出

λ=2−2α0+α1±8α1+4α0−4α0α1−4−3α122. \begin{equation} \lambda =\frac{2-2{\alpha}_0+{\alpha}_1\pm \sqrt{8{\alpha}_1+4{\alpha}_0-4{\alpha}_0{\alpha}_1-4-3{\alpha}_1^2}}{2}. \tag {13} \end{equation}

当以下不等式组满足时，CES 是稳定的

{(α0−2.5)2+α12>1.25(α0−0.5)2+(α1−1)2>0.25(α0−1.5)2+(α1−0.5)2<1.5. \begin{equation} \left\{\begin{array}{l}{\left({\alpha}_0-2.5\right)}^2+{\alpha}_1^2>1.25\\ {}{\left({\alpha}_0-0.5\right)}^2+{\left({\alpha}_1-1\right)}^2>0.25\\ {}{\left({\alpha}_0-1.5\right)}^2+{\left({\alpha}_1-0.5\right)}^2<1.5\end{array}.\right. \tag {14} \end{equation} (14)​

平稳区域和稳定区域均显示在 图 1 中。平稳区域 (10) 对应于三角形。三角形内曲线下方所有平滑参数的组合将产生平稳谐波轨迹，而其余组合则产生指数轨迹。稳定性条件 (14) 对应于暗区域。稳定区域与平稳区域相交，但一般来说，稳定的 CES 可以产生平稳和非平稳预测。

​
CES 的条件均值和方差

使用状态空间模型 (6) 可以计算已知 lt{l}_t 和 ct{c}_t 时，CES 向前 hh 步的条件均值

E(yt+h∣vt)=w′Fh−1vt, \begin{equation} \mathrm{E}\left({y}_{t+h}\mid {\mathbf{v}}_t\right)={\mathbf{w}}^{\prime }{\mathbf{F}}^{h-1}{\mathbf{v}}_t, \tag {15} \end{equation}

其中

E(yt+h∣vt)=y^t+h\mathrm{E}\left({y}_{t+h}\mid {\mathbf{v}}_t\right)={\hat{y}}_{t+h}

而 F\mathbf{F} 和 w\mathbf{w} 是来自 (7) 的矩阵。

(15) 的预测轨迹将根据 lt,ct{l}_t, {c}_t 的值以及复平滑参数而有所不同。对平稳性条件的分析表明，CES 的预测轨迹有几种类型，具体取决于复平滑参数的特定值。

  1. 当 α1=1{\alpha}_1=1 时，所有预测值都将等于最后获得的预测值，这对应于一条平线。该轨迹显示在 图 2A 中。
  2. 当 α1>1{\alpha}_1>1 时，模型产生的轨迹呈指数增长，如 图 2B 所示。
  3. 当 4−α024<α1<1\frac{4-{\alpha}_0^2}{4}<{\alpha}_1<1 时，轨迹变为平稳，并且 CES 产生指数衰减，如 图 2C 所示。
  4. 当 1−α0<α1<4−α0241-{\alpha}_0<{\alpha}_1<\frac{4-{\alpha}_0^2}{4} 轨迹变为谐波并收敛到零 图 2D。
  5. 最后，当 0<α1<1−α00<{\alpha}_1<1-{\alpha}_0 产生发散的谐波轨迹时，模型变为非平稳。此轨迹在预测中无用，因此我们在图表上不显示它。

使用 (7) 计算 CES 的条件方差，用于 hh 步长预测，已知 lt{l}_t 和 ct{c}_t，可以与纯加性 ETS 模型类似地计算 (Hyndman et al., 2008, p. 96)。

​
加载库和数据

提示

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

接下来，我们导入绘图库并配置绘图样式。

import pandas as pd

import scipy.stats as stats

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()
年值
01960-01-0169.123902
11961-01-0169.760244
21962-01-0169.149756
31963-01-0169.248049
41964-01-0170.311707

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()
dsyunique_id
01960-01-0169.1239021
11961-01-0169.7602441
21962-01-0169.1497561
31963-01-0169.2480491
41964-01-0170.3117071

现在，我们使用 .tail() 函数查看时间序列的最后几行。

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");

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

#plt.savefig("Gráfico de Densidad y qq")
plt.show();

​
时间序列分解

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

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

  • 水平 (Level): 这是随时间平均的主要值。
  • 趋势 (Trend): 趋势是导致时间序列中出现增加或减少模式的值。
  • 季节性 (Seasonality): 这是时间序列中短期发生的周期性事件，导致时间序列中出现短期增加或减少模式。
  • 残差/噪声 (Residual/Noise): 这些是时间序列中的随机变动。

随着时间的推移，这些分量的组合形成了时间序列。大多数时间序列包含水平和噪声/残差，趋势或季节性是可选值。

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

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

​
加法时间序列

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

​
乘法时间序列

如果时间序列的分量相乘，则该时间序列称为乘法时间序列。在可视化方面，如果时间序列随时间呈现指数增长或下降，则可以将其视为乘法时间序列。乘法时间序列的数学函数可以表示为。

y(t)=Level∗Trend∗seasonality∗Noisey(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. 用于测试模型的数据。

对于测试数据，我们将使用最近的 12 个月来测试和评估我们模型的性能。

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 实现 AutoCES

​
加载库

from statsforecast import StatsForecast
from statsforecast.models import AutoCES

​
实例化模型

导入并实例化模型。设置 参数有时很棘手。这篇关于 季节周期）由大师 Rob Hyndmann 撰写，可能会有所帮助 season_length

注意

使用信息准则自动选择最佳的 Complex Exponential Smoothing model。默认为 Akaike Information Criterion (AICc)，而特定模型使用 maximum likelihood 进行估计。状态空间方程可以根据其 SS 简单、PP 部分、ZZ 优化或 NN 忽略的分量来确定。模型字符串参数定义了 CES 模型的类型：NN 用于简单 CES（无季节性）、SS 用于简单季节性（滞后 CES）、PP 用于部分季节性（无复杂部分）、FF 用于完整季节性（包含实部和复部季节性分量的滞后 CES）。

如果将分量选择为 ZZ，它充当一个占位符，要求 AutoCES 模型找出最佳参数。

season_length = 1 # year data 
horizon = len(test) # number of predictions

# We call the model that we are going to use
models = [AutoCES(season_length=season_length)]

我们通过实例化新的 StatsForecast 对象并使用以下参数来拟合模型

models: 模型列表。从模型中选择所需的模型并导入它们。

  • freq: 一个字符串，指示数据的频率。（请参阅 pandas 可用的频率。）

  • n_jobs: n_jobs: int, 并行处理中使用的作业数，使用 -1 表示所有核心。

  • fallback_model: 当模型失败时使用的模型。

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

sf = StatsForecast(models=models, freq='YS')

​
拟合模型

sf.fit(df=train)

StatsForecast(models=[CES])

result=sf.fitted_[0,0].model_
print(result.keys())
print(result['fit'])

dict_keys(['loglik', 'aic', 'bic', 'aicc', 'mse', 'amse', 'fit', 'fitted', 'residuals', 'm', 'states', 'par', 'n', 'seasontype', 'sigma2', 'actual_residuals'])
results(x=array([1.63706552, 1.00511519]), fn=76.78049826760919, nit=27, simplex=array([[1.63400329, 1.00510199],
       [1.63706552, 1.00511519],
       [1.63638944, 1.00512037]]))

现在让我们可视化模型的残差。

如上所示，获得的结果以字典形式输出，为了从字典中提取每个元素，我们将使用 .get() 函数提取元素，然后将其保存在 pd.DataFrame() 中。

residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual
残差模型
0-0.727729
10.144552
2-0.762086
……
51-0.073258
52-0.234578
530.109990
fig, axs = plt.subplots(nrows=2, ncols=2)

# plot[1,1]
residual.plot(ax=axs[0,0])
axs[0,0].set_title("Residuals");

# plot
sns.distplot(residual, ax=axs[0,1]);
axs[0,1].set_title("Density plot - Residual");

# plot
stats.probplot(residual["residual Model"], dist="norm", plot=axs[1,0])
axs[1,0].set_title('Plot Q-Q')

# plot
plot_acf(residual,  lags=35, ax=axs[1,1],color="fuchsia")
axs[1,1].set_title("Autocorrelation");

plt.show();

​
预测方法

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

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

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

  • h (int): 表示未来 h 步的预测。在本例中，提前 12 个月。

  • level (list of floats): 这个可选参数用于概率预测。设置预测区间的置信水平（或置信百分位数）。例如，level=[90] 意味着模型预计实际值有 90% 的时间会落在此区间内。

这里的 forecast 对象是一个新的数据框，包含一列用于存放模型名称和 y hat 值，以及用于不确定性区间的列。根据您的计算机性能，此步骤大约需要 1 分钟。（如果想将速度提高到几秒，请移除 AutoModels，例如 ARIMA 和 Theta）

# Prediction
Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat
unique_iddsCES
012014-01-0182.906075
112015-01-0183.166687
212016-01-0183.424744
312017-01-0183.685760
412018-01-0183.946213
512019-01-0184.208359
values=sf.forecast_fitted_values()
values.head()
unique_iddsyCES
011960-01-0169.12390269.851631
111961-01-0169.76024469.615692
211962-01-0169.14975669.911842
311963-01-0169.24804969.657822
411964-01-0170.31170769.601196
StatsForecast.plot(values)

使用 forecast 方法添加 95% 置信区间

sf.forecast(df=train, h=horizon, level=[95])
unique_iddsCESCES-lo-95CES-hi-95
012014-01-0182.90607582.34248383.454016
112015-01-0183.16668782.60402983.717271
212016-01-0183.42474482.85857383.975870
312017-01-0183.68576083.11894684.239582
412018-01-0183.94621383.37690584.501133
512019-01-0184.20835983.63773884.765408
# Merge the forecasts with the true values
Y_hat = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat
dsyunique_idCES
02014-01-0183.090244182.906075
12015-01-0182.543902183.166687
22016-01-0183.243902183.424744
32017-01-0182.946341183.685760
42018-01-0183.346341183.946213
52019-01-0183.197561184.208359
sf.plot(train, Y_hat)

​
带置信区间的 predict 方法

要生成预测，请使用 predict 方法。

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

  • h (int): 表示未来 h 步的预测。在本例中，提前 12 个月。

  • level (list of floats): 这个可选参数用于概率预测。设置预测区间的置信水平（或置信百分位数）。例如，level=[95] 意味着模型预计实际值有 95% 的时间会落在此区间内。

这里的 forecast 对象是一个新的数据框，包含一列用于存放模型名称和 y hat 值，以及用于不确定性区间的列。

此步骤应耗时不到 1 秒。

sf.predict(h=horizon)
unique_iddsCES
012014-01-0182.906075
112015-01-0183.166687
212016-01-0183.424744
312017-01-0183.685760
412018-01-0183.946213
512019-01-0184.208359
forecast_df = sf.predict(h=horizon, level=[95]) 
forecast_df
unique_iddsCESCES-lo-95CES-hi-95
012014-01-0182.90607582.34248383.454016
112015-01-0183.16668782.60402983.717271
212016-01-0183.42474482.85857383.975870
312017-01-0183.68576083.11894684.239582
412018-01-0183.94621383.37690584.501133
512019-01-0184.20835983.63773884.765408

现在让我们可视化预测结果和时间序列的历史数据，同时绘制在进行 95% 置信度预测时获得的置信区间。

sf.plot(train, test.merge(forecast_df), level=[95])

​
交叉验证

在之前的步骤中，我们使用历史数据预测未来。然而，为了评估其准确性，我们还想知道模型在过去的表现如何。为了评估模型在数据上的准确性和鲁棒性，请执行交叉验证。

对于时间序列数据，通过定义一个滑动窗口在历史数据上预测其后的时期来完成交叉验证。这种形式的交叉验证使我们能够在更广泛的时间实例范围内更好地估计模型的预测能力，同时保持训练集中的数据连续，这是我们的模型所要求的。

下图描绘了这种交叉验证策略

​
执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现都非常慢。statsforecast 库将交叉验证实现为分布式操作，从而减少了执行过程的时间消耗。如果您有大数据集，还可以使用 Ray、Dask 或 Spark 在分布式集群中执行交叉验证。

在本例中，我们希望评估每个模型在过去 5 个月内的性能 (n_windows=5)，以 12 个月为步长进行预测 (step_size=12)。根据您的计算机性能，此步骤大约需要 1 分钟。

来自 StatsForecast 类中的 cross_validation 方法接受以下参数。

  • df: 训练数据框

  • h (int): 表示未来 h 步的预测。在本例中，提前 12 个月。

  • step_size (int): 每个窗口之间的步长。换句话说：您希望多久运行一次预测过程。

  • n_windows(int): 用于交叉验证的窗口数量。换句话说：您想评估过去多少个预测过程。

crossvalidation_df = sf.cross_validation(df=train,
                                         h=horizon,
                                         step_size=12,
                                         n_windows=3)

crossvaldation_df 对象是一个新的数据框，包含以下列

  • unique_id: 序列标识符
  • ds: 日期戳或时间索引
  • cutoff: n_windows 的最后一个日期戳或时间索引。
  • y: 真实值
  • "model": 包含模型名称和拟合值的列。
crossvalidation_df.head()
unique_iddscutoffyCES
011984-01-011983-01-0175.38951274.952705
111985-01-011983-01-0175.47073275.161736
211986-01-011983-01-0175.77073275.377945
311987-01-011983-01-0176.21951275.590378
411988-01-011983-01-0176.37073275.806343

​
模型评估

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

from functools import partial

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

evaluate(
    Y_hat,
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=season_length), ufl.rmse, ufl.smape],
    train_df=train,
)
unique_id指标CES
01mae0.556314
11mape0.006699
21mase1.770512
31rmse0.630183
41smape0.003336

​
参考资料

  1. Nixtla-AutoCES
  2. Rob J. Hyndman and George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”.
  3. Ivan Svetunkov, Nikolaos Kourentzes, John Keith Ord, “Complex exponential smoothing”