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

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


1. 变量说明

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

我们用变量X_{t}表示一个时间序列,x_{t}表示序列中的第t个点,t=1,2,3\dots,NN表示序列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)模型,一般会有如下假设:
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}的相关性,所以,

序列的偏相关系数PACFpacf(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}\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
  4. 因为噪声项默认是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)
  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(\begin{matrix} r_{1}\\ r_{2}\\ .\\ r_{k-1}\\ r_{k} \end{matrix}\right) = \left(\begin{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}\right)\times \left(\begin{matrix} \phi_{1}\\ \phi_{2}\\ .\\ \phi_{k-1}\\ \phi_{k} \end{matrix}\right)
这里的r_{k}便是之前提到的序列自相关系数ACF,同时,可以注意到r_{0}=1,因此有
\left(\begin{matrix} r_{1}\\ r_{2}\\ .\\ r_{k-1}\\ r_{k} \end{matrix}\right) = \left( \begin{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}\right)\times \left(\begin{matrix} \phi_{1}\\ \phi_{2}\\ .\\ \phi_{k-1}\\ \phi_{k} \end{matrix}\right)
简化起见,我们令
r=\left(\begin{matrix} r_{1}\\ r_{2}\\ .\\ r_{k-1}\\ r_{k} \end{matrix}\right), R= \left(\begin{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}\right), \Phi= \left(\begin{matrix} \phi_{1}\\ \phi_{2}\\ .\\ \phi_{k-1}\\ \phi_{k} \end{matrix}\right)
则有等式如下,Rtoeplitz矩阵,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)的理解所总结的,存在错误,烦请指出,本文仅作为各位学习相关算法的参考。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,743评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,296评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,285评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,485评论 1 283
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,581评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,821评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,960评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,719评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,186评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,516评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,650评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,329评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,936评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,757评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,991评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,370评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,527评论 2 349