ARIMA模型預測

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

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

ARIMA（Autoregressive Integrated Moving Average model，差分整合移動平均自回歸模型）不僅可以很好地對已有數(shù)據(jù)進行擬合，并且還可以在此基礎上對未來數(shù)據(jù)進行預測。

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

ARIMA可以拆分為AR、I、MA，分別對應p、d、q三個參數(shù)。

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

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

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

確定參數(shù)值

根據(jù)BIC

因此要使用ARIMA模型進行預測首先需要確定參數(shù)p、q的值。因為一般階數(shù)不超過整體數(shù)據(jù)的十分之一，因此這里采用窮舉法。

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

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

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

def get_order(data):
    pmax = int(len(data) / 10)    #一般階數(shù)不超過 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值最小對應的索引
    return p,q

根據(jù)圖像

還可以根據(jù)自相關系數(shù)ACF與偏自相關系數(shù)PACF圖像來確定階數(shù)p、q的值，若下所示使用庫函數(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)遞減或震蕩衰減，如下左圖所示。

截尾指序列從某個時點變得非常小，如下右圖所示。

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

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

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

進行預測

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

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

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

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

特別注意的是在使用pandas的Datetime時間索引的數(shù)據(jù)進行訓練和預測時，要保證時間間隔是有規(guī)律的，否則在預測時索引會出錯而報錯如下：

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']    # 設置字符集顯示中文
plt.rcParams['axes.unicode_minus'] = False      # 設置負號正確顯示
 
# 顯示數(shù)據(jù)趨勢圖
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個作為樣本
    draw_data(sub_data)
    
    # 確定模型階數(shù)
    # p,q=get_order(sub_data)
    # print('階數(shù)：',p,q) 
    
    # 利用ARIMA模型進行預測
    model=ARIMA(sub_data,order=(5,1,5)).fit()        # 傳入?yún)?shù)，構(gòu)建并擬合模型
    predict_data=model.predict(0,150)                # 預測數(shù)據(jù)
    forecast=model.forecast(21)                      # 預測未來數(shù)據(jù)
    
    # 繪制原數(shù)據(jù)和預測數(shù)據(jù)對比圖
    plt.plot(sub_data,label='原數(shù)據(jù)')
    plt.plot(predict_data,label='預測數(shù)據(jù)')
    plt.legend()
    plt.show()
    

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

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

穩(wěn)定性

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

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

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

如下通過statsmodels包中的adfuller來進行檢測：

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

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

移動平均

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

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

差分

通過將相同間隔的數(shù)據(jù)做差來消除周期性的因素。pandas的diff可以完成對序列相隔steps作差，例如當steps=7時，前7個數(shù)據(jù)是沒法和前面的數(shù)據(jù)作差的，因此前七個數(shù)據(jù)為Nan，通過dropna()刪除。

在預測操作之后，還需要將差分的數(shù)據(jù)還原回來。差分操作只是簡單的相減，因此還原就是相加，將原數(shù)據(jù)通過shift向后移動7個位置然后相加，即可實現(xiàn)還原操作。

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

總結(jié)

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

最新評論