[python] 時間序列分析之ARIMA

1 時間序列與時間序列分析

在生產和科學研究中，對某一個或者一組變量  進行觀察測量，將在一系列時刻  所得到的離散數字組成的序列集合，稱之為時間序列。 
時間序列分析是根據系統觀察得到的時間序列數據，通過曲線擬合和參數估計來建立數學模型的理論和方法。時間序列分析常用於國民宏觀經濟控制、市場潛力預測、氣象預測、農作物害蟲災害預報等各個方面。

2 時間序列建模基本步驟

1. 獲取被觀測系統時間序列數據；
2. 對數據繪圖，觀測是否為平穩時間序列；對於非平穩時間序列要先進行d階差分運算，化為平穩時間序列；
3. 經過第二步處理，已經得到平穩時間序列。要對平穩時間序列分別求得其自相關系數ACF 和偏自相關系數PACF ，通過對自相關圖和偏自相關圖的分析，得到最佳的階層 p 和階數 q
4. 由以上得到的 ，得到ARIMA模型。然后開始對得到的模型進行模型檢驗。

3 ARIMA實戰解剖

原理大概清楚，實踐卻還是會有諸多問題。相比較R語言，Python在做時間序列分析的資料相對少很多。下面就通過Python語言詳細解析后三個步驟的實現過程。 
文中使用到這些基礎庫： 。 對其調用如下

1.  
   from __future__ import print_function
2.  
   import pandas as pd
3.  
   import numpy as np
4.  
   from scipy import stats
5.  
   import matplotlib.pyplot as plt
6.  
   import statsmodels.api as sm
7.  
   from statsmodels.graphics.api import qqplot

3.1 獲取數據

這里我們使用一個具有周期性的測試數據，進行分析。 
數據如下：

dta=[10930,10318,10595,10972,7706,6756,9092,10551,9722,10913,11151,8186,6422, 
6337,11649,11652,10310,12043,7937,6476,9662,9570,9981,9331,9449,6773,6304,9355, 
10477,10148,10395,11261,8713,7299,10424,10795,11069,11602,11427,9095,7707,10767, 
12136,12812,12006,12528,10329,7818,11719,11683,12603,11495,13670,11337,10232, 
13261,13230,15535,16837,19598,14823,11622,19391,18177,19994,14723,15694,13248, 
9543,12872,13101,15053,12619,13749,10228,9725,14729,12518,14564,15085,14722, 
11999,9390,13481,14795,15845,15271,14686,11054,10395]

1.  
   dta=pd.Series(dta)
2.  
   dta. index = pd.Index(sm.tsa.datetools.dates_from_range('2001','2100'))
3.  
   dta.plot(figsize=( 12,8))

3.2 時間序列的差分

ARIMA 模型對時間序列的要求是平穩型。因此，當你得到一個非平穩的時間序列時，首先要做的即是做時間序列的差分，直到得到一個平穩時間序列。如果你對時間序列做d次差分才能得到一個平穩序列，那么可以使用ARIMA(p,d,q)模型，其中d是差分次數。

1.  
   fig = plt.figure(figsize=(12,8))
2.  
   ax1= fig.add_subplot(111)
3.  
   diff1 = dta.diff(1)
4.  
   diff1.plot(ax=ax1)

 
一階差分的時間序列的均值和方差已經基本平穩，不過我們還是可以比較一下二階差分的效果

1.  
   fig = plt.figure(figsize=(12,8))
2.  
   ax2= fig.add_subplot(111)
3.  
   diff2 = dta.diff(2)
4.  
   diff2.plot(ax=ax2)

 
可以看出二階差分后的時間序列與一階差分相差不大，並且二者隨着時間推移，時間序列的均值和方差保持不變。因此可以將差分次數d設置為1。 
其實還有針對平穩的檢驗，叫“ADF單位根平穩型檢驗”，以后再更。

3.3 合適的

現在我們已經得到一個平穩的時間序列，接來下就是選擇合適的ARIMA模型，即ARIMA模型中合適的。 
第一步我們要先檢查平穩時間序列的自相關圖和偏自相關圖。

1.  
   dta= dta.diff(1) #我們已經知道要使用一階差分的時間序列，之前判斷差分的程序可以注釋掉
2.  
   fig = plt.figure(figsize=(12,8))
3.  
   ax1=fig.add_subplot(211)
4.  
   fig = sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax1)
5.  
   ax2 = fig.add_subplot(212)
6.  
   fig = sm.graphics.tsa.plot_pacf(dta,lags=40,ax=ax2)

其中lags 表示滯后的階數，以上分別得到acf 圖和pacf 圖 
 
通過兩圖觀察得到： 
* 自相關圖顯示滯后有三個階超出了置信邊界； 
* 偏相關圖顯示在滯后1至7階（lags 1,2,…，7）時的偏自相關系數超出了置信邊界，從lag 7之后偏自相關系數值縮小至0 
則有以下模型可以供選擇： 
1. ARMA(0,1)模型：即自相關圖在滯后1階之后縮小為0，且偏自相關縮小至0，則是一個階數q=1的移動平均模型； 
2. ARMA(7,0)模型：即偏自相關圖在滯后7階之后縮小為0，且自相關縮小至0，則是一個階層p=3的自回歸模型； 
3. ARMA(7,1)模型：即使得自相關和偏自相關都縮小至零。則是一個混合模型。 
4. …還可以有其他供選擇的模型 
現在有以上這么多可供選擇的模型，我們通常采用ARMA模型的AIC法則。我們知道：增加自由參數的數目提高了擬合的優良性，AIC鼓勵數據擬合的優良性但是盡量避免出現過度擬合(Overfitting)的情況。所以優先考慮的模型應是AIC值最小的那一個。赤池信息准則的方法是尋找可以最好地解釋數據但包含最少自由參數的模型。不僅僅包括AIC准則，目前選擇模型常用如下准則： 
* AIC=-2 ln(L) + 2 k 中文名字：赤池信息量 akaike information criterion 
* BIC=-2 ln(L) + ln(n)*k 中文名字：貝葉斯信息量 bayesian information criterion 
* HQ=-2 ln(L) + ln(ln(n))*k hannan-quinn criterion 
構造這些統計量所遵循的統計思想是一致的，就是在考慮擬合殘差的同時，依自變量個數施加“懲罰”。但要注意的是，這些准則不能說明某一個模型的精確度，也即是說，對於三個模型Ａ，Ｂ，Ｃ，我們能夠判斷出Ｃ模型是最好的，但不能保證Ｃ模型能夠很好地刻畫數據，因為有可能三個模型都是糟糕的。

1.  
   arma_mod20 = sm.tsa.ARMA(dta,(7,0)).fit()
2.  
   print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)
3.  
   arma_mod30 = sm.tsa.ARMA(dta,(0,1)).fit()
4.  
   print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)
5.  
   arma_mod40 = sm.tsa.ARMA(dta,(7,1)).fit()
6.  
   print(arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic)
7.  
   arma_mod50 = sm.tsa.ARMA(dta,(8,0)).fit()
8.  
   print(arma_mod50.aic,arma_mod50.bic,arma_mod50.hqic)

 
可以看到ARMA(7,0)的aic，bic，hqic均最小，因此是最佳模型。

3.4 模型檢驗

在指數平滑模型下，觀察ARIMA模型的殘差是否是平均值為0且方差為常數的正態分布（服從零均值、方差不變的正態分布），同時也要觀察連續殘差是否（自）相關。

3.4.1 我們對ARMA(7,0)模型所產生的殘差做自相關圖

1.  
   fig = plt.figure(figsize=(12,8))
2.  
   ax1 = fig.add_subplot(211)
3.  
   fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
4.  
   ax2 = fig.add_subplot(212)
5.  
   fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)

3.4.2 做D-W檢驗

德賓-沃森（Durbin-Watson）檢驗。德賓-沃森檢驗,簡稱D-W檢驗，是目前檢驗自相關性最常用的方法，但它只使用於檢驗一階自相關性。因為自相關系數ρ的值介於-1和1之間，所以 0≤DW≤４。並且DW＝O＝＞ρ＝１　　 即存在正自相關性 
DW＝４＜＝＞ρ＝－１　即存在負自相關性 
DW＝２＜＝＞ρ＝０　　即不存在（一階）自相關性 
因此，當DW值顯著的接近於O或４時，則存在自相關性，而接近於２時，則不存在（一階）自相關性。這樣只要知道ＤＷ統計量的概率分布，在給定的顯著水平下，根據臨界值的位置就可以對原假設Ｈ０進行檢驗。

print(sm.stats.durbin_watson(arma_mod20.resid.values)) 

檢驗結果是2.02424743723，說明不存在自相關性。

3.4.3 觀察是否符合正態分布

這里使用QQ圖，它用於直觀驗證一組數據是否來自某個分布，或者驗證某兩組數據是否來自同一（族）分布。在教學和軟件中常用的是檢驗數據是否來自於正態分布。QQ圖細節，下次再更。

1.  
   resid = arma_mod20.resid #殘差
2.  
   fig = plt.figure(figsize=( 12,8))
3.  
   ax = fig.add_subplot( 111)
4.  
   fig = qqplot(resid, line= 'q', ax=ax, fit=True)

3.4.4 Ljung-Box檢驗

Ljung-Box test是對randomness的檢驗,或者說是對時間序列是否存在滯后相關的一種統計檢驗。對於滯后相關的檢驗，我們常常采用的方法還包括計算ACF和PCAF並觀察其圖像，但是無論是ACF還是PACF都僅僅考慮是否存在某一特定滯后階數的相關。LB檢驗則是基於一系列滯后階數，判斷序列總體的相關性或者說隨機性是否存在。 
時間序列中一個最基本的模型就是高斯白噪聲序列。而對於ARIMA模型，其殘差被假定為高斯白噪聲序列，所以當我們用ARIMA模型去擬合數據時，擬合后我們要對殘差的估計序列進行LB檢驗，判斷其是否是高斯白噪聲，如果不是，那么就說明ARIMA模型也許並不是一個適合樣本的模型。

1.  
   r,q,p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
2.  
   data = np.c_[range(1,41), r[1:], q, p]
3.  
   table = pd.DataFrame( data, columns=['lag', "AC", "Q", "Prob(>Q)"])
4.  
   print(table.set_index( 'lag'))

 
檢驗的結果就是看最后一列前十二行的檢驗概率（一般觀察滯后1~12階），如果檢驗概率小於給定的顯著性水平，比如0.05、0.10等就拒絕原假設，其原假設是相關系數為零。就結果來看，如果取顯著性水平為0.05，那么相關系數與零沒有顯著差異，即為白噪聲序列。

3.5 模型預測

模型確定之后，就可以開始進行預測了，我們對未來十年的數據進行預測。

1.  
   predict_sunspots = arma_mod20.predict( '2090', '2100', dynamic=True)
2.  
   print(predict_sunspots)
3.  
   fig, ax = plt.subplots(figsize=( 12, 8))
4.  
   ax = dta.ix[ '2001':].plot(ax=ax)
5.  
   predict_sunspots.plot(ax=ax)

 
前面90個數據為測試數據，最后10個為預測數據；從圖形來，預測結果較為合理。至此，本案例的時間序列分析也就結束了。

參考文獻與推薦閱讀

1. statsmodels–statistics in python
2. 時間序列分析—(ARIMA模型)
3. Arima預測模型（R語言）
4. 介紹QQplot
5. LBQ檢驗
6. 經管之家