呂秀麗，虞棟杰，鄭 靜

(杭州電子科技大學(xué)經(jīng)濟(jì)學(xué)院，浙江 杭州 310018)

0 引 言

研究變量間相互依賴關(guān)系的時(shí)間序列是較為經(jīng)典的統(tǒng)計(jì)分析方法，其中向量自回歸模型(Vector Autoregressive Models, VAR)是常用的計(jì)量經(jīng)濟(jì)模型。使用VAR模型對時(shí)間序列進(jìn)行分析時(shí)，需先假設(shè)時(shí)間序列線性同質(zhì)，但在實(shí)際應(yīng)用中，各種復(fù)雜因素致使線性假定成立模型的分析結(jié)果不夠準(zhǔn)確。為了解決這個(gè)問題，學(xué)者們將研究方向轉(zhuǎn)向非線性模型，龍如銀等[1]運(yùn)用區(qū)制為3、滯后階數(shù)為1的Markov區(qū)制轉(zhuǎn)換模型(Markov Switching Models, MS)研究我國1984年至2003年的通貨膨脹率，與傳統(tǒng)VAR模型相比，更好地?cái)M合了我國經(jīng)濟(jì)的通貨膨脹率。劉金全等[2]運(yùn)用Markov轉(zhuǎn)換向量自回歸(Markov Switching-Vector Autoregressive Models, MS-VAR)模型將我國通貨膨脹率分為3個(gè)區(qū)制，分析經(jīng)濟(jì)政策與通貨膨脹率在各區(qū)制下的關(guān)系，合理解釋了當(dāng)時(shí)我國社會(huì)經(jīng)濟(jì)及通貨膨脹風(fēng)險(xiǎn)狀況，指出治理國家通貨膨脹應(yīng)先判斷其所處區(qū)制。朱慧明等[3]引入貝葉斯MS-VAR模型，2次調(diào)用Gibbs抽樣方法，研究我國股市與通貨膨脹率的波動(dòng)關(guān)系，指出可通過機(jī)制轉(zhuǎn)換模型來刻畫該波動(dòng)關(guān)系。以上模型的參數(shù)較多，需事先進(jìn)行參數(shù)估計(jì)。目前常用的模型參數(shù)估計(jì)方法有2類，一是確定性方法，如極大似然估計(jì)、EM算法等，在估計(jì)過程中未考慮模型參數(shù)的先驗(yàn)信息；二是貝葉斯方法、蒙特卡羅法等，計(jì)算成本較高，易產(chǎn)生過擬合或欠擬合。變分貝葉斯(Variational Bayes,VB)是貝葉斯推理的近似方法。Adami[4]將變分方法應(yīng)用到貝葉斯推理中，介紹了如何用變分方法來逼近這類推理中出現(xiàn)的難以處理的后驗(yàn)分布；Alan等[5]采用參數(shù)擴(kuò)展變分貝葉斯方法解決了VB收斂到近似解較慢的問題；Ormerod等[6]將變分近似引入統(tǒng)計(jì)學(xué)，用統(tǒng)計(jì)學(xué)術(shù)語解釋變分近似概念；Koop等[7]提出一種均值場變分貝葉斯算法，分析美國宏觀經(jīng)濟(jì)數(shù)據(jù)，將其拓展至計(jì)量經(jīng)濟(jì)學(xué)領(lǐng)域；馬占宇等[8]結(jié)合變分推斷，提出一種學(xué)術(shù)研究熱點(diǎn)關(guān)鍵詞提取方法，高效、準(zhǔn)確地提取學(xué)術(shù)研究的熱點(diǎn)關(guān)鍵詞；沈圓圓等[9]提出一種Logistic組稀疏性回歸的Bayes概率模型的變分貝葉斯算法，提高了模型的預(yù)測性能；萬志成等[10]建立了一種無限變量高斯混合模型，解決了無限變量高斯混合模型的參數(shù)估計(jì)問題；付維明等[11]運(yùn)用擴(kuò)散方法對分布式隨機(jī)變分推斷算法進(jìn)行改進(jìn)，并將其應(yīng)用至主題模型，速度快，可擴(kuò)展強(qiáng)。本文提出一種Markov轉(zhuǎn)換向量自回歸模型的變分貝葉斯推斷算法，利用參數(shù)的先驗(yàn)信息對MS-VAR模型的參數(shù)進(jìn)行估計(jì)，并給出完整的變分推理過程。

1 MS-VAR模型

1.1 馬爾可夫轉(zhuǎn)移模型

給定所有先前的狀態(tài)，馬爾可夫鏈(本文簡稱馬氏鏈)的當(dāng)前狀態(tài)依賴且僅依賴其前一個(gè)狀態(tài)，這一性質(zhì)被稱為馬爾可夫特性。馬爾可夫區(qū)制轉(zhuǎn)移模型依賴具有馬爾可夫性的馬氏鏈生成的隨機(jī)過程。

設(shè){yt,t=1,…,T}為具有T個(gè)單變量觀測值的時(shí)間序列，是隨機(jī)過程{yt}的子集。隱藏的離散隨機(jī)過程{zt}的實(shí)現(xiàn)決定了隨機(jī)過程{yt}的概率分布。{yt}是直接可觀測的，而{zt}是一個(gè)潛在的隨機(jī)變量，只能通過觀測{yt}來間接觀測過程{zt}。假設(shè)隱過程{zt}是一個(gè)具有有限狀態(tài)空間{0,…,K-1}的不可約、非周期的馬氏鏈，其隨機(jī)性質(zhì)用K×K階轉(zhuǎn)移矩陣A描述，其中aij是轉(zhuǎn)移矩陣A的第i+1行第j+1列元素，表示在t-1時(shí)刻處于狀態(tài)i，下一時(shí)刻轉(zhuǎn)移到狀態(tài)j的概率。轉(zhuǎn)移矩陣A中每一行的所有元素之和為1且所有元素非負(fù)，即：

設(shè)參數(shù)θi∈Θ是由轉(zhuǎn)移區(qū)制決定的參數(shù)，而ψ∈Ψ是不依賴于區(qū)制的參數(shù)，Θ,Ψ為參數(shù)空間，θi和ψ共同決定了參數(shù)條件密度p(·|θi,ψ|)。本文假設(shè)在給定z1,…,zΤ時(shí)的隨機(jī)變量y1,…,yΤ是條件獨(dú)立的，其密度為：

yt|(zt=i)|～p(·|θi,ψ|)

對于聯(lián)合隨機(jī)過程{zt,yt}，yt的條件密度為：

其中，It-1={yt-1,yt-2,…}為t-1時(shí)刻可獲得的觀測信息。

1.2 馬爾可夫區(qū)制轉(zhuǎn)移的向量自回歸模型

假設(shè)yt=(y1t,y2t,…,ymt)′為m×1的觀察值向量，隨機(jī)變量st=1,2,…,J(J≥2)表示在t時(shí)刻的狀態(tài)，與誤差向量εst獨(dú)立，誤差向量εst服從均值為0且協(xié)方差陣為Σst的正態(tài)分布，存在一個(gè)初始狀態(tài)概率向量為Π、狀態(tài)轉(zhuǎn)移概率矩陣為A的馬氏鏈，St={s1,s2,…,st}是由該馬氏鏈生成的不可觀測的隨機(jī)狀態(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

進(jìn)一步轉(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í)刻所有參數(shù)的信息。

2 變分貝葉斯推斷

貝葉斯估計(jì)方法通過設(shè)置參數(shù)的先驗(yàn)分布，根據(jù)貝葉斯規(guī)則，利用樣本信息推導(dǎo)出參數(shù)后驗(yàn)分布，得到參數(shù)估計(jì)值。但是，由于大部分真實(shí)后驗(yàn)概率的積分表達(dá)式難以簡化計(jì)算，估計(jì)工作無法順利進(jìn)行。本文運(yùn)用變分法來逼近后驗(yàn)概率，其關(guān)鍵思想是選擇一系列易于操作的分布族來逼近真實(shí)的后驗(yàn)概率。

假設(shè)Ω=(Ω1,Ω2,…,ΩJ)是待估參數(shù)向量，對于一個(gè)給定的模型設(shè)計(jì)參數(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ù)都被視為隨機(jī)變量，假設(shè)隱變量和模型參數(shù)的變分后驗(yàn)分布為q(S,Φ1,Φ2)，結(jié)合式(4)可得對數(shù)邊際似然函數(shù)，

(5)

進(jìn)一步整理式(5)，可得：

其中，

(6)

為了計(jì)算變分的下界L(q)，需要指定變分后驗(yàn)概率的形式和模型參數(shù)的先驗(yàn)分布。本文中，變分后驗(yàn)概率的因式分解形式如下：

q(S,Φ1,Φ2)=q(S)q(Π)q(A)q(Ω,Σ)

(7)

為了簡便計(jì)算，需要先驗(yàn)分布和后驗(yàn)分布的形式保持高度一致。設(shè)Π和Α矩陣中每一行的先驗(yàn)分布為狄利克雷(Dirichlet Distribution, Dir)分布，即

(8)

(9)

文獻(xiàn)[12-13]認(rèn)為，參數(shù)Ω的先驗(yàn)分布為正態(tài)-逆Wishart分布，即

(10)

式中，

(11)

(12)

式中，W-1表示逆Wishart分布。

模型參數(shù)的先驗(yàn)分布表示為：

p(Φ1,Φ2)=p(Ω,Σ)p(Π)p(Α)

(13)

此時(shí)，變分下界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ù)，當(dāng)s(t-1)i與stj的乘積為1時(shí)，δ(·)為1，否則為0。

在變分貝葉斯的框架下，變分后驗(yàn)概率的近似最優(yōu)解為：

lnqg(Zg)=Eh≠g[lnp(Y,Z)]+C

(17)

式中，Z={S,Φ1,Φ2}表示狀態(tài)變量和參數(shù)的集合，Zg表示集合Z中第g項(xiàng)參數(shù)，不涉及Zg參數(shù)的項(xiàng)均視為常數(shù)，記作C。求解參數(shù)Zg的變分后驗(yàn)時(shí)，將除Zg外的參數(shù)視為隨機(jī)變量，參數(shù)Zg視為固定變量，對式(17)中與Zg參數(shù)相關(guān)、但不等于參數(shù)Zg的項(xiàng)進(jìn)行計(jì)算，得到隨機(jī)變量的數(shù)學(xué)期望Eh≠g[lnp(Y,Z)]。

3 變分后驗(yàn)分布

運(yùn)用式(17)對模型中的各個(gè)參數(shù)進(jìn)行變分推斷，求解參數(shù)的變分后驗(yàn)分布，并利用后驗(yàn)分布計(jì)算變分下界。

3.1 q(Π),q(A)變分推斷

求解變分后驗(yàn)概率q(Π)時(shí)，將Π視為一個(gè)參數(shù)，其余的S,A,Ω,Σ視為隨機(jī)變量，因此與Π無關(guān)的項(xiàng)即為常數(shù)。將所有與Π相關(guān)的項(xiàng)合并在一起，計(jì)算結(jié)果如下：

lnq*(Π)=ES,A[lnp(S|Π,A)]+lnp(Π)+C=

(18)

同理，可以計(jì)算出q(A)，

lnq*(A)=ES,Π[lnp(S|Π,A)]+lnp(A)+C=

(19)

3.2 q(Ω,Σ)變分推斷

將所有與Ω,Σ相關(guān)的項(xiàng)整合在一起，得到：

(20)

3.3 q(S)變分推斷

將所有與隱變量S相關(guān)的項(xiàng)整合，得到：

lnq*(S)=EΠ,A[lnp(S|Π,A)]+EΩ,Σ[lnp(Y|S,Ω,Σ)]+C=

(21)

對于隱馬爾可夫模型，隱變量q(S)變分后驗(yàn)的主要形式為：

(22)

比對式(21)和式(22)，可得：

(23)

(24)

(25)

在得到q(S)的分布后，計(jì)算得出隱變量的單點(diǎn)均值ES(stj)和兩點(diǎn)均值ES(s(t-1)i,stj)，

(26)

(27)

式中，fvar(·)為前向概率，bvar(·)為后向概率。

3.4 變分下界L(q)

結(jié)合隱變量和參數(shù)的先驗(yàn)分布，運(yùn)用變分貝葉斯方法推導(dǎo)出隱變量和參數(shù)的變分后驗(yàn)分布后，將式(18)—式(27)代入式(14)，計(jì)算得出的變分下界如下：

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ì)，計(jì)算得出的相應(yīng)期望如下：

4 實(shí)證分析

4.1 數(shù)據(jù)選取與處理

本文從我國國家外匯管理局門戶網(wǎng)站選擇人民幣兌美元匯率中間價(jià)歷史數(shù)據(jù)，從網(wǎng)易財(cái)經(jīng)官網(wǎng)選擇上海證券綜合指數(shù)(簡稱上證綜指或滬指)、深圳證券成份指數(shù)(簡稱深證成指)的歷史交易日度數(shù)據(jù)，組成一個(gè)3變量原始樣本，數(shù)據(jù)采集時(shí)間從2015-01-05日至2021-06-30，根據(jù)金融行業(yè)交易特點(diǎn)即有交易日與非交易日之分，剔除非交易日后共收集到1 580組數(shù)據(jù)。人民幣兌美元匯率的收益率可由匯率中間價(jià)時(shí)間序列數(shù)據(jù)經(jīng)過計(jì)算而得，計(jì)算公式如下：

Rt=100×(lnrt-lnrt-1)

式中，rt為當(dāng)日匯率的中間價(jià)，Rt為處理后的匯率收益率，自然對數(shù)用于消除系列的不穩(wěn)定性，為了更好地減少誤差，其百分比形式需要乘以100。

用相同原理處理滬指和深證成指數(shù)據(jù)，得到相應(yīng)的收益率。本文用Rratio表示人民幣兌美元匯率中間價(jià)的日對數(shù)收益率、Rsz表示上證綜合指數(shù)的日對數(shù)收益率、Rsc表示深證成分指數(shù)的日對數(shù)收益率。

對模型設(shè)計(jì)參數(shù)進(jìn)行分析估計(jì)前，先繪制3組日對數(shù)收益率數(shù)據(jù)的時(shí)間序列圖，判斷其經(jīng)濟(jì)平穩(wěn)性。時(shí)間序列圖如圖1所示。

圖1 3組日對數(shù)收益率的時(shí)間序列圖

從圖1可以看出，3組對數(shù)收益率曲線均圍繞0值上下波動(dòng)，除圖1(a)中出現(xiàn)1個(gè)異常點(diǎn)外，其余部分波動(dòng)幅度前后、上下一致，故認(rèn)為是平穩(wěn)序列。由于圖檢驗(yàn)法具有較強(qiáng)的主觀性，為了確保檢驗(yàn)結(jié)果的穩(wěn)健性，使用Eviews10.0軟件進(jìn)行增廣迪基-富勒(Augmented Dickey-Fuller, ADF)檢驗(yàn)和菲利普斯-珀森(Phillips-Perron, PP)檢驗(yàn)。ADF檢驗(yàn)和PP檢驗(yàn)的零假設(shè)H0均為存在單位根，區(qū)別在于兩者的檢驗(yàn)統(tǒng)計(jì)量不同，如果序列平穩(wěn)，則該序列不存在單位根。當(dāng)顯著性水平為α?xí)r，檢驗(yàn)統(tǒng)計(jì)量小于該水平下的迪基-富勒(Dickey-Fuller, DF)臨界值，此時(shí)應(yīng)該拒絕原假設(shè)，在α顯著性水平下認(rèn)為該序列平穩(wěn)。3組日對數(shù)收益率的檢驗(yàn)結(jié)果如表1所示。

表1 不同對數(shù)收益率的平穩(wěn)性檢驗(yàn)

表1中，變量Rratio的ADF值為-36.520 76，顯著性水平1%DF臨界值為-3.434 157，其ADF值小于顯著性水平1%DF臨界值，且Rratio相應(yīng)的PP統(tǒng)計(jì)量值為-36.794 94，同樣小于1%DF臨界值，表明拒絕原假設(shè)，故認(rèn)為變量Rratio為平穩(wěn)序列。同理可得，變量Rsz和Rsc拒絕原假設(shè)，因此，3組參數(shù)變量均為平穩(wěn)序列。

4.2 參數(shù)估計(jì)

在建立模型之前，需要選擇合適的模型參數(shù)。在選擇時(shí)滯順序時(shí)，應(yīng)考慮模型的整體動(dòng)態(tài)性和自由度。階數(shù)較大時(shí)，可以更好地反映一個(gè)模型研究動(dòng)態(tài)性，但自由度會(huì)減少，需要估計(jì)的參數(shù)會(huì)增加。通常情況下，滯后階數(shù)的確定根據(jù)赤池信息準(zhǔn)則(Akaike Information Criterion , AIC)、施瓦茨準(zhǔn)則(Schwarz Criterion, SC)、漢南-奎因準(zhǔn)則(Hannan-Quinn Criterion, HQ)和最終預(yù)測誤差準(zhǔn)則(Final Prediction Error, FPE)來進(jìn)行綜合判斷，準(zhǔn)則的值越小，表明模型擬合效果越好。為了選取合適的滯后階數(shù)，運(yùn)用Eviews10.0軟件對3組對數(shù)收益率數(shù)據(jù)做VAR分析，結(jié)果如表2所示。

表2 不同準(zhǔn)則下的滯后階數(shù)

從表2可以看出，在滯后階數(shù)為1時(shí)，其FPE，AIC與HQ準(zhǔn)則的值均達(dá)到最小，表明該準(zhǔn)則下的最優(yōu)滯后期為1。本文經(jīng)過綜合分析，選取MS-VAR模型的滯后階數(shù)為1。

由于本文研究的參數(shù)是人民幣兌美元匯率及滬深指數(shù)，分析時(shí)結(jié)合其經(jīng)濟(jì)發(fā)展意義，將馬爾可夫的區(qū)制即狀態(tài)變量設(shè)為2，即劃分為匯率上升或匯率下降這2種不同體制。假設(shè)d=1,j=2，狀態(tài)變量為st，MS-VAR模型表示為：

(38)

由式(38)可以看出，模型中有28個(gè)待估計(jì)參數(shù)。進(jìn)行模型估計(jì)時(shí)，若待估參數(shù)過多，一方面會(huì)增加估計(jì)成本，另一方面預(yù)測效果并不一定能達(dá)到預(yù)期，甚至可能偏離經(jīng)濟(jì)理論。本文提出的Markov轉(zhuǎn)換向量自回歸模型的變分貝葉斯推斷算法使得參數(shù)估計(jì)變得更加精確，且收斂速度較快，節(jié)約了時(shí)間成本。

取T=1 580，運(yùn)行MATLAB R2018a軟件，使用本文提出的Markov轉(zhuǎn)換向量自回歸模型的變分貝葉斯推斷算法迭代1 000次，并針對迭代結(jié)果，計(jì)算樣本均值的標(biāo)準(zhǔn)差，即標(biāo)準(zhǔn)誤，得到模型轉(zhuǎn)移概率的估計(jì)結(jié)果如表3所示。

表3 轉(zhuǎn)移概率的估計(jì)結(jié)果

迭代后得到模型參數(shù)矩陣的估計(jì)結(jié)果是：

5 結(jié)束語

針對Markov轉(zhuǎn)換向量自回歸模型應(yīng)用傳統(tǒng)參數(shù)估計(jì)方法成本高、易過擬合等缺陷，本文在變分貝葉斯的理論框架下，提出一種Markov轉(zhuǎn)換向量自回歸模型的變分貝葉斯推斷算法。詳細(xì)演算MS-VAR模型的變分推導(dǎo)過程，給出參數(shù)估計(jì)算法，并結(jié)合實(shí)證分析驗(yàn)證了本文算法在高維數(shù)據(jù)結(jié)構(gòu)情形下的高收斂性。但是，本文模型中存在多個(gè)系統(tǒng)，難以將觀測值劃分為各狀態(tài)參數(shù)，無法估計(jì)所有系統(tǒng)的參數(shù)，下一步計(jì)劃將支持向量機(jī)、邏輯回歸等機(jī)器學(xué)習(xí)分類方法與MS-VAR模型變分貝葉斯方法結(jié)合起來，研究模型中不同區(qū)制的區(qū)分及其參數(shù)估計(jì)。