計(jì)算自相關(guān)系數(shù)acf和偏相關(guān)系數(shù)pacf

時(shí)間序列分析中侧啼，自相關(guān)系數(shù)ACF和偏相關(guān)系數(shù)PACF是兩個(gè)比較重要的統(tǒng)計(jì)指標(biāo)钥屈，在使用arima模型做序列預(yù)測時(shí)，我們可以根據(jù)這兩個(gè)統(tǒng)計(jì)值來判斷模型類型（ar還是ma）以及選擇參數(shù)队萤。目前網(wǎng)上（百度）關(guān)于這兩個(gè)系數(shù)的資料已經(jīng)相當(dāng)豐富了下愈，不過大部分的內(nèi)容都只著重于介紹它們的含義以及使用方式击蹲，并沒有對計(jì)算方法有詳細(xì)的說明。所以雖然這兩個(gè)系數(shù)的計(jì)算比較簡單婉宰，但是我認(rèn)為還是有必要做一下總結(jié)說明歌豺，以便于其他人作為參考。本文的內(nèi)容將主要集中于如何計(jì)算ACF和PACF心包，關(guān)于這兩個(gè)系數(shù)的介紹和說明类咧，大家可以參考網(wǎng)上的博客。

1. 變量說明

首先對基本變量做一下說明蟹腾，后續(xù)的公式和計(jì)算都將以這些變量為準(zhǔn)痕惋。

我們用變量 $X_{t}$ 表示一個(gè)時(shí)間序列， $x_{t}$ 表示序列中的第 $t$ 個(gè)點(diǎn)娃殖， $t=1,2,3\dots,N$ 值戳， $N$ 表示序列 $X_{t}$ 的長度。

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

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

序列的標(biāo)準(zhǔn)差： $\sigma$

對于長度一樣的兩條不同序列 $X_{t}$ 和 $Y_{t}$ 炉爆，可以使用協(xié)方差來刻畫它們的相關(guān)性堕虹。

序列的協(xié)方差： $cov(X_{t},Y_{t})=E((X_{t}-\mu_{x})(Y_{t}-\mu_{y}))$

協(xié)方差的值 $|cov(X_{t},Y_{t})|$ 越大，說明序列 $X_{t}$ 和 $Y_{t}$ 的相關(guān)性越強(qiáng)（大于0時(shí)為正相關(guān)叶洞，小于0時(shí)為負(fù)相關(guān)）鲫凶。類似地禀崖，對于序列 $X_{t}$ 衩辟，我們根據(jù)序列的滯后次數(shù) $k$ 來計(jì)算對應(yīng)的序列自協(xié)方差，

序列的自協(xié)方差（有偏）： $\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}$ 波附，我們有兩種估計(jì)值艺晴，有偏估計(jì)（上式）和無偏估計(jì)，

序列的自協(xié)方差（無偏）： $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}$ 掸屡，進(jìn)一步地封寞，我們根據(jù)序列的自協(xié)方差來定義序列的自相關(guān)系數(shù)：

序列的自相關(guān)系數(shù)（有偏）： $\hat{r}_{k}=\frac{\hat{c}_{k}}{\hat{c}_{0}}$

序列的自相關(guān)系數(shù)（無偏）： $r_{k}=\frac{c_{k}}{c_{0}}$

后續(xù)關(guān)于PACF的計(jì)算將以無偏估計(jì)值（ $c_{k}$ 和 $r_{k}$ ）為代表，大家可自行替換為有偏估計(jì)（ $\hat{c}_{k}$ 和 $\hat{r}_{k}$ ）仅财。

2. 序列的自相關(guān)系數(shù)ACF

由前文易得ACF的計(jì)算公式：

自相關(guān)系數(shù)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)}$

自相關(guān)系數(shù)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. 序列的偏相關(guān)系數(shù)PACF

偏相關(guān)系數(shù)PACF的計(jì)算相較于自相關(guān)系數(shù)ACF要復(fù)雜一些狈究。網(wǎng)上大部分資料都只給出了PACF的公式和理論說明，對于PACF的值則沒有具體的介紹盏求，所以我們首先需要說明一下PACF指的是什么抖锥。這里我們借助AR模型來說明，對于AR(p)模型碎罚，一般會(huì)有如下假設(shè)：
$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$ 是線性相關(guān)系數(shù)， $\xi_{i+1}$ 是噪聲荆烈，即我們假設(shè)點(diǎn) $x_{i+1}$ 與前 $p$ 個(gè)點(diǎn) $x_{i-p+1},x_{i-p+2},\dots,x_{i}$ 是線性相關(guān)的拯勉。而PACF所要表示的就是點(diǎn) $x_{i}$ 與點(diǎn) $x_{i-p}$ 的相關(guān)性竟趾，所以，

序列的偏相關(guān)系數(shù)PACF： $pacf(p)=\phi_{p}$

我們沒有辦法單獨(dú)求解 $\phi_{p}$ 宫峦。對于一般的線性回歸問題岔帽，可以使用最小二乘法（MLS）來求解相關(guān)系數(shù)，而這里可以使用Yule-Walker等式來求解斗遏，詳情可以參考YW-Eshel山卦。對于滯后次數(shù) $k$ ，我們依照如下過程來構(gòu)建Yule-Walker等式：

首先诵次，我們有假設(shè)的原等式:
$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}$
接著，對等式的兩邊同時(shí)求期望逾一，可以得到：
$\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$
因?yàn)樵肼曧?xiàng)默認(rèn)是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$ （無偏，有偏估計(jì)時(shí)遵堵，除以 $N-1$ ）箱玷，同時(shí)我們又有 $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}$

根據(jù)最后的等式，我們將所有相關(guān)項(xiàng)列出來后壳坪，可以得到：
$\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}$ 便是之前提到的序列自相關(guān)系數(shù)ACF舶得，同時(shí)，可以注意到 $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)$
則有等式如下， $R$ 為toeplitz矩陣蝎亚， $R\Phi=r$
由于矩陣 $R$ 是對稱滿秩矩陣九孩，所以存在可逆矩陣 $R^{-1}$ ，使得 $\hat{\Phi}=R^{-1}r$ 发框，則可以求得滯后次數(shù)為 $k$ 的偏相關(guān)系數(shù)（上標(biāo) $(k)$ 表示第 $k$ 個(gè)解向量）： $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. 總結(jié)

對如何計(jì)算自相關(guān)系數(shù)ACF和偏相關(guān)系數(shù)PACF的介紹就到這里了，由于本人并非統(tǒng)計(jì)學(xué)相關(guān)專業(yè)梅惯，上述內(nèi)容是基于個(gè)人對網(wǎng)上資料和開源代碼（python：statsmodels）的理解所總結(jié)的宪拥，存在錯(cuò)誤，煩請指出个唧，本文僅作為各位學(xué)習(xí)相關(guān)算法的參考江解。