霍尔特模型


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

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 的逐步指南。

在这个过程中特别介绍 StatsForecast 类及一些相关的方法，比如 StatsForecast.plot、StatsForecast.forecast 和 StatsForecast.cross_validation 等等。

介绍

霍尔特模型，也称为双指数平滑法，是一种广泛用于时间序列分析的预测技术。它由查尔斯·霍尔特在1957年开发，作为对布朗简单指数平滑法的改进。

霍尔特模型用于预测表现出趋势的时间序列的未来值。该模型使用两个平滑参数，一个用于估计趋势，另一个用于估计时间序列的水平或基线。这些参数分别称为\(\alpha\)和\(\beta\)。

霍尔特模型是布朗简单指数平滑法的扩展，后者仅使用一个平滑参数来估计时间序列的趋势和基线水平。霍尔特模型通过添加第二个平滑参数来提高预测的准确性。

霍尔特模型的主要优点之一是易于实施，并且不需要大量历史数据就能生成准确的预测。此外，该模型高度适应性强，可以定制以适应各种时间序列。

然而，霍尔特模型也存在一些局限性。例如，该模型假设时间序列是平稳的，并且趋势是线性的。如果时间序列不是平稳的或具有非线性趋势，则霍尔特模型可能不是最合适的。

总体而言，霍尔特模型是时间序列分析中一种有用且广泛使用的技术，尤其是在预计系列会表现出线性趋势的情况下。

霍尔特方法

简单指数平滑在数据有趋势时效果不佳。在这种情况下，我们可以使用双指数平滑。与其他方法相比，这是一种处理不含季节性趋势数据的更可靠的方法。此方法在公式中添加了一个时间趋势方程。使用两个不同的权重或平滑参数同时更新这两个组件。

霍尔特的指数平滑有时也称为双指数平滑。这里的主要思想是使用SES并将其提升以捕捉趋势组件。

霍尔特（1957）将简单指数平滑扩展到可以预测具有趋势的数据。该方法包括一个预测方程和两个平滑方程（一个用于水平，一个用于趋势）：

假设一个序列具有以下特征：

水平
趋势
无季节性
噪声

\[\begin{align*} \text{预测方程}&& \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} \\ \text{水平方程} && \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ \text{趋势方程} && b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)b_{t-1}, \end{align*}\]

其中 \(\ell_{t}\) 表示在时间 \(t\) 对序列水平的估计，\(b_t\) 表示在时间 \(t\) 对序列趋势（斜率）的估计，\(\alpha\) 是水平的平滑参数，\(0\le\alpha\le1\)，而 \(\beta^{*}\) 是趋势的平滑参数，\(0\le\beta^*\le1\)。

与简单指数平滑一样，这里的水平方程显示 \(\ell_{t}\) 是观测值 \(y_{t}\) 和时间 \(t\) 的一步预测训练值 \(\ell_{t-1} + b_{t-1}\) 的加权平均。趋势方程显示 \(b_t\) 是基于 \(\ell_{t} - \ell_{t-1}\) 和趋势的先前估计 \(b_{t-1}\) 对时间 \(t\) 估计趋势的加权平均。

预测函数不再是平坦的，而是有趋势的。\(h\) 步提前的预测等于最后估计的水平加上 \(h\) 倍的最后估计的趋势值。因此，预测是 \(h\) 的线性函数。

创新状态空间模型用于指数平滑

表7.6中介绍的指数平滑方法是生成点预测的算法。本教程中的统计模型生成相同的点预测，但也可以生成预测（或预报）区间。统计模型是一种随机（或随机的）数据生成过程，可以产生整个预测分布。

每个模型由一个测量方程组成，用于描述观察到的数据，以及一些状态方程，用于描述未观察到的组件或状态（水平、趋势、季节性）随时间的变化。因此，这些被称为状态空间模型。

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

为了区分加性误差模型和乘性误差模型，我们将每个状态空间模型标记为ETS( .,.,.)表示（误差、趋势、季节性）。这个标签也可以被视为指数平滑。使用与表7.5中相同的符号，组件的可能性为：\(Error=\{A,M\}\)，\(Trend=\{N,A,A_d\}\) 和 \(Seasonal=\{N,A,M\}\)。

在我们的案例中，即具有趋势的线性霍尔特模型，我们将看到两个情况，包括加性和乘性。

ETS(A,A,N): 霍尔特线性方法与加性误差

对于该模型，我们假设一步预测的训练误差为 \(\varepsilon_t=y_t-\ell_{t-1}-b_{t-1} \sim \text{NID}(0,\sigma^2)\)。将其代入霍尔特线性方法的误差修正方程，得到

\[\begin{align*} y_t&=\ell_{t-1}+b_{t-1}+\varepsilon_t\\ \ell_t&=\ell_{t-1}+b_{t-1}+\alpha \varepsilon_t\\ b_t&=b_{t-1}+\beta \varepsilon_t, \end{align*}\]

为了简单起见，我们设定 \(\beta=\alpha \beta^*\)。

ETS(M,A,N): 霍尔特线性方法与乘性误差

将一步预测的训练误差指定为相对误差，形式为

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

并采用与上述类似的方法，霍尔特线性方法的乘性误差基础上的创新状态空间模型可指定为

\[\begin{align*} y_t&=(\ell_{t-1}+b_{t-1})(1+\varepsilon_t)\\ \ell_t&=(\ell_{t-1}+b_{t-1})(1+\alpha \varepsilon_t)\\ b_t&=b_{t-1}+\beta(\ell_{t-1}+b_{t-1}) \varepsilon_t \end{align*}\]

其中同样有 \(\beta=\alpha \beta^*\) 和 \(\varepsilon_t \sim \text{NID}(0,\sigma^2)\)。

指数平滑方法的分类

基于时间序列成分的思想，我们可以转向ETS分类法。ETS代表“误差-趋势-季节性”，定义了成分之间的具体互动方式。根据误差、趋势和季节性的类型，Pegels（1969年）提出了一种分类方法，随后由Hyndman等人（2002年）进一步发展，并在Hyndman等人（2008年）中进行了修订。根据这一分类方法，误差、趋势和季节性可以是：

误差：“加性”（A）或“乘性”（M）；
趋势：“无”（N）、“加性”（A）、“加性衰减”（Ad）、“乘性”（M）或“乘性衰减”（Md）；
季节性：“无”（N）、“加性”（A）或“乘性”（M）。

ETS分类法中的成分具有明确的解释：水平表示每个时间段的平均值，趋势反映了值的变化，而季节性对应于周期性波动（例如，每年一月销售额的增长）。根据上述成分类型，理论上可以设计出30种不同误差、趋势和季节性的ETS模型。图1展示了具有确定性（随时间不变）水平、趋势、季节性以及加性误差项的不同时间序列示例。

“图1：对应于加性误差ETS模型的时间序列” 图4.1：对应于加性误差ETS模型的时间序列

从图1中的图形中需要注意的事项：

当季节性为乘性时，其幅度随着数据水平的提高而增加，而当季节性为加性时，幅度保持不变。例如，比较ETS(A,A,A)与ETS(A,A,M)：对于前者，第一年和最后年之间的最高点和最低点之间的距离大致相同。而在ETS(A,A,M)的情况下，这个距离随着系列水平的提高而增加；
当趋势为乘性时，数据表现出指数增长/衰减；
衰减趋势减缓了加性和乘性趋势；
如果系列水平不变，几乎不可能区分加性和乘性季节性，因为在两种情况下，季节性的幅度都将保持不变（比较ETS(A,N,A)和ETS(A,N,M)）。

ETS 分类法中的数学模型

我希望读者可以更加清楚地了解到 ETS 框架是建立在时间序列分解的理念之上的。通过引入不同的组件，定义它们的类型，并添加它们更新的方程式，我们可以构建更能捕捉时间序列关键特征的模型。但我们也应该考虑组件随时间可能发生的变化。“转移”或“状态”方程应反映这一变化：它们解释了水平、趋势或季节性组件是如何演变的。

正如第 2.2 节所讨论的，考虑到不同类型的组件及其相互作用，我们最终得到了分类法中的 30 个模型。表 1 和表 2 数学上总结了图 1 和图 2 中图示的所有 30 个 ETS 模型，展示了测量和转移方程的公式。

表 1：加性误差 ETS 模型 | | 非季节性 |加性 |乘性| |—-|———–|———–|————–| |无趋势|\(\begin{aligned} &y_{t} = l_{t-1} + \epsilon_t \\ &l_t = l_{t-1} + \alpha \epsilon_t \end{aligned}\) |\(\begin{aligned} &y_{t} = l_{t-1} + s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} + \alpha \epsilon_t \\ &s_t = s_{t-m} + \gamma \epsilon_t \end{aligned}\) |\(\begin{aligned} &y_{t} = l_{t-1} s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} + \alpha \frac{\epsilon_t}{s_{t-m}} \\ &s_t = s_{t-m} + \gamma \frac{\epsilon_t}{l_{t-1}} \end{aligned}\)| |加性| \(\begin{aligned} &y_{t} = l_{t-1} + b_{t-1} + \epsilon_t \\ &l_t = l_{t-1} + b_{t-1} + \alpha \epsilon_t \\ &b_t = b_{t-1} + \beta \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} + b_{t-1} + s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} + b_{t-1} + \alpha \epsilon_t \\ &b_t = b_{t-1} + \beta \epsilon_t \\ &s_t = s_{t-m} + \gamma \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} + b_{t-1}) s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} + b_{t-1} + \alpha \frac{\epsilon_t}{s_{t-m}} \\ &b_t = b_{t-1} + \beta \frac{\epsilon_t}{s_{t-m}} \\ &s_t = s_{t-m} + \gamma \frac{\epsilon_t}{l_{t-1} + b_{t-1}} \end{aligned}\)| |加性衰减| \(\begin{aligned} &y_{t} = l_{t-1} + \phi b_{t-1} + \epsilon_t \\ &l_t = l_{t-1} + \phi b_{t-1} + \alpha \epsilon_t \\ &b_t = \phi b_{t-1} + \beta \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} + \phi b_{t-1} + s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} + \phi b_{t-1} + \alpha \epsilon_t \\ &b_t = \phi b_{t-1} + \beta \epsilon_t \\ &s_t = s_{t-m} + \gamma \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} + \phi b_{t-1}) s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} + \phi b_{t-1} + \alpha \frac{\epsilon_t}{s_{t-m}} \\ &b_t = \phi b_{t-1} + \beta \frac{\epsilon_t}{s_{t-m}} \\ &s_t = s_{t-m} + \gamma \frac{\epsilon_t}{l_{t-1} + \phi b_{t-1}} \end{aligned}\)| |乘性| \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1} + \epsilon_t \\ &l_t = l_{t-1} b_{t-1} + \alpha \epsilon_t \\ &b_t = b_{t-1} + \beta \frac{\epsilon_t}{l_{t-1}} \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1} + s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} b_{t-1} + \alpha \epsilon_t \\ &b_t = b_{t-1} + \beta \frac{\epsilon_t}{l_{t-1}} \\ &s_t = s_{t-m} + \gamma \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1} s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} b_{t-1} + \alpha \frac{\epsilon_t}{s_{t-m}} \\ &b_t = b_{t-1} + \beta \frac{\epsilon_t}{l_{t-1}s_{t-m}} \\ &s_t = s_{t-m} + \gamma \frac{\epsilon_t}{l_{t-1} b_{t-1}} \end{aligned}\)| |乘性衰减| \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1}^\phi + \epsilon_t \\ &l_t = l_{t-1} b_{t-1}^\phi + \alpha \epsilon_t \\ &b_t = b_{t-1}^\phi + \beta \frac{\epsilon_t}{l_{t-1}} \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1}^\phi + s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} b_{t-1}^\phi + \alpha \epsilon_t \\ &b_t = b_{t-1}^\phi + \beta \frac{\epsilon_t}{l_{t-1}} \\ &s_t = s_{t-m} + \gamma \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1}^\phi s_{t-m} + \epsilon_t \\ &l_t = l_{t-1} b_{t-1}^\phi + \alpha \frac{\epsilon_t}{s_{t-m}} \\ &b_t = b_{t-1}^\phi + \beta \frac{\epsilon_t}{l_{t-1}s_{t-m}} \\ &s_t = s_{t-m} + \gamma \frac{\epsilon_t}{l_{t-1} b_{t-1}} \end{aligned}\)|

表 2：乘性误差 ETS 模型 | |非季节性 |加性 |乘性| |——|————-|———-|————–| |无趋势| \(\begin{aligned} &y_{t} = l_{t-1}(1 + \epsilon_t) \\ &l_t = l_{t-1}(1 + \alpha \epsilon_t) \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} + s_{t-m})(1 + \epsilon_t) \\ &l_t = l_{t-1} + \alpha \mu_{y,t} \epsilon_t \\ &s_t = s_{t-m} + \gamma \mu_{y,t} \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} s_{t-m}(1 + \epsilon_t) \\ &l_t = l_{t-1}(1 + \alpha \epsilon_t) \\ &s_t = s_{t-m}(1 + \gamma \epsilon_t) \end{aligned}\)| |加性| \(\begin{aligned} &y_{t} = (l_{t-1} + b_{t-1})(1 + \epsilon_t) \\ &l_t = (l_{t-1} + b_{t-1})(1 + \alpha \epsilon_t) \\ &b_t = b_{t-1} + \beta \mu_{y,t} \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} + b_{t-1} + s_{t-m})(1 + \epsilon_t) \\ &l_t = l_{t-1} + b_{t-1} + \alpha \mu_{y,t} \epsilon_t \\ &b_t = b_{t-1} + \beta \mu_{y,t} \epsilon_t \\ &s_t = s_{t-m} + \gamma \mu_{y,t} \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} + b_{t-1}) s_{t-m}(1 + \epsilon_t) \\ &l_t = (l_{t-1} + b_{t-1})(1 + \alpha \epsilon_t) \\ &b_t = b_{t-1} + \beta (l_{t-1} + b_{t-1}) \epsilon_t \\ &s_t = s_{t-m} (1 + \gamma \epsilon_t) \end{aligned}\)| |加性衰减| \(\begin{aligned} &y_{t} = (l_{t-1} + \phi b_{t-1})(1 + \epsilon_t) \\ &l_t = (l_{t-1} + \phi b_{t-1})(1 + \alpha \epsilon_t) \\ &b_t = \phi b_{t-1} + \beta \mu_{y,t} \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} + \phi b_{t-1} + s_{t-m})(1 + \epsilon_t) \\ &l_t = l_{t-1} + \phi b_{t-1} + \alpha \mu_{y,t} \epsilon_t \\ &b_t = \phi b_{t-1} + \beta \mu_{y,t} \epsilon_t \\ &s_t = s_{t-m} + \gamma \mu_{y,t} \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} + \phi b_{t-1}) s_{t-m}(1 + \epsilon_t) \\ &l_t = l_{t-1} + \phi b_{t-1} (1 + \alpha \epsilon_t) \\ &b_t = \phi b_{t-1} + \beta (l_{t-1} + \phi b_{t-1}) \epsilon_t \\ &s_t = s_{t-m}(1 + \gamma \epsilon_t) \end{aligned}\)| |乘性| \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1} (1 + \epsilon_t) \\ &l_t = l_{t-1} b_{t-1} (1 + \alpha \epsilon_t) \\ &b_t = b_{t-1} (1 + \beta \epsilon_t) \end{aligned}\) | \(\begin{aligned} &y_{t} = (l_{t-1} b_{t-1} + s_{t-m})(1 + \epsilon_t) \\ &l_t = l_{t-1} b_{t-1} + \alpha \mu_{y,t} \epsilon_t \\ &b_t = b_{t-1} + \beta \frac{\mu_{y,t}}{l_{t-1}} \epsilon_t \\ &s_t = s_{t-m} + \gamma \mu_{y,t} \epsilon_t \end{aligned}\)| \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1} s_{t-m} (1 + \epsilon_t) \\ &l_t = l_{t-1} b_{t-1} (1 + \alpha \epsilon_t) \\ &b_t = b_{t-1} (1 + \beta \epsilon_t) \\ &s_t = s_{t-m} (1 + \gamma \epsilon_t) \end{aligned}\)| |乘性衰减| \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1}^\phi (1 + \epsilon_t) \\ &l_t = l_{t-1} b_{t-1}^\phi (1 + \alpha \epsilon_t) \\ &b_t = b_{t-1}^\phi (1 + \beta \epsilon_t) \end{aligned}\)| \(\begin{aligned} &y_{t} = (l_{t-1} b_{t-1}^\phi + s_{t-m})(1 + \epsilon_t) \\ &l_t = l_{t-1} b_{t-1}^\phi + \alpha \mu_{y,t} \epsilon_t \\ &b_t = b_{t-1}^\phi + \beta \frac{\mu_{y,t}}{l_{t-1}} \epsilon_t \\ &s_t = s_{t-m} + \gamma \mu_{y,t} \epsilon_t \end{aligned}\) | \(\begin{aligned} &y_{t} = l_{t-1} b_{t-1}^\phi s_{t-m} (1 + \epsilon_t) \\ &l_t = l_{t-1} b_{t-1}^\phi \left(1 + \alpha \epsilon_t\right) \\ &b_t = b_{t-1}^\phi \left(1 + \beta \epsilon_t\right) \\ &s_t = s_{t-m} \left(1 + \gamma \epsilon_t\right) \end{aligned}\)|

从统计学的角度来看，表 1 和表 2 中的公式对应于“真实模型”，它们解释了潜在数据背后的模型，但在构建和估计时，\(\epsilon_t\) 被估计的 \(e_t\) 替代（依赖于误差类型的不同而计算方式不同），时间序列组件和光滑参数也被它们的估计值替代（例如用 \(\hat \alpha\) 替代 \(\alpha\)）。然而，如果这些模型参数的值是已知的，那么就可以通过这些模型生成点预测和条件 h 步的期望。

霍尔特线性趋势方法的属性

霍尔特线性趋势方法是一种时间序列预测技术，利用指数平滑来估计时间序列的水平和趋势成分。该方法具有几个属性，包括：

加性模型：霍尔特线性趋势方法假设时间序列可以分解为一个加性模型，其中观察值是水平、趋势和误差成分的总和。
平滑参数：该方法使用两个平滑参数，α 和 β，用于估计时间序列的水平和趋势成分。这些参数分别控制施加于水平和趋势成分的平滑量。
线性趋势：霍尔特线性趋势方法假设时间序列的趋势成分遵循一条直线。这意味着该方法适用于显示出恒定线性趋势的时间序列数据。
预测：该方法使用估计的水平和趋势成分来预测时间序列的未来值。下一个时期的预测值由水平和趋势成分的总和给出。
优化：平滑参数α和β通过优化过程来估计，目的是最小化预测值与观察值之间的平方误差和。这涉及对不同平滑参数值的迭代，直到找到最佳值。
季节性：霍尔特线性趋势方法可以扩展以纳入季节性成分。这涉及在模型中添加一个季节性成分，以捕捉时间序列中定期发生的任何系统性变化。

总体而言，霍尔特线性趋势方法是一种强大且广泛使用的预测技术，可以用来为具有恒定线性趋势的时间序列数据生成准确的预测。该方法易于实现，并且可以扩展以处理具有季节性变化的时间序列数据。

加载库和数据

Tip

需要使用Statsforecast。有关安装，请参见说明。

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

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.style.use('grayscale') # 五三八 灰度 经典
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#008080',  # #212946
    'axes.facecolor': '#008080',
    'savefig.facecolor': '#008080',
    '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': '#000000',  #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)

读取数据

import pandas as pd

df=pd.read_csv("https://raw.githubusercontent.com/Naren8520/Serie-de-tiempo-con-Machine-Learning/main/Data/ads.csv")
df.head()

	Time	Ads
0	2017-09-13T00:00:00	80115
1	2017-09-13T01:00:00	79885
2	2017-09-13T02:00:00	89325
3	2017-09-13T03:00:00	101930
4	2017-09-13T04:00:00	121630

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	2017-09-13T00:00:00	80115	1
1	2017-09-13T01:00:00	79885	1
2	2017-09-13T02:00:00	89325	1
3	2017-09-13T03:00:00	101930	1
4	2017-09-13T04:00:00	121630	1

print(df.dtypes)

ds           object
y             int64
unique_id    object
dtype: object

我们可以看到我们的时间变量(ds)是以对象格式表示的，我们需要将其转换为日期格式。

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

使用 plot 方法探索数据

使用 StatsForecast 类中的 plot 方法绘制一些系列。该方法从数据集中打印随机系列，适用于基本的探索性数据分析 (EDA)。

from statsforecast import StatsForecast

StatsForecast.plot(df)

增强型迪基-福勒检验

增强型迪基-福勒（ADF）检验是一种统计检验，旨在确定时间序列数据中是否存在单位根。单位根会导致时间序列分析中出现不可预测的结果。在单位根检验中形成的原假设用于确定时间序列数据受趋势影响的程度。接受原假设意味着我们接纳时间序列数据不是平稳的证据。拒绝原假设或接受备择假设，则意味着我们接纳时间序列数据是由平稳过程生成的证据。这个过程也被称为平稳趋势。ADF检验统计量的值是负数。较低的ADF值表示对原假设的更强拒绝。

增强型迪基-福勒检验是一种常见的统计检验，用于检验给定的时间序列是否平稳。我们可以通过定义原假设和备择假设来实现这一点。

原假设：时间序列是非平稳的。它表现出时间依赖性趋势。
备择假设：时间序列是平稳的。换句话说，该序列不依赖于时间。
ADF或t统计量 < 临界值：拒绝原假设，时间序列是平稳的。
ADF或t统计量 > 临界值：未能拒绝原假设，时间序列是非平稳的。

from statsmodels.tsa.stattools import adfuller

def Augmented_Dickey_Fuller_Test_func(series , column_name):
    print (f'Dickey-Fuller test results for columns: {column_name}')
    dftest = adfuller(series, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','No Lags Used','Number of observations used'])
    for key,value in dftest[4].items():
       dfoutput['Critical Value (%s)'%key] = value
    print (dfoutput)
    if dftest[1] <= 0.05:
        print("Conclusion:====>")
        print("Reject the null hypothesis")
        print("The data is stationary")
    else:
        print("Conclusion:====>")
        print("The null hypothesis cannot be rejected")
        print("The data is not stationary")

Augmented_Dickey_Fuller_Test_func(df["y"],'Ads')

Dickey-Fuller test results for columns: Ads
Test Statistic         -7.089634e+00
p-value                 4.444804e-10
No Lags Used            9.000000e+00
                            ...     
Critical Value (1%)    -3.462499e+00
Critical Value (5%)    -2.875675e+00
Critical Value (10%)   -2.574304e+00
Length: 7, dtype: float64
Conclusion:====>
Reject the null hypothesis
The data is stationary

自相关图

自相关函数

定义 1. 设 \(\{x_t;1 ≤ t ≤ n\}\) 为来自 \(\{X_t\}\) 的大小为 n 的时间序列样本。 1. \(\bar x = \sum_{t=1}^n \frac{x_t}{n}\) 称为 \(\{X_t\}\) 的样本均值。 2. \(c_k =\sum_{t=1}^{n−k} (x_{t+k}- \bar x)(x_t−\bar x)/n\) 被称为 \(\{X_t\}\) 的样本自协方差函数。 3. \(r_k = c_k /c_0\) 被称为 \(\{X_t\}\) 的样本自相关函数。

关于这个定义，有以下几点说明：

像大多数文献一样，本指南使用 ACF 来表示样本自相关函数和自相关函数。ACF 所表示的内容可以在上下文中很容易识别。
显然，c0 是 \(\{X_t\}\) 的样本方差。此外，\(r_0 = c_0/c_0 = 1\)，对于任何整数 \(k\)，都有 \(|r_k| ≤ 1\)。
当我们计算任何固定长度为 \(n\) 的样本序列的 ACF 时，对于较大的 \(k\)，我们不能对 \(r_k\) 的值过于自信，因为随着 \(k\) 的增大，用于计算 \(r_k\) 的 \((x_{t +k }, x_t )\) 的对数减少。一个经验法则是不估计 \(r_k\) 对于 \(k > n/3\)，另一个是 \(n ≥ 50, k ≤ n/4\)。在任何情况下，谨慎总是一个好主意。
我们也可以根据定义 1 计算非平稳时间序列样本的 ACF。然而，在这种情况下，随着 \(k\) 增加，ACF 或 \(r_k\) 的衰减速度非常慢或几乎不衰减。
将 ACF \((r_k)\) 绘制与滞后 \(k\) 相对容易，但在分析时间序列样本时非常有帮助。这样的 ACF 图称为自相关图。
如果 \(\{X_t\}\) 是平稳的并且 \(E(X_t)=0\) 且对于所有 \(k \neq 0\)，\(\rho_k =0\)，即，它是一个白噪声序列，那么 \(r_k\) 的采样分布渐近于正态分布，均值为 0，方差为 \(1/n\)。因此，约有 95% 的机会 \(r_k\) 落在区间 \([−1.96/\sqrt{n}, 1.96/\sqrt{n}]\)。

现在我们可以总结如下： (1) 如果时间序列图清楚地显示出趋势或/和季节性，那么它肯定是非平稳的； (2) 如果 ACF \(r_k\) 随着滞后 \(k\) 的增加非常缓慢或几乎不衰减，则该时间序列也应该是非平稳的。

偏自相关

设 \(\{X_t\}\) 为一个平稳时间序列，且 \(E(X_t) = 0\)。这里假设 \(E(X_t) = 0\) 只是为了简洁。如果 \(E(X_t) = \mu \neq 0\)，可以将 \(\{X_t\}\) 替换为 \(\{X_t −\mu \}\)。现在考虑 \(X_t\) 关于 \(\{X_{t−k+1:t−1}\}\) 的线性回归（预测），对于任何整数 \(k ≥ 2\)。我们用 \(\hat X_t\) 表示这个回归（预测）： \[\hat X_t =\alpha_1 X_{t−1}+···+\alpha_{k−1} X_{t−k+1}\]

其中 \(\{\alpha_1, · · · , \alpha_{k−1} \}\) 满足

\[\{\alpha_1, · · · , \alpha_{k−1} \}=\argmin_{\beta_1,···,\beta_{k−1}} E[X_t −(\beta_1 X_{t−1} +···+\beta_{k−1} X_{t−k+1})]^2\]

也就是说，\(\{\alpha_1, · · · , \alpha_{k−1} \}\) 是通过最小化预测的均方误差而选择的。同样，令 \(\hat X_{t −k}\) 表示 \(X_{t −k}\) 关于 \(\{X_{t −k+1:t −1}\}\) 的回归（预测）：

\[\hat X_{t−k} =\eta_1 X_{t−1}+···+\eta_{k−1} X_{t−k+1}\]

注意如果 \(\{X_t\}\) 是平稳的，则 \(\{\alpha_{1:k−1} \} = \{\eta_{1:k−1} \}\). 现在令 \(\hat Z_{t−k} = X_{t−k} − \hat X_{t−k}\) 和 \(\hat Z_t = X_t − \hat X_t\)。那么 \(\hat Z_{t−k}\) 是去除中介变量 \(\{X_{t−k+1:t−1} \}\) 对 \(X_{t−k}\) 影响的残差，而 \(\hat Z_t\) 是去除 \(\{X_{t −k+1:t −1} \}\) 对 \(X_t\) 影响的残差。

定义 2. 平稳时间序列 \(\{X_t \}\) 的滞后 \(k\) 的偏自相关函数 (PACF) 为：

\[\phi_{11} = Corr(X_{t−1}, X_t ) = \frac{Cov(X_{t−1}, X_t )} {[Var(X_{t−1})Var(X_t)]^{1/2}} = \rho_1\] 和

\[\phi_{kk} = Corr(\hat Z_{t−k},\hat Z_t) = \frac{Cov(\hat Z_{t−k},\hat Z_t)} {[Var(\hat Z_{t −k} )Var(\hat Z_t )]^{1/2}}, \ k ≥ 2\]

另一方面，以下定理为估计平稳时间序列的PACF铺平了道路，其证明可以在 Fan 和 Yao (2003) 中找到。

定理 1. 设 \(\{X_t \}\) 为一个平稳时间序列，且 \(E(X_t ) = 0\)，并且 \(\{a_{1k},··· ,a_{kk}\}\) 满足

\[\{a_{1k},··· ,a_{kk}\}= \argmin_{a_1 ,··· ,a_k} E(X_t − a_1 X_{t−1}−···−a_k X_{t−k})^2\]

那么，对于 \(k≥1\)，有 \(\phi_{kk} =a_{kk}\)。

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

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

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

plt.show();

时间序列的分解

如何分解时间序列，以及为什么要分解？

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

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

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

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

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

加法型时间序列

如果时间序列的组件相加以形成时间序列，则该时间序列称为加法型时间序列。通过可视化，我们可以说，如果时间序列的增加或减少模式在整个序列中相似，则时间序列是加法型的。任何加法型时间序列的数学函数可以表示为： \[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 = "additive", period=12)
a.plot();

乘法性

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

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

让我们将数据分成两部分： 1. 用于训练我们的 Holt 模型 的数据。 2. 用于测试我们的模型的数据。

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

train = df[df.ds<='2017-09-20 17:00:00'] 
test = df[df.ds>'2017-09-20 17:00:00']

train.shape, test.shape

((186, 3), (30, 3))

现在让我们绘制训练数据和测试数据。

sns.lineplot(train,x="ds", y="y", label="Train", linestyle="--")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.title("Ads watched (hourly data)");
plt.show()

Holt 方法的实现与 StatsForecast

为了更好地了解 Holt 模型 函数的参数，以下列出了相关信息。有关更多信息，请访问文档。

season_length : int
    每单位时间内的观察次数。例如：12 个月的数据。
error_type : str
    ETS 模型的误差类型。可以是加性 (A) 或乘性 (M)。
alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    用于计算符合预测区间的信息。
    默认情况下，模型将计算本地预测
    区间。

加载库

from statsforecast import StatsForecast
from statsforecast.models import Holt

实例化模型

导入并实例化模型。设置参数有时会很棘手。大师 Rob Hyndmann 的文章季节性周期对于 season_length 可能会有所帮助。

season_length = 24 # 每小时数据 
horizon = len(test) # 预测数量

models = [Holt(season_length=season_length, error_type="A", alias="Add"),
          Holt(season_length=season_length, error_type="M", alias="Multi")]

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

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

freq: 一个字符串，表示数据的频率。（请参阅 pandas 提供的频率。）
n_jobs: n_jobs：整数，表示并行处理使用的工作数量，使用 -1 表示所有核心。
fallback_model: 在模型失败时使用的模型。

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

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

适配模型

sf.fit()

StatsForecast(models=[Add,Multi])

让我们看看我们的 Holt 模型 的结果。我们可以通过以下指令来观察它：

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

dict_keys(['loglik', 'aic', 'bic', 'aicc', 'mse', 'amse', 'fit', 'residuals', 'components', 'm', 'nstate', 'fitted', 'states', 'par', 'sigma2', 'n_params', 'method', 'actual_residuals'])
results(x=array([ 9.99900000e-01,  1.00000000e-04,  7.99488495e+04, -4.74727943e+00]), fn=5213.773789117043, nit=72, simplex=array([[ 9.99900000e-01,  1.00000000e-04,  7.99083958e+04,
         3.48822281e+00],
       [ 9.99900000e-01,  1.00000000e-04,  7.99488495e+04,
        -4.74727943e+00],
       [ 9.99900000e-01,  1.00000000e-04,  8.02890214e+04,
         7.41636110e+00],
       [ 9.99900000e-01,  1.00000000e-04,  8.01679532e+04,
        -1.60232330e+01],
       [ 9.99900000e-01,  1.00000000e-04,  8.03170190e+04,
        -3.36503173e+00]]))

现在让我们可视化我们模型的拟合值。

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

residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual

	residual Model
0	170.897743
1	-225.252721
2	9444.730190
...	...
213	-20317.858920
214	-7924.609136
215	-14867.577350

216 rows × 1 columns

import scipy.stats as stats

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

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 积值列，以及不确定区间的列。根据你的计算机，这一步大约需要 1 分钟。

Y_hat = sf.forecast(horizon, fitted=True)
Y_hat

	ds	Add	Multi
unique_id
1	2017-09-22 00:00:00	80281.781250	81537.476562
1	2017-09-22 01:00:00	80277.085938	82788.335938
1	2017-09-22 02:00:00	80272.382812	84039.195312
...	...	...	...
1	2017-09-23 03:00:00	80154.835938	115310.718750
1	2017-09-23 04:00:00	80150.132812	116561.578125
1	2017-09-23 05:00:00	80145.429688	117812.437500

30 rows × 3 columns

values=sf.forecast_fitted_values()
values.head()

	ds	y	Add	Multi
unique_id
1	2017-09-13 00:00:00	80115.0	79944.101562	79406.656250
1	2017-09-13 01:00:00	79885.0	80110.250000	81393.054688
1	2017-09-13 02:00:00	89325.0	79880.273438	81163.125000
1	2017-09-13 03:00:00	101930.0	89320.250000	90602.976562
1	2017-09-13 04:00:00	121630.0	101926.195312	103208.789062

StatsForecast.plot(values)

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

sf.forecast(h=horizon, level=[95])

	ds	Add	Add-lo-95	Add-hi-95	Multi	Multi-lo-95	Multi-hi-95
unique_id
1	2017-09-22 00:00:00	80281.781250	56698.507812	103865.062500	81537.476562	65208.687500	97866.257812
1	2017-09-22 01:00:00	80277.085938	46925.292969	113628.875000	82788.335938	59458.406250	106118.265625
1	2017-09-22 02:00:00	80272.382812	39423.585938	121121.179688	84039.195312	55172.894531	112905.500000
...	...	...	...	...	...	...	...
1	2017-09-23 03:00:00	80154.835938	-44792.613281	205102.281250	115310.718750	3485.836914	227135.593750
1	2017-09-23 04:00:00	80150.132812	-47015.281250	207315.546875	116561.578125	1745.947998	231377.203125
1	2017-09-23 05:00:00	80145.429688	-49200.355469	209491.218750	117812.437500	1.485114	235623.390625

30 rows × 7 columns

Y_hat=Y_hat.reset_index()
Y_hat

	unique_id	ds	Add	Multi
0	1	2017-09-22 00:00:00	80281.781250	81537.476562
1	1	2017-09-22 01:00:00	80277.085938	82788.335938
2	1	2017-09-22 02:00:00	80272.382812	84039.195312
...	...	...	...	...
27	1	2017-09-23 03:00:00	80154.835938	115310.718750
28	1	2017-09-23 04:00:00	80150.132812	116561.578125
29	1	2017-09-23 05:00:00	80145.429688	117812.437500

30 rows × 4 columns

# 将预测结果与真实值合并
Y_hat1 = pd.concat([df,Y_hat],  keys=['unique_id', 'ds'])
Y_hat1

		ds	y	unique_id	Add	Multi
unique_id	0	2017-09-13 00:00:00	80115.0	1	NaN	NaN
	1	2017-09-13 01:00:00	79885.0	1	NaN	NaN
	2	2017-09-13 02:00:00	89325.0	1	NaN	NaN
...	...	...	...	...	...	...
ds	27	2017-09-23 03:00:00	NaN	1	80154.835938	115310.718750
	28	2017-09-23 04:00:00	NaN	1	80150.132812	116561.578125
	29	2017-09-23 05:00:00	NaN	1	80145.429688	117812.437500

246 rows × 5 columns

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat1]).set_index('ds').tail(300)
plot_df['y'].plot(ax=ax, linewidth=2)
plot_df[ "Add"].plot(ax=ax, linewidth=2)
plot_df[ "Multi"].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Ads watched (hourly data)', fontsize=20)
ax.set_xlabel('Monthly [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	Add	Multi
unique_id
1	2017-09-22 00:00:00	80281.781250	81537.476562
1	2017-09-22 01:00:00	80277.085938	82788.335938
1	2017-09-22 02:00:00	80272.382812	84039.195312
...	...	...	...
1	2017-09-23 03:00:00	80154.835938	115310.718750
1	2017-09-23 04:00:00	80150.132812	116561.578125
1	2017-09-23 05:00:00	80145.429688	117812.437500

30 rows × 3 columns

forecast_df = sf.predict(h=horizon, level=[80,95]) 

forecast_df

	ds	Add	Add-lo-95	Add-lo-80	Add-hi-80	Add-hi-95	Multi	Multi-lo-95	Multi-lo-80	Multi-hi-80	Multi-hi-95
unique_id
1	2017-09-22 00:00:00	80281.781250	56698.507812	64861.507812	95702.062500	103865.062500	81537.476562	65208.687500	70860.656250	92214.289062	97866.257812
1	2017-09-22 01:00:00	80277.085938	46925.292969	58469.519531	102084.648438	113628.875000	82788.335938	59458.406250	67533.710938	98042.953125	106118.265625
1	2017-09-22 02:00:00	80272.382812	39423.585938	53562.789062	106981.976562	121121.179688	84039.195312	55172.894531	65164.535156	102913.851562	112905.500000
...	...	...	...	...	...	...	...	...	...	...	...
1	2017-09-23 03:00:00	80154.835938	-44792.613281	-1543.909790	161853.578125	205102.281250	115310.718750	3485.836914	42192.359375	188429.078125	227135.593750
1	2017-09-23 04:00:00	80150.132812	-47015.281250	-2998.863037	163299.125000	207315.546875	116561.578125	1745.947998	41487.671875	191635.484375	231377.203125
1	2017-09-23 05:00:00	80145.429688	-49200.355469	-4429.233887	164720.093750	209491.218750	117812.437500	1.485114	40779.996094	194844.875000	235623.390625

30 rows × 11 columns

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

pd.concat([df, forecast_df]).set_index('ds')

	y	unique_id	Add	Add-lo-95	Add-lo-80	Add-hi-80	Add-hi-95	Multi	Multi-lo-95	Multi-lo-80	Multi-hi-80	Multi-hi-95
ds
2017-09-13 00:00:00	80115.0	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2017-09-13 01:00:00	79885.0	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2017-09-13 02:00:00	89325.0	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...
2017-09-23 03:00:00	NaN	NaN	80154.835938	-44792.613281	-1543.909790	161853.578125	205102.281250	115310.718750	3485.836914	42192.359375	188429.078125	227135.593750
2017-09-23 04:00:00	NaN	NaN	80150.132812	-47015.281250	-2998.863037	163299.125000	207315.546875	116561.578125	1745.947998	41487.671875	191635.484375	231377.203125
2017-09-23 05:00:00	NaN	NaN	80145.429688	-49200.355469	-4429.233887	164720.093750	209491.218750	117812.437500	1.485114	40779.996094	194844.875000	235623.390625

246 rows × 12 columns

df_plot= pd.concat([df, forecast_df]).set_index('ds')
df_plot

	y	unique_id	Add	Add-lo-95	Add-lo-80	Add-hi-80	Add-hi-95	Multi	Multi-lo-95	Multi-lo-80	Multi-hi-80	Multi-hi-95
ds
2017-09-13 00:00:00	80115.0	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2017-09-13 01:00:00	79885.0	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2017-09-13 02:00:00	89325.0	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...
2017-09-23 03:00:00	NaN	NaN	80154.835938	-44792.613281	-1543.909790	161853.578125	205102.281250	115310.718750	3485.836914	42192.359375	188429.078125	227135.593750
2017-09-23 04:00:00	NaN	NaN	80150.132812	-47015.281250	-2998.863037	163299.125000	207315.546875	116561.578125	1745.947998	41487.671875	191635.484375	231377.203125
2017-09-23 05:00:00	NaN	NaN	80145.429688	-49200.355469	-4429.233887	164720.093750	209491.218750	117812.437500	1.485114	40779.996094	194844.875000	235623.390625

246 rows × 12 columns

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

def plot_forecasts(y_hist, y_pred, models):
    _, ax = plt.subplots(1, 1, figsize = (20, 7))
    df_plot = pd.concat([y_hist, y_pred]).set_index('ds').tail(12*10)
    df_plot[['y'] + models].plot(ax=ax, linewidth=3 , )
    colors = ['green', "lime"]
    ax.fill_between(df_plot.index, 
                df_plot['Add-lo-80'], 
                df_plot['Add-hi-80'],
                alpha=.20,
                color='fuchsia',
                label='Holt_level_80')
    ax.fill_between(df_plot.index, 
                df_plot['Add-lo-95'], 
                df_plot['Add-hi-95'],
                alpha=.3,
                color='lime',
                label='Holt_level_95')
    ax.set_title('', fontsize=22)
    ax.set_ylabel("Ads watched (hourly data)", fontsize=20)
    ax.set_xlabel('Hourly', fontsize=20)
    ax.legend(prop={'size': 20})
    ax.grid(True)
    plt.show()

plot_forecasts(df, forecast_df, models=["Add"])

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

sf.plot(df, forecast_df)

交叉验证

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

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

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

执行时间序列交叉验证

时间序列模型的交叉验证被视为最佳实践，但大多数实现非常慢。statsforecast库将交叉验证实现为分布式操作，从而使得这一过程降低耗时。如果您有大型数据集，也可以使用Ray、Dask或Spark在分布式集群中进行交叉验证。

在这种情况下，我们希望评估每个模型在过去5个月的表现(n_windows=)，并每隔两个月进行一次预测(step_size=12)。根据您的计算机，这一步骤大约需要1分钟。

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

df: 训练数据框
h (int): 表示向未来预测的步数。在这种情况下，提前30小时。
step_size (int): 每个窗口之间的步长。换句话说：您希望多长时间运行一次预测过程。
n_windows(int): 用于交叉验证的窗口数量。换句话说：您希望评估过去多少个预测过程。

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

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

unique_id: 索引。如果您不喜欢使用索引，只需运行 crossvalidation_df.resetindex()。
ds: 日期戳或时间索引
cutoff: n_windows的最后日期戳或时间索引。
y: 实际值
model: 包含模型名称和拟合值的列。

crossvalidation_df

	ds	cutoff	y	Add	Multi
unique_id
1	2017-09-18 06:00:00	2017-09-18 05:00:00	99440.0	111573.328125	112901.664062
1	2017-09-18 07:00:00	2017-09-18 05:00:00	97655.0	111820.390625	114476.921875
1	2017-09-18 08:00:00	2017-09-18 05:00:00	97655.0	112067.453125	116052.179688
...	...	...	...	...	...
1	2017-09-21 21:00:00	2017-09-20 17:00:00	103080.0	148230.671875	183085.953125
1	2017-09-21 22:00:00	2017-09-20 17:00:00	95155.0	148541.937500	184642.046875
1	2017-09-21 23:00:00	2017-09-20 17:00:00	80285.0	148853.203125	186198.140625

90 rows × 5 columns

模型评估

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

%%capture
!pip install datasetsforecast

计算 RMSE 的函数需要两个参数：

实际值。
预测值，在这种情况下是 Holt 模型。

from datasetsforecast.losses import rmse

def evaluate_cross_validation(df, metric):
    models = df.drop(columns=['ds', 'cutoff', 'y']).columns.tolist()
    evals = []
    for model in models:
        eval_ = df.groupby(['unique_id', 'cutoff']).apply(lambda x: metric(x['y'].values, x[model].values)).to_frame() # 计算每个唯一标识符、模型和截止时间的损失。
        eval_.columns = [model]
        evals.append(eval_)
    evals = pd.concat(evals, axis=1)
    evals = evals.groupby(['unique_id']).mean(numeric_only=True) # 对每种模型和unique_id组合的所有截止值的误差指标取平均值。
    evals['best_model'] = evals.idxmin(axis=1)
    return evals

evaluation_df = evaluate_cross_validation(crossvalidation_df, rmse)
evaluation_df

	Add	Multi	best_model
unique_id
1	31882.099609	46589.265625	Add

致谢

我们感谢Naren Castellon撰写本教程。

参考文献

Changquan Huang • Alla Petukhina. Springer系列（2022）。使用Python的应用时间序列分析与预测。
Ivan Svetunkov. 使用增强动态自适应模型（ADAM）进行预测和分析。
James D. Hamilton. 时间序列分析，普林斯顿大学出版社，新泽西州普林斯顿，第一版，1994年。
Nixtla参数。
Pandas可用频率。
Rob J. Hyndman 和 George Athanasopoulos（2018）。 “预测原则与实践，时间序列交叉验证”。。
季节性周期 - Rob J Hyndman。

Give us a ⭐ on Github

目录

介绍

霍尔特方法

创新状态空间模型用于指数平滑

ETS(A,A,N): 霍尔特线性方法与加性误差

ETS(M,A,N): 霍尔特线性方法与乘性误差

指数平滑方法的分类

ETS 分类法中的数学模型

霍尔特线性趋势方法的属性

加载库和数据

读取数据

使用 plot 方法探索数据

增强型迪基-福勒检验

自相关图

时间序列的分解

加法型时间序列

乘法型时间序列

加法性

乘法性

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

Holt 方法的实现与 StatsForecast

加载库

实例化模型

适配模型

预测方法

使用置信区间的预测方法

交叉验证

执行时间序列交叉验证

模型评估

致谢

参考文献