平稳时间序列分析与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 + ⋯ + a p z t − p = h ( t ) z _ { t } + a _ { 1 } z _ { t - 1 } + a _ { 2 } z _ { t - 2 } + \cdots + a _ { p } z _ { t - p } = h ( t ) zt+a1zt1+a2zt2++apztp=h(t)其中, p ⩾ 1 ; a 1 , a 2 , ⋯   , a p p \geqslant 1 ; a _ { 1 } , a _ { 2 } , \cdots , a _ { p } p1;a1,a2,,ap为实数, h ( t ) h ( t ) h(t) t t t的已知函数。当 h ( t ) = 0 h ( t )=0 h(t)=0时,则称差分方程为齐次线性差分方程,否则称为非齐次线性差分方程

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

    齐次线性差分方程的求解需要借助特征方程和特征根,其特征方程为:
λ p + a 1 λ p − 1 + a 2 λ p − 2 + ⋯ + a p = 0 \lambda ^ { p } + a _ { 1 } \lambda ^ { p - 1 } + a _ { 2 } \lambda ^ { p - 2 } + \cdots + a _ { p } = 0 λp+a1λp1+a2λp2++ap=0,特征根为该特征方程对应的根,不妨记作
λ 1 , λ 2 , ⋯   , λ p \lambda _ { 1 } , \lambda _ { 2 } , \cdots , \lambda _ { p } λ1,λ2,,λp根据特征根的取值情况,齐次线性差分方程的通解会有不同的表达方式,下面分三种情况列举:

  • λ 1 , λ 2 , ⋯   , λ p \lambda _ { 1 } , \lambda _ { 2 } , \cdots , \lambda _ { p } λ1,λ2,,λp无重根时,此时齐次线性差分方程的通解为:
    z t = c 1 λ 1 t + c 2 λ 2 t + ⋯ + c p λ p t z _ { t } = c _ { 1 } \lambda _ { 1 } ^ { t } + c _ { 2 } \lambda _ { 2 } ^ { t } + \cdots + c _ { p } \lambda _ { p } ^ { t } zt=c1λ1t+c2λ2t++cpλpt其中, c 1 , c 2 , ⋯   , c p c _ { 1 } , c _ { 2 } , \cdots , c _ { p } c1,c2,,cp为任意的实数。
  • λ 1 , λ 2 , ⋯   , λ p \lambda _ { 1 } , \lambda _ { 2 } , \cdots , \lambda _ { p } λ1,λ2,,λp有重根时,不妨假定 λ 1 = λ 2 = ⋯ = λ d \lambda _ { 1 } = \lambda _ { 2 } = \dots = \lambda _ { d } λ1=λ2==λd d d d个相同的实根, λ d + 1 , λ d + 2 , ⋯   , λ p \lambda _ { d + 1 } , \lambda _ { d + 2 } , \cdots , \lambda _ { p } λd+1,λd+2,,λp为其他不相同的实根,则此时齐次线性差分方程的通解为:
    z t = ( c 1 + c 2 t + ⋯ + c d t d − 1 ) λ 1 t + c d + 1 λ d + 1 t + ⋯ + c p λ p t 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 } zt=(c1+c2t++cdtd1)λ1t+cd+1λd+1t++cpλpt其中, c 1 , c 2 , ⋯   , c p c _ { 1 } , c _ { 2 } , \cdots , c _ { p } c1,c2,,cp为任意的实数。
  • λ 1 , λ 2 , ⋯   , λ p \lambda _ { 1 } , \lambda _ { 2 } , \cdots , \lambda _ { p } λ1,λ2,,λp有复根时,不妨设 λ 1 = a + i b = r e i ω , λ 2 = a − i b = r e − i ω \lambda _ { 1 } = a + i b = r \mathrm { e } ^ { i \omega } , \lambda _ { 2 } = a - i b = r \mathrm { e } ^ { - i \omega } λ1=a+ib=reiω,λ2=aib=reiω为其中的一对复根,由欧拉公式易得 r = a 2 + b 2 , ω = arccos ⁡ a r r = \sqrt { a ^ { 2 } + b ^ { 2 } } , \omega = \arccos \frac { a } { r } r=a2+b2 ,ω=arccosra,而 λ 3 , λ 4 , ⋯   , λ p \lambda _ { 3 } , \lambda _ { 4 } , \cdots , \lambda _ { p } λ3,λ4,,λp为其中不相同的实根,则此时齐次线性差分方程的通解为:
    z t = c 1 λ 1 t + c 2 λ 2 t + ⋯ + c p λ p t = r t ( c 1 e i t ω + c 2 e − i t ω ) + c 3 λ 3 t + ⋯ + c p λ p t \begin{aligned} 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 } ^ { i t \omega } + c _ { 2 } \mathrm { e } ^ { - i t \omega } \right) + c _ { 3 } \lambda _ { 3 } ^ { t } + \cdots + c _ { p } \lambda _ { p } ^ { t } \end{aligned} zt=c1λ1t+c2λ2t++cpλpt=rt(c1eitω+c2eitω)+c3λ3t++cpλpt其中, c 1 , c 2 , ⋯   , c p c _ { 1 } , c _ { 2 } , \cdots , c _ { p } c1,c2,,cp为任意的实数。
2.1.3 非齐次线性差分方程的通解

    非齐次线性差分方程的通解为对应的齐次线性差分方程的通解 z t ′ z _ { t } ^ { \prime } zt和一个特解 z t ′ ′ z _ { t } ^ { \prime \prime } zt的加和,即:
z t = z t ′ + z t ′ ′ z _ { t } = z _ { t } ^ { \prime } + z _ { t } ^ { \prime \prime } zt=zt+zt其中,特解 z t ′ ′ z _ { t } ^ { \prime \prime } zt就是任意一个使得非齐次线性差分方程成立的值,即:
z t ′ ′ + a 1 z t − 1 ′ ′ + a 2 z t − 2 ′ ′ + ⋯ + a p z t − p ′ ′ = h ( t ) 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 ) zt+a1zt1+a2zt2++apztp=h(t)

2.2 AR模型
2.2.1 AR模型的定义

    AR模型也称为自回归模型,其定义如下:
{ x t = ϕ 0 + ϕ 1 x t − 1 + ϕ 2 x t − 2 + ⋯ + ϕ p x t − p + ε t ϕ p ≠ 0 E ( ε t ) = 0 , Var ⁡ ( ε t ) = σ ε 2 , E ( ε t ε s ) = 0 , s ≠ t E ( x s ε t ) = 0 , ∀ s &lt; t \left\{ \begin{array} { l } { x _ { t } = \phi _ { 0 } + \phi _ { 1 } x _ { t - 1 } + \phi _ { 2 } x _ { t - 2 } + \cdots + \phi _ { p } x _ { t - p } + \varepsilon _ { t } } \\ { \phi _ { p } \neq 0 } \\ { E \left( \varepsilon _ { t } \right) = 0 , \operatorname { Var } \left( \varepsilon _ { t } \right) = \sigma _ { \varepsilon } ^ { 2 } , E \left( \varepsilon _ { t } \varepsilon _ { s } \right) = 0 , \quad s \neq t } \\ { E \left( x _ { s } \varepsilon _ { t } \right) = 0 , \quad \forall s &lt; t } \end{array} \right. xt=ϕ0+ϕ1xt1+ϕ2xt2++ϕpxtp+εtϕp̸=0E(εt)=0,Var(εt)=σε2,E(εtεs)=0,s̸=tE(xsεt)=0,s<t上述的结构也称为 p p p阶自回归模型,简记为 A R ( p ) AR(p) AR(p)
    当对序列进行中心化后,并将AR模型用延迟算子表示时,此时AR模型可以表示为:
Φ ( B ) x t = ε t \Phi ( B ) x _ { t } = \varepsilon _ { t } Φ(B)xt=εt其中, Φ ( B ) = 1 − ϕ 1 B − ϕ 2 B 2 − ⋯ − ϕ p B p \Phi ( B ) = 1 - \phi _ { 1 } B - \phi _ { 2 } B ^ { 2 } - \cdots - \phi _ { p } B ^ { p } Φ(B)=1ϕ1Bϕ2B2ϕpBp,称为 p p p阶自回归系数多项式。

2.2.2 AR模型平稳性的判别

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

  1. 特征根判别
    对于任一中心化 A R ( p ) \mathrm { AR } ( p ) AR(p)模型 Φ ( B ) x t = ε t \Phi ( B ) x _ { t } = \varepsilon _ { t } Φ(B)xt=εt都可以被视为一个非齐性线性差分方程:
    x t − ϕ 1 x t − 1 − ϕ 2 x t − 2 − ⋯ − ϕ p x t − p = ε t x _ { t } - \phi _ { 1 } x _ { t - 1 } - \phi _ { 2 } x _ { t - 2 } - \cdots - \phi _ { p } x _ { t - p } = \varepsilon _ { t } xtϕ1xt1ϕ2xt2ϕpxtp=εt其通解为:
    x t = x t ′ + x t ′ ′ x _ { t } = x _ { t } ^ { \prime } + x _ { t } ^ { \prime \prime } xt=xt+xt其中, x t ′ x _ { t } ^ { \prime } xt为齐次线性差分方程 Φ ( B ) x t = 0 \Phi ( B ) x _ { t } = 0 Φ(B)xt=0的通解, x t ′ ′ x _ { t } ^ { \prime \prime } xt为特解。
    对于齐次线性差分方程 Φ ( B ) x t = 0 \Phi ( B ) x _ { t } = 0 Φ(B)xt=0的通解,记其特征根为 λ 1 , λ 2 , ⋯ &ThinSpace; , λ p \lambda _ { 1 } , \lambda _ { 2 } , \cdots , \lambda _ { p } λ1,λ2,,λp,为了不失一般性,我们假设有 d d d个相同的实根, p − d − 2 m p-d-2m pd2m个不相等的实根, m m m对复根,则由前面介绍的线性差分方程我们可以得到通解的计算公式如下:
    x t ′ = ∑ j = 1 d c j t j − 1 λ 1 t + ∑ j = d + 1 p − 2 m c j λ j t + ∑ j = 1 m r j t ( c 1 j e i t w j + c 2 j e − i t w j ) x _ { t } ^ { \prime } = \sum _ { j = 1 } ^ { d } c _ { j } t ^ { j - 1 } \lambda _ { 1 } ^ { t } + \sum _ { j = d + 1 } ^ { p - 2 m } c _ { j } \lambda _ { j } ^ { t } + \sum _ { j = 1 } ^ { m } r _ { j } ^ { t } \left( c _ { 1 j } \mathrm { e } ^ { i t w _ { j } } + c _ { 2 j } \mathrm { e } ^ { - i t w _ { j } } \right) xt=j=1dcjtj1λ1t+j=d+1p2mcjλjt+j=1mrjt(c1jeitwj+c2jeitwj)其中, c 1 , ⋯ &ThinSpace; , c p − 2 m , c 1 j , c 2 j ( j = 1 , ⋯ &ThinSpace; , m ) c _ { 1 } , \cdots , c _ { p - 2 m } , c _ { 1 j } , c _ { 2 j } ( j = 1 , \cdots , m ) c1,,cp2m,c1j,c2j(j=1,,m)为任意实数。
    对于特解,首先,可以证明自回归系数多项式方程 Φ ( u ) = 0 \Phi ( u ) = 0 Φ(u)=0的根是齐次线性差分方程 Φ ( B ) x t = 0 \Phi ( B ) x _ { t } = 0 Φ(B)xt=0的特征根的倒数。因此, Φ ( B ) \Phi ( B ) Φ(B)可以因子分解为:
    Φ ( B ) = ∏ i = 1 p ( 1 − λ i B ) \Phi ( B ) = \prod _ { i = 1 } ^ { p } \left( 1 - \lambda _ { i } B \right) Φ(B)=i=1p(1λiB)这样一来,特解 x t ′ ′ x _ { t } ^ { \prime \prime } xt的计算公式可以表示为:
    x t ′ ′ = ε t Φ ( B ) = ε t ∏ i = 1 p ( 1 − λ i B ) = ∑ i = 1 p k i 1 − λ i B ε t x _ { t } ^ { \prime \prime } = \frac { \varepsilon _ { t } } { \Phi ( B ) } = \frac { \varepsilon _ { t } } { \prod _ { i = 1 } ^ { p } \left( 1 - \lambda _ { i } B \right) } = \sum _ { i = 1 } ^ { p } \frac { k _ { i } } { 1 - \lambda _ { i } B } \varepsilon _ { t } xt=Φ(B)εt=i=1p(1λiB)εt=i=1p1λiBkiεt其中, k i ( i = 1 , 2 , ⋯ &ThinSpace; , p ) k _ { i } ( i = 1,2 , \cdots , p ) ki(i=1,2,,p)为常数。
    因此,最终 A R ( p ) \mathrm { AR } ( p ) AR(p)模型的通解的计算公式为:
    x t = x t ′ + x t ′ ′ = ∑ j = 1 d c j t j − 1 λ 1 t + ∑ j = d + 1 p − 2 m c j λ j t + ∑ j = 1 m r j t ( c 1 j e i t w j + c 2 j e − i t w j ) + ∑ i = 1 p k i 1 − λ i B ε t \begin{aligned} x _ { t } &amp; = x _ { t } ^ { \prime } + x _ { t } ^ { \prime \prime } \\ &amp; = \sum _ { j = 1 } ^ { d } c _ { j } t ^ { j - 1 } \lambda _ { 1 } ^ { t } + \sum _ { j = d + 1 } ^ { p - 2 m } c _ { j } \lambda _ { j } ^ { t } + \sum _ { j = 1 } ^ { m } r _ { j } ^ { t } \left( c _ { 1 j } \mathrm { e } ^ { i t w _ { j } } + c _ { 2 j } \mathrm { e } ^ { - i t w _ { j } } \right) + \sum _ { i = 1 } ^ { p } \frac { k _ { i } } { 1 - \lambda _ { i } B } \varepsilon _ { t } \end{aligned} xt=xt+xt=j=1dcjtj1λ1t+j=d+1p2mcjλjt+j=1mrjt(c1jeitwj+c2jeitwj)+i=1p1λiBkiεt
    要使得 A R ( p ) \mathrm { AR } ( p ) AR(p)模型平稳,即要求对于任意的实数 c 1 , ⋯ &ThinSpace; , c p − 2 m , c 1 j , c 2 j ( j = 1 , 2 , ⋯ &ThinSpace; , m ) c _ { 1 } , \cdots , c _ { p - 2 m } , c _ { 1 j } , c _ { 2 j } ( j = 1,2 , \cdots,m) c1,,cp2m,c1j,c2j(j=1,2,,m)有:
    lim ⁡ t → ∞ x t = 0 \lim _ { t \rightarrow \infty } x _ { t } = 0 tlimxt=0要使得上式成立,其充要条件为:
    ∣ λ i ∣ &lt; 1 , i = 1 , 2 , ⋯ &ThinSpace; , p − 2 m ∣ r i ∣ &lt; 1 , i = 1 , 2 , ⋯ &ThinSpace; , m \begin{array} { l } { \left| \lambda _ { i } \right| &lt; 1 , i = 1,2 , \cdots , p - 2 m } \\ { \left| r _ { i } \right| &lt; 1 , i = 1,2 , \cdots , m } \end{array} λi<1,i=1,2,,p2mri<1,i=1,2,,m A R ( p ) \mathrm { AR } ( p ) AR(p) p p p个特征根都在单位圆内,或者自回归系数多项式的根都在单位圆外。
  2. 平稳域判别
    平稳域判别则直接根据自回归系数进行判断,即对参数向量 ( ϕ 1 , ϕ 2 , ⋯ &ThinSpace; , ϕ p ) ′ \left( \phi _ { 1 } , \phi _ { 2 } , \cdots , \phi _ { p } \right) ^ { \prime } (ϕ1,ϕ2,,ϕp)的取值范围进行限制,使得特征根都在单位圆内,即:
    { ϕ 1 , ϕ 2 , ⋯ &ThinSpace; , ϕ p ∣ 特 征 根 都 在 单 位 圆 内 } \left\{ \phi _ { 1 } , \phi _ { 2 } , \cdots , \phi _ { p } |\right.特征根都在单位圆内\} {ϕ1,ϕ2,,ϕp}该集合即称为 A R ( p ) \mathrm { AR } ( p ) AR(p)模型的平稳域,一般这种方法对于低阶 A R AR AR模型比较方便。
2.2.3 平稳AR模型的统计性质
  1. 均值
    对于 A R ( p ) AR(p) AR(p)模型,在等式两边取期望得:
    E x t = E ( ϕ 0 + ϕ 1 x t − 1 + ϕ 2 x t − 2 + ⋯ + ϕ p x t − p + ε t ) E x _ { t } = E \left( \phi _ { 0 } + \phi _ { 1 } x _ { t - 1 } + \phi _ { 2 } x _ { t - 2 } + \cdots + \phi _ { p } x _ { t - p } + \varepsilon _ { t } \right) Ext=E(ϕ0+ϕ1xt1+ϕ2xt2++ϕpxtp+εt),由于平稳序列的均值为常数 E x t = μ ( ∀ t ∈ T ) E x _ { t } = \mu ( \forall t \in T ) Ext=μ(tT),并且 { ε t } \left\{ \varepsilon _ { t } \right\} {εt}为白噪声序列,其均值 E ε t = 0 E \varepsilon _ { t } = 0 Eεt=0,所以有:
    ( 1 − ϕ 1 − … − ϕ p ) μ = ϕ 0 ⇒ μ = ϕ 0 1 − ϕ 1 − ⋯ − ϕ p \begin{aligned} &amp; \left( 1 - \phi _ { 1 } - \ldots - \phi _ { p } \right) \mu = \phi _ { 0 } \\ \Rightarrow \mu &amp; = \frac { \phi _ { 0 } } { 1 - \phi _ { 1 } - \cdots - \phi _ { p } } \end{aligned} μ(1ϕ1ϕp)μ=ϕ0=1ϕ1ϕpϕ0
  2. 方差
    要计算 A R ( p ) AR(p) AR(p)模型的方差,需要引入Green函数,其定义如下,设 λ 1 , λ 2 , ⋯ &ThinSpace; , λ p \lambda _ { 1 } , \lambda _ { 2 } , \cdots , \lambda _ { p } λ1,λ2,,λp为平稳 A R ( p ) AR(p) AR(p)模型的特征根,则有:
    x t = ε t Φ ( B ) = ∑ i = 1 p k i 1 − λ i B ε t = ∑ i = 1 ∞ ∑ j = 0 ∞ k i ( λ i B ) j ε t = ∑ j = 0 ∞ ∑ i = 1 p k i λ i j ε t − j ≡ ∑ j = 0 ∞ G j ε t − j \begin{aligned} x _ { t } &amp; = \frac { \varepsilon _ { t } } { \Phi ( B ) } \\ &amp; = \sum _ { i = 1 } ^ { p } \frac { k _ { i } } { 1 - \lambda _ { i } B } \varepsilon _ { t } \\ &amp; = \sum _ { i = 1 } ^ { \infty } \sum _ { j = 0 } ^ { \infty } k _ { i } \left( \lambda _ { i } B \right) ^ { j } \varepsilon _ { t } \\ &amp; = \sum _ { j = 0 } ^ { \infty } \sum _ { i = 1 } ^ { p } k _ { i } \lambda _ { i } ^ { j } \varepsilon _ { t - j } \\ &amp; \equiv \sum _ { j = 0 } ^ { \infty } G _ { j } \varepsilon _ { t - j } \end{aligned} xt=Φ(B)εt=i=1p1λiBkiεt=i=1j=0ki(λiB)jεt=j=0i=1pkiλijεtjj=0Gjεtj其中, G 0 = 1 , G j = ∑ i = 1 p k i λ i j ( j = 1 , 2 , ⋯ &ThinSpace; ) G _ { 0 } = 1 , G _ { j } = \sum _ { i = 1 } ^ { p } k _ { i } \lambda _ { i } ^ { j } ( j = 1,2 , \cdots ) G0=1,Gj=i=1pkiλij(j=1,2,)
    因此, A R ( p ) AR(p) AR(p)模型的方差为:
    Var ⁡ ( x t ) = ∑ j = 0 ∞ G j 2 Var ⁡ ( ε t ) \operatorname { Var } \left( x _ { t } \right) = \sum _ { j = 0 } ^ { \infty } G _ { j } ^ { 2 } \operatorname { Var } \left( \varepsilon _ { t } \right) Var(xt)=j=0Gj2Var(εt)
  3. 协方差函数
    在平稳模型 x t = ϕ 1 x t − 1 + ϕ 2 x t − 2 + ⋯ + ϕ p x t − p + ε t x _ { t } = \phi _ { 1 } x _ { t - 1 } + \phi _ { 2 } x _ { t - 2 } + \cdots + \phi _ { p } x _ { t - p } + \varepsilon _ { t } xt=ϕ1xt1+ϕ2xt2++ϕpxtp+εt两边同乘以 x t − k ( ∀ k ⩾ 1 ) x _ { t - k } ( \forall k \geqslant 1 ) xtk(k1),并求期望得:
    E ( x t x t − k ) = ϕ 1 E ( x t − 1 x t − k ) + ⋯ + ϕ p E ( x t − p x t − k ) + E ( ε t x t − k ) , ∀ k ⩾ 1 E \left( x _ { t } x _ { t - k } \right) = \phi _ { 1 } E \left( x _ { t - 1 } x _ { t - k } \right) + \cdots + \phi _ { p } E \left( x _ { t - p } x _ { t - k } \right) + E \left( \varepsilon _ { t } x _ { t - k } \right) , \quad \forall k \geqslant 1 E(xtxtk)=ϕ1E(xt1xtk)++ϕpE(xtpxtk)+E(εtxtk),k1由于 E ( ε t x t − k ) = 0 , ∀ k ⩾ 1 E \left( \varepsilon _ { t } x _ { t - k } \right) = 0 , \quad \forall k \geqslant 1 E(εtxtk)=0,k1,因此,可得自协方差函数的递推公式:
    γ k = ϕ 1 γ k − 1 + ϕ 2 γ k − 2 + ⋯ + ϕ p γ k − p \gamma _ { k } = \phi _ { 1 } \gamma _ { k - 1 } + \phi _ { 2 } \gamma _ { k - 2 } + \cdots + \phi _ { p } \gamma _ { k - p } γk=ϕ1γk1+ϕ2γk2++ϕpγkp
  4. 自相关系数具有拖尾性
    由自相关系数的计算公式 ρ k = γ k γ 0 \rho _ { k } = \frac { \gamma _ { k } } { \gamma _ { 0 } } ρk=γ0γk易得 A R ( p ) AR(p) AR(p)模型自相关系数的递推公式:
    ρ k = ϕ 1 ρ k − 1 + ϕ 2 ρ k − 2 + ⋯ + ϕ p ρ k − p \rho _ { k } = \phi _ { 1 } \rho _ { k - 1 } + \phi _ { 2 } \rho _ { k - 2 } + \cdots + \phi _ { p } \rho _ { k - p } ρk=ϕ1ρk1+ϕ2ρk2++ϕpρkp由自相关系数的递推公式可以发现,其实该表达式是一个 p p p阶齐次线性差分方程,其滞后任意 k k k阶的自相关系数的通解为:
    ρ k = ∑ i = 1 p c i λ i k \rho _ { k } = \sum _ { i = 1 } ^ { p } c _ { i } \lambda _ { i } ^ { k } ρk=i=1pciλik式中, ∣ λ i ∣ &lt; 1 ( i = 1 , 2 , ⋯ &ThinSpace; , p ) \left| \lambda _ { i } \right| &lt; 1 ( i = 1,2 , \cdots , p ) λi<1(i=1,2,,p)为该差分方程的特征根, c 1 , c 2 , ⋯ &ThinSpace; , c p c _ { 1 } , c _ { 2 } , \cdots , c _ { p } c1,c2,,cp为任意实数,如果已知任意 p p p个自相关系数,那么就可以求出 c 1 , c 2 , ⋯ &ThinSpace; , c p c _ { 1 } , c _ { 2 } , \cdots , c _ { p } c1,c2,,cp对应的值,因此, c 1 , c 2 , ⋯ &ThinSpace; , c p c _ { 1 } , c _ { 2 } , \cdots , c _ { p } c1,c2,,cp不可能全为0。因此,由该式也容易得到 ρ k \rho _ { k } ρk始终有非0取值,即不会在 k k k大于某个值后就恒等于0,这也是 A R AR AR模型自相关系数的一个重要特性,即拖尾性
    另外,随着 k k k的增大,由于 ∣ λ i ∣ &lt; 1 ( i = 1 , 2 , ⋯ &ThinSpace; , p ) \left| \lambda _ { i } \right| &lt; 1 ( i = 1,2 , \cdots , p ) λi<1(i=1,2,,p),所以当 k → ∞ k \rightarrow \infty k时, ρ k = ∑ i = 1 p c i λ i k → 0 \rho _ { k } = \sum _ { i = 1 } ^ { p } c _ { i } \lambda _ { i } ^ { k } \rightarrow 0 ρk=i=1pciλik0,因此,可以发现 A R AR AR模型自相关系数的另一个重要特性,即指数衰减性,这也是我们前面讲的为什么平稳时间序列的自相关系数很快会衰减到0附近。
  5. 偏自相关系数具有p步截尾性
    对于一个平稳 A R ( p ) AR(p) AR(p)模型,前面所计算出的滞后 k k k阶的自相关系数 ρ k \rho_{k} ρk实际上并不是 x t x_t xt x t − k x_{t-k} xtk之间单纯的相关关系,因为 x t x_t xt还会受到中间 k − 1 k-1 k1个随机变量 x t − 1 , x t − 2 , ⋯ &ThinSpace; , x t − k + 1 x_{t-1}, x_{t-2}, \cdots, x_{t-k+1} xt1,xt2,,xtk+1的影响,因此,为了能单纯地测度 x t − k x_{t-k} xtk x t x_t xt的影响,引进了偏自相关系数的概念,其定义如下:
    ρ x t , x t − k ∣ x t − 1 , ⋯ &ThinSpace; , x t − k + 1 = E [ ( x t − E ^ x t ) ( x t − k − E ^ x t − k ) ] E [ ( x t − k − E ^ x t − k ) 2 ] \rho_{x_{t}, x_{t-k}\left|x_{t-1}, \cdots, x_{t-k+1}\right.}=\frac{E\left[\left(x_{t}-\widehat{E} x_{t}\right)\left(x_{t-k}-\widehat{E} x_{t-k}\right)\right]}{E\left[\left(x_{t-k}-\widehat{E} x_{t-k}\right)^{2}\right]} ρxt,xtkxt1,,xtk+1=E[(xtkE xtk)2]E[(xtE xt)(xtkE xtk)]
    其中, E ^ x t = E [ x t ∣ x t − 1 , ⋯ &ThinSpace; , x t − k + 1 ] \widehat{E} x_{t}=E\left[x_{t} | x_{t-1}, \cdots, x_{t-k+1}\right] E xt=E[xtxt1,,xtk+1] E ^ x t − k = E [ x t − k ∣ x t − 1 , ⋯ &ThinSpace; , x t − k + 1 ] \widehat{E} x_{t-k}=E\left[x_{t-k} | x_{t-1}, \cdots, x_{t-k+1}\right] E xtk=E[xtkxt1,,xtk+1]。在具体计算时,可以借鉴线性回归中的偏相关系数的计算,假定 { x t } \left\{x_{t}\right\} {xt}为一个中心化之后的平稳序列,用过去的 k k k期序列值 x t − 1 , x t − 2 , ⋯ &ThinSpace; , x t − k x_{t-1}, x_{t-2}, \cdots, x_{t-k} xt1,xt2,,xtk x t x_t xt k k k阶自回归拟合,即:
    x t = ϕ k 1 x t − 1 + ϕ k 2 x t − 2 + ⋯ + ϕ k k x t − k + ε t x_{t}=\phi_{k 1} x_{t-1}+\phi_{k 2} x_{t-2}+\cdots+\phi_{k k} x_{t-k}+\varepsilon_{t} xt=ϕk1xt1+ϕk2xt2++ϕkkxtk+εt x t − 1 , x t − 2 , ⋯ &ThinSpace; , x t − k + 1 x_{t-1}, x_{t-2}, \cdots, x_{t-k+1} xt1,xt2,,xtk+1取条件,并对方程两边取期望得:
    E ^ x t = ϕ k 1 x t − 1 + ϕ k 2 x t − 2 + ⋯ + ϕ k ( k − 1 ) x t − k + 1 + ϕ k k E ^ ( x t − k ) + E ( ε t ∣ x t − 1 , ⋯ &ThinSpace; , x t − k + 1 ) \begin{aligned} \widehat{E} x_{t}=&amp; \phi_{k 1} x_{t-1}+\phi_{k 2} x_{t-2}+\cdots+\phi_{k(k-1)} x_{t-k+1}+\phi_{k k} \widehat{E}\left(x_{t-k}\right) +E\left(\varepsilon_{t} | x_{t-1}, \cdots, x_{t-k+1}\right) \end{aligned} E xt=ϕk1xt1+ϕk2xt2++ϕk(k1)xtk+1+ϕkkE (xtk)+E(εtxt1,,xtk+1)
    由于 E ( ε t ∣ x t − 1 , ⋯ &ThinSpace; , x t − k + 1 ) = E ( ε t ) = 0 E\left(\varepsilon_{t} | x_{t-1}, \cdots, x_{t-k+1}\right)=E\left(\varepsilon_{t}\right)=0 E(εtxt1,,xtk+1)=E(εt)=0,因此,直接将原式子减去上式得:
    x t − E ^ x t = ϕ k k ( x t − k − E ^ x t − k ) + ε t x_{t}-\widehat{E} x_{t}=\phi_{k k}\left(x_{t-k}-\widehat{E} x_{t-k}\right)+\varepsilon_{t} xtE xt=ϕkk(xtkE xtk)+εt
    此时,再对上式两边同时乘以 x t − k − E ^ x t − k x_{t-k}-\widehat{E} x_{t-k} xtkE xtk,并取期望,可得:
    E [ ( x t − E ^ x t ) ( x t − k − E ^ x t − k ) ] = ϕ k k E [ ( x t − k − E ^ x t − k ) 2 ] E\left[\left(x_{t}-\widehat{E} x_{t}\right)\left(x_{t-k}-\widehat{E} x_{t-k}\right)\right]=\phi_{k k} E\left[\left(x_{t-k}-\widehat{E} x_{t-k}\right)^{2}\right] E[(xtE xt)(xtkE xtk)]=ϕkkE[(xtkE xtk)2]
    因此,可以推出:
    ϕ k k = E [ ( x t − E ^ x t ) ( x t − k − E ^ x t − k ) ] E [ ( x t − k − E ^ x t − k ) 2 ] \phi_{k k}=\frac{E\left[\left(x_{t}-\widehat{E} x_{t}\right)\left(x_{t-k}-\widehat{E} x_{t-k}\right)\right]}{E\left[\left(x_{t-k}-\widehat{E} x_{t-k}\right)^{2}\right]} ϕkk=E[(xtkE xtk)2]E[(xtE xt)(xtkE xtk)]
    即滞后 k k k阶的偏自相关系数其实就等于自回归的系数。对于 ϕ k k \phi_{k k} ϕkk的求解,一般可以这样计算,即对 k 阶 k阶 k自回归方程两边同乘以 x t − l x_{t-l} xtl,得:
    ρ l = ϕ k 1 ρ l − 1 + ϕ k 2 ρ l − 2 + ⋯ + ϕ k k ρ l − k , ∀ l ⩾ 1 \rho_{l}=\phi_{k 1} \rho_{l-1}+\phi_{k 2} \rho_{l-2}+\cdots+\phi_{k k} \rho_{l-k}, \quad \forall l \geqslant 1 ρl=ϕk1ρl1+ϕk2ρl2++ϕkkρlk,l1
    然后取前 k k k个方程构成方程组:
    { ρ 1 = ϕ k 1 ρ 0 + ϕ k 2 ρ 1 + ⋯ + ϕ k k ρ k − 1 ρ 2 = ϕ k 1 ρ 1 + ϕ k 2 ρ 0 + ⋯ + ϕ k k ρ k − 2 ⋮ ρ k = ϕ k 1 ρ k − 1 + ϕ k 2 ρ k − 2 + ⋯ + ϕ k k ρ 0 \left\{\begin{aligned} \rho_{1} &amp;=\phi_{k 1} \rho_{0}+\phi_{k 2} \rho_{1}+\cdots+\phi_{k k} \rho_{k-1} \\ \rho_{2} &amp;=\phi_{k 1} \rho_{1}+\phi_{k 2} \rho_{0}+\cdots+\phi_{k k} \rho_{k-2} \\ &amp; \vdots \\ \rho_{k} &amp;=\phi_{k 1} \rho_{k-1}+\phi_{k 2} \rho_{k-2}+\cdots+\phi_{k k} \rho_{0} \end{aligned}\right. ρ1ρ2ρk=ϕk1ρ0+ϕk2ρ1++ϕkkρk1=ϕk1ρ1+ϕk2ρ0++ϕkkρk2=ϕk1ρk1+ϕk2ρk2++ϕkkρ0该方程组称为Yule-Walker方程,用矩阵形式表示如下:
    ( 1 ρ 1 ⋯ ρ k − 1 ρ 1 1 ⋯ ρ k − 2 ⋮ ⋮ ⋮ ρ k − 1 ρ k − 2 ⋯ 1 ) ( ϕ k 1 ϕ k 2 ⋮ ϕ k k ) = ( ρ 1 ρ 2 ⋮ ρ k ) \left( \begin{array}{cccc}{1} &amp; {\rho_{1}} &amp; {\cdots} &amp; {\rho_{k-1}} \\ {\rho_{1}} &amp; {1} &amp; {\cdots} &amp; {\rho_{k-2}} \\ {\vdots} &amp; {\vdots} &amp; { } &amp; {\vdots} \\ {\rho_{k-1}} &amp; {\rho_{k-2}} &amp; {\cdots} &amp; {1}\end{array}\right) \left( \begin{array}{c}{\phi_{k 1}} \\ {\phi_{k 2}} \\ {\vdots} \\ {\phi_{k k}}\end{array}\right)=\left( \begin{array}{c}{\rho_{1}} \\ {\rho_{2}} \\ {\vdots} \\ {\rho_{k}}\end{array}\right) 1ρ1ρk1ρ11ρk2ρk1ρk21ϕk1ϕk2ϕkk=ρ1ρ2ρk
    根据线性方程组求解的Gramer法则,有:
    ϕ k k = D k D \phi_{k k}=\frac{D_{k}}{D} ϕkk=DDk其中, D D D为系数矩阵的行列式, D k D_k Dk则是将 D D D中第 k k k列换成等号右边的系数。
    D = ∣ 1 ρ 1 ⋯ ρ k − 1 ρ 1 1 ⋯ ρ k − 2 ⋮ ⋮ ⋮ ρ k − 1 ρ k − 2 ⋯ 1 ∣ D=\left| \begin{array}{cccc}{1} &amp; {\rho_{1}} &amp; {\cdots} &amp; {\rho_{k-1}} \\ {\rho_{1}} &amp; {1} &amp; {\cdots} &amp; {\rho_{k-2}} \\ {\vdots} &amp; {\vdots} &amp; { } &amp; {\vdots} \\ {\rho_{k-1}} &amp; {\rho_{k-2}} &amp; {\cdots} &amp; {1}\end{array}\right| D=1ρ1ρk1ρ11ρk2ρk1ρk21
    D k = ∣ 1 ρ 1 ⋯ ρ 1 ρ 1 1 ⋯ ρ 2 ⋮ ⋮ ⋮ ρ k − 1 ρ k − 2 ⋯ ρ k ∣ D_{k}=\left| \begin{array}{cccc}{1} &amp; {\rho_{1}} &amp; {\cdots} &amp; {\rho_{1}} \\ {\rho_{1}} &amp; {1} &amp; {\cdots} &amp; {\rho_{2}} \\ {\vdots} &amp; {\vdots} &amp; {} &amp; {\vdots} \\ {\rho_{k-1}} &amp; {\rho_{k-2}} &amp; {\cdots} &amp; {\rho_{k}}\end{array}\right| Dk=1ρ1ρk1ρ11ρk2ρ1ρ2ρk
    对于 A R ( p ) AR(p) AR(p)模型,其偏自相关系数具有一个特性,即 p p p截尾性,即对于 ∀ k &gt; p \forall k&gt;p k>p,有 ϕ k k = 0 \phi_{k k}=0 ϕkk=0。由于篇幅原因,这里证明略。
2.2 MA模型
2.2.1 MA模型的定义

    一般称具有如下结构的模型为 q q q阶移动平均模型,简称 MA ⁡ ( q ) \operatorname{MA}(q) MA(q)
{ x t = μ + ε t − θ 1 ε t − 1 − θ 2 ε t − 2 − ⋯ − θ q ε t − q θ q ≠ 0 E ( ε t ) = 0 , Var ⁡ ( ε t ) = σ ε 2 , E ( ε t ε s ) = 0 , s ≠ t \left\{\begin{array}{l}{x_{t}=\mu+\varepsilon_{t}-\theta_{1} \varepsilon_{t-1}-\theta_{2} \varepsilon_{t-2}-\cdots-\theta_{q} \varepsilon_{t-q}} \\ {\theta_{q} \neq 0} \\ {E\left(\varepsilon_{t}\right)=0, \operatorname{Var}\left(\varepsilon_{t}\right)=\sigma_{\varepsilon}^{2}, E\left(\varepsilon_{t} \varepsilon_{s}\right)=0, s \neq t}\end{array}\right. xt=μ+εtθ1εt1θ2εt2θqεtqθq̸=0E(εt)=0,Var(εt)=σε2,E(εtεs)=0,s̸=t
μ = 0 \mu=0 μ=0,则上式称为中心化 M A ( q ) \mathrm{MA}(q) MA(q)模型,采用延迟算子表示,可以表示为:
x t = Θ ( B ) ε t x_{t}=\Theta(B) \varepsilon_{t} xt=Θ(B)εt
其中, Θ ( B ) = 1 − θ 1 B − θ 2 B 2 − ⋯ − θ q B q \Theta(B)=1-\theta_{1} B-\theta_{2} B^{2}-\cdots-\theta_{q} B^{q} Θ(B)=1θ1Bθ2B2θqBq,称为 q q q阶移动平均系数多项式。

2.2.2 MA模型的统计性质
  1. 常数均值
    q &lt; ∞ q&lt;\infty q<时, M A ( q ) \mathrm{MA}(q) MA(q)的均值为:
    E x t = E ( μ + ε t − θ 1 ε t − 1 − θ 2 ε t − 2 − ⋯ − θ q ε t − q ) = μ E x_{t}=E\left(\mu+\varepsilon_{t}-\theta_{1} \varepsilon_{t-1}-\theta_{2} \varepsilon_{t-2}-\cdots-\theta_{q} \varepsilon_{t-q}\right)=\mu Ext=E(μ+εtθ1εt1θ2εt2θqεtq)=μ
    特别地,对于中心化 M A ( q ) \mathrm{MA}(q) MA(q)模型,其均值为0。
  2. 常数方差
    Var ⁡ ( x t ) = Var ⁡ ( μ + ε t − θ 1 ε t − 1 − θ 2 ε t − 2 − ⋯ − θ q ε t − q ) = ( 1 + θ 1 2 + ⋯ + θ q 2 ) σ ε 2 \operatorname{Var}\left(x_{t}\right)=\operatorname{Var}\left(\mu+\varepsilon_{t}-\theta_{1} \varepsilon_{t-1}-\theta_{2} \varepsilon_{t-2}-\cdots-\theta_{q} \varepsilon_{t-q}\right)=\left(1+\theta_{1}^{2}+\cdots+\theta_{q}^{2}\right) \sigma_{\varepsilon}^{2} Var(xt)=Var(μ+εtθ1εt1θ2εt2θqεtq)=(1+θ12++θq2)σε2
  3. 自协方差函数
    γ k = E ( x t x t − k ) = E [ ( ε t − θ 1 ε t − 1 − ⋯ − θ q ε t − q ) ( ε t − k − θ 1 ε t − k − 1 − ⋯ − θ q ε t − k − q ) ] = { ( 1 + θ 1 2 + ⋯ + θ q 2 ) σ ε 2 , k = 0 ( − θ k + ∑ i = 1 q − k θ i θ k + i ) σ ε 2 , 1 ⩽ k ⩽ q 0 , k &gt; q \begin{aligned} \gamma_{k} &amp;=E\left(x_{t} x_{t-k}\right) \\ &amp;=E\left[\left(\varepsilon_{t}-\theta_{1} \varepsilon_{t-1}-\cdots-\theta_{q} \varepsilon_{t-q}\right)\left(\varepsilon_{t-k}-\theta_{1} \varepsilon_{t-k-1}-\cdots-\theta_{q} \varepsilon_{t-k-q}\right)\right] \\ &amp;=\left\{\begin{array}{ll}{\left(1+\theta_{1}^{2}+\cdots+\theta_{q}^{2}\right) \sigma_{\varepsilon}^{2},} &amp; {k=0} \\ {\left(-\theta_{k}+\sum_{i=1}^{q-k} \theta_{i} \theta_{k+i}\right) \sigma_{\varepsilon}^{2},} &amp; {1 \leqslant k \leqslant q} \\ {0,} &amp; {k&gt;q}\end{array}\right.\end{aligned} γk=E(xtxtk)=E[(εtθ1εt1θqεtq)(εtkθ1εtk1θqεtkq)]=(1+θ12++θq2)σε2,(θk+i=1qkθiθk+i)σε2,0,k=01kqk>q
  4. 自相关系数q阶截尾
    ρ k = γ k γ 0 = { 1 , k = 0 − θ k + ∑ i = 1 q − k θ i θ k + i 1 + θ 1 2 + ⋯ + θ q 2 , 1 ⩽ k ⩽ q 0 , k &gt; q \rho_{k}=\frac{\gamma_{k}}{\gamma_{0}}=\left\{\begin{array}{ll}{1,} &amp; {k=0} \\ {\frac{ -\theta_{k}+\sum_{i=1}^{q-k} \theta_{i} \theta_{k+i}}{1+\theta_{1}^{2}+\cdots+\theta_{q}^{2}},} &amp; {1 \leqslant k \leqslant q} \\ {0,} &amp; {k&gt;q}\end{array}\right. ρk=γ0γk=1,1+θ12++θq2θk+i=1qkθiθk+i,0,k=01kqk>q
  5. 偏自相关系数拖尾
    由于不同的MA模型可能会有对应同一个自相关系数,因此,为了保证一个给定的自相关系数可以唯一对应一个MA模型,需要对模型添加一个约束条件,这个约束条件称为可逆性条件,即对于 M A ( q ) \mathrm{MA}(q) MA(q)模型,将其表达成 A R ( p ) \mathrm{AR}(p) AR(p)的形式:
    ε t = x t Θ ( B ) \varepsilon_{t}=\frac{x_{t}}{\Theta(B)} εt=Θ(B)xt
    其中, Θ ( B ) = 1 − θ 1 B − ⋯ − θ q B q \Theta(B)=1-\theta_{1} B-\cdots-\theta_{q} B^{q} Θ(B)=1θ1BθqBq称为移动平均系数多项式,假定 1 λ 1 , ⋯ &ThinSpace; , 1 λ q \frac{1}{\lambda_{1}}, \cdots, \frac{1}{\lambda_{q}} λ11,,λq1为该系数多项式的 q q q个根,则 Θ ( B ) \Theta(B) Θ(B)可以分解为:
    Θ ( B ) = ∏ k = 1 q ( 1 − λ k B ) \Theta(B)=\prod_{k=1}^{q}\left(1-\lambda_{k} B\right) Θ(B)=k=1q(1λkB)
    因此,可以得:
    ε t = x t ( 1 − λ 1 B ) ⋯ ( 1 − λ q B ) \varepsilon_{t}=\frac{x_{t}}{\left(1-\lambda_{1} B\right) \cdots\left(1-\lambda_{q} B\right)} εt=(1λ1B)(1λqB)xt
    要使得MA模型可逆,则要求 ∣ λ i ∣ &lt; 1 \left|\lambda_{i}\right|&lt;1 λi<1,即 M A ( q ) \mathrm{MA}(q) MA(q)模型的系数多项式的根都在单位圆外,这个条件记为可逆性条件。
    这样一来,如果一个 M A ( q ) \mathrm{MA}(q) MA(q)模型满足可逆性条件,那么,它就可以有以下两种等价表达方式:
    { Θ ( B ) ε t = x t ε t = I ( B ) x t \left\{\begin{array}{l}{\Theta(B) \varepsilon_{t}=x_{t}} \\ {\varepsilon_{t}=I(B) x_{t}}\end{array}\right. {Θ(B)εt=xtεt=I(B)xt
    有上式可得:
    Θ ( B ) I ( B ) x t = x t \Theta(B) I(B) x_{t}=x_{t} Θ(B)I(B)xt=xt
    展开得:
    ( 1 − ∑ k = 1 q θ k B k ) ( 1 + ∑ j = 0 ∞ I j B j ) 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} (1k=1qθkBk)(1+j=0IjBj)xt=xt
    其中,
    { I 0 = 1 I j = ∑ k = 1 j θ k ′ I j − k , j ⩾ 1 \left\{\begin{array}{l}{I_{0}=1} \\ {I_{j}=\sum_{k=1}^{j} \theta_{k}^{\prime} I_{j-k}, \quad j \geqslant 1}\end{array}\right. {I0=1Ij=k=1jθkIjk,j1
    θ k ′ = { θ k , k ⩽ q 0 , k &gt; q \theta_{k}^{\prime}=\left\{\begin{array}{ll}{\theta_{k},} &amp; {k \leqslant q} \\ {0,} &amp; {k&gt;q}\end{array}\right. θk={θk,0,kqk>q
    因此,对于一个可逆 M A ( q ) \mathrm{MA}(q) MA(q)模型,可以将其等价写成 A R ( ∞ ) \mathrm{AR}(\infty) AR()形式:
    I ( B ) x t = ε t I(B) x_{t}=\varepsilon_{t} I(B)xt=εt
    而由前面我们知道, A R ( p ) \mathrm{AR}(p) AR(p)模型具有 p p p阶截尾,因此, M A ( q ) \mathrm{MA}(q) MA(q)模型的偏自相关系数 ∞ \infty 截尾,即具有拖尾性
2.3 ARMA模型
2.3.1 ARMA模型的定义

    把具有如下结构的模型称为自回归移动平均模型,简记为 A R M A ( p , q ) \mathrm{ARMA}(p, q) ARMA(p,q)
{ x t = ϕ 0 + ϕ 1 x t − 1 + ⋯ + ϕ p x t − p + ε t − θ 1 ε t − 1 − ⋯ − θ q ε t − q ϕ p ≠ 0 , θ q ≠ 0 E ( ε t ) = 0 , Var ⁡ ( ε t ) = σ ε 2 , E ( ε t ε s ) = 0 , s ≠ t E ( x s ε t ) = 0 , ∀ s &lt; t \left\{\begin{array}{l}{x_{t}=\phi_{0}+\phi_{1} x_{t-1}+\cdots+\phi_{p} x_{t-p}+\varepsilon_{t}-\theta_{1} \varepsilon_{t-1}-\cdots-\theta_{q} \varepsilon_{t-q}} \\ {\phi_{p} \neq 0, \theta_{q} \neq 0} \\ {E\left(\varepsilon_{t}\right)=0, \operatorname{Var}\left(\varepsilon_{t}\right)=\sigma_{\varepsilon}^{2}, E\left(\varepsilon_{t} \varepsilon_{s}\right)=0, s \neq t} \\ {E\left(x_{s} \varepsilon_{t}\right)=0, \forall s&lt;t}\end{array}\right. xt=ϕ0+ϕ1xt1++ϕpxtp+εtθ1εt1θqεtqϕp̸=0,θq̸=0E(εt)=0,Var(εt)=σε2,E(εtεs)=0,s̸=tE(xsεt)=0,s<t
ϕ 0 = 0 \phi_{0}=0 ϕ0=0时,即为中心化 A R M A ( p , q ) \mathrm{ARMA}(p, q) ARMA(p,q)模型,用延迟算子可以表示为:
Φ ( B ) x t = Θ ( B ) ε t \Phi(B) x_{t}=\Theta(B) \varepsilon_{t} Φ(B)xt=Θ(B)εt
其中, Φ ( B ) = 1 − ϕ 1 B − ⋯ − ϕ p B p \Phi(B)=1-\phi_{1} B-\cdots-\phi_{p} B^{p} Φ(B)=1ϕ1BϕpBp p p p阶自回归系数多项式, Θ ( B ) = 1 − θ 1 B − ⋯ − θ q B q \Theta(B)=1-\theta_{1} B-\cdots-\theta_{q} B^{q} Θ(B)=1θ1BθqBq q q q阶移动平均系数多项式。当 q = 0 q=0 q=0 p = 0 p=0 p=0时,则 A R M A ( p , q ) \mathrm{ARMA}(p, q) ARMA(p,q)则相当于 A R ( p ) AR(p) AR(p) M A ( q ) MA(q) MA(q)模型。

2.3.2 ARMA模型的统计性质
  1. 均值
    对模型两边同取均值得:
    E x t = ϕ 0 1 − ϕ 1 − ⋯ − ϕ p E x_{t}=\frac{\phi_{0}}{1-\phi_{1}-\cdots-\phi_{p}} Ext=1ϕ1ϕpϕ0
  2. 自协方差函数
    γ ( k ) = E ( x t x t + k ) = E [ ( ∑ i = 0 ∞ G i ε t − i ) ( ∑ j = 0 ∞ G j ε t + k − j ) ] = E [ ∑ i = 0 ∞ G i ∑ j = 0 ∞ G j ε t − i ε t + k − j ] = σ ε 2 ∑ i = 0 G i G i + k \begin{aligned} \gamma(k) &amp;=E\left(x_{t} x_{t+k}\right) \\ &amp;=E\left[\left(\sum_{i=0}^{\infty} G_{i} \varepsilon_{t-i}\right)\left(\sum_{j=0}^{\infty} G_{j} \varepsilon_{t+k-j}\right)\right] \\ &amp;=E\left[\sum_{i=0}^{\infty} G_{i} \sum_{j=0}^{\infty} G_{j} \varepsilon_{t-i} \varepsilon_{t+k-j}\right] \\ &amp;=\sigma_{\varepsilon}^{2} \sum_{i=0} G_{i} G_{i+k} \end{aligned} γ(k)=E(xtxt+k)=E[(i=0Giεti)(j=0Gjεt+kj)]=E[i=0Gij=0Gjεtiεt+kj]=σε2i=0GiGi+k
  3. 自相关系数
    ρ ( k ) = γ ( k ) γ ( 0 ) = ∑ j = 0 ∞ G j G j + k ∑ j = 0 ∞ G j 2 \rho(k)=\frac{\gamma(k)}{\gamma(0)}=\frac{\sum_{j=0}^{\infty} G_{j} G_{j+k}}{\sum_{j=0}^{\infty} G_{j}^{2}} ρ(k)=γ(0)γ(k)=j=0Gj2j=0GjGj+k
    根据 ARMA ⁡ ( p , q ) \operatorname{ARMA}(p, q) ARMA(p,q)模型的自相关系数计算公式可以知道, ARMA ⁡ ( p , q ) \operatorname{ARMA}(p, q) ARMA(p,q)的自相关系数具有拖尾性, 并且由于 ARMA ⁡ ( p , q ) \operatorname{ARMA}(p, q) 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,因此,我们可以认为该时间序列不是白噪声序列。接下来,对于该时间序列,我们采用 A R M A ( p , q ) ARMA(p,q) ARMA(p,q)进行拟合,但是需要确定 p p p q q 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 p p定为2,而 q q q定为0,即采用 A R ( 2 ) AR(2) AR(2)模型进行拟合。不过,通过自相关图和偏自相关图有时比较难做出判断,此时,也可以通过模型的BIC值来确定 p p p q q 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

     定阶后,我们根据阶数拟合 A R M A ARMA 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 = ε t 1 − 0.7185 B + 0.5294 B 2 x_{t}-11.0227=\frac{\varepsilon_{t}}{1-0.7185 B+0.5294 B^{2}} xt11.0227=10.7185B+0.5294B2εt

4. 总结

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值