基于Gevers-Wouters算法的ARMA模型改進(jìn)

唐欣欣，鄧光明,b

（桂林理工大學(xué)a.理學(xué)院；b.應(yīng)用統(tǒng)計(jì)研究所，廣西桂林541006）

0 引言

在時(shí)間序列分析中，一般會(huì)假設(shè)認(rèn)為某一變量的未來值能夠完全由該變量的現(xiàn)在和過去值預(yù)測(cè)，而不受到其他變量的影響。由于這個(gè)假設(shè)很強(qiáng)，通常不能完全滿足，所以誤差是不可避免的。例如對(duì)于平穩(wěn)的時(shí)間序列會(huì)構(gòu)造ARMA模型對(duì)數(shù)據(jù)進(jìn)行擬合，但是往往擬合的效果并不能達(dá)到預(yù)期，為了能提高預(yù)測(cè)的準(zhǔn)確性，這個(gè)時(shí)候就要通過一些特殊的統(tǒng)計(jì)方法對(duì)ARMA模型進(jìn)行進(jìn)一步的改進(jìn)。

目前在改進(jìn)ARMA模型方面已經(jīng)有很多研究，有學(xué)者將ARMA模型與其他的模型相結(jié)合，吳朝陽[1]基于灰色GM(1,1)模型和ARMA模型得到組合模型GM-ARMA模型減小了誤差，邵龍鋒[2]把小波變換的思想引入到ARMA模型中對(duì)其進(jìn)行改進(jìn)而后預(yù)測(cè)，任強(qiáng)等[3]用Leslie矩陣和ARMA模型得到人口隨機(jī)預(yù)測(cè)方法。同時(shí)也有學(xué)者在已構(gòu)造的ARMA模型基礎(chǔ)上，通過調(diào)整各項(xiàng)參數(shù)達(dá)到提高預(yù)測(cè)精度的目的，何永沛[4]結(jié)合優(yōu)化理論中的阻尼最小二乘法求解ARMA模型參數(shù)并在預(yù)測(cè)性能上有了較大的提高，黃雁勇等[5]通過DFP算法構(gòu)造具有遺傳對(duì)稱正定性的矩陣來近似Hesse矩陣的逆從而實(shí)現(xiàn)ARMA模型的參數(shù)估計(jì)，孫汝儒等[6]提出用改進(jìn)的PSO算法對(duì)ARMA(r,m)模型定階的新方法。

本文以我國(guó)國(guó)內(nèi)航線的旅客周轉(zhuǎn)量為例，考慮到航空客運(yùn)行業(yè)會(huì)因?yàn)楣?jié)假、氣候等因素出現(xiàn)市場(chǎng)波動(dòng)，需要一種精度較高的時(shí)間序列短期預(yù)測(cè)模型進(jìn)行擬合，所以在剔除了季節(jié)趨勢(shì)之后用ARMA模型解決這類時(shí)間序列建模的問題。但是由于預(yù)測(cè)的誤差較大，為提高精度引入了Gevers-Wouters算法調(diào)整模型的參數(shù)，以此達(dá)到改進(jìn)模型的目的。

1 理論基礎(chǔ)

觀測(cè)到的時(shí)間序列用{Yt}表示，并用{et}表示未觀測(cè)到的白噪聲序列（即一列均值為零的獨(dú)立同分布的隨機(jī)變量）。

一般線性過程{Yt}可以表示為現(xiàn)在與過去白噪聲變量的一種加權(quán)線性組合，如公式(1)所示：

如果公式(1)的右側(cè)是一個(gè)無窮級(jí)數(shù)，為了讓表達(dá)式具有數(shù)學(xué)意義，則需要對(duì)權(quán)數(shù)ψ做出如下假設(shè)：

因?yàn)閧et}是無法觀測(cè)到的白噪聲序列，所以假設(shè)et的系數(shù)為1不會(huì)導(dǎo)致公式(2)失去一般性，即ψ0=1[7]。

當(dāng)有限個(gè)系數(shù)ψ不為零的時(shí)候，得到移動(dòng)平均過程，寫成公式(3)的形式：

公式(3)被稱之為q階移動(dòng)平均過程，簡(jiǎn)單記為MA(q)。

以自身做為回歸變量的過程稱之為自回歸過程，記為AR(p)，即序列Yt的當(dāng)期值是自身最近p階滯后項(xiàng)和新信息項(xiàng)et的線性組合，由此得到p階自回歸過程寫為公式(4)的形式：

其中，et包括了序列在t期無法用過去數(shù)值來解釋的所有新信息。對(duì)于每一個(gè)時(shí)期，假設(shè)et獨(dú)立于Yt-1,Yt-2,Yt-3,…。

如果假定序列中部分為自回歸過程，部分為滑動(dòng)平均過程，由此就可以得到一個(gè)很普遍的時(shí)間序列模型，自回歸移動(dòng)平均過程，用公式(5)表示如下：

可以簡(jiǎn)單記為ARMA(p,q)。這種情況下，時(shí)刻t系統(tǒng)值Yt不僅和以前時(shí)刻的序列值相關(guān)，而且同以前時(shí)刻進(jìn)入系統(tǒng)的擾動(dòng)項(xiàng)也存在著一定的依存關(guān)系[8]。

2 數(shù)據(jù)來源及預(yù)處理

本文從2010年1月至2016年3月的中國(guó)民航主要運(yùn)輸生產(chǎn)指標(biāo)統(tǒng)計(jì)獲取國(guó)內(nèi)航線的旅客周轉(zhuǎn)量x1（單位：萬人·公里），繪制出時(shí)序圖，如圖1所示。

圖1 國(guó)內(nèi)航線旅客周轉(zhuǎn)量時(shí)序圖

從圖1可以看出，我國(guó)國(guó)內(nèi)航線的旅客周轉(zhuǎn)量自2010年以來呈現(xiàn)出穩(wěn)步上升的趨勢(shì)，不可忽視的是，其中伴隨著一定的季節(jié)波動(dòng)。導(dǎo)致季節(jié)變動(dòng)的原因是多種多樣的，譬如氣候、節(jié)假日等都是使時(shí)間序列發(fā)生規(guī)律性變動(dòng)的因素。通過季節(jié)調(diào)整，消除序列中的季節(jié)性影響，能夠更清晰地揭示趨勢(shì)。為了可以更加準(zhǔn)確地反映客觀經(jīng)濟(jì)的本質(zhì)屬性，就需要對(duì)時(shí)間序列中的季節(jié)變動(dòng)因素采取一定消除和調(diào)整的方法。為了能夠準(zhǔn)確地構(gòu)建時(shí)間序列模型，需要剔除掉數(shù)據(jù)中的季節(jié)趨勢(shì)部分?；贚oess方法對(duì)旅客周轉(zhuǎn)量數(shù)據(jù)做季節(jié)趨勢(shì)分解，得到三個(gè)部分，如圖2所示，自上而下分別是趨勢(shì)、季節(jié)以及冗余。

圖2 國(guó)內(nèi)航線旅客周轉(zhuǎn)量的趨勢(shì)分解圖

輸出季節(jié)部分的數(shù)據(jù)，各月季節(jié)波動(dòng)的數(shù)據(jù)如表1所示。從表1中可以看出，在每年的8月份左右達(dá)到一個(gè)明顯的高峰，在11月份左右產(chǎn)生低谷。了解了國(guó)內(nèi)航線旅客周轉(zhuǎn)量的基本趨勢(shì)之后，需要進(jìn)一步擬合曲線才能有效預(yù)測(cè)未來發(fā)展趨勢(shì)。

表1 各月季節(jié)變化情況

首先需要去除季節(jié)變動(dòng)成分，將剩余部分的數(shù)據(jù)取對(duì)數(shù)后進(jìn)行ADF檢驗(yàn)。ADF檢驗(yàn)顯著性水平為0.01，0.05和0.1的臨界值分別是-4.04，-3.45，-3.15（值越小越顯著）。該數(shù)據(jù)的檢驗(yàn)統(tǒng)計(jì)量為-6.5804，在顯著性水平為0.01、0.05和0.1時(shí)均小于臨界值，說明這組對(duì)數(shù)序列是平穩(wěn)的。所以可以直接對(duì)這組對(duì)數(shù)數(shù)據(jù)構(gòu)建ARMA模型。

3 構(gòu)建ARMA模型

構(gòu)建ARMA模型一般情況下是輸出acf和pacf這兩個(gè)條形圖用以判斷階數(shù)，但是這樣的判斷帶有一定的主觀性而且很多時(shí)候階數(shù)并不容易看出。所以，本文用R軟件程序包TSA中的函數(shù)armasubsets()根據(jù)BIC準(zhǔn)則來判斷階數(shù)。

本文將最大的p和q均設(shè)置為14，對(duì)于不同的p和q會(huì)顯示出對(duì)應(yīng)的BIC值。BIC是從貝葉斯的角度發(fā)展而來，近似于要求后驗(yàn)?zāi)Ｐ偷母怕首畲?，而先?yàn)?zāi)Ｐ驮谒心Ｐ椭杏芯鶆虻姆植糩9]。

BIC（貝葉斯信息量）的表達(dá)公式如下：

其中，L表示的是最大似然，n表示的是數(shù)據(jù)數(shù)量，k表示模型中的變量個(gè)數(shù)。

去除季節(jié)變動(dòng)的國(guó)內(nèi)航線旅客周轉(zhuǎn)量At擬合的時(shí)間序列模型如公式(7)所示：

確定時(shí)間序列模型的階數(shù)后，可通過Ljung-Box檢驗(yàn)了解模型的擬合情況。

Ljung-Box檢驗(yàn)的零假設(shè)是對(duì)于某個(gè)滯后序列獨(dú)立，其p值小說明可能有相關(guān)性，而對(duì)于不相關(guān)的觀測(cè)值（如：純隨機(jī)過程）p值則應(yīng)該很大。對(duì)數(shù)據(jù)的殘差做Ljung-Box檢驗(yàn)，并將檢驗(yàn)得到的p值做出點(diǎn)圖，如圖3所示。

圖3 數(shù)據(jù)殘差的Ljung-Box檢驗(yàn)的p值

從圖3可以看出，在滯后1～5的Ljung-Box檢驗(yàn)的p值較小，這個(gè)滯后范圍看不出顯著的不相關(guān)性，但是對(duì)較長(zhǎng)滯后（5～30）做檢驗(yàn)p值呈現(xiàn)出不斷增大的趨勢(shì)并接近于1。對(duì)于所有的滯后階數(shù)p值均大于0.05這一水平，因此不能夠拒絕原假設(shè)，即數(shù)據(jù)的殘差不存在相關(guān)性[10]。所以，可以用公式(7)這個(gè)模型對(duì)國(guó)內(nèi)航線旅客周轉(zhuǎn)量的情況進(jìn)行預(yù)測(cè)。

通過構(gòu)造的時(shí)間序列模型預(yù)測(cè)2016年4月至2017年3月lnAt的值，并繪制出時(shí)序圖，如圖4所示。

圖4 ARMA模型預(yù)測(cè)圖

圖4標(biāo)出了誤差允許范圍內(nèi)的旅客運(yùn)輸周轉(zhuǎn)量預(yù)測(cè)區(qū)間。在此基礎(chǔ)上，加上各月季節(jié)變動(dòng)量（見表1）求出旅客周轉(zhuǎn)量預(yù)測(cè)值，即得到未來十二個(gè)月我國(guó)民航旅客運(yùn)輸周轉(zhuǎn)量，如表2所示。目前中國(guó)民用航空局最新的數(shù)據(jù)公布為2016年4月的中國(guó)民航主要生產(chǎn)指標(biāo)，通過與真實(shí)值比較可以發(fā)現(xiàn)，誤差值為1.82%，該模型預(yù)測(cè)的誤差還是較大的。為提高預(yù)測(cè)的精度，還需要對(duì)模型進(jìn)行進(jìn)一步的改進(jìn)。

表2 旅客周轉(zhuǎn)量預(yù)測(cè)值和真實(shí)值之間的比較

4 模型的改進(jìn)

模型改進(jìn)的目標(biāo)主要是對(duì)國(guó)內(nèi)航線構(gòu)造的ARMA模型的參數(shù)進(jìn)行合理的調(diào)整，由此引入ARMA新息模型的概念。ARMA新息模型是現(xiàn)代時(shí)間序列分析方法的基本工具，它是觀測(cè)信號(hào)的ARMA模型?；贏RMA新息模型，可以把最優(yōu)濾波問題轉(zhuǎn)化為由新息序列生成的線性空間中的射影問題[11]。構(gòu)造ARMA新息模型的關(guān)鍵就是構(gòu)造其MA（移動(dòng)平均過程）部分。

估計(jì)MA部分的參數(shù)用到了Gevers-Wouters算法[12]。

用r(t)表示m維的移動(dòng)平均模型MA(n)：

該試驗(yàn)選擇在水稻成熟期（10月末）展開植株采樣和水稻經(jīng)濟(jì)性狀測(cè)試，專門對(duì)水稻的株高、有效穗、每穗實(shí)粒數(shù)、結(jié)實(shí)率以及千粒重進(jìn)行對(duì)比分析。實(shí)驗(yàn)中采用了烘干法進(jìn)行產(chǎn)量折算，最后對(duì)水稻生產(chǎn)的經(jīng)濟(jì)效益進(jìn)行了精確評(píng)估。

設(shè)初始時(shí)刻t0=-∞，則r(t)是一個(gè)平穩(wěn)隨機(jī)序列。已知r(t)的階次nd和相關(guān)函數(shù)：

由Rr(k)求MA參數(shù)d1,…,dnd和[13]。

考慮可逆的純量MA(nd)過程即公式(8)，已知其相關(guān)函數(shù)Rr(k)，k=0,1,…,nd，則可用Gevers-Wouters迭代算法求MA參數(shù)di和：

其中t=0,1,2,…，i=t,t-1,…,0且規(guī)定：

在較弱的情況下，上述的極限關(guān)系是成立的，并且可以保證MA(nd)過程是可逆的。應(yīng)用上述Gevers-Wouters算法，取迭代次數(shù)t充分大，則在時(shí)刻t時(shí)有估值：

通過迭代次數(shù)t=50～100便可以達(dá)到參數(shù)估計(jì)誤差不超過10-3的滿意精度。

將上述單變量MA參數(shù)估計(jì)的Gevers-Wouters算法推廣到多變量的情形，考慮可逆的多變量MA過程：

其中r(t)∈Rm，ε(t)∈Rm是零均值、方差陣為Qε的白噪聲，q-1為單位滯后算子，且多項(xiàng)式矩陣D(q-1)=Im+D1q-1+…+Dndq-nd是穩(wěn)定的（即detD(x)的零點(diǎn)全位于單位圓外）。記它的相關(guān)函數(shù)為Rr(i)=E[r(t) rT(t-i)]，i=1,…,nd；Rr(i)=0(i＞nd)。

將構(gòu)造的國(guó)內(nèi)航線旅客周轉(zhuǎn)量ARMA(6,7)模型中移動(dòng)平均過程單獨(dú)提取出來，通過Gevers-Wouters算法重新估計(jì)參數(shù)。畫出t=1～100時(shí)MA的各項(xiàng)參數(shù)估值，如圖5所示。

圖5 用Gevers-Wouters算法進(jìn)行MA參數(shù)估計(jì)

圖5中分別顯示的是在迭代次數(shù)t=1～100時(shí)MA的各項(xiàng)參數(shù)的估值（即），可以看出前期參數(shù)的估值波動(dòng)很頻繁，在t=80～100時(shí)參數(shù)的估值趨向于穩(wěn)定。由于用于擬合模型的原數(shù)據(jù)（2010年1月至2016年3月）為75組，為了能提高預(yù)測(cè)的準(zhǔn)確性，于是根據(jù)Gevers-Wouters算法得到的參數(shù)估值對(duì)之后每個(gè)時(shí)期分別構(gòu)造模型。

根據(jù)結(jié)果調(diào)整ARMA(6,7)模型中MA過程的各個(gè)參數(shù)，由此計(jì)算出=15.461932，最后得到2016年4月國(guó)內(nèi)航線旅客周轉(zhuǎn)量的預(yù)測(cè)值為5032135.06。

以此類推，分別可以得到2016年5月至2017年3月的國(guó)內(nèi)航線旅客周轉(zhuǎn)量模型，并按照各個(gè)模型計(jì)算出預(yù)測(cè)值。通過整理匯總，與ARMA模型參數(shù)調(diào)整前的預(yù)測(cè)值相比較，如表3所示。從表3可以看出，基于Gevers-Wouters算法改進(jìn)ARMA模型的參數(shù)后，誤差值由原來的1.82%降到了0.74%，預(yù)測(cè)的精度明顯提高。

表3 模型改進(jìn)前后預(yù)測(cè)值和真實(shí)值之間的比較

5 結(jié)束語

本文引入Gevers-Wouters算法，對(duì)ARMA模型中MA部分的各項(xiàng)參數(shù)進(jìn)行調(diào)整。本文以國(guó)內(nèi)航線的旅客周轉(zhuǎn)量為例做實(shí)證分析，將訓(xùn)練集用于擬合ARMA模型，通過貝葉斯準(zhǔn)則定階，通過G-W算法對(duì)模型參數(shù)進(jìn)行調(diào)整，最后得到的預(yù)測(cè)結(jié)果有了明顯的改善。試驗(yàn)證明基于Gevers-Wouters算法改進(jìn)ARMA模型可以有效提高預(yù)測(cè)的精度，且這種方法適用于其他的時(shí)間序列數(shù)據(jù)的預(yù)測(cè)研究，具有普遍性。

[1]吳朝陽.改進(jìn)的灰色模型與ARMA模型的股指預(yù)測(cè)[J].智能系統(tǒng)學(xué)報(bào),2010,5(3).

[2]邵龍鋒.小波變換下ARMA改進(jìn)模型預(yù)測(cè)話務(wù)總量的研究[D].重慶:重慶大學(xué)碩士論文,2015.

[3]任強(qiáng),侯大道.人口預(yù)測(cè)的隨機(jī)方法：基于Leslie矩陣和ARMA模型[J].人口研究,2011,35(2).

[4]何永沛.ARMA模型參數(shù)估計(jì)算法改進(jìn)及在股票預(yù)測(cè)中的應(yīng)用[J].重慶工學(xué)院學(xué)報(bào),2009,23(2).

[5]黃雁勇,王沁,李裕奇.ARMA模型參數(shù)估計(jì)算法的改進(jìn)[J].統(tǒng)計(jì)與決策,2009,(16).

[6]孫汝儒,肖迪.基于改進(jìn)PSO算法對(duì)ARMA模型定階新方法[J].計(jì)算機(jī)應(yīng)用與軟件,2013,30(12).

[7]Jonathan D.Cryer,Kung-Sik Chan.時(shí)間序列分析及應(yīng)用：R語言[M].北京:機(jī)械工業(yè)出版社,2011.

[8]王婷.民航客運(yùn)量的ARIMA模型與預(yù)測(cè)[J].五邑大學(xué)學(xué)報(bào),2007,21(1).

[9]宋勇林.拓展的貝葉斯信息準(zhǔn)則的一些性質(zhì)[D].武漢:華中師范大學(xué)碩士論文，2014.

[10]Jafar Nouri Hajwal Hashem.時(shí)間序列模型的應(yīng)用研究[D].武漢:華中師范大學(xué)碩士論文，2014.

[11]鄧自立,王欣,高媛.建模與估計(jì)[M].北京:科學(xué)出版社,2007.

[12]De Keyser R M C,Van Cauwenberghe A R.A Self-Tuning Multistep Preditor Application[J].Automatic,1981,17(1).

[13]Zhang G S,Xu Y.Application of Gevers-Wouters Algorithm to Intelligent Dew Point Measurer Control Scheme[C].Switzerland:Trans Tech Publications,2014,(643).