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+ϕ2xi−1+⋯+ϕpxi−p+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
⎝⎜⎜⎜⎛x2x3⋮xN⎠⎟⎟⎟⎞=A
⎝⎜⎜⎜⎛x1x2⋮xN−1⎠⎟⎟⎟⎞ϕ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=1N−1xi2∑i=1N−1xixi+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+ϕ2xi−1+ϵ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 ⎝⎜⎜⎜⎛x3x4⋮xN⎠⎟⎟⎟⎞=A ⎝⎜⎜⎜⎛x2x3⋮xN−1x1x2⋮xN−2⎠⎟⎟⎟⎞Φ (ϕ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=⎣⎡(x2x1x3x2⋯⋯xN−1xN−2)⎝⎛x2x3xN−1x1x2xN−2⎠⎞⎦⎤−1=(∑i=2N−1xi2∑i=2N−1xixi−1∑i=2N−1xixi−1∑i=1N−2xi2)−1=∑i=2N−1xi2∑i=1N−2xi2−∑i=2N−1xixi−1∑i=2N−1xixi−11(∑i=1N−2xi2−∑i=2N−1xixi−1−∑i=2N−1xixi−1∑i=2N−1xi2)
接下来,假设时间序列是平稳的,所以说自协方差只和滞后多少项相关。利用这个性质,在这个案例里面可以得到这样的等式:
(
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=co2−c121(co−c1−c1co)(ATA)−1=co2(1−r12)1(co−c1−c1co)(ATA)−1=co(1−r12)1(ro−r1−r1ro)
而且:
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=(x2x1x3x2⋯⋯xN−1xN−2)⎝⎜⎜⎜⎛x3x4⋮xN⎠⎟⎟⎟⎞=(∑i=3Nxixi−1∑i=3Nxixi−2,)
又因为这个时间序列是平稳的,所以:
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(1−r12)1(ro−r1−r1ro)(c1c2)=1−r121(1−r1−r11)(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=1−r12r1(1−r2)
以及
ϕ
^
2
=
r
2
−
r
1
2
1
−
r
1
2
\hat{\phi}_{2}=\frac{r_{2}-r_{1}^{2}}{1-r_{1}^{2}}
ϕ^2=1−r12r2−r12
虽然可以按照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+ϕ2xi−1+⋯+ϕpxi−p+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=1∑p(ϕjxixi−j+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=1∑p(ϕj⟨xixi−j+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=1∑p(ϕj⟨xixi−j+1⟩)
再除以N-1,再利用自协方差的平衡性:
c
−
l
=
c
l
c_{-l} = c_{l}
c−l=cl,得:
c
1
=
∑
j
=
1
p
ϕ
j
c
j
−
1
c_{1}=\sum_{j=1}^{p} \phi_{j} c_{j-1}
c1=j=1∑pϕjcj−1
再除以
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=1∑pϕjrj−1
2.2 lag=2(滞后2)
两边乘以 x i − 1 x_{i-1} xi−1, 得到:
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} xi−1xi+1=j=1∑p(ϕjxi−1xi−j+1)+xi−1ξ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
⟨xi−1xi+1⟩=j=1∑p(ϕj⟨xi−1xi−j+1⟩)+⟨xi−1ξ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)
⟨xi−1xi+1⟩=j=1∑p(ϕj⟨xi−1xi−j+1⟩)
然后:
c
2
=
∑
j
=
1
p
ϕ
j
c
j
−
2
c_{2}=\sum_{j=1}^{p} \phi_{j} c_{j-2}
c2=j=1∑pϕjcj−2
然后:
r
2
=
∑
j
=
1
p
ϕ
j
r
j
−
2
r_{2}=\sum_{j=1}^{p} \phi_{j} r_{j-2}
r2=j=1∑pϕjrj−2
2.3 lag=k(滞后k)
两边乘以 x i − k − 1 x_{i-k-1} xi−k−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 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} xi−k+1xi+1=j=1∑p(ϕjxi−k+1xi−j+1)+xi−k+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
⟨xi−k+1xi+1⟩=j=1∑p(ϕj⟨xi−k+1xi−j+1⟩)+⟨xi−k+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)
⟨xi−k+1xi+1⟩=j=1∑p(ϕj⟨xi−k+1xi−j+1⟩)
然后:
c
k
=
∑
j
=
1
p
ϕ
j
c
j
−
k
c_{k}=\sum_{j=1}^{p} \phi_{j} c_{j-k}
ck=j=1∑pϕjcj−k
然后:
r
k
=
∑
j
=
1
p
ϕ
j
r
j
−
k
r_{k}=\sum_{j=1}^{p} \phi_{j} r_{j-k}
rk=j=1∑pϕjrj−k
2.4 lag=p(滞后p)
两边乘以 x i − p − 1 x_{i-p-1} xi−p−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 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} xi−p+1xi+1=j=1∑p(ϕjxi−p+1xi−j+1)+xi−p+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
⟨xi−p+1xi+1⟩=j=1∑p(ϕj⟨xi−p+1xi−j+1⟩)+⟨xi−p+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)
⟨xi−p+1xi+1⟩=j=1∑p(ϕj⟨xi−p+1xi−j+1⟩)
然后:
c
p
=
∑
j
=
1
p
ϕ
j
c
j
−
p
c_{p}=\sum_{j=1}^{p} \phi_{j} c_{j-p}
cp=j=1∑pϕjcj−p
然后:
r
p
=
∑
j
=
1
p
ϕ
j
r
j
−
p
r_{p}=\sum_{j=1}^{p} \phi_{j} r_{j-p}
rp=j=1∑pϕjrj−p
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+⋯+ϕp−1rp−2+ϕprp−1r2=ϕ1r1+ϕ2ro+ϕ3r1+⋯+ϕp−1rp−3+ϕprp−2⋮rp−1=ϕ1rp−2+ϕ2rp−3+ϕ3rp−4+⋯+ϕp−1ro+ϕpr1rp=ϕ1rp−1+ϕ2rp−2+ϕ3rp−3+⋯+ϕp−1r1+ϕ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)
⎝⎜⎜⎜⎜⎜⎛r1r2⋮rp−1rp⎠⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎛ror1rp−2rp−1r1ro⋮rp−3rp−2r2r1rp−4rp−3⋯⋯⋯⋯rp−2rp−3⋮ror1rp−1rp−2r1ro⎠⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛ϕ1ϕ2⋮ϕp−1ϕ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
⎝⎜⎜⎜⎜⎜⎛r1r2⋮rp−1rp⎠⎟⎟⎟⎟⎟⎞R
⎝⎜⎜⎜⎜⎜⎛1r1rp−2rp−1r11⋮rp−3rp−2r2r1rp−4rp−3⋯⋯⋯⋯rp−2rp−3⋮1r1rp−1rp−2r11⎠⎟⎟⎟⎟⎟⎞Φ
⎝⎜⎜⎜⎜⎜⎛ϕ1ϕ2⋮ϕp−1ϕ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}
Φ^=R−1r
3. Yule-Walker公式求解过程
-
循环i, 1 ≤ i ≤ p 1 \leq i \leq p 1≤i≤p
- 计算 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} 1≤j≤i−1范围内的 ϕ ^ 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)
参考资料:
- http://www.stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YuleWalkerAndMore.htm
- https://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.pacf.html
- https://gitee.com/yuanzhoulvpi/time_series