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?=?1?xi?+?2?xi?1?+?+?p?xi?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?=?1?xi?+?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
??????x2?x3??xN??????????=A
??????x1?x2??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?1?xi2?∑i=1N?1?xi?xi+1??=co?c1??=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?=?1?xi?+?2?xi?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 ??????x3?x4??xN??????????=A ??????x2?x3??xN?1??x1?x2??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=???(x2?x1??x3?x2?????xN?1?xN?2??)???x2?x3?xN?1??x1?x2?xN?2?????????1=(∑i=2N?1?xi2?∑i=2N?1?xi?xi?1??∑i=2N?1?xi?xi?1?∑i=1N?2?xi2??)?1=∑i=2N?1?xi2?∑i=1N?2?xi2??∑i=2N?1?xi?xi?1?∑i=2N?1?xi?xi?1?1?(∑i=1N?2?xi2??∑i=2N?1?xi?xi?1???∑i=2N?1?xi?xi?1?∑i=2N?1?xi2??)?
接下来，假设时间序列是平稳的，所以说自协方差只和滞后多少项相关。利用这个性质，在这个案例里面可以得到这样的等式：
(
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??c12?1?(co??c1???c1?co??)(ATA)?1=co2?(1?r12?)1?(co??c1???c1?co??)(ATA)?1=co?(1?r12?)1?(ro??r1???r1?ro??)?
而且：
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=(x2?x1??x3?x2?????xN?1?xN?2??)??????x3?x4??xN????????=(∑i=3N?xi?xi?1?∑i=3N?xi?xi?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=(c1?c2??)
结合上面的公式，可以有：
(
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???r1?ro??)(c1?c2??)=1?r12?1?(1?r1???r1?1?)(r1?r2??).?
将上面的矩阵分为两个部分，可以得到：
?
^
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?r12?r1?(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?r12?r2??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?=?1?xi?+?2?xi?1?+?+?p?xi?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}
xi?xi+1?=j=1∑p?(?j?xi?xi?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
?xi?xi+1??=j=1∑p?(?j??xi?xi?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)
?xi?xi+1??=j=1∑p?(?j??xi?xi?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??j?cj?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??j?rj?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?1?xi+1?=j=1∑p?(?j?xi?1?xi?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?1?xi+1??=j=1∑p?(?j??xi?1?xi?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?1?xi+1??=j=1∑p?(?j??xi?1?xi?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??j?cj?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??j?rj?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+1?xi+1?=j=1∑p?(?j?xi?k+1?xi?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+1?xi+1??=j=1∑p?(?j??xi?k+1?xi?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+1?xi+1??=j=1∑p?(?j??xi?k+1?xi?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??j?cj?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??j?rj?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+1?xi+1?=j=1∑p?(?j?xi?p+1?xi?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+1?xi+1??=j=1∑p?(?j??xi?p+1?xi?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+1?xi+1??=j=1∑p?(?j??xi?p+1?xi?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??j?cj?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??j?rj?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?=?1?ro?+?2?r1?+?3?r2?+?+?p?1?rp?2?+?p?rp?1?r2?=?1?r1?+?2?ro?+?3?r1?+?+?p?1?rp?3?+?p?rp?2??rp?1?=?1?rp?2?+?2?rp?3?+?3?rp?4?+?+?p?1?ro?+?p?r1?rp?=?1?rp?1?+?2?rp?2?+?3?rp?3?+?+?p?1?r1?+?p?ro??
可以被写成：
(
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)
????????r1?r2??rp?1?rp??????????=????????ro?r1?rp?2?rp?1??r1?ro??rp?3?rp?2??r2?r1?rp?4?rp?3???????rp?2?rp?3??ro?r1??rp?1?rp?2?r1?ro???????????????????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
????????r1?r2??rp?1?rp????????????R
????????1r1?rp?2?rp?1??r1?1?rp?3?rp?2??r2?r1?rp?4?rp?3???????rp?2?rp?3??1r1??rp?1?rp?2?r1?1???????????Φ
?????????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
阅读更多
