NGINAR(1)模型的Bayes估計及預測

楊艷秋,王德輝

(1.吉林大學 數(shù)學學院,長春 130012;2.吉林師范大學 數(shù)學學院,吉林 四平 136000)

現(xiàn)實生活中的許多計數(shù)過程,如某醫(yī)院某時刻住院的病患人數(shù)、某時刻某地區(qū)的生物種群數(shù)量、某地區(qū)某階段的犯罪數(shù)量等,這些數(shù)據(jù)通常取非負整數(shù)值,用一般的時間序列模型擬合這些數(shù)據(jù)通常會產(chǎn)生異常預測,因此需要引入整值自回歸模型.目前對整值時間序列數(shù)據(jù)的建模主要有如下兩種方法: 借助于潛過程的狀態(tài)空間建模過程[1]；借助于稀疏算子的建模過程,這類方法是建模整值時間序列數(shù)據(jù)的主要方法.二項稀疏算子“°”[2]是整值時間序列發(fā)展的基礎,定義為

其中:α∈(0,1);X為非負整值隨機變量;{Bi}為獨立同分布(i.i.d.)的Bernoulli隨機變量序列且與X互相獨立,滿足

P(Bi=1)=1-P(Bi=0)=α.

但二項稀疏算子有一個局限,即求和序列{Bi}為i.i.d.的Bernoulli隨機變量序列,只能取0或1值,因此α°X的取值總是小于或等于X的值.但在實際問題中,每一個事件都可能關(guān)聯(lián)更多相關(guān)的計數(shù)事件,因此用幾何隨機變量來描述這些事件更合適.文獻[3]引入了負二項稀疏算子“*”,定義為

k∈0.由于負二項稀疏算子中求和序列{Wj}的取值是非負整數(shù),使得β*X的取值可能大于也可能小于X的取值,從而很好地突破了二項稀疏算子的局限.本文給出基于負二項稀疏算子的一階整值自回歸模型參數(shù)的Bayes估計,先進行數(shù)值模擬,再與條件最小二乘估計和Yule-Walker估計進行均方誤差的比較,最后對新西蘭牛皮膚病數(shù)據(jù)進行實例分析及模型預測.

1 模型簡介

基于二項稀疏算子“°”的一階整值自回歸模型[4-5]為

Xt=α°Xt-1+Zt.

若Xt表示t時刻住院的患者數(shù),α°Xt-1表示上個月以概率α仍繼續(xù)住院的患者數(shù),Zt為t時刻新住院的患者數(shù),則Xt為α°Xt-1與Zt之和.但對于過度分散的計數(shù)過程,負二項稀疏算子“*”則更合適.文獻[6]基于負二項稀疏算子“*”,提出了一階整數(shù)值自回歸模型(new geometric first-order integer-valued autoregressive process,NGINAR(1))如下:

Xt=α*Xt-1+εt,

其中: 負二項稀疏算子“*”定義為

2 Bayes估計

其中

i=1,2,…,T.

考慮α的先驗分布取(0,1)上的均勻分布,即π(α)=1,0<α<1.根據(jù)Bayes原理,參數(shù)α的后驗分布為

綜上可得：

定理1設樣本x0,x1,…,xT來自于NGINAR(1)模型,在二次損失函數(shù)和均勻先驗分布下,參數(shù)α的Bayes估計形如式(1).

3 數(shù)值模擬

下面通過數(shù)值模擬給出NGINAR(1)模型參數(shù)Bayes估計的優(yōu)良性,將NGINAR(1)模型參數(shù)的Bayes估計與條件最小二乘(CLS)估計和Yule-Walker估計(Y-W)進行均方誤差的比較.先給出Bayes估計的算法:

1) 選擇迭代初值α(0),并令i=1;

5) 令i=i+1,返回步驟2),直到算法達到事先約定的收斂標準.

下面進行數(shù)值模擬,樣本容量分別取T=100,500.取μ=0.5,5,10,表1列出了Bayes估計值的偏差(Bias)和均方誤差(MSE),表2列出了參數(shù)的Bayes估計與條件最小二乘估計和Yule-Walker估計的均方誤差比.模擬運行過程中先進行2 000次預迭代,以確保參數(shù)的收斂性,然后再進行500次迭代,得到模擬結(jié)果.

表1 Bayes估計的偏差和均方誤差

表2 3種估計方法的均方誤差比

由表1可見,Bayes估計的偏差和均方誤差都比較小.以均方誤差為準則,表2中3種估計方法的均方誤差比可見,Bayes估計優(yōu)于條件最小二乘估計和Yule-Walker估計.

4 實例分析及模型預測

圖1 牛皮膚病數(shù)據(jù)樣本路徑Fig.1 Sample path of skin-lesions data

下面將NGINAR(1)模型應用到一組牛皮膚病患病數(shù)據(jù)集中[9],并進行分析及預測.該數(shù)據(jù)集來源于新西蘭農(nóng)林部,記錄了新西蘭某地區(qū)動物衛(wèi)生實驗室記錄的2003—2009年間每月患皮膚病的牛數(shù)量.將該數(shù)據(jù)集分為兩部分：2003-01—2009-08的數(shù)據(jù)用于估計參數(shù)值,2009-09—2009-12的數(shù)據(jù)作為樣本外待預測值.

將數(shù)據(jù)集的前80個觀測數(shù)據(jù)記作X1,X2,…,X80,統(tǒng)計結(jié)果表明,牛皮膚病數(shù)據(jù)的樣本均值為1.5,樣本方差為3.417 72,方差比均值為Id=2.278 5.牛皮膚病數(shù)據(jù)樣本路徑如圖1所示,自相關(guān)函數(shù)圖像及偏自相關(guān)函數(shù)圖像分別如圖2和圖3所示.由圖2和圖3可見,該組數(shù)據(jù)為一階相關(guān),因此可以用NGINAR(1)模型對其進行擬合.

圖2 牛皮膚病數(shù)據(jù)自相關(guān)函數(shù)圖像Fig.2 Autocorrelation function plot of skin-lesions data

圖3 牛皮膚病數(shù)據(jù)偏自相關(guān)函數(shù)圖像Fig.3 Partial autocorrelation function plot of skin-lesions data

下面通過序列{Xt}的近似h步(h∈0)預測條件分布方法對NGINAR(1)模型進行預測(簡稱條件分布預測).條件分布預測方法較傳統(tǒng)條件期望預測方法更適用于整數(shù)值時間序列.雖然條件期望預測方法可以使預測值的均方誤差最小,但當觀察值和待預測值為整數(shù)值時,利用條件期望預測方法得到的預測值卻很少取到整數(shù)值點.為了解決上述問題,文獻[10]提出通過條件分布預測方法對整數(shù)值模型進行預測,用這種預測方法得到的預測值和整數(shù)值時間序列本身的狀態(tài)空間一致,而且利用條件分布預測方法來計算條件中位數(shù)、條件均值及條件眾數(shù)等點的預測,甚至預測值的置信區(qū)間都比較容易,能得到較理想的預測值.

由于NGINAR(1)過程具有Markov性,在給定Xn的條件下,Xn+h的條件分布(即Xn的條件預測分布)為

P(Xn+h=xn+h|Xn=xn)=[Ph]xn+h,xn,

其中轉(zhuǎn)移概率為

下面考察NGINAR(1)模型的預測效果.利用牛皮膚病數(shù)據(jù)集對模型進行預測.將條件分布預測和條件期望預測方法進行比較,結(jié)果列于表3.由表3可見,當h=1,2,3,4時,使用條件分布預測方法得到的預測值均為0,與實際值相符,而條件期望預測方法的預測結(jié)果分別為1.412 7,1.470 4,1.478 5,1.479 7,與實際值0有一定的偏差,而且條件分布預測方法的預測均值絕對偏差為0,條件期望預測方法的均值絕對偏差為1.460 3,因此,條件分布預測的方法更適用于整數(shù)值時間序列.

表3 牛皮膚病數(shù)據(jù)的條件期望與條件分布預測結(jié)果比較

圖4為牛皮膚病數(shù)據(jù)的h步條件預測分布圖像.用NGINAR(1)模型、ZIPINAR(1)模型、PINAR(1)模型來擬合該組數(shù)據(jù),并用AIC(Akaike信息準則)、BIC(Bayes信息準則)、均方根值(RMS)和方差比均值對上述3個模型進行比較,結(jié)果列于表4.由表4可見,NGINAR(1)模型具有最小的AIC值和BIC值,3個模型的均方根值相差不大,而NGINAR(1)模型的方差比均值為2.479 9,更接近于數(shù)據(jù)集自身的方差比均值2.278 5,這些數(shù)據(jù)均表明用NGINAR(1)模型擬合該組數(shù)據(jù)集較合適.

圖4 牛皮膚病數(shù)據(jù)的h步條件預測分布Fig.4 h-Step conditional prediction distribution of skin-lesions data

模型估計AICBIC均方根值IdZIPINAR(1)^α=0.164 1276.875 6284.021 71.806 81.674 0^λ=2.031 1^ρ=0.386 3PINAR(1)^α=0.157 3293.339298.103 11.807 51.000 0^λ=1.256 7NGINAR(1)^α=0.140 0271.103275.8671.809 72.479 9^μ=1.479 9

通過上述模擬結(jié)果可知,NGINAR(1)模型參數(shù)的Bayes估計效果優(yōu)于條件最小二乘估計和Yule-Walker估計,且條件分布預測方法比條件期望預測方法更適用于整數(shù)值時間序列.