舒 婷， 羅幼喜， 李翰芳

(湖北工業(yè)大學理學院， 湖北武漢430068)

在貝葉斯框架下，利用Lasso懲罰的方法對面板數據中固定效應與隨機效應進行降維，從而進行參數估計與變量選擇已經是學術界討論的熱門問題之一。在統計模型參數估計中，除經典的損失框架、極大似然框架、經驗似然框架之外就是基于貝葉斯框架。英國學者Thomas Bayes于1764年提出貝葉斯統計[1]，之后Laplace進一步完善貝葉斯統計理論。與似然估計框架不同的是，貝葉斯將所有待估參數都看作某種分布，通過給定條件先驗分布，利用貝葉斯估計原理推導出各參數對應的條件后驗分布，以系統的統計推斷方法進行參數估計，即使是在觀測時間減少或者樣本數據量較小的面板數據下，也能得到精確的估計值。Jaynes[2-3]也相繼展開對貝葉斯統計理論的研究。如今，貝葉斯統計方法不斷被廣大學者應用于理論研究與解決實際問題之中。

待估計變量的懲罰函數不同，即以不同的懲罰函數對其進行壓縮時得到的估計值也有所不同，從而有不同的變量選擇結果。Lasso懲罰由Tibshirani[4]于1996年提出，與常用的SCAD懲罰相比，Lasso懲罰降低了模型的計算復雜程度。由于Lasso懲罰能對變量進行自動選擇，Park等[5]將其引入貝葉斯分層模型中，提出待估參數的有條件Laplace先驗，為變量選擇開辟了一條新的理論道路。楊磊等[6]通過這種貝葉斯分層模型，利用Lasso懲罰實現自動化參數的調整，避免了繁瑣的調參過程。在對固定效應與隨機效應系數施加Lasso懲罰后，待估參數的壓縮系數是相同的，而在實際數據分析中往往不同的變量其重要程度也不盡相同，這就表明對應的壓縮程度不能都一致。因此，Zou[7]提出了Adaptive Lasso懲罰函數，對每一個待估參數系數添加不同的權重，從變量選擇和估計的效果來看，該方法能對待估參數進行不同程度的壓縮，大大提高模型估計的準確率。Li等[8]利用貝葉斯方法與Adaptive Lasso懲罰，分別對線性和logistic回歸模型進行參數后驗分布近似并證明其具有漸近性和Oracle性質。李翰芳等[9]研究含固定效應的一般線性模型的變量選擇問題，在貝葉斯框架下引入Lasso懲罰函數對固定效應系數進行估計，李子強等[10]、王小燕等[11]、張巧琦[12]在此基礎上用Adaptive Lasso懲罰對固定效應的系數進行不同程度的壓縮。目前，無論是基于貝葉斯框架，還是似然函數框架，利用不同懲罰函數對固定效應的變量選擇與估計問題已被廣泛討論，但是針對含隨機效應的面板數據模型的變量選擇與參數估計問題研究還相對較少。

對于含隨機效應的復雜面板數據模型而言，由于解釋變量受隨機效應不同程度的影響，其變量選擇與估計的難度遠遠大于只含固定效應的一般線性模型。本文提出一種基于貝葉斯框架下新的雙懲罰分位回歸方法，在參數估計的同時進行變量選擇。首先，建立貝葉斯分層分位回歸模型，對固定效應與隨機效應同時添加Adaptive Lasso懲罰，使不同待估參數得到不同程度的壓縮，從而建立貝葉斯雙Adaptive Lasso分位回歸模型(BALQR)；其次，在給定各待估參數的條件先驗分布后，構造固定效應與隨機效應參數估計的切片Gibbs抽樣算法；最后，利用R語言進行編程實現蒙特卡羅模擬比較研究，對比了不同隨機效應、不同誤差分布情形下變量選擇和參數估計的效果。

1 模型與方法

1.1 建立含隨機效應的面板數據貝葉斯分層分位回歸模型

含混合效應的線性模型定義如下

(1)

式中:xij為第i個個體在第j次時間點下的觀測值；β是解釋變量的固定效應向量；αi為第i個個體的隨機效應向量；zij是與隨機效應相對應的協變量,擾動項εij分布未知。

對于式(1)中的響應變量yij，其條件分位回歸函數為：

(2)

(3)

(4)

由于非對稱Laplace分布沒有共軛先驗，本文利用Neal[13]提出的正態(tài)指數分解法將ALD分布表示為正態(tài)和指數2個常見分布的混合，在后面構造的Gibbs抽樣算法中，由于常見分布利于數據的抽取，降低參數估計的計算復雜度，同時也減小混合效應模型的估計誤差。

(5)

(6)

假設模型(6)中β～π(β)，αi～π(αi)，σ～π(σ)均有先驗分布，且ξij=ξi1,ξi2,…,ξiMT,i=1,…,N，IG為逆Gamma分布。

由式(4)～(6)可推導出如下貝葉斯分層模型：

(7)

1.2 建立貝葉斯雙Adaptive Lasso分位回歸模型

在Alhamzawi等[14]研究的Lasso懲罰下，由于固定效應β和隨機效應αi都分別具有依賴的條件參數η1和η2，但對于每一個不同的β和αi而言，其壓縮系數應該是不同的。因此用Lasso懲罰就會限制固定效應與隨機效應變動的靈活性，為了得到更好的參數估計結果，本文建立貝葉斯雙Adaptive Lasso分位回歸模型。

首先假設固定效應β與隨機效應αi都有條件Laplace先驗如下：

(8)

(9)

由式(3)、(8)、(9)，可以推導出固定效應β與隨機效應αi的后驗分布為

(10)

(11)

要使式(11)中固定效應與隨機效應聯合后驗對數密度函數極大化，也即使式(12)雙Adaptive Lasso貝葉斯分位回歸函數極小化。

(12)

1.3 參數估計的切片Gibbs抽樣算法

在上述建立的雙Adaptive Lasso懲罰模型之中，直接抽取樣本是很難實現的，為了更加有效地進行參數估計的Gibbs抽樣[15]，需要進一步的理論推導。

對于固定效應β和隨機效應αi，式(8)、(9)可改寫為：

(13)

(14)

利用文獻[16]中的積分恒等式

(15)

可將式(13)固定效應的先驗分布改寫為

(16)

式中R=r1,…,rlT，通過引入輔助變量R,將式(16)的聯合條件先驗分布轉換成常見分布的混合，即正態(tài)分布與指數分布的組合。

(17)

其他未知參數的后驗條件密度函數可由式(13)～(17)進一步推導得出。

對于rl,l=1,…,k，

(18)

式中GIGρ,m,n為廣義逆高斯分布(generalized inverse Gaussian, GIG)。

(19)

對于固定效應β,

(20)

同理，對于隨機效應αi，由積分恒等式(15)將式(14)改寫為

(21)

(22)

(23)

對于隨機效應αi，

(24)

對于σ的先驗密度函數為π(σ)～IG(e0,f0)，有

π(σyij,β,αi,ξij,τ)∝L(yijβ,αi,ξij,σ,τ)π(ξijσ)πσ～IG(ν,ω),

(25)

對于ξij：

(26)

因此，可構造出的切片Gibbs抽樣算法為：

③ 從條件后驗分布πσyij,β,αi,ξij,τ～IGν,ω中生成σ；

⑥ 從條件后驗分布πβyij,αi,ξ,σ,R∝NkBb,B中生成β；

⑨ 從條件后驗分布παi|yij,β,σ,ξ,S～Nk(Aa,A)中生成αi；

⑩ 從②開始往后直至迭代收斂。

2 蒙特卡羅模擬

本文首先在不同分位點、不同類型解釋變量系數下，模擬貝葉斯分位回歸(BQR)[17]、貝葉斯雙Lasso分位回歸(BLQR)[9]、貝葉斯雙Adaptive Lasso分位回歸(BALQR) 3種方法的估計結果，并進行對比研究；然后在不同方法、不同強度的隨機效應干擾下進行對比；最后研究不同方法、不同分布的隨機誤差對不同類型面板數據模型估計的影響。為了評價模型估計的準確性，本文選取多次模擬均方誤差(MSE)的均值作為評價指標，MMSE為不同估計方法下的均方誤差，定義為

(27)

2.1 不同方法對估計的影響

模擬一保持ρ=0.5、D=diag1,1,0、εij～N0,1不變，表1和表2分別為兩組解釋變量系數在不同分位點、不同方法下的估計結果，其中包含MMSE與各解釋變量系數估計的平均值、標準差及其對應的區(qū)間估計。

表1 稀疏面板數據下3種方法的估計結果

表2 稠密面板數據下3種方法的估計結果

雙Adaptive Lasso分位回歸方法與已有研究不同之處在于固定效應與隨機效應都有對應的壓縮系數，從而更大程度地對無關變量進行剔除。本文突破已有文獻中僅針對固定效應的變量選擇進行研究，考慮不同隨機效應影響下變量的估計與選擇。

根據表1，針對稀疏型面板數據展開討論，在τ=0.25、0.5、0.75處，BALQR與BLQR的MMSE均值相差不大，前者比后者略低一些，意味著2種方法都能得到精確的估計結果。在高分位點τ=0.95處，3種方法的MMSE均值明顯增大，但BALQR在100次模擬中得到的MMSE均值、標準差最小，區(qū)間估計的范圍更加精確。

根據表2，在稠密型面板數據的條件下，在τ=0.25、0.5、0.75處，BALQR與BLQR的MMSE均值比BQR要小，在該條件下BQR的參數估計與變量選擇效果沒有BLQR與BALQR好。從MMSE的均值與標準差來看，BALQR與BLQR的估計效果相差不大，與表1所得結果類似。在高分位點τ=0.95處，3種方法的均方誤差都相對變大，BQR與BLQR的MMSE均值分別為0.057、0.055，兩者差距不大，而此時BALQR的MMSE均值為0.038，說明其在高分位點處的估計效果最好。

綜上，本文可以得出BALQR方法無論在低分位點還是高分位點處，其參數估計與變量選擇效果都最優(yōu)。

2.2 隨機效應對估計的影響

模擬二保持ρ=0.5、εij～N0,1不變，通過不斷增加隨機效應的干擾強度來比較不同分位點下3種方法的估計結果，將其協方差分別設為：D1=diag(1,1,0)、D2=diag(2,2,0)、D3=diag(3,3,0)展開討論。由于各分位點的估計結果類似，表3、表4僅展示0.5分位點下的模擬結果。

表3 不同隨機效應影響下稀疏面板數據的估計結果

表4 不同隨機效應影響下稠密面板數據的估計結果

根據表3，在稀疏型面板數據模型中，隨機效應的干擾強度由弱到強，對比D1=diag(1,1,0)、D2=diag(2,2,0)、D3=diag(3,3,0) 3種條件下，BQR的MMSE均值分別為0.043、0.065、0.083；BLQR的MMSE均值分別為0.035、0.053、0.072；BALQR的MMSE均值分別為0.023、0.042、0.065。3種方法對應的MMSE均值逐步降低，這是因為Lasso分位回歸相比分位回歸而言，對每一個解釋變量施加Lasso懲罰，從而提高模型計算的速度降低估計的偏差，而Lasso懲罰對每一個β都有相同的壓縮系數，本文采用的Adaptive Lasso懲罰函數突破了這一限制，從模擬結果來看，在不同干擾強度下、BALQR方法對稀疏型面板數據模型能進行有效的參數估計，從而進行變量選擇。

根據表4，針對稠密型面板數據而言，在D1=diag(1,1,0)條件下，BLQR與BALQR的MMSE均值分別為0.027、0.021。隨著協方差增加到D2=diag(2,2,0)時，BLQR與BALQR都能得到相對精確的估計效果。但當協方差增加到D3=diag(3,3,0)時，BQR與BLQR的MMSE均值分別為0.082、0.069，其標準差也顯著上升，說明在隨機效應的干擾強度不斷增強時，BQR與BLQR的估計偏差較大。此時BALQR的MMSE均值為0.042，標準差的均值為0.016，表明BALQR在強隨機效應干擾狀態(tài)下對稠密型面板數據的估計效果比對稀疏型面板數據的估計效果要好。

就MMSE均值這一指標而言，隨著隨機效應的干擾強度增加，其均值和標準差也在增加，意味著固定效應的估計精度在下降，尤其是β1與β2，這是由于本文在模擬時假定前面2個變量受隨機效應影響的結果。綜合兩組解釋變量系數的模擬結果來看，BALQR的估計效果最好，其變量選擇與估計的優(yōu)勢比BQR、BLQR方法更加突出。

2.3 不同分布的隨機誤差對估計的影響

模擬三保持ρ=0.5、D=diag(1,1,0)不變，本文將考慮隨機誤差分別服從標準正態(tài)分布、t(3)分布、ALD(0, 0.5, 1)分布下，模擬雙Adaptive Lasso分位回歸的面板數據變量選擇與估計結果?；谪惾~斯框架的優(yōu)勢在于未知參數都可以看成服從某一先驗條件分布，在給定不同待估參數的先驗條件分布之后，利用Gibbs抽樣算法進行變量選擇與估計，表5、6為0.5分位點下的估計結果。

根據表5，從不同分布不同方法對稀疏型面板數據的估計結果來看，當隨機誤差服從標準正態(tài)分布時，BQR、BLQR、BALQR 3種方法估計的MMSE均值分別為0.043、0.035、0.023，3種方法都能得到偏差較小的參數估計結果，但BALQR更勝一籌，且BALQR對應的MMSE標準差最小，為0.017，說明100次模擬結果的波動較小，其得到精確估計的效果較穩(wěn)定。當隨機誤差服從t(3)分布時，BQR與BLQR的MMSE均值顯著增加，分別為0.060、0.052，說明在t(3)分布下，這兩種估計方法沒有在標準正態(tài)分布下精確，而此時BALQR對應的MMSE均值為0.025，充分說明BALQR估計效果具有巨大的優(yōu)勢。當隨機誤差服從ALD(0,0.5,1)時，BQR與BLQR的估計效果幾乎一樣，均方誤差都比BALQR高，此時BALQR的MMSE均值為0.020，即在稀疏型面板數據中，隨機誤差分別服從3種不同分布下，BALQR都能保持穩(wěn)定的精確的估計結果。

根據表6，在稠密型面板數據模型中，當隨機誤差服從標準正態(tài)分布時，BLQR與BALQR都能得到較為精確的估計結果，相反BQR的誤差項偏高。當隨機誤差服從t(3)分布時，與表5的結果相似，BQR與BLQR的MMSE均值也顯著增加，表明這2種方法的估計效果不佳，而BALQR的MMSE均值為0.032，在該方法下能得到偏差較小的參數估計。對比稀疏型面板數據而言(表5)，BALQR在t(3)分布下對應的MMSE均值為0.025，說明雙Adaptive Lasso分位回歸在處理稀疏型面板數據方面更有優(yōu)勢。當隨機誤差服從ALD(0,0.5,1)時，BQR與BLQR的MMSE均值都為0.037，而BALQR的為0.024，即BALQR的估計效果最優(yōu)。

表 5 不同分布下稀疏面板數據的估計結果

表 6 不同分布下稠密面板數據的估計結果

綜上，不論隨機誤差是服從標準正態(tài)分布還是t(3)分布、ALD(0,0.5,1)分布，貝葉斯雙Adaptive Lasso分位回歸方法的估計誤差最小。針對不同類型的面板數據而言，BALQR對含隨機效應的混合線性模型能得到更好的變量選擇結果。

3 實例分析

使用本文提出的貝葉斯雙Adaptive Lasso分位回歸方法對來自29家民營上市公司15年(2003-2017年)的凈利潤數據進行研究分析。凈利潤是上市公司的稅后利潤，是衡量上市公司經營效益的重要指標之一。民營上市公司相關數據來源于CSMAR數據庫，本文初步選取來自上市公司的股本結構、資產價值、盈利能力、潛在價值這4個方面共計12個指標作為解釋變量，如表7所示，凈利潤為響應變量。接下來本文研究各個解釋變量在不同分位點上對凈利潤的貢獻程度，以及利用雙Adaptive Lasso懲罰對不同重要程度的變量進行選擇。

表 7 解釋變量描述

考慮到各個解釋變量之間存在多重共線性，想要分析清楚總股數、股東總數、流動資產、營業(yè)收入等到底如何影響凈利潤的大小，一般的線性回歸模型難以解決這一問題，因此我們建立如下的分層隨機效應面板數據模型：

(28)

下面將在4個分位點上，對29家民營上市公司的凈利潤值建立貝葉斯雙Adaptive Lasso分位回歸模型，以此來分析各個解釋變量的貢獻率。表8為各分位點下的估計結果。

表 8 BALQR在4個分位點處的估計結果

表8中各分位點下，對上市公司凈利潤影響較大的變量有所不同，按權重系數大小將低于0.1的變量自動剔除，重要變量的估計均值加粗表示。在低分位點τ=0.25處，雙Adaptive Lasso變量選擇出5個貢獻率較大的變量，重要程度按絕對值從大到小依次為X8、X6、X5、X12、X4；在中分位點τ=0.5處，變量選擇結果按重要程度依次為X10、X8、X6、X11；在較高分位點τ=0.75處，變量選擇結果按重要程度依次為X10、X8、X1、X6、X7；在高分位點τ=0.95處，變量選擇結果按重要程度依次為X8、X12、X7、X6、X2、X1。

不難看出扣除非經常性損益后的凈利潤、股東權益在每個分位點上對凈利潤的影響力較強，說明股東權益是影響凈利潤的關鍵指標。總股數在低分位點和中分位點都被剔除，但在高分位點變成重要變量之一，說明隨著總股數的增加對提高上市公司凈利潤有巨大的潛力。股東總數、年度股東大會出席率、每股凈資產在各分位點大多為負值，說明要保持上市公司凈利潤值穩(wěn)定上升，需要減少凈利潤值對其的依賴程度，這一結果為公司調整內部結構提供重要的理論依據。

4 結論

本文提出的貝葉斯雙Adaptive Lasso分位回歸方法對固定效應與隨機效應參數共同給定非條件Laplace先驗，使得不同的參數壓縮程度不同，從而達到更加精確的估計結果。無論是稀疏型面板數據還是稠密型面板數據，貝葉斯雙Adaptive Lasso分位回歸方法都能得到精確的參數估計并同時進行變量選擇。

通過對比研究BQR、BLQR、BALQR 3種分位回歸方法，在變量選擇方面，3種懲罰分位回歸方法對冗余變量的系數幾乎都壓縮到零,從而進行精確的變量選擇，但是對于固定效應系數的估計卻不盡相同，BALQR的估計誤差明顯小于BQR與BLQR的估計誤差。在模擬過程中各參數的核密度圖顯示誤差基本集中在零左右，表明BALQR下Gibbs抽樣算法已達到收斂，限于文章篇幅，在此不做展示。在不同隨機效應的干擾強度下，該方法仍舊可以得到誤差最小的參數估計值。

貝葉斯分層模型的極大優(yōu)勢在于對來源于不同先驗分布的參數進行估計，由于隨機誤差的產生人為不可控，其分布也是不明確的，一般都是人為設定。本文從不同分布的角度出發(fā)，對比BQR、BLQR、BALQR的估計結果，對于不同的面板數據類型，在不同分位點上貝葉斯雙Adaptive Lasso方法的預測結果最佳。