计算自相关系数acf和偏相关系数pacf

Posted 2021-03-17 nkuhyx

时间序列分析中，自相关系数ACF和偏相关系数PACF是两个比较重要的统计指标，在使用arma模型做序列分析时，我们可以根据这两个统计值来判断模型类型（ar还是ma）以及选择参数。目前网上关于这两个系数的资料已经相当丰富了，不过大部分内容都着重于介绍它们的含义以及使用方式，而没有对计算方法有详细的说明。所以虽然这两个系数的计算并不复杂，但是我认为还是有必要做一下总结，以便于其他人参考。本文的内容将主要集中于如何计算ACF和PACF，关于这两个系数的详细描述，大家可以参考网上的其它博客。

---

1. 变量说明

首先对基本变量做一下说明，后续的公式和计算都将以这些变量为准。

我们用变量(X_{t})表示一个时间序列，(x_{t})表示序列中的第(t)个点，(t=1,2,3dots,N)，(N)表示序列(X_{t})的长度。

序列的均值：(mu=E(X_{t}))

序列的方差：(sigma^{2}=D(X_{t})=E((X_{t}-mu)^{2}))

序列的标准差：(sigma)

对于长度一样的两条不同序列(X_{t})和(Y_{t})，可以使用协方差来刻画它们的相关性。

序列的协方差：(cov(X_{t},Y_{t})=E((X_{t}-mu_{x})(Y_{t}-mu_{y})))

协方差的值(|cov(X_{t},Y_{t})|)越大，说明序列(X_{t})和(Y_{t})的相关性越强（大于0时为正相关，小于0时为负相关）。类似地，对于序列(X_{t})，我们根据序列的滞后次数(k)来计算对应的序列自协方差，

序列的自协方差（有偏）：(hat{c}_{k}=E((X_{t}-mu)(X_{t-k}-mu))=frac{1}{N}sum_{t=k+1}^{N}(x_{t}-mu)(x_{t-k}-mu))

对于(c_{k})，我们有两种估计值，有偏估计（上式）和无偏估计，

序列的自协方差（无偏）：(c_{k}=frac{1}{N-k}sum_{t=k+1}^{N}(x_{t}-mu)(x_{t-k}-mu))

可以注意到(c_{0}(hat{c}_{0})=sigma^{2})，进一步地，我们根据序列的自协方差来定义序列的自相关系数：

序列的自相关系数（有偏）：(hat{r}_{k}=frac{hat{c}_{k}}{hat{c}_{0}})

序列的自相关系数（无偏）：(r_{k}=frac{c_{k}}{c_{0}})

后续关于PACF的计算将以无偏估计值（(c_{k})和(r_{k})）为代表，大家可自行替换为有偏估计（(hat{c}_{k})和(hat{r}_{k})）。

---

2. 序列的自相关系数ACF

由前文易得ACF的计算公式：

自相关系数ACF（无偏）：(acf(k) = r_{k}=frac{c_{k}}{c_{0}}=frac{N}{N-k} imes frac{sum_{t=k+1}^{N}(x_{t}-mu)(x_{t-k}-mu)}{sum_{t=1}^{N}(x_{t}-mu)(x_{t}-mu)})

自相关系数ACF（有偏）：(hat{acf}(k) = hat{r}_{k}=frac{(N-k)c_{k}}{Nc_{0}}=frac{sum_{t=k+1}^{N}(x_{t}-mu)(x_{t-k}-mu)}{sum_{t=1}^{N}(x_{t}-mu)(x_{t}-mu)})

代码（python）如下

import numpy as np

# acf method in statsmodels
#import statsmodels.tsa.stattools as stattools
#def default_acf(ts, k):
#    return statools.acf(ts, nlags=k, unbiased=False)

def acf(ts, k):
    """ Compute autocorrelation coefficient, biased
    """
    x = np.array(ts) - np.mean(ts)
    coeff = np.zeros(k+1, np.float64) # to store acf
    coeff[0] = x.dot(x) # N*c(0)

    for i in range(1, k+1):
        coeff[i] = x[:-i].dot(x[i:]) # (N-k)*c(i)
        
    return coeff / coeff[0]

---

3. 序列的偏相关系数PACF

偏相关系数PACF的计算相较于自相关系数ACF要复杂一些。网上大部分资料都只给出了PACF的公式和理论说明，对于PACF的值则没有具体的介绍，所以我们首先需要说明一下PACF指的是什么。这里我们借助AR模型来说明，对于AR(p)模型，一般会有如下假设：
[x_{i+1} = phi_{1}x_{i}+phi_{2}x_{i-1}+...phi_{p}x_{i-p+1}+xi_{i+1}]
其中，(phi_{j},j=1,2,dots,P)是线性相关系数，(xi_{i+1})是噪声，即我们假设点(x_{i+1})与前(p)个点(x_{i-p+1},x_{i-p+2},dots,x_{i})是线性相关的，而PACF所要表示的就是点(x_{i})与点(x_{i-p})的相关性。所以，

序列的偏相关系数PACF：(pacf(p)=phi_{p})

我们没有办法单独求解(phi_{p})。对于一般的线性回归问题，可以使用最小二乘法（MLS）来求解相关系数，而这里可以使用Yule-Walker等式来求解，详情可以参考YW-Eshel。对于滞后次数(k)，我们依照如下过程来构建Yule-Walker等式：

1. 首先，我们有假设的原等式:
   [x_{i+1}=sum_{j=1}^{k}(phi_{j}x_{i-j+1})+xi_{i+1}]
2. 将等式的两边同乘以(x_{i-k+1})，可以得到：
   [x_{i-k+1}.x_{i+1}=sum_{j=1}^{k}(phi_{j}x_{i-j+1}.x_{i-k+1})+x_{i-k+1}.xi_{i+1}]
3. 接着，对等式的两边同时求期望，可以得到：
   [langle x_{i-k+1}.x_{i+1} angle=sum_{j=1}^{k}(phi_{j}langle x_{i-j+1}.x_{i-k+1} angle)+langle x_{i-k+1}.xi_{i+1} angle]
4. 因为噪声项默认是0均值的，所以可以去掉噪声：
   [langle x_{i-k+1}.x_{i+1} angle=sum_{j=1}^{k}(phi_{j}langle x_{i-j+1}.x_{i-k+1} angle)]
5. 等式两边同除以(N-k)（无偏，有偏估计时，除以(N-1)），同时我们又有(c_{l} = c_{-l})，因此可以得到：
   [c_{k}=sum_{j=1}^{k}phi_{j}c_{j-k}]
6. 最后，将等式两边同除以(c_{0})，可以得到：
   [r_{k}=sum_{j=1}^{k}phi_{j}r_{j-k}]

根据最后的等式，我们将所有相关项列出来后，可以得到：
[left(egin{matrix} r_{1} r_{2} . r_{k-1} r_{k} end{matrix} ight) = left(egin{matrix} r_{0} & r_{1} & r_{2} & . & r_{k-2} & r_{k-1} r_{1} & r_{0} & r_{1} & . & r_{k-3} & r_{k-2} . & . & . & . & . & . r_{k-2} & r_{k-3} & r_{k-4} & . & r_{0} & r_{1} r_{k-1} & r_{k-2} & r_{k-3} & . & r_{1} & r_{0} end{matrix} ight) imes left(egin{matrix} phi_{1} phi_{2} . phi_{k-1} phi_{k} end{matrix} ight)]
这里的(r_{k})便是之前提到的序列自相关系数ACF，同时，可以注意到(r_{0}=1)，因此有
[left(egin{matrix} r_{1} r_{2} . r_{k-1} r_{k} end{matrix} ight) = left( egin{matrix} 1 & r_{1} & r_{2} & . & r_{k-2} & r_{k-1} r_{1} &1 & r_{1} & . & r_{k-3} & r_{k-2} . & . & . & . & . & . r_{k-2} & r_{k-3} & r_{k-4} & . &1 & r_{1} r_{k-1} & r_{k-2} & r_{k-3} & . & r_{1} &1 end{matrix} ight) imes left(egin{matrix} phi_{1} phi_{2} . phi_{k-1} phi_{k} end{matrix} ight)]
简化起见，我们令
[r=left(egin{matrix} r_{1} r_{2} . r_{k-1} r_{k} end{matrix} ight), R= left(egin{matrix} 1 & r_{1} & r_{2} & . & r_{k-2} & r_{k-1} r_{1} &1 & r_{1} & . & r_{k-3} & r_{k-2} . & . & . & . & . & . r_{k-2} & r_{k-3} & r_{k-4} & . &1 & r_{1} r_{k-1} & r_{k-2} & r_{k-3} & . & r_{1} &1 end{matrix} ight), Phi= left(egin{matrix} phi_{1} phi_{2} . phi_{k-1} phi_{k} end{matrix} ight)]
则有等式如下，(R)为toeplitz矩阵，[RPhi=r]
由于矩阵(R)是对称满秩矩阵，所以存在可逆矩阵(R^{-1})，使得(hat{Phi}=R^{-1}r)，则可以求得滞后次数为(k)的偏相关系数（上标((k))表示第(k)个解向量）：[pacf(k)=hat{Phi}_{k}^{(k)}]
代码（python）如下

import numpy as np
from scipy.linalg import toeplitz

# pacf method in statsmodels
#import statsmodels.tsa.stattools as stattools
#def default_pacf(ts, k):
#    return statools.pacf(ts, nlags=k, unbiased=True)

def yule_walker(ts, order):
    ''' Solve yule walker equation
    '''
    x = np.array(ts) - np.mean(ts)
    n = x.shape[0]

    r = np.zeros(order+1, np.float64) # to store acf
    r[0] = x.dot(x) / n # r(0)
    for k in range(1, order+1):
        r[k] = x[:-k].dot(x[k:]) / (n - k) # r(k)

    R = toeplitz(r[:-1])

    return np.linalg.solve(R, r[1:]) # solve `Rb = r` to get `b`

def pacf(ts, k):
    ''' Compute partial autocorrelation coefficients for given time series，unbiased
    '''
    res = [1.]
    for i in range(1, k+1):
        res.append(yule_walker(ts, i)[-1])
    return np.array(res)

---

4. 总结

对如何计算自相关系数ACF和偏相关系数PACF的介绍就到这里了，由于本人并非统计学相关专业，上述内容是基于个人对网上资料和开源代码（python：statsmodels）的理解所总结的，存在错误，烦请指出，本文仅作为各位学习相关算法的参考。