基于ABC-SMC回歸算法具有觀測噪聲的AR(p)模型的參數(shù)估計(jì)

許淑婷，鄭斌斌，李安水，張慧增

(杭州師范大學(xué)數(shù)學(xué)學(xué)院,浙江 杭州 311121)

時(shí)間序列自回歸(AR)模型參數(shù)估計(jì)是統(tǒng)計(jì)學(xué)中一個(gè)重要的研究課題,其估計(jì)方法已經(jīng)日趨成熟,主要有格點(diǎn)搜索、Newton-Raphson方法、最陡爬坡法、極大似然估計(jì)等[1]方法以及貝葉斯統(tǒng)計(jì)框架下的馬爾科夫鏈蒙特卡洛(MCMC)方法和近似貝葉斯計(jì)算(ABC)算法[2].然而當(dāng)觀測數(shù)據(jù)存在噪聲時(shí)，如何對模型參數(shù)進(jìn)行有效的估計(jì)更是具有理論意義和應(yīng)用價(jià)值.

對于具有觀測噪聲的AR模型,1979年Sakai和Arase[3]使用補(bǔ)償最小二乘法(CBLS)對其參數(shù)進(jìn)行估計(jì),在此基礎(chǔ)上,鄭衛(wèi)星[4]提出了一種具有直接結(jié)構(gòu)的(ILSD)改進(jìn)最小二乘法, Mahmoudi和Karimi[5]提出了一種將逆濾波方程與Yule-Walker方程結(jié)合的方法(IFILS).2008年，Diversi等[6]利用誤差不變(EIV)的概念來估計(jì)具有觀測噪聲的AR模型.2020年，Esfandiari等[7]提出了4種方法來估計(jì)具有觀測噪聲的AR模型：第一種方法旨在減少低階Yule-Walker方程中觀測噪聲方差的破壞性影響,并以遞歸的方式估計(jì)參數(shù);第二種方法將其看成約束最小二乘優(yōu)化問題并找到最佳值,通過迭代算法確定觀測噪聲方差,然后估計(jì)參數(shù);第三種是使用近似值,將未知參數(shù)的數(shù)量減少到僅兩個(gè)參數(shù),并以遞歸方式估計(jì)參數(shù);最后一種是將觀測噪聲的方差估計(jì)為充分放大的自相關(guān)矩陣的最小特征值,并使用它來估計(jì)參數(shù).以上方法或計(jì)算成本大、不易實(shí)現(xiàn)，或是在高階情況下估計(jì)效果不佳.本文在貝葉斯統(tǒng)計(jì)的框架下,采用簡單又易于實(shí)現(xiàn)的近似貝葉斯序貫蒙特卡洛(ABC-SMC)回歸算法來估計(jì)具有觀測噪聲的AR模型的參數(shù).

近似貝葉斯計(jì)算(ABC)[8]是近年來貝葉斯統(tǒng)計(jì)推斷中常用的方法.當(dāng)模型比較復(fù)雜,似然函數(shù)難以寫出具體表達(dá)形式時(shí),利用ABC算法模擬參數(shù)近似后驗(yàn)分布的樣本進(jìn)而得到參數(shù)估計(jì).1984 年Rubin[9]提出了ABC基本思想,2002 年Beaumont等[10]提出了最基本的ABC拒絕采樣算法,基本ABC算法的缺點(diǎn)是當(dāng)待估參數(shù)的維數(shù)比較大時(shí)計(jì)算成本很大,并且當(dāng)模型具有噪聲時(shí),可能導(dǎo)致估計(jì)效率低下.2007年Sisson等[11]首次將粒子濾波器融入ABC方法,粒子濾波器[12]能從不完整的或者是含有噪聲的觀測序列中估計(jì)出動態(tài)系統(tǒng)的狀態(tài),從而提升估計(jì)效率.2009年Toni等[13]在處理動態(tài)系統(tǒng)問題時(shí)引入了序列重要性采樣(SIS),提出ABC序貫蒙特卡洛算法.然而對于不同的采樣粒子而言,生成的模擬數(shù)據(jù)集與觀測數(shù)據(jù)集的差異程度是不同的,為此我們引入回歸模型對采樣結(jié)果進(jìn)行調(diào)整,形成ABC-SMC回歸算法.

本文內(nèi)容組織如下:第1節(jié)介紹了具有觀測噪聲的AR模型及其自協(xié)方差函數(shù);第2節(jié)闡述了具有觀測噪聲的AR模型的ABC-SMC算法.最后通過數(shù)值模擬,并與ABC算法、ABC-SMC算法進(jìn)行了比較.模擬結(jié)果表明,所提出的ABC-SMC回歸算法顯著提高了具有觀測噪聲的AR模型參數(shù)估計(jì)精度.

1 具有觀測噪聲的AR(p)模型

定義1設(shè){Yt:t∈T}為一時(shí)間序列,其中T=0,±1,±2,…,如果對于?t∈T,有

其中{εt:t∈T}和{νt:t∈T}為均值為0，方差分別為σ和ν的高斯白噪聲序列,且{Xt}和{νt}相互獨(dú)立，則稱{Yt}為具有觀測噪聲的AR(p)模型.

1.1 具有觀測噪聲的AR(p)模型的自協(xié)方差函數(shù)

命題1假設(shè){Xt:t∈T} 為寬平穩(wěn)的中心化AR(p)過程,令γX(0)=Var(Xt),γX(h)=Cov(Xt,Xt-h),h=±1,±2,…,則

1){Yt:t∈T}為寬平穩(wěn)過程;

2)令γY(0)=Var(Yt),γY(h)=Cov(Yt,Yt-h),h=±1,±2,…,則γY(0)=γX(0)+ν,γY(h)=γX(h),h=±1,±2,….

證明由{Xt:t∈T}是寬平穩(wěn)的中心化AR(p)過程,可知E(Xt)=0,則

E(Yt)=E(Xt+vt)=E(Xt)+E(vt)=0.

(1)

由{Xt}和{νt}相互獨(dú)立,則對于?t,τ∈T,t≠τ,有E(Xtvτ)=E(Xt)E(ντ)=0.由{νt}是高斯白噪聲序列,可知

因此，{Yt:t∈T}的方差為

Var(Xt)+ν=γX(0)+ν,

(2)

并且當(dāng)h=±1,±2,…時(shí),有

Cov(Yt,Yt-h)=E(YtYt-h)=E(Xt+vt)(Xt-h+vt-h)=E(XtXt-h+Xtvt-h+Xt-hvt+vtvt-h)=

E(XtXt-h)=Cov(Xt,Xt-h)=γX(h).

(3)

由式(1)、(2)、(3)可知,{Yt:t∈T}是寬平穩(wěn)的.由式(2)可知,γY(0)=γX(0)+ν;由式(3)可知,γY(h)=γX(h)，h=±1,±2,….

注以下討論的 {Xt:t∈T} 均為寬平穩(wěn)的中心化AR(p)過程.

命題2對于具有觀測噪聲的AR(p)模型 {Yt:t∈T},自協(xié)方差函數(shù)滿足:

1)γY(0)=φ1γY(1)+φ2γY(2)+…+φpγY(p)+σ+ν;

2)γY(h)=φ1γY(h-1)+φ2γY(h-2)+…+φpγY(h-p)-φhν,h=1,2,…,p;

3)γY(h)=φ1γY(h-1)+φ2γY(h-2)+…+φpγY(h-p),h=p+1,p+2,….

證明由{Xt:t∈T}為寬平穩(wěn)的中心化AR(p)模型可知,{Xt:t∈T}的自協(xié)方差函數(shù)滿足以下關(guān)系:

由命題1可得

γY(0)=γX(0)+ν=φ1γX(1)+φ2γX(2)+…+φpγX(p)+σ+ν=

φ1γY(1)+φ2γY(2)+…+φpγY(p)+σ+ν.

當(dāng)h=1,2,…,p時(shí),

γY(h)=γX(h)=φ1γX(h-1)+φ2γX(h-2)+…+φpγX(h-p)=

φ1γX(h-1)+φ2γX(h-2)+…+φhγX(0)+φh+1γX(-1)+…+φpγX(h-p)=

φ1γY(h-1)+φ2γY(h-2)+…+φh(γY(0)-ν)+φh+1γY(-1)+…+φpγY(h-p)=

φ1γY(h-1)+φ2γY(h-2)+…+φpγY(h-p)-φhν.

當(dāng)h=p+1,p+2,… 時(shí),

γY(h)=γX(h)=φ1γX(h-1)+φ2γX(h-2)+…+φpγX(h-p)=

φ1γY(h-1)+φ2γY(h-2)+…+φpγY(h-p).

1.2 具有觀測噪聲的貝葉斯AR(p)模型

1.2.1 具有觀測噪聲的AR(p)模型的似然函數(shù)

設(shè){Yt:t∈T}為一具有觀測噪聲的AR(p)模型,記 Φ=(φ1,…,φp) 為模型系數(shù).令y1:n為{Yt:t∈T}的可觀測樣本,x1:n為{Xt:t∈T}的不可觀測樣本,其中1:n表示1,2,…,n.那么y1:n,x1:n的似然函數(shù)為

p(y1:n,x1:n|Φ,σ,ν)=p(y1:n|x1:n,Φ,σ,ν)p(x1:n|Φ,σ,ν),

樣本y1:n的似然函數(shù)為

1.2.2 具有觀測噪聲的AR(p)的參數(shù)先驗(yàn)分布

記I(λ)=λp-φ1λp-1-…-φp=0 是AR(p)模型 {Xt:t∈T} 的特征方程, Λ=(λ1,…,λp) 為該特征方程的特征根.眾所周知, AR(p) 模型平穩(wěn)當(dāng)且僅當(dāng) |λi|<1,i=1,2,…,p.

令特征多項(xiàng)式因式分解如下:

防治適期掌握在發(fā)生高峰期前，且田間若蟲占總蟲量80%以上。無公害茶園用高效、低毒、低殘留的2.5%高效氯氟氰菊酯微乳劑 1 500 克/公頃、25%噻蟲嗪水分散粒劑50克/公頃進(jìn)行防治；有機(jī)茶園則選用生物藥劑0.3%苦參堿水劑4.5～9.0克/公頃[4]、30%茶皂素水劑30%茶皂素水劑60毫升/公頃[5]。采用低容量蓬面噴霧施藥。

λp-φ1λp-1-φ2λp-2-…-φp=(λ-λ1)(λ-λ2)…(λ-λp).

(4)

式(4)右邊展開得到

與式(4)左邊比較有如下關(guān)系:

由此可知, AR(p)模型的模型系數(shù)φ1,φ2,…,φp由特征根λ1,λ2,…,λp唯一確定.

為了書寫方便,我們將模型系數(shù)與特征根的關(guān)系記為 Φ=f(Λ).

由于在高階情況下,難以刻畫 AR(p) 模型的平穩(wěn)域,可由根與系數(shù)的關(guān)系,將對模型系數(shù)Φ的估計(jì)轉(zhuǎn)換成對特征根 Λ 的估計(jì).因此在貝葉斯統(tǒng)計(jì)的框架下,具有觀測噪聲的AR(p)模型的待估參數(shù)為Λ,σ和ν.假設(shè)Λ,σ,ν相互獨(dú)立,通常我們選取σ的先驗(yàn)分布為逆Gamma分布IG(ασ,βσ),ν的先驗(yàn)分布為逆Gamma分布IG(αν,βν),其中超參數(shù)ασ,βσ,αν,βν均為已知.

特征方程I(Λ)=0既有實(shí)根也有復(fù)根,且復(fù)根以共軛的形式成對出現(xiàn).若{Xt:t∈T}平穩(wěn),則特征根的模小于1.若特征根的先驗(yàn)分布為復(fù)平面單位圓上連續(xù)型分布，則特征根為實(shí)根的概率為0,顯然這是不合理的.本文對于特征根先驗(yàn)分布構(gòu)造思想如下:首先定義共軛復(fù)根的對數(shù)為均勻分布,在共軛復(fù)數(shù)對數(shù)已知的條件下,遵循先驗(yàn)分布簡單的原則,先在上半單位圓內(nèi)按照均勻分布取復(fù)根,并確定其共軛復(fù)根,然后在區(qū)間(-1,1)上按照均勻分布對實(shí)根進(jìn)行抽樣.下面我們給出Λ先驗(yàn)分布的表示.

當(dāng)k=0時(shí),令λi～U(-1,1),i=1,2,…,p,且λ1,λ2,…,λp相互獨(dú)立.則有Λ=λ1:p的條件分布如下:

故參數(shù) Λ 的先驗(yàn)分布為

1.2.3 具有觀測噪聲的AR(p)模型的參數(shù)后驗(yàn)分布

假設(shè)具有觀測噪聲的AR(p)模型的參數(shù)先驗(yàn)分布為π(Λ),π(σ),π(ν).由1.2.1有y1:n,x1:n的似然函數(shù)為

p(y1:n,x1:n|Λ,σ,ν)=p(y1:n|x1:n,f(Λ),σ,ν)p(x1:n|f(Λ),σ,ν),

樣本y1:n的似然函數(shù)為

則有模型的聯(lián)合分布為

p(y1:n,Λ,σ,ν)=p(y1:n|Λ,σ,ν)π(Λ)π(σ)π(ν).

由貝葉斯公式可知,具有觀測噪聲的AR(p)模型的參數(shù)后驗(yàn)分布為

由于模型參數(shù)的后驗(yàn)分布無法給出具體的表達(dá)形式,下面我們通過給出參數(shù)后驗(yàn)分布的ABC-SMC抽樣算法和ABC-SMC回歸抽樣算法,給出后驗(yàn)分布的估計(jì).

2 具有觀測噪聲的AR(p)模型的貝葉斯估計(jì)

在本節(jié)中我們首先使用ABC-SMC算法對參數(shù)后驗(yàn)分布進(jìn)行抽樣,然后引入回歸模型對抽樣進(jìn)行調(diào)整,本文把該算法稱為ABC-SMC回歸算法.

2.1 具有觀測噪聲的AR(p)模型的ABC-SMC算法

算法1Λ 的先驗(yàn)分布抽樣算法

輸入:模型的階數(shù)p.

輸出:生成Λ的隨機(jī)向量.

算法步驟:

1.生成k～U{0,1,2,…,p2}.2.if k=0 then3.生成λi～U(-1,1),i=1,2,…,p.4.end if5.if k>0 then 6.生成λ'i～U(D),i=1,2,…,k.7.令λ2i-1=λ'i,i=1,2,…,k.8.令λ2i= λ'i,i=1,2,…,k.9. end if 10.生成λi～U(-1,1),i=2k+1,…,p. 11.輸出Λ=λ1:p.

ABC-SMC算法基于蒙特卡洛方法,引入了序貫重要性采樣方法.在ABC-SMC算法中,每次采樣不再是一個(gè)樣本，而是預(yù)定數(shù)量的樣本N,每個(gè)樣本稱為粒子,一次采樣的樣本序列稱為“粒子池”.對于粒子池中的每個(gè)粒子,設(shè)定擾動核q(·|·)進(jìn)行擾動,再判斷其是否滿足所設(shè)定的條件,以此選擇接受或拒絕.對接受的粒子進(jìn)行賦權(quán),以此權(quán)重進(jìn)行下一次采樣,生成新一代的粒子池.

對于具有觀測噪聲的AR(p)模型,由于Λ,σ,ν,相互獨(dú)立,因此模型參數(shù)的先驗(yàn)分布為

π(θ)=π(Λ)π(σ)π(ν).

算法2Λ的擾動算法

輸入:當(dāng)前參數(shù)Λ(t)=λ1:p;當(dāng)前共軛復(fù)數(shù)對數(shù)k.

輸出:擾動后的參數(shù)Λ*.

算法步驟:

1.ifk=0 then

2.k*=1.

5.生成λ*～U(D).

7.end if

9.生成k*～U{k-1,k+1}.

10.ifk*=k-1 then

11.抽取I～U{1,2,…,k}.

14.end if

15.ifk*=k+1 then

18.生成λ*～U(D).

20.end if

21.end if

24.抽取I～U{1,2,…,k}.

26.ifp為奇數(shù) then

28.else

30.end if

31.end if

令k為Λ的共軛復(fù)數(shù)的對數(shù),k*為Λ*的共軛復(fù)數(shù)的對數(shù).由算法2可知,當(dāng)k*=k-1 時(shí),

當(dāng)k*=k+1時(shí)，

設(shè)參數(shù)σ的擾動核為q(σ*|σ,ησ)=N(σ,ησ).參數(shù)ν的擾動核為q(ν*|ν,ην)=N(ν,ην).為了保證擾動得到的σ,ν均大于等于0,先擾動得到σ*和ν*,直至滿足條件為止.下面是參數(shù)σ,ν的擾動算法.

算法3σ,ν的擾動算法

輸入:當(dāng)前參數(shù)σ(t),ν(t).

輸出:擾動得到的下一代參數(shù)σ(t+1),ν(t+1).

算法步驟:

1.repeat 2.σ*～q(·|σ(t),ησ).3.until σ*>04.repeat5.ν*～q(·|ν(t),ην). 6.until ν*>07.輸出 σ(t+1)=σ*,ν(t+1)=ν*.

令θ*=(Λ*,σ*,ν*),θ=(Λ,σ,ν).由算法2和算法3可得

q(θ*|θ)=q(Λ*|Λ)q(σ*|σ,ησ)q(ν*|ν,ην).

(5)

下面給出具有觀測噪聲的AR(p)模型的ABC-SMC算法.

算法4ABC-SMC算法

輸入:觀測數(shù)據(jù)集y∈Rn;模型階數(shù)p;先驗(yàn)分布π(θ);超參數(shù)ασ,βσ,αν,βν;距離函數(shù)d(·,·);擾動核函數(shù)q(·|·);容忍度δ>0;粒子池?cái)?shù)目N;迭代次數(shù)T.

算法步驟:

2.fori=1 toNdo

4.end for

6.fort=1 toTdo

7.fori=1 toNdo

8.repeat

10.ifπ(θ**)≠0 then

12.end if

13.untild(s,S*)≤δ

15.end for

17.end for

2.2 具有觀測噪聲的AR(p)模型的ABC-SMC回歸算法

(6)

算法5ABC-SMC回歸算法

輸入:觀測數(shù)據(jù)集y∈Rn;模型階數(shù)p;先驗(yàn)分布π(θ);超參數(shù)ασ,βσ,αν,βν;Rp+1上距離函數(shù)d(·,·);擾動核函數(shù)q(·|·);容忍度δ>0;粒子池?cái)?shù)目N;迭代次數(shù)T.

算法步驟:

3.由式(6)定義X,Θ.

7.fori=1 toNdo

9.end for

3 數(shù)值模擬

為了說明以上方法的有效性,下面給出如下數(shù)值模擬.

假定模型為

先由模型隨機(jī)地生成觀測數(shù)據(jù)集y1,…,yn;取樣本量n=1 000,超參數(shù)ασ=1,βσ=0.3,αν=1.2,βν=0.5,粒子池?cái)?shù)目為N=100,迭代次數(shù)T=1 000,容忍度δ=3,參數(shù)的擾動核方差ησ=ην=0.1;分別使用ABC-SMC算法以及ABC-SMC回歸算法對參數(shù)后驗(yàn)分布進(jìn)行抽樣,計(jì)算樣本均值作為參數(shù)的估計(jì)值,并與ABC算法進(jìn)行比較.為了觀察先驗(yàn)分布選取對模擬結(jié)果的影響,再取先驗(yàn)分布σ～Exp(0.5),ν～Exp(0.8)進(jìn)行模擬,并與前面的估計(jì)結(jié)果比較.

以下圖1、 圖2和 圖3分別為ABC算法、 ABC-SMC算法以及ABC-SMC回歸算法抽樣所得的參數(shù)近似后驗(yàn)分布結(jié)果,實(shí)曲線為核密度估計(jì),實(shí)直線表示真實(shí)值,虛線表示估計(jì)值.其估計(jì)值以及估計(jì)相對誤差的比較見表1.由表1可以看出,ABC-SMC算法的估計(jì)精度比ABC算法高很多,而引入回歸之后的ABC-SMC回歸算法估計(jì)效果有顯著提升.表2為不同先驗(yàn)分布的模擬結(jié)果.由表2可以看出,不同的先驗(yàn)分布均能將估計(jì)誤差控制在10%以內(nèi),說明先驗(yàn)分布的選取對模擬結(jié)果的影響不大.

φ1 φ2 σ υ

φ1 φ2 σ υ

φ1 φ2 σ υ

表1 ABC算法、ABC-SMC算法和ABC-SMC回歸算法估計(jì)結(jié)果對比

表2 ABC-SMC回歸算法選取不同先驗(yàn)分布的估計(jì)結(jié)果對比

4 結(jié)論

對于具有觀測噪聲的AR(p)模型,采用前p階樣本自協(xié)方差函數(shù)作為統(tǒng)計(jì)量,使用ABC-SMC算法對其參數(shù)進(jìn)行估計(jì)是有效的;在此基礎(chǔ)上引入回歸模型調(diào)整采樣結(jié)果后,估計(jì)精度有顯著提升.