楊艷秋,王德輝

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

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

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

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

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

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

1 模型簡介

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

Xt=α°Xt-1+Zt.

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

Xt=α*Xt-1+εt,

其中: 負(fù)二項(xiàng)稀疏算子“*”定義為

2 Bayes估計(jì)

其中

i=1,2,…,T.

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

綜上可得：

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

3 數(shù)值模擬

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

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

5) 令i=i+1,返回步驟2),直到算法達(dá)到事先約定的收斂標(biāo)準(zhǔn).

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

表1 Bayes估計(jì)的偏差和均方誤差

表2 3種估計(jì)方法的均方誤差比

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

4 實(shí)例分析及模型預(yù)測

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

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

將數(shù)據(jù)集的前80個觀測數(shù)據(jù)記作X1,X2,…,X80,統(tǒng)計(jì)結(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)模型對其進(jìn)行擬合.

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

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

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

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

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

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

圖4為牛皮膚病數(shù)據(jù)的h步條件預(yù)測分布圖像.用NGINAR(1)模型、ZIPINAR(1)模型、PINAR(1)模型來擬合該組數(shù)據(jù),并用AIC(Akaike信息準(zhǔn)則)、BIC(Bayes信息準(zhǔn)則)、均方根值(RMS)和方差比均值對上述3個模型進(jìn)行比較,結(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步條件預(yù)測分布Fig.4 h-Step conditional prediction distribution of skin-lesions data

模型估計(jì)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估計(jì)效果優(yōu)于條件最小二乘估計(jì)和Yule-Walker估計(jì),且條件分布預(yù)測方法比條件期望預(yù)測方法更適用于整數(shù)值時間序列.