基于Gibbs-DA算法的貝葉斯分位回歸模型研究*

朱慧明，翁 元，曾昭法，任英華，李招來

(1．湖南大學(xué) 工商管理學(xué)院，湖南 長(zhǎng)沙 410082；2．湖南大學(xué) 金融與統(tǒng)計(jì)學(xué)院，湖南 長(zhǎng)沙 410079)

回歸分析模型是實(shí)踐中應(yīng)用廣泛的一類統(tǒng)計(jì)模型，如果模型中的隨機(jī)擾動(dòng)項(xiàng)來自均值為零而且同方差的分布，回歸系數(shù)的OLS估計(jì)具有最佳線性無偏性,并且，如果隨機(jī)擾動(dòng)項(xiàng)服從正態(tài)分布，回歸系數(shù)的OLS估計(jì)量和ML估計(jì)為最小方差無偏估計(jì)．但是，在現(xiàn)實(shí)社會(huì)經(jīng)濟(jì)系統(tǒng)中，這些假設(shè)條件通常難以得到滿足，例如，數(shù)據(jù)可能出現(xiàn)尖峰或厚尾分布和異方差等問題，此時(shí)OLS估計(jì)量不再具有上述優(yōu)良性．為了彌補(bǔ)OLS估計(jì)方法的不足之處， Koenker和Bassett[1]提出了分位回歸模型的建模思想．與OLS方法相比，分位回歸通過加權(quán)誤差絕對(duì)值之和最小化方法得到參數(shù)的估計(jì)，這些估計(jì)量不容易受到異常值的影響，穩(wěn)健性強(qiáng)．分位回歸模型在金融風(fēng)險(xiǎn)度量、時(shí)間序列和生存分析等領(lǐng)域中獲得了廣泛應(yīng)用[2-5]．

分位回歸模型參數(shù)估計(jì)方法很多，例如，對(duì)偶單純形估計(jì)算法、內(nèi)點(diǎn)法和外點(diǎn)法等．但是，它們是建立在經(jīng)典統(tǒng)計(jì)理論基礎(chǔ)之上，參數(shù)是固定不變常數(shù)．為了考慮參數(shù)不確定性問題，Yu和Moyeed[6]，Tony和Sung[7]，Tsionas[8]和曾惠芳[9]等利用MCMC抽樣方法構(gòu)建貝葉斯分位模型．但是，它們沒有研究模型參數(shù)的隨機(jī)抽樣問題．本文利用Gibbs抽樣和數(shù)據(jù)擴(kuò)展(Data Augmentation，DA)算法，構(gòu)建基于Gibbs-DA抽樣算法的貝葉斯分位回歸模型，解決模型參數(shù)不確定性風(fēng)險(xiǎn)問題．

1 非對(duì)稱Laplace分布的混合表示

定義1如果隨機(jī)變量U具有如下密度函數(shù)：

(1)

則稱U服從位置參數(shù)μ，尺度參數(shù)σ和偏度參數(shù)τ的非對(duì)稱Laplace分布，記作U～ALD(μ,σ,τ)．特別地，當(dāng)μ=0,σ=1時(shí)，記作U～ALD(τ)．如果隨機(jī)變量U服從非對(duì)稱Laplace分布ALD(μ,σ,τ)，則它的數(shù)學(xué)期望和方差如下：

U的特征函數(shù)為:

(2)

非對(duì)稱Laplace分布與指數(shù)分布、正態(tài)分布之間存在密切關(guān)系：服從非對(duì)稱Laplace分布的隨機(jī)變量可以用指數(shù)分布、正態(tài)分布的隨機(jī)變量的函數(shù)來表示．如果U1,U2分別服從指數(shù)分布Exp(1)和正態(tài)分布N(0,1)，并且兩者相互獨(dú)立，記:

(3)

(4)

這就是μ=0,σ=1時(shí)非對(duì)陣Laplace分布的特征函數(shù)．根據(jù)特征函數(shù)的唯一性定理可知：U～ALD(τ)．

2 貝葉斯分位回歸算法

假設(shè)y與p維變量x之間具有線性函數(shù)關(guān)系：

y=x'β+ε

(5)

此處β=(β1,β2,…,βp)'，隨機(jī)誤差項(xiàng)ε的概率分布函數(shù)為F．顯然，對(duì)于給定的x，隨機(jī)變量y的τ分位數(shù)等于x'β+F-1(τ)，即:

Qy(τ|x,β)=x'β+F-1(τ)

(6)

其等價(jià)形式為:

Qy(τ|x,β(τ))=x'β(τ)

(7)

此處β(τ)=(β1(τ),β2(τ),…,βm(τ))'，也就是說，β(τ)與分位數(shù)τ的取值大小有關(guān)．

假設(shè)(x1,y1),(x2,y2),…,(xn,yn)是(x,y)的一組觀察值，則β(τ)的估計(jì)為:

(8)

其中ρτ(u)=u(τ-I(u<0))為損失函數(shù)，I(u<0)為示性函數(shù)．由于

(9)

等價(jià)于

(10)

而后者恰好是非對(duì)稱Laplace變量的聯(lián)合分布密度函數(shù)的核，因此，根據(jù)非對(duì)陣Laplace分布的性質(zhì)，分位回歸建模問題可以進(jìn)一步轉(zhuǎn)化為

i=1,2,…,n

(11)

中的參數(shù)β(τ)；此處zi,ui相互獨(dú)立，zi服從標(biāo)準(zhǔn)指數(shù)分布Exp(1)，ui服從標(biāo)準(zhǔn)正態(tài)分布；模型(8)稱為數(shù)據(jù)擴(kuò)展模型，zi稱為潛變量(latent variable)．

則模型(8)可以簡(jiǎn)化為如下矩陣形式：

(12)

L(Y|X;Z,β(τ))

(13)

由于Z的分量相互獨(dú)立，并且均服從標(biāo)準(zhǔn)指數(shù)分布Exp(1)，其概率分布密度函數(shù)

(14)

為了進(jìn)行貝葉斯分析，需要進(jìn)一步設(shè)置分位回歸模型參數(shù)的先驗(yàn)分布．根據(jù)文獻(xiàn)[10]的觀點(diǎn)，選擇如下均值向量β0，協(xié)方差陣Σ0的多元正態(tài)分布作為模型參數(shù)β(τ)的先驗(yàn)分布，即

β(τ)～Np(β0,Σ0)，β0∈Rp,Σ0>0

(15)

此處β0,Σ0可以由Koenker的經(jīng)典分位回歸方法確定．根據(jù)貝葉斯定理，潛變量和參數(shù)的后驗(yàn)分布密度與似然函數(shù)、先驗(yàn)分布密度的乘積成正比，即:

(16)

為了能夠進(jìn)行MCMC抽樣算法，研究潛變量和參數(shù)完全條件后驗(yàn)分布．

1)參數(shù)β(τ)的完全條件后驗(yàn)分布．通過π(β(τ),Z|Y,X)關(guān)于β(τ)的積分運(yùn)算，獲得Z的邊緣后驗(yàn)分布密度函數(shù)π(Z|Y,X)，據(jù)此計(jì)算β(τ)的完全條件后驗(yàn)分布密度函數(shù)π(β(τ)|Y,X'Z)．不難驗(yàn)證：參數(shù)β(τ)的完全條件后驗(yàn)分布密度函數(shù)為

(17)

顯然，β(τ)的完全條件后驗(yàn)分布是均值μβ(τ)，協(xié)方差陣Στ的正態(tài)分布，即

(β(τ)|Y,X;Z)～Np(μβ(τ),Στ)

(18)

2)潛變量U的完全條件后驗(yàn)分布．同樣地，通過π(β(τ),Z|Y,X)關(guān)于Z的積分運(yùn)算，求出參數(shù)β(τ)的邊緣后驗(yàn)分布密度函數(shù)π(β(τ)|Y,X)，據(jù)此進(jìn)一步計(jì)算Z的完全條件后驗(yàn)分布密度函數(shù)π(Z|Y,X;β(τ))．對(duì)于給定的Y,X和β(τ)，Z的完全條件后驗(yàn)分布密度函數(shù)為:

(19)

(zi|Y,X;β(τ))～GIG(1/2,χi,ψ),

i=1,2,…,n

(20)

其分布密度函數(shù)為

π(Z|Y,X;β(τ))

其中

根據(jù)潛變量Z和模型參數(shù)β(τ)的完全后驗(yàn)條件分布，確定貝葉斯分位回歸模型的MCMC抽樣步驟：

Step 1．設(shè)置參數(shù)β(τ)的初始值β(0)(τ)，以及樣本量lN．并且，假設(shè)(Z(j-1),β(j-1)(τ))為(Z,β(τ))的第j-1次抽樣結(jié)果;

以上是基于Gibbs-DA抽樣的單鏈條算法．當(dāng)然，為了便于分析Markov鏈的平穩(wěn)性，也利用Gibbs-DA抽樣同時(shí)生成β(τ)的多條鏈．例如，給定β(τ)的m個(gè)初始值{β(i0)(τ),i=1,2,…,m}，生成如下m個(gè)長(zhǎng)度為lN的Markov鏈：

然后，根據(jù)Gelman-Rubin收斂診斷檢驗(yàn)統(tǒng)計(jì)量分析Markov鏈{β(ij)(τ),i=1,…,m；j=1,…,lN}的收斂性問題．如果Gelman-Rubin檢驗(yàn)統(tǒng)計(jì)量的取值介于1.0到1.2之間，則鏈條達(dá)到平穩(wěn)狀態(tài)，此時(shí)利用樣本{(Z(i, lj),β(i, lj)(τ)),i=1,2,…,m；j=M+1,M+2,…,N}，對(duì)模型參數(shù)進(jìn)行統(tǒng)計(jì)推斷，求出模型參數(shù)點(diǎn)估計(jì)及置信水平1-α的區(qū)間估計(jì)，并對(duì)模型進(jìn)行擬合分析．

3 實(shí)證研究

作為上述方法的應(yīng)用，我們考慮我國(guó)歷年能源消耗彈性的貝葉斯分位回歸問題．模型變量包括能源消耗總量(y，單位：萬噸標(biāo)準(zhǔn)煤)、國(guó)內(nèi)生產(chǎn)總值(GDP，單位：億元)，以及第二產(chǎn)業(yè)所占比重(P，單位：%)數(shù)據(jù)研究．?dāng)?shù)據(jù)來源于歷年 《中國(guó)統(tǒng)計(jì)年鑒》、《中國(guó)能源統(tǒng)計(jì)年鑒》，涉及全國(guó)30個(gè)省市自治區(qū)(不含西藏)，同時(shí)，為了便于比較分析，GDP數(shù)據(jù)按2005年價(jià)格計(jì)算．能源消耗總量、國(guó)內(nèi)生產(chǎn)總值和第二產(chǎn)業(yè)所占比重之間存在如下關(guān)系：

ln (y)=β0+β1ln (GDP)+β2ln (P)+ε

此處P=GDP2/GDP×100；β1表示能源消耗的GDP彈性系數(shù)，而β2為能源消耗的產(chǎn)業(yè)結(jié)構(gòu)彈性系數(shù)．

首先，利用2010年全國(guó)30個(gè)地區(qū)的數(shù)據(jù)，分別建立3個(gè)分位點(diǎn)(τ=0.25，0.50和0.75)的貝葉斯分位回歸模型．為此，運(yùn)用Quantreg初步估計(jì)分位回歸模型系數(shù)，并將估計(jì)結(jié)果作為貝葉斯分位回歸模型參數(shù)先驗(yàn)分布的均值向量；同時(shí)，參數(shù)先驗(yàn)分布的協(xié)方差矩陣如下：Σ0=diag(1 000,1 000,1 000)．為獲得模型參數(shù)的貝葉斯估計(jì)，通過運(yùn)用R2WinBUGS軟件，對(duì)每個(gè)參數(shù)模擬生成兩條Markov鏈，每條Markov鏈的迭代次數(shù)均為100 000次，以確保Markov鏈的收斂性．在MCMC模擬初始階段，Markov鏈取值可能與參數(shù)的初始值選擇有關(guān)，為此，舍棄每條鏈的前50 000次抽樣結(jié)果，利用余下的50 000次數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析．表1給出了分位回歸模型系數(shù)的貝葉斯估計(jì)結(jié)果．

表1 分位回歸模型參數(shù)的貝葉斯估計(jì)結(jié)果

根據(jù)表1所列出的計(jì)算結(jié)果，0.25分位點(diǎn)的GDP彈性系數(shù)為0.681 2，也就是說，在產(chǎn)業(yè)結(jié)構(gòu)保持不變的條件下，GDP增加大于1%，能源消耗總量增加0.681 2%；0.75分位點(diǎn)的GDP彈性系數(shù)為0.491 2，低于0.25分位點(diǎn)的估計(jì)值；與能源消耗的GDP彈性系數(shù)情況不同，0.75分位點(diǎn)的產(chǎn)業(yè)結(jié)構(gòu)彈性系數(shù)的取值為1.016 0，大于0.25分位點(diǎn)的產(chǎn)業(yè)結(jié)構(gòu)彈性系數(shù)．為了能夠從時(shí)間維度和不同分位點(diǎn)兩個(gè)方面分析能源消耗總量與GDP、產(chǎn)業(yè)結(jié)構(gòu)之間的關(guān)系，利用2001-2009年的統(tǒng)計(jì)數(shù)據(jù)，根據(jù)前面的建模思路，進(jìn)一步構(gòu)建0.25，0.50和0.75三個(gè)分位點(diǎn)的貝葉斯分位回歸分析模型，合計(jì)27個(gè)模型．

表2匯總了2001-2010年能源消耗的GPD彈性系數(shù)，即：30個(gè)分位回歸模型系數(shù) 的估計(jì)，該表最后一列等于(β1(0.25)-β1(0.75))/β1(0.75)×100. 從時(shí)間維度來看，2001-2010年3個(gè)分位點(diǎn)的GDP彈性系數(shù)變化不大，但是，0.75分位點(diǎn)的GDP彈性系數(shù)明顯小于0.25分位點(diǎn)的彈性系數(shù)．由于0.75分位點(diǎn)代表能源消費(fèi)量較高的地區(qū)，主要是東部地區(qū)；而0.25分位點(diǎn)代表能源消費(fèi)量低的地區(qū)，主要是中、西部地區(qū)．這說明我國(guó)東部整體單位經(jīng)濟(jì)增長(zhǎng)所消耗的能源少于中、西部整體單位經(jīng)濟(jì)增長(zhǎng)所消耗的能源，東部的能源利用率和節(jié)能力度優(yōu)于中、西部.

表2 2001-2010年能源消耗的GDP彈性系數(shù)

表3匯總了2001-2010年能源消耗的產(chǎn)業(yè)結(jié)構(gòu)彈性系數(shù)，即：30個(gè)分位回歸模型系數(shù)β2的估計(jì)，該表的最后一列等于(β2(0.75)-β2(0.25))/β2(0.25)×100. 從時(shí)間維度來看，2001-2010年三個(gè)分位點(diǎn)的產(chǎn)業(yè)結(jié)構(gòu)彈性系數(shù)下降趨勢(shì)明顯，但是，與GDP彈性系數(shù)情況不同，0.75分位點(diǎn)的產(chǎn)業(yè)結(jié)構(gòu)彈性系數(shù)大于0.25分位點(diǎn)的彈性系數(shù)．這說明2001-2010年?yáng)|部的高產(chǎn)出相對(duì)西部來說是依靠高能耗實(shí)現(xiàn)的．

表3 2001-2010年能源消耗的產(chǎn)業(yè)結(jié)構(gòu)彈性系數(shù)

4 結(jié) 論

本文討論了貝葉斯線性分位回歸模型的構(gòu)建、MCMC仿真分析及其應(yīng)用問題．通過模型的統(tǒng)計(jì)結(jié)構(gòu)分析，根據(jù)非對(duì)稱Laplace分布的正態(tài)—指數(shù)分布的混合表示性質(zhì)，利用服從指數(shù)分布的潛變量，據(jù)此獲得模型的似然函數(shù)，推斷了多元正態(tài)先驗(yàn)分布條件下分位回歸模型參數(shù)的后驗(yàn)分布，證明了潛變量的完全條件分布為廣義逆高斯分布．在此基礎(chǔ)上，根據(jù)MCMC仿真分析思想，設(shè)計(jì)了貝葉斯分位回歸模型的Gibbs-DA隨機(jī)抽樣步驟及統(tǒng)計(jì)分析．利用我國(guó)2001-2010年期間各地區(qū)的能源消耗總量，GDP和產(chǎn)業(yè)結(jié)構(gòu)數(shù)據(jù)進(jìn)行實(shí)證分析，研究結(jié)果表明，貝葉斯方法可以有效地應(yīng)用于分位回歸建模問題．

[1]KOENKER R, BASSETT G. Regression quantiles[J]. Econometrica, 1978, 46(1): 33-50.

[2]WOLTERS M H. Estimating monetary policy reaction functions using quantile regressions[J]． Journal of Macroeconomics, 2012, 34(2): 342-361.

[3]BAUR D G , DIMPFL T, JUNG R C. Stock return autocorrelations revisited: A quantile regression approach[J]. Journal of Empirical Finance, 2012, 19(2): 254-265.

[4]FITZENBERGER B F, KOENKER R, MACHADO J A F. Economic applications of quantile regression[M]. Heidelberg: Sprink, 2002：135-148.

[5]羅幼喜,田茂再.面板數(shù)據(jù)的分位回歸方法及其模擬研究[J]. 統(tǒng)計(jì)研究, 2010, 27(10): 81-87.

LUO You-xi,TIAN Mao-zai. Quantile regression for panel data and its simulation study [J]. Statistical Research, 2010, 27(10): 81-87.(In Chinese)

[6]YU K M, MOYEED R A. Bayesian quantile regression[J]. Statistics & Probability Letters, 2001, 54(4): 437-447.

[7]TONY L, SUNG J J. Bayesian quantile regression methods[J]. Journal of Applied Econometrics, 2010, 25(2): 287-307.

[8]TSIONAS E G. Bayesian quantile inference[J]. Journal of Statistical Computation and Simulation，2003，73(9): 659 - 674.

[9]曾惠芳，朱慧明，李素芳，等．基于MH算法的貝葉斯分位自回歸模型[J]．湖南大學(xué)學(xué)報(bào):自然科學(xué)版，2010，37(2)：88-92.

ZENG Hui-fang, ZHU Hui-ming, LI Su-fang,etal．Bayesian inference on the quantile autoregressive models with MH algorithm[J]. Journal of Hunan University: Natural Sciences, 2010，37(2)：88-92. (In Chinese)

[10]ALHAMZAWI R, YU K M. Power prior elicitation in bayesian quantile regression[J]. Journal of Probability and Statistics, 2011, 3(1): 1-16.