唐欣欣,鄧光明,b
(桂林理工大學a.理學院;b.應用統計研究所,廣西桂林541006)
在時間序列分析中,一般會假設認為某一變量的未來值能夠完全由該變量的現在和過去值預測,而不受到其他變量的影響。由于這個假設很強,通常不能完全滿足,所以誤差是不可避免的。例如對于平穩的時間序列會構造ARMA模型對數據進行擬合,但是往往擬合的效果并不能達到預期,為了能提高預測的準確性,這個時候就要通過一些特殊的統計方法對ARMA模型進行進一步的改進。
目前在改進ARMA模型方面已經有很多研究,有學者將ARMA模型與其他的模型相結合,吳朝陽[1]基于灰色GM(1,1)模型和ARMA模型得到組合模型GM-ARMA模型減小了誤差,邵龍鋒[2]把小波變換的思想引入到ARMA模型中對其進行改進而后預測,任強等[3]用Leslie矩陣和ARMA模型得到人口隨機預測方法。同時也有學者在已構造的ARMA模型基礎上,通過調整各項參數達到提高預測精度的目的,何永沛[4]結合優化理論中的阻尼最小二乘法求解ARMA模型參數并在預測性能上有了較大的提高,黃雁勇等[5]通過DFP算法構造具有遺傳對稱正定性的矩陣來近似Hesse矩陣的逆從而實現ARMA模型的參數估計,孫汝儒等[6]提出用改進的PSO算法對ARMA(r,m)模型定階的新方法。
本文以我國國內航線的旅客周轉量為例,考慮到航空客運行業會因為節假、氣候等因素出現市場波動,需要一種精度較高的時間序列短期預測模型進行擬合,所以在剔除了季節趨勢之后用ARMA模型解決這類時間序列建模的問題。但是由于預測的誤差較大,為提高精度引入了Gevers-Wouters算法調整模型的參數,以此達到改進模型的目的。
觀測到的時間序列用{Yt}表示,并用{et}表示未觀測到的白噪聲序列(即一列均值為零的獨立同分布的隨機變量)。
一般線性過程{Yt}可以表示為現在與過去白噪聲變量的一種加權線性組合,如公式(1)所示:

如果公式(1)的右側是一個無窮級數,為了讓表達式具有數學意義,則需要對權數ψ做出如下假設:

因為{et}是無法觀測到的白噪聲序列,所以假設et的系數為1不會導致公式(2)失去一般性,即ψ0=1[7]。
當有限個系數ψ不為零的時候,得到移動平均過程,寫成公式(3)的形式:

公式(3)被稱之為q階移動平均過程,簡單記為MA(q)。
以自身做為回歸變量的過程稱之為自回歸過程,記為AR(p),即序列Yt的當期值是自身最近p階滯后項和新信息項et的線性組合,由此得到p階自回歸過程寫為公式(4)的形式:

其中,et包括了序列在t期無法用過去數值來解釋的所有新信息。對于每一個時期,假設et獨立于Yt-1,Yt-2,Yt-3,…。
如果假定序列中部分為自回歸過程,部分為滑動平均過程,由此就可以得到一個很普遍的時間序列模型,自回歸移動平均過程,用公式(5)表示如下:

可以簡單記為ARMA(p,q)。這種情況下,時刻t系統值Yt不僅和以前時刻的序列值相關,而且同以前時刻進入系統的擾動項也存在著一定的依存關系[8]。
本文從2010年1月至2016年3月的中國民航主要運輸生產指標統計獲取國內航線的旅客周轉量x1(單位:萬人·公里),繪制出時序圖,如圖1所示。

圖1 國內航線旅客周轉量時序圖
從圖1可以看出,我國國內航線的旅客周轉量自2010年以來呈現出穩步上升的趨勢,不可忽視的是,其中伴隨著一定的季節波動。導致季節變動的原因是多種多樣的,譬如氣候、節假日等都是使時間序列發生規律性變動的因素。通過季節調整,消除序列中的季節性影響,能夠更清晰地揭示趨勢。為了可以更加準確地反映客觀經濟的本質屬性,就需要對時間序列中的季節變動因素采取一定消除和調整的方法。為了能夠準確地構建時間序列模型,需要剔除掉數據中的季節趨勢部分。基于Loess方法對旅客周轉量數據做季節趨勢分解,得到三個部分,如圖2所示,自上而下分別是趨勢、季節以及冗余。

圖2 國內航線旅客周轉量的趨勢分解圖
輸出季節部分的數據,各月季節波動的數據如表1所示。從表1中可以看出,在每年的8月份左右達到一個明顯的高峰,在11月份左右產生低谷。了解了國內航線旅客周轉量的基本趨勢之后,需要進一步擬合曲線才能有效預測未來發展趨勢。

表1 各月季節變化情況
首先需要去除季節變動成分,將剩余部分的數據取對數后進行ADF檢驗。ADF檢驗顯著性水平為0.01,0.05和0.1的臨界值分別是-4.04,-3.45,-3.15(值越小越顯著)。該數據的檢驗統計量為-6.5804,在顯著性水平為0.01、0.05和0.1時均小于臨界值,說明這組對數序列是平穩的。所以可以直接對這組對數數據構建ARMA模型。
構建ARMA模型一般情況下是輸出acf和pacf這兩個條形圖用以判斷階數,但是這樣的判斷帶有一定的主觀性而且很多時候階數并不容易看出。所以,本文用R軟件程序包TSA中的函數armasubsets()根據BIC準則來判斷階數。
本文將最大的p和q均設置為14,對于不同的p和q會顯示出對應的BIC值。BIC是從貝葉斯的角度發展而來,近似于要求后驗模型的概率最大,而先驗模型在所有模型中有均勻的分布[9]。
BIC(貝葉斯信息量)的表達公式如下:

其中,L表示的是最大似然,n表示的是數據數量,k表示模型中的變量個數。
去除季節變動的國內航線旅客周轉量At擬合的時間序列模型如公式(7)所示:

確定時間序列模型的階數后,可通過Ljung-Box檢驗了解模型的擬合情況。
Ljung-Box檢驗的零假設是對于某個滯后序列獨立,其p值小說明可能有相關性,而對于不相關的觀測值(如:純隨機過程)p值則應該很大。對數據的殘差做Ljung-Box檢驗,并將檢驗得到的p值做出點圖,如圖3所示。

圖3 數據殘差的Ljung-Box檢驗的p值
從圖3可以看出,在滯后1~5的Ljung-Box檢驗的p值較小,這個滯后范圍看不出顯著的不相關性,但是對較長滯后(5~30)做檢驗p值呈現出不斷增大的趨勢并接近于1。對于所有的滯后階數p值均大于0.05這一水平,因此不能夠拒絕原假設,即數據的殘差不存在相關性[10]。所以,可以用公式(7)這個模型對國內航線旅客周轉量的情況進行預測。
通過構造的時間序列模型預測2016年4月至2017年3月lnAt的值,并繪制出時序圖,如圖4所示。

圖4 ARMA模型預測圖
圖4標出了誤差允許范圍內的旅客運輸周轉量預測區間。在此基礎上,加上各月季節變動量(見表1)求出旅客周轉量預測值,即得到未來十二個月我國民航旅客運輸周轉量,如表2所示。目前中國民用航空局最新的數據公布為2016年4月的中國民航主要生產指標,通過與真實值比較可以發現,誤差值為1.82%,該模型預測的誤差還是較大的。為提高預測的精度,還需要對模型進行進一步的改進。

表2 旅客周轉量預測值和真實值之間的比較
模型改進的目標主要是對國內航線構造的ARMA模型的參數進行合理的調整,由此引入ARMA新息模型的概念。ARMA新息模型是現代時間序列分析方法的基本工具,它是觀測信號的ARMA模型。基于ARMA新息模型,可以把最優濾波問題轉化為由新息序列生成的線性空間中的射影問題[11]。構造ARMA新息模型的關鍵就是構造其MA(移動平均過程)部分。
估計MA部分的參數用到了Gevers-Wouters算法[12]。
用r(t)表示m維的移動平均模型MA(n):
該試驗選擇在水稻成熟期(10月末)展開植株采樣和水稻經濟性狀測試,專門對水稻的株高、有效穗、每穗實粒數、結實率以及千粒重進行對比分析。實驗中采用了烘干法進行產量折算,最后對水稻生產的經濟效益進行了精確評估。


設初始時刻t0=-∞,則r(t)是一個平穩隨機序列。已知r(t)的階次nd和相關函數:

由Rr(k)求MA參數d1,…,dnd和[13]。
考慮可逆的純量MA(nd)過程即公式(8),已知其相關函數Rr(k),k=0,1,…,nd,則可用Gevers-Wouters迭代算法求MA參數di和:

其中t=0,1,2,…,i=t,t-1,…,0且規定:

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

通過迭代次數t=50~100便可以達到參數估計誤差不超過10-3的滿意精度。
將上述單變量MA參數估計的Gevers-Wouters算法推廣到多變量的情形,考慮可逆的多變量MA過程:

其中r(t)∈Rm,ε(t)∈Rm是零均值、方差陣為Qε的白噪聲,q-1為單位滯后算子,且多項式矩陣D(q-1)=Im+D1q-1+…+Dndq-nd是穩定的(即detD(x)的零點全位于單位圓外)。記它的相關函數為Rr(i)=E[r(t) rT(t-i)],i=1,…,nd;Rr(i)=0(i>nd)。
將構造的國內航線旅客周轉量ARMA(6,7)模型中移動平均過程單獨提取出來,通過Gevers-Wouters算法重新估計參數。畫出t=1~100時MA的各項參數估值,如圖5所示。

圖5 用Gevers-Wouters算法進行MA參數估計
圖5中分別顯示的是在迭代次數t=1~100時MA的各項參數的估值(即),可以看出前期參數的估值波動很頻繁,在t=80~100時參數的估值趨向于穩定。由于用于擬合模型的原數據(2010年1月至2016年3月)為75組,為了能提高預測的準確性,于是根據Gevers-Wouters算法得到的參數估值對之后每個時期分別構造模型。

根據結果調整ARMA(6,7)模型中MA過程的各個參數,由此計算出=15.461932,最后得到2016年4月國內航線旅客周轉量的預測值為5032135.06。
以此類推,分別可以得到2016年5月至2017年3月的國內航線旅客周轉量模型,并按照各個模型計算出預測值。通過整理匯總,與ARMA模型參數調整前的預測值相比較,如表3所示。從表3可以看出,基于Gevers-Wouters算法改進ARMA模型的參數后,誤差值由原來的1.82%降到了0.74%,預測的精度明顯提高。

表3 模型改進前后預測值和真實值之間的比較
本文引入Gevers-Wouters算法,對ARMA模型中MA部分的各項參數進行調整。本文以國內航線的旅客周轉量為例做實證分析,將訓練集用于擬合ARMA模型,通過貝葉斯準則定階,通過G-W算法對模型參數進行調整,最后得到的預測結果有了明顯的改善。試驗證明基于Gevers-Wouters算法改進ARMA模型可以有效提高預測的精度,且這種方法適用于其他的時間序列數據的預測研究,具有普遍性。
[1]吳朝陽.改進的灰色模型與ARMA模型的股指預測[J].智能系統學報,2010,5(3).
[2]邵龍鋒.小波變換下ARMA改進模型預測話務總量的研究[D].重慶:重慶大學碩士論文,2015.
[3]任強,侯大道.人口預測的隨機方法:基于Leslie矩陣和ARMA模型[J].人口研究,2011,35(2).
[4]何永沛.ARMA模型參數估計算法改進及在股票預測中的應用[J].重慶工學院學報,2009,23(2).
[5]黃雁勇,王沁,李裕奇.ARMA模型參數估計算法的改進[J].統計與決策,2009,(16).
[6]孫汝儒,肖迪.基于改進PSO算法對ARMA模型定階新方法[J].計算機應用與軟件,2013,30(12).
[7]Jonathan D.Cryer,Kung-Sik Chan.時間序列分析及應用:R語言[M].北京:機械工業出版社,2011.
[8]王婷.民航客運量的ARIMA模型與預測[J].五邑大學學報,2007,21(1).
[9]宋勇林.拓展的貝葉斯信息準則的一些性質[D].武漢:華中師范大學碩士論文,2014.
[10]Jafar Nouri Hajwal Hashem.時間序列模型的應用研究[D].武漢:華中師范大學碩士論文,2014.
[11]鄧自立,王欣,高媛.建模與估計[M].北京:科學出版社,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).