​
目录

​
引言

自回归 (autoregressive) 时间序列模型 (AutoRegressive) 是一种用于分析和预测单变量时间序列的 统计 技术。本质上，自回归模型 基于时间序列的先前值可用于预测未来值的思想。

在此模型中，因变量（时间序列）在不同时间点回归自身，从而在过去值和当前值之间建立依赖关系。其思想是过去值可以帮助我们理解和预测序列的未来值。

自回归模型 可以拟合不同的阶数，这表明使用了多少个过去值来预测当前值。例如，1 阶 自回归模型 (AR(1))(AR(1)) 仅使用紧邻的前一个值来预测当前值，而 p(AR(p))p (AR(p)) 阶 自回归模型 则使用 pp 个先前值。

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

在 多元回归模型 中，我们使用预测变量的线性组合来预测感兴趣的变量。在 自回归模型 中，我们使用变量过去值的 线性组合 来预测感兴趣的变量。术语 自回归 表明它是变量对自身的回归。

​
自回归模型的定义

在给出 ARCH 模型的正式定义之前，我们先概括地定义 ARCH 模型的组成部分

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

AR 模型是单变量时间序列最基本的构建模块。正如您之前所见，单变量时间序列是一系列模型，它们仅使用目标变量的过去信息来预测其未来，而不依赖于其他解释变量。

定义 1. (1) 以下方程称为 pp 阶自回归模型，记为 AR(p)\text{AR(p)}

Xt=φ0+φ1Xt−1+φ2Xt−2+⋯+φpXt−p+εt \begin{equation} Xt =\varphi_0 +\varphi_1X_{t−1}+\varphi_2X_{t−2}+\cdots+\varphi_p X_{t−p}+\varepsilon_t \tag 1 \end{equation}

其中 {εt}∼WN(0,σϵ2)\{\varepsilon_t \} \sim WN(0,\sigma_{\epsilon}^2)，E(Xsεt)=0E(X_s \varepsilon_t) = 0（当 s<ts < t 时）且 φ0,φ1,⋯ ,φp\varphi_0,\varphi_1,\cdots ,\varphi_p 是实值参数（系数），其中 φp≠0\varphi_p \neq 0。

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

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

  • 为简单起见，我们通常假设截距（常数项） φ0=0\varphi_0 = 0；否则，我们可以考虑 {Xt−μ}\{X_t −\mu \}，其中 μ=φ0/(1−φ1−⋯−φp)\mu =\varphi_0 /(1−\varphi_1 − \cdots −\varphi_p)。

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

  • E(Xsεt)=0(s<t)E(X_s \varepsilon_t) = 0(s < t) 表示过去时刻的 XsX_s 与当前时刻 tt 的 εt\varepsilon_t 无关。

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

此外，使用后移运算符 (见) BB，AR(p)\text{AR(p)} 模型可以改写为

φ(B)Xt=εt\varphi(B)X_t = \varepsilon_t

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

​
PACF 的定义

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

X^t=α1Xt−1+⋯+αk−1Xt−k+1\hat X_t =\alpha_1 X_{t−1}+ \cdots +\alpha_{k−1} X_{t−k+1}

其中 {α1,⋯ ,αk−1}\{\alpha_1, \cdots , \alpha_{k−1} \} 满足

{α1,⋯ ,αk−1}=arg min⁡β1,⋅⋅⋅,βk−1E[Xt−(β1Xt−1+⋯+βk−1Xt−k+1)]2\{\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

也就是说，选择 {α1,⋯ ,αk−1}\{\alpha_1, \cdots , \alpha_{k−1} \} 以最小化预测的均方误差。类似地，设 X^t−k\hat X_{t −k} 表示对 Xt−kX_{t −k} 在 {Xt−k+1:t−1}\{X_{t −k+1:t −1} \} 上的回归（预测）。

X^t−k=η1Xt−1+⋯+ηk−1Xt−k+1\hat X_{t−k} =\eta_1 X_{t−1}+ \cdots +\eta_{k−1} X_{t−k+1}

请注意，如果 {Xt}\{X_t \} 是平稳的，则 {α1:k−1}={η1:k−1}\{ \alpha_{1:k−1}\} = \{\eta_{1:k−1} \}。现在设 Z^t−k=Xt−k−X^t−k\hat Z_{t−k} = X_{t−k} − \hat X_{t−k} 和 Z^t=Xt−X^t\hat Z_t = X_t − \hat X_t。那么 Z^t−k\hat Z_{t−k} 是移除中间变量 {Xt−k+1:t−1}\{X_{t−k+1:t−1} \} 对 Xt−kX_{t−k} 影响后的残差，而 Z^t\hat Z_t 是移除 {Xt−k+1:t−1}\{X_{t −k+1:t −1} \} 对 XtX_t 影响后的残差。

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

ϕ11=Corr(Xt−1,Xt)=Cov(Xt−1,Xt)[Var(Xt−1)Var(Xt)]1/2=ρ1\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

以及

ϕkk=Corr(Z^t−k,Z^t)=Cov(Z^t−k,Z^t)[Var(Z^t−k)Var(Z^t)]1/2\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 年出版物第 172 页），∣ϕkk∣≤1|\phi_{kk}| \leq 1。另一方面，以下定理为估计平稳时间序列的偏自相关函数（PACF）提供了途径，其证明可在 Fan 和 Yao (2003) 中找到。

另一方面，以下定理为估计平稳时间序列的偏自相关函数（PACF）提供了途径，其证明可在 Fan 和 Yao (2003) 中找到。

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

{a1k,⋯ ,akk}=arg min⁡a1,⋯ ,akE(Xt−a1Xt−1−⋯−akXt−k)2\{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

则 ϕkk=akk\phi_{kk}=a_{kk} 对于 k≥1k≥1 成立。

​
自回归模型的性质

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

{X(t−p):(t−1)}={x(t−p):(t−1)}\{X_{(t−p):(t−1)} \} = \{x_{(t−p):(t−1)} \}

则有 E(Xt∣X(t−p):(t−1))=φ0+φ1xt−1+φ2xt−2+⋯+φpxt−pE(X_t |X_{(t−p):(t−1)}) = \varphi_0 + \varphi_1x_{t−1} + \varphi_2 x_{t−2} + \cdots + \varphi_p x_{t−p}

这表明，给定过去的值，该等式的右侧是对 XtX_t 的一个很好的估计。此外

Var(Xt∣X(t−p):(t−1))=Var(εt)=σε2Var(X_t |X_{(t −p):(t −1)}) = Var(\varepsilon_t ) = \sigma_{\varepsilon}^2

现在假设 AR(p)\text{AR(p)} 模型（即，公式 (1)）是平稳的；则我们有

  1. 模型均值 E(Xt)=μ=φ0/(1−φ1−⋅⋅⋅−φp)E(_Xt)=\mu =\varphi_0 / (1−\varphi_1−···−\varphi_p)。因此，模型均值 μ=0\mu=0 当且仅当 φ0=0\varphi_0 =0。

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

Var(Xt)=γ0=φ1γ1+φ2γ2+⋯+φpγp+σε2\text {Var} (X_t) = \gamma_0 = \varphi_1 \gamma_1 + \varphi_2 \gamma2 + \cdots + \varphi_p \gamma_p + \sigma_{\varepsilon}^2

此外

γ0=σε2/(1−φ1ρ1−φ2ρ2−⋯−φpρp)\gamma_0 = \sigma_{\varepsilon}^2 / ( 1 − \varphi_1 \rho_1 − \varphi_2 \rho_2 − \cdots − \varphi_p \rho_p )

  1. 对于所有k>pk > p，偏自相关函数ϕkk=0\phi_{kk} = 0，也就是说AR(p)模型的偏自相关函数在滞后阶数pp之后截尾，这对于识别AR\text{AR}模型非常有帮助。事实上，此时，对XtX_t基于{Xt−k+1:t−1}\{X_{t−k+1:t−1} \}的预测或回归方程是

X^t=φ1Xt−1+⋯+φk−1Xt−k+1\hat X_t =\varphi_1 X_{t−1}+\cdots +\varphi_{k−1} X_{t−k+1}

因此，Xt−X^t=εtX_t − \hat X_t = \varepsilon_t。此外，Xt−k−X^t−kX_{t−k} − \hat X_{t−k}是{Xt−k:t−1}\{ X_{t−k:t−1} \}的函数，并且εt\varepsilon_t与{Xt−k:t−1}\{X_{t−k:t−1} \}中的所有元素不相关。因此

Cov(Xt−k−X^t−k,Xt−X^t)=Cov(Xt−k−X^t−k,εt)=0.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，ϕkk=0\phi_{kk} = 0。

  1. 我们将等式(1)乘以Xt−kX_{t−k}，取期望，然后除以γ0\gamma_0，得到自相关系数之间的递推关系
for k≥1,ρk=φ1ρk−1+φ2ρk−2+⋯+φpρk−p \begin{equation} for \ 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,⋅⋅⋅,pk = 1,2,··· ,p。然后我们得到一组差分方程，这就是著名的Yule-Walker方程。如果已知ACF{ρ1:p}\text{ACF} \{\rho_{1:p} \}，则可以求解Yule-Walker方程得到{φ1:p}\{\varphi_{1:p} \}的估计值，这些解称为Yule-Walker估计。

  1. 由于现在模型是一个平稳的AR(p)\text{AR(p)}模型，它自然满足Xt=φ1Xt−1+φ2Xt−2+⋯+φpXt−p+εtX_t =\varphi_1 X_{t−1}+ \varphi_2 X_{t−2} + \cdots + \varphi_p X_{t−p} + \varepsilon_t。因此ϕpp=φp\phi_{pp} = \varphi_p。如果该AR(p)\text{AR(p)}模型进一步是高斯模型且样本大小为T\text{T}已知，则 (a) ϕ^pp→φp\hat \phi_{pp} → \varphi_p当T→∞T → ∞时； (b) 根据Quenouille (1949)的结论，对于k>p,Tϕ^kkk > p, \sqrt{T} \hat \phi_{kk}渐近服从标准正态（高斯）分布N(0,1)\text{N(0,1)}，或者说ϕkk\phi_{kk}渐近分布为N(0, 1/T )\text{N(0, 1/T )}。

​
AR模型的平稳性和因果性

考虑AR(1)模型

Xt=φXt−1+εt,εt∼WN(0,σε2) \begin{equation} X_t = \varphi X_{t − 1} + \varepsilon_t , \varepsilon_t \sim W N( 0 , \sigma_{\varepsilon}^2 ) \tag 3 \end{equation}

对于 ∣φ∣<1|\varphi|<1，令 X1t=∑j=0∞φjεt−jX_{1t} =\sum_{j=0}^{\infty} \varphi^j \varepsilon_{t−j}；对于 ∣φ∣>1|\varphi|>1，令 X2t=−∑j=1∞φ−jεt+jX_{2t} =− \sum_{j=1}^{\infty} \varphi^{-j} \varepsilon_{t+j}。很容易看出， {X1t}\{X_{1t } \} 和 {X2t}\{X_{2t } \} 都是平稳的并满足方程 (3)。也就是说，两者都是方程 (3) 的平稳解。这就产生了一个问题：两者中哪个更可取？显然，{X2t}\{X_{2t } \} 取决于不可观测的 {εt}\{\varepsilon_t \} 的未来值，因此不自然。因此，我们选择 {X1t}\{X_{1t } \} 并放弃 {X2t}\{X_{2t } \}。换句话说，我们要求方程 (3) 中的系数 φ\varphi 的绝对值小于 1。此时，AR(1)\text{AR}(1) 模型被称为因果模型，其因果表达式为 Xt=∑j=0∞φjεt−jX_t = \sum_{j=0}^{\infty} \varphi^j \varepsilon_{t−j}。一般来说，因果性的定义如下。

定义 3 (1) 时间序列 {Xt}\{X_t \} 是因果的，如果存在系数 ψj\psi_j 使得

Xt=∑j=0∞ψjεt−j,  ∑j=0∞∣ψj∣<∞X_t =\sum_{j=0}^{\infty} \psi_j \varepsilon_{t-j}, \ \ \sum_{j=0}^{\infty} |\psi_j |< \infty

其中 ψ0=1,{εt}∼WN(0,σε2)\psi_0 = 1, \{\varepsilon_t \} \sim WN(0, \sigma_{\varepsilon}^2 )。此时，我们称时间序列 {Xt}\{X_t \} 具有 MA(∞)\text{MA}(\infty) 表示。

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

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

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

注意以下几点：* 根据 Brockwell 和 Davis (2016) 第 75 页的存在性和唯一性，由方程 (1) 定义的 AR\text{AR} 模型是平稳的当且仅当其 AR\text{AR} 多项式 φ(z)=1−φ1z−⋯−φpzp≠0\varphi(z)=1−\varphi_1 z− \cdots − \varphi_p z^p \neq 0 对于所有 ∣z∣=1|z|=1，即 AR\text{AR} 多项式的所有根都不在单位圆上。因此，对于由方程 (1) 定义的 AR 模型，其平稳性条件弱于其因果性条件。

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

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

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

  • 根据二次多项式 φ(z)=1−φ1z−φ2z2\varphi(z) = 1 − \varphi_1 z − \varphi_2 z^2 的根和系数之间的关系，可以证明多项式的两个根模都大于 1 当且仅当

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

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

ψj=∑k=1jφk′ψj−k,  j≥1 where  φk′=φk  if  k≤p  and  φk′=0  if  k>p\psi_j=\sum_{k=1}^{j} \varphi '_k \psi_{j-k}, \ \ j \geq 1 \ where \ \ \varphi '_k =\varphi_k \ \ if \ \ k \leq p \ \ and \ \ \varphi '_k =0 \ \ if \ \ k>p

​
自相关：过去影响现在

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

举一个直观的例子，考虑医生诊室外的等候队列。想象一下医生有一个计划，每个病人可以与他交流 20 分钟。如果每个病人都恰好花费 20 分钟，那么这就会进行得很顺利。但如果一个病人花费的时间稍长呢？如果一次就诊的时长对下一次就诊的时长产生影响，那么可能存在自相关。所以，如果医生需要加快一次预约的速度，因为上一次预约花了太长时间，那么这就是过去和现在之间存在相关性的体现。过去的值会影响未来的值。

​
正自相关和负自相关

与“常规”相关类似，自相关可以是正的或负的。正自相关意味着当前较高的值很可能在下一期产生较高的值。这可以在股票交易中观察到：一旦很多人想买一只股票，它的价格就会上涨。这种积极趋势使得人们更想买这只股票，因为它带来了正回报。买入这只股票的人越多，它的价格就越高，想买入的人就越多。

正相关在下跌趋势中也适用。如果今天的股票价值很低，明天的值可能更低，因为人们开始抛售。当许多人抛售时，价值就会下跌，更多的人会想卖出。这也是一个正自相关的例子，因为过去和现在朝同一个方向发展。如果过去很低，现在就低；如果过去很高，现在就高。

如果两个趋势相反，则存在负自相关。医生就诊时长示例就是这种情况。如果一次就诊时间较长，下一次就诊时间就会较短。如果一次就诊花费的时间较少，医生可能会为下一次就诊多花一点时间。

​
平稳性和 ADF 检验

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

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

一个随时间显著增长或下降的时间序列很容易发现。但有时很难判断时间序列是否平稳。这时，增强迪基-福勒 (ADF) 检验就派上用场了。

​
加载库和数据

提示

需要 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()
日期总计
01986-1-019034
11986-2-019596
21986-3-0110558
31986-4-019002
41986-5-019239

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
01986-1-0190341
11986-2-0195961
21986-3-01105581
31986-4-0190021
41986-5-0192391
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 模型中是很有趣的。

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

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

阶数取决于你。理论上说，你可以根据 PACF 图来确定阶数。理论告诉你，你应该选取自相关性变为 0 之前的滞后数。所有其他滞后应为 0。

理论上，你经常会看到很棒的图表，其中第一个峰值非常高，而其余的都等于零。在这种情况下，选择很容易：你正在处理一个非常“纯粹”的 AR(1) 示例。另一种常见情况是自相关开始很高并缓慢下降到零。在这种情况下，你应该使用 PACF 尚未为零的所有延迟。

然而，在实践中，事情并不总是那么简单。记住那句名言：“所有模型都是错误的，但有些是有用的”。要找到完全符合 AR 模型的情况是非常罕见的。总的来说，自回归过程可以帮助解释变量变异的一部分，但不是全部。

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

​
将数据分割成训练集和测试集

让我们将数据分成两组：1. 用于训练 AutoRegressive 模型的数据 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()

​
使用 StatsForecast 实现 AutoRegressive

​
加载库

from statsforecast import StatsForecast
from statsforecast.models import AutoRegressive

​
实例化模型

导入并实例化模型。设置参数有时会很棘手。主教练 Rob Hyndmann 关于季节周期的这篇文章可能很有用。season_length。

方法 1：我们使用整数格式的 lags 参数，也就是说，我们将希望在模型中评估的 lags 以整数形式输入。

season_length = 12 # Monthly data 
horizon = len(test) # number of predictions biasadj=True, include_drift=True,

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

方法 2：我们使用列表格式的 lags 参数，也就是说，我们将希望在模型中评估的 lags 以列表的形式输入，如下所示。

season_length = 12 # Monthly data 
horizon = len(test) # number of predictions

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

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

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

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

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

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

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

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

​
拟合模型

sf.fit(df=train)

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
残差模型
0-11998.537347
1NaN
2NaN
......
309-2718.312961
310-1306.795172
311-2713.284999
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 不存储拟合值，并且在分布式环境中具有高度的可扩展性。

预测方法接受两个参数：预测未来 h（horizon）和 level。

  • h (int): 表示预测未来 h 个步骤。在此例中为未来 12 个月。

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

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

# Prediction
Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat
unique_idds自回归模型
012012-01-0115905.582031
112012-02-0113597.894531
212012-03-0115488.883789
............
912012-10-0114087.901367
1012012-11-0113274.105469
1112012-12-0112498.226562
values=sf.forecast_fitted_values()
values.head()
unique_iddsy自回归模型
011986-01-019034.021032.537109
111986-02-019596.0NaN
211986-03-0110558.0NaN
311986-04-019002.0126172.937500
411986-05-019239.010020.040039

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

sf.forecast(df=train, h=horizon, level=[95])
unique_idds自回归模型自回归模型-低-95自回归模型-高-95
012012-01-0115905.5820312119.58642629691.578125
112012-02-0113597.894531-188.10113527383.890625
212012-03-0115488.8837891702.88806229274.878906
..................
912012-10-0114087.901367-1050.06835929225.871094
1012012-11-0113274.105469-1886.97314528435.183594
1112012-12-0112498.226562-2675.54760727672.001953
# Merge the forecasts with the true values
Y_hat1 = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat1
dsyunique_id自回归模型
02012-01-0113427115905.582031
12012-02-0114447113597.894531
22012-03-0114717115488.883789
...............
92012-10-0113795114087.901367
102012-11-0113352113274.105469
112012-12-0112716112498.226562
sf.plot(train, Y_hat1)

​
带置信区间的预测方法

使用 predict 方法生成预测。

predict 方法接受两个参数：预测未来 h（horizon）和 level。

  • h (int): 表示预测未来 h 个步骤。在此例中为未来 12 个月。

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

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

此步骤应耗时少于 1 秒。

sf.predict(h=horizon)
unique_idds自回归模型
012012-01-0115905.582031
112012-02-0113597.894531
212012-03-0115488.883789
............
912012-10-0114087.901367
1012012-11-0113274.105469
1112012-12-0112498.226562
forecast_df = sf.predict(h=horizon, level=[95]) 
forecast_df
unique_idds自回归模型自回归模型-低-95自回归模型-高-95
012012-01-0115905.5820312119.58642629691.578125
112012-02-0113597.894531-188.10113527383.890625
212012-03-0115488.8837891702.88806229274.878906
..................
912012-10-0114087.901367-1050.06835929225.871094
1012012-11-0113274.105469-1886.97314528435.183594
1112012-12-0112498.226562-2675.54760727672.001953
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=6,
                                         n_windows=5)

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

  • unique_id：序列标识符
  • ds：日期戳或时间索引
  • cutoff：n_windows 的最后一个日期戳或时间索引。
  • y：真实值
  • "model"：包含模型名称和拟合值的列。
crossvalidation_df
unique_idds截止日期y自回归模型
012009-01-012008-12-0119262.024295.837891
112009-02-012008-12-0120658.023993.947266
212009-03-012008-12-0122660.021201.121094
..................
5712011-10-012010-12-0112893.019349.708984
5812011-11-012010-12-0111843.016899.849609
5912011-12-012010-12-0111321.018159.574219

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

from IPython.display import display

crossvalidation_df.rename(columns = {'y' : 'actual'}, inplace = True) # rename actual values 

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

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

​
模型评估

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

from functools import partial

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

evaluate(
    test.merge(Y_hat),
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=season_length), ufl.rmse, ufl.smape],
    train_df=train,
    agg_fn='mean',
)
指标自回归模型
0平均绝对误差（MAE）962.023763
1平均绝对百分比误差（MAPE）0.072733
2平均绝对比例误差（MASE）0.601808
3均方根误差（RMSE）1195.013050
4对称平均绝对百分比误差（SMAPE）0.034858

​
参考

  1. Changquan Huang • Alla Petukhina. Springer series (2022). Applied Time Series Analysis and Forecasting with Python.
  2. Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). “Models for optimising the theta method and their relationship to state space models”. International Journal of Forecasting.
  3. Nixtla 参数.
  4. Pandas 可用频率.
  5. Rob J. Hyndman and George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”
  6. 季节周期 - Rob J Hyndman.