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

楊艷秋,王德輝

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

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

其中:α∈(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的值.但在實際問題中,每一個事件都可能關聯更多相關的計數事件,因此用幾何隨機變量來描述這些事件更合適.文獻[3]引入了負二項稀疏算子“*”,定義為

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

1 模型簡介

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

Xt=α°Xt-1+Zt.

若Xt表示t時刻住院的患者數,α°Xt-1表示上個月以概率α仍繼續住院的患者數,Zt為t時刻新住院的患者數,則Xt為α°Xt-1與Zt之和.但對于過度分散的計數過程,負二項稀疏算子“*”則更合適.文獻[6]基于負二項稀疏算子“*”,提出了一階整數值自回歸模型(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.根據Bayes原理,參數α的后驗分布為

綜上可得：

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

3 數值模擬

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

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

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

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

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

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

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

4 實例分析及模型預測

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

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

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

圖2 牛皮膚病數據自相關函數圖像Fig.2 Autocorrelation function plot of skin-lesions data

圖3 牛皮膚病數據偏自相關函數圖像Fig.3 Partial autocorrelation function plot of skin-lesions data

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

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

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

其中轉移概率為

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

表3 牛皮膚病數據的條件期望與條件分布預測結果比較

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

圖4 牛皮膚病數據的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

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