GJR-CAViaR模型的貝葉斯分位數(shù)回歸
——基于Gibbs抽樣的MCMC算法實現(xiàn)

張 穎 傅 強

一、引言

目前,常用的分位數(shù)回歸模型的估計方法分為兩類。一類是直接進行優(yōu)化求解,如單純形法和內點法。另一類是借助于貝葉斯原理進行參數(shù)估計。

直接優(yōu)化求解屬于頻率學派的范疇,是傳統(tǒng)的經典統(tǒng)計學方法。經典估計方法將參數(shù)視為固定常數(shù),然后利用最小二乘或極大似然等方法計算參數(shù)的估計值,得到參數(shù)的漸近分布和統(tǒng)計性質,并進行假設檢驗。貝葉斯學派與經典統(tǒng)計法在參數(shù)估計的原理上存在不同。貝葉斯學派將待估參數(shù)視為隨機變量,利用貝葉斯原理和觀測樣本得到參數(shù)的后驗分布。在無法得到參數(shù)后驗分布的具體表達形式時,采用重復抽樣技術解決參數(shù)的估計問題。因此,相對于傳統(tǒng)統(tǒng)計對樣本量的敏感,貝葉斯統(tǒng)計在小樣本情形下也能得到可靠的參數(shù)信息。與經典統(tǒng)計方法的另一區(qū)別是貝葉斯方法在估計中加入了參數(shù)的先驗信息。合適的先驗分布可提高估計的精度。而當無先驗信息可提供時,參數(shù)后驗分布的估計會和最小二乘估計量一致,并不會影響到估計量的良好的統(tǒng)計性質。

總之,貝葉斯估計可以充分利用樣本信息以及參數(shù)的先驗信息,是簡單有效地獲得參數(shù)后驗分布的方法。從損失函數(shù)最小化出發(fā),Koenker和 Machado(1999)［1］將分位數(shù)回歸與貝葉斯統(tǒng)計推斷相結合,首次建立了服從非對稱拉普拉斯分布ALD(Asymmetric Laplace Distribution)的誤差項和目標參數(shù)的似然函數(shù)之間的關系。Yu和Moyeed(2001)［2］完善了這一方法,并詳細證明了求損失函數(shù)最小等價于誤差項服從ALD時最大化似然函數(shù)。

在分位數(shù)回歸的CAViaR模型中,貝葉斯方法的應用也越來越普遍?；诓粚ΨQ拉普拉斯分布,Hsu(2010)［3］討論了CAViaR模型中參數(shù)的先驗選擇,并利用MCMC方法完成了對模型的貝葉斯分析。曾惠芳(2011)［4］利用三階傅里葉級數(shù)來擬合CAViaR模型中的非參部分,研究了上證綜指周數(shù)據(jù)的在險價值VaR?；贑hou(2005)［5］提出的日內極差(CARR)模型,Chen等(2012)［6］建立了極差以及閥值范圍值(TRV,Threshold Range Value)的CAViaR模型。他們采用亞太經濟合作體的五國的市場指數(shù)數(shù)據(jù)(S＆P500,Nikkei225,TAIEX,HIS,KOSPI) 以及兩國匯率數(shù)據(jù)(歐元對美元和日元對美元),發(fā)現(xiàn)貝葉斯CAViaR模型比傳統(tǒng)的CAViaR能更準確地預測市場風險VaR,且優(yōu)于歷史模擬法及參數(shù)模型(GARCH和RiskMetrics) 方法。陳磊等(2013)［7］對布倫特原油價格的實證表明,貝葉斯CAViaR模型比傳統(tǒng)的CAViaR模型對油價風險VaR的預測更好,根據(jù)平均相對偏差和市場風險資本要求發(fā)現(xiàn),不對稱CAViaR模型可以有效刻畫油價VaR的動態(tài)變化。王新宇等(2013)［8］構建了有變點的CAViaR模型,提出變點的SAV和IGARCH模型的貝葉斯推理及MCMC實現(xiàn),并檢驗了創(chuàng)業(yè)板指數(shù)市場風險因首次限售解禁而發(fā)生的結構改變。這些研究表明,與傳統(tǒng)的優(yōu)化求解相比,貝葉斯方法能更好地估計CAViaR模型的參數(shù)。

張穎等(2012)［9］提出GJR-CAViaR模型研究收益對風險的不對稱影響,彌補了AS-CAViaR和TARCH-CAViaR模型設定中在參數(shù)限制上的不足,并利用實證指出在GJR、AS及TARCH三種CAViaR模型中,GJR-CAViaR模型的樣本內預測結果總體效果較好。但是,尚未有文獻提出如何進行GJR-CAViaR模型參數(shù)的貝葉斯推斷。GJR-CAViaR模型屬于遞歸的分位數(shù)回歸方程,模型的非線性及半?yún)⑿栽斐韶惾~斯估計比較復雜,而且,不對稱拉普拉斯分布屬于非標準的分布,只有將模型參數(shù)的滿條件分布轉化為常用已知分布,才能借助于Gibbs抽樣來估計參數(shù)。因此,有必要對如何解決這類模型的貝葉斯估計提出分析框架。本文提出GJR-CAViaR模型的貝葉斯估計的算法實現(xiàn),并基于Gibbs抽樣建立上證綜指日收益率的GJR-CAViaR模型參數(shù)估計。

二、分位數(shù)回歸的貝葉斯估計理論基礎

（一）分位數(shù)回歸理論

Koenker和 Bassett(1978)［10］提出了分位數(shù)回歸理論,該理論是對中位數(shù)回歸的擴展。假設線性回歸方程為:

定義誤差項u的損失函數(shù)ρ(u)如式(2)所示。其中,I(u)為指示函數(shù),0＜τ＜1。

分位數(shù)回歸估計的目標是損失函數(shù)最小,其目標函數(shù)表示如式(4)。

可以看出,分位數(shù)回歸對正負誤差項賦予了不同的權重。當τ＝0.5時,式(4)的目標函數(shù)將簡化為誤差項的絕對值之和最小(least absolute deviation,LAD),也就是傳統(tǒng)上的中位數(shù)回歸。

目標函數(shù)式(4)不可微,因此無法得到模型的解析解。這個時候只能通過線性規(guī)劃或者貝葉斯的方法來求解。

（二）貝葉斯分位數(shù)回歸

Yu和Moyeed(2001)［2］的研究是該領域最重要的文獻之一。他們給出貝葉斯分位數(shù)回歸的原理,并基于不對稱拉普拉斯分布建立分位數(shù)回歸模型的貝葉斯估計框架。

拉普拉斯分布以皮埃爾—西蒙·拉普拉斯的名字命名。標準的對稱的拉普拉斯分布是由兩個具有同樣的位置參數(shù)的獨立同分布的指數(shù)分布背靠背拼接在一起,也被稱為雙指數(shù)分布。一般化的拉普拉斯分布考慮了分布的非對稱性,Yu和Moyeed(2001)［2］給出服從不對稱拉普拉斯分布的隨機變量x的概率密度函數(shù)表示,如式(5)。

式(5)中,σ＞0是尺度參數(shù)(scale parameter),μ是位置參數(shù)(location parameter),0＜p＜1是偏度參數(shù)(asymmetrc parameter),ρp(x－μ)是損失函數(shù),定義如式(2)。x滿足不對稱拉普拉斯分布,記作x～ALD(μ,σ,p) 。

隨機變量x的分位數(shù)函數(shù)為

可以看出,

式(7)說明隨機變量x的p分位數(shù)等于其位置參數(shù)μ,這個就是ALD分布可作為分位數(shù)回歸模型的誤差分布的依據(jù)。而偏度參數(shù)p就是式(2)中的分位數(shù)水平τ。

對于線性回歸模型yi＝x＇iβ＋ui,假定模型的殘差項u服從不對稱拉普拉斯分布,ui～ALD(0,σ,τ)。根據(jù)拉普拉斯分布的線性變換性質,yi～ALD(x＇iβ,σ,τ) 。則因變量y的概率密度函數(shù)可表示成式(8)。

則樣本y＝(y1,y2,…,yn)的似然函數(shù)可表示成式(9)。

其中,

可以看出,求分位數(shù)回歸的損失函數(shù)最小化等價于式(9)的似然函數(shù)的最大化。貝葉斯分位數(shù)回歸的關鍵就是利用誤差項的非對稱的拉普拉斯分布,將對損失函數(shù)的求解轉化為對似然函數(shù)的求解。

Yu和Moyeed(2001)［2］還證明了即使待估參數(shù)β的先驗分布是不真實的分布(Improper distribution),只要β滿足p(β)∝1,那么β的后驗分布就是真實分布(Proper distribution)。也就是說,β的后驗密度函數(shù)p(β｜y)存在,y是觀測值。即:

MCMC模擬所必須遵循的原則就是參數(shù)的后驗分布是真實分布,但對參數(shù)的先驗分布并沒有這個要求。所以,在貝葉斯分位數(shù)回歸中,即使參數(shù)的先驗分布設定為均勻分布,得到的后驗分布也是真實的,不會對MCMC的估計結果造成影響。

三、GJR-CAViaR模型的貝葉斯估計理論

張穎等(2012)［9］提出的GJR-CAViaR模型是基于CAViaR模型的不對稱收益模型,就正負收益對市場風險VaR的不對稱性建模,并定量描述此影響的不對稱程度。模型的設定如下:

當利用如式(15)的目標函數(shù)估計參數(shù)β的時候,也就是執(zhí)行分位數(shù)回歸的時候,gt(β)就是yt的τ條件分位數(shù)(τ是分位數(shù)水平),即在險價值VaR。但是,我們并不知道gt(β)的具體形式,只知道gt(β)的演化方程如式(13)。

下面討論GJR-CAViaR模型的Gibbs抽樣。Gibbs抽樣是將多變量抽樣轉化為對滿條件分布進行單變量抽樣的算法,重點是如何推斷出條件分布。以不對稱拉普拉斯分布為基礎的半?yún)⒌姆治粩?shù)回歸模型CAViaR模型的貝葉斯推斷很難得到參數(shù)的解析后驗分布,也就是說,無法得到參數(shù)解析的滿條件分布,不能直接使用Gibbs抽樣來實現(xiàn)對模型的模擬。因此,需要引入其他標準的分布來得到模型參數(shù)解析的滿條件分布,并實現(xiàn)模型的Gibbs抽樣。文獻中主要通過兩種方式解決這個問題。其一,利用混合正態(tài)分布與有偏的拉普拉斯分布(skewed-laplace)的等價關系(Tsionas,2003［11］),推導出不同參數(shù)的后驗邊緣分布,解決拉普拉斯分布對應的貝葉斯分位數(shù)回歸問題。其二,把不對稱拉普拉斯分布表示成標準指數(shù)分布和標準正態(tài)分布的組合,實現(xiàn)貝葉斯分位數(shù)回歸(Kozumi和 Kobayashi,2011［12］)。不對稱拉普拉斯分布和常用的正態(tài)分布及指數(shù)分布存在如下關系。假設隨機變量ε服從于不對稱拉普拉斯分布,z服從標準指數(shù)分布z～E(1),u是標準正態(tài)分布,u～N(0,1)。則有式(16)成立。

其中,

式(17)與(18)中的τ就是式(15)分位數(shù)回歸中的分位數(shù)水平τ。Kozumi和 Kobayashi(2011)［12］指出 Tsionas(2003)［11］的方法有兩點不足。其一,算法會產生高度相關的抽樣點從而降低有效性；其二,當樣本量很大的時候運算速度慢。故本文選擇 Kozumi和 Kobayashi(2011)［12］的方法來處理不對稱拉普拉斯分布的問題。根據(jù)式(12),假設隨機變量ε服從不對稱拉普拉斯分布,則GJR-CAViaR模型可等價表示為

當引入尺度參數(shù)σ時,上式可改寫為

令vt＝σ zt,則vt～E(σ),上式重新寫作

在給定vt的時候,yt服從均值為gtβ()＋θ vt,方差為φ2σ vt的條件正態(tài)分布。此時,可求得y＝(y1,…,yn)的聯(lián)合密度函數(shù)為

假定參數(shù)β的先驗分布為正態(tài)分布,σ的先驗分布為逆伽馬分布,σ～IG(n0/2,s0/2),其中,n0,s0為固定常數(shù),則β的滿條件分布也是正態(tài)分布,并可推導出下面的結論:vt的滿條件分布為逆高斯分布。

其中,

同時,σ的滿條件分布為逆伽馬分布。

其中,

從上述推導結果可以看到,模型各個參數(shù)的滿條件分布都服從常用的已知分布。這時,就可以利用Gibbs方法模擬并構造馬爾科夫鏈,具體步驟如下:

第一步,選擇參數(shù)的初始值,β0＝

第二步,從條件正態(tài)密度函數(shù)f(β0｜y,σ,v)中抽取隨機數(shù),記作β1。

第三步,從條件逆伽馬密度函數(shù)f(σ｜y,β1,v)中抽取隨機數(shù),記作σ1。

第四步,從條件逆高斯密度函數(shù)f(v｜y,β1,σ1)中抽取隨機數(shù),記作v1。

由此成功完成一次狀態(tài)轉移,參數(shù)變?yōu)?β1,σ1,v1)。利用這個新值,重復步驟1到步驟4,完成新的參數(shù)的迭代。如此重復N次,直到得到穩(wěn)定的馬爾科夫鏈,就可以得到參數(shù)的后驗分布及估計結果。

四、GJR-CAViaR模型的參數(shù)估計結果

王新宇等(2010)［13］以及邸俊鵬(2013)［14］都利用模擬數(shù)據(jù)證明了在分位數(shù)回歸貝葉斯推理中,估計量依賴于尺度參數(shù)的設置。將尺度參數(shù)σ參數(shù)化,可以降低參數(shù)估計值的標準差,提高被估系數(shù)的估計精度。這里,我們不再構造模型進行模擬分析,直接將尺寸參數(shù)σ看作待估變量。在實證部分,利用同花順軟件下載上證綜指的日收盤價數(shù)據(jù),計算日收盤指數(shù)的對數(shù)收益率。

自1996年12月16日開始,中國股市開始實行漲停板制度。樣本區(qū)間選擇為:1996年12月16日至2014年12月31日。因為方差很大的正態(tài)分布可近似為均勻分布,設定參數(shù)的先驗分布為方差大的正態(tài)分布,β0～N(0,100),尺度參數(shù)σ的先驗分布服從逆伽馬分布,σ～IG(10－5,10－5)。Gibbs抽樣共模擬20 000次,去掉前面10 000個抽樣值,利用后面的10 000個點計算參數(shù)估計值。本文涉及的抽樣方法通過R軟件和WinBUGS軟件實現(xiàn)。

（一）不同置信水平下GJR-CAViaR模型的參數(shù)估計

在樣本區(qū)間內,選擇最近的300個樣本①對于貝葉斯推斷,估計結果對樣本量的大小不太敏感。在邸俊鵬(2013)［14］的論文中,當樣本量大于200時Gibbs抽樣估計的參數(shù)的偏誤和標準差都比較小。因此,基于充足樣本量的考慮,本文選取最近的300個上證綜指日收益率的數(shù)據(jù)。,即:2013年7月25日至2014年12月31日的收益率數(shù)據(jù),計算99%置信水平下GJR-CAViaR模型的參數(shù)估計值。表1列出了估計參數(shù)的后驗均值、標準差、蒙特卡洛誤差以及95%的貝葉斯置信區(qū)間。

對Gibbs抽樣算法的結果進行評價有三種方式:接受率、收斂性以及抽樣的精度。我們的實證以收斂性作為判斷的標準。若馬爾科夫鏈收斂,則我們認為Gibbs模擬值是真實后驗分布的合理近似。簡單地,可以從抽樣值軌跡直觀地對抽樣分布進行檢驗。如果抽樣值估計在某值附近上下波動,則該馬爾科夫鏈收斂。另外,抽樣值的相關性太強將導致算法無法到達參數(shù)的整個空間。若隨著滯后階數(shù)的增加,自相關系數(shù)趨于零,該馬爾科夫鏈收斂。

表1 99%置信水平下模型參數(shù)的Gibbs抽樣估計值

圖1至圖5是樣本容量為300時,不同參數(shù)在0.01分位點下的抽樣軌跡和核密度圖。圖中顯示不同參數(shù)的抽樣序列都比較集中,在小范圍內上下波動,因此抽樣值構成的馬爾科夫鏈均收斂。

圖1 β0的MCMC抽樣值軌跡和核密度圖

圖2 β1的MCMC抽樣值軌跡和核密度圖

圖3 β2的MCMC抽樣值軌跡和核密度圖

圖4 β3的MCMC抽樣值軌跡和核密度圖

圖5 σ的MCMC抽樣值軌跡和核密度圖

圖6和圖7分別給出了樣本容量為300時,在0.01分位點下,各參數(shù)抽樣值的自相關圖。從圖中可以看出,隨著滯后期的增加,自相關系數(shù)逐漸趨近于零。

圖6 β0，β1的自相關圖

圖7 β2，β3的自相關圖

由于馬爾科夫鏈收斂,因此各參數(shù)的估計結果可靠?；谕瑯拥臉颖?我們還計算了0.05分位點下各參數(shù)的Gibbs抽樣算法的估計值。實驗結果表明,0.05分位點的各參數(shù)的抽樣值軌跡以及自相關圖與0.01分位點的結果類似,馬爾科夫鏈均收斂。參數(shù)的抽樣值估計圖及自相關圖略。表2給出了0.05分位點下各參數(shù)的Gibbs抽樣的估計值。

從表1和表2的估計結果看,1%VaR以及5%VaR的GJR-CAViaR模型的所有參數(shù)的后驗均值都位于95%的后驗置信區(qū)間內,說明在5%顯著性水平下各參數(shù)是顯著的。也就是說,中國證券市場VaR具有自回歸的特征,即當前風險的大小會受到前一期風險的影響。且存在收益對風險的不對稱效應,負收益對VaR的影響更大。

表2 95%置信水平下模型參數(shù)的Gibbs抽樣估計值

（二）不同樣本容量的估計結果

為了考察不同樣本容量的貝葉斯估計結果,對于全樣本區(qū)間內的數(shù)據(jù),分別選擇最近的500個、800個、1 000個以及2 000個樣本估計GJR-CAViaR模型的參數(shù),并計算99%置信水平以及95%置信水平下的參數(shù)估計值。

表3 99%置信水平下不同樣本容量的參數(shù)的Gibbs抽樣估計值

續(xù)前表

參數(shù)的MCMC抽樣軌跡以及抽樣值的自相關性顯示,所有的參數(shù)抽樣構成的馬爾科夫鏈都是收斂的,參數(shù)的估計結果都是可靠的。

從表3和表4的MCMC參數(shù)估計結果看,所有參數(shù)的后驗均值都落在95%的后驗置信區(qū)間內,因此,5%顯著性水平下各參數(shù)是顯著的。β1的值都比較高,β2是正數(shù),且β3顯著大于零。說明在小樣本的情況下,仍然可以得出與大樣本相同的結論,即中國證券市場VaR具有自回歸的特征,當前風險的大小受前一期風險值的影響較大。且存在收益對風險的不對稱效應,負收益對VaR的影響更大。而且,95%和99%置信水平下模型參數(shù)的估計結果顯示,所有估計參數(shù)的標準差都比較小。由于參數(shù)采用的是近乎無信息先驗的正態(tài)分布,計算出的GJR-CAViaR模型的參數(shù)標準差都比較小,也說明了貝葉斯方法的有效性。

表4 95%置信水平下不同樣本容量的參數(shù)的Gibbs抽樣估計值

五、結論

通過不對稱拉普拉斯分布,本文探討了基于Gibbs抽樣的GJR-CAViaR模型的MCMC的算法實現(xiàn),并利用上證綜指的日收益率數(shù)據(jù)進行實證分析,發(fā)現(xiàn)中國證券市場VaR具有自回歸性質,且呈現(xiàn)出收益對風險的不對稱特征,且該性質不受樣本容量和置信水平的影響。

本文關于不同樣本容量的Gibbs抽樣估計值的討論為參數(shù)的結構變化提供了思路。在樣本分布不變的情況下,隨著樣本量的增加,參數(shù)的估計精度會隨之提高。當我們的樣本從500、800、1 000增加到2 000的時候,參數(shù)的抽樣分布的標準差并未出現(xiàn)隨著樣本量增加而越來越小的情況。而且,不同的樣本容量下,參數(shù)的估計值也有不同。如樣本量從1 000增加到2 000的時候,5%VaR的參數(shù)估計值的差異比較大。因此,在這個區(qū)間樣本的分布并不是保持不變的。此時,選擇基于樣本容量為1 000的GJR-CAViaR模型去估計第1 001個交易日的VaR,或者選擇基于樣本容量為2 000的GJR-CAViaR模型去估計第1 001個交易日的VaR是值得思考的問題。在樣本分布沒有發(fā)生改變的情況下,樣本容量2 000的模型應該是更好的選擇。如果樣本在這個區(qū)間內發(fā)生了比較大的結構變化,也許樣本容量較小而日期更近的模型會有更不錯的結果。在CAViaR模型的估計樣本區(qū)間內是否需要加入結構參數(shù),如何確定區(qū)間內樣本是否產生結構變化的問題,將是我們未來的研究方向。

［1］Koenker R,Machado J A.Goodness of fit and related inference processes for quantile regression［J］.Journal of The American Statistical Association,1999,94(448):1296－1310.

［2］Yu K,Moyeed R A.Bayesian quantile regression［J］.Statistics＆ Probability Letters,2001,54(4):437 －447.

［3］Hsu Y H.Applications of quantile regression to estimation and detection of some tail characteristics ［D］.University Of Illinois At Urbana-Champaign,2010.

［4］曾惠芳.基于MCMC算法的貝葉斯分位回歸計量模型及應用研究 ［D］.湖南大學,2011.

［5］Chou R Y T.Forecasting financial volatilities with extreme values:the conditional autoregressive range(CARR)model ［J］.Journal of Money,Credit And Banking,2005,37(3),561－582.

［6］Chen C W S,Gerlach R,Hwang B B K,Mcaleer M.Forecasting Value-at-Risk using nonlinear regression quantiles and the intra-day range ［J］.International Journal of Forecasting,2012(3),557－574.

［7］陳磊,杜化宇,曾勇.基于貝葉斯CAViaR模型的油價風險研究 ［J］.系統(tǒng)工程理論與實踐,2013(11),2757－2765.

［8］王新宇,吳祝武,宋學鋒.變點CAViaR市場風險測量模型及創(chuàng)業(yè)板應用 ［J］.中國礦業(yè)大學學報,2013(03),506－512.

［9］張穎,孫和風,張富祥.中美股票市場風險差異的新解釋——收益對市場風險不對稱效應的CAViaR模型與實證 ［J］,南開經濟研究,2012(05),111－120.

［10］Koenker R,Bassett,G.Regression Quantiles ［J］.Econometrica,1978,46(1),33 －50.

［11］Tsionas E G.Bayesian quantile inference ［J］.Journal of Statistical Computation And Simulation,2003,73(9),659 －674.

［12］Kozumi H,Kobayashi G.Gibbs sampling methods for Bayesian quantile regression ［J］.Journal of Statistical Computation And Simulation,2011,81(11),1565－1578.

［13］王新宇,宋學鋒,吳瑞明.基于AAVS-CAViaR模型的股市風險測量研究 ［J］.系統(tǒng)工程學報,2010(03),326－333.

［14］邸俊鵬.分位數(shù)回歸的貝葉斯估計與應用研究 ［D］,南開大學,2013.