基于Gibbs抽样方法的空间滞后随机前沿模型Bayesian估计

窦剑军，王 媛，张辉国，胡锡健

（1.兰州财经大学陇桥学院，兰州 730101；2.新疆大学 数学与系统科学学院，乌鲁木齐 830046）

0 引言

每个生产单元追求的是低投入高产出的利润最大化，这就要求生产单元将技术效率发挥到最大化，对于每个生产单元面临的问题是如何测算它们的技术效率。作为技术效率测算的参数方法的随机前沿模型应运而生，该模型于 1977 年被 Meeusen[1]，Aigner[2]与 Battese[3]提出后被广泛地应用于诸多行业技术效率的测算中。经过四十多年的发展，一系列关于随机前沿模型的设定和模型的参数估计方法、技术效率估计推断以及将随机前沿模型应用到具体行业测算技术效率的成果不断出现，为每个行业技术效率的测算提供的一种全新的技术方法。

在随机前沿模型发展的初期主要研究经典的随机前沿模型，其假定不同的生产单元相互独立。但是经济、知识、技术效率往往具有空间上的溢出效应，技术效率高的生产单元将会集聚到一起，技术效率低的生产单元将会集聚到一起，形成“马太效应”。所以，地理位置上的邻接和经济上的关联是产生空间效应的重要因素。所以经典的随机前沿模型在技术效率的测算过程中存在很大的缺陷。

近年来，随着产业结构的调整，产业出现了集聚化的发展，生产单元与生产单元之间的联系越来越紧密，空间交互效应发挥着重要的作用。将空间交互效应引入到经典随机前沿模型中形成空间随机前沿模型，其理论研究层出不穷并且将该方法应用具体行业的技术效率测算中。在模型参数估计方面，张进峰（2014）[4]应用二阶段最小二乘法和矩估计方法对空间滞后随机前沿模型作了估计。在应用方面，陈关聚（2014）[5]应用随机前沿技术测算了中国制造业30个行业的全要素能源效率，张月玲等（2015）[6]应用随机前沿模型研究了人力资本结构、适宜技术选择与全要素生产率变动的分解，陈洁（2016）[7]应用随机前沿模型测算了我国各省份用电技术效率。

上述估计方法可以得到模型参数的一致估计，但对参数的显著性进行假设检验时检验统计量的构造比较困难。鉴于此，目前应用贝叶斯方法对空间滞后随机前沿模型的参数做估计还没有出现，本文用贝叶斯方法对模型参数的估计进行了研究。贝叶斯方法的最大优点在于得到了模型参数的后验分布，在进行参数的显著性假设检验时检验统计量比较容易得到。

1 模型介绍

考虑如下随机前沿模型：

其中y是n个生产单元的产出构成的一个n×1维的向量，X是一个n×(m+1)的设计矩阵，由每个生产单元的m个投入构成的矩阵，W是n个单元的空间权重矩阵。v是不可控的随机因素，设v～N(0，)，u是技术无效率项，设u～N+(0，)，v与u相互独立，在u已知的条件下，y服从正态分布。

在模型(1)中待估计的参数为ρ，β，，和技术无效率项u,为了对参数进行贝叶斯估计，需要选择参数的先验分布。与的先验分布选择为倒Gamma分布，即取：

取ρ和β的先验分布为：

π(β)∝const

其中λmin，λmax为空间邻接矩阵的最小特征值和最大特征值。

由y=ρWy+Xβ+v-u和u已知的条件下得v=(I-ρW)y+u-Xβ，根 据v～N(0，) 有所以

由贝叶斯公式，经计算得参数联合后验分布：

其中R+表示正实数。

从上可以得到，模型参数的联合后验分布是一个高维分布，得到每个参数的后验分布需要对联合后验作积分，该积分有可能不存在或没有办法得到。基于此，将采用马尔科夫链蒙特卡罗(MCMC)方法中的Gibbs抽样方法对模型参数的后验均值进行推断。Gibbs抽样是一种解决高维Bayesian模型的后验积分困难问题的迭代Monte Carlo方法，它避免了对复杂表达式进行积分的问题。Gibbs抽样得到的样本{ρ(i)，β(i)，，，u(i);i=1，2，…n}随着i的增大将会收敛到已知核的参数后验分布。因此所有关于参数的向量函数的后验期望可以通过Gibbs抽样得到的样本平均值来近似。即：

2 参数的条件分布及其抽样方法

为了进行抽样，首先需要推导出每个参数的后验条件分布，在给定u的条件下，由公式(2)可以得到每个参数后验条件分布的核。这里记D={y，W，X}。β的后验条件分布为：

由上式可以看出β后验条件分布为多元正态分布。

的后验条件分布为：

恰好为倒Gamma分布密度函数的核，即：

的后验条件分布为：

这恰好也为倒Gamma分布密度函数的核，即：

参数ρ的条件后验分布为：

这个后验条件分布不再是一个常见的分布形式，对参数ρ的抽样需要采用Metropolis-Hastings方法(M-H)。对于M-H抽样，需要一个建议分布，从中产生参数ρ的一个候选点，记为ρ*，将这个候选值与当前值ρc代入上式中去计算马尔科夫链转移概率：

假设ρ*从一个对称的提议分布g(ρ*|ρc)=g(ρc|ρ*)中产生，则在每一次迭代中，从g(·)中产生一个增量Z，然后ρ*=ρc+Z,比如增量Z可以从标准正态分布中产生，此时候选点ρ*|ρc～N(ρc，σ2)，σ2＞ 0

产生ρ的算法如下：

（1）构造合适的建议分布g(·|ρc)；

（2）从 建 议 分 布g(·) 中 产 生 一 个 增 量 Z，令ρ*=ρc+Z；

（3）重复上述步骤(直到马氏链达到平稳状态）

(a)从 U(0,1)中产生 u;

(c)增加t，返回到(a)。

参数u的条件分布为：

可见u的后验条件分布服是n元截尾正态分布N+(μ，Ω)，可直接对其进行抽样。

3 技术效率的估计

Jondrow等(1982)[8]研究得到各个生产单元的技术无效率项ui可以通过E(ui|εi)或Mode(ui|εi)来测算，第i个生产单 元 的 技 术 效 率TEi=exp{-ûi} ，其 中=E(ui|εi) 或Mode(ui|εi)。因此，只要求得u的后验条件分布的均值或众数即可。因为u的后验条件分布为截尾正态分布，所以

每个生产单元的技术效率为：

4 模特卡罗模拟

为了验证解决Bayesian估计中高维积分问题的Gibbs抽样方法的效果，本文进行了数值模拟，考虑如下数据生成过程：

表1 参数后验均值模拟结果

从表1可以看出，参数的模拟值非常接近真实值（其中第一列中括号里的为真实值），该方法在模拟过程中取得了比较好的效果，因此Bayesian方法能够很好地对空间滞后随机前沿模型作出估计。

5 结论与展望

本文基于Gibbs抽样方法对空间滞后随机前沿模型进行了Bayesian推断，得到了模型参数的后验条件分布和技术效率的估计。Gibbs抽样方法是一种解决高维Bayesian模型积分问题的后验迭代Monte Carlo方法，它避免了对联合后验分布进行数值积分计算的问题。模拟结果表明，参数的估计值不受空间权重矩阵的影响，并且估计值十分接近真实值。

本文仅考虑因变量的空间滞后，没有考虑误差的空间自相关。沿用本文的方法和思想，可以把该方法应用到到空间滞后误差自相关随机前沿模型的参数估计中。同时可以把技术无效率的分布扩展到指数分布、截尾正态分布等其他形式的分布。

参考文献：

[1]Meeusen W,Broeck J V D.Effciency Estimation from Cobb-Douglas Production Functions With Composed Error[J].International Econom⁃ic Review,1977,（18）.

[2]Aigner D,Lovell K，Schmidt P.Formulation and Estimation of Stochas⁃tic Frontier Production Function Models[J].Journal of Economet⁃rics,1977,（6）.

[3]Battese G，Corra G.Estimation of a Production Frontier Model:With Application to the Pastoral Zone of Eastern Australia[J].Australian Journal of Agricultural Economics,1977,（21）.

[4]张进峰.空间滞后随机前沿模型的估计研究[J].商业经济与管理，2014,（8）.

[5]陈关聚.中国制造业全要素能源效率及影响因素研究——基于面板数据的随机前沿分析[J].中国软科学，2014，(1).

[6]张月玲，叶阿忠，陈泓.人力资本结构、适宜技术选择与全要素变动分解——基于区域异质性随机前沿生产函数的经验分析[J].财经研究，2015,41(6).

[7]陈洁.我国用电技术的随机前沿分析[J].生态经济，2016,32(8).

[8]Jondrow J,Lovell K,Materov I，et al.On the Estimation of Technical Ineffciency in the Stochastic Frontier Production Function Model[J].Journal of Econometrics,1982,（19）.