平稳时间序列分析与python实现

1. 引言

    前面我们介绍了时间序列的定义以及如何对一个时间序列进行平稳性检验和随机性检验，那么，当一个序列被检测出来是平稳时间序列，并且是非白噪声序列时，我们应该如何对他进行进一步分析，以提取出句子中的规律信息呢？目前，对平稳时间序列进行拟合的模型主要有AR模型、MA模型和ARMA模型，本文将分别对这三者进行具体介绍，并用python来实现它。

2. 平稳时间序列分析

2.1 线性差分方程

2.1.1 线性差分方程的定义

    在介绍这三个模型之前，先介绍一个知识点，即线性差分方程。一般来说，我们称具有如下形式的方程为序列 { z t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{ z _ { t } , t = 0 , \pm 1 , \pm 2 , \cdots \right\} {zt​,t=0,±1,±2,⋯}的线性差分方程：
$$z_t+a_1{z}_{t-1}+a_2{z}_{t-2}+\cdots +a_p{z}_{t-p}=h(t)$$其中，$p⩾1;a_1,a_2,\cdots ,a_p$为实数，$h(t)$为$t$的已知函数。当$h(t)=0$时，则称差分方程为齐次线性差分方程，否则称为非齐次线性差分方程。

2.1.2 齐次线性差分方程的通解

    齐次线性差分方程的求解需要借助特征方程和特征根，其特征方程为：
$${\lambda }^p+a_1{\lambda }^{p-1}+a_2{\lambda }^{p-2}+\cdots +a_p=0$$，特征根为该特征方程对应的根，不妨记作
$${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$$根据特征根的取值情况，齐次线性差分方程的通解会有不同的表达方式，下面分三种情况列举：

• 当${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$无重根时，此时齐次线性差分方程的通解为：
  $$z_t=c_1{\lambda }_1^t+c_2{\lambda }_2^t+\cdots +c_p{\lambda }_p^t$$其中，$c_1,c_2,\cdots ,c_p$为任意的实数。
• 当${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$有重根时，不妨假定${\lambda }_1={\lambda }_2=\cdots ={\lambda }_d$为$d$个相同的实根，${\lambda }_{d+1},{\lambda }_{d+2},\cdots ,{\lambda }_p$为其他不相同的实根，则此时齐次线性差分方程的通解为：
  $$z_t=\left(c_1+c_{2}t+\cdots +c_d{t}^{d-1}\right){\lambda }_1^t+c_{d+1}{\lambda }_{d+1}^t+\cdots +c_p{\lambda }_p^t$$其中，$c_1,c_2,\cdots ,c_p$为任意的实数。
• 当${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$有复根时，不妨设${\lambda }_1=a+ib=r{\mathrm{e}}^{i\omega },{\lambda }_2=a-ib=r{\mathrm{e}}^{-i\omega }$为其中的一对复根，由欧拉公式易得$r=\sqrt{a^2+b^2},\omega =\arccos \frac{a}{r}$，而${\lambda }_3,{\lambda }_4,\cdots ,{\lambda }_p$为其中不相同的实根，则此时齐次线性差分方程的通解为：
  $$\begin{array}{rl}{z}_t& =c_1{\lambda }_1^t+c_2{\lambda }_2^t+\cdots +c_p{\lambda }_p^t\\ & =r^t\left(c_1{\mathrm{e}}^{it\omega }+c_2{\mathrm{e}}^{-it\omega }\right)+c_3{\lambda }_3^t+\cdots +c_p{\lambda }_p^t\end{array}$$其中，$c_1,c_2,\cdots ,c_p$为任意的实数。

2.1.3 非齐次线性差分方程的通解

    非齐次线性差分方程的通解为对应的齐次线性差分方程的通解$z_t^{\prime }$和一个特解$z_t^{\prime \prime }$的加和，即：
$$z_t=z_t^{\prime }+z_t^{\prime \prime }$$其中，特解$z_t^{\prime \prime }$就是任意一个使得非齐次线性差分方程成立的值，即：
$$z_t^{\prime \prime }+a_1{z}_{t-1}^{\prime \prime }+a_2{z}_{t-2}^{\prime \prime }+\cdots +a_p{z}_{t-p}^{\prime \prime }=h(t)$$

2.2 AR模型

2.2.1 AR模型的定义

    AR模型也称为自回归模型，其定义如下：
$$\{\begin{array}{l}{x}_t={\varphi }_0+{\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+\cdots +{\varphi }_p{x}_{t-p}+{\epsilon }_t\\ {\varphi }_p̸=0\\ E\left({\epsilon }_t\right)=0,\mathrm{Var}\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s̸=t\\ E\left(x_s{\epsilon }_t\right)=0,\forall s<t\end{array}$$上述的结构也称为$p$阶自回归模型，简记为$AR(p)$。
    当对序列进行中心化后，并将AR模型用延迟算子表示时，此时AR模型可以表示为：
$$\Phi (B)x_t={\epsilon }_t$$其中，$\Phi (B)=1-{\varphi }_{1}B-{\varphi }_2{B}^2-\cdots -{\varphi }_p{B}^p$，称为$p$阶自回归系数多项式。

2.2.2 AR模型平稳性的判别

    虽然AR模型可以拟合平稳时间序列，但是并不是所有的AR模型都是平稳的，一般来讲，可以通过特征根判别或平稳域判别来判断AR模型的平稳性。

1. 特征根判别
   对于任一中心化$\mathrm{A}\mathrm{R}(p)$模型$\Phi (B)x_t={\epsilon }_t$都可以被视为一个非齐性线性差分方程：
   $$x_t-{\varphi }_1{x}_{t-1}-{\varphi }_2{x}_{t-2}-\cdots -{\varphi }_p{x}_{t-p}={\epsilon }_t$$其通解为：
   $$x_t=x_t^{\prime }+x_t^{\prime \prime }$$其中，$x_t^{\prime }$为齐次线性差分方程$\Phi (B)x_t=0$的通解，$x_t^{\prime \prime }$为特解。
   对于齐次线性差分方程$\Phi (B)x_t=0$的通解，记其特征根为${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$，为了不失一般性，我们假设有$d$个相同的实根，$p-d-2m$个不相等的实根，$m$对复根，则由前面介绍的线性差分方程我们可以得到通解的计算公式如下：
   $$x_t^{\prime }=\sum _{j=1}^d{c}_j{t}^{j-1}{\lambda }_1^t+\sum _{j=d+1}^{p-2m}{c}_j{\lambda }_j^t+\sum _{j=1}^m{r}_j^t\left(c_{1j}{\mathrm{e}}^{it{w}_j}+c_{2j}{\mathrm{e}}^{-it{w}_j}\right)$$其中，$c_1,\cdots ,c_{p-2m},c_{1j},c_{2j}(j=1,\cdots ,m)$为任意实数。
   对于特解，首先，可以证明自回归系数多项式方程$\Phi (u)=0$的根是齐次线性差分方程$\Phi (B)x_t=0$的特征根的倒数。因此，$\Phi (B)$可以因子分解为：
   $$\Phi (B)=\prod _{i=1}^p\left(1-{\lambda }_{i}B\right)$$这样一来，特解$x_t^{\prime \prime }$的计算公式可以表示为：
   $$x_t^{\prime \prime }=\frac{{\epsilon }_t}{\Phi (B)}=\frac{{\epsilon }_t}{{\prod }_{i=1}^p\left(1-{\lambda }_{i}B\right)}=\sum _{i=1}^p\frac{k_i}{1-{\lambda }_{i}B}{\epsilon }_t$$其中，$k_i(i=1,2,\cdots ,p)$为常数。
   因此，最终$\mathrm{A}\mathrm{R}(p)$模型的通解的计算公式为：
   $$\begin{array}{rl}{x}_t& =x_t^{\prime }+x_t^{\prime \prime }\\ & =\sum _{j=1}^d{c}_j{t}^{j-1}{\lambda }_1^t+\sum _{j=d+1}^{p-2m}{c}_j{\lambda }_j^t+\sum _{j=1}^m{r}_j^t\left(c_{1j}{\mathrm{e}}^{it{w}_j}+c_{2j}{\mathrm{e}}^{-it{w}_j}\right)+\sum _{i=1}^p\frac{k_i}{1-{\lambda }_{i}B}{\epsilon }_t\end{array}$$
   要使得$\mathrm{A}\mathrm{R}(p)$模型平稳，即要求对于任意的实数$c_1,\cdots ,c_{p-2m},c_{1j},c_{2j}(j=1,2,\cdots ,m)$有：
   $$\underset{t\to \infty }{\lim }{x}_t=0$$要使得上式成立，其充要条件为：
   $$\begin{array}{l}\left|{\lambda }_i\right|<1,i=1,2,\cdots ,p-2m\\ \left|r_i\right|<1,i=1,2,\cdots ,m\end{array}$$即$\mathrm{A}\mathrm{R}(p)$的$p$个特征根都在单位圆内，或者自回归系数多项式的根都在单位圆外。
2. 平稳域判别
   平稳域判别则直接根据自回归系数进行判断，即对参数向量${\left({\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_p\right)}^{\prime }$的取值范围进行限制，使得特征根都在单位圆内，即：
   { ϕ 1 , ϕ 2 , ⋯ &ThinSpace; , ϕ p ∣ 特 征 根 都 在 单 位 圆 内 } \left\{ \phi _ { 1 } , \phi _ { 2 } , \cdots , \phi _ { p } |\right.特征根都在单位圆内\} {ϕ1​,ϕ2​,⋯,ϕp​∣特征根都在单位圆内}该集合即称为$\mathrm{A}\mathrm{R}(p)$模型的平稳域，一般这种方法对于低阶$AR$模型比较方便。

2.2.3 平稳AR模型的统计性质

1. 均值
   对于$AR(p)$模型，在等式两边取期望得：
   $$E{x}_t=E\left({\varphi }_0+{\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+\cdots +{\varphi }_p{x}_{t-p}+{\epsilon }_t\right)$$，由于平稳序列的均值为常数$E{x}_t=\mu (\forall t\in T)$，并且 { ε t } \left\{ \varepsilon _ { t } \right\} {εt​}为白噪声序列，其均值$E{\epsilon }_t=0$，所以有：
   $$\begin{array}{rl}& \left(1-{\varphi }_1-\dots -{\varphi }_p\right)\mu ={\varphi }_0\\ \Rightarrow \mu & =\frac{{\varphi }_0}{1-{\varphi }_1-\cdots -{\varphi }_p}\end{array}$$
2. 方差
   要计算$AR(p)$模型的方差，需要引入Green函数，其定义如下，设${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$为平稳$AR(p)$模型的特征根，则有：
   $$\begin{array}{rl}{x}_t& =\frac{{\epsilon }_t}{\Phi (B)}\\ & =\sum _{i=1}^p\frac{k_i}{1-{\lambda }_{i}B}{\epsilon }_t\\ & =\sum _{i=1}^{\infty }\sum _{j=0}^{\infty }{k}_i{\left({\lambda }_{i}B\right)}^j{\epsilon }_t\\ & =\sum _{j=0}^{\infty }\sum _{i=1}^p{k}_i{\lambda }_i^j{\epsilon }_{t-j}\\ & \equiv \sum _{j=0}^{\infty }{G}_j{\epsilon }_{t-j}\end{array}$$其中，$G_0=1,G_j={\sum }_{i=1}^p{k}_i{\lambda }_i^j(j=1,2,\cdots )$。
   因此，$AR(p)$模型的方差为：
   $$\mathrm{Var}\left(x_t\right)=\sum _{j=0}^{\infty }{G}_j^2\mathrm{Var}\left({\epsilon }_t\right)$$
3. 协方差函数
   在平稳模型$x_t={\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+\cdots +{\varphi }_p{x}_{t-p}+{\epsilon }_t$两边同乘以$x_{t-k}(\forall k⩾1)$，并求期望得：
   $$E\left(x_t{x}_{t-k}\right)={\varphi }_{1}E\left(x_{t-1}{x}_{t-k}\right)+\cdots +{\varphi }_{p}E\left(x_{t-p}{x}_{t-k}\right)+E\left({\epsilon }_t{x}_{t-k}\right),\forall k⩾1$$由于$E\left({\epsilon }_t{x}_{t-k}\right)=0,\forall k⩾1$，因此，可得自协方差函数的递推公式：
   $${\gamma }_k={\varphi }_1{\gamma }_{k-1}+{\varphi }_2{\gamma }_{k-2}+\cdots +{\varphi }_p{\gamma }_{k-p}$$
4. 自相关系数具有拖尾性
   由自相关系数的计算公式${\rho }_k=\frac{{\gamma }_k}{{\gamma }_0}$易得$AR(p)$模型自相关系数的递推公式：
   $${\rho }_k={\varphi }_1{\rho }_{k-1}+{\varphi }_2{\rho }_{k-2}+\cdots +{\varphi }_p{\rho }_{k-p}$$由自相关系数的递推公式可以发现，其实该表达式是一个$p$阶齐次线性差分方程，其滞后任意$k$阶的自相关系数的通解为：
   $${\rho }_k=\sum _{i=1}^p{c}_i{\lambda }_i^k$$式中，$\left|{\lambda }_i\right|<1(i=1,2,\cdots ,p)$为该差分方程的特征根，$c_1,c_2,\cdots ,c_p$为任意实数，如果已知任意$p$个自相关系数，那么就可以求出$c_1,c_2,\cdots ,c_p$对应的值，因此，$c_1,c_2,\cdots ,c_p$不可能全为0。因此，由该式也容易得到${\rho }_k$始终有非0取值，即不会在$k$大于某个值后就恒等于0，这也是$AR$模型自相关系数的一个重要特性，即拖尾性。
   另外，随着$k$的增大，由于$\left|{\lambda }_i\right|<1(i=1,2,\cdots ,p)$，所以当$k\to \infty$时，${\rho }_k={\sum }_{i=1}^p{c}_i{\lambda }_i^k\to 0$，因此，可以发现$AR$模型自相关系数的另一个重要特性，即指数衰减性，这也是我们前面讲的为什么平稳时间序列的自相关系数很快会衰减到0附近。
5. 偏自相关系数具有p步截尾性
   对于一个平稳$AR(p)$模型，前面所计算出的滞后$k$阶的自相关系数${\rho }_k$实际上并不是$x_t$与$x_{t-k}$之间单纯的相关关系，因为$x_t$还会受到中间$k-1$个随机变量$x_{t-1},x_{t-2},\cdots ,x_{t-k+1}$的影响，因此，为了能单纯地测度$x_{t-k}$对$x_t$的影响，引进了偏自相关系数的概念，其定义如下：
   $${\rho }_{x_t,x_{t-k}|x_{t-1},\cdots ,x_{t-k+1}}=\frac{E\left[\left(x_t-\hat{E}{x}_t\right)\left(x_{t-k}-\hat{E}{x}_{t-k}\right)\right]}{E\left[{\left(x_{t-k}-\hat{E}{x}_{t-k}\right)}^2\right]}$$
   其中，$\hat{E}{x}_t=E\left[x_t|x_{t-1},\cdots ,x_{t-k+1}\right]$，$\hat{E}{x}_{t-k}=E\left[x_{t-k}|x_{t-1},\cdots ,x_{t-k+1}\right]$。在具体计算时，可以借鉴线性回归中的偏相关系数的计算，假定 { x t } \left\{x_{t}\right\} {xt​}为一个中心化之后的平稳序列，用过去的$k$期序列值$x_{t-1},x_{t-2},\cdots ,x_{t-k}$对$x_t$坐$k$阶自回归拟合，即：
   $$x_t={\varphi }_{k1}{x}_{t-1}+{\varphi }_{k2}{x}_{t-2}+\cdots +{\varphi }_{kk}{x}_{t-k}+{\epsilon }_t$$对$x_{t-1},x_{t-2},\cdots ,x_{t-k+1}$取条件，并对方程两边取期望得：
   $$\begin{array}{rl}\hat{E}{x}_t=& {\varphi }_{k1}{x}_{t-1}+{\varphi }_{k2}{x}_{t-2}+\cdots +{\varphi }_{k(k-1)}{x}_{t-k+1}+{\varphi }_{kk}\hat{E}\left(x_{t-k}\right)+E\left({\epsilon }_t|x_{t-1},\cdots ,x_{t-k+1}\right)\end{array}$$
   由于$E\left({\epsilon }_t|x_{t-1},\cdots ,x_{t-k+1}\right)=E\left({\epsilon }_t\right)=0$，因此，直接将原式子减去上式得：
   $$x_t-\hat{E}{x}_t={\varphi }_{kk}\left(x_{t-k}-\hat{E}{x}_{t-k}\right)+{\epsilon }_t$$
   此时，再对上式两边同时乘以$x_{t-k}-\hat{E}{x}_{t-k}$，并取期望，可得：
   $$E\left[\left(x_t-\hat{E}{x}_t\right)\left(x_{t-k}-\hat{E}{x}_{t-k}\right)\right]={\varphi }_{kk}E\left[{\left(x_{t-k}-\hat{E}{x}_{t-k}\right)}^2\right]$$
   因此，可以推出：
   $${\varphi }_{kk}=\frac{E\left[\left(x_t-\hat{E}{x}_t\right)\left(x_{t-k}-\hat{E}{x}_{t-k}\right)\right]}{E\left[{\left(x_{t-k}-\hat{E}{x}_{t-k}\right)}^2\right]}$$
   即滞后$k$阶的偏自相关系数其实就等于自回归的系数。对于${\varphi }_{kk}$的求解，一般可以这样计算，即对$k阶$自回归方程两边同乘以$x_{t-l}$，得：
   $${\rho }_l={\varphi }_{k1}{\rho }_{l-1}+{\varphi }_{k2}{\rho }_{l-2}+\cdots +{\varphi }_{kk}{\rho }_{l-k},\forall l⩾1$$
   然后取前$k$个方程构成方程组：
   $$\{\begin{array}{rl}{\rho }_1& ={\varphi }_{k1}{\rho }_0+{\varphi }_{k2}{\rho }_1+\cdots +{\varphi }_{kk}{\rho }_{k-1}\\ {\rho }_2& ={\varphi }_{k1}{\rho }_1+{\varphi }_{k2}{\rho }_0+\cdots +{\varphi }_{kk}{\rho }_{k-2}\\ & ⋮\\ {\rho }_k& ={\varphi }_{k1}{\rho }_{k-1}+{\varphi }_{k2}{\rho }_{k-2}+\cdots +{\varphi }_{kk}{\rho }_0\end{array}$$该方程组称为Yule-Walker方程，用矩阵形式表示如下：
   $$\left(\begin{array}{cccc}1& {\rho }_1& \cdots & {\rho }_{k-1}\\ {\rho }_1& 1& \cdots & {\rho }_{k-2}\\ ⋮& ⋮& & ⋮\\ {\rho }_{k-1}& {\rho }_{k-2}& \cdots & 1\end{array}\right)\left(\begin{array}{c}{\varphi }_{k1}\\ {\varphi }_{k2}\\ ⋮\\ {\varphi }_{kk}\end{array}\right)=\left(\begin{array}{c}{\rho }_1\\ {\rho }_2\\ ⋮\\ {\rho }_k\end{array}\right)$$
   根据线性方程组求解的Gramer法则，有：
   $${\varphi }_{kk}=\frac{D_k}{D}$$其中，$D$为系数矩阵的行列式，$D_k$则是将$D$中第$k$列换成等号右边的系数。
   $$D=\left|\begin{array}{cccc}1& {\rho }_1& \cdots & {\rho }_{k-1}\\ {\rho }_1& 1& \cdots & {\rho }_{k-2}\\ ⋮& ⋮& & ⋮\\ {\rho }_{k-1}& {\rho }_{k-2}& \cdots & 1\end{array}\right|$$
   $$D_k=\left|\begin{array}{cccc}1& {\rho }_1& \cdots & {\rho }_1\\ {\rho }_1& 1& \cdots & {\rho }_2\\ ⋮& ⋮& & ⋮\\ {\rho }_{k-1}& {\rho }_{k-2}& \cdots & {\rho }_k\end{array}\right|$$
   对于$AR(p)$模型，其偏自相关系数具有一个特性，即$p$步截尾性，即对于$\forall k>p$，有${\varphi }_{kk}=0$。由于篇幅原因，这里证明略。

2.2 MA模型

2.2.1 MA模型的定义

    一般称具有如下结构的模型为$q$阶移动平均模型，简称$\mathrm{MA}(q)$：
$$\{\begin{array}{l}{x}_t=\mu +{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-{\theta }_2{\epsilon }_{t-2}-\cdots -{\theta }_q{\epsilon }_{t-q}\\ {\theta }_q̸=0\\ E\left({\epsilon }_t\right)=0,\mathrm{Var}\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s̸=t\end{array}$$
当$\mu =0$，则上式称为中心化$\mathrm{M}\mathrm{A}(q)$模型，采用延迟算子表示，可以表示为：
$$x_t=\Theta (B){\epsilon }_t$$
其中，$\Theta (B)=1-{\theta }_{1}B-{\theta }_2{B}^2-\cdots -{\theta }_q{B}^q$，称为$q$阶移动平均系数多项式。

2.2.2 MA模型的统计性质

1. 常数均值
   当$q<\infty$时，$\mathrm{M}\mathrm{A}(q)$的均值为：
   $$E{x}_t=E\left(\mu +{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-{\theta }_2{\epsilon }_{t-2}-\cdots -{\theta }_q{\epsilon }_{t-q}\right)=\mu$$
   特别地，对于中心化$\mathrm{M}\mathrm{A}(q)$模型，其均值为0。
2. 常数方差
   $$\mathrm{Var}\left(x_t\right)=\mathrm{Var}\left(\mu +{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-{\theta }_2{\epsilon }_{t-2}-\cdots -{\theta }_q{\epsilon }_{t-q}\right)=\left(1+{\theta }_1^2+\cdots +{\theta }_q^2\right){\sigma }_{\epsilon }^2$$
3. 自协方差函数
   $$\begin{array}{rl}{\gamma }_k& =E\left(x_t{x}_{t-k}\right)\\ & =E\left[\left({\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-\cdots -{\theta }_q{\epsilon }_{t-q}\right)\left({\epsilon }_{t-k}-{\theta }_1{\epsilon }_{t-k-1}-\cdots -{\theta }_q{\epsilon }_{t-k-q}\right)\right]\\ & =\{\begin{array}{ll}\left(1+{\theta }_1^2+\cdots +{\theta }_q^2\right){\sigma }_{\epsilon }^2,& k=0\\ \left(-{\theta }_k+{\sum }_{i=1}^{q-k}{\theta }_i{\theta }_{k+i}\right){\sigma }_{\epsilon }^2,& 1⩽k⩽q\\ 0,& k>q\end{array}\end{array}$$
4. 自相关系数q阶截尾
   $${\rho }_k=\frac{{\gamma }_k}{{\gamma }_0}=\{\begin{array}{ll}1,& k=0\\ \frac{-{\theta }_k+{\sum }_{i=1}^{q-k}{\theta }_i{\theta }_{k+i}}{1+{\theta }_1^2+\cdots +{\theta }_q^2},& 1⩽k⩽q\\ 0,& k>q\end{array}$$
5. 偏自相关系数拖尾
   由于不同的MA模型可能会有对应同一个自相关系数，因此，为了保证一个给定的自相关系数可以唯一对应一个MA模型，需要对模型添加一个约束条件，这个约束条件称为可逆性条件，即对于$\mathrm{M}\mathrm{A}(q)$模型，将其表达成$\mathrm{A}\mathrm{R}(p)$的形式：
   $${\epsilon }_t=\frac{x_t}{\Theta (B)}$$
   其中，$\Theta (B)=1-{\theta }_{1}B-\cdots -{\theta }_q{B}^q$称为移动平均系数多项式，假定$\frac{1}{{\lambda }_1},\cdots ,\frac{1}{{\lambda }_q}$为该系数多项式的$q$个根，则$\Theta (B)$可以分解为：
   $$\Theta (B)=\prod _{k=1}^q\left(1-{\lambda }_{k}B\right)$$
   因此，可以得：
   $${\epsilon }_t=\frac{x_t}{\left(1-{\lambda }_{1}B\right)\cdots \left(1-{\lambda }_{q}B\right)}$$
   要使得MA模型可逆，则要求$\left|{\lambda }_i\right|<1$，即$\mathrm{M}\mathrm{A}(q)$模型的系数多项式的根都在单位圆外，这个条件记为可逆性条件。
   这样一来，如果一个$\mathrm{M}\mathrm{A}(q)$模型满足可逆性条件，那么，它就可以有以下两种等价表达方式：
   $$\{\begin{array}{l}\Theta (B){\epsilon }_t=x_t\\ {\epsilon }_t=I(B)x_t\end{array}$$
   有上式可得：
   $$\Theta (B)I(B)x_t=x_t$$
   展开得：
   $$\left(1-\sum _{k=1}^q{\theta }_k{B}^k\right)\left(1+\sum _{j=0}^{\infty }{I}_j{B}^j\right)x_t=x_t$$
   其中，
   $$\{\begin{array}{l}{I}_0=1\\ I_j={\sum }_{k=1}^j{\theta }_k^{\prime }{I}_{j-k},j⩾1\end{array}$$
   $${\theta }_k^{\prime }=\{\begin{array}{ll}{\theta }_k,& k⩽q\\ 0,& k>q\end{array}$$
   因此，对于一个可逆$\mathrm{M}\mathrm{A}(q)$模型，可以将其等价写成$\mathrm{A}\mathrm{R}(\infty )$形式：
   $$I(B)x_t={\epsilon }_t$$
   而由前面我们知道，$\mathrm{A}\mathrm{R}(p)$模型具有$p$阶截尾，因此，$\mathrm{M}\mathrm{A}(q)$模型的偏自相关系数$\infty$截尾，即具有拖尾性。

2.3 ARMA模型

2.3.1 ARMA模型的定义

    把具有如下结构的模型称为自回归移动平均模型，简记为$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(p,q)$：
$$\{\begin{array}{l}{x}_t={\varphi }_0+{\varphi }_1{x}_{t-1}+\cdots +{\varphi }_p{x}_{t-p}+{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-\cdots -{\theta }_q{\epsilon }_{t-q}\\ {\varphi }_p̸=0,{\theta }_q̸=0\\ E\left({\epsilon }_t\right)=0,\mathrm{Var}\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s̸=t\\ E\left(x_s{\epsilon }_t\right)=0,\forall s<t\end{array}$$
当${\varphi }_0=0$时，即为中心化$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(p,q)$模型，用延迟算子可以表示为：
$$\Phi (B)x_t=\Theta (B){\epsilon }_t$$
其中，$\Phi (B)=1-{\varphi }_{1}B-\cdots -{\varphi }_p{B}^p$为$p$阶自回归系数多项式，$\Theta (B)=1-{\theta }_{1}B-\cdots -{\theta }_q{B}^q$为$q$阶移动平均系数多项式。当$q=0$或$p=0$时，则$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{A}(p,q)$则相当于$AR(p)$或$MA(q)$模型。

2.3.2 ARMA模型的统计性质

1. 均值
   对模型两边同取均值得：
   $$E{x}_t=\frac{{\varphi }_0}{1-{\varphi }_1-\cdots -{\varphi }_p}$$
2. 自协方差函数
   $$\begin{array}{rl}\gamma (k)& =E\left(x_t{x}_{t+k}\right)\\ & =E\left[\left(\sum _{i=0}^{\infty }{G}_i{\epsilon }_{t-i}\right)\left(\sum _{j=0}^{\infty }{G}_j{\epsilon }_{t+k-j}\right)\right]\\ & =E\left[\sum _{i=0}^{\infty }{G}_i\sum _{j=0}^{\infty }{G}_j{\epsilon }_{t-i}{\epsilon }_{t+k-j}\right]\\ & ={\sigma }_{\epsilon }^2\sum _{i=0}{G}_i{G}_{i+k}\end{array}$$
3. 自相关系数
   $$\rho (k)=\frac{\gamma (k)}{\gamma (0)}=\frac{{\sum }_{j=0}^{\infty }{G}_j{G}_{j+k}}{{\sum }_{j=0}^{\infty }{G}_j^2}$$
   根据$\mathrm{ARMA}(p,q)$模型的自相关系数计算公式可以知道，$\mathrm{ARMA}(p,q)$的自相关系数具有拖尾性， 并且由于$\mathrm{ARMA}(p,q)$可以转化为无穷阶的自回归模型，因此，其偏自相关系数也具有拖尾性。
   最后，我们可以总结出如下规律：
   

3. ARMA模型的python实现

     接下来，本文将以我国1950-2008年邮路及农村投递线路每年新增里程数为例，来介绍一下平稳时间序列分析的python实现。根据我们在上一篇文章《时间序列的平稳性检验与随机性检验》所介绍的，当拿到一个时间序列时，首先我们应该先观察它的时序图，根据时序图的走势做一个初步的分析，分析其是否是平稳时间序列。在python中，其代码实现如下：

import pandas as pd
def load_data(path,index_col):
   """
   加载时序数据，数据存储形式为csv格式
   :param path: 数据存储路径.[str]
   :param index_col: 时间日期所在列名. [str]
   :return:
   """
   data = pd.read_csv(path, index_col=index_col, encoding='gbk')
   return data

def plot_ts(data):
   """
   绘制时序图
   :param data:时序数据集. [DataFrame]
   :return:
   """
   plt.rcParams['font.sans-serif'] = ['SimHei']
   plt.rcParams['axes.unicode_minus'] = False
   data.plot()
   plt.show()

data = load_data('./data/arma_data.csv',index_col=u'year')
plot_ts(data)

     从上图可以发现，该时间序列的分布接近平稳时间序列的分布，接下来，我们通过绘制其自相关图来做进一步判断，代码如下：

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf as acf


def plot_acf(data):
   """
   绘制自相关图
   :param data:时序数据集.[DataFrame]
   :return:
   """
   acf(data)
   plt.show()

     从该时间序列的自相关图可以发现，在第三阶之后，其自相关系数基本收敛到0附近，并随着阶数的增加，自相关系数一直在0附近上下波动，因此，从自相关图也可以初步判定在时间序列是一个平稳时间序列。那么，当判定时间序列为平稳时间序列后，我们还要进一步判断该时间序列是否是白噪声序列，如果是白噪声序列，则此时对时间序列的分析将到此为此。采用LB统计量进行白噪声检验：

from statsmodels.stats.diagnostic import acorr_ljungbox


def LB_test(data,lags=10):
    """
    白噪声检验，采用LB统计量检验
    :param data: 时序数据集.[DataFrame]
    :param lags: 阶数.[int]
    :return:
    """
    result = acorr_ljungbox(data, lags=lags)
    return result

     可以发现，在0.05的显著性水平下，此时前10阶的LB统计量对应的p值均小于0.05，因此，我们可以认为该时间序列不是白噪声序列。接下来，对于该时间序列，我们采用$ARMA(p,q)$进行拟合，但是需要确定$p$和$q$的阶数到底取多少比较合适。我们可以直接通过自相关图和偏自相关图来进行判定：

from statsmodels.graphics.tsaplots import plot_pacf as pacf

def plot_pacf(data,lags=20):
    """
    绘制偏自相关图
    :param data: 时序数据集.[DataFrame]
    :param lags: 阶数.[int]
    :return:
    """
    pacf(data,lags=lags)
    plt.show()

     从该时间序列的自相关图和偏自相关图可以发现，该时间序列的自相关图在3阶后逐渐收敛至0，但是一直围绕0上下波动，而不是跳高式地迅速缩减到0附近，因此是我们所说的拖尾的状态，而偏自相关图则在2阶之后迅速收敛到0附近，因此，可以认为是截尾状态，因此，我们可以将$p$定为2，而$q$定为0，即采用$AR(2)$模型进行拟合。不过，通过自相关图和偏自相关图有时比较难做出判断，此时，也可以通过模型的BIC值来确定$p$和$q$的最优取值，最终选取出来的阶数与通过自相关图和偏自相关图得出的分析结果一致。

from statsmodels.api import tsa

def best_lags(data,max_p=None,max_q=None):
    """
    采用BIC值确定arma模型的最佳阶数
    :param data: 时序数据集. [DataFrame]
    :param max_p: 最大的p值. [int]
    :param max_q: 最大的q值. [int]
    :return:
    """
    # 一般阶数不超过length/10
    if max_p == None:
        max_p = int(len(data) / 10)
    if max_q == None:
        max_q = int(len(data) / 10)
    bic_matrix = []
    for p in range(max_p + 1):
        tmp = []
        for q in range(max_q + 1):
            try:
                tmp.append(tsa.ARMA(data, (p, q)).fit().bic)
            except:
                tmp.append(None)
        bic_matrix.append(tmp)
    bic_matrix = pd.DataFrame(bic_matrix)
    p, q = bic_matrix.stack().idxmin()
    return p,q

     定阶后，我们根据阶数拟合$ARMA$模型，最终结果如下：

from statsmodels.api import tsa

best_p,best_q = best_lags(data)
arma = tsa.ARMA(data,(best_p,best_q)).fit()
arma.summary2()

最终拟合出来的模型如下：
$$x_t-11.0227=\frac{{\epsilon }_t}{1-0.7185B+0.5294B^2}$$

4. 总结

     以上就是对平稳时间序列模型$AR、MA、ARMA$模型的介绍以及python实现，整体涉及到的知识点还是蛮多的，不过需要注意的是，在现实生活中，平稳时间序列其实是很少出现的，因此，对于非平稳时间序列，这三种模型都不能直接使用，关于非平稳时间序列的分析我们将在后续的文章中进行介绍。