计算自相关系数acf和偏相关系数pacf
时间序列分析中,自相关系数ACF和偏相关系数PACF是两个比较重要的统计指标,在使用arma模型做序列分析时,我们可以根据这两个统计值来判断模型类型(ar还是ma)以及选择参数。目前网上关于这两个系数的资料已经相当丰富了,不过大部分内容都着重于介绍它们的含义以及使用方式,而没有对计算方法有详细的说明。所以虽然这两个系数的计算并不复杂,但是我认为还是有必要做一下总结,以便于其他人参考。本文的内容将主要集中于如何计算ACF和PACF,关于这两个系数的详细描述,大家可以参考网上的其它博客。
1. 变量说明
首先对基本变量做一下说明,后续的公式和计算都将以这些变量为准。
我们用变量\(X_{t}\)表示一个时间序列,\(x_{t}\)表示序列中的第\(t\)个点,\(t=1,2,3\dots,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}\times \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)模型,一般会有如下假设:
其中,\(\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等式:
- 首先,我们有假设的原等式:\[x_{i+1}=\sum_{j=1}^{k}(\phi_{j}x_{i-j+1})+\xi_{i+1} \]
- 将等式的两边同乘以\(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} \]
- 接着,对等式的两边同时求期望,可以得到:\[\langle x_{i-k+1}.x_{i+1}\rangle=\sum_{j=1}^{k}(\phi_{j}\langle x_{i-j+1}.x_{i-k+1}\rangle)+\langle x_{i-k+1}.\xi_{i+1}\rangle \]
- 因为噪声项默认是0均值的,所以可以去掉噪声:\[\langle x_{i-k+1}.x_{i+1}\rangle=\sum_{j=1}^{k}(\phi_{j}\langle x_{i-j+1}.x_{i-k+1}\rangle) \]
- 等式两边同除以\(N-k\)(无偏,有偏估计时,除以\(N-1\)),同时我们又有\(c_{l} = c_{-l}\),因此可以得到:\[c_{k}=\sum_{j=1}^{k}\phi_{j}c_{j-k} \]
- 最后,将等式两边同除以\(c_{0}\),可以得到:\[r_{k}=\sum_{j=1}^{k}\phi_{j}r_{j-k} \]
根据最后的等式,我们将所有相关项列出来后,可以得到:
这里的\(r_{k}\)便是之前提到的序列自相关系数ACF,同时,可以注意到\(r_{0}=1\),因此有
简化起见,我们令
则有等式如下,\(R\)为toeplitz矩阵,$$R\Phi=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)的理解所总结的,存在错误,烦请指出,本文仅作为各位学习相关算法的参考。