計算自相關系數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)模型，一般會有如下假設：

\[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}\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)\]

則有等式如下，\(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）的理解所總結的，存在錯誤，煩請指出，本文僅作為各位學習相關算法的參考。