介绍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=1p=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个自协方差和自协方差系数。
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公式。
这是一个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
在左右两边乘上
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
⟩
leftlangle x_{i} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i} x_{i-j+1}rightrangleright)+leftlangle x_{i} xi_{i+1}rightrangle
⟨xixi+1⟩=j=1∑p(ϕj⟨xixi−j+1⟩)+⟨xiξi+1⟩
注意到
⟨
x
i
ξ
i
+
1
⟩
=
0
leftlangle x_{i} xi_{i+1}rightrangle = 0
⟨xiξi+1⟩=0 因为截距
ξ
xi
ξ和当前时间是无关的，因此：
⟨
x
i
x
i
+
1
⟩
=
∑
j
=
1
p
(
ϕ
j
⟨
x
i
x
i
−
j
+
1
⟩
)
leftlangle x_{i} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i} x_{i-j+1}rightrangleright)
⟨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
两边乘以 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
⟩
leftlangle x_{i-1} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i-1} x_{i-j+1}rightrangleright)+leftlangle x_{i-1} xi_{i+1}rightrangle
⟨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
⟩
)
leftlangle x_{i-1} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i-1} x_{i-j+1}rightrangleright)
⟨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
两边乘以 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
⟩
leftlangle x_{i-k+1} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i-k+1} x_{i-j+1}rightrangleright)+leftlangle x_{i-k+1} xi_{i+1}rightrangle
⟨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
⟩
)
leftlangle x_{i-k+1} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i-k+1} x_{i-j+1}rightrangleright)
⟨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
两边乘以 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
⟩
leftlangle x_{i-p+1} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i-p+1} x_{i-j+1}rightrangleright)+leftlangle x_{i-p+1} xi_{i+1}rightrangle
⟨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
⟩
)
leftlangle x_{i-p+1} x_{i+1}rightrangle=sum_{j=1}^{p}left(phi_{j}leftlangle x_{i-p+1} x_{i-j+1}rightrangleright)
⟨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
有：
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
-
循环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
代码如下：
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
