王 安

(桐城師范高等?？茖W(xué)校 商貿(mào)與電子信息系,安徽 桐城 231400)

貝葉斯時間序列建模與分析是統(tǒng)計學(xué)中一個重要的理論與實踐課題。早期的學(xué)術(shù)研究開始了貝葉斯時間序列建模[1-2]。在歐美國家,貝葉斯模型首先取得了良好的效果,對于預(yù)測分析的效果改善明顯,而使得其得到了廣泛的使用。相較于非貝葉斯模型,貝葉斯模型更加靈活,正因為這個特點,更多的學(xué)者、專家以及相關(guān)行業(yè)的從業(yè)者,逐漸對貝葉斯模型持認可的態(tài)度[3-4]。尤其是近些年來,貝葉斯方法在社會經(jīng)濟預(yù)測分析中的應(yīng)用逐年增加,并取得了很大的成功。貝葉斯統(tǒng)計推理可以應(yīng)用的范圍十分廣泛,上至核電廠的工程質(zhì)量控制、國家經(jīng)濟民生相關(guān)的數(shù)據(jù)統(tǒng)計分析,下至對社會犯罪情況的分析與預(yù)測,包羅萬象。貝葉斯方法完全可以說是人們生產(chǎn)生活中的一個有力的工具[5]。但就目前的研究狀況來說，貝葉斯方法在我國的應(yīng)用和發(fā)展還處于起步階段:一是與西方國家相比,我國目前的觀點受到頻率學(xué)派的限制,大多數(shù)統(tǒng)計專著和學(xué)術(shù)論文缺乏嚴格、系統(tǒng)的討論，使國內(nèi)學(xué)者對貝葉斯方法感到陌生[6-7]。二是由于缺乏相應(yīng)的軟件支持,高維數(shù)據(jù)分析中的復(fù)雜問題難以計算。該問題的求解阻礙了貝葉斯方法在實際問題中的應(yīng)用。

1 ARMA模型

1.1 模型

作為最常用的擬合平穩(wěn)時間序列的模型,ARMA模型共有三類,分別為AR模型、MA模型和ARMA模型。ARMA模型本身是一種自回歸的移動平均模型,具有自回歸模型和移動平均模型二者的特性。

AR(p)模型——p階自回歸模型

xt=φ0+φ1xt-1+…φpxt-p+εt

由于該序列是平穩(wěn)序列,據(jù)此可推導(dǎo)得出均值為如下形式:

若φ0=0,則該模型稱為中心化的AR(p)模型;若想對非中心化的平穩(wěn)時間序列,進行中心化操作,則可做如下的變換:

Φp(B)=I-φ1B-…-φpBp稱為p階自回歸多項式,則AR(p)模型可表示為:Φp(B)xt=εt,其中為延遲算子。

1.2 格林函數(shù)

格林函數(shù)反映了影響效應(yīng)衰減的快慢程度(回到平衡位置的速度)是用來描述系統(tǒng)記憶擾動程度的函數(shù),Gj表示擾動εt-j對系統(tǒng)現(xiàn)在行為影響的權(quán)數(shù)。

1.3 模型的方差

1.4 模型的自協(xié)方差

對中心化的平穩(wěn)模型,可推得自協(xié)方差函數(shù)的遞推公式:

γ(K)=φ1γ(k-1)+φ2γ(k-2)+…+φpγ(k-p)

用格林函數(shù)顯示表示:

對于AR(1)模型,

1.5 模型的自相關(guān)函數(shù)

遞推公式:

ρ(k)=φ1ρ(k-1)+φ2ρ(k-2)+…+φpρ(k-p)

如圖1所示,平穩(wěn)AR(p)模型的自相關(guān)函數(shù)有兩個顯著的性質(zhì)[8]:

1)拖尾性。自相關(guān)函數(shù)總是有一個非零值,當k大于某個常數(shù)時,ρ(k)不等于零。

2)負指數(shù)衰減。隨著時間的推移,自相關(guān)函數(shù)ρ(k)會迅速衰減,且以負指數(shù)(其中λi為自相關(guān)函數(shù)差分方程的特征根)的速度在減小。

(a)自相關(guān)函數(shù)ACF(k)xt=0.8xt-1+εtFig(a) Autocorrelation function of ACF(k)xt=0.8xt-1+εt

(b)自相關(guān)函數(shù)ACF(k)xt=-0.8xt-1+εtFig(b) Autocorrelation function of ACF(k)xt=-0.8xt-1+εt

圖1AR(1)模型的自相關(guān)函數(shù)
Fig.1AutocorrelationfunctionofAR(1)model

1.6 模型的偏自相關(guān)函數(shù)

如圖2(a-b)所示,自相關(guān)函數(shù)ρ(k)實際上并不只是xt與xt-k之間的相關(guān)關(guān)系,它還會受到中間k-1個隨機變量xt-1, …,xt-k+1的影響[9]。為了消除隨機變量在中間造成的干擾,我們將單純測量xt與xt-k之間的相關(guān)關(guān)系,引入了滯后k偏自相關(guān)函數(shù)(PACF),計算公式為:

其中,

實際上k階自回歸模型第k個回歸系數(shù)φkk就是滯后k偏自相關(guān)函數(shù):

xt=Φk1xt-1+Φk2xt-2+…+Φkkxt-k+εt

兩邊同乘以xt-k,求期望再除以γ(0)得到:

ρt=Φk1ρt-1+Φk2ρt-2+…+Φkkρt-k,t=1,2,…,k,…,n

取前K個方程構(gòu)成的方程組:

…

ρ1=Φk1ρk-1+Φk2ρk-2+…+Φkkρ0

稱為Yule-Walker方程,可以解出Φkk。

(a)偏自相關(guān)函數(shù)PACF(k)xt=0.8xt-1+εt(a) Partial autocorrelation function of PACF(k)xt=0.8xt-1+εt

(b)偏自相關(guān)函數(shù)PACF(k)xt=-0.8xt-1+εt(b) Partial autocorrelation function of PACF(k)xt=-0.8xt-1+εt

圖2AR(1)模型的偏自相關(guān)函數(shù)
Fig.2PartialautocorrelationfunctionofAR(1)model

2 ARIMA模型——混合自回歸移動平均模型

ARIMA 是一種單變量同方差的線性模型,對于滿足參數(shù)模型的平穩(wěn)時間序列,主要有以下三種基本形式:自回歸模型(AR:Auto-Regression),移動平均模型 (MA:Moving-Average)和混合模型 (ARIMA:Auto-Regression-Moving-Average)[9]。

簡短地說, AR模型是建立當前值和歷史值之間的聯(lián)系,MA模型是計算AR部分的誤差的累計,此模型是分析兩者與差分的和。

2.1 原理

它用于處理單變量同方差的非平穩(wěn)時間序列,并通過差分法或適當?shù)淖儞Q轉(zhuǎn)化為平穩(wěn)時間序列,再使用ARMA模型,這種模型也稱Box-Jenkins模型。

殘差的條件方差是異方差的時間序列,適合用GARCH模型。ARIMA(p,d,q)模型的形式如下:

其中,Δd=(I-B)d為d階差分,

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

Θ(B)=1-θ1B-θ2B2-…-θpBp

d階差分后的序列可表示為:

2.2 建模流程

2.2.1識別階段

在識別后面ARIMA模型之前,首先需要指定響應(yīng)變量序列,計算自相關(guān)系數(shù)(ACF)、逆自相關(guān)系數(shù)(IACF)、偏自相關(guān)系數(shù)(PACF)和相關(guān)數(shù)。在這個工作階段,通常有幾個ARIMA模型可供選擇,要區(qū)分模型的類型需要通過檢驗樣本自相關(guān)系數(shù)(SACF)和樣本偏自相關(guān)系數(shù)(SPACF)。

2.2.2估計階段

指定的響應(yīng)變量需要通過指定的ARIMA模型來擬合,同時進行該模型的參數(shù)估計。模型的適用性需要通過診斷統(tǒng)計量來判斷。

考量參數(shù)估計值的好壞有幾個重要的指標,從多個維度檢驗估計結(jié)果的合理性。

1)模型里各項是否必要通過顯著性校驗來進行;

2)各模型的優(yōu)劣通過R2來決定;

3)殘差序列所蘊含額外信息通過白噪聲殘差檢驗來指明。

經(jīng)過以上各項檢查后,若發(fā)現(xiàn)所選模型并不合適,可嘗試其他模型,重新開始估計診斷。

3 預(yù)測階段

預(yù)測時間序列的未來值,并根據(jù)ARIMA模型的預(yù)測值生成置信區(qū)間如圖3所示。

圖3 xlog殘差相關(guān)診斷Fig.3 xlog residual-related diagnosis

圖4 xlog殘差正態(tài)診斷Fig.4 xlog residual normal diagnosis

按AIC和SBC標準計算的統(tǒng)計量,值較小,表明對模型擬合得較好。

另外,從-441.319 5<-438.081 43,-435.569 11<-426.580 65可看出,本次擬合的ARIMA(0,1,1)×(0,1,1)12模型比上次擬合的模型要好。

估計值之間的相關(guān)系數(shù)為-0.064,這是一個較小的相關(guān)系數(shù),如果這個相關(guān)系數(shù)較大時,就需要考慮是否刪除其中一個參數(shù)。

對應(yīng)的QLB統(tǒng)計量分別為6.54和5.60。其所對應(yīng)的p值分別為0.768和0.231,其所蘊含的信息均已被提取。

最后,如圖5所示將表中移動平均MA的兩個因子(1-0.50275B)和(1-0.52993B12)代入ARIMA(0,1,1)×(0,1,1)12模型得到:

(1-B)(1-B12)Zt=(1-0.50275B)(1-0.52993B12)εt

其中,Z=log(xt)

圖5 xlog預(yù)測Fig.5 xlog prediction

圖6 xlog置信度曲線Fig.6 xlog confidence curve

預(yù)測今后一年(12期)國際航線各個月度的旅客人數(shù),結(jié)果存在數(shù)據(jù)集forxlog中。

最后,繪制置信度曲線,在圖中標注預(yù)測數(shù)據(jù)線和原始數(shù)據(jù)線的時序圖,置信度為95%如圖6所示。

4 結(jié)語

先驗概率的概念一直是貝葉斯統(tǒng)計和非貝葉斯統(tǒng)計的研究熱點之一。然而,隨著科學(xué)的進步,貝葉斯統(tǒng)計在實際應(yīng)用中的成功已經(jīng)慢慢改變了人們的看法。貝葉斯統(tǒng)計已經(jīng)受到人們的重視,成為統(tǒng)計研究的熱點。貝葉斯統(tǒng)計推斷方法體系是近幾十年來迅速發(fā)展起來的一種系統(tǒng)建模方法。本文對時間序列預(yù)測模型的常用模型進行了深入闡述,提供了標準的建模、估計以及優(yōu)化思路和流程。在今后的研究工作中,在考慮貝葉斯方法高精度的優(yōu)勢的同時,需要解決貝葉斯估計計算復(fù)雜、開銷較大的問題。