自回归模型


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

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 等。

目录

引言

自回归时间序列模型(AutoRegressive)是一种统计技术，主要用于分析和预测单变量时间序列。从本质上讲，自回归模型基于这样一个理念：时间序列的历史值可以用来预测未来值。

在该模型中，因变量（时间序列）在不同的时间点返回到自身，形成过去值与现在值之间的依赖关系。其核心思想是，以往的值可以帮助我们理解和预测该序列的未来值。

自回归模型可以拟合不同的阶数，这表示用于预测当前值的过去值的数量。例如，阶数为1的自回归模型 \((AR(1))\) 仅使用前一个值来预测当前值，而阶数为 \(p\) 的自回归模型 \((AR(p))\) 则使用 \(p\) 个以前的值。

自回归模型是时间序列分析的基本模型之一，广泛应用于多个领域，从金融和经济到气象学和社会科学。该模型能够捕捉时间序列数据中的非线性依赖性，这使其在预测和长期趋势分析中尤为有用。

在 多元回归模型中，我们使用预测变量的线性组合来预测感兴趣的变量。在自回归模型中，我们使用该变量的历史值的线性组合来预测感兴趣的变量。术语自回归表明这是一个变量对自身的回归分析。

自回归模型的定义

在给出ARCH模型的正式定义之前，我们先一般性地定义ARCH模型的组成部分：

  • 自回归，一个我们已经了解的概念，是使用统计方法构建单变量时间序列模型，这意味着一个变量的当前值受到其在不同时间段过去值的影响。
  • 异方差性意味着模型在不同时间点可以具有不同的幅度或变异性（方差随时间变化）。
  • 条件性，由于波动性不是固定的，这里的参考是我们在模型中引入的常数，以限制异方差性并使其在条件上依赖于变量的先前值或多个值。

AR模型是单变量时间序列的最基本构建块。正如你之前所看到的，单变量时间序列是一个仅使用目标变量的过去信息来预测其未来，并且不依赖于其他解释变量的模型家族。

定义 1. (1) 下面的方程称为阶数为 \(p\) 的自回归模型，记作 \(\text{AR(p)}\)：

\[ \begin{equation} X_t = \varphi_0 + \varphi_1 X_{t-1} + \varphi_2 X_{t-2} + \cdots + \varphi_p X_{t-p} + \varepsilon_t \tag 1 \end{equation} \]

其中 \(\{\varepsilon_t \} \sim WN(0, \sigma_{\epsilon}^2)\)，\(E(X_s \varepsilon_t) = 0\) 当 \(s < t\) 时，且 \(\varphi_0, \varphi_1, \cdots, \varphi_p\) 是实值参数（系数），且 \(\varphi_p \neq 0\)。

  1. 如果时间序列 \(\{X_t \}\) 是平稳的，并满足如 (1) 所示的方程，我们称其为 \(\text{AR(p)}\) 过程。

注意以下关于此定义的说明：

  • 为了简化，我们通常假设截距（常数项） \(\varphi_0 = 0\)；否则，我们可以考虑 \(\{X_t - \mu \}\)，其中 \(\mu = \varphi_0 /(1−\varphi_1 − \cdots −\varphi_p)\)。

  • 我们将 \(\text{AR}\) 模型的概念与 \(\text{AR}\) 过程的概念区分开。 \(\text{AR}\) 模型可能是平稳的，也可能不是，而 \(\text{AR}\) 过程必须是平稳的。

  • \(E(X_s \varepsilon_t) = 0 (s < t)\) 意味着过去的 \(X_s\) 与当前时间 \(t\) 的 \(\varepsilon_t\) 无关。

  • 类似于MA模型的定义，有时方程 (1) 中的 \(\varepsilon_t\) 被称为创新或冲击项。

此外，使用反向滞后（见）运算符 \(B\)，\(\text{AR(p)}\) 模型可以重写为

\[\varphi(B)X_t = \varepsilon_t\]

其中 \(\varphi(z) = 1 − \varphi_1 z − \cdots − \varphi_p z^p\) 被称为（相应的）\(\text{AR}\) 多项式。此外，在Python包 |StatsModels| 中，\(\varphi(B)\) 被称为 \(\text{AR}\) 滞后多项式。

PACF的定义

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

\[\hat X_t =\alpha_1 X_{t−1}+ \cdots +\alpha_{k−1} X_{t−k+1}\]

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

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

也就是说，\(\{\alpha_1, \cdots , \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}+ \cdots +\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}\) 中去除干扰变量 \(\{X_{t−k+1:t−1} \}\) 效应后的残差，而 \(\hat Z_t\) 是从 \(X_t\) 中去除 \(\{X_{t −k+1:t −1} \}\) 效应后的残差。

定义 2. 对于平稳时间序列 \(\{X_t \}\)，其 \(E(X_t ) = 0\) 的部分自相关函数(PACF) 在滞后 \(k\) 的定义为

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

根据相关系数的性质（参见，例如，Casella 和 Berger 2002，P172），|φkk| ≤ 1。另一方面，以下定理为估计平稳时间序列的PACF提供了依据，其证明可见于 Fan 和 Yao（2003）。

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

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

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

自回归模型的性质

从 \(\text{AR(p)}\) 模型，即公式 (1)，我们可以看出它与多个线性回归模型的形式相同。然而，它用自己的过去来解释当前值。给定过去的值

\[\{X_{(t−p):(t−1)} \} = \{x_{(t−p):(t−1)} \}\]

我们有 \[E(X_t |X_{(t−p):(t−1)}) = \varphi_0 + \varphi_1x_{t−1} + \varphi_2 x_{t−2} + \cdots + \varphi_p x_{t−p}\]

这表明，给定过去，这个方程右侧是 \(X_t\) 的一个很好估计。此外

\[Var(X_t |X_{(t −p):(t −1)}) = Var(\varepsilon_t ) = \sigma_{\varepsilon}^2\]

现在我们假设 AR(p) 模型，即公式 (1)，是平稳的；那么我们有

  1. 模型的均值 \(E(X_t)=\mu =\varphi_0 / (1−\varphi_1−···−\varphi_p)\) 。因此，当且仅当 \(\varphi_0 =0\) 时，模型均值 \(\mu=0\)。

  2. 如果均值为零或 \(\varphi_0 = 0\) （下面的 (3) 和 (4) 有相同的假设），注意到 \(E(X_t \varepsilon_t ) = \sigma_{\varepsilon}^2\) ，我们将公式 (1) 乘以 \(X_t\) ，取期望，得到

\[\text {Var} (X_t) = \gamma_0 = \varphi_1 \gamma_1 + \varphi_2 \gamma_2 + \cdots + \varphi_p \gamma_p + \sigma_{\varepsilon}^2\]

此外

\[\gamma_0 = \sigma_{\varepsilon}^2 / ( 1 − \varphi_1 \rho_1 − \varphi_2 \rho_2 − \cdots − \varphi_p \rho_p )\]

  1. 对于所有 \(k > p\)，偏自相关 \(\phi_{kk} = 0\)，即 \(\text{AR(p)}\) 模型的 PACF 在滞后 \(p\) 之后截断，这在识别 \(\text{AR}\) 模型时非常有用。实际上，在这一点上，\(X_t\) 对 \(\{X_{t−k+1:t−1} \}\) 的预测或回归为

\[\hat X_t =\varphi_1 X_{t−1}+\cdots +\varphi_{k−1} X_{t−k+1}\]

因此，\(X_t − \hat X_t = \varepsilon_t\)。此外，\(X_{t−k} − \hat X_{t−k}\) 是 \(\{ X_{t−k:t−1} \}\) 的一个函数，并且 \(\varepsilon_t\) 与 \(\{X_{t−k:t−1} \}\) 中的每个元素不相关。因此

\[Cov(X_{t−k} −\hat X_{t−k},X_t −\hat X_t)=Cov(X_{t−k} −\hat X_{t−k},\varepsilon_t)=0.\]

根据定义 2，\(\phi_{kk} = 0\).

  1. 我们将方程（1）两边乘以 \(X_{t−k}\)，取期望，除以 \(\gamma_0\)，然后得到自相关之间的递归关系：

\[ \begin{equation} 对于 \ k ≥ 1, \rho_k = \varphi_1 \rho_{k−1} + \varphi_2 \rho_{k−2} + \cdots + \varphi_p \rho_{k−p} \tag 2 \end{equation} \]

对于方程（2），令 \(k = 1,2,··· ,p\)。然后我们得到一组差分方程，这被称为尤尔-沃克方程。如果给定 \(\text{ACF} \{\rho_{1:p} \}\)，那么我们可以解尤尔-沃克方程以获得对 \(\{\varphi_{1:p} \}\) 的估计，这些解称为尤尔-沃克估计。

  1. 由于模型现在是平稳的 \(\text{AR(p)}\)，它自然满足 \(X_t = \varphi_1 X_{t−1}+ \varphi_2 X_{t−2} + \cdots + \varphi_p X_{t−p} + \varepsilon_t\)。因此 \(\phi_{pp} = \varphi_p\)。如果 \(\text{AR(p)}\) 模型进一步符合高斯分布，并且给定样本大小为 \(\text{T}\)，那么 (a) \(\hat \phi_{pp} → \varphi_p\) 当 \(T → ∞\)；(b) 根据Quenouille (1949)，对于 \(k > p, \sqrt{T} \hat \phi_{kk}\) 渐近地遵循标准正态（高斯）分布 \(\text{N(0,1)}\)，或者 \(\phi_{kk}\) 渐近分布为 \(\text{N(0, 1/T )}\)。

AR模型的平稳性和因果性

考虑AR(1)模型：

\[ \begin{equation} X_t = \varphi X_{t − 1} + \varepsilon_t , \varepsilon_t \sim W N( 0 , \sigma_{\varepsilon}^2 ) \tag 3 \end{equation} \]

对于\(|\varphi|<1\)，设\(X_{1t} =\sum_{j=0}^{\infty} \varphi^j \varepsilon_{t−j}\)；对于\(|\varphi|>1\)，设\(X_{2t} =− \sum_{j=1}^{\infty} \varphi^{-j} \varepsilon_{t+j}\)。很容易证明\(\{X_{1t } \}\)和\(\{X_{2t } \}\)都是平稳的，并且满足式(3)。也就是说，二者都是方程(3)的平稳解。这引发了一个问题：二者中哪个更可取？显然，\(\{X_{2t } \}\)依赖于不可观测的\(\{\varepsilon_t \}\)的未来值，因此这是不自然的。因此，我们选择\(\{X_{1t } \}\)而放弃\(\{X_{2t } \}\)。换句话说，我们要求方程(3)中的系数\(\varphi\)的绝对值小于1。在这一点上，\(\text{AR}(1)\)模型被称为因果的，其因果表达式为\(X_t = \sum_{j=0}^{\infty} \varphi^j \varepsilon_{t−j}\)。一般而言，因果性的定义如下。

定义 3 (1) 如果存在系数\(\psi_j\)使得

\[X_t =\sum_{j=0}^{\infty} \psi_j \varepsilon_{t-j}, \ \ \sum_{j=0}^{\infty} |\psi_j |< \infty\]

其中\(\psi_0 = 1, \{\varepsilon_t \} \sim WN(0, \sigma_{\varepsilon}^2 )\)。在这一点上，我们说时间序列\(\{X_t \}\)具有\(\text{MA}(\infty)\)表示。

  1. 如果由模型生成的时间序列是因果的，我们就说该模型是因果的。

因果性暗示时间序列\(\{X_t\}\)是由过去到时间t的白噪声（或创新）造成的。此外，时间序列\(\{X_{2t } \}\)是一个平稳但不因果的例子。为了确定一个\(\text{AR}\)模型是否是因果的，类似于\(\text{MA}\)模型的可逆性，我们有以下定理。

定理 2（因果性定理） 由 Eq.(1) 定义的 \(\text{AR}\) 模型是因果的，当且仅当其 \(\text{AR}\) 多项式 \(\varphi(z)=1−\varphi_1 z− \cdots − \varphi_p z^p\) 的根的模长大于 1 或位于复平面单位圆外。

请注意以下几点： * 在 Brockwell 和 Davis（2016）第 75 页有关存在性和唯一性的论述中，由 Eq.(1) 定义的 \(\text{AR}\) 模型是平稳的，当且仅当其 \(\text{AR}\) 多项式 \(\varphi(z)=1−\varphi_1 z− \cdots − \varphi_p z^p \neq 0\) 对于所有 \(|z|=1\) 或所有的 \(\text{AR}\) 多项式的根不位于单位圆上。因此，对于由 Eq. (1) 定义的 AR 模型，其平稳条件弱于因果性条件。

  • 因果时间序列必然是平稳的。因此，满足因果条件的 \(\text{AR}\) 模型自然是平稳的。但平稳的 \(\text{AR}\) 模型不一定是因果的。

  • 如果由 Eq. (1) 生成的时间序列 \(\{X_t \}\) 不是来自遥远的过去，即 \[t \in T = {\cdots ,−n,\cdots ,−1,0,1,\cdots ,n,\cdots}\]

而是从初始值 \(X_0\) 开始，那么它可能是非平稳的，更不用说因果性了。

  • 根据二次多项式 \(\varphi(z) = 1 − \varphi_1 z − \varphi_2 z^2\) 的根与系数之间的关系，可以证明当且仅当以下条件成立时，多项式的两个根的模长都超过 1：

\[\begin{equation} \left\{ \begin{array}{ll} |\varphi_2| < 1, \\ \varphi_2 + \varphi_1 < 1, \\ \varphi_2 - \varphi_1 < 1,\\ \end{array} \right. \end{equation}\]

因此，我们可以方便地使用这三个不等式来决定一个 \(\text{AR(2)}\) 模型是否因果。

  • 可以证明，对于由 Eq. (1) 定义的 \(\text{AR(p)}\) 模型，定义 3 中的系数 \(\{\psi_j \}\) 满足 \(\psi_0=1\)，并且

\[\psi_j=\sum_{k=1}^{j} \varphi '_k \psi_{j-k}, \ \ j \geq 1 \text{ 其中 } \ \ \varphi '_k =\varphi_k \ \text{如果} \ \ k \leq p \ \text{且} \ \ \varphi '_k =0 \ \text{如果} \ \ k>p\]

自相关：过去影响现在

自回归模型描述了变量的现在与其过去之间的关系。因此，它适合于过去和现在值相关的变量。

作为一个直观的例子，考虑一下医生的等待队伍。假设医生的计划是每位患者与他有20分钟的时间。如果每位患者正好花费20分钟，这项计划运作良好。但是，如果一位患者花费的时间稍长呢？如果一个查询的持续时间对下一个查询的持续时间有影响，可能会存在自相关。因此，如果医生需要加快一个预约，因为上一个预约花费了太长时间，那么就要看看过去与现在之间的相关性。过去的值影响未来的值。

正自相关与负自相关

与“常规”相关性一样，自相关可以是正的或负的。正自相关意味着现在的高值可能在下一个时期也会导致高值。这可以在股票交易中观察到：一旦很多人想要购买一只股票，它的价格就会上涨。这种正趋势使得人们更想购买这只股票，因为它有正回报。购买股票的人越多，价格就越高，越多的人会想要购买。

正相关在下行趋势中也同样有效。如果今天的股票价值较低，那么明天的价值很可能会更低，因为人们开始抛售。当很多人抛售时，价值下跌，甚至更多的人会想要抛售。这也是正自相关的一个例子，因为过去和现在的趋势朝着同一个方向发展。如果过去低，现在就低；如果过去高，现在就高。

如果两个趋势相反，则存在负自相关。这在医生就诊时间的例子中就是这种情况。如果一个查询花费的时间较长，那么下一个查询将会更短。如果一次就诊花费的时间较少，医生可能会在下一次花费更多时间。

平稳性和ADF检验

数据中存在趋势的问题在单变量时间序列建模中是普遍存在的。时间序列的平稳性意味着时间序列没有（长期）趋势：它在相同的平均值周围是稳定的。否则，时间序列被称为非平稳的。

理论上，AR模型可以在模型中具有趋势系数，但由于平稳性是一般时间序列理论中的一个重要概念，因此最好立即学习如何处理它。许多模型只能在平稳时间序列上运行。

随着时间的推移，明显持续增长或下降的时间序列容易被发现。但有时很难判断一个时间序列是否是平稳的。这就是增强型迪基-福勒（ADF）检验派上用场的地方。

加载库和数据

Tip

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

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

import pandas as pd

import scipy.stats as stats

import statsmodels.api as sm
import statsmodels.tsa.api as smt
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, 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/catfish.csv")
df.head()
Date Total
0 1986-1-01 9034
1 1986-2-01 9596
2 1986-3-01 10558
3 1986-4-01 9002
4 1986-5-01 9239

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 1986-1-01 9034 1
1 1986-2-01 9596 1
2 1986-3-01 10558 1
3 1986-4-01 9002 1
4 1986-5-01 9239 1 
print(df.dtypes)
ds           object
y             int64
unique_id    object
dtype: object

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

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

使用plot方法探索数据

使用StatsForecast类的plot方法绘制一些序列。此方法从数据集中随机打印8个序列，适用于基本的探索性数据分析（EDA）。

from statsforecast import StatsForecast

StatsForecast.plot(df)

增强型迪基-福勒检验

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

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

原假设：时间序列是非平稳的。它给出了一个时间依赖的趋势。 备择假设：时间序列是平稳的。换句话说，该序列不依赖于时间。

ADF或t统计量 < 临界值：拒绝原假设，时间序列是平稳的。 ADF或t统计量 > 临界值：未能拒绝原假设，时间序列是非平稳的。

让我们检查一下我们正在分析的序列是否是平稳序列。让我们创建一个函数来进行检查，使用Dickey Fuller检验。

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"],'Sales')
Dickey-Fuller test results for columns: Sales
Test Statistic          -1.589903
p-value                  0.488664
No Lags Used            14.000000
                          ...    
Critical Value (1%)     -3.451691
Critical Value (5%)     -2.870939
Critical Value (10%)    -2.571778
Length: 7, dtype: float64
Conclusion:====>
The null hypothesis cannot be rejected
The data is not stationary

在之前的结果中，我们可以看到“增强型迪基-福勒（Augmented Dickey-Fuller）”检验给出的“p值”为0.488664，这告诉我们不能拒绝原假设，另一方面，我们序列的数据是不平稳的。

我们需要对我们的时间序列进行差分，以便将数据转换为平稳序列。

Augmented_Dickey_Fuller_Test_func(df["y"].diff().dropna(),"Sales")
Dickey-Fuller test results for columns: Sales
Test Statistic          -4.310935
p-value                  0.000425
No Lags Used            17.000000
                          ...    
Critical Value (1%)     -3.451974
Critical Value (5%)     -2.871063
Critical Value (10%)    -2.571844
Length: 7, dtype: float64
Conclusion:====>
Reject the null hypothesis
The data is stationary

通过应用微分，我们的时间序列现在是平稳的。

def tsplot(y, lags=None, figsize=(12, 7), style='bmh'): # [3]
    if not isinstance(y, pd.Series):
        y = pd.Series(y)
        
    with plt.style.context(style):    
        fig = plt.figure(figsize=figsize)
        layout = (2, 2)
        ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
        acf_ax = plt.subplot2grid(layout, (1, 0))
        pacf_ax = plt.subplot2grid(layout, (1, 1))
        
        y.plot(ax=ts_ax)
        p_value = sm.tsa.stattools.adfuller(y)[1]
        ts_ax.set_title('Time Series Analysis plot\n Dickey-Fuller: p={0:.5f}'.format(p_value))
        smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
        smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
        plt.tight_layout()
tsplot(df["y"].diff().dropna(), lags=20);

正如您所见，从图表中蓝色背景阴影区域可以看出，PACF显示了第一、第二、第三、第四、第六、第七、第九和第十等延迟超出了阴影区域。这意味着将这些滞后值纳入AR模型也会很有趣。

我们应该包含多少个滞后？

现在，时间序列分析中的一个大问题始终是包含多少个滞后。这被称为时间序列的阶数。阶数的符号为AR(1)表示阶数为1，AR(p)表示阶数为p。

阶数由您决定。从理论上讲，您可以根据PACF图表确定阶数。理论告诉您在自相关为0之前应取滞后数。其他所有滞后应为0。

在理论上，您经常会看到很好的图表，其中第一个峰值非常高，其他全部为零。在这种情况下，选择很简单：您正在处理一个非常“纯粹”的AR(1)示例。另一个常见情况是当您的自相关值较高并缓慢降低到零。在这种情况下，您应该使用PACF尚未为零的所有滞后。

然而，在实践中，这并不总是那么简单。请记住那句著名的格言“所有模型都是错误的，但有些是有用的”。很少有案例完美符合AR模型。一般而言，自回归过程可以帮助解释变量的一部分变异，但不是全部。

在实践中，您将尝试选择能够为您的模型提供最佳预测性能的滞后数。最佳预测性能通常不是通过查看自相关图来定义的：那些图给您提供的是理论估计。然而，预测性能最佳的定义是通过模型评估和基准测试来实现的，使用您在模块2中见过的技术。在本模块的后面，我们将看到如何使用模型评估为AR模型选择性能阶数。但在我们深入探讨之前，是时候仔细研究AR模型的确切定义了。

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

让我们将数据分为多个集合 1. 用于训练我们的 自回归 模型的数据 2. 用于测试我们模型的数据

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

train = df[df.ds<='2011-12-01'] 
test = df[df.ds>'2011-12-01']
train.shape, test.shape
((312, 3), (12, 3))

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

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

自动回归模型的实现

要了解 自动回归模型 函数的参数，下面列出了相关信息。更多详情请访问 文档

lags : int 或 list
    包含在模型中的滞后数量。
    如果传入一个整数，则考虑所有滞后至 `lags`。
    如果是列表，则仅考虑列表中的元素作为滞后。
include_mean : bool (默认为 True)
    自动回归模型是否应包含均值项？
    对于未差分的序列，默认值为 True，对于差分序列默认为 False （均值不会影响拟合或预测）。
include_drift : bool (默认为 False)
    自动回归模型是否应包含线性漂移项？
    （即，拟合一个带有自动回归误差的线性回归。）
blambda : float，可选（默认为 None）
    Box-Cox 变换参数。
biasadj : bool (默认为 False)
    使用调整后的反变换均值 Box-Cox。
method : str (默认为 'CSS-ML')
    拟合方法：最大似然或最小条件平方和。
    默认（除非有缺失值）是使用条件平方和来寻找起始值，然后进行最大似然估计。
fixed : dict，可选（默认为 None）
    包含自动回归模型的固定系数的字典。示例： `{'ar1': 0.5, 'ar5': 0.75}`。
    对于自回归项使用 `ar{i}` 键。
alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    用于计算符合预测区间的信息。
    默认情况下，模型将计算原生预测区间。

加载库

from statsforecast import StatsForecast
from statsforecast.models import AutoRegressive

实例化模型

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

方法 1： 我们使用整数格式的滞后参数，也就是说，我们将要评估的滞后放入模型中。

season_length = 12 # 月度数据 
horizon = len(test) # 预测数量偏差调整=True，包含漂移=True，

models2 = [AutoRegressive(lags=[14], include_mean=True)]

方法 2: 我们以列表格式使用 lags 参数，也就是说，我们将想要在模型中评估的滞后值以列表的形式放置，如下所示。

season_length = 12 # 月度数据 
horizon = len(test) # 预测的数量

models = [AutoRegressive(lags=[3,4,6,7,9,10,11,12,13,14], include_mean=True)]

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

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

  • freq: 字符串，表示数据的频率。 （请参见 pandas 的可用频率。）

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

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

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

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

拟合模型

sf.fit()
StatsForecast(models=[AutoRegressive])

让我们看看我们的Theta模型的结果。我们可以通过以下指令进行观察：

result=sf.fitted_[0,0].model_
print(result.keys())
dict_keys(['coef', 'sigma2', 'var_coef', 'mask', 'loglik', 'aic', 'arma', 'residuals', 'code', 'n_cond', 'nobs', 'model', 'aicc', 'bic', 'xreg', 'lambda', 'x'])

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

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

residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual
residual Model
0 -11990.435671
1 NaN
2 NaN
... ...
309 -2717.643794
310 -1306.165240
311 -2712.663024

312 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 AutoRegressive
unique_id
1 2012-01-01 15904.990234
1 2012-02-01 13597.334961
1 2012-03-01 15488.281250
... ... ...
1 2012-10-01 14086.938477
1 2012-11-01 13273.106445
1 2012-12-01 12497.278320

12 rows × 2 columns

values=sf.forecast_fitted_values()
values.head()
ds y AutoRegressive
unique_id
1 1986-01-01 9034.0 21024.435547
1 1986-02-01 9596.0 NaN
1 1986-03-01 10558.0 NaN
1 1986-04-01 9002.0 129223.968750
1 1986-05-01 9239.0 10019.207031

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

sf.forecast(h=horizon, level=[95])
ds AutoRegressive AutoRegressive-lo-95 AutoRegressive-hi-95
unique_id
1 2012-01-01 15904.990234 1759.574219 30050.406250
1 2012-02-01 13597.334961 -548.080872 27742.751953
1 2012-03-01 15488.281250 1342.864990 29633.697266
... ... ... ... ...
1 2012-10-01 14086.938477 -1445.719727 29619.597656
1 2012-11-01 13273.106445 -2283.256104 28829.470703
1 2012-12-01 12497.278320 -3072.115234 28066.671875

12 rows × 4 columns

Y_hat=Y_hat.reset_index()
Y_hat
unique_id ds AutoRegressive
0 1 2012-01-01 15904.990234
1 1 2012-02-01 13597.334961
2 1 2012-03-01 15488.281250
... ... ... ...
9 1 2012-10-01 14086.938477
10 1 2012-11-01 13273.106445
11 1 2012-12-01 12497.278320

12 rows × 3 columns

# 将预测结果与真实值合并
test['unique_id'] = test['unique_id'].astype(int)
Y_hat1 = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat1
ds y unique_id AutoRegressive
0 2012-01-01 13427 1 15904.990234
1 2012-02-01 14447 1 13597.334961
2 2012-03-01 14717 1 15488.281250
... ... ... ... ...
9 2012-10-01 13795 1 14086.938477
10 2012-11-01 13352 1 13273.106445
11 2012-12-01 12716 1 12497.278320

12 rows × 4 columns

# 将预测结果与真实值合并

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat1]).set_index('ds')
plot_df[['y', "AutoRegressive"]].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)
plt.show()

带置信区间的预测方法

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

预测方法接受两个参数：预测下一个 h（代表时间跨度）和 level。

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

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

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

此步骤应少于 1 秒。

sf.predict(h=horizon)
ds AutoRegressive
unique_id
1 2012-01-01 15904.990234
1 2012-02-01 13597.334961
1 2012-03-01 15488.281250
... ... ...
1 2012-10-01 14086.938477
1 2012-11-01 13273.106445
1 2012-12-01 12497.278320

12 rows × 2 columns

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

forecast_df
ds AutoRegressive AutoRegressive-lo-95 AutoRegressive-hi-95
unique_id
1 2012-01-01 15904.990234 1759.574219 30050.406250
1 2012-02-01 13597.334961 -548.080872 27742.751953
1 2012-03-01 15488.281250 1342.864990 29633.697266
... ... ... ... ...
1 2012-10-01 14086.938477 -1445.719727 29619.597656
1 2012-11-01 13273.106445 -2283.256104 28829.470703
1 2012-12-01 12497.278320 -3072.115234 28066.671875

12 rows × 4 columns

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

pd.concat([df, forecast_df]).set_index('ds')
y unique_id AutoRegressive AutoRegressive-lo-95 AutoRegressive-hi-95
ds
1986-01-01 9034.0 1 NaN NaN NaN
1986-02-01 9596.0 1 NaN NaN NaN
1986-03-01 10558.0 1 NaN NaN NaN
... ... ... ... ... ...
2012-10-01 NaN NaN 14086.938477 -1445.719727 29619.597656
2012-11-01 NaN NaN 13273.106445 -2283.256104 28829.470703
2012-12-01 NaN NaN 12497.278320 -3072.115234 28066.671875

336 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)
    plt.show()
plot_forecasts(train, test, forecast_df, models=["AutoRegressive"])

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

交叉验证

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

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

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

执行时间序列交叉验证

时间序列模型的交叉验证被视为最佳实践，但大多数实现速度非常慢。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=6,
                                         n_windows=5)

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

  • unique_id: 索引。如果您不喜欢使用索引，只需运行crossvalidation_df.resetindex()
  • ds: 日期戳或时间索引
  • cutoff: n_windows的最后日期戳或时间索引。
  • y: 实际值
  • "model": 包含模型名称和拟合值的列。
crossvalidation_df
ds cutoff y AutoRegressive
unique_id
1 2009-01-01 2008-12-01 19262.0 24295.751953
1 2009-02-01 2008-12-01 20658.0 23993.804688
1 2009-03-01 2008-12-01 22660.0 21200.978516
... ... ... ... ...
1 2011-10-01 2010-12-01 12893.0 19350.546875
1 2011-11-01 2010-12-01 11843.0 16900.646484
1 2011-12-01 2010-12-01 11321.0 18160.392578

60 rows × 4 columns

现在我们将绘制每个截止期的预测。为了使图表更清晰，我们将重新命名每个期间的实际值。

crossvalidation_df.rename(columns = {'y' : 'actual'}, inplace = True) # 重命名实际值 

cutoff = crossvalidation_df['cutoff'].unique()

for k in range(len(cutoff)): 
    cv = crossvalidation_df[crossvalidation_df['cutoff'] == cutoff[k]]
    StatsForecast.plot(df, cv.loc[:, cv.columns != 'cutoff'])

模型评估

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

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

计算RMSE的函数有两个参数：

  1. 实际值。
  2. 预测值，在这种情况下是自回归。
rmse = rmse(crossvalidation_df['actual'], crossvalidation_df["AutoRegressive"])
print("RMSE using cross-validation: ", rmse)
RMSE using cross-validation:  3566.5479

正如您所注意到的，我们使用交叉验证结果对我们的模型进行了评估。

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

from datasetsforecast.losses import (mae, mape, mase, rmse, smape)
def evaluate_performace(y_hist, y_true, y_pred, model):
    y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
    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_true['y'].values, 
                                                y_true[model].values, 
                                                y_hist['y'].values, seasonality=12)
        else:
            evaluation[model][metric_name] = metric(y_true['y'].values, y_true[model].values)
    return pd.DataFrame(evaluation).T
evaluate_performace(train, test, Y_hat, model="AutoRegressive")
mae mape mase rmse smape
AutoRegressive 961.773193 7.271437 0.601651 1194.660588 6.970027

致谢

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

参考文献

  1. Changquan Huang • Alla Petukhina. Springer 系列 (2022). 使用 Python 进行应用时间序列分析与预测。
  2. Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). “优化 theta 方法的模型及其与状态空间模型的关系”. 国际预测学杂志。
  3. Nixtla 参数。
  4. Pandas 可用频率。
  5. Rob J. Hyndman 和 George Athanasopoulos (2018). “预测原则与实践，时间序列交叉验证”.。
  6. 季节性周期 - Rob J Hyndman。

Give us a ⭐ on Github