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

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

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

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

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

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

1　引言

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

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

2　模型介紹

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

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

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

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

3　Bayesian分析及抽樣方法

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

由于參數(shù)的先驗信息很難得到，而不恰當?shù)南闰灂δM結果產(chǎn)生影響.因此，參數(shù)β，γ，α，λ的先驗皆選為無信息先驗.另外，為確保方差hit為正和模型的平穩(wěn)性，系數(shù)λ和α必須滿足λi≥0，αi≥0和

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

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

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

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

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

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

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

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

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

參數(shù)β，γ和σ2的后驗條件分布都為標準分布，可直接進行抽樣.

隱含變量uit的后驗分布隨著時間變化，第i個公司t時刻無效率項uit的后驗分布為

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

文中令

其中

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

式

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

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

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

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

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

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

其中，

那么，λ建議分布為

這里，

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

上式中α?是α的加權非線性最小二乘估計值，即

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

則有

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

令

則

4　模擬實驗

為驗證Gibbs抽樣方法的效果，這里進行了模擬實驗.為方便計算，令p=q=1并對兩個公司的情況進行模擬，考慮如下模型：

取k=3，m=3，xit，zit的第一列均取1，為表現(xiàn)出相關性，xit，zit的剩余k-1列的元素從二元標準正態(tài)分布中隨機產(chǎ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.對于先驗分布的參數(shù)分別取p0=4，q0=0.01.分別取樣本量100，500和1000三種情況來進行模擬比較，模擬軌道長度60000并拋去軌道的前10000個點.各參數(shù)后驗估計值和后驗標準差見表1.

表1　參數(shù)后驗均值結果

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

5　實證研究

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

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

表2　參數(shù)后驗均值結果

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

［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隨機前沿模型的Bayesian推斷［J］.高校應用數(shù)學報，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

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