AutoCES模型


# 此单元格将不会被渲染，但用于隐藏警告并限制显示的行数。

import warnings
warnings.filterwarnings("ignore")

import logging
logging.getLogger('statsforecast').setLevel(logging.ERROR)

import pandas as pd
pd.set_option('display.max_rows', 6)

使用Statsforecast的AutoCES模型的逐步指南。

本文的目的是提供一个关于如何使用AutoCES和Statsforecast构建CES模型的逐步指南。

在这个过程中，我们将熟悉主要的StatsForecast类以及一些相关的方法，如StatsForecast.plot、StatsForecast.forecast和StatsForecast.cross_validation等。

目录

引言

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

复合指数平滑

方法与模型

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

\[ \begin{equation} {\hat{y}}_t=\alpha {y}_{t-1}+\left(1-\alpha \right){\hat{y}}_{t-1} \tag{1} \end{equation} \]

其中，\(\alpha\)是平滑参数，\({\hat{y}}_t\)是系列的估计值。如果我们将\({\hat{y}}_{t-1}\)用公式（1）替代，且使用一个索引而不是（Brown, 1956），同样的方法可以表示为过去实际观察值的加权平均：

\[ \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} \]

这种表示法的想法是展示权重\(\alpha {\left(1-\alpha \right)}^{j-1}\)在我们样本中是如何随时间分配的。如果平滑参数\(\alpha \in \left(0,1\right)\)，那么权重会随着的增加而指数下降。如果它位于所谓的“可接受范围”内（Brenner等，1968），即\(\alpha \in \left(0,2\right)\)，则权重会以振荡方式下降。在实践和学术文献中，传统和可接受范围都得到了有效的使用（对于后者的应用，参见例如Gardner & Diaz-Saiz，2008；Snyder等，2017）。然而，在现实生活中，权重的分配可能更加复杂，呈现出谐波而非指数下降，这意味着一些过去的观察可能比最近的观察更重要。为了实现这种权重分配，我们在（2）的基础上，通过将（2）中的实变量替换为复变量，引入复杂的动态交互。首先，我们将\({y}_{t-j}\)替换为复变量\({y}_{t-j}+{ie}_{t-j}\)，其中\({e}_t\)是模型的误差项，\(i\)是满足方程\({i}^2=-1\)的虚数单位。这样做的想法是考虑过去每个观察值的实际值和误差对最终预测的影响。其次，我们将\(\alpha\)替换为复变量\({\alpha}_0+i{\alpha}_1\)，并将1替换为\(1+i\)以引入谐波下降的权重。根据复平滑参数的值，权重分配将在时间上展现出多种路径，包括指数型、振荡型和谐波型。最后，两个复数相乘的结果将是另一个复数，因此我们将\({\hat{y}}_{t-j}\)替换为\({\hat{y}}_{t-j}+i{\hat{e}}_{t-j}\)，其中\({\hat{e}}_{t-j}\)是误差项的代理。得到的CES可以写为：

\[ \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} \]

在得到具有谐波分布的权重的模型后，我们现在可以通过替换简化形式

\[{\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}}\]

在（3）中得到：

\[ \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} \]

请注意，\({\hat{e}}_t\) 在时间序列分析和预测目的上并不是很有趣，但它被用作承载方法先前错误信息的载体。在（4）中使用复数变量而不是实数变量，允许对实际值和预测误差进行指数加权。通过改变\({\alpha}_0+i{\alpha}_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背后所依据的统计模型。这种模型可以写成以下状态空间形式：

\[ \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} \]

其中，\({\epsilon}_t\) 是白噪声误差项，\({l}_t\) 是水平成分，\({c}_t\) 是观察值的非线性趋势成分。请注意，时间序列中的依赖关系具有交互结构，并且该时间序列中没有显式的趋势成分，因为此模型不需要像ETS那样将序列人为分解为水平和趋势。尽管我们将\({c}_t\)成分称为“非线性趋势”，但它并不对应于传统的趋势成分，因为它包含了前一个\({c}_{t-1}\)和水平\({l}_{t-1}\)的相关信息。此外，请注意，在（6）中我们使用\({\epsilon}_t\)而不是\({e}_t\)，这意味着只有在没有错误指定的情况下，CES才将（6）视为基础统计模型。在估计该模型的情况下，\({\epsilon}_t\)将被\({e}_t\)替代，随后将引导我们回到原始的公式（4）。

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

\[ \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} \]

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

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

\[\mathbf{g}=\left(\begin{array}{c}{\alpha}_0-{\alpha}_1\\ {}{\alpha}_0+{\alpha}_1\end{array}\right)\] 是持久性向量和

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

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

（6）中的误差项是可加的，因此CES的似然函数是平凡的，类似于加法指数平滑模型中的一种（Hyndman等，2008，第68页）：

\[ \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} \]

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

CES的平稳性和稳定性条件

为了理解CES的性质，我们需要研究其平稳性和稳定性条件。当所有特征值位于单位圆内时，平稳性对状态空间形式的通用指数平滑（8）成立（Hyndman等，2008，第38页）。CES可以是平稳的，也可以不是，这取决于复杂平滑参数的值，与总是非平稳的ETS模型形成对比。计算CES的特征值可以得到以下根：

\[ \begin{equation} \lambda =\frac{2-{\alpha}_0\pm \sqrt{\alpha_0^2+4{\alpha}_1-4}}{2}. \tag 9 \end{equation} \]

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

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

\[ \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）来限制CES，我们只是展示了模型将如何依据复杂平滑参数的值而表现。这一CES性质意味着它能够建模平稳或非平稳过程，而无需在两者之间切换。CES对每个独立时间序列的特性取决于平滑参数的值。

另一个重要属性是来自（7）的稳定性条件对于CES。设\(\epsilon_t={y}_t-{l}_{t-1}\)，可以得到以下结果：

\[ \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} \]

矩阵 \[\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)\] 被称为折现矩阵，可以写成一般形式：

\[ \begin{equation} \mathbf{D}=\mathbf{F}-\mathbf{g}{\mathbf{w}}^{\prime }. \tag {12} \end{equation} \]

如果（12）的所有特征值都位于单位圆内部，则模型被称为稳定。这比模型的平稳性条件更重要，因为它确保复杂权重随着时间的推移而下降，并且较旧的观察值的权重小于新的观察值，这也是传统ETS模型的主要特征之一。特征值由以下公式给出：

\[ \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将是稳定的：

\[ \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} \]

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

条件均值和方差的CES

已知\({l}_t\)和\({c}_t\)的情况下，CES的条件均值可以通过状态空间模型（6）计算：

\[ \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} \]

其中

\[\mathrm{E}\left({y}_{t+h}\mid {\mathbf{v}}_t\right)={\hat{y}}_{t+h}\]

而\(\mathbf{F}\)和\(\mathbf{w}\)是（7）中的矩阵。

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

  1. 当\({\alpha}_1=1\)时，所有预测值将等于最后获得的预测值，这对应于一条平坦的线。这条轨迹在图2A中展示。
  2. 当\({\alpha}_1>1\)时，模型产生的轨迹具有指数增长，如图2B所示。
  3. 当\(\frac{4-{\alpha}_0^2}{4}<{\alpha}_1<1\)时，轨迹变得平稳，CES产生的指数下降如图2C所示。
  4. 当\(1-{\alpha}_0<{\alpha}_1<\frac{4-{\alpha}_0^2}{4}\)时，轨迹变得谐波并将收敛于零，如图2D所示。
  5. 最后，当\(0<{\alpha}_1<1-{\alpha}_0\)时，产生发散的谐波轨迹，模型变得非平稳。这条轨迹在预测中没有用处，因此我们不在图表中显示它。

通过（7），已知\({l}_t\)和\({c}_t\)的情况下，CES的条件方差可以类似于纯加性ETS模型（Hyndman等，2008年，第96页）进行计算。

加载库和数据

Tip

将需要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()
year value
0 1960-01-01 69.123902
1 1961-01-01 69.760244
2 1962-01-01 69.149756
3 1963-01-01 69.248049
4 1964-01-01 70.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()
ds y unique_id
0 1960-01-01 69.123902 1
1 1961-01-01 69.760244 1
2 1962-01-01 69.149756 1
3 1963-01-01 69.248049 1
4 1964-01-01 70.311707 1

现在，让我们使用 .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, engine="matplotlib")

自相关图

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_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();

时间序列的分解

如何分解时间序列以及为什么要这样做？

在时间序列分析中，为了预测新值，了解过去的数据是非常重要的。更正式地说，了解数值随时间变化所遵循的模式是非常重要的。导致我们预测值偏离正确方向的原因可能有很多。基本上，时间序列由四个组成部分组成。这些组成部分的变化导致时间序列模式的变化。这些组成部分包括：

  • 水平： 这是随着时间变化的主要平均值。
  • 趋势： 趋势是导致时间序列中出现增长或下降模式的值。
  • 季节性： 这是在时间序列中短时间内发生的周期性事件，导致时间序列中出现短期增长或下降模式。
  • 残差/噪声： 这些是时间序列中的随机变化。

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

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

时间序列组件的组合可以分为两种类型： * 加法 * 乘法

加法时间序列

如果时间序列的组件是相加形成时间序列的，那么该时间序列称为加法时间序列。通过可视化，我们可以说如果时间序列的增加或减少模式在整个序列中相似，则该时间序列是加法的。任何加法时间序列的数学函数可以表示为： \[y(t) = 水平 + 趋势 + 季节性 + 噪声\]

乘法时间序列

如果时间序列的组件是彼此相乘的，那么该时间序列称为乘法时间序列。通过可视化，如果时间序列随着时间的推移呈现指数增长或下降，则该时间序列可以视为乘法时间序列。乘法时间序列的数学函数可以表示为： \[y(t) = 水平 * 趋势 * 季节性 * 噪声\]

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

要了解AutoCES模型的函数参数，下面列出了相关信息。有关更多信息，请访问文档

model : str
    控制状态空间方程。
season_length : int
    每单位时间的观察次数。例如：24小时数据。
alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    计算符合预测区间的信息。
    默认情况下，模型将计算本地预测区间。

导入库

from statsforecast import StatsForecast
from statsforecast.models import AutoCES

from statsforecast.arima import arima_string

实例化模型

导入并实例化模型。设置参数有时会很棘手。Rob Hyndmann大师的这篇关于季节周期的文章对于season_length的理解会很有帮助。

注意

自动选择最佳的 复杂指数平滑模型，使用信息准则。默认使用赤池信息准则（AICc），而特定模型则采用最大似然法估计。状态空间方程可以根据它们的 \(S\) 简单、\(P\) 部分、\(Z\) 优化或 \(N\) 省略的组件来确定。模型字符串参数定义了CES模型的类型：\(N\) 表示简单CES（无季节性），\(S\) 表示简单季节性（滞后CES），\(P\) 表示部分季节性（无复杂部分），\(F\) 表示完整季节性（滞后CES，包含真实和复杂的季节部分）。

如果选定的组件为 \(Z\)，则作为占位符，要求AutoCES模型找出最佳参数。

season_length = 1 # 年度数据 
horizon = len(test) # 预测数量

# 我们称之为即将使用的模型
models = [AutoCES(season_length=season_length)]

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

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

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

  • n_jobs: n_jobs: int，用于并行处理的工作数量，使用 -1 表示使用所有核心。

  • fallback_model: 如果模型失败，将使用的模型。

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

sf = StatsForecast(df=train,
                   models=models,
                   freq='YS', 
                   n_jobs=-1)

拟合模型

sf.fit()
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
residual Model
0 -0.727729
1 0.144552
2 -0.762086
... ...
51 -0.073258
52 -0.234578
53 0.109990

54 rows × 1 columns

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

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

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

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

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

plt.show();

预测方法

如果您希望在生产环境中提高速度，尤其是在处理多个序列或模型时，我们建议使用 StatsForecast.forecast 方法，而不是 .fit 和 .predict。

主要区别在于 .forecast 不会存储拟合值，并且在分布式环境中具有较高的可扩展性。

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

  • h (int): 表示预测未来 h 步。在这种情况下，指的是未来12个月。

  • level (list of floats): 这个可选参数用于概率预测。设置预测区间的水平（或置信百分位数）。例如，level=[90] 表示模型预期真实值在该区间内的概率为90%。

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

# 预测
Y_hat = sf.forecast(horizon, fitted=True)

Y_hat
ds CES
unique_id
1 2014-01-01 82.906075
1 2015-01-01 83.166687
1 2016-01-01 83.424744
1 2017-01-01 83.685760
1 2018-01-01 83.946213
1 2019-01-01 84.208359 
values=sf.forecast_fitted_values()
values.head()
ds y CES
unique_id
1 1960-01-01 69.123901 69.851631
1 1961-01-01 69.760246 69.615692
1 1962-01-01 69.149757 69.911842
1 1963-01-01 69.248047 69.657822
1 1964-01-01 70.311707 69.601196 
StatsForecast.plot(values)

添加95%的置信区间与预测方法

sf.forecast(h=horizon, level=[95])
ds CES CES-lo-95 CES-hi-95
unique_id
1 2014-01-01 82.906075 82.342484 83.454018
1 2015-01-01 83.166687 82.604027 83.717270
1 2016-01-01 83.424744 82.858574 83.975868
1 2017-01-01 83.685760 83.118942 84.239578
1 2018-01-01 83.946213 83.376907 84.501137
1 2019-01-01 84.208359 83.637741 84.765411 
Y_hat=Y_hat.reset_index()
Y_hat
unique_id ds CES
0 1 2014-01-01 82.906075
1 1 2015-01-01 83.166687
2 1 2016-01-01 83.424744
3 1 2017-01-01 83.685760
4 1 2018-01-01 83.946213
5 1 2019-01-01 84.208359 
# 将预测结果与真实值合并
test['unique_id'] = test['unique_id'].astype(int)
Y_hat = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat
ds y unique_id CES
0 2014-01-01 83.090244 1 82.906075
1 2015-01-01 82.543902 1 83.166687
2 2016-01-01 83.243902 1 83.424744
3 2017-01-01 82.946341 1 83.685760
4 2018-01-01 83.346341 1 83.946213
5 2019-01-01 83.197561 1 84.208359 
fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat]).set_index('ds')
plot_df[['y', "CES"]].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Year ', fontsize=20)
ax.set_xlabel('Timestamp [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)

预测方法及置信区间

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

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

  • h (int): 表示预测在未来 h 步的值。在此情况下，代表提前 12 个月。

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

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

此步骤应在 1 秒内完成。

sf.predict(h=horizon)
ds CES
unique_id
1 2014-01-01 82.906075
1 2015-01-01 83.166687
1 2016-01-01 83.424744
1 2017-01-01 83.685760
1 2018-01-01 83.946213
1 2019-01-01 84.208359 
forecast_df = sf.predict(h=horizon, level=[95]) 

forecast_df
ds CES CES-lo-95 CES-hi-95
unique_id
1 2014-01-01 82.906075 82.342484 83.454018
1 2015-01-01 83.166687 82.604027 83.717270
1 2016-01-01 83.424744 82.858574 83.975868
1 2017-01-01 83.685760 83.118942 84.239578
1 2018-01-01 83.946213 83.376907 84.501137
1 2019-01-01 84.208359 83.637741 84.765411

我们可以使用pandas函数pd.concat()将预测结果与历史数据连接起来，然后可以使用这个结果进行绘图。

pd.concat([df, forecast_df]).set_index('ds')
y unique_id CES CES-lo-95 CES-hi-95
ds
1960-01-01 69.123902 1 NaN NaN NaN
1961-01-01 69.760244 1 NaN NaN NaN
1962-01-01 69.149756 1 NaN NaN NaN
... ... ... ... ... ...
2017-01-01 NaN NaN 83.685760 83.118942 84.239578
2018-01-01 NaN NaN 83.946213 83.376907 84.501137
2019-01-01 NaN NaN 84.208359 83.637741 84.765411

66 rows × 5 columns

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

def plot_forecasts(y_hist, y_true, y_pred, models):
    _, ax = plt.subplots(1, 1, figsize = (20, 7))
    y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
    df_plot = pd.concat([y_hist, y_true]).set_index('ds').tail(12*10)
    df_plot[['y'] + models].plot(ax=ax, linewidth=2)
    colors = ['orange', 'black', 'green']
    for model, color in zip(models, colors):
        ax.fill_between(df_plot.index, 
                        df_plot[f'{model}-lo-95'], 
                        df_plot[f'{model}-hi-95'],
                        alpha=.35,
                        color=color,
                        label=f'{model}-level-90')
    ax.set_title('', fontsize=22)
    ax.set_ylabel('', fontsize=20)
    ax.set_xlabel('Timestamp [t]', fontsize=20)
    ax.legend(prop={'size': 15})
    ax.grid(True)
plot_forecasts(train, test, forecast_df, models=['CES'])

让我们使用Statsforecast中提供的plot函数绘制相同的图，如下所示。

sf.plot(df, forecast_df, level=[95])

交叉验证

在之前的步骤中，我们利用历史数据来预测未来。然而，为了评估预测的准确性，我们也希望了解模型在过去的表现。为了评估模型在数据上的准确性和鲁棒性，进行交叉验证是必要的。

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

以下图表展示了这种交叉验证策略：

执行时间序列交叉验证

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

在这种情况下，我们希望评估每个模型在最后5个月的表现（n_windows=5），每隔一个月进行一次预测（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: 索引。如果您不喜欢使用索引，可以运行crossvalidation_df.resetindex()
  • ds: 日期戳或时间索引
  • cutoff: n_windows的最后一个日期戳或时间索引。
  • y: 真实值
  • "model": 包含模型名称和拟合值的列。
crossvalidation_df.head()
ds cutoff y CES
unique_id
1 1984-01-01 1983-01-01 75.389511 74.952705
1 1985-01-01 1983-01-01 75.470734 75.161736
1 1986-01-01 1983-01-01 75.770729 75.377945
1 1987-01-01 1983-01-01 76.219513 75.590378
1 1988-01-01 1983-01-01 76.370735 75.806343

模型评估

现在我们可以使用适当的准确性指标来计算预测的准确性。在这里，我们将使用均方根误差（RMSE）。为此，我们首先需要安装datasetsforecast，这是一个由Nixtla开发的Python库，其中包含计算RMSE的功能。

%%capture
!pip install datasetsforecast
from datasetsforecast.losses import rmse

计算RMSE的函数接受两个参数：

  1. 实际值。
  2. 预测值，在本例中是AutoCES。
rmse = rmse(crossvalidation_df['y'], crossvalidation_df["CES"])
print("RMSE using cross-validation: ", rmse)
RMSE using cross-validation:  0.40604722

如您所注意到的，我们已经使用交叉验证结果来评估我们的模型。

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

from datasetsforecast.losses import (mae, mape, mase, rmse, smape)
def evaluate_performace(y_hist, y_hat, model):
    evaluation = {}
    evaluation[model] = {}
    for metric in [mase, mae, mape, rmse, smape]:
        metric_name = metric.__name__
        if metric_name == 'mase':
            evaluation[model][metric_name] = metric(y_hat['y'].values, 
                                                y_hat[model].values, 
                                                y_hist['y'].values, seasonality=24)
        else:
            evaluation[model][metric_name] = metric(y_hat['y'].values, y_hat[model].values)
    return pd.DataFrame(evaluation).T
evaluate_performace(train, Y_hat, model="CES")
mae mape mase rmse smape
CES 0.556314 0.669916 0.087037 0.630183 0.667133

鸣谢

我们要感谢 Naren Castellon 为编写本教程所做的贡献。

参考文献

  1. Nixtla-AutoCES

  2. Rob J. Hyndman 和 George Athanasopoulos (2018)。“预测原则与实践，时间序列交叉验证”。

  3. 复杂指数平滑

Give us a ⭐ on Github