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

王英偉 馬樹才

(遼寧大學(xué)經(jīng)濟學(xué)院 遼寧 沈陽 110036)

0 引 言

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

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

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

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

1 ARIMA模型

ARIMA模型是由Box和Jenkins提出的用于時間序列分析和預(yù)測的線性模型。在時間序列分析中，需要設(shè)置3個參數(shù)，分別為自回歸階數(shù)(p)、差分階數(shù)(d)和移動平均階數(shù)(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)分別為自回歸參數(shù)和移動平均參數(shù)；εt為符合N(0,σ2)正態(tài)分布的誤差項。式(2)和式(3)分別為自回歸相關(guān)系數(shù)多項式和移動平均系數(shù)多項式。

時間序列yt應(yīng)先采用ADF檢驗選取差分階數(shù)，將時間序列變?yōu)槠椒€(wěn)序列，再根據(jù)赤池信息準則(Akaike Information Criterion，AIC)選擇最佳模型。

2 支持向量回歸

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

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

(4)

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

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

(5)

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

(6)

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

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

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

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

(7)

3 LSTM模型

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

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

圖1 LSTM記憶單元圖

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

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

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

(8)

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

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

(9)

(10)

3) 將舊單元狀態(tài)St-1更新為新單元狀態(tài)St：

(11)

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

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

(12)

ht=Ot?tanh(St)

(13)

LSTM神經(jīng)網(wǎng)絡(luò)要求輸入連續(xù)時間步上的特征向量值以便進行訓(xùn)練。假設(shè)每個序列共有N個時間步，并且SN表示第N個時間步的特征向量,所以第M個序列可以表示為{SM,SM+1,…,SM+N-1},輸入特征向量序列的構(gòu)造方式如圖2所示。

圖2 LSTM神經(jīng)網(wǎng)絡(luò)的輸入特征向量序列構(gòu)造方式

4 深度LSTM神經(jīng)網(wǎng)絡(luò)

深度LSTM神經(jīng)網(wǎng)絡(luò)是由多個基本LSTM模塊組成，核心思想是在較低LSTM層建立輸入數(shù)據(jù)的局部特征，并在較高LSTM層進行整合。實驗證明，深度LSTM神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的學(xué)習(xí)能力更強[12]。

圖3 深度LSTM網(wǎng)絡(luò)結(jié)構(gòu)

5 貝葉斯優(yōu)化算法

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

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

(14)

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

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

6 ARIMA_DLSTM模型

圖4 ARIMA_DLSTM模型流程

7 實 驗

7.1 實驗數(shù)據(jù)與預(yù)處理

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

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

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

(15)

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

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

7.2 衡量指標

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

(16)

(17)

(18)

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

7.3 實驗結(jié)果

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

表1 模型參數(shù)及取值范圍

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

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

圖5 貝葉斯優(yōu)化算法在Canadian Lynx的驗證誤差迭代圖

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

表2 對Canadian Lynx時間序列預(yù)測結(jié)果

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

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

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

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

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

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

表3 對Colorado River flow時間序列預(yù)測結(jié)果

圖8給出ARIMA模型和ARIMA_DLSTM模型在Colorado River flow測試集上的預(yù)測結(jié)果。可以看出，ARIMA_DLSTM模型預(yù)測值存在滯后現(xiàn)象，但和真實值更接近，優(yōu)于ARIMA模型。

圖8 ARIMA模型和ARIMA_DLSTM模型在Colorado River flow的時序預(yù)測結(jié)果

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

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

圖9 貝葉斯優(yōu)化算法在Airline Passengers的驗證誤差迭代圖

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

表4 對Airline Passengers時間序列預(yù)測結(jié)果

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

圖10 ARIMA模型和ARIMA_DLSTM模型在Airline Passengers 時間序列預(yù)測結(jié)果

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

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

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

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

表5 對Bitcoin匯率時間序列預(yù)測結(jié)果

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

圖12 ARIMA模型和ARIMA_DLSTM模型在Bitcoin匯率時間序列預(yù)測結(jié)果

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

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

圖13 貝葉斯優(yōu)化算法在糖價格指數(shù)的驗證誤差迭代圖

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

表6 對糖價格指數(shù)時間序列預(yù)測結(jié)果

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

圖14 ARIMA模型和ARIMA_DLSTM模型在該時間序列預(yù)測結(jié)果

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

圖15 5個數(shù)據(jù)集下的MSE誤差比雷達圖

8 結(jié) 語

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