前段時間整理了一個預測的基本思考架構和常見的方法，其中提到了ARIMA模型，在《大資料預測》那本書裡，ARIMA是單獨開辟一章講的，比較複雜和難了解的一個模型，自己最近找了點資料粗淺學習了一下理論，并嘗試用Python實作一下。
一、時間序列資料及其預處理
1.時序資料
時序資料顧名思義就是随着時間而變動的資料，是指某個個體在不同時間點上收集到的資料。已經被收集(或者叫觀察)到的資料其實是時間序列變量的一個觀察值，但由于時間的不可逆性，每一個時間點的變量有且僅能有一個觀察值，我們用這些觀察值拟合預測模型，用來預測未來時刻時間序列變量的值(此時，未發生的時間序列變量是一個随機變量)。
2.時序資料預處理
時序資料的預處理主要是平穩性和随機性檢驗，根據檢驗結果的不同，用不同的模型對資料進行模組化。 純随機序列，又稱為白噪聲序列，序列的變量之間無任何關系，是以沒有分析價值；平穩非白噪聲序列，其均值和方差是常數，目前有比較成熟的模型對其進行模拟，主要是AR(Autoregressive model)自回歸模型，MA(Movingaverage model)移動平均模型、以及ARMA-一個AR和MA結合的模型。
(1)平穩性檢驗
平穩性檢驗的基本原理是檢驗不同時間點時序變量之間的相關關系，由于是一個時間序列資料的不同時刻之間的互相比較，是以取名為自協方差函數和自相關系數，但基本原理和協方差函數、相關系數類似，當一個時間序列有為常數的均值和方差，且其自協方差函數隻跟時間間隔有關，即γ(t,s) =γ(k,k+s-t) ，其中γ是自協方差函數。 平穩性檢驗的方法分為兩類，一類是偏主觀和模糊的圖示檢驗，一類是比較精确的數學檢驗。 第一類主要是用時序圖和自相關圖進行觀察，如果一個時間序列的時序圖圍繞某一個常數上下波動，且波動範圍有界，可以認為是平穩的，除此之外，如果自相關圖中自相關系數比較快速衰減到0，并且在零附近震蕩，可以認為序列是平穩序列，這是因為平穩序列通常具有短期相關性，随着延遲期數的增加，變量之間的相關關系會指數降低。
第二類資料的檢驗方法是機關根檢驗。機關根又稱為機關圓，指模型的特征方程的根和1的關系，如果模型的特征根小于1，則序列非平穩，如果特征根大于1，則序列平穩。
(2)白噪聲檢驗
白噪聲序列的各個序列值之間沒有任何關系，即γ(k)=0，在實際中，隻要序列變量之間的協方差系數在0值附近波動，就可以認為是白噪聲序列了。
二、常用模型
平穩時間序列的模組化通常是AR、MA和ARMA模型，它們都是一種多元線性回歸模型。
1.AR模型
這個模型的定義是t時刻的序列取值和前p期的序列值有關，且是前p期的線性組合。這個比較容易了解，在通常的印象中，我們大都覺得這次發生的事和最近幾次發生的事有關聯。 并非所有的AR模型都是平穩的，AR模型滿足平穩性的條件是系數|φ|<1； 同時自相關系數ACF呈指數的速度衰減，在0附近波動，表現出拖尾性；其次，由于AR模型裡在計算t和t-k期之間的相關性時會受到中間k期的影響，剔除k期影響後的自相關系數稱為偏自相關系數PACF，平穩AR模型的PACF具有截尾性(尾巴突然斷掉)，這是因為yt隻跟y(t-1)和y(t-p)之間的序列值有關，當時間間隔k>p的時候，yt和y(t-k)之間就沒有直接關系了。
2.MA模型
這個模型剛開始不太好了解，某個時刻的序列值跟過去曆史q期的白噪聲有關？are you kidding me？後來看網上有人說MA模型關注的是AR模型中的随機擾動項的累加，通過擾動項的移動平均能夠有效消除預測中的随機波動。 MA模型的自相關系數在q階後全部等于0，具有截尾性，偏自相關系數在q階後拖尾，這是因為任何一個可逆的MA模型都是可以轉換為無窮階的AR模型(各種嵌套回去)，而無窮階的AR模型的PACF是拖尾的。
3.ARMA模型
ARMA模型是AR和MA模型的線性組合，特别的，當q=0，它是一個AR(p)模型，如果p=0，它是一個MA(q)模型。 平穩條件下的三類模型的自相關系數ACF和偏自相關系數PACF特性如下：
利用這些特性可以通過圖示的方法判斷平穩性。
三、模組化步驟
時間序列經過資料預處理被判定為平穩非白噪聲序列後，就可以利用ARMA進行模組化。 1.計算ACF、PACF； 2.ARMA模型定階，即确定最優p、q值； 3.估計模型參數、進行參數檢驗； 4.模型檢驗； 5.模型預測及評價 在這裡提及一下ARIMA模型，ARIMA是指時間序列經過I階差分之後可以滿足平穩非白噪聲序列的條件，并利用ARMA為其進行模組化。是以，ARIMA和ARMA最大的差異就是前者對原資料序列進行了差分，差分次數一般為1次，2次都不太推薦，因為2次差分之後很可能就成為白噪聲序列，别問我怎麼知道的(╯︵╰)。
四、Python試一試
會列出一些重點步驟和代碼。 1.把資料分成了data_time的模型内資料和用于判斷模型效果的data_time_test模型外資料
data_time_origin = pd.read_excel(r'C:\*\arima_data_gmv.xlsx',index_col = '日期')data_time = data_time_origin[:'2017-05-20'] #将5-20号之前的資料作為模型資料data_time_test = data_time_origin['2017-05-21':]
2.看一下時序圖和自相關函數圖判斷是否平穩序列
data_time.plot(figsize=(20,8),fontsize=15)#導入相應的統計模型包from statsmodels.graphics.tsaplots import plot_acf #acf計算from statsmodels.tsa.stattools import adfuller as ADF #機關根檢測plot_acf(data_time).show()
時序圖上能看到還是有比較明顯的趨勢的，初步判斷不是平穩序列。
ACF上看自相關系數長期大于0，且後期又長期小于0，呈現一種對稱模式，可以判斷不是平穩序列。機關根檢驗也證明不是平穩序列。
adf_test = ADF(data_time['gmv'])print(adf_test)'''統計值=臨界值不能拒絕零假設，存在機關根，序列非平穩。可以看到p-value=0.6377顯著大于0.05，另外統計量計算值-1.28顯著大于1%、5%、10%的臨界值，是以序列存在機關根，序列是非平穩的'''result:(-1.281137399884022, 0.6377571909319921, 5, 75, {'1%': -3.520713130074074, '5%': -2.9009249540740742, '10%': -2.5877813777777776}, 573.5874067738569)
3.進行一階差分運算看其是否能達到平穩
data_time2 = data_time.diff(periods=1)data_time2 = data_time2.dropna()data_time2.columns = ['gmv差分'] data_time2.plot(figsize=(20,8),fontsize = 15)
可以看到基本消除了趨勢，資料有穩定的均值。
#自相關圖plot_acf(data_time2).show()#偏自相關圖from statsmodels.graphics.tsaplots import plot_pacfplot_pacf(data_time2).show() #機關根檢驗adf_test_d = ADF(data_time2['gmv差分'])print(adf_test_d)result:(-8.458458355369716, 1.591087857743347e-13, 4, 75, {'1%': -3.520713130074074, '5%': -2.9009249540740742, '10%': -2.5877813777777776}, 565.8248717064124)
自相關圖能看出有一些短期的相關性，機關根檢驗的結果是統計量-8.45小于臨界值，p-value也遠遠小于0.05，差分之後是平穩序列。
#白噪聲檢驗from statsmodels.stats.diagnostic import acorr_ljungboxprint(acorr_ljungbox(data_time2, lags=[3]))result:(array([13.84640456]), array([0.00312186]))
白噪聲檢驗的p-value 0.003小于0.05，是以拒絕原假設，差分後的時序不是白噪聲序列。是以上面的檢驗過程說明差分後的時間序列是平穩非白噪聲序列，可以進行平穩時序模組化。
4.模型定階
上面的PACF圖，能看到詭異的拖尾性，前半段截尾，後半段拖尾...看着像是1階AR模型，具體是不是、行不行我也不知道。 暴力用BIC資訊值算一下，找最優的pq值。
from statsmodels.tsa.arima_model import ARIMAdata_time['gmv'] = data_time['gmv'].astype(float)#定階pmax = int(len(data_time)/10) #一般階數不超過length/10qmax = int(len(data_time)/10) #一般階數不超過length/10bic_matrix = [] #bic矩陣for p in range(pmax + 1): temp = [] for q in range(qmax + 1): try: temp.append(ARIMA(data_time,(p,1,q)).fit().bic) except: temp.append(None) bic_matrix.append(temp)bic_matrix = pd.DataFrame(bic_matrix)bic_matrix = bic_matrix.stack()p,q = bic_matrix.idxmin()print('BIC值最小的p和q分别為：p=%s, q=%s' %(p,q))result:BIC值最小的p和q分别為：p=3, q=3
是以定階結果竟然是3和3，那就試着用這個模拟一下看看。
5.建構最終模型
arima_model = ARIMA(data_time,(p,1,q)).fit() arima_model.summary()
強烈懷疑這個模型的有效性... wierd，誰能給指條明路啊。用了80天的資料做訓練，然後在p=3，q=3的時候，模型結果隻給了系數，沒給有效性的資訊，攤手。 6.既然模型參數沒法檢驗，隻能用預測結果看看效果
forecast_n = pd.DataFrame(arima_model.forecast(11)[0],index = pd.date_range('2017-05-21',periods = 11))forecast_n.columns = ['gmv_fore']forecast_nerror_test = pd.merge(data_time_test,forecast_n,left_index = True, right_index = True)error_test['誤差'] = error_test['gmv'] - error_test['gmv_fore']plt.figure(figsize=(15,6))plt.plot(error_test.index,error_test['gmv'],'-*',error_test.index,error_test['gmv_fore'],'-go')plt.legend(loc = 'best', fontsize=20)plt.show()print('預測值與真實值之間的絕對誤差為：', np.abs(sum(error_test['gmv']) - sum(error_test['gmv_fore'])))result:預測值與真實值之間的絕對誤差為：17.937944428429546
其中藍線是原值，綠線是預測值，看結果還行，把上升和下降的趨勢都已經預測出來了，用原值的和同預測值的和做差取絕對值，絕對值17.9，在10天的預測區間内，200多的基礎數量級上，這個誤差我個人覺得還挺nice啊，比我之前做的那些勞什子玩意強多了，而且還是到天啊，到天！以前都是月度的。
7.做一點模型評價
計量經濟學裡進行模型評價的名額有三種：一是均方誤差MSE、一個是絕對平均誤差MAE，另外一個是符号正确百分率。符号正确百分率基本用于經濟領域，或者說對增長和下降的趨勢預測，隻是方向對了就算預測不錯了，是以此處用前兩種。
'''1.均方誤差2.絕對平均誤差'''#均方根誤差&均方誤差rmse = np.sqrt(sum(error_test['誤差']**2)/len(error_test['誤差']))mse = sum(error_test['誤差']**2)/len(error_test['誤差'])print('均方根誤差為：%.2f' % rmse,'\n','均方誤差為：%.2f' % mse)#絕對平均誤差mae = sum(abs(error_test['誤差']))/len(error_test['誤差'])print('絕對平均誤差為：%.2f' % mae)result:均方根誤差為：18.69 均方誤差為：349.30絕對平均誤差為：13.10
當然這幾個名額是需要在不同模型之間比較才能有優劣關系的，也就是說選用不同的p，q，或者甚至是非ARIMA模型的其它模型預測出來的結果之間的互相比較，但鑒于我今天已經不想搞了，是以就不比較了。但無論怎樣，我覺得這個預測的結果還真是挺不錯的，對于天的預測竟然趨勢全對，而且在200的基數量級上最多是5-30日那天預測誤差49，預測多了49萬，其它時候模型的效果還是蠻好的，我很滿意，嗯，就是這樣！
參考資料： 1.中國大學MOOC-中南财經政法大學《時間序列分析》. 2.中國大學MOOC-對外經濟貿易大學《計量經濟學導論》. 3.張良均 等. Python資料分析與挖掘實戰[M].北京:機械工業出版社,2015:119-134. 4.網上衆多牛人文章.