以前都是用SPSS/EVIEWS等來做的，現在用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

執行的時候會報錯：Undefined variable from import:，不過好像不影響結果

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)
dta.index = pd.Index(sm.tsa.datetools.dates_from_range('2001','2090')) //應該是2090，不是2100

dta.plot(figsize=(12,8))

plt.show() // 在Scala IDE要輸入這個命令才能顯示圖片

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)

已經平穩，二階效果也差不多，具體見原文

3.合適的p,q 


現在我們已經得到一個平穩的時間序列，接來下就是選擇合適的ARIMA模型，即ARIMA模型中合適的p,q 。 

 3.1 第一步我們要先檢查平穩時間序列的自相關圖和偏自相關圖。

diff1= dta.diff(1)#我們已經知道要使用一階差分的時間序列，之前判斷差分的程式可以註釋掉 //原文有錯誤應該是diff1= dta.diff(1)，而非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 

3.2模型選擇
根據上圖，猜測有以下模型可以供選擇：
1）ARMA(0,1)模型：即自相關圖在滯後1階之後縮小為0，且偏自相關縮小至0，則是一個階數q=1的移動平均模型；
2.）ARMA(7,0)模型：即偏自相關圖在滯後7階之後縮小為0，且自相關縮小至0，則是一個階層p=7的自迴歸模型； //原文錯寫為3
3.）ARMA(7,1)模型：即使得自相關和偏自相關都縮小至零。則是一個混合模型。
4） …還可以有其他供選擇的模型 （實際上下文選了ARMA(8,0)）
現在有以上這麼多可供選擇的模型，我們通常採用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()
print(arma_mod70.aic,arma_mod70.bic,arma_mod70.hqic)
arma_mod30 = sm.tsa.ARMA(dta,(0,1)).fit()
print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)
arma_mod71 = sm.tsa.ARMA(dta,(7,1)).fit()
print(arma_mod71.aic,arma_mod71.bic,arma_mod71.hqic)
arma_mod80 = sm.tsa.ARMA(dta,(8,0)).fit()
print(arma_mod80.aic,arma_mod80.bic,arma_mod80.hqic)

把原文的mod名字改了，以qp階數命名更直觀，不過我跑出來的結果和原文不太一樣

1619.19185018 1641.69013721 1628.26448199
1657.21729729 1664.71672631 1660.2415079
1605.68656094 1630.68465765 1615.76726295
1597.93598102 1622.93407772 1608.01668303

這樣的話應該是ARMA(8，0)模型擬合效果最好。

3.3 接下來檢驗下殘差序列：

resid = arma_mod80.resid //原文把這個變數賦值語句漏了，所以老是出錯


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)
plt.show()

殘差的ACF和PACF圖，可以看到序列殘差基本為白噪聲

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

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

結果=2.02332930932，所以殘差序列不存在自相關性。

3.4 觀察是否符合正態分佈

這裡使用QQ圖，它用於直觀驗證一組資料是否來自某個分佈，或者驗證某兩組資料是否來自同一（族）分佈。

print(stats.normaltest(resid))
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
#plt.show()

結果表明基本符合正態分佈

3.5 殘差序列Ljung-Box檢驗，也叫Q檢驗

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.014746   0.020229  0.886899
2.0  -0.045590   0.215795  0.897720
3.0   0.101588   1.197979  0.753489
4.0   0.063346   1.584316  0.811608
5.0   0.176146   4.606765  0.465727
6.0   0.004994   4.609224  0.594816
7.0  -0.209083   8.970252  0.254799
8.0   0.115683  10.321552  0.243178
9.0   0.020911  10.366252  0.321657
10.0 -0.220114  15.380865  0.118781
11.0 -0.050668  15.649938  0.154632
12.0 -0.031545  15.755573  0.202690
13.0 -0.055370  16.085244  0.244557
14.0  0.195340  20.242428  0.122682
15.0 -0.204417  24.855630  0.051916
16.0  0.034558  24.989255  0.070015
17.0  0.181463  28.724206  0.037155
18.0 -0.043332  28.940139  0.049116
19.0 -0.045327  29.179739  0.063210
20.0  0.044198  29.410802  0.079980
21.0  0.141917  31.827664  0.060944
22.0  0.118171  33.528034  0.054822
23.0  0.004488  33.530524  0.072258
24.0 -0.133612  35.770166  0.057769
25.0  0.061816  36.256927  0.067791
26.0 -0.021829  36.318576  0.085964
27.0  0.047270  36.612251  0.102536
28.0  0.131287  38.914114  0.082314
29.0 -0.080915  39.802824  0.087213
30.0 -0.026180  39.897408  0.106867
31.0  0.011463  39.915849  0.130949
32.0 -0.015224  39.948936  0.157836
33.0 -0.042675  40.213480  0.181101
34.0  0.006060  40.218910  0.214085
35.0 -0.015980  40.257355  0.248807
36.0  0.000413  40.257381  0.287349
37.0 -0.083807  41.354656  0.286211
38.0 -0.091967  42.701428  0.276142
39.0  0.002435  42.702391  0.315032
40.0 -0.071601  43.551379  0.322770

prob值均大於0.05，所以殘差序列不存在自相關性

4. 預測

predict_dta = arma_mod80.predict('2090', '2100', dynamic=True)
print(predict_dta)


fig, ax = plt.subplots(figsize=(12, 8))
ax = dta.ix['2000':].plot(ax=ax)
fig = arma_mod80.plot_predict('2090', '2100', dynamic=True, ax=ax, plot_insample=False)
plt.show()

結果：

2090-12-31     9541.934170
2091-12-31    12906.497016
2092-12-31    13979.448587
2093-12-31    14499.448647
2094-12-31    13892.217206
2095-12-31    13247.336277
2096-12-31    10958.837240
2097-12-31    10070.285452
2098-12-31    12680.448582
2099-12-31    13472.910370
2100-12-31    13611.974609

結尾

1.模型的好壞暫時還沒查到如何檢驗，比如APE,MAPE等等

2.系統還有2個報錯，雖然不影響結果，原因待查

D:\Python27\lib\site-packages\statsmodels\tsa\arima_model.py:1724: FutureWarning: TimeSeries is deprecated. Please use Series
  forecast = TimeSeries(forecast, index=self.data.predict_dates)
D:\Python27\lib\site-packages\matplotlib\collections.py:446: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == 'face':

3. 另外還參考了官方的example，http://statsmodels.sourceforge.net/devel/examples/notebooks/generated/tsa_arma_0.html