【时间序列】python与时间序列3(Yule-Walker方程:计算时间序列神器最详细的公式推导 python代码复现)...

客官°小女子只卖身不卖艺 2022-09-16 11:18 148阅读 0赞

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

## 介绍 ##

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

**在找pacf如何计算的时候，发现现在计算pacf都在用Yule-Walker方法，包括在计算AR(p)的时候也是这个方程（不是不用别的方法了）。**

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

最新的代码在gitee仓库里面：**【点击文末左下角的阅读原文】**

# 使用Yule-Walker方法计算AR(p) #

计算一个均值为0的离散随机时间序列的AR(p)参数，怎么做？我们要估计的也就是下面公式中的

# 1.转置 #

## 1.1 p=1 ##

p=1的时候，公式变成这样

上面的形式可以进一步写成这样：

可以使用最小二乘法求解：

上面的、分别是第i个自协方差和自协方差系数。

## 1.2 p=2 ##

p=2的时候，公式变成这样

上面的形式进一步写成这样：

和第一个不一样，这个时候上面的等式写成这样：

是这么计算的：

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

而且：

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

结合上面的公式，可以有：

将上面的矩阵分为两个部分，可以得到：

以及

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

# 2. Yule-Waler公式 #

这是一个AR(p)公式：

## 2.1 lag=1（滞后1） ##

在左右两边乘上，得

i和j是各自的时间索引。把都拎出来，都写到一起，得到这样的公式：

注意到 因为截距和当前时间是无关的，因此：

再除以N-1，再利用自协方差的平衡性:，得：

再除以，得：

## 2.2 lag=2（滞后2） ##

两边乘以, 得到：

然后：

## 2.3 lag=k(滞后k) ##

两边乘以, 得到：

然后：

## 2.4 lag=p(滞后p) ##

两边乘以, 得到：

然后：

## 2.5 将上面的公式放一起 ##

有：

可以被写成：

又因为，所以可以写成：

整理成：

最终可以写成这样

# 3. Yule-Walker公式求解过程 #

*  循环i，
    
     *  计算 和
     *  然后计算，公式为：
     *  只保留，在范围内的都不要。
     *  第i个pacf系数这个时候就等于
 *  结束循环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/tree/master

公众号：AI蜗牛车
    
    保持谦逊、保持自律、保持进步