以前都是用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