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

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

2 時間序列建模基本步驟

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

3 ARIMA實戰解剖

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

from __future__ import print_function
import pandas as pd
import numpy as np
from scipy import  stats
import matplotlib.pyplot as plt
import statsmodels.api as sm
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]

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

3.2 時間序列的差分d

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

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

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

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

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

3.3 合適的p,q

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

dta= dta.diff(1)#我們已經知道要使用一階差分的時間序列，之前判斷差分的程式可以註釋掉
fig = plt.figure(figsize=(12,8))
ax1=fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax1)
ax2 = fig.add_subplot(212)
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
構造這些統計量所遵循的統計思想是一致的，就是在考慮擬合殘差的同時，依自變數個數施加“懲罰”。但要注意的是，這些準則不能說明某一個模型的精確度，也即是說，對於三個模型Ａ，Ｂ，Ｃ，我們能夠判斷出Ｃ模型是最好的，但不能保證Ｃ模型能夠很好地刻畫資料，因為有可能三個模型都是糟糕的。

arma_mod20 = sm.tsa.ARMA(dta,(7,0)).fit()
print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)
arma_mod30 = sm.tsa.ARMA(dta,(0,1)).fit()
print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)
arma_mod40 = sm.tsa.ARMA(dta,(7,1)).fit()
print(arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic)
arma_mod50 = sm.tsa.ARMA(dta,(8,0)).fit()
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)模型所產生的殘差做自相關圖

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
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圖細節，下次再更。

resid = arma_mod20.resid#殘差
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
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模型也許並不是一個適合樣本的模型。

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

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

3.5 模型預測

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

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

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

參考文獻與推薦閱讀