基于ARIMA和LSTM混合模型的時間序列預測

王英偉 馬樹才

(遼寧大學經濟學院 遼寧 沈陽 110036)

0 引 言

時間序列預測在眾多領域有廣泛應用，如金融、經濟、工程和航空等，并成為機器學習領域的重要研究課題[1]。現實中時間序列通常同時具有線性和非線性特征，因此對不同時間序列現象建立模型并準確預測已成為最具挑戰的應用之一。

ARIMA等統計模型以簡單和靈活性被大量用于時間序列預測[2-3]，但現實中的時間序列通常具有非線性特征，而傳統時間序列預測方法是線性模型，在對時間序列建模中表現出一定局限性[4]。因此，支持向量機和神經網絡等非線性模型，被廣泛用于時間序列預測領域[5]。其中神經網絡具有強大的非線性、映射性和自適應性等特征，可以有效提高預測精度，以LSTM為代表的深度神經網絡成為近年來的研究熱點。但單一非線性模型對同時具有線性和非線性特征的時間序列并不能獲得最優結果[6]。

由于單一模型的局限性，眾多學者提出由線性和非線性模型組成的混合模型。Zhang[7]和Oliveira等[1]假設時間序列同時具有線性和非線性特征，并分別提出ARIMA-ANN混合模型和ARIMA-SVR混合模型，ARIMA用于提取線性特征，而ANN和SVR用于提取殘差特征，最后對兩者進行組合，后者應用PSO選擇模型參數。文獻[8]提出ARIMA-SVR_s混合模型，首先將時間序列分解為高波動率和低波動率兩部分，并分別使用ARIMA和AR_SVR對兩部分建模，最后將兩者結果組合。實驗結果證明，混合模型可以有效提高預測精度。以上模型均假設時間序列是線性特征和非線性特征的線性組合,但在實際應用中，兩者可能存在非線性關系，從而影響混合模型的性能。

針對上述問題，本文對線性和非線性模型預測結果的組合方式提出新的策略。首先采用ARIMA模型提取時間序列線性特征，然后用SVR模型對ARIMA模型的預測值和實際值間的誤差序列進行預測，最后將前兩者預測結果作為LSTM模型的輸入對時間序列進行預測，并應用貝葉斯優化算法選擇深度LSTM模型的超參數。

1 ARIMA模型

ARIMA模型是由Box和Jenkins提出的用于時間序列分析和預測的線性模型。在時間序列分析中，需要設置3個參數，分別為自回歸階數(p)、差分階數(d)和移動平均階數(q)，ARIMA(p,d,q)的一般形式如下：

(1)

φ(B)=1-φ1B-φ2B2-…-φpBp

(2)

θ(B)=1-θ1B-θ2B2-…-θqBq

(3)

式中：B為后移算子；▽d=(1-B)d為高階差分；φi(i=1,2,…,p)和θj(j=1,2,…,q)分別為自回歸參數和移動平均參數；εt為符合N(0,σ2)正態分布的誤差項。式(2)和式(3)分別為自回歸相關系數多項式和移動平均系數多項式。

時間序列yt應先采用ADF檢驗選取差分階數，將時間序列變為平穩序列，再根據赤池信息準則(Akaike Information Criterion，AIC)選擇最佳模型。

2 支持向量回歸

支持向量回歸是基于結構風險最小化(SRM)原則進行泛化誤差上界最小化，最初由文獻[9]提出的機器學習方法。

f(x)=ωφ(x)+b

(4)

式中：ω為權值向量;φ是將樣本空間映射到高維空間的非線性變換函數;b為偏差。

為保證支持向量回歸具有稀疏性，引入ε不敏感損失函數，允許樣本值落入不敏感損失帶，即允許最大為ε的誤差，公式如下：

(5)

參數ω和b通過式(6)求得。

(6)

s.t.yi-ωφ(x)-b≤ε+ξ

ωφ(x)+b-yi≤ε+ξ*

ξ,ξ*≥0,i=0,1,…,l

式中：‖w‖2為模型復雜度項;C為懲罰因子，用于平衡模型復雜度和最小訓練誤差;ξ和ξ*表示松弛變量，用以度量不敏感損失帶外的訓練樣本偏離程度。將式(6)引入拉格朗日函數和核函數，則非線性回歸函數表示為：

(7)

3 LSTM模型

循環神經網絡(RNN)是用于處理序列數據的神經網絡。其最大特點是神經元在某一時刻的輸出能夠作為輸入再次輸入到神經元，從而使網絡具有記憶能力，但RNN無法解決長期時間序列依賴問題。因此，Hochreiter等[10]提出長短期記憶(LSTM)神經網絡。

LSTM模型是標準RNN的一種變體，通過引入記憶單元(Memory Cell)解決長期依賴問題，即梯度消失和梯度爆炸[11]，記憶單元結構如圖1所示。

圖1 LSTM記憶單元圖

每個記憶單元包含三個門控結構：遺忘門(ft)、輸入門(it)和輸出門(ot)，用于對類似傳送帶的單元狀態(Cell State)進行移除和添加信息。每個門控結構由一個sigmoid神經網絡層和一個點乘操作組成。假設xt表示t時刻的輸入向量，Wf、Ws、Wi和Wo表示循環層權重矩陣，Uf、Us、Ui和Uo表示輸入層權重矩陣，bf、bs、bi和bo表示偏置向量，ht表示LSTM層的輸出向量，計算過程如下：

1) 通過遺忘門移除單元狀態St-1中的無用歷史信息，并計算遺忘門激活值：

ft=sigmoid(Ufxt+Wfht-1+bf)

(8)

2) 通過輸入門決定將哪些新信息存儲到單元狀態：

it=sigmoid(Uixt+Wiht-1+bi)

(9)

(10)

3) 將舊單元狀態St-1更新為新單元狀態St：

(11)

4) 輸出門通過sigmoid層決定需要輸出的單元狀態，并將輸出結果Ot和經過tanh層的新單元狀態tanh(St)相乘計算記憶單元的輸出ht：

Ot=sigmoid(Uoxt+Woht-1+bo)

(12)

ht=Ot?tanh(St)

(13)

LSTM神經網絡要求輸入連續時間步上的特征向量值以便進行訓練。假設每個序列共有N個時間步，并且SN表示第N個時間步的特征向量,所以第M個序列可以表示為{SM,SM+1,…,SM+N-1},輸入特征向量序列的構造方式如圖2所示。

圖2 LSTM神經網絡的輸入特征向量序列構造方式

4 深度LSTM神經網絡

深度LSTM神經網絡是由多個基本LSTM模塊組成，核心思想是在較低LSTM層建立輸入數據的局部特征，并在較高LSTM層進行整合。實驗證明，深度LSTM神經網絡結構的學習能力更強[12]。

圖3 深度LSTM網絡結構

5 貝葉斯優化算法

貝葉斯優化算法是一種高效的全局優化算法[13]。它能夠在每次迭代中，根據代理模型擬合實際目標函數的結果選擇最優評估點，減少目標函數的迭代次數，尤其適合對評估代價高昂的黑盒目標函數進行優化，因此被廣泛用于機器學習和深度學習的超參數優化。

假設目標函數f:X→R，貝葉斯優化采用基于模型的序貫優化法對式(14)中的最優化問題求解。

(14)

貝葉斯優化算法的迭代過程可以分為三個部分。(1) 根據最大化采集函數選擇最優評估點，即xn+1=argmaxap(x)，其中ap(f):X→R為采集函數，p(f)為f的先驗概率分布；(2) 評估目標函數yn+1=f(xn+1)+ε，并將新的點(xn+1,yn+1)加入觀測數據Dn=(xj,yj),j=1，2，…，n；(3) 更新f的后驗概率分布和采集函數，分別表示為p(f|Dn+1)和ap(f|Dn+1)。

本文采用高斯過程(GP)作為代理模型，而采集函數則使用期望提高(Expected Improvement, EI)函數。

6 ARIMA_DLSTM模型

圖4 ARIMA_DLSTM模型流程

7 實 驗

7.1 實驗數據與預處理

本文采用5組來自不同領域的時間序列數據集進行實證研究，分別為Canadian Lynx、Colorado River flow、Airline Passengers、比特幣兌美元匯率和國內糖期貨價格指數。以上具有不同統計特征的數據集在時間序列預測研究的文獻中被廣泛應用[7,14-15]。

首先對Canadian Lynx數據集進行對數化處理，即將所有值取常用對數，再對每個數據集按4∶1的比例劃分為訓練集和測試集，訓練集用于訓練模型，測試集用于評估模型性能。Canadian Lynx數據集共有114個樣本，其中前100個樣本作為訓練集，后14個樣本作為測試集。Colorado River flow數據集共有744個樣本，前595個樣本作為訓練集，后149個樣本作為測試集。Airline Passengers數據集共有144個樣本，前115個樣本作為訓練集，后29個樣本作為測試集。比特幣兌美元匯率采用2014年12月至2019年7月間的數據，共有1 672個樣本，前1 472個樣本作為訓練集，后200個樣本作為測試集。國內糖期貨價格指數采用2006年3月至2019年4月間共3 190個樣本，前3 000個樣本作為訓練集，后190個樣本作為測試集。

在神經網絡訓練中，數據間的量綱差別對網絡訓練是否能收斂以及預測準確性起到關鍵作用[16]。因此，在建模前需要對輸入數據進行預處理，本文利用式(15)將樣本數據歸一化至[0,1]區間。

(15)

式中：minI(t)和maxI(t)分別為訓練數據集的最小值和最大值。訓練后的輸出數據進行反歸一化處理，以便得到預測值。

為提高模型泛化能力，采用如下兩種方式避免過擬合。首先，采用失活(dropout)正則化，在網絡訓練的每次更新中，以一定失活率隨機選擇一部分單元失活，包括輸入連接和遞歸連接，可以有效防止過擬合。如果采用深度LSTM網絡，可以同時在每層間采用失活正則化，所以每層網絡共有三個失活參數。其次采用早停法，將訓練樣本劃分為訓練集和驗證集，在每次迭代中，分別計算訓練集和驗證集的損失值，如果驗證集損失值在步數k內不再減少，則停止訓練，并返回最低驗證損失值的模型參數。本文中使用兩層LSTM神經網絡，因此共有六個失活參數。在訓練集中，80%樣本作為訓練集，20%樣本作為驗證集，步數k設置為50。

7.2 衡量指標

為評估ARIMA_DLSTM模型的預測性能，本文采用均方誤差(MSE)、平均絕對百分比誤差(MAPE)和平均絕對誤差(MAE)作為衡量預測精度的指標。計算公式分別為：

(16)

(17)

(18)

式中：Xt表示實際值;Ft表示預測值;N為時間序列數據集的樣本數量。

7.3 實驗結果

表1中列出ARIMA_DLSTM混合模型中SVR模型和DLSTM模型的參數及其選擇范圍。SVR模型采用RBF核函數，超參數包括懲罰系數C、不敏感損失系數ε、寬度系數γ和時間步timestep，可以通過網格搜索進行選擇，時間步數為輸入序列的長度。DLSTM模型的層數為2，超參數包括每層神經元數量和dropout參數，其中dropout參數包括輸入連接、遞歸連接和每層連接之間失活參數，所以2層網絡共有6個參數，并采用貝葉斯優化算法[17]進行超參數選擇，算法迭代次數為50。為評估ARIMA_DLSTM混合模型的性能，將ARIMA[7]、ARIMAMLP[18]、ARIMA_MLP[7]、ARIMA_SVR[8]和ARIMA_SVR_s[1]預測模型作為對比模型。

表1 模型參數及取值范圍

(1) Canadian Lynx數據集。在SVR模型中，懲罰系數C為1 000，不敏感損失系數ε為0.1，寬度系數γ為1.0，時間步(輸入序列長度)為10。LSTM模型的每層神經元個數分別為26和21，dropout為[0.10,0.26,0.16,0.12,0.27,0.26]，時間步為5，迭代次數為2 000次。

圖5給出利用貝葉斯優化算法對LSTM模型進行超參數優化的最優驗證誤差迭代圖。可以看出，算法在第42次迭代后誤差值最優。

圖5 貝葉斯優化算法在Canadian Lynx的驗證誤差迭代圖

表2給出了不同模型在Canadian Lynx時間序列測試集中的預測結果。可見，ARIMA_SVR模型在比較模型中最優，而ARIMA_DLSTM模型的測試結果比ARIMA_SVR模型的MSE、MAPE和MAE值分別降低53.47%、24.84%和26.80%。

表2 對Canadian Lynx時間序列預測結果

圖6給出ARIMA模型和ARIMA_DLSTM模型在Canadian Lynx測試集上的預測結果。由圖6可知，ARIMA_DLSTM模型的預測結果和真實值更加接近，有效提高了ARIMA模型的性能。

圖6 ARIMA模型和ARIMA_DLSTM模型在Canadian Lynx的時序列測結果

(2) Colorado River flow數據集。在SVR模型中，懲罰系數C為1 000，不敏感損失系數ε為0.001，寬度系數γ為0.01，時間步(輸入序列長度)為24。LSTM模型的每層神經元個數分別為24和34，dropout為[0.14,0.32,0.22,0.18,0.28,0.18]，時間步為43，迭代次數為2 000次。

圖7給出利用貝葉斯優化算法對LSTM模型進行超參數優化的最優驗證誤差迭代圖。可見，算法在第30次迭代后誤差值最優。

圖7 貝葉斯優化算法在Colorado River flow的驗證誤差迭代圖

表3給出了不同模型在Colorado River flow時間序列測試集中的預測結果。可見，ARIMA_MLP模型的MSE和MAE值在比較模型中最優，ARIMA_SVR_s模型的MAPE值最優，而ARIMA_DLSTM模型的測試結果比最優MSE、MAPE和MAE值分別降低82.15%、41.58%和38.26%。

表3 對Colorado River flow時間序列預測結果

圖8給出ARIMA模型和ARIMA_DLSTM模型在Colorado River flow測試集上的預測結果。可以看出，ARIMA_DLSTM模型預測值存在滯后現象，但和真實值更接近，優于ARIMA模型。

圖8 ARIMA模型和ARIMA_DLSTM模型在Colorado River flow的時序預測結果

(3) Airline Passengers數據集。在SVR模型中，懲罰系數C為1 000，不敏感損失系數ε為0.01，寬度系數γ為0.01，時間步(輸入序列長度)為16。LSTM模型的每層神經元個數分別為35和22，dropout為[0.14,0.27,0.11,0.10,0.27,0.25]，時間步為47，迭代次數為2 000次。

圖9給出利用貝葉斯優化算法對LSTM模型進行超參數優化的最優驗證誤差迭代圖。可以看出，算法在第41次迭代后誤差值最優。

圖9 貝葉斯優化算法在Airline Passengers的驗證誤差迭代圖

表4給出了不同模型在Airline Passengers時間序列測試集中的預測結果。可以看出，ARIMA_SVR_s模型預測結果在比較模型中最優，而ARIMA_DLSTM模型的測試結果比最優MSE、MAPE和MAE值分別降低23.61%、8.09%和5.72%。

表4 對Airline Passengers時間序列預測結果

圖10給出ARIMA模型和ARIMA_DLSTM模型在Airline Passengers測試集上的預測結果。可以看出，ARIMA_DLSTM模型預測值和真實值間趨勢基本吻合，有效提高了ARIMA模型的預測精度。

圖10 ARIMA模型和ARIMA_DLSTM模型在Airline Passengers 時間序列預測結果

(4) 比特幣匯率數據集。在SVR模型中，懲罰系數C為1 000，不敏感損失系數ε為0.001，寬度系數γ為0.01，時間步(輸入序列長度)為7。LSTM模型的每層神經元個數分別為23和22，dropout為[0.23,0.11,0.17,0.16,0.32,0.16]，時間步為40，迭代次數為2 000次。

圖11給出利用貝葉斯優化算法對LSTM模型進行超參數優化的最優驗證誤差迭代圖。可以看出，算法在第16次迭代后誤差值最優。

圖11 貝葉斯優化算法在Bitcoin匯率的驗證誤差迭代圖

表5給出了不同模型在Bitcoin匯率時間序列測試集中的預測結果。可以看出，ARIMA_SVR_s模型預測結果在比較模型中最優，而ARIMA_DLSTM模型的測試結果比最優MSE、MAPE和MAE值分別降低17.8%、17.13%和10.4%。

表5 對Bitcoin匯率時間序列預測結果

圖12給出ARIMA模型和ARIMA_DLSTM模型在Bitcoin匯率測試集上后20組數據的預測結果。可以看出，ARIMA_DLSTM模型預測值和真實值之間雖有滯后現象，但趨勢吻合，能夠提高ARIMA模型的預測精度。

圖12 ARIMA模型和ARIMA_DLSTM模型在Bitcoin匯率時間序列預測結果

(5) 糖價格指數數據集。在SVR模型中，懲罰系數C為1 000，不敏感損失系數ε為0.001，寬度系數γ為1，時間步(輸入序列長度)為3。LSTM模型的每層神經元個數分別為27和5，dropout為[0.03,0.18,0.19,0.06,0.01,0.3]，時間步為20，迭代次數為2 000次。

圖13給出利用貝葉斯優化算法對LSTM模型進行超參數優化的最優驗證誤差迭代圖。可以看出，算法在第25次迭代后誤差值最優。

圖13 貝葉斯優化算法在糖價格指數的驗證誤差迭代圖

表6給出了不同模型在糖價格指數時間序列測試集中的預測結果。可以看出，ARIMA_SVR_s模型預測結果在比較模型中最優，而ARIMA_DLSTM模型的測試結果比最優MSE、MAPE和MAE值分別降低12.48%、3.17%和3.88%。

表6 對糖價格指數時間序列預測結果

圖14給出ARIMA模型和ARIMA_DLSTM模型在糖價格指數測試集上后20組數據的預測結果。可以看出，ARIMA_DLSTM模型預測值和真實值更接近，精度更高。

圖14 ARIMA模型和ARIMA_DLSTM模型在該時間序列預測結果

圖15給出了5個數據集下的ARIMA模型和5種混合模型的MSE誤差比雷達圖，圖中曲線的頂點位置越靠近邊緣則顯示MSE誤差比越大，說明相應的混合模型相比于ARIMA模型的MSE誤差更低。可以看出，ARIMA_DLSTM模型在所有時間序列中MSE誤差最低，其中在Colorado River flow時間序列中的誤差比最大，優勢最為明顯。

圖15 5個數據集下的MSE誤差比雷達圖

8 結 語

本文提出了一種基于ARIMA模型、SVR模型和深度LSTM模型的ARIMA_DLSTM時間序列混合預測模型。ARIMA模型和SVR模型能夠分別提取時間序列的線性特征和誤差序列的非線性特征，深度LSTM模型對線性和非線性預測結果進行非線性組合。對來自不同領域的時間序列進行實證分析，實驗結果表明，提出的ARIMA_DLSTM模型和其他混合模型相比，能夠有效提高預測精度，有一定的實際應用價值。