双边随机前沿模型的最大似然估计法拟合【Use R】

QQ461157910
双边随机前沿模型，也称双边随机边界模型，近年来在计量经济中常用于对市场上供求双方的议价能力进行测度。例如 Kumbhakar等在一篇发表于2009年的论文中运用双边随机边界模型分析劳动力市场工资的定价行为和劳资双方的议价能力。
双边随机前沿模型基本结构如下图。

其中ω为正向随机边界，μ为负向随机边界，分别代表市场活动中博弈或议价的双方。ν是随机误差。一般令ω、μ服从参数分别为Σω和Σμ的指数分布，ν则为高斯分布。而μ(X)=Xβ，X为解释变量，y为被解释变量。
随机前沿模型一般用MLE方法来拟合，双边随机前沿模型也不能免俗。下面的代码示例中，我们先根据双边随机前沿模型生成一组虚拟数据。

x=data.frame(x1=rnorm(200,mean=3,sd=1),x2=rnorm(200,mean=8,sd=1.02),x3=rnorm(200,mean=4,sd=1.5))
beta=c(0.37,-0.12,0.25)
muta=as.matrix(x)%*%as.matrix(beta)
y=muta+rexp(200,2)-rexp(200,1.5)+rnorm(200,0,1)

然后写出该模型的对数似然函数。

ll=function(p){
  sigmu=p[1]
  sigom=p[2]
  sigmv=p[3]
  
  if(sigmu<=0|sigom<=0|sigmv<=0)
    return(-Inf)
  
  beta=p[4:6]
  
  eps=y-(as.matrix(x)%*%as.matrix(beta))
  
  a=sigmv^2/(2*(sigmu^2))+(eps/sigmu)
  b=sigmv^2/(2*(sigom^2))-(eps/sigom)
  h=(eps/sigmv)-(sigmv/sigom)
  c=-(eps/sigmv)-(sigmv/sigmu)
  
  n=length(eps)
  
  return(-n*log(sigmu+sigom)+sum(log(exp(a)*pnorm(c)+exp(b)*pnorm(h))))

}

最后用maxLik包内的maxLik函数对模型参数进行最大似然估计。这里参数有6个，3个β，3个Σ。

library(maxLik)

a <- maxLik(ll, start=c(0.2,1,0.5,0.3,0.1,0.2),method = 'BFGSR')

--------------------------------------------
Maximum Likelihood estimation
BFGSR maximization, 36 iterations
Return code 3: Last step could not find a value above the current.
Boundary of parameter space?  
Consider switching to a more robust optimisation method temporarily.
Log-Likelihood: -377.3712 
6  free parameters
Estimates:
     Estimate Std. error t value  Pr(> t)    
[1,]  1.87428    0.77887   2.406  0.01611 *  
[2,]  1.40025    0.35599   3.933 8.38e-05 ***
[3,]  0.87413    0.19707   4.436 9.18e-06 ***
[4,]  0.34218    0.12443   2.750  0.00596 ** 
[5,] -0.09927    0.05973  -1.662  0.09652 .  
[6,]  0.25861    0.08741   2.958  0.00309 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

估计值十分接近真实值，表明程序正确。
最大似然估计本质上是一个最优化问题，因此除了maxLik外，R中其他用于最优化的包也能完成MLE。
此外，可以尝试先用差分进化等全局最优化算法求得一个初始参数集S，再以S为起点，用maxLik对S进一步优化。对于示例中的简单数据而言，这样做毫无必要，但对于研究中遇到的真实数据，用上述两步法常能得到更好的参数组合。
除了在R中手工编程外，stata也有一个第三方命令SFA2tier可以拟合双边随机前沿模型，但据说该命令易用性比较一般。