Yule-Walker方程 完整推导过程以及python代码复现

Every cell phone call solves the Yule-Walker equations every ten microseconds. ——Thierry Dutoit

介绍

代码、数据全部免费,都放在我的gitee仓库里面https://gitee.com/yuanzhoulvpi/time_series,想要使用的可以直接到我这个仓库下载。

本文是继python与时间序列(开篇)的第3篇。在找pacf如何计算的时候,发现现在计算pacf都在用Yule-Walker方法,包括在计算AR§的时候也是这个方程。(不是不用别的方法了)。

我在这里把Yule-Walker方法的公式推导写出来,并且把我写的python代码也整理出来。方便参考。

最新的代码在gitee仓库里面:https://gitee.com/yuanzhoulvpi/time_series/blob/master/%E6%A1%88%E4%BE%8B/03python%E5%A4%8D%E7%8E%B0Yule-Walker%E6%96%B9%E7%A8%8B.ipynb

使用Yule-Walker方法计算AR§

计算一个均值为0的离散随机时间序列 { X i } N \{X_i\} ^ N {Xi}N的AR§参数,怎么做?
我们要估计的也就是下面公式中的 ξ i + 1 \xi_{i+1} ξi+1

x i + 1 = ϕ 1 x i + ϕ 2 x i − 1 + ⋯ + ϕ p x i − p + 1 + ξ i + 1 (1) x_{i+1}=\phi_{1} x_{i}+\phi_{2} x_{i-1}+\cdots+\phi_{p} x_{i-p+1}+\xi_{i+1} \tag 1 xi+1=ϕ1xi+ϕ2xi1++ϕpxip+1+ξi+1(1)

1.转置

1.1 p=1

p=1的时候,公式变成这样 X i + 1 = ϕ 1 x i + ϵ i + 1 X_{i+1} = \phi_1 x_i + \epsilon_{i+1} Xi+1=ϕ1xi+ϵi+1
上面的形式可以进一步写成这样:
( x 2 x 3 ⋮ x N ) ⏟ b = ( x 1 x 2 ⋮ x N − 1 ) ⏟ A ϕ 1 \underbrace{\left(\begin{array}{c} x_{2} \\ x_{3} \\ \vdots \\ x_{N} \end{array}\right)}_{\mathbf{b}}=\underbrace{\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{N-1} \end{array}\right)}_{\mathbf{A}} \phi_{1} b x2x3xN=A x1x2xN1ϕ1

可以使用最小二乘法求解:
ϕ ^ 1 = ( A T A ) − 1 A T b = ∑ i = 1 N − 1 x i x i + 1 ∑ i = 1 N − 1 x i 2 = c 1 c o = r 1 \hat{\phi}_{1}=\left(\mathbf{A}^{T} \mathbf{A}\right)^{-1} \mathbf{A}^{T} \mathbf{b}=\frac{\sum_{i=1}^{N-1} x_{i} x_{i+1}}{\sum_{i=1}^{N-1} x_{i}^{2}}=\frac{c_{1}}{c_{o}}=r_{1} ϕ^1=(ATA)1ATb=i=1N1xi2i=1N1xixi+1=coc1=r1
上面的 c i c_{i} ci r i r_{i} ri分别是第i个自协方差和自协方差系数。

1.2 p=2

p=2的时候,公式变成这样 X i + 1 = ϕ 1 x i + ϕ 2 x i − 1 + ϵ i + 1 X_{i+1} = \phi_1 x_i +\phi_2 x_{i-1} + \epsilon_{i+1} Xi+1=ϕ1xi+ϕ2xi1+ϵi+1
上面的形式进一步写成这样:

( x 3 x 4 ⋮ x N ) ⏟ b = ( x 2 x 1 x 3 x 2 ⋮ ⋮ x N − 1 x N − 2 ) ⏟ A ( ϕ 1 ϕ 2 ) ⏟ Φ . \underbrace{\left(\begin{array}{c} x_{3} \\ x_{4} \\ \vdots \\ x_{N} \end{array}\right)}_{\mathbf{b}}=\underbrace{\left(\begin{array}{cc} x_{2} & x_{1} \\ x_{3} & x_{2} \\ \vdots & \vdots \\ x_{N-1} & x_{N-2} \end{array}\right)}_{\mathbf{A}} \underbrace{\left(\begin{array}{c} \phi_{1} \\ \phi_{2} \end{array}\right)}_{\Phi} . b x3x4xN=A x2x3xN1x1x2xN2Φ (ϕ1ϕ2).

和第一个不一样,这个时候上面的等式写成这样:
Φ ^ = ( A T A ) − 1 A T b \hat{\boldsymbol{\Phi}}=\left(\mathbf{A}^{T} \mathbf{A}\right)^{-1} \mathbf{A}^{T} \mathbf{b} Φ^=(ATA)1ATb
是这么计算的:
( A T A ) − 1 = [ ( x 2 x 3 ⋯ x N − 1 x 1 x 2 ⋯ x N − 2 ) ( x 2 x 1 x 3 x 2 x N − 1 x N − 2 ) ] − 1 = ( ∑ i = 2 N − 1 x i 2 ∑ i = 2 N − 1 x i x i − 1 ∑ i = 2 N − 1 x i x i − 1 ∑ i = 1 N − 2 x i 2 ) − 1 = 1 ∑ i = 2 N − 1 x i 2 ∑ i = 1 N − 2 x i 2 − ∑ i = 2 N − 1 x i x i − 1 ∑ i = 2 N − 1 x i x i − 1 ( ∑ i = 1 N − 2 x i 2 − ∑ i = 2 N − 1 x i x i − 1 − ∑ i = 2 N − 1 x i x i − 1 ∑ i = 2 N − 1 x i 2 ) \begin{gathered} \left(\mathbf{A}^{T} \mathbf{A}\right)^{-1}=\left[\left(\begin{array}{llll} x_{2} & x_{3} & \cdots & x_{N-1} \\ x_{1} & x_{2} & \cdots & x_{N-2} \end{array}\right)\left(\begin{array}{cc} x_{2} & x_{1} \\ x_{3} & x_{2} \\ x_{N-1} & x_{N-2} \end{array}\right)\right]^{-1} \\ =\left(\begin{array}{cc} \sum_{i=2}^{N-1} x_{i}^{2} & \sum_{i=2}^{N-1} x_{i} x_{i-1} \\ \sum_{i=2}^{N-1} x_{i} x_{i-1} & \sum_{i=1}^{N-2} x_{i}^{2} \end{array}\right)^{-1} \\ =\frac{1}{\sum_{i=2}^{N-1} x_{i}^{2} \sum_{i=1}^{N-2} x_{i}^{2}-\sum_{i=2}^{N-1} x_{i} x_{i-1} \sum_{i=2}^{N-1} x_{i} x_{i-1}}\left(\begin{array}{cc} \sum_{i=1}^{N-2} x_{i}^{2} & -\sum_{i=2}^{N-1} x_{i} x_{i-1} \\ -\sum_{i=2}^{N-1} x_{i} x_{i-1} & \sum_{i=2}^{N-1} x_{i}^{2} \end{array}\right) \end{gathered} (ATA)1=(x2x1x3x2xN1xN2)x2x3xN1x1x2xN21=(i=2N1xi2i=2N1xixi1i=2N1xixi1i=1N2xi2)1=i=2N1xi2i=1N2xi2i=2N1xixi1i=2N1xixi11(i=1N2xi2i=2N1xixi1i=2N1xixi1i=2N1xi2)

接下来,假设时间序列是平稳的,所以说自协方差只和滞后多少项相关。利用这个性质,在这个案例里面可以得到这样的等式:

( A T A ) − 1 = 1 c o 2 − c 1 2 ( c o − c 1 − c 1 c o ) ( A T A ) − 1 = 1 c o 2 ( 1 − r 1 2 ) ( c o − c 1 − c 1 c o ) ( A T A ) − 1 = 1 c o ( 1 − r 1 2 ) ( r o − r 1 − r 1 r o ) \begin{gathered} \left(\mathbf{A}^{T} \mathbf{A}\right)^{-1}=\frac{1}{c_{o}^{2}-c_{1}^{2}}\left(\begin{array}{rr} c_{o} & -c_{1} \\ -c_{1} & c_{o} \end{array}\right) \\ \left(\mathbf{A}^{T} \mathbf{A}\right)^{-1}=\frac{1}{c_{o}^{2}\left(1-r_{1}^{2}\right)}\left(\begin{array}{rr} c_{o} & -c_{1} \\ -c_{1} & c_{o} \end{array}\right) \\ \left(\mathbf{A}^{T} \mathbf{A}\right)^{-1}=\frac{1}{c_{o}\left(1-r_{1}^{2}\right)}\left(\begin{array}{rr} r_{o} & -r_{1} \\ -r_{1} & r_{o} \end{array}\right) \end{gathered} (ATA)1=co2c121(coc1c1co)(ATA)1=co2(1r12)1(coc1c1co)(ATA)1=co(1r12)1(ror1r1ro)
而且:
A T b = ( x 2 x 3 ⋯ x N − 1 x 1 x 2 ⋯ x N − 2 ) ( x 3 x 4 ⋮ x N ) = ( ∑ i = 3 N x i x i − 1 ∑ i = 3 N x i x i − 2 , ) \mathbf{A}^{T} \mathbf{b}=\left(\begin{array}{llll} x_{2} & x_{3} & \cdots & x_{N-1} \\ x_{1} & x_{2} & \cdots & x_{N-2} \end{array}\right)\left(\begin{array}{l} x_{3} \\ x_{4} \\ \vdots \\ x_{N} \end{array}\right)=\left(\begin{array}{l} \sum_{i=3}^{N} x_{i} x_{i-1} \\ \sum_{i=3}^{N} x_{i} x_{i-2}, \end{array}\right) ATb=(x2x1x3x2xN1xN2)x3x4xN=(i=3Nxixi1i=3Nxixi2,)

又因为这个时间序列是平稳的,所以:

A T b = ( c 1 c 2 ) \mathbf{A}^{T} \mathbf{b}=\left(\begin{array}{l} c_{1} \\ c_{2} \end{array}\right) ATb=(c1c2)

结合上面的公式,可以有:
( A T A ) − 1 A T b = 1 c o ( 1 − r 1 2 ) ( r o − r 1 − r 1 r o ) ( c 1 c 2 ) = 1 1 − r 1 2 ( 1 − r 1 − r 1 1 ) ( r 1 r 2 ) . \begin{gathered} \left(\mathbf{A}^{T} \mathbf{A}\right)^{-1} \mathbf{A}^{T} \mathbf{b}=\frac{1}{c_{o}\left(1-r_{1}^{2}\right)}\left(\begin{array}{rr} r_{o} & -r_{1} \\ -r_{1} & r_{o} \end{array}\right)\left(\begin{array}{l} c_{1} \\ c_{2} \end{array}\right) \\ =\frac{1}{1-r_{1}^{2}}\left(\begin{array}{rr} 1 & -r_{1} \\ -r_{1} & 1 \end{array}\right)\left(\begin{array}{l} r_{1} \\ r_{2} \end{array}\right) . \end{gathered} (ATA)1ATb=co(1r12)1(ror1r1ro)(c1c2)=1r121(1r1r11)(r1r2).
将上面的矩阵分为两个部分,可以得到:
ϕ ^ 1 = r 1 ( 1 − r 2 ) 1 − r 1 2 \hat{\phi}_{1}=\frac{r_{1}\left(1-r_{2}\right)}{1-r_{1}^{2}} ϕ^1=1r12r1(1r2)

以及
ϕ ^ 2 = r 2 − r 1 2 1 − r 1 2 \hat{\phi}_{2}=\frac{r_{2}-r_{1}^{2}}{1-r_{1}^{2}} ϕ^2=1r12r2r12

虽然可以按照p=2继续计算p=3的情况,但是那样的话,代数计算会变得很复杂。
幸运的是,有一种简单的方法来计算AR§的系数,这个方法就叫Yule-Waler公式。

2. Yule-Waler公式

这是一个AR§公式:
x i + 1 = ϕ 1 x i + ϕ 2 x i − 1 + ⋯ + ϕ p x i − p + 1 + ξ i + 1 x_{i+1}=\phi_{1} x_{i}+\phi_{2} x_{i-1}+\cdots+\phi_{p} x_{i-p+1}+\xi_{i+1} xi+1=ϕ1xi+ϕ2xi1++ϕpxip+1+ξi+1

2.1 lag=1(滞后1)

在左右两边乘上 x i x_{i} xi,得
x i x i + 1 = ∑ j = 1 p ( ϕ j x i x i − j + 1 ) + x i ξ i + 1 x_{i} x_{i+1}=\sum_{j=1}^{p}\left(\phi_{j} x_{i} x_{i-j+1}\right)+x_{i} \xi_{i+1} xixi+1=j=1p(ϕjxixij+1)+xiξi+1

i和j是各自的时间索引。把 { ϕ j } \{\phi_j\} {ϕj}都拎出来, x j x_j xj都写到一起,得到这样的公式:
⟨ x i x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i x i − j + 1 ⟩ ) + ⟨ x i ξ i + 1 ⟩ \left\langle x_{i} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i} x_{i-j+1}\right\rangle\right)+\left\langle x_{i} \xi_{i+1}\right\rangle xixi+1=j=1p(ϕjxixij+1)+xiξi+1

注意到 ⟨ x i ξ i + 1 ⟩ = 0 \left\langle x_{i} \xi_{i+1}\right\rangle = 0 xiξi+1=0 因为截距 ξ \xi ξ和当前时间是无关的,因此:
⟨ x i x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i x i − j + 1 ⟩ ) \left\langle x_{i} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i} x_{i-j+1}\right\rangle\right) xixi+1=j=1p(ϕjxixij+1)

再除以N-1,再利用自协方差的平衡性: c − l = c l c_{-l} = c_{l} cl=cl,得:
c 1 = ∑ j = 1 p ϕ j c j − 1 c_{1}=\sum_{j=1}^{p} \phi_{j} c_{j-1} c1=j=1pϕjcj1

再除以 c 0 c_0 c0,得:
r 1 = ∑ j = 1 p ϕ j r j − 1 r_{1}=\sum_{j=1}^{p} \phi_{j} r_{j-1} r1=j=1pϕjrj1

2.2 lag=2(滞后2)

两边乘以 x i − 1 x_{i-1} xi1, 得到:

x i − 1 x i + 1 = ∑ j = 1 p ( ϕ j x i − 1 x i − j + 1 ) + x i − 1 ξ i + 1 x_{i-1} x_{i+1}=\sum_{j=1}^{p}\left(\phi_{j} x_{i-1} x_{i-j+1}\right)+x_{i-1} \xi_{i+1} xi1xi+1=j=1p(ϕjxi1xij+1)+xi1ξi+1

然后:
⟨ x i − 1 x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i − 1 x i − j + 1 ⟩ ) + ⟨ x i − 1 ξ i + 1 ⟩ \left\langle x_{i-1} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i-1} x_{i-j+1}\right\rangle\right)+\left\langle x_{i-1} \xi_{i+1}\right\rangle xi1xi+1=j=1p(ϕjxi1xij+1)+xi1ξi+1

然后:
⟨ x i − 1 x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i − 1 x i − j + 1 ⟩ ) \left\langle x_{i-1} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i-1} x_{i-j+1}\right\rangle\right) xi1xi+1=j=1p(ϕjxi1xij+1)

然后:
c 2 = ∑ j = 1 p ϕ j c j − 2 c_{2}=\sum_{j=1}^{p} \phi_{j} c_{j-2} c2=j=1pϕjcj2

然后:
r 2 = ∑ j = 1 p ϕ j r j − 2 r_{2}=\sum_{j=1}^{p} \phi_{j} r_{j-2} r2=j=1pϕjrj2

2.3 lag=k(滞后k)

两边乘以 x i − k − 1 x_{i-k-1} xik1, 得到:

x i − k + 1 x i + 1 = ∑ j = 1 p ( ϕ j x i − k + 1 x i − j + 1 ) + x i − k + 1 ξ i + 1 x_{i-k+1} x_{i+1}=\sum_{j=1}^{p}\left(\phi_{j} x_{i-k+1} x_{i-j+1}\right)+x_{i-k+1} \xi_{i+1} xik+1xi+1=j=1p(ϕjxik+1xij+1)+xik+1ξi+1

然后:
⟨ x i − k + 1 x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i − k + 1 x i − j + 1 ⟩ ) + ⟨ x i − k + 1 ξ i + 1 ⟩ \left\langle x_{i-k+1} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i-k+1} x_{i-j+1}\right\rangle\right)+\left\langle x_{i-k+1} \xi_{i+1}\right\rangle xik+1xi+1=j=1p(ϕjxik+1xij+1)+xik+1ξi+1

然后:
⟨ x i − k + 1 x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i − k + 1 x i − j + 1 ⟩ ) \left\langle x_{i-k+1} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i-k+1} x_{i-j+1}\right\rangle\right) xik+1xi+1=j=1p(ϕjxik+1xij+1)

然后:
c k = ∑ j = 1 p ϕ j c j − k c_{k}=\sum_{j=1}^{p} \phi_{j} c_{j-k} ck=j=1pϕjcjk

然后:
r k = ∑ j = 1 p ϕ j r j − k r_{k}=\sum_{j=1}^{p} \phi_{j} r_{j-k} rk=j=1pϕjrjk

2.4 lag=p(滞后p)

两边乘以 x i − p − 1 x_{i-p-1} xip1, 得到:

x i − p + 1 x i + 1 = ∑ j = 1 p ( ϕ j x i − p + 1 x i − j + 1 ) + x i − p + 1 ξ i + 1 x_{i-p+1} x_{i+1}=\sum_{j=1}^{p}\left(\phi_{j} x_{i-p+1} x_{i-j+1}\right)+x_{i-p+1} \xi_{i+1} xip+1xi+1=j=1p(ϕjxip+1xij+1)+xip+1ξi+1

然后:
⟨ x i − p + 1 x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i − p + 1 x i − j + 1 ⟩ ) + ⟨ x i − p + 1 ξ i + 1 ⟩ \left\langle x_{i-p+1} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i-p+1} x_{i-j+1}\right\rangle\right)+\left\langle x_{i-p+1} \xi_{i+1}\right\rangle xip+1xi+1=j=1p(ϕjxip+1xij+1)+xip+1ξi+1

然后:
⟨ x i − p + 1 x i + 1 ⟩ = ∑ j = 1 p ( ϕ j ⟨ x i − p + 1 x i − j + 1 ⟩ ) \left\langle x_{i-p+1} x_{i+1}\right\rangle=\sum_{j=1}^{p}\left(\phi_{j}\left\langle x_{i-p+1} x_{i-j+1}\right\rangle\right) xip+1xi+1=j=1p(ϕjxip+1xij+1)

然后:
c p = ∑ j = 1 p ϕ j c j − p c_{p}=\sum_{j=1}^{p} \phi_{j} c_{j-p} cp=j=1pϕjcjp

然后:
r p = ∑ j = 1 p ϕ j r j − p r_{p}=\sum_{j=1}^{p} \phi_{j} r_{j-p} rp=j=1pϕjrjp

2.5 将上面的公式放一起

有:
r 1 = ϕ 1 r o + ϕ 2 r 1 + ϕ 3 r 2 + ⋯ + ϕ p − 1 r p − 2 + ϕ p r p − 1 r 2 = ϕ 1 r 1 + ϕ 2 r o + ϕ 3 r 1 + ⋯ + ϕ p − 1 r p − 3 + ϕ p r p − 2 ⋮ r p − 1 = ϕ 1 r p − 2 + ϕ 2 r p − 3 + ϕ 3 r p − 4 + ⋯ + ϕ p − 1 r o + ϕ p r 1 r p = ϕ 1 r p − 1 + ϕ 2 r p − 2 + ϕ 3 r p − 3 + ⋯ + ϕ p − 1 r 1 + ϕ p r o \begin{gathered} r_{1}=\phi_{1} r_{o}+\phi_{2} r_{1}+\phi_{3} r_{2}+\cdots+\phi_{p-1} r_{p-2}+\phi_{p} r_{p-1} \\ r_{2}=\phi_{1} r_{1}+\phi_{2} r_{o}+\phi_{3} r_{1}+\cdots+\phi_{p-1} r_{p-3}+\phi_{p} r_{p-2} \\ \vdots \\ r_{p-1}=\phi_{1} r_{p-2}+\phi_{2} r_{p-3}+\phi_{3} r_{p-4}+\cdots+\phi_{p-1} r_{o}+\phi_{p} r_{1} \\ r_{p}=\phi_{1} r_{p-1}+\phi_{2} r_{p-2}+\phi_{3} r_{p-3}+\cdots+\phi_{p-1} r_{1}+\phi_{p} r_{o} \end{gathered} r1=ϕ1ro+ϕ2r1+ϕ3r2++ϕp1rp2+ϕprp1r2=ϕ1r1+ϕ2ro+ϕ3r1++ϕp1rp3+ϕprp2rp1=ϕ1rp2+ϕ2rp3+ϕ3rp4++ϕp1ro+ϕpr1rp=ϕ1rp1+ϕ2rp2+ϕ3rp3++ϕp1r1+ϕpro
可以被写成:
( r 1 r 2 ⋮ r p − 1 r p ) = ( r o r 1 r 2 ⋯ r p − 2 r p − 1 r 1 r o r 1 ⋯ r p − 3 r p − 2 ⋮ ⋮ r p − 2 r p − 3 r p − 4 ⋯ r o r 1 r p − 1 r p − 2 r p − 3 ⋯ r 1 r o ) ( ϕ 1 ϕ 2 ⋮ ϕ p − 1 ϕ p ) \left(\begin{array}{c} r_{1} \\ r_{2} \\ \vdots \\ r_{p-1} \\ r_{p} \end{array}\right)=\left(\begin{array}{cccccc} r_{o} & r_{1} & r_{2} & \cdots & r_{p-2} & r_{p-1} \\ r_{1} & r_{o} & r_{1} & \cdots & r_{p-3} & r_{p-2} \\ & \vdots & & & \vdots & \\ r_{p-2} & r_{p-3} & r_{p-4} & \cdots & r_{o} & r_{1} \\ r_{p-1} & r_{p-2} & r_{p-3} & \cdots & r_{1} & r_{o} \end{array}\right)\left(\begin{array}{c} \phi_{1} \\ \phi_{2} \\ \vdots \\ \phi_{p-1} \\ \phi_{p} \end{array}\right) r1r2rp1rp=ror1rp2rp1r1rorp3rp2r2r1rp4rp3rp2rp3ror1rp1rp2r1roϕ1ϕ2ϕp1ϕp

又因为 r 0 = 1 r_0 = 1 r0=1,所以可以写成:
( r 1 r 2 ⋮ r p − 1 r p ) ⏟ r ( 1 r 1 r 2 ⋯ r p − 2 r p − 1 r 1 1 r 1 ⋯ r p − 3 r p − 2 ⋮ ⋮ r p − 2 r p − 3 r p − 4 ⋯ 1 r 1 r p − 1 r p − 2 r p − 3 ⋯ r 1 1 ) ⏟ R ( ϕ 1 ϕ 2 ⋮ ϕ p − 1 ϕ p ) ⏟ Φ \underbrace{\left(\begin{array}{c} r_{1} \\ r_{2} \\ \vdots \\ r_{p-1} \\ r_{p} \end{array}\right)}_{\mathbf{r}} \underbrace{\left(\begin{array}{cccccc} 1 & r_{1} & r_{2} & \cdots & r_{p-2} & r_{p-1} \\ r_{1} & 1 & r_{1} & \cdots & r_{p-3} & r_{p-2} \\ & \vdots & & & \vdots & \\ r_{p-2} & r_{p-3} & r_{p-4} & \cdots & 1 & r_{1} \\ r_{p-1} & r_{p-2} & r_{p-3} & \cdots & r_{1} & 1 \end{array}\right)}_{\mathbf{R}} \underbrace{\left(\begin{array}{c} \phi_{1} \\ \phi_{2} \\ \vdots \\ \phi_{p-1} \\ \phi_{p} \end{array}\right)}_{\boldsymbol{\Phi}} r r1r2rp1rpR 1r1rp2rp1r11rp3rp2r2r1rp4rp3rp2rp31r1rp1rp2r11Φ ϕ1ϕ2ϕp1ϕp
整理成:
R Φ = r (2) \mathbf{R} \boldsymbol{\Phi}=\mathbf{r} \tag 2 RΦ=r(2)

最终可以写成这样
Φ ^ = R − 1 r \hat{\boldsymbol{\Phi}}=\mathbf{R}^{-1} \mathbf{r} Φ^=R1r

3. Yule-Walker公式求解过程

  • 循环i, 1 ≤ i ≤ p 1 \leq i \leq p 1ip

    • 计算 R ( i ) \mathbf{R} ^ {(i)} R(i) r ( i ) \mathbf{r} ^ {(i)} r(i)
    • 然后计算 Φ ^ ( i ) \hat{\boldsymbol{\Phi}}^{(i)} Φ^(i),公式为:
      Φ ^ ( i ) = ( R ( i ) ) − 1 r ( i ) = ( ϕ ^ 1 ϕ ^ 2 ⋮ ϕ ^ i ) \hat{\boldsymbol{\Phi}}^{(i)}=\left(\mathbf{R}^{(i)}\right)^{-1} \mathbf{r}^{(i)}=\left(\begin{array}{c} \hat{\phi}_{1} \\ \hat{\phi}_{2} \\ \vdots \\ \hat{\phi}_{i} \end{array}\right) Φ^(i)=(R(i))1r(i)=ϕ^1ϕ^2ϕ^i
    • 只保留 ϕ ^ i \hat{\phi}_{i} ϕ^i,在 1 ≤ j ≤ i − 1 1 \leq j \leq {i-1} 1ji1范围内的 ϕ ^ j \hat{\phi}_{j} ϕ^j都不要
    • 第i个pacf系数这个时候就等于 p a c f ( i ) = ϕ ^ i pacf(i) = \hat{\phi}_{i} pacf(i)=ϕ^i
  • 结束循环i

4. python计算Yule-Walker公式

代码如下:

from scipy.linalg import toeplitz
import numpy as np


def cal_my_yule_walker(x, nlags=5):
    """
    自己实现yule_walker理论
    :param x:
    :param nlags:
    :return:
    """
    x = np.array(x, dtype=np.float64)
    x -= x.mean()
    n = x.shape[0]

    r = np.zeros(shape=nlags+1, dtype=np.float64)
    r[0] = (x ** 2).sum()/n

    for k in range(1, nlags+1):
        r[k] = (x[0:-k] * x[k:]).sum() / (n-k*1)


    R = toeplitz(c=r[:-1])
    result = np.linalg.solve(R, r[1:])
    return result

def cal_my_pacf_yw(x, nlags=5):
    """
    自己通过yule_walker方法求出pacf的值
    :param x:
    :param nlags:
    :return:
    """
    pacf = np.empty(nlags+1) * 0
    pacf[0] = 1.0
    for k in range(1, nlags+1):
        pacf[k] = cal_my_yule_walker(x,nlags=k)[-1]

    return pacf


# 测试函数
data_x = np.random.randint(low=10, high=20, size=20)

# 使用yelu_walker方法计算pacf
cal_my_pacf_yw(data_x)
    

参考资料:

  1. http://www.stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YuleWalkerAndMore.htm
  2. https://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.pacf.html
  3. https://gitee.com/yuanzhoulvpi/time_series

阅读更多

在这里插入图片描述

  • 5
    点赞
  • 37
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

yuanzhoulvpi

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值