ARIMA模型預(yù)測(cè)

時(shí)間序列分析預(yù)測(cè)就是在已有的和時(shí)間有關(guān)的數(shù)據(jù)序列的基礎(chǔ)上構(gòu)建其數(shù)據(jù)模型并預(yù)測(cè)其未來(lái)的數(shù)據(jù)，例如航空公司的一年內(nèi)每日乘客數(shù)量、某個(gè)地區(qū)的人流量，這些數(shù)據(jù)往往具有周期性的規(guī)律。

如下圖所示，有的數(shù)據(jù)呈現(xiàn)出簡(jiǎn)單的周期性循環(huán)，有的呈現(xiàn)出周期性循環(huán)變化。

ARIMA（Autoregressive Integrated Moving Average model，差分整合移動(dòng)平均自回歸模型）不僅可以很好地對(duì)已有數(shù)據(jù)進(jìn)行擬合，并且還可以在此基礎(chǔ)上對(duì)未來(lái)數(shù)據(jù)進(jìn)行預(yù)測(cè)。

其函數(shù)原型為：ARIMA(data,order=(p,d,q))，因此除了需要傳入原始數(shù)據(jù)data之外，還需要三個(gè)參數(shù)p、d、q，這三個(gè)參數(shù)與ARIMA的原理有關(guān)。

ARIMA可以拆分為AR、I、MA，分別對(duì)應(yīng)p、d、q三個(gè)參數(shù)。

1.AR（AutoRegressive，自回歸模型）描述當(dāng)前值與歷史值之間的關(guān)系，用變量自身的歷史時(shí)間數(shù)據(jù)對(duì)自身進(jìn)行預(yù)測(cè)，自回歸過(guò)程的公式定義：yt是當(dāng)前值 u是常數(shù)項(xiàng) P是階數(shù) ri是自相關(guān)系數(shù) et是誤差。因此參數(shù)p為自回歸模型的階數(shù)。

2.MA（Moving Average，移動(dòng)平均值）利用一定時(shí)間間隔內(nèi)的平均值作為某一期的估計(jì)值，這樣可以消除數(shù)據(jù)的周期性影響。其公式為。因此參數(shù)q為移動(dòng)平均的階數(shù)。

3.I代表差分操作，時(shí)間序列最常用來(lái)剔除周期性因素的方法，它主要是對(duì)等周期間隔的數(shù)據(jù)進(jìn)行線性求減。從而使數(shù)據(jù)變得平穩(wěn)，ARIMA一般進(jìn)行一次差分即可穩(wěn)定，因此d一般取值為0、1、2，這里取1，進(jìn)行一次差分。

確定參數(shù)值

根據(jù)BIC

因此要使用ARIMA模型進(jìn)行預(yù)測(cè)首先需要確定參數(shù)p、q的值。因?yàn)橐话汶A數(shù)不超過(guò)整體數(shù)據(jù)的十分之一，因此這里采用窮舉法。

分別從0~10取p、q，并且得到該模型的BIC值（貝葉斯信息量 bayesian information criterion）。

BIC準(zhǔn)則綜合考慮了殘差大小和自變量的個(gè)數(shù)，殘差越小、自變量個(gè)數(shù)越少BIC值越小，模型越優(yōu)。

因此從所有的p、q取值中找出bic值最小的即為理想的模型參數(shù)。經(jīng)過(guò)如下函數(shù)，我得到p、q為5、5

def get_order(data):
    pmax = int(len(data) / 10)    #一般階數(shù)不超過(guò) length /10
    qmax = int(len(data) / 10)
    bic_matrix = []
    for p in range(pmax +1):
        temp= []
        for q in range(qmax+1):
            try:
                temp.append(ARIMA(data, order=(p, 1, q)).fit().bic)    # 將bic值存入二維數(shù)組
            except:
                temp.append(None)
        bic_matrix.append(temp)    
    bic_matrix = pd.DataFrame(bic_matrix)   #將其轉(zhuǎn)換成Dataframe 數(shù)據(jù)結(jié)構(gòu)
    p,q = bic_matrix.stack().astype('float64').idxmin()        # 找出bic值最小對(duì)應(yīng)的索引
    return p,q

根據(jù)圖像

還可以根據(jù)自相關(guān)系數(shù)ACF與偏自相關(guān)系數(shù)PACF圖像來(lái)確定階數(shù)p、q的值，若下所示使用庫(kù)函數(shù)輸出數(shù)據(jù)的圖像

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
 
def draw_acf_pacf(data,lags):
    f = plt.figure(facecolor='white')
    ax1 = f.add_subplot(211)
    plot_acf(data,ax=ax1,lags=lags)
    ax2 = f.add_subplot(212)
    plot_pacf(data,ax=ax2,lags=lags)
    plt.subplots_adjust(hspace=0.5)
    plt.show()

繪制圖像有如下兩種類型，拖尾是指序列以指數(shù)率單調(diào)遞減或震蕩衰減，如下左圖所示。

截尾指序列從某個(gè)時(shí)點(diǎn)變得非常小，如下右圖所示。

根據(jù)ACF圖像與PACF圖像的拖尾結(jié)尾選取不同的模型：

• 自相關(guān)系數(shù)拖尾，則q=0，偏自相關(guān)系數(shù)p階截尾：AR(p)模型，
• 自相關(guān)系數(shù)q階截尾，偏自相關(guān)函數(shù)拖尾，p=0：MA(q)模型
• 自相關(guān)函數(shù)和偏自相關(guān)函數(shù)均拖尾，階數(shù)分別為q、p：ARMA(p,q)模型

 那么圖像如何確定階數(shù)呢？即觀察圖片從第幾個(gè)數(shù)據(jù)之后收斂到置信區(qū)間，例如上面右圖從24之后再無(wú)超出置信區(qū)間的數(shù)據(jù)，所以認(rèn)為其24階截尾。

進(jìn)行預(yù)測(cè)

接下來(lái)就可以使用arima模型進(jìn)行模型擬合與預(yù)測(cè)了，這里使用的是python第三方包statsmodels.tsa.arima.model中的ARIMA模型。這是Statsmodels自從0.11版本新獨(dú)立的模塊，其原來(lái)的模塊為statsmodels.tsa.arima_model.ARIMA，二者在功能上都實(shí)現(xiàn)了arima模型，并且具有相同的屬性和方法名，其返回值均為ARIMAResults對(duì)象，通過(guò)該對(duì)象的predict()、forecast()得到預(yù)測(cè)值。值得注意的是新舊模塊的方法的傳入?yún)?shù)和返回值類型不太相同，因此使用時(shí)需要注意。

其官方文檔連接為：statsmodels.tsa.arima.model.ARIMAResults、

這里使用predict()函數(shù)進(jìn)行預(yù)測(cè)，其可以傳入start、end參數(shù)代表預(yù)測(cè)數(shù)據(jù)的開(kāi)始和結(jié)束坐標(biāo)，坐標(biāo)值不僅可以是int、str，還可以是 datetime時(shí)間類型，如果start、end介于原有數(shù)據(jù)區(qū)域內(nèi)，即為對(duì)原數(shù)據(jù)的預(yù)測(cè)擬合，如果end超過(guò)了原有數(shù)據(jù)長(zhǎng)度，即代表對(duì)未來(lái)數(shù)據(jù)進(jìn)行預(yù)測(cè)。

forecast()函數(shù)接收一個(gè)steps參數(shù)，代表對(duì)未來(lái)多少個(gè)數(shù)據(jù)進(jìn)行預(yù)測(cè)，也可以傳入datetime時(shí)間，代表從樣本數(shù)據(jù)結(jié)束一直預(yù)測(cè)到該時(shí)間。

特別注意的是在使用pandas的Datetime時(shí)間索引的數(shù)據(jù)進(jìn)行訓(xùn)練和預(yù)測(cè)時(shí)，要保證時(shí)間間隔是有規(guī)律的，否則在預(yù)測(cè)時(shí)索引會(huì)出錯(cuò)而報(bào)錯(cuò)如下：

KeyError: 'The `end` argument could not be matched to a location related to the index of the data.'

import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
 
%matplotlib inline  
plt.rcParams['font.sans-serif'] = ['SimHei']    # 設(shè)置字符集顯示中文
plt.rcParams['axes.unicode_minus'] = False      # 設(shè)置負(fù)號(hào)正確顯示
 
# 顯示數(shù)據(jù)趨勢(shì)圖
def draw_data(data):
    plt.figure()
    plt.plot(data)
    
with h5py.File('./taxi.h5', 'r') as hf:
    # 讀取數(shù)據(jù)
    in_data=np.array(hf['in_data'])
    in_area = in_data.transpose(2, 0, 1)  # 調(diào)整維度
    in_mean_day = in_area.mean(axis=2)
    sub_data=pd.Series(in_mean_day[116][0:100])      # 截取數(shù)據(jù)集中的100個(gè)作為樣本
    draw_data(sub_data)
    
    # 確定模型階數(shù)
    # p,q=get_order(sub_data)
    # print('階數(shù)：',p,q) 
    
    # 利用ARIMA模型進(jìn)行預(yù)測(cè)
    model=ARIMA(sub_data,order=(5,1,5)).fit()        # 傳入?yún)?shù)，構(gòu)建并擬合模型
    predict_data=model.predict(0,150)                # 預(yù)測(cè)數(shù)據(jù)
    forecast=model.forecast(21)                      # 預(yù)測(cè)未來(lái)數(shù)據(jù)
    
    # 繪制原數(shù)據(jù)和預(yù)測(cè)數(shù)據(jù)對(duì)比圖
    plt.plot(sub_data,label='原數(shù)據(jù)')
    plt.plot(predict_data,label='預(yù)測(cè)數(shù)據(jù)')
    plt.legend()
    plt.show()
    

如下所示為擬合預(yù)測(cè)結(jié)果：

數(shù)據(jù)預(yù)處理

穩(wěn)定性

ARIMA所需要的數(shù)據(jù)是穩(wěn)定的，這樣才能進(jìn)行規(guī)律的發(fā)掘和擬合，穩(wěn)定的數(shù)據(jù)是指其均值和方差為一個(gè)常量、自協(xié)方差與時(shí)間無(wú)關(guān)，如下所示分別為穩(wěn)定的數(shù)據(jù)和均值不穩(wěn)定、方差不穩(wěn)定、自協(xié)方差隨時(shí)間變化

數(shù)據(jù)的穩(wěn)定性可以利用單位根檢驗(yàn)Dickey-Fuller Test 進(jìn)行判斷，即在一定置信水平下，假設(shè)時(shí)序數(shù)據(jù)Null hypothesis: 非穩(wěn)定、序列具有單位根。

因此對(duì)于一個(gè)平穩(wěn)的時(shí)序數(shù)據(jù)，就需要在給定的置信水平上顯著，拒絕原假設(shè)。

如下通過(guò)statsmodels包中的adfuller來(lái)進(jìn)行檢測(cè)：

from statsmodels.tsa.stattools import adfuller
 
# 利用ADF檢測(cè)數(shù)據(jù)是否平穩(wěn)
def check_stable(data):
    adf_res=adfuller(data)
    print('t統(tǒng)計(jì)量',adf_res[0])
    print('t統(tǒng)計(jì)量的P值',adf_res[1])
    print('延遲階數(shù)',adf_res[2])
    print('觀測(cè)值的個(gè)數(shù)',adf_res[3])
    for key,value in adf_res[4].items():
        print('臨界區(qū)間：',key,'，值：',value)

統(tǒng)計(jì)結(jié)果如下：其中P值低于0.05，可以拒絕原假設(shè)，即數(shù)據(jù)是平穩(wěn)的

移動(dòng)平均

那么如何使數(shù)據(jù)變得平穩(wěn)，首先可以使用數(shù)據(jù)平滑技術(shù)消除數(shù)據(jù)的周期性起伏變化，常用的平滑技術(shù)有移動(dòng)平均、加權(quán)移動(dòng)平均，移動(dòng)平均即利用一定時(shí)間間隔內(nèi)的平均值作為某一期的估計(jì)值，而指數(shù)平均則是用變權(quán)的方法來(lái)計(jì)算均值。通過(guò)pandas的rolling()方法可以實(shí)現(xiàn)數(shù)據(jù)移動(dòng)，ewa()實(shí)現(xiàn)加權(quán)移動(dòng)，mean()實(shí)現(xiàn)平均。step平均的間隔。

# 對(duì)數(shù)據(jù)進(jìn)行移動(dòng)平均來(lái)平滑數(shù)據(jù)
def move_average(data,step):
    series=pd.Series(data)
    rol_mean=series.rolling(window=step).mean()     # 移動(dòng)平均
    rol_weight_mean=pd.DataFrame.ewm(series,span=step).mean()       # 加權(quán)移動(dòng)平均
    
    return rol_mean,rol_weight_mean

差分

通過(guò)將相同間隔的數(shù)據(jù)做差來(lái)消除周期性的因素。pandas的diff可以完成對(duì)序列相隔steps作差，例如當(dāng)steps=7時(shí)，前7個(gè)數(shù)據(jù)是沒(méi)法和前面的數(shù)據(jù)作差的，因此前七個(gè)數(shù)據(jù)為Nan，通過(guò)dropna()刪除。

在預(yù)測(cè)操作之后，還需要將差分的數(shù)據(jù)還原回來(lái)。差分操作只是簡(jiǎn)單的相減，因此還原就是相加，將原數(shù)據(jù)通過(guò)shift向后移動(dòng)7個(gè)位置然后相加，即可實(shí)現(xiàn)還原操作。

# 進(jìn)行差分
diff=series.diff(steps)
diff=diff.dropna()
 
# 差分還原
diff_shift=sub_data.shift(7)
diff_recover=predict_data.add(diff_shift)

總結(jié)

以上為個(gè)人經(jīng)驗(yàn)，希望能給大家一個(gè)參考，也希望大家多多支持腳本之家。