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§参数，怎么做？
我们要估计的也就是下面公式中的${\xi }_{i+1}$

$$x_{i+1}={\varphi }_1{x}_i+{\varphi }_2{x}_{i-1}+\cdots +{\varphi }_p{x}_{i-p+1}+{\xi }_{i+1}\text{\hspace{1em}(1)}$$

1.转置

1.1 p=1

p=1的时候，公式变成这样$$X_{i+1}={\varphi }_1{x}_i+{ϵ}_{i+1}$$
上面的形式可以进一步写成这样：
$$\underset{\mathbf{b}}{\underbrace{\left(\begin{array}{c}{x}_2\\ x_3\\ ⋮\\ x_N\end{array}\right)}}=\underset{\mathbf{A}}{\underbrace{\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{N-1}\end{array}\right)}}{\varphi }_1$$

可以使用最小二乘法求解：
$${\hat{\varphi }}_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$$
上面的$c_i$、$r_i$分别是第i个自协方差和自协方差系数。

1.2 p=2

p=2的时候，公式变成这样$$X_{i+1}={\varphi }_1{x}_i+{\varphi }_2{x}_{i-1}+{ϵ}_{i+1}$$
上面的形式进一步写成这样：

$$\underset{\mathbf{b}}{\underbrace{\left(\begin{array}{c}{x}_3\\ x_4\\ ⋮\\ x_N\end{array}\right)}}=\underset{\mathbf{A}}{\underbrace{\left(\begin{array}{cc}{x}_2& x_1\\ x_3& x_2\\ ⋮& ⋮\\ x_{N-1}& x_{N-2}\end{array}\right)}}\underset{\Phi }{\underbrace{\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right)}}.$$

和第一个不一样，这个时候上面的等式写成这样：
$$\hat{\mathbf{\Phi }}={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathbf{b}$$
是这么计算的：
$$\begin{array}{c}{\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{array}$$

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

$$\begin{array}{c}{\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{array}$$
而且：
$${\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\\ ⋮\\ 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)$$

又因为这个时间序列是平稳的，所以：

$${\mathbf{A}}^T\mathbf{b}=\left(\begin{array}{l}{c}_1\\ c_2\end{array}\right)$$

结合上面的公式，可以有：
$$\begin{array}{c}{\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{array}$$
将上面的矩阵分为两个部分，可以得到：
$${\hat{\varphi }}_1=\frac{r_1\left(1-r_2\right)}{1-r_1^2}$$

以及
$${\hat{\varphi }}_2=\frac{r_2-r_1^2}{1-r_1^2}$$

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

2. Yule-Waler公式

这是一个AR§公式：
$$x_{i+1}={\varphi }_1{x}_i+{\varphi }_2{x}_{i-1}+\cdots +{\varphi }_p{x}_{i-p+1}+{\xi }_{i+1}$$

2.1 lag=1（滞后1）

在左右两边乘上$x_i$，得
$$x_i{x}_{i+1}=\sum _{j=1}^p\left({\varphi }_j{x}_i{x}_{i-j+1}\right)+x_i{\xi }_{i+1}$$

i和j是各自的时间索引。把 { ϕ j } \{\phi_j\} {ϕj​}都拎出来，$x_j$都写到一起，得到这样的公式：
$$⟨x_i{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_i{x}_{i-j+1}⟩\right)+⟨x_i{\xi }_{i+1}⟩$$

注意到$⟨x_i{\xi }_{i+1}⟩=0$ 因为截距$\xi$和当前时间是无关的，因此：
$$⟨x_i{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_i{x}_{i-j+1}⟩\right)$$

再除以N-1，再利用自协方差的平衡性:$c_{-l}=c_l$，得：
$$c_1=\sum _{j=1}^p{\varphi }_j{c}_{j-1}$$

再除以$c_0$，得：
$$r_1=\sum _{j=1}^p{\varphi }_j{r}_{j-1}$$

2.2 lag=2（滞后2）

两边乘以$x_{i-1}$, 得到：

$$x_{i-1}{x}_{i+1}=\sum _{j=1}^p\left({\varphi }_j{x}_{i-1}{x}_{i-j+1}\right)+x_{i-1}{\xi }_{i+1}$$

然后：
$$⟨x_{i-1}{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_{i-1}{x}_{i-j+1}⟩\right)+⟨x_{i-1}{\xi }_{i+1}⟩$$

然后：
$$⟨x_{i-1}{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_{i-1}{x}_{i-j+1}⟩\right)$$

然后：
$$c_2=\sum _{j=1}^p{\varphi }_j{c}_{j-2}$$

然后：
$$r_2=\sum _{j=1}^p{\varphi }_j{r}_{j-2}$$

2.3 lag=k(滞后k)

两边乘以$x_{i-k-1}$, 得到：

$$x_{i-k+1}{x}_{i+1}=\sum _{j=1}^p\left({\varphi }_j{x}_{i-k+1}{x}_{i-j+1}\right)+x_{i-k+1}{\xi }_{i+1}$$

然后：
$$⟨x_{i-k+1}{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_{i-k+1}{x}_{i-j+1}⟩\right)+⟨x_{i-k+1}{\xi }_{i+1}⟩$$

然后：
$$⟨x_{i-k+1}{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_{i-k+1}{x}_{i-j+1}⟩\right)$$

然后：
$$c_k=\sum _{j=1}^p{\varphi }_j{c}_{j-k}$$

然后：
$$r_k=\sum _{j=1}^p{\varphi }_j{r}_{j-k}$$

2.4 lag=p(滞后p)

两边乘以$x_{i-p-1}$, 得到：

$$x_{i-p+1}{x}_{i+1}=\sum _{j=1}^p\left({\varphi }_j{x}_{i-p+1}{x}_{i-j+1}\right)+x_{i-p+1}{\xi }_{i+1}$$

然后：
$$⟨x_{i-p+1}{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_{i-p+1}{x}_{i-j+1}⟩\right)+⟨x_{i-p+1}{\xi }_{i+1}⟩$$

然后：
$$⟨x_{i-p+1}{x}_{i+1}⟩=\sum _{j=1}^p\left({\varphi }_j⟨x_{i-p+1}{x}_{i-j+1}⟩\right)$$

然后：
$$c_p=\sum _{j=1}^p{\varphi }_j{c}_{j-p}$$

然后：
$$r_p=\sum _{j=1}^p{\varphi }_j{r}_{j-p}$$

2.5 将上面的公式放一起

有：
$$\begin{array}{c}{r}_1={\varphi }_1{r}_o+{\varphi }_2{r}_1+{\varphi }_3{r}_2+\cdots +{\varphi }_{p-1}{r}_{p-2}+{\varphi }_p{r}_{p-1}\\ r_2={\varphi }_1{r}_1+{\varphi }_2{r}_o+{\varphi }_3{r}_1+\cdots +{\varphi }_{p-1}{r}_{p-3}+{\varphi }_p{r}_{p-2}\\ ⋮\\ r_{p-1}={\varphi }_1{r}_{p-2}+{\varphi }_2{r}_{p-3}+{\varphi }_3{r}_{p-4}+\cdots +{\varphi }_{p-1}{r}_o+{\varphi }_p{r}_1\\ r_p={\varphi }_1{r}_{p-1}+{\varphi }_2{r}_{p-2}+{\varphi }_3{r}_{p-3}+\cdots +{\varphi }_{p-1}{r}_1+{\varphi }_p{r}_o\end{array}$$
可以被写成：
$$\left(\begin{array}{c}{r}_1\\ r_2\\ ⋮\\ 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}\\ & ⋮& & & ⋮& \\ 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}{\varphi }_1\\ {\varphi }_2\\ ⋮\\ {\varphi }_{p-1}\\ {\varphi }_p\end{array}\right)$$

又因为$r_0=1$，所以可以写成：
$$\underset{\mathbf{r}}{\underbrace{\left(\begin{array}{c}{r}_1\\ r_2\\ ⋮\\ r_{p-1}\\ r_p\end{array}\right)}}\underset{\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}\\ & ⋮& & & ⋮& \\ 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)}}\underset{\mathbf{\Phi }}{\underbrace{\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ ⋮\\ {\varphi }_{p-1}\\ {\varphi }_p\end{array}\right)}}$$
整理成：
$$\mathbf{R}\mathbf{\Phi }=\mathbf{r}\text{\hspace{1em}(2)}$$

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

3. Yule-Walker公式求解过程

• 循环i，$1\le i\le p$

  • 计算${\mathbf{R}}^{(i)}$ 和 ${\mathbf{r}}^{(i)}$
  • 然后计算${\hat{\mathbf{\Phi }}}^{(i)}$，公式为：
    $${\hat{\mathbf{\Phi }}}^{(i)}={\left({\mathbf{R}}^{(i)}\right)}^{-1}{\mathbf{r}}^{(i)}=\left(\begin{array}{c}{\hat{\varphi }}_1\\ {\hat{\varphi }}_2\\ ⋮\\ {\hat{\varphi }}_i\end{array}\right)$$
  • 只保留${\hat{\varphi }}_i$，在$1\le j\le i-1$范围内的${\hat{\varphi }}_j$都不要
  • 第i个pacf系数这个时候就等于$pacf(i)={\hat{\varphi }}_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

阅读更多