程　迪，張世斌（.內(nèi)蒙古大學(xué)數(shù)學(xué)科學(xué)學(xué)院，內(nèi)蒙古呼和浩特000；.上海海事大學(xué)數(shù)學(xué)系，上海浦東0306）

動(dòng)態(tài)異方差隨機(jī)前沿模型的Bayesian推斷

程迪1，張世斌2
（1.內(nèi)蒙古大學(xué)數(shù)學(xué)科學(xué)學(xué)院，內(nèi)蒙古呼和浩特010021；2.上海海事大學(xué)數(shù)學(xué)系，上海浦東201306）

隨機(jī)前沿模型中如果忽略單邊干擾項(xiàng)的異質(zhì)性（heterogeneity）往往導(dǎo)致錯(cuò)誤的效率估計(jì).從個(gè)體特征的影響和方差的時(shí)變性?xún)煞矫鎸?duì)單邊干擾項(xiàng)進(jìn)行考慮，提出異方差動(dòng)態(tài)隨機(jī)前沿模型.利用Gibbs抽樣方法對(duì)動(dòng)態(tài)異方差隨機(jī)前沿模型進(jìn)行Bayesian分析.導(dǎo)出了模型參數(shù)的后驗(yàn)條件分布，對(duì)中小樣本的模擬實(shí)驗(yàn)顯示在最小后驗(yàn)均方誤差準(zhǔn)則下得到的參數(shù)估計(jì)值非常接近真值.對(duì)電力公司的實(shí)際數(shù)據(jù)進(jìn)行分析顯示對(duì)數(shù)無(wú)效率項(xiàng)的方差有一定的時(shí)變性.

隨機(jī)前沿模型；Bayesian分析；異方差；Gibbs抽樣；Metropolis-Hastings抽樣

1　引言

隨機(jī)前沿模型由Aigner等（1977）［1］及Meeusen和Broeck（1977）［2］提出，被廣泛用于衡量經(jīng)濟(jì)效率.隨機(jī)前沿模型將誤差項(xiàng)分為對(duì)稱(chēng)誤差項(xiàng)和單邊干擾項(xiàng)兩部分，分別反映了測(cè)量誤差和模型中的技術(shù)無(wú)效率對(duì)產(chǎn)出或成本造成的負(fù)面影響，即單邊干擾項(xiàng)又為無(wú)效率項(xiàng).通常假定無(wú)效率項(xiàng)在所有決策單元中是獨(dú)立同分布的，無(wú)效率項(xiàng)的分布一般假設(shè)為半正態(tài)分布［1］，指數(shù)分布［2］，截?cái)嗾龖B(tài)分布［3］，Gamma分布［4］或倒Gamma分布［5］等.然而，影響一個(gè)公司效率的因素除了成本函數(shù)模型中的輸入輸出變量，還包括一些公司特征，如公司的所有權(quán)，公司政策的改變，行業(yè)的競(jìng)爭(zhēng)等.這些因素將對(duì)公司的效率產(chǎn)生一定的影響.如果公司的這些特征產(chǎn)生異質(zhì)性（heterogeneity）或者對(duì)單邊干擾項(xiàng)產(chǎn)生影響，這些特征就可以在前沿模型中表現(xiàn)出來(lái).文［6］將協(xié)變量直接以函數(shù)形式表現(xiàn)并研究了單邊干擾項(xiàng)中未被觀(guān)察到的異質(zhì)性.文［7］針對(duì)多輸出隨機(jī)前沿模型，提出用copula函數(shù)來(lái)刻畫(huà)無(wú)效率項(xiàng)分布.文［8］提出了動(dòng)態(tài)隨機(jī)前沿模型，無(wú)效率項(xiàng)隨著時(shí)間的變化而變化，并用Bayesian分析和MCMC方法對(duì)參數(shù)進(jìn)行估計(jì).文［9］將隨機(jī)前沿模型進(jìn)行轉(zhuǎn)換，消去個(gè)體的固有影響，并對(duì)其參數(shù)估計(jì)的一致性進(jìn)行分析.文［10］對(duì)無(wú)效率項(xiàng)服從廣義Gamm a分布和廣義混合Gamma分布的隨機(jī)前沿模型進(jìn)行Bayesian分析，考慮了無(wú)效率項(xiàng)中可觀(guān)察到的異方差.文［11］提出用參數(shù)來(lái)刻畫(huà)未被觀(guān)測(cè)到的異質(zhì)性，指出了異質(zhì)性可能影響的變量并通過(guò)實(shí)例對(duì)變量進(jìn)行分析論證.

現(xiàn)有隨機(jī)前沿模型較少涉及無(wú)效率項(xiàng)中未觀(guān)測(cè)到的異質(zhì)性，本文提出用GARCH（p，q）來(lái)刻畫(huà)對(duì)數(shù)無(wú)效率項(xiàng).其中，均值方程刻畫(huà)可以觀(guān)測(cè)到的異質(zhì)性，條件方差方程刻畫(huà)未觀(guān)測(cè)到的異質(zhì)性.當(dāng)p=q=0時(shí)對(duì)數(shù)無(wú)效率項(xiàng)方差為常數(shù)，無(wú)效率項(xiàng)不存在異方差.否則，對(duì)數(shù)無(wú)效率項(xiàng)方差是過(guò)去誤差項(xiàng)和滯后條件方差的函數(shù)，隨時(shí)間變化方差的值不同.本文利用Gibbs抽樣對(duì)動(dòng)態(tài)異方差隨機(jī)前沿模型進(jìn)行Bayesian分析.通過(guò)對(duì)模型參數(shù)選取合理的先驗(yàn)分布，導(dǎo)出模型參數(shù)的后驗(yàn)條件分布，并給出了Gibbs抽樣的具體策略.由于參數(shù)λ和α的后驗(yàn)條件分布不是標(biāo)準(zhǔn)密度函數(shù)，且極難找到標(biāo)準(zhǔn)分布的密度函數(shù)作為控制函數(shù)進(jìn)行舍選，常用的隨機(jī)游走抽樣和G riddy-Gibbs抽樣方法在該分層隨機(jī)前沿模型中的舍選效率非常低.本文參照Nakatsum a（參見(jiàn)［12］）提出的近似GARCH模型方法，將條件方差方程構(gòu)造成新的ARMA方程，得到參數(shù)的建議分布，最后通過(guò)Metropolis-Hastings抽樣實(shí)現(xiàn)對(duì)參數(shù)λ和α的抽樣.此方法大大提高了動(dòng)態(tài)異方差隨機(jī)前沿模型的計(jì)算效率.

2　模型介紹

與截面數(shù)據(jù)相比，面板數(shù)據(jù)更能提供對(duì)公司技術(shù)效率的可靠估計(jì).考慮隨機(jī)前沿模型

其中yit表示第i個(gè)個(gè)體t時(shí)間的產(chǎn)出，xit為k×1維向量，表示第i個(gè)個(gè)體t時(shí)間的投?入，β￠是k×1的參數(shù)向量.vit是統(tǒng)計(jì)誤差，設(shè)vit關(guān)于i和t相互獨(dú)立并且有相同的分布vit～N 0，σ2.uit是無(wú)效率項(xiàng)，為非負(fù)的隨機(jī)變量.假設(shè)vit和uit相互獨(dú)立，且無(wú)效率項(xiàng)uit滿(mǎn)足

對(duì)所有的i，t都有ξit～N（0，1）且相互獨(dú)立，zit為m×1維的協(xié)變量，γ為m×1維的參數(shù)向量.這里γ用來(lái)描述無(wú)效率項(xiàng)對(duì)數(shù)中可以觀(guān)測(cè)到的異質(zhì)性，未能觀(guān)測(cè)到的異質(zhì)性服從N（0，hit）.方差hit為GARCH（p，q）模型，當(dāng)λ1=···=λp=α1=···=αq=0時(shí)，hit恒為常數(shù)λ0，即無(wú)效率項(xiàng)uit不存在異方差，否則方差hit隨時(shí)間變化，依賴(lài)于其以前時(shí)刻的方差和式（1）中的擾動(dòng)項(xiàng)平方.現(xiàn)有文獻(xiàn)中較多的涉及無(wú)效率項(xiàng)中已觀(guān)測(cè)到的異質(zhì)性，而對(duì)未觀(guān)察到的異質(zhì)性研究較少，文［8］和［11］中假設(shè)未觀(guān)測(cè)到的異質(zhì)性服從一定分布.本文在已有文獻(xiàn)的基礎(chǔ)上，用協(xié)變量表示無(wú)效率中可以觀(guān)測(cè)到的異質(zhì)性.對(duì)未觀(guān)察到的異質(zhì)性，通過(guò)假設(shè)對(duì)數(shù)無(wú)效率項(xiàng)方差為以前時(shí)刻方差和擾動(dòng)平方項(xiàng)的線(xiàn)性組合，討論無(wú)效率項(xiàng)中未觀(guān)測(cè)到的異質(zhì)性及其變化規(guī)律.

y和u的聯(lián)合分布為

3　Bayesian分析及抽樣方法

本節(jié)對(duì)動(dòng)態(tài)異方差隨機(jī)前沿模型進(jìn)行Bayesian分析并采用Gibbs抽樣方法對(duì)模擬參數(shù)的后驗(yàn)均值進(jìn)行推斷.首先確定參數(shù)的先驗(yàn)分布，先驗(yàn)分布的選取應(yīng)滿(mǎn)足相容性原則，即隨著樣本量的增大，先驗(yàn)信息的影響減弱.假設(shè)參數(shù)β，γ，α，λ和尺度參數(shù)σ2相互獨(dú)立，則參數(shù)β，γ，α，λ的聯(lián)合先驗(yàn)具有形式

由于參數(shù)的先驗(yàn)信息很難得到，而不恰當(dāng)?shù)南闰?yàn)會(huì)對(duì)模擬結(jié)果產(chǎn)生影響.因此，參數(shù)β，γ，α，λ的先驗(yàn)皆選為無(wú)信息先驗(yàn).另外，為確保方差hit為正和模型的平穩(wěn)性，系數(shù)λ和α必須滿(mǎn)足λi≥0，αi≥0和

那么，參數(shù)β，γ，α，λ的聯(lián)合先驗(yàn)可以寫(xiě)為

其中，I（·）為示性函數(shù).

尺度參數(shù)σ2的先驗(yàn)分布可取倒Gamma分布（參見(jiàn)［13］），即

其中，p0＞0和q0＞0為先驗(yàn)分布的參數(shù).

由Bayesian理論，將參數(shù)的先驗(yàn)分布和樣本信息綜合可得到參數(shù)后驗(yàn)分布.θ和u的聯(lián)合后驗(yàn)分布為

為對(duì)模型進(jìn)行預(yù)測(cè)或參數(shù)推斷，Bayesian模型往往需要對(duì)高維概率密度分布進(jìn)行積分. MCMC是一種解決高維積分問(wèn)題的迭代Monte Carlo方法，它解決了復(fù)雜表達(dá)式難以進(jìn)行高維積分的問(wèn)題，應(yīng)用非常廣泛.MCMC通過(guò)利用M arkov鏈進(jìn)行積分得到一系列的M arkov鏈模擬值，該鏈的平穩(wěn)分布收斂到已知參數(shù)的后驗(yàn)分布，即可看作后驗(yàn)分布的獨(dú)立樣本.Gibbs抽樣是一種最常見(jiàn)且最簡(jiǎn)單的MCMC方法，這里使用Gibbs抽樣方法對(duì)參數(shù)后驗(yàn)均值進(jìn)行推斷.為了使MCMC方法產(chǎn)生的Markov鏈的平穩(wěn)分布是（2），下面討論Gibbs抽樣方法的具體策略.設(shè)Gibbs抽樣得到的Markov鏈

收斂到已知的參數(shù)后驗(yàn)分布.因此，參數(shù)后驗(yàn)期望可以通過(guò)Gibbs抽樣得到的平均值來(lái)近似，即

下面對(duì)參數(shù)抽樣的具體策略進(jìn)行分析.由（2）式可得到每個(gè)參數(shù)后驗(yàn)分布的核.

其中，x為（N·T）×k矩陣，y為（N·T）×1矩陣.由（3）式可知，

參數(shù)β，γ和σ2的后驗(yàn)條件分布都為標(biāo)準(zhǔn)分布，可直接進(jìn)行抽樣.

隱含變量uit的后驗(yàn)分布隨著時(shí)間變化，第i個(gè)公司t時(shí)刻無(wú)效率項(xiàng)uit的后驗(yàn)分布為

其后驗(yàn)分布不是標(biāo)準(zhǔn)分布，不能進(jìn)行直接抽樣，可以通過(guò)切片抽樣或舍選抽樣得到.文中使用的是切片抽樣方法，產(chǎn)生uit的具體步驟見(jiàn)算法1.

文中令

其中

從后驗(yàn)分布可知，在給定參數(shù)u，γ的前提下，參數(shù)λ和α中給定任意一個(gè)參數(shù)的值，另一個(gè)參數(shù)的條件后驗(yàn)核密度函數(shù)都不是標(biāo)準(zhǔn)分布密度函數(shù).在這里取樣方法參照Nakatsum a（參見(jiàn)［12］）所提議的近似GARCH模型，用Metropolis-Hastings抽樣方法對(duì)其進(jìn)行抽樣.令

式

其中l(wèi)=m ax｛p，q｝，且有

以λ的取樣為例，Metropolis-Hastings抽樣方法的具體步驟為：

算法2（產(chǎn)生λ）

步1從建議分布π（λ）中抽取λ′.

步2從U（0，1）中抽取e.

算法2中λ的建議分布需要主觀(guān)設(shè)定，近似GARCH的似然函數(shù)可以寫(xiě)為：

其中，

那么，λ建議分布為

這里，

為生成α的建議分布，將wit（α）在α?處用Tay lor展式展至一次項(xiàng)，則近似GARCH模型的似然函數(shù)為

上式中α?是α的加權(quán)非線(xiàn)性最小二乘估計(jì)值，即

其中ζitj是wit在點(diǎn)處的一階導(dǎo)數(shù)值的相反數(shù).又

則有

那么，參數(shù)α的建議分布為

令

則

4　模擬實(shí)驗(yàn)

為驗(yàn)證Gibbs抽樣方法的效果，這里進(jìn)行了模擬實(shí)驗(yàn).為方便計(jì)算，令p=q=1并對(duì)兩個(gè)公司的情況進(jìn)行模擬，考慮如下模型：

取k=3，m=3，xit，zit的第一列均取1，為表現(xiàn)出相關(guān)性，xit，zit的剩余k-1列的元素從二元標(biāo)準(zhǔn)正態(tài)分布中隨機(jī)產(chǎn)生，兩者的相關(guān)系數(shù)為0.5，如文［14］.取參數(shù)β真值為［1，1，1］′，γ真值為［-0.3，1，1］′.取σ2=0.05，λ0=0.05，λ1=0.2，α1=0.7.對(duì)于先驗(yàn)分布的參數(shù)分別取p0=4，q0=0.01.分別取樣本量100，500和1000三種情況來(lái)進(jìn)行模擬比較，模擬軌道長(zhǎng)度60000并拋去軌道的前10000個(gè)點(diǎn).各參數(shù)后驗(yàn)估計(jì)值和后驗(yàn)標(biāo)準(zhǔn)差見(jiàn)表1.

表1　參數(shù)后驗(yàn)均值結(jié)果

由表1，對(duì)比各個(gè)參數(shù)的真值和估計(jì)值，參數(shù)β的估計(jì)值在樣本量很小時(shí)就非常接近真實(shí)值，參數(shù)z的估計(jì)值樣本量很小時(shí)就很接近真值.參數(shù)λ0，λ1，α1和σ2在樣本量小時(shí)接近真值，參數(shù)的估計(jì)隨著樣本量的增大，模擬效果變好.參數(shù)的后驗(yàn)標(biāo)準(zhǔn)差隨著樣本量的增大明顯減小，特別是參數(shù)λ0，λ1和α1，在樣本量較小時(shí)標(biāo)準(zhǔn)差較大，隨著樣本量的增大明顯減小.可見(jiàn)，用Gibbs抽樣對(duì)動(dòng)態(tài)異方差隨機(jī)前沿模型進(jìn)行Bayesian分析有比較好的效果.在模擬實(shí)驗(yàn)過(guò)程中發(fā)現(xiàn)，當(dāng)樣本量較小特別是公司多時(shí)間短的情況下，λ0的模擬值會(huì)出現(xiàn)均非常接近于零的情況，使hit的后驗(yàn)均值幾乎為零，即參數(shù)估計(jì)結(jié)果中無(wú)效率項(xiàng)方差估計(jì)值非常小.這是因?yàn)?，?的抽樣采用的是M etropolis-Hasting算法.在抽樣過(guò)程中，一旦在一步中出現(xiàn)λ0的模擬值很接近于零時(shí)，λ0的建議分布的方差將很小，使得Metropolis-Hasting抽樣更新效率極低.而當(dāng)樣本抽樣時(shí)間區(qū)間很短時(shí)，極易出現(xiàn)在某一步迭代中λ0的模擬值非常接近于零的情況.將λ0的先驗(yàn)由無(wú)信息先驗(yàn)變?yōu)榫鶆蚍植糢［δ，1］，其中δ為一很小的正數(shù)，可避免λ0的模擬值幾乎為零的情況.

5　實(shí)證研究

將本文模型應(yīng)用于哥倫比亞電力改革.樣本數(shù)據(jù)來(lái)自文［15］中，取自1998-2012年哥倫比亞的電力分銷(xiāo)公司，本文對(duì)原始數(shù)據(jù)進(jìn)行篩選，取出其中時(shí)間序列較長(zhǎng)的12個(gè)公司共174個(gè)觀(guān)察值.考慮如下隨機(jī)前沿模型，

其中，yit，x1it，x2it分別表示第i個(gè)公司t時(shí)刻的輸入總費(fèi)用，輸出總能量和消費(fèi)者數(shù)量.無(wú)效率項(xiàng)異方差變量zit為第i個(gè)公司t時(shí)刻的用戶(hù)密度.數(shù)據(jù)和模型的進(jìn)一步解釋可參見(jiàn)文［16］.先驗(yàn)參數(shù)取值與模擬實(shí)驗(yàn)部分的值一致，p和q的值均取1.G ibbs抽樣模擬軌道長(zhǎng)60000并拋去模擬軌道的前10000個(gè)點(diǎn).主要參數(shù)的后驗(yàn)均值和標(biāo)準(zhǔn)差見(jiàn)表2，其中括號(hào)中數(shù)字為參數(shù)后驗(yàn)標(biāo)準(zhǔn)差.

表2　參數(shù)后驗(yàn)均值結(jié)果

由表2可知，從γ的估計(jì)值可以看出，公司的用戶(hù)密度會(huì)對(duì)效率產(chǎn)生一定的影響.參數(shù)λ0的估計(jì)值為0.011，而λ1和α1的估計(jì)值分別為0.192和0.160，說(shuō)明對(duì)數(shù)無(wú)效率項(xiàng)的方差依賴(lài)于過(guò)去誤差項(xiàng)和滯后條件誤差.可見(jiàn)，隨機(jī)前沿模型中不能忽略公司特征引起的異方差和方差的時(shí)變性.實(shí)例主要研究了條件方差hit為GARCH（1，1）的情況，未對(duì)高階的情況進(jìn)行深入研究.另外，可以將本文與文［8］和［11］進(jìn)行結(jié)合，對(duì)無(wú)效率項(xiàng)進(jìn)行更全面的詮釋.

［1］Aigner D，Knox-Lovell C A，Schm dt P.Formu lation and estimation of stochastic frontier production function m odels［J］.J Econom etrics，1977，6:21-37.

［2］M eeusen W，van der B roeck J.E fficiency estim ation from Cobb-Douglas p roduction functions w ith Com posed Error［J］.Internat Econom Rev，1977，18:435-444.

［3］Stevenson R E.Likelihood functions for generalized stochastic frontier model［J］.J Econom etrics，1980，13:57-66.

［4］G reene W H.A gamm a-d istributed stochastic frontier m odel［J］.J Econom etrics，1994，61: 273-303.

［5］劉曉君，張世斌.正態(tài)倒Gamm a隨機(jī)前沿模型的Bayesian推斷［J］.高校應(yīng)用數(shù)學(xué)報(bào)，2013，28（4）: 488-496.

［6］Greene W.Reconsidering heterogeneity in panel data estimators of the stochastic frontier model［J］.J Econometrics，2005，126:269-303.

［7］Carta A，SteelM F J.M odellingmu lti-output stochastic frontiers using copu las［J］.Com put Statist Data Anal，2012，56:3757-3773.

［8］Tsionas E.In ference in dynam ic stochastic frontier models［J］.J App l Econom，2006，21: 669-676.

［9］Chen Y i-Y i，Schm id t P，W ang Hung-Jen.Consistent estim ation of the fixed effects stochastic frontiermodel［J］.J Econometrics，2014，181:65-76.

［10］Griffin J，Steel M.Flexiblem ixturemodelling of stochastic frontiers［J］.J Prod Anal，2008，29:33-50.

［11］Galán J，Veiga H，W iper M.Bayesian estimation of inefficiency heterogeneity in stochastic frontiermodels［J］.J Prod Anal，2014，42:85-101.

［12］Nakatsuma T.Bayesian Analysis of ARMA-GARCH Models：A Markov chain Sam pling Approach［J］.J Econometrics，2000，95:57-69.

［13］Fernandez C，Osiewalski J，Steel M F J.On the use of panel data in stochastic frontier models［J］.J Econometrics，1997，79:169-193.

［14］Wang Hung-Jen，Schm id t P.One step and two step estimation of the effects of exogenous variables on technical efficiency levels［J］.Prod Anal，2002，18:129-144.

［15］Galán J E，Pollitt M.Inefficiency persistence and heterogeneity in Colombian electricity utilities［J］.Energ Econ，2014，46:31-44.

［16］Grow itsch T，Jamasb T，Pollitt M.Quality of service，efficiency and scale in network industries:an analysis of European electricity d istribu tion［J］.App l Econ，2005，41:2555-2570.

M R Sub ject C lassification:62F15；65C60

Bayesian in ference for dynam ic heterogeneity stochastic fron tier m odel

CHENG D i1，ZHANG Shi-bin2
（1.School of M ath.Sci.，Inner M ongolian Univ.，Hohhot 010021，China；2.Dept.of M ath.，Shanghai M aritim e Univ.，Shanghai201306，China）

If heterogeneity of the“inefficiency”term is disregarded，it w ill resu lt in the incorrect estim ate of this term in the stochastic frontier m odel.By com bining the in fluence from characteristic differencesof individualsw ith the time-varying property of variance，a dynam ic heterogeneity stochastic frontier m odel is p roposed.By the G ibbs sam p ling，the m ethodology for Bayesian analysis of the dynam ic heterogeneity stochastic frontier model is given.For each model parameter，the posterior distribu tion is derived.A sim ulation study show s that under the criterion ofm inim izing the posterior mean square error，the Bayesian estimate is close to its true value for small and medium sized sam p les. From the Bayesian analysis based on the real electric pow er com pany generation data，it is evidenced that there exists the time-varying p roperty for the variance of the logarithm“inefficiency”term.

stochastic frontiermodel；Bayesian inference；heterogeneity；Gibbs sam p ling；Metropolis-Hastings sam p ling

O212.8

A

1000-4424（2016）02-0127-09

2015-11-13

2016-04-19

上海市自然科學(xué)基金（13ZR1419100）；上海市教委科研創(chuàng)新項(xiàng)目（14YZ115）