Markov轉(zhuǎn)換向量自回歸模型的變分貝葉斯推斷

呂秀麗，虞棟杰，鄭靜(杭州電子科技大學經(jīng)濟學院，浙江杭州310018)0引言研究變量間相互依賴關(guān)系的時間序列是較為經(jīng)典的統(tǒng)計分析方法，其中向量自回歸模型(VectorAutoregressiveModels,VAR)是常用的計量經(jīng)濟模型。使用VAR模型對時間序列進行分析時，需先假設(shè)時間序列線性同質(zhì)，但在實際應(yīng)用中，各種復雜因素致使線性假定成立模型的分析結(jié)果不夠準確。為了解決這個問題，學者們將研究方向轉(zhuǎn)向非線性模型，龍如銀等[1]運用區(qū)制為3、滯后階數(shù)為1的Markov區(qū)制轉(zhuǎn)換模型(MarkovSwitchingModels,MS)研究我國1984年至2003年的通貨膨脹率，與傳統(tǒng)VAR模型相比，更好地擬合了我國經(jīng)濟的通貨膨脹率。劉金全等[2]運用Markov轉(zhuǎn)換向量自回歸(MarkovSwitching-VectorAutoregressiveModels,MS-VAR)模型將我國通貨膨脹率分為3個區(qū)制，分析經(jīng)濟政策與通貨膨脹率在各區(qū)制下的關(guān)系，合理解釋了當時我國社會經(jīng)濟及通貨膨脹風險狀況，指出治理國家通貨膨脹應(yīng)先判斷其所處區(qū)制。朱慧明等[3]引入貝葉斯MS-VAR模型，2次調(diào)用Gibbs抽樣方法，研究我國股市與通貨膨脹率的波動關(guān)系，指出可通過機制轉(zhuǎn)換模型來刻畫該波動關(guān)系。以上模型的參數(shù)較多，需事先進行參數(shù)估計。目前常用的模型參數(shù)估計方法有2類，一是確定性方法，如極大似然估計、EM算法等，在估計過程中未考慮模型參數(shù)的先驗信息；二是貝葉斯方法、蒙特卡羅法等，計算成本較高，易產(chǎn)生過擬合或欠擬合。變分貝葉斯(VariationalBayes,VB)是貝葉斯推理的近似方法。Adami[4]將變分方法應(yīng)用到貝葉斯推理中，介紹了如何用變分方法來逼近這類推理中出現(xiàn)的難以處理的后驗分布；Alan等[5]采用參數(shù)擴展變分貝葉斯方法解決了VB收斂到近似解較慢的問題；Ormerod等[6]將變分近似引入統(tǒng)計學，用統(tǒng)計學術(shù)語解釋變分近似概念；Koop等[7]提出一種均值場變分貝葉斯算法，分析美國宏觀經(jīng)濟數(shù)據(jù)，將其拓展至計量經(jīng)濟學領(lǐng)域；馬占宇等[8]結(jié)合變分推斷，提出一種學術(shù)研究熱點Key提取方法，高效、準確地提取學術(shù)研究的熱點Key；沈圓圓等[9]提出一種Logistic組稀疏性回歸的Bayes概率模型的變分貝葉斯算法，提高了模型的預測性能；萬志成等[10]建立了一種無限變量高斯混合模型，解決了無限變量高斯混合模型的參數(shù)估計問題；付維明等[11]運用擴散方法對分布式隨機變分推斷算法進行改進，并將其應(yīng)用至主題模型，速度快，可擴展強。本文提出一種Markov轉(zhuǎn)換向量自回歸模型的變分貝葉斯推斷算法，利用參數(shù)的先驗信息對MS-VAR模型的參數(shù)進行估計，并給出完整的變分推理過程。1MS-VAR模型1.1馬爾可夫轉(zhuǎn)移模型給定所有先前的狀態(tài)，馬爾可夫鏈(本文簡稱馬氏鏈)的當前狀態(tài)依賴且僅依賴其前一個狀態(tài)，這一性質(zhì)被稱為馬爾可夫特性。馬爾可夫區(qū)制轉(zhuǎn)移模型依賴具有馬爾可夫性的馬氏鏈生成的隨機過程。設(shè){yt,t=1,…,T}為具有T個單變量觀測值的時間序列，是隨機過程{yt}的子集。隱藏的離散隨機過程{zt}的實現(xiàn)決定了隨機過程{yt}的概率分布。{yt}是直接可觀測的，而{zt}是一個潛在的隨機變量，只能通過觀測{yt}來間接觀測過程{zt}。假設(shè)隱過程{zt}是一個具有有限狀態(tài)空間{0,…,K-1}的不可約、非周期的馬氏鏈，其隨機性質(zhì)用K×K階轉(zhuǎn)移矩陣A描述，其中aij是轉(zhuǎn)移矩陣A的第i+1行第j+1列元素，表示在t-1時刻處于狀態(tài)i，下一時刻轉(zhuǎn)移到狀態(tài)j的概率。轉(zhuǎn)移矩陣A中每一行的所有元素之和為1且所有元素非負，即：設(shè)參數(shù)θi∈Θ是由轉(zhuǎn)移區(qū)制決定的參數(shù)，而ψ∈Ψ是不依賴于區(qū)制的參數(shù)，Θ,Ψ為參數(shù)空間，θi和ψ共同決定了參數(shù)條件密度p(·|θi,ψ|)。本文假設(shè)在給定z1,…,zΤ時的隨機變量y1,…,yΤ是條件獨立的，其密度為：yt|(zt=i)|～p(·|θi,ψ|)對于聯(lián)合隨機過程{zt,yt}，yt的條件密度為：其中，It-1={yt-1,yt-2,…}為t-1時刻可獲得的觀測信息。1.2馬爾可夫區(qū)制轉(zhuǎn)移的向量自回歸模型假設(shè)yt=(y1t,y2t,…,ymt)′為m×1的觀察值向量，隨機變量st=1,2,…,J(J≥2)表示在t時刻的狀態(tài)，與誤差向量εst獨立，誤差向量εst服從均值為0且協(xié)方差陣為Σst的正態(tài)分布，存在一個初始狀態(tài)概率向量為Π、狀態(tài)轉(zhuǎn)移概率矩陣為A的馬氏鏈，St={s1,s2,…,st}是由該馬氏鏈生成的不可觀測的隨機狀態(tài)序列，Π=(πj),πj=p(s1=j),1≤j≤J；A=(aij),aij=p(st=j|st-1=i|),1≤i,j≤J,1≤t≤T其中，πj為初始狀態(tài)為j的概率。Ust為st狀態(tài)下的均值向量，故MS-VAR模型的向量形式為：其中，yt∈Rm，α1,st,…,αd,st∈Rm×m表示自回歸系數(shù)，d為模型的滯后階數(shù)，εst∈Rm。設(shè)Ωst=[Ust,α1,st,…,αd,st]′∈R(md+1)×m，xt=[1,yt-1,…,yt-d]′∈R(md+1)×1，故MS-VAR模型改寫成為：yt=Ω′stxt+εst進一步轉(zhuǎn)化為：Y≈XΩst根據(jù)以上參數(shù)的相關(guān)假設(shè)可以得到y(tǒng)t～Nm×1(Ω′stxt,Σst)，則向量yt的條件概率密度分布函數(shù)為：其中，Φt-1為變量參數(shù)集合，包含t-1時刻所有參數(shù)的信息。2變分貝葉斯推斷貝葉斯估計方法通過設(shè)置參數(shù)的先驗分布，根據(jù)貝葉斯規(guī)則，利用樣本信息推導出參數(shù)后驗分布，得到參數(shù)估計值。但是，由于大部分真實后驗概率的積分表達式難以簡化計算，估計工作無法順利進行。本文運用變分法來逼近后驗概率，其關(guān)鍵思想是選擇一系列易于操作的分布族來逼近真實的后驗概率。假設(shè)Ω=(Ω1,Ω2,…,ΩJ)是待估參數(shù)向量，對于一個給定的模型設(shè)計參數(shù)Φ1={Ω,Σ}，Φ2={Π,A}，完整數(shù)據(jù)的似然函數(shù)為：p(Y,S|Φ1,Φ2|)=p(Y|S,Φ1|)p(S|Φ2|)(1)其中，(2)(3)根據(jù)貝葉斯規(guī)則，邊際似然函數(shù)為：(4)此時，隱式變量和模型參數(shù)都被視為隨機變量，假設(shè)隱變量和模型參數(shù)的變分后驗分布為q(S,Φ1,Φ2)，結(jié)合式(4)可得對數(shù)邊際似然函數(shù)，(5)進一步整理式(5)，可得：其中，(6)為了計算變分的下界L(q)，需要指定變分后驗概率的形式和模型參數(shù)的先驗分布。本文中，變分后驗概率的因式分解形式如下：q(S,Φ1,Φ2)=q(S)q(Π)q(A)q(Ω,Σ)(7)為了簡便計算，需要先驗分布和后驗分布的形式保持高度一致。設(shè)Π和Α矩陣中每一行的先驗分布為狄利克雷(DirichletDistribution,Dir)分布，即(8)(9)文獻[12-13]認為，參數(shù)Ω的先驗分布為正態(tài)-逆Wishart分布，即(10)式中，(11)(12)式中，W-1表示逆Wishart分布。模型參數(shù)的先驗分布表示為：p(Φ1,Φ2)=p(Ω,Σ)p(Π)p(Α)(13)此時，變分下界L(q)可通過式(6)求解，將式(6)改寫為：L(q)=Eq[lnp(Y,S,Φ1,Φ2)]-Eq[lnq(S,Φ1,Φ2)](14)式中，E(·)表示期望，將式(1)—式(3)、式(7)—式(13)代入式(14)，可得：Eq[lnp(Y,S,Φ1,Φ2)]=Eq[lnp(Y|S,Φ1)+lnp(S|Φ2)+lnp(Φ1)+lnp(Φ2)]=(15)Eq[lnq(S,Φ1,Φ2)]=E[lnq(S)+lnq(A)+lnq(Π)+lnq(Ω,Σ)](16)式(15)中，δ(·)表示δ函數(shù)，當s(t-1)i與stj的乘積為1時，δ(·)為1，否則為0。在變分貝葉斯的框架下，變分后驗概率的近似最優(yōu)解為：lnqg(Zg)=Eh≠g[lnp(Y,Z)]+C(17)式中，Z={S,Φ1,Φ2}表示狀態(tài)變量和參數(shù)的集合，Zg表示集合Z中第g項參數(shù)，不涉及Zg參數(shù)的項均視為常數(shù)，記作C。求解參數(shù)Zg的變分后驗時，將除Zg外的參數(shù)視為隨機變量，參數(shù)Zg視為固定變量，對式(17)中與Zg參數(shù)相關(guān)、但不等于參數(shù)Zg的項進行計算，得到隨機變量的數(shù)學期望Eh≠g[lnp(Y,Z)]。3變分后驗分布運用式(17)對模型中的各個參數(shù)進行變分推斷，求解參數(shù)的變分后驗分布，并利用后驗分布計算變分下界。3.1q(Π),q(A)變分推斷求解變分后驗概率q(Π)時，將Π視為一個參數(shù)，其余的S,A,Ω,Σ視為隨機變量，因此與Π無關(guān)的項即為常數(shù)。將所有與Π相關(guān)的項合并在一起，計算結(jié)果如下：lnq*(Π)=ES,A[lnp(S|Π,A)]+lnp(Π)+C=(18)同理，可以計算出q(A)，lnq*(A)=ES,Π[lnp(S|Π,A)]+lnp(A)+C=(19)3.2q(Ω,Σ)變分推斷將所有與Ω,Σ相關(guān)的項整合在一起，得到：(20)3.3q(S)變分推斷將所有與隱變量S相關(guān)的項整合，得到：lnq*(S)=EΠ,A[lnp(S|Π,A)]+EΩ,Σ[lnp(Y|S,Ω,Σ)]+C=(21)對于隱馬爾可夫模型，隱變量q(S)變分后驗的主要形式為：(22)比對式(21)和式(22)，可得：(23)(24)(25)在得到q(S)的分布后，計算得出隱變量的單點均值ES(stj)和兩點均值ES(s(t-1)i,stj)，(26)(27)式中，fvar(·)為前向概率，bvar(·)為后向概率。3.4變分下界L(q)結(jié)合隱變量和參數(shù)的先驗分布，運用變分貝葉斯方法推導出隱變量和參數(shù)的變分后驗分布后，將式(18)—式(27)代入式(14)，計算得出的變分下界如下：L(q)=E[lnp(Π)]+E[lnp(A)]+EΩ,Σ[lnp(Ω,Σ)]+ES,Ω,Σ[lnp(Y|S,Ω,Σ)]+ES,Π,A[lnp(S|A,Π)]-E[lnq(S)]-E[lnq(A)]-E[lnq(Π)]-E[lnq(Ω,Σ)](28)式中，(29)(30)(31)(32)(33)(34)(35)(36)(37)將式(29)—式(37)代入式(28)，可得：根據(jù)狄利克雷分布的性質(zhì)，計算得出的相應(yīng)期望如下：4實證分析4.1數(shù)據(jù)選取與處理本文從我國國家外匯管理局門戶網(wǎng)站選擇人民幣兌美元匯率中間價歷史數(shù)據(jù)，從網(wǎng)易財經(jīng)官網(wǎng)選擇上海證券綜合指數(shù)(簡稱上證綜指或滬指)、深圳證券成份指數(shù)(簡稱深證成指)的歷史交易日度數(shù)據(jù)，組成一個3變量原始樣本，數(shù)據(jù)采集時間從2015-01-05日至2021-06-30，根據(jù)金融行業(yè)交易特點即有交易日與非交易日之分，剔除非交易日后共收集到1580組數(shù)據(jù)。人民幣兌美元匯率的收益率可由匯率中間價時間序列數(shù)據(jù)經(jīng)過計算而得，計算公式如下：Rt=100×(lnrt-lnrt-1)式中，rt為當日匯率的中間價，Rt為處理后的匯率收益率，自然對數(shù)用于消除系列的不穩(wěn)定性，為了更好地減少誤差，其百分比形式需要乘以100。用相同原理處理滬指和深證成指數(shù)據(jù)，得到相應(yīng)的收益率。本文用Rratio表示人民幣兌美元匯率中間價的日對數(shù)收益率、Rsz表示上證綜合指數(shù)的日對數(shù)收益率、Rsc表示深證成分指數(shù)的日對數(shù)收益率。對模型設(shè)計參數(shù)進行分析估計前，先繪制3組日對數(shù)收益率數(shù)據(jù)的時間序列圖，判斷其經(jīng)濟平穩(wěn)性。時間序列圖如圖1所示。圖13組日對數(shù)收益率的時間序列圖從圖1可以看出，3組對數(shù)收益率曲線均圍繞0值上下波動，除圖1(a)中出現(xiàn)1個異常點外，其余部分波動幅度前后、上下一致，故認為是平穩(wěn)序列。由于圖檢驗法具有較強的主觀性，為了確保檢驗結(jié)果的穩(wěn)健性，使用Eviews10.0軟件進行增廣迪基-富勒(AugmentedDickey-Fuller,ADF)檢驗和菲利普斯-珀森(Phillips-Perron,PP)檢驗。ADF檢驗和PP檢驗的零假設(shè)H0均為存在單位根，區(qū)別在于兩者的檢驗統(tǒng)計量不同，如果序列平穩(wěn)，則該序列不存在單位根。當顯著性水平為α時，檢驗統(tǒng)計量小于該水平下的迪基-富勒(Dickey-Fuller,DF)臨界值，此時應(yīng)該拒絕原假設(shè)，在α顯著性水平下認為該序列平穩(wěn)。3組日對數(shù)收益率的檢驗結(jié)果如表1所示。表1不同對數(shù)收益率的平穩(wěn)性檢驗表1中，變量Rratio的ADF值為-36.52076，顯著性水平1%DF臨界值為-3.434157，其ADF值小于顯著性水平1%DF臨界值，且Rratio相應(yīng)的PP統(tǒng)計量值為-36.79494，同樣小于1%DF臨界值，表明拒絕原假設(shè)，故認為變量Rratio為平穩(wěn)序列。同理可得，變量Rsz和Rsc拒絕原假設(shè)，因此，3組參數(shù)變量均為平穩(wěn)序列。4.2參數(shù)估計在建立模型之前，需要選擇合適的模型參數(shù)。在選擇時滯順序時，應(yīng)考慮模型的整體動態(tài)性和自由度。階數(shù)較大時，可以更好地反映一個模型研究動態(tài)性，但自由度會減少，需要估計的參數(shù)會增加。通常情況下，滯后階數(shù)的確定根據(jù)赤池信息準則(AkaikeInformationCriterion,AIC)、施瓦茨準則(SchwarzCriterion,SC)、漢南-奎因準則(Hannan-QuinnCriterion,HQ)和最終預測誤差準則(FinalPredictionError,FPE)來進行綜合判斷，準則的值越小，表明模型擬合效果越好。為了選取合適的滯后階數(shù)，運用Eviews10.0軟件對3組對數(shù)收益率數(shù)據(jù)做VAR分析，結(jié)果如表2所示。表2不同準則下的滯后階數(shù)從表2可以看出，在滯后階數(shù)為1時，其FPE，AIC與HQ準則的值均達到最小，表明該準則下的最優(yōu)滯后期為1。本文經(jīng)過綜合分析，選取MS-VAR模型的滯后階數(shù)為1。由于本文研究的參數(shù)是人民幣兌美元匯率及