廣義二項自回歸模型的貝葉斯統計推斷

邵思雨, 龍 琪, 張 潔

(長春工業大學 數學與統計學院, 吉林 長春 130012)

0 引 言

時間序列數據是指把具有某一個相同統計指標的數值依據它發生的日期先后順序排列而成的數列。而取值為非負整數的時間序列數據作為一種計數數據,其研究范圍更為廣泛。因為這類時間序列數據只能取整值,而傳統實數范圍上的時間序列數據能取到小數,所以這兩種數據的研究存在較大差別。一般情況下,依據非負整數值時間序列數據有無上界,可以將此數據劃分成有上限的數據(0,1,…,n)和無上限的數據(0,1,…)。例如:某一地區一周七天下雨的天數數據,以及某一地區患新型肺炎的人數等。這些年,研究者為了擬合這類數據建立了不少模型。迄今為止,最常見的兩種建模方法分別為基于潛過程構造的狀態空間模型[1],以及基于稀疏算子提出的整數值自回歸模型[2-6]。

針對取值范圍有上限的整數值時間序列數據,例如長春市每天有多少個城區發生交通事故數據、全國每個月有多少城市發生地震數據等。擬合這類時間序列數據較為常見的方法是McKenzie E[7]提出的First-order binomial autoregre-ssive(BAR(1))模型。Weiβ C H[8]研究BAR(1)模型的相應性質,并將其推廣到高階情形。此后,越來越多BAR模型的衍生模型涌現[9-11]。然而,二項稀疏算子是基于獨立同分布的伯努利隨機變量序列提出的,所以,如果將BAR(1)模型應用到個體之間具有相依結構的數據中,可能會使該模型的擬合效果不好。基于此,Ristic M M等[12]提出一個新的算子——廣義二項稀疏算子,它的提出為具有相依結構的整數值時間序列數據的研究提供了參考價值。近期,Kang Y等[13]和康堯[14]將廣義二項稀疏算子應用到二項自回歸模型中,得出一類一階廣義二項自回歸(GBAR(1))模型,并且對該模型進行了相應的統計推斷。由于此算子放松了對二項稀疏算子的獨立性假設,從而對于包含相依結構的整數值時間序列數據的研究具有較大意義,并且能夠有效地解決古典BAR(1)模型解釋實際背景問題時的部分局限性。

對于上述模型的統計推斷研究,模型相應的參數估計問題是研究的重點內容之一。除了常用的Yule-Walker 估計、條件最小二乘估計、條件最大似然估計和經驗似然估計以外,貝葉斯估計也是一種對模型參數進行統計推斷的方法,不少學者對貝葉斯方法進行了研究。Lazar N A等[15]針對貝葉斯經驗似然檢驗了后驗推斷的有效性。張哲等[16]利用貝葉斯方法研究了First-order integer-valued autoregressive(INAR(1))模型參數的估計問題。Xi R等[17]針對spike-and-slab先驗的經驗似然的貝葉斯分位數回歸,研究了其中非參數貝葉斯變量選擇等問題。朱復康等[18]利用矩方法和貝葉斯方法研究Integer-valued generalized autoregressive conditional heteroscedastic(INGARCH (1,1))模型的參數估計問題。近期部分研究詳見文獻[19-20] 。

基于此,考慮針對GBAR(1)模型,研究相應的貝葉斯統計推斷問題,并將該模型應用于實際數據的分析研究中。通過貝葉斯估計方法得到有效的擬合結果,為具有相依結構的實際問題研究提供理論支撐。

1 GBAR(1)模型的貝葉斯估計

Kang Y等[13]將廣義二項稀疏算子 “ °θ”引入到BAR (1)模型中,提出一類一階廣義二項自回歸(GBAR (1))模型,并且給出模型的Yule-Walker估計、條件最小二乘估計、條件最大似然估計。文中通過進一步研究給出了GBAR (1)模型的貝葉斯估計。首先給出GBAR (1)模型的定義及相關性質。

定義1如果{Xt}是有上限的非負整數值時間序列數據,并且{Xt}t∈N滿足方程:

Xt=α°θΧt-1+β°θ(n-Χt-1),t≥1,

(1)

其中n∈N為數據上限,α、θ∈(0,1),各稀疏運算之間相互獨立,則稱{Xt}為廣義二項自回歸模型,記作GBAR(1)模型。

“α°θ”是廣義二項稀疏算子, 其定義為:

Ui=(1-Vi)Wi+ViZ,

式中:{Wi}i∈N,{Vi}i∈N----兩列獨立同分布的Bernoulli(α)和Bernoulli(θ)隨機變量序列;

Z----一個Bernoulli(α)隨機變量。

{Wi}i∈N,{Vi}i∈N和Z之間是相互獨立的,α、θ∈(0,1)。

對于給定的n∈N,GBAR(1)模型的轉移概率為:

P(Xt=k|Xt-1=l)=

(2)

在轉移概率已知的條件下進一步考慮GBAR(1)模型的貝葉斯估計。假設X1,X2,…,XT是來自GBAR(1)模型的T個樣本觀測值,α,β,θ為待估參數。由式(2)可得樣本的聯合似然函數為

P(x|α,β,θ)=L(x|α,β,θ)=

(3)

取參數α,β,θ的先驗分布分別為Beta(e,f)、Beta(p,q)和Beta(h,g),這里我們分別記作π(α),π(β),π(θ),即

(4)

(5)

(6)

由此得出聯合密度函數

P(x,α,β,θ)=L(x|α,β,θ)×π(α)×

π(β)×π(θ)。

(7)

在給定聯合似然函數以及先驗分布后,由貝葉斯原理可知,參數α,β,θ的聯合后驗為

(8)

(9)

(10)

(11)

2 數值模擬

對GBAR(1)模型進行數值模擬,通過參數估計的偏差和均方誤差來比較Conditional Least Squares (CLS) 估計、Conditional Maximum Likelihood (CML) 估計和Bayesian (BYE)估計的估計效果。根據參數α,β,θ的取值空間,選取以下四組不同的參數組合進行數值模擬:

1)n=5,α=0.25,β=0.75,θ=0.65;

2)n=5,α=0.40,β=0.60,θ=0.55;

3)n=8,α=0.40,β=0.20,θ=0.45;

4)n=8,α=0.80,β=0.20,θ=0.35。

選取參數α,β,θ的先驗分布為貝塔分布。令貝塔分布的均值等于參數真值的方法來選取超參數,并令幾組參數組合的先驗分布分別為:

1)α～Beta(1,3),β～Beta(6,2),θ～Beta(2,1);

2)α～Beta(2,3),β～Beta(6,4),θ～Beta(2,2);

3)α～Beta(2,3),β～Beta(1,4),θ～Beta(1,1);

4)α～Beta(4,1),β～Beta(1,4),θ～Beta(1,2)。

鑒于在進行貝葉斯估計時,需要計算多重積分,文中采用樣本平均值法求解數值積分。此方法是利用樣本的均值來估計所求的定積分

若Y服從(c,d)上的均勻分布,則有

得出定積分

I=E(g(Y))×(d-c)。

因此,當序列{Yt,t=1,2,…,N}獨立同分布于(c,d)上均勻分布,則

Gi=g(Yi),i=1,2,…,N

同為一列獨立同分布的隨機變量序列,根據強大數定律,則有

所以

是I的強相合無偏估計,如此求解定積分的方法稱為樣本平均值法。

運用此法求積分

能將積分求解問題轉化為求g(Y)期望的問題。在數值模擬時,樣本平均值法對應的樣本量N取值為105。

模擬考慮在觀測數據樣本量T分別為20,50,100 和200的情況下進行,最終結果是在R軟件環境下取 500 次估計結果的平均值。四個參數組合的樣本路徑和ACF如圖 1所示。

從圖1可以看出, 序列為平穩時間序列,并且存在短期自相關性。

(a) 樣本路徑 (b) ACF

GBAR(1)模型在四組參數下CLS 估計、CML估計和BYE估計的模擬結果分別見表1～表4。

表1 第1組參數的估計結果

表2 第2組參數的估計結果

表3 第3組參數的估計結果

表4 第4組參數的估計結果

為了展現三種估計效果的好壞,我們用偏差 (bias) 和均方誤差(MSE)兩種度量指標。通過表中模擬結果可以看出,CLS、CML和BYE三種估計的偏差和均方誤差都比較小,這表明以上三種估計方法的效果都比較合適。并且隨著樣本量的增大, 估計量的偏差和均方誤差都在減小,估計效果越來越好。其次,在樣本量相對比較小時,貝葉斯估計的偏差和均方誤差比條件最大似然估計的略小一些,而當樣本量增加后,貝葉斯估計的效果和條件最大似然估計的效果相差也不大,且貝葉斯估計的變化幅度相對較小, 后續我們通過對數據污染情況下參數估計的比較更是驗證了貝葉斯估計的穩健性。

為了對三種估計方法做進一步的比較,接下來考慮三種估計方法對帶污染的GBAR(1)模型的估計效果。

定義2若一個隨機過程{Mt}t∈N滿足方程

Mt=(1-δt)Xt+δtξt,t∈N,

式中:{Xt}t∈N----一個沒有被污染的GBAR(1) 過程;

{ξt}t∈N,{δt}t∈N----一列i.i.d.的隨機變量序列,分別服從B(n,p)和B(1,r);

r----污染比例,該模型為帶污染的GBAR(1)模型。

樣本量為100和200時,污染比例r=0.3,0.5時GBAR(1)模型的三種估計方法模擬結果分別見表5～表6。

從表中可知,對于有污染的數據,貝葉斯估計的效果要優于CML和CLS估計。由此可知,當數據存在污染時,貝葉斯估計方法更為可靠。進而說明,貝葉斯估計更加穩健,估計結果合理有效。

3 實例分析

實際生活中,有上限范圍的整數值時間序列數據是很常見的,例如:某一地區一周七天下雨天數數據、全國每個月有多少城市發生地震數據等。文中從天氣網中查詢了2020年1月1日至2021年12月31日湖南省長沙市的歷史天氣,然后統計其每周下雨天數的數據情況,通過整理得到了有上限范圍的時間序列數據,其中上限范圍n=7,樣本量為105,數據的相應統計信息見表7 。

表7 基本統計量

周降雨天數數據的樣本路徑和自相關圖如圖2所示。

(a) 樣本路徑 (b) ACF

從圖2(a)可以看出,數據沒有向上或向下的趨勢,而是在 3 附近波動,說明數據是平穩的。觀察圖2(b)可以看出,自相關系數在2階后衰減到了標準差以內,也能確定序列是平穩時間序列。并且序列具有自相關性,因此用上限n=7的GBAR(1)模型來刻畫這組數據是合理的。

考慮將GBAR(1)模型應用于該數據,并與文獻[7] 提出的BAR(1)模型、文獻[21]提出的衍生模型、BAR(1)、(IZ-BAR(1))模型和(ZIB-AR(1))模型進行比較。其中,GBAR(1)模型的貝葉斯估計中超參數設置為e=5,f=5,p=2,q=3,h=6,g=4,對應樣本平均值法的N取106。

通過將GBAR(1)模型與各種其他模型參數的赤池信息準則(AIC)和貝葉斯信息準則(BIC)進行對比,用來說明GBAR(1)模型擬合的優良性。

一方面,無論是GBAR(1)模型的貝葉斯估計還是條件最大似然估計,其AIC、BIC值都小于BAR(1)模型、IZ-BAR(1)模型和ZIB-AR(1)模型,從而得出GBAR(1)模型的擬合效果優于BAR(1)模型以及其衍生模型;另一方面,GBAR(1)模型的貝葉斯估計對于樣本量接近100的數據估計結果和AIC、BIC值與條件最大似然估計的相差很小,說明應用貝葉斯方法估計模型參數也能得到一個比較好的擬合效果,GBAR(1) 模型與其他模型AIC、BIC比較見表8。

表8 GBAR(1) 模型與其他模型AIC、BIC比較

進一步對GBAR(1)模型的皮爾遜殘差進行分析,其中皮爾遜殘差公式為

根據GBAR(1)模型的條件期望和條件方差公式,我們分別計算了GBAR(1)模型在貝葉斯和條件最大似然兩種估計方法下殘差的均值和方差。認為擬合殘差的均值和方差越接近0和1,模型的擬合效果越好。在本節實例數據中,所擬合的貝葉斯方法下GBAR(1)模型殘差的均值和方差分別為0.025 6和1.293 3。該結果也再次說明GBAR(1)模型擬合這組數據是合適的。

4 結 語

為了進一步研究具有相依結構的廣義二項自回歸模型的參數估計問題,充分利用樣本信息、總體信息,以及通過歷史經驗或者規律得到的參數的先驗信息,求得GBAR(1)模型參數的貝葉斯估計表示式。證明了對于GBAR(1)模型在小樣本情況下,貝葉斯估計的效果要優于條件最大似然估計和條件最小二乘估計。通過對帶污染的GBAR(1)模型進行數值模擬,驗證了貝葉斯估計方法的穩健性。最后將GBAR(1)模型應用于長沙市下雨天數據集進行擬合,與幾類二項自回歸模型進行比較,得到最小的AIC與BIC值。表明GBAR(1)模型的擬合效果優于BAR(1)模型以及其衍生模型。