前言：我用的是python3.6版本。這篇部落格參考了以下兩篇博文，根據自己的了解重新整合，稍作修改代碼，希望能夠幫到大家。如果在過程中發現錯誤或者有更好的方式實作，歡迎大家發表評論。
參考1: https://blog.csdn.net/u010414589/article/details/49622625
參考2: https://blog.csdn.net/hal_sakai/article/details/51965657
正文：
一、時間序列模組化基本步驟
1、擷取被觀測系統時間序列資料；
2、對資料繪圖，觀測是否為平穩時間序列；對于非平穩時間序列要先進行d階差分運算，化為平穩時間序列；
3、經過第二步處理，已經得到平穩時間序列。要對平穩時間序列分别求得其自相關系數ACF 和偏自相關系數PACF ，通過對自相關圖和偏自相關圖的分析，得到最佳的階層 p 和階數 q
4、由以上得到的d、q、p得到ARIMA模型。然後開始對得到的模型進行模型檢驗。
二、實戰解析
基礎庫： pandas；numpy；scipy；matplotlib；statsmodels ：
#-*- coding:utf-8 -*-
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
plt.rcParams['font.sans-serif']=['SimHei'] #用來畫圖時正常顯示中文标簽
plt.rcParams['axes.unicode_minus']=False #用來畫圖時正常顯示坐标軸負數
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf 2.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=np.array(dta,dtype=np.float) #轉換資料類型
dta=pd.Series(dta)
#print(type(dta))
dta.index=pd.Index(sm.tsa.datetools.dates_from_range('2001','2090'))
plt.plot(dta,label='原始資料')
plt.legend() #讓label生效 #plt.show() #PyCharm IDE要輸入這個指令才能顯示圖檔，由于下面還要畫圖， 是以我隻要在程式最後附上這句話即可 
2.2 差分d
ARIMA 模型對時間序列的要求是平穩型。是以，當你得到一個非平穩的時間序列時，首先要做的即是做時間序列的差分，直到得到一個平穩時間序列。如果你對時間序列做d次差分才能得到一個平穩序列，那麼可以使用ARIMA(p,d,q)模型，其中d是差分次數。
plt.figure() #新開一個圖窗
diff1=dta.diff(1)
plt.plot(diff1,label='一階')
diff2=dta.diff(2)
plt.plot(diff2,label='二階')
plt.legend() 
一階差分的時間序列的均值和方差已經基本平穩，不過我們還是可以比較一下二階差分的效果。可以看出二階差分後的時間序列與一階差分相差不大，并且二者随着時間推移，時間序列的均值和方差保持不變。是以可以将差分次數d設定為1。
其實還有針對平穩的檢驗，叫“ADF機關根平穩型檢驗”(還沒有學習這部分，等我學習了會更這部分内容的)。
2.3 合适的p,q
現在我們已經得到一個平穩的時間序列，接來下就是選擇合适的ARIMA模型，即ARIMA模型中合适的p,qp,q。
2.3.1 第一步我們要先檢查平穩時間序列的自相關圖和偏自相關圖。
其中lags 表示滞後的階數，以下分别得到差分後的平穩序列的acf 圖和pacf 圖。
因為差分後的序列第一個是NAN，是以要去掉它，直接對diff1運作圖形會運作不出來，有兩種方法，#後是第二種
fig=plt.figure()
ax1=fig.add_subplot(211) # 畫子圖
plot_acf(diff1[1:],lags=40,ax=ax1)
#diff11=np.diff(dta)
#plot_acf(diff11,lags=40,ax=ax1)
ax2=fig.add_subplot(212)
plot_pacf(diff1[1:],lags=40,ax=ax2) 通過兩圖觀察得到：
* 自相關圖顯示滞後有三個階超出了置信邊界；
* 偏相關圖顯示在滞後1至7階（lags 1,2,…，7）時的偏自相關系數超出了置信邊界，從lag 7之後偏自相關系數值縮小至0
2.3.2 模型選擇（p,q選擇）
這部分文字來自參考的第一篇博文，根據acf和pacf圖得到以下模型，但是我并不了解這部分，我不明白這些模型是如何通過觀察得到的，如果有知道的童鞋希望可以告知。
根據上圖，猜測有以下模型可以供選擇：
1. ARMA(0,1)模型：即自相關圖在滞後1階之後縮小為0，且偏自相關縮小至0，則是一個階數q=1的移動平均模型；
2. ARMA(7,0)模型：即偏自相關圖在滞後7階之後縮小為0，且自相關縮小至0，則是一個階層p=7的自回歸模型；
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_mod70 = sm.tsa.ARMA(dta,(7,0)).fit()
arma_mod01 = sm.tsa.ARMA(dta,(0,1)).fit()
arma_mod71 = sm.tsa.ARMA(dta,(7,1)).fit()
arma_mod80 = sm.tsa.ARMA(dta,(8,0)).fit()
print(arma_mod70.aic,arma_mod70.bic,arma_mod70.hqic)
print(arma_mod01.aic,arma_mod01.bic,arma_mod01.hqic)
print(arma_mod71.aic,arma_mod71.bic,arma_mod71.hqic)
print(arma_mod80.aic,arma_mod80.bic,arma_mod80.hqic) 1619.1917731179842 1641.69006015 1628.26440493
1657.2172636424214 1664.71669265 1660.24147424
1605.686576696957 1630.6846734 1615.7672787
1597.935980998883 1622.9340777 1608.01668301
可以看到ARMA(8,0)的aic，bic，hqic均最小，是以是最佳模型。
2.4 模型檢驗
在指數平滑模型下，觀察ARIMA模型的殘差是否是平均值為0且方差為常數的正态分布（服從零均值、方差不變的正态分布），同時也要觀察連續殘差是否（自）相關。
2.4.1 對ARMA(8,0)模型所産生的殘差做自相關圖
殘差的ACF和PACF圖，可以看到序列殘差基本為白噪聲。
2.4.2 做D-W檢驗
德賓-沃森（Durbin-Watson）檢驗。德賓-沃森檢驗,簡稱D-W檢驗，是目前檢驗自相關性最常用的方法，但它隻使用于檢驗一階自相關性。因為自相關系數ρ的值介于-1和1之間，是以 0≤DW≤４。并且DW＝O＝＞ρ＝１ 即存在正自相關性
DW＝４＜＝＞ρ＝－１ 即存在負自相關性
DW＝２＜＝＞ρ＝０ 即不存在（一階）自相關性
是以，當DW值顯著的接近于0或４時，則存在自相關性，而接近于２時，則不存在（一階）自相關性。這樣隻要知道ＤＷ統計量的機率分布，在給定的顯著水準下，根據臨界值的位置就可以對原假設Ｈ０進行檢驗。
print(sm.stats.durbin_watson(arma_mod80.resid.values)) 結果是 2.02323317584，說明不存在自相關性。
2.4.3 是否符合正态分布
這裡使用QQ圖，它用于直覺驗證一組資料是否來自某個分布，或者驗證某兩組資料是否來自同一（族）分布。在教學和軟體中常用的是檢驗資料是否來自于正态分布。（QQ圖還未學習，待更）
resid = arma_mod80.resid#殘差
qqplot(resid, line='q',fit=True) 2.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')) 結果：
AC Q Prob(>Q)
lag
1.0 -0.014699 0.020101 0.887255
2.0 -0.045553 0.215348 0.897920
3.0 0.101611 1.197985 0.753488
4.0 0.063415 1.585168 0.811455
5.0 0.176177 4.608680 0.465475
6.0 0.004990 4.611135 0.594563
7.0 -0.209066 8.971443 0.254713
8.0 0.115697 10.323082 0.243078
9.0 0.020907 10.367765 0.321541
10.0 -0.220087 15.381123 0.118772
11.0 -0.050643 15.649936 0.154632
12.0 -0.031516 15.755374 0.202699
13.0 -0.055326 16.084532 0.244595
14.0 0.195279 20.239100 0.122782
15.0 -0.204407 24.851871 0.051968
16.0 0.034586 24.985716 0.070078
17.0 0.181442 28.719772 0.037199
18.0 -0.043290 28.935286 0.049176
19.0 -0.045298 29.174578 0.063288
20.0 0.044176 29.405417 0.080078
21.0 0.141981 31.824456 0.060989
22.0 0.118182 33.525145 0.054859
23.0 0.004506 33.527654 0.072303
24.0 -0.133570 35.765894 0.057824
25.0 0.061832 36.252912 0.067849
26.0 -0.021823 36.314524 0.086035
27.0 0.047269 36.608186 0.102616
28.0 0.131266 38.909330 0.082392
29.0 -0.080900 39.797706 0.087299
30.0 -0.026185 39.892329 0.106967
31.0 0.011452 39.910736 0.131066
32.0 -0.015213 39.943774 0.157970
33.0 -0.042655 40.208071 0.181254
34.0 0.006051 40.213485 0.214256
35.0 -0.015970 40.251881 0.248995
36.0 0.000407 40.251907 0.287553
37.0 -0.083808 41.349205 0.286411
38.0 -0.091955 42.695609 0.276349
39.0 0.002442 42.696577 0.315252
40.0 -0.071594 43.545399 0.322996
檢驗的結果就是看最後一列前十二行的檢驗機率（一般觀察滞後1~12階），如果檢驗機率小于給定的顯著性水準，比如0.05、0.10等就拒絕原假設，其原假設是相關系數為零。就結果來看，如果取顯著性水準為0.05，那麼相關系數與零沒有顯著差異，即為白噪聲序列。即：prob值均大于0.05，是以殘差序列不存在自相關性。
2.5 預測結果
模型确定之後，就可以開始進行預測了，我們對未來十年的資料進行預測。
predict_sunspots = arma_mod80.predict('2090', '2100', dynamic=True)
print(predict_sunspots)
fig, ax = plt.subplots()
ax = dta.ix['2001':].plot(ax=ax)
predict_sunspots.plot(ax=ax) plt.show() 2090-12-31 9542.423714
2091-12-31 12907.339818
2092-12-31 13980.500757
2093-12-31 14500.625514
2094-12-31 13893.455315
2095-12-31 13248.652060
2096-12-31 10960.214775
2097-12-31 10072.055950
2098-12-31 12682.506171
2099-12-31 13475.109548
2100-12-31 13614.243276
Freq: A-DEC, dtype: float64
前面90個資料為測試資料，最後10個為預測資料。（還沒有預測模型好壞的标準）
三、
在模型選擇方面我用的是自動選擇BIC值最小的值：
from statsmodels.tsa.arima_model import ARIMA #定階:導入ARIMA
pmax=int(len(diff1)/10) #一般階數不超過length/10
qmax=int(len(diff1)/10)
bic_matrix=[] #BIC
for p in range(pmax+1):
tmp=[]
for q in range(qmax+1):
try: #處理報錯
tmp.append(ARIMA(dta,(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("BIC最小的p值和q值為：%s、%s"%(p,q))
model=ARIMA(dta,(p,1,q)).fit()
predict_sunspots = model.predict('2090', '2100', dynamic=True)
print(predict_sunspots)
fig, ax = plt.subplots()
ax = dta.ix['2001':].plot(ax=ax)
predict_sunspots.plot(ax=ax)
plt.show() 結果：BIC最小的p值和q值為：6、6
2090-12-31 -867.898294
2091-12-31 4643.628545
2092-12-31 775.860782
2093-12-31 357.981031
2094-12-31 -1137.189339
2095-12-31 -406.305865
2096-12-31 -2943.157268
2097-12-31 -654.484380
2098-12-31 4519.613997
2099-12-31 727.625300
2100-12-31 119.724101
Freq: A-DEC, dtype: float64
問題：這樣做出的模型 （6,6）效果還不如（8,0）的效果好，而且預測結果還出現了負數，希望有人能夠解答這個問題。