利用python進行時間序列分析——數據挖掘

09-05

來自專欄 Eureka機器學習讀書筆記7 人贊了文章

1.時間序列常用模型

1 平滑法：用於趨勢分析與預測，削弱短息隨機波動對序列的影響，平滑技術有移動平均法和指數平滑法

2 趨勢擬合法：把時間作為自變數，觀察序列值為因變數，建立回歸模型，具體可分為線性和曲線擬合

3 組合模型：時間序列趨勢主要收到長期趨勢（T）、季節變動（S）、周期波動（C）和不規則波動（ $varepsilon$ ）四個因素影響，可構建

加法模型： $x_{t}=T_{t}+S_{t}+C_{t}+varepsilon_ {t}$
乘法模型： $x_{t}=T_{t}*S_{t}*C_{t}*varepsilon _{t}$

4 AR模型： $x_{t}=phi _{0}+phi _{1}x_{t-1}+phi_ {2}x_{t-2}+...+phi_ {p}x_{t-p}+varepsilon _{t}$

以前p期的序列值 $x_{t-1}、x_{t-2}、 x_{t-p}$ 為自變數， $x_{t}$ 為其線性模型

5 MA模型： $x_{t}=mu +varepsilon _{t}- heta_{1}varepsilon _{t-1}- heta_{2}varepsilon_ {t-2}-...- heta_{q}varepsilon _{t-q}$

$x{t}$與以前各期的序列值無關，是前q期的隨機擾動項的線性回歸模型

6 ARMA模型： $x_{t}=phi _{0}+phi _{1}x_{t-1}+phi _{2}x_{t-2}+...+phi _{p}x_{t-p}+varepsilon _{t}- heta_{1}varepsilon _{t-1}- heta_{2}varepsilon _{t-2}-...- heta_{q}varepsilon _{t-q}$ $x_{t}$ 不僅與前各期序列值有關，還與隨機擾動項有關

7 ARIMA模型：許多非平穩序列差分後會顯示平穩序列的性質，稱這個非平穩為差分平穩序列。對查分平穩用ARIMA模型進行擬合

8 ARCH模型：ARCH模型能準確地模擬時間序列變數的波動性的變化，適用於序列具有異方差性質且異方差函數短期自相關

9 GARCH模型及衍生模型：GARCH模型稱為廣義ARCH模型，是ARCH模型的拓展。相比ARCH模型，GARCH模型及其衍生模型更能反映實際序列中的長期記憶性、信息的非對稱性等性質

2. 時間序列預處理

純隨機性和平穩性檢驗（時間序列預處理）

1 純隨機序列又稱白雜訊序列，序列的各項沒有任何關係。白雜訊序列是沒有信息可提取的平穩序列

2 對於平穩的非白雜訊序列，它的均值和方差是常數，ARMA是最常用的平穩序列擬合模型

3 非平穩序列，均值方差不穩定，一般將其轉為平穩序列，這樣可以用如ARMA。如果一個序列差分後具有平穩性，則該序列為差分平穩序列，可用ARIMA

2.1 平穩性檢驗

（1）平穩時間序列定義

如果時間序列 $x_{t}$ 在某一常數附件波動且範圍有限，即有常數均值方差，且延遲$k$期的序列變數的子自協方差和自相關係數是相等的，則 $x_{t}$ 是平穩序列。

自協方差和自相關係數衡量同一個事件在兩個不同時期時期（t,s）之間的相關程度，形象可以說是度量自己過去行為對自己現在行為的影響。

任取 $t,sin T$ 定義時間序列 $left { X_{t} ight }$ 的自協方差函數 $gamma (t,s))=E[(X_{t}-mu _{t})(X_{s}-mu _{s})]$ 和自相關係數 $ho (t,s))=frac{cov(X_{t},X_{s})}{sigma _{t}sigma _{s}}$（$gamma (t,t)=gamma (0)=1, ho _{0}=1$ ）

（2）平穩性檢驗

1）時序圖檢驗：根據平穩時間序列的均值和方差都為常數的性質，平穩時間序列的時序圖顯示該序列始終在一個常數值附近隨機波動，而且波動有界。

2）自相關檢驗：；平穩時間序列具有短期相關性，這表明平穩時間序列通常只有近期的序列值對現在的值影響較明顯。隨著延遲期數$k$的增加，平穩序列的自相關係數會衰減趨向於零，並在零附近隨機波動，而非平穩序列的的自相關係數衰減速度較慢。

3）單位根檢驗：如果存在單位根就是非平穩時間序列了

2.2 純隨機性檢驗

如果一個序列式是純隨機序列，滿足 $gamma (k)=0$ 。實際上純隨機序列的樣本自相關係數不會絕對為零，但很接近零並在零附近隨機擾動。

純隨機性檢驗又稱白雜訊檢驗，一般是構造檢驗統計量：Q統計量和LB統計量

3 平穩時間序列分析

3.1 AR模型

$AR(p)$ : $x_{t}=phi _{0}+phi _{1}x_{t-1}+phi _{2}x_{t-2}+...+phi _{p}x_{t-p}+varepsilon _{t}$ $

隨機干擾項 $varepsilon _{t}$ 為均值為零的白雜訊序列

均值方差為常數，ACF拖尾，PACF p階截尾

3.2 MA模型

$MA(q)$ : $x_{t}=mu +varepsilon _{t}- heta_{1}varepsilon _{t-1}- heta_{2}varepsilon _{t-2}-...- heta_{q}varepsilon _{t-q}$

隨機干擾項 $varepsilon _{t}$ 為均值為零的白雜訊序列， $mu$ 時序列均值

均值方差為常數，ACF q階截尾，PACF拖尾

3.3 ARMA模型

$ARMA(p,q)$ : $x_{t}=phi _{0}+phi _{1}x_{t-1}+phi _{2}x_{t-2}+...+phi _{p}x_{t-p}+varepsilon _{t}- heta_{1}varepsilon _{t-1}- heta_{2}varepsilon _{t-2}-...- heta_{q}varepsilon _{t-q}$

隨機干擾項 $varepsilon _{t}$ 為均值為零的白雜訊序列

均值方差為常數，ACF拖尾，PACF拖尾

3.4 平穩時間序列建模

平穩非白雜訊序列
計算ACF,PACF
ARMA模型識別【模型定階，根據ACF,PACF的拖尾與截尾性質】
估計模型中未知參數的值
模型檢驗（Y進入第六步，N返回第三部）
模型優化
預測將來趨勢

4. 非平穩時間序列分析

對非平穩時間的序列的分析方法可以分為確定性因素分解的時間序列分析和隨機時間序列分析的兩大類

確定性因素分解把時間序列變化歸結為4個因素（長期趨勢，季節變動，循環變動和隨機波動）的綜合影響。

隨機時序分析可以建立的模型有ARIMA模型、殘差自回歸模型、季節模型、異方差模型

4.1 差分運算

1） p階差分運算

相距一期的兩個序列值之間的減法運算稱為1階差分運算

2）k步差分

相距k期的兩個序列值之間的減法運算稱為k步差分運算

4.2 ARIMA模型

許多非平穩序列差分後會顯示出平穩序列的性質，這時這個非平穩序列為差分平穩序列。對差分平穩序列用ARIMA模型進行擬合。ARIMA模型的實質就是差分運算和ARMA模型的組合。

步驟：

1. 獲取觀察值序列

2. 平穩性檢驗（平穩則跳過第三步）

3. 差分運算

4. 白雜訊檢驗（純隨機性檢驗，Y退出分析，N進入下一步）

5. ARIMA擬合

案例：2015/1/1~2015/2/6某餐廳銷售數據進行建模

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as sns #seaborn畫出的圖更好看，且代碼更簡單，缺點是可塑性差from statsmodels.graphics.tsaplots import plot_acf #自相關圖from statsmodels.tsa.stattools import adfuller as ADF #平穩性檢測from statsmodels.graphics.tsaplots import plot_pacf #偏自相關圖from statsmodels.stats.diagnostic import acorr_ljungbox #白雜訊檢驗from statsmodels.tsa.arima_model import ARIMA#seaborn 是建立在matplotlib之上的%matplotlib inline %pylab inline#jupyter中文顯示是方框，加入下面兩行即可顯示中文，若嫌麻煩，也可去網上搜索如何永久顯示中文plt.rcParams[font.sans-serif] = [SimHei] #用來正常顯示中文標籤plt.rcParams[axes.unicode_minus] = False #用來正常顯示負號pylab.rcParams[figure.figsize] = (10, 6) #設置輸出圖片大小sns.set(color_codes=True) #seaborn設置背景#讀取數據，指定日期列為指標，Pandas自動將「日期」列識別為Datetime格式data = pd.read_excel(arima_data.xls, index_col = u日期)#時序圖data.plot()<matplotlib.axes._subplots.AxesSubplot at 0x1cea2ebceb8>

有明顯的遞增趨勢，可以判斷是非平穩的

#自相關圖plot_acf(data).show()

自相關圖顯示自相關係數長期大於零，說明時間序列有很強的相關性

#平穩性檢測print(u原始序列的ADF檢驗結果為：, ADF(data[u銷量]))#返回值依次為adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore原始序列的ADF檢驗結果為： (1.8137710150945274, 0.9983759421514264, 10, 26, {1%: -3.7112123008648155, 5%: -2.981246804733728, 10%: -2.6300945562130176}, 299.46989866024177)

單位根統計量對應的p的值顯著大於0.05，最終判斷該序列是非平穩序列的（非平穩不一定不是白雜訊）

#1階差分後的結果D_data = data.diff().dropna() #1階差分，丟棄na值D_data.columns = [u銷量差分] #更改列名D_data.plot() #時序圖<matplotlib.axes._subplots.AxesSubplot at 0x1cea2e92630>

plot_acf(D_data).show() #自相關圖

print(u差分序列的ADF檢驗結果為：, ADF(D_data[u銷量差分])) #平穩性檢測差分序列的ADF檢驗結果為： (-3.1560562366723532, 0.02267343544004886, 0, 35, {1%: -3.6327426647230316, 5%: -2.9485102040816327, 10%: -2.6130173469387756}, 287.5909090780334)

1階差分後時序圖在均值附近平穩的波動、自相關圖有較強的短期相關性、單位根檢驗p值小於0.05，所以1階差分後的序列是平穩序列

print(u差分序列的白雜訊檢驗結果為：, acorr_ljungbox(D_data, lags=1)) #返回統計量和p值差分序列的白雜訊檢驗結果為： (array([11.30402222]), array([0.00077339]))

輸出p值遠小於0.05，所以1階差分之後是平穩非白雜訊序列

擬合ARMA模型：一種人為識別，看差分後的ACFPACF，第二種方法是相對最有模型

plot_pacf(D_data).show()

一階差分後的自相關圖顯示出一階截尾，偏自相關顯示拖尾性，所以考慮MA(1)模型進行擬合1階差分後的序列，即對原序列建立ARIMA(0,1,1)模型

#相對最優模型data[u銷量] = data[u銷量].astype(float) #銷量轉為float類型#定階pmax = int(len(D_data)/10) #一般階數不超過length/10qmax = int(len(D_data)/10) #一般階數不超過length/10bic_matrix = [] #bic矩陣for p in range(pmax+1): tmp = [] for q in range(qmax+1): try: #存在部分報錯，所以用try來跳過報錯。 tmp.append(ARIMA(data, (p,1,q)).fit().bic) except: tmp.append(None) bic_matrix.append(tmp)bic_matrix = pd.DataFrame(bic_matrix) #從中可以找出最小值p,q = bic_matrix.stack().idxmin() #先用stack展平，然後用idxmin找出最小值位置。print(uBIC最小的p值和q值為：%s、%s %(p,q))BIC最小的p值和q值為：0、1model = ARIMA(data, (p,1,q)).fit() #建立ARIMA(0, 1, 1)模型model.summary2() #給出一份模型報告

model.forecast(5) #作為期5天的預測，返回預測結果、標準誤差、置信區間。(array([4873.9665477 , 4923.92261622, 4973.87868474, 5023.83475326, 5073.79082178]), array([ 73.08574262, 142.32680643, 187.54283213, 223.80283273, 254.95705912]), array([[4730.72112437, 5017.21197102], [4644.96720157, 5202.87803086], [4606.3014882 , 5341.45588128], [4585.18926146, 5462.48024505], [4574.0841683 , 5573.49747526]]))

python主要的時間序列模式演算法函數

Time Series analysis 官方文檔

Time Series analysis?

www.statsmodels.org

acf() 計算自相關係數 statsmodels.tsa.stattols

plt_acf() 畫自相關係數 statsmodels.graphics.tsaplots

pacf() 計算片相關係數 statsmodels.tsa.stattols

plot_acf() 畫偏相關係數 statsmodels.graphics.tsaplots

adfuller() 對觀測值序列進行單位根檢驗 statsmodels.tsa.stattols//ADF test

diff() 對觀測值進行查分計算 pandas對象自帶的方法

ARIMA() 創建一個ARIMA序列模型 statsmodels.tsa.arima_model

summary()&summaty2 給出一分ARIMA模型的報告 ARIMA模型對象自帶方法

aic/bic/hqic 計算ARIMA模型的AIC/BIC/HQIC指標值 ARIMA模型對象自帶方法

froecast() 應用構建的時間序列進行預測 ARIMA模型對象自帶方法

acorr_ljungbox() Ljung-Box檢驗，檢驗是否為白雜訊 statsmodels.stats.diagnostic

參考資料：

《python數據分析與挖掘實戰》張良均、王路、譚立雲、蘇劍林著