前段時(shí)間整理了一個(gè)預(yù)測(cè)的基本思考框架和常見(jiàn)的方法，其中提到了ARIMA模型，在《大數(shù)據(jù)預(yù)測(cè)》那本書(shū)里，ARIMA是單獨(dú)開(kāi)辟一章講的，比較復(fù)雜和難理解的一個(gè)模型，自己最近找了點(diǎn)資料粗淺學(xué)習(xí)了一下理論，并嘗試用實(shí)現(xiàn)一下。
一、時(shí)間序列數(shù)據(jù)及其預(yù)處理
1.時(shí)序數(shù)據(jù)
時(shí)序數(shù)據(jù)顧名思義就是隨著時(shí)間而變動(dòng)的數(shù)據(jù)，是指某個(gè)個(gè)體在不同時(shí)間點(diǎn)上收集到的數(shù)據(jù)。已經(jīng)被收集（或者叫觀察）到的數(shù)據(jù)其實(shí)是時(shí)間序列變量的一個(gè)觀察值，但由于時(shí)間的不可逆性，每一個(gè)時(shí)間點(diǎn)的變量有且僅能有一個(gè)觀察值，我們用這些觀察值擬合預(yù)測(cè)模型，用來(lái)預(yù)測(cè)未來(lái)時(shí)刻時(shí)間序列變量的值（此時(shí)，未發(fā)生的時(shí)間序列變量是一個(gè)隨機(jī)變量）。
2.時(shí)序數(shù)據(jù)預(yù)處理
時(shí)序數(shù)據(jù)的預(yù)處理主要是平穩(wěn)性和隨機(jī)性檢驗(yàn)，根據(jù)檢驗(yàn)結(jié)果的不同，用不同的模型對(duì)數(shù)據(jù)進(jìn)行建模。
純隨機(jī)序列，又稱為白噪聲序列，序列的變量之間無(wú)任何關(guān)系，因此沒(méi)有分析價(jià)值；平穩(wěn)非白噪聲序列，其均值和方差是常數(shù)，目前有比較成熟的模型對(duì)其進(jìn)行模擬，主要是AR( model)自回歸模型，MA( model)移動(dòng)平均模型、以及ARMA-一個(gè)AR和MA結(jié)合的模型。
（1）平穩(wěn)性檢驗(yàn)
平穩(wěn)性檢驗(yàn)的基本原理是檢驗(yàn)不同時(shí)間點(diǎn)時(shí)序變量之間的相關(guān)關(guān)系，由于是一個(gè)時(shí)間序列數(shù)據(jù)的不同時(shí)刻之間的相互比較，因此取名為自協(xié)方差函數(shù)和自相關(guān)系數(shù)，但基本原理和協(xié)方差函數(shù)、相關(guān)系數(shù)類似，當(dāng)一個(gè)時(shí)間序列有為常數(shù)的均值和方差，且其自協(xié)方差函數(shù)只跟時(shí)間間隔有關(guān)，即γ(t,s) =γ(k,k+s-t)，其中γ是自協(xié)方差函數(shù)。平穩(wěn)性檢驗(yàn)的方法分為兩類，一類是偏主觀和模糊的圖示檢驗(yàn)，一類是比較精確的數(shù)學(xué)檢驗(yàn)。第一類主要是用時(shí)序圖和自相關(guān)圖進(jìn)行觀察，如果一個(gè)時(shí)間序列的時(shí)序圖圍繞某一個(gè)常數(shù)上下波動(dòng)，且波動(dòng)范圍有界，可以認(rèn)為是平穩(wěn)的，除此之外，如果自相關(guān)圖中自相關(guān)系數(shù)比較快速衰減到0，并且在零附近震蕩，可以認(rèn)為序列是平穩(wěn)序列，這是因?yàn)槠椒€(wěn)序列通常具有短期相關(guān)性，隨著延遲期數(shù)的增加，變量之間的相關(guān)關(guān)系會(huì)指數(shù)降低。
第二類數(shù)據(jù)的檢驗(yàn)方法是單位根檢驗(yàn)。單位根又稱為單位圓，指模型的特征方程的根和1的關(guān)系，如果模型的特征根小于1，則序列非平穩(wěn)，如果特征根大于1，則序列平穩(wěn)。
（2）白噪聲檢驗(yàn)
白噪聲序列的各個(gè)序列值之間沒(méi)有任何關(guān)系，即γ(k)=0，在實(shí)際中，只要序列變量之間的協(xié)方差系數(shù)在0值附近波動(dòng)，就可以認(rèn)為是白噪聲序列了。
二、常用模型
平穩(wěn)時(shí)間序列的建模通常是AR、MA和ARMA模型，它們都是一種多元線性回歸模型。
1.AR模型
這個(gè)模型的定義是t時(shí)刻的序列取值和前p期的序列值有關(guān)，且是前p期的線性組合。這個(gè)比較容易理解，在通常的印象中，我們大都覺(jué)得這次發(fā)生的事和最近幾次發(fā)生的事有關(guān)聯(lián)。
并非所有的AR模型都是平穩(wěn)的，AR模型滿足平穩(wěn)性的條件是系數(shù)|φ|同時(shí)自相關(guān)系數(shù)ACF呈指數(shù)的速度衰減，在0附近波動(dòng)時(shí)間序列模型是什么，表現(xiàn)出拖尾性；其次，由于AR模型里在計(jì)算t和t-k期之間的相關(guān)性時(shí)會(huì)受到中間k期的影響，剔除k期影響后的自相關(guān)系數(shù)稱為偏自相關(guān)系數(shù)PACF，平穩(wěn)AR模型的PACF具有截尾性（尾巴突然斷掉），這是因?yàn)閥t只跟y(t-1)和y(t-p)之間的序列值有關(guān)，當(dāng)時(shí)間間隔k>p的時(shí)候，yt和y(t-k)之間就沒(méi)有直接關(guān)系了。
2.MA模型
這個(gè)模型剛開(kāi)始不太好理解，某個(gè)時(shí)刻的序列值跟過(guò)去歷史q期的白噪聲有關(guān)？are you me？后來(lái)看網(wǎng)上有人說(shuō)MA模型關(guān)注的是AR模型中的隨機(jī)擾動(dòng)項(xiàng)的累加，通過(guò)擾動(dòng)項(xiàng)的移動(dòng)平均能夠有效消除預(yù)測(cè)中的隨機(jī)波動(dòng)。MA模型的自相關(guān)系數(shù)在q階后全部等于0，具有截尾性，偏自相關(guān)系數(shù)在q階后拖尾，這是因?yàn)槿魏我粋€(gè)可逆的MA模型都是可以轉(zhuǎn)換為無(wú)窮階的AR模型（各種嵌套回去），而無(wú)窮階的AR模型的PACF是拖尾的。
3.ARMA模型
ARMA模型是AR和MA模型的線性組合，特別的，當(dāng)q=0，它是一個(gè)AR(p)模型，如果p=0，它是一個(gè)MA(q)模型。平穩(wěn)條件下的三類模型的自相關(guān)系數(shù)ACF和偏自相關(guān)系數(shù)PACF特性如下：
利用這些特性可以通過(guò)圖示的方法判斷平穩(wěn)性。
三、建模步驟
時(shí)間序列經(jīng)過(guò)數(shù)據(jù)預(yù)處理被判定為平穩(wěn)非白噪聲序列后，就可以利用ARMA進(jìn)行建模。1.計(jì)算ACF、PACF；
2.ARMA模型定階，即確定最優(yōu)p、q值；3.估計(jì)模型參數(shù)、進(jìn)行參數(shù)檢驗(yàn)；
4.模型檢驗(yàn)；
5.模型預(yù)測(cè)及評(píng)價(jià)
在這里提及一下ARIMA模型，ARIMA是指時(shí)間序列經(jīng)過(guò)I階差分之后可以滿足平穩(wěn)非白噪聲序列的條件，并利用ARMA為其進(jìn)行建模。因此，ARIMA和ARMA最大的差異就是前者對(duì)原數(shù)據(jù)序列進(jìn)行了差分，差分次數(shù)一般為1次，2次都不太推薦，因?yàn)?次差分之后很可能就成為白噪聲序列，別問(wèn)我怎么知道的(╯︵╰)。
四、試一試
會(huì)列出一些重點(diǎn)步驟和代碼。
1.把數(shù)據(jù)分成了的模型內(nèi)數(shù)據(jù)和用于判斷模型效果的模型外數(shù)據(jù)
data_time_origin?=?pd.read_excel(r'C:\*\arima_data_gmv.xlsx',index_col?=?'日期')data_time?=?data_time_origin[:'2017-05-20']??#將5-20號(hào)之前的數(shù)據(jù)作為模型數(shù)據(jù)data_time_test?=?data_time_origin['2017-05-21':]
2.看一下時(shí)序圖和自相關(guān)函數(shù)圖判斷是否平穩(wěn)序列
data_time.plot(figsize=(20,8),fontsize=15)#導(dǎo)入相應(yīng)的統(tǒng)計(jì)模型包from statsmodels.graphics.tsaplots import plot_acf #acf計(jì)算from?statsmodels.tsa.stattools?import?adfuller?as?ADF?#單位根檢測(cè)plot_acf(data_time).show()
時(shí)序圖上能看到還是有比較明顯的趨勢(shì)的，初步判斷不是平穩(wěn)序列。
ACF上看自相關(guān)系數(shù)長(zhǎng)期大于0，且后期又長(zhǎng)期小于0，呈現(xiàn)一種對(duì)稱模式，可以判斷不是平穩(wěn)序列。單位根檢驗(yàn)也證實(shí)不是平穩(wěn)序列。
adf_test = ADF(data_time['gmv'])print(adf_test)'''統(tǒng)計(jì)值<臨界值，拒絕零假設(shè)，平穩(wěn)或趨勢(shì)平穩(wěn)，計(jì)算值>=臨界值不能拒絕零假設(shè)，存在單位根，序列非平穩(wěn)。可以看到p-value=0.6377顯著大于0.05，另外統(tǒng)計(jì)量計(jì)算值-1.28顯著大于1%、5%、10%的臨界值，因此序列存在單位根，序列是非平穩(wěn)的'''result:(-1.281137399884022,?0.6377571909319921,?5,?75,?{'1%':?-3.520713130074074,?'5%':?-2.9009249540740742,?'10%':?-2.5877813777777776},?573.5874067738569)
3.進(jìn)行一階差分運(yùn)算看其是否能達(dá)到平穩(wěn)
data_time2 = data_time.diff(periods=1)data_time2 = data_time2.dropna()data_time2.columns = ['gmv差分']data_time2.plot(figsize=(20,8),fontsize = 15)
可以看到基本消除了趨勢(shì)，數(shù)據(jù)有穩(wěn)定的均值。
#自相關(guān)圖plot_acf(data_time2).show()#偏自相關(guān)圖from statsmodels.graphics.tsaplots import plot_pacfplot_pacf(data_time2).show()?#單位根檢驗(yàn)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)
自相關(guān)圖能看出有一些短期的相關(guān)性，單位根檢驗(yàn)的結(jié)果是統(tǒng)計(jì)量-8.45小于臨界值，p-value也遠(yuǎn)遠(yuǎn)小于0.05，差分之后是平穩(wěn)序列。
#白噪聲檢驗(yàn)from statsmodels.stats.diagnostic import acorr_ljungboxprint(acorr_ljungbox(data_time2, lags=[3]))result:(array([13.84640456]),?array([0.00312186]))
白噪聲檢驗(yàn)的p-value 0.003小于0.05，因此拒絕原假設(shè)，差分后的時(shí)序不是白噪聲序列。因此上面的檢驗(yàn)過(guò)程說(shuō)明差分后的時(shí)間序列是平穩(wěn)非白噪聲序列，可以進(jìn)行平穩(wěn)時(shí)序建模。
4.模型定階
上面的PACF圖，能看到詭異的拖尾性，前半段截尾，后半段拖尾...看著像是1階AR模型，具體是不是、行不行我也不知道。
暴力用BIC信息值算一下，找最優(yōu)的pq值。
from statsmodels.tsa.arima_model import ARIMAdata_time['gmv'] = data_time['gmv'].astype(float)#定階pmax = int(len(data_time)/10) #一般階數(shù)不超過(guò)length/10qmax = int(len(data_time)/10) #一般階數(shù)不超過(guò)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
所以定階結(jié)果竟然是3和3，那就試著用這個(gè)模擬一下看看。
5.構(gòu)建最終模型
arima_model = ARIMA(data_time,(p,1,q)).fit()arima_model.summary()
強(qiáng)烈懷疑這個(gè)模型的有效性... wierd，誰(shuí)能給指條明路啊。用了80天的數(shù)據(jù)做訓(xùn)練，然后在p=3，q=3的時(shí)候，模型結(jié)果只給了系數(shù)，沒(méi)給有效性的信息，攤手。
6.既然模型參數(shù)沒(méi)法檢驗(yàn)，只能用預(yù)測(cè)結(jié)果看看效果
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('預(yù)測(cè)值與真實(shí)值之間的絕對(duì)誤差為：', np.abs(sum(error_test['gmv']) - sum(error_test['gmv_fore'])))result:預(yù)測(cè)值與真實(shí)值之間的絕對(duì)誤差為：17.937944428429546
其中藍(lán)線是原值，綠線是預(yù)測(cè)值，看結(jié)果還行，把上升和下降的趨勢(shì)都已經(jīng)預(yù)測(cè)出來(lái)了，用原值的和同預(yù)測(cè)值的和做差取絕對(duì)值，絕對(duì)值17.9時(shí)間序列模型是什么，在10天的預(yù)測(cè)區(qū)間內(nèi)，200多的基礎(chǔ)數(shù)量級(jí)上，這個(gè)誤差我個(gè)人覺(jué)得還挺nice啊，比我之前做的那些勞什子玩意強(qiáng)多了，而且還是到天啊，到天！以前都是月度的。
7.做一點(diǎn)模型評(píng)價(jià)
計(jì)量經(jīng)濟(jì)學(xué)里進(jìn)行模型評(píng)價(jià)的指標(biāo)有三種：一是均方誤差MSE、一個(gè)是絕對(duì)平均誤差MAE，另外一個(gè)是符號(hào)正確百分率。符號(hào)正確百分率基本用于經(jīng)濟(jì)領(lǐng)域，或者說(shuō)對(duì)增長(zhǎng)和下降的趨勢(shì)預(yù)測(cè)，只是方向?qū)α司退泐A(yù)測(cè)不錯(cuò)了，所以此處用前兩種。
'''1.均方誤差2.絕對(duì)平均誤差'''#均方根誤差&均方誤差rmse = np.sqrt(sum(error_test['誤差']**2)/len(error_test['誤差']))mse = sum(error_test['誤差']**2)/len(error_test['誤差'])print('均方根誤差為：%.2f'?% rmse,'\n','均方誤差為：%.2f'?% mse)#絕對(duì)平均誤差mae = sum(abs(error_test['誤差']))/len(error_test['誤差'])print('絕對(duì)平均誤差為：%.2f' % mae)result:均方根誤差為：18.69均方誤差為：349.30絕對(duì)平均誤差為：13.10
當(dāng)然這幾個(gè)指標(biāo)是需要在不同模型之間比較才能有優(yōu)劣關(guān)系的，也就是說(shuō)選用不同的p，q，或者甚至是非ARIMA模型的其它模型預(yù)測(cè)出來(lái)的結(jié)果之間的相互比較，但鑒于我今天已經(jīng)不想搞了，所以就不比較了。但無(wú)論怎樣，我覺(jué)得這個(gè)預(yù)測(cè)的結(jié)果還真是挺不錯(cuò)的，對(duì)于天的預(yù)測(cè)竟然趨勢(shì)全對(duì)，而且在200的基數(shù)量級(jí)上最多是5-30日那天預(yù)測(cè)誤差49，預(yù)測(cè)多了49萬(wàn)，其它時(shí)候模型的效果還是蠻好的，我很滿意，嗯，就是這樣！
參考資料：
1.中國(guó)大學(xué)MOOC-中南財(cái)經(jīng)政法大學(xué)《時(shí)間序列分析》.
2.中國(guó)大學(xué)MOOC-對(duì)外經(jīng)濟(jì)貿(mào)易大學(xué)《計(jì)量經(jīng)濟(jì)學(xué)導(dǎo)論》.3.張良均 等. 數(shù)據(jù)分析與挖掘?qū)崙?zhàn)[M].北京:機(jī)械工業(yè)出版社,2015:119-134.4.網(wǎng)上眾多牛人文章.