時間序列 ACF 和 PACF 理解、程式碼、視覺化

時間序列系列的相關介紹，從零梳理時序概念、相關技術、和實戰案例，歡迎訂閱 ?時間序列專欄 跟蹤全部內容。

我們說時間序列可以被預測，主要基於以下事實：我們可以部分掌握影響該時間序列的因素的變化情況。換句話說，對時間序列進行預測，其實就是利用各種理論和工具，對觀察到的時間序列進行“抽絲剝繭”，以試圖掌握其變化的本質，從而對未來的表現進行預測。

而自相關性是時序預測的基礎，對於時序的平穩性、白噪聲檢測、確定 模型中的階數（p/q）有著重要的作用。本篇將著重介紹自相關的概念 和 。

ACF 自相關函式

概念理解

ACF（Autocorrelation Function）就是用來計算時間序列自身的相關性的函式。

對於同一時間 的計算，，這個很好理解。

如果是不同的時間，比如 ，該如何計算呢？

實際上，在應用自相關函式時，其輸入分別為原始的時間序列 及其 階滯後序列 ，於是 就變成了：

，這裡兩個序列的長度是一致的，如下圖所示：

計算和程式碼

ACF的公式定義為：

Python程式碼實現可以直接使用statsmodels包進行計算，當然也可以自己透過Numpy復現一遍公式，結果是一樣的。

import numpy as np

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

# 手擼公式計算acf，有偏
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]

視覺化

透過視覺化可以更清楚的看出不同lag的係數值和趨勢變化，透過statsmodels函式的直接繪製，以下是示例。

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
df = pd.read_excel('data.xlsx',engine='openpyxl')

df['price'].plot()
plot_acf(df['price'], lags=40, adjusted=False)
plt.show()

第一個圖是一組時間序列的資料。第二個圖是計算的ACF相關係數圖。

ACF圖的橫座標表示滯後的階數，縱座標表示對應的滯後序列與原始序列的相關係數。可以看出，隨著滯後階數的增加，滯後序列與原始序列的相關性也在不斷地降低。圖中的藍色區域表示置信區間，用來標識相關係數是否具有統計顯著性。簡單來說，如果相關係數落在置信區間內，表明對應的兩個序列的相關係數並不能代表其真實相關性。

即使是兩個完全不相干的白噪聲序列，由於隨機性的影響，其相關係數也不可能全都為0，因此，需要使用置信區間來過濾掉那些由於隨機性造成的“偽相關”。

PACF 偏自相關函式

概念理解

我們知道求導是對所有項都求導，求偏導只對某一個求導忽略其他項。 和 也可以理解為這樣的關係。

前面我們計算 自相關函式時，得到的並不是 與 之間單純的相關關係。因為 同時還會受到中間 個隨機變數 、... 的影響，而這 個隨機變數又都和 具有相關關係，所以自相關係數裡面實際摻雜了其他變數對 與 的影響。

為了得到  對 的直接影響，引入了偏自相關係數 的概念。滯後 偏自相關係數是指，對於平穩時間序列 ，在剔除了中間 個隨機變數 、... 的干擾之後， 對 影響的純相關程度。

計算和程式碼

的計算比 要複雜很多。這裡我們藉助AR模型來說明，對於AR(p)模型，一般會有如下假設：

其中， 是線性相關係數， 是噪聲，即我們假設點 與前 個點 是線性相關的。而 所要表示的就是點 與點 的相關性，所以，

序列的偏相關係數PACF：

有幾種方法可以求解相關係數，方法包括最小二乘法（MLS）、尤爾-沃克方程（Yule-Walker equation）、伯格演算法（Burg"s method。由於公式推導內容較多，本篇對求解方法不做詳細介紹。

Python計算程式碼如下：

import numpy as np
from scipy.linalg import toeplitz

# 使用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)

視覺化

使用statsmodels函式直接繪製，methond可以選擇合適的方法求解，下面使用最小二乘法ols進行求解。

# 使用最小二乘法ols求解
plot_pacf(df['price'], lags=40,  method='ols')
plt.show()

以上就是對 和 的介紹，理解自相關的概念對於學習時間序列非常重要，下一篇將介紹它們的應用場景。

參考連結

[1].https://blog.csdn.net/SunJW_2017/article/details/126993853，芳樽裡的歌
[2].，洪於祥 [3].，gwave [4]. [5].https://mp.weixin.qq.com/s/llMZaMkhoXLRDlFxoFlXiw，seriesc