极大似然估计 (MLE) 及 Stata 实现

发布时间：2020-09-18 阅读 7200

Stata 连享会主页 || 视频 || 推文

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者：郭李鹏 (山西晋城市财政局)
邮箱： 497724596@qq.com

1. MLE 简介

最大似然估计 ( MLE ) 在计量经济学中有广泛的应用。本文主要介绍 MLE 的基本原理和应用，并演示如何在 Stata 软件中进行最大似然估计。

1.1 MLE 的基本原理

计量经济学中，最小二乘估计、最大似然估计和广义矩估计是构造统计量的三种基本方法。

最小二乘估计，最合理的参数估计量应该使得模型能最好地拟合样本数据，也就是估计值和观测值之差的平方和最小。

矩估计的基本思想是，简单随机样本的原点矩依概率收敛到相应的总体原点矩，用样本矩替换总体矩，进而找出未知参数的估计。

最大似然估计，最合理的参数估计量应该使得从模型中抽取该 n 组样本观测值的概率最大，也就是概率分布函数或者说是似然函数最大。

最小二乘估计和广义矩估计进行统计推断时并不需要设定密度函数，即无需对干扰项分布进行假设；而最大似然估计的前提条件是，能够写出密度函数的完整设定，其中包含一系列待估参数。

1.2 似然函数

假设针对一个包含 $m$ 个随机变量的随机向量 $y_{i} = [y_{1 i}, y_{2 i}, \dots, y_{m i}]^{'}$ ，我们收集了 n 个样本值： $y_{1}$ , $y_{2}$ , $\dots$ , $y_{n}$ ，则这 $n$ 个随机向量的联合密度函数 $L (θ)$ 可以写成

L (θ) = f (y_{1}, y_{2}, \dots, y_{n} | θ),

其中 $θ$ 是一个参数向量，包含了联合密度函数中所涉及到的所有参数。联合密度函数 $L (θ)$ 又称作概似函数，极大化概似函数所得到的解通常以 $\hat{θ}$ 表示，即最大似然估计值。

1.3 对数似然函数

对于分布独立的样本 $x_{1}, x_{2}, \dots, x_{n}$ ，假设其密度函数为 $f_{i} (x_{i} | θ)$ ，经对数转换后的联合密度函数即对数似然函数：

\ln L_{n} (θ) = \sum_{i = 1}^{n} \ln f_{i} (x_{i} | θ)

根据最大似然估计的原理，我们可以通过极大化对数似然函数获得未知参数 $θ$ 的估计值。

1.4 MLE 估计的基本步骤

在 Stata 中进行最大似然估计的基本步骤如下：

推导最大似然函数
编写似然函数的 Stata 程序
设定解释变量和被解释变量，完整设定 ml model 命令
估计最大似然函数：执行 ml maximize 命令

1.5 何时使用 MLE ？

MLE 的主要运用在离散数据模型和一些分布比较复杂的模型，主要有：Logit 模型，Probit 模型，泊松模型，负二项模型，非似然函数模型以及随机边界模型等。

2. MLE 的 Stata 实现

2.1 范例 1：线性回归模型的 MLE 估计

对于线性回归模型，若假设其满足所有基本假设条件，且干扰项服从正态分布， $y$ 服从均值为 $μ$ ，方差为 $σ^{2}$ 的正态分布，可以将 $μ$ 设定为某些变量的线性函数，即 $μ = x β$ 。具体完整形式如下：

y_{i} \sim independent N (x_{i}^{'} β, σ^{2}), i = 1, \dots, n,

矩阵形式为：

y \sim N (X β, σ^{2} I_{n}) .

我们可以通过 MLE 获得该模型中两组参数 $β$ 和 $σ^{2}$ 的估计值。

在用 Stata 进行似然估计前，需要先写出对数似然函数基本设定形式：

\ln L (θ) = \sum_{i = 1}^{n} \ln \frac{1}{σ \sqrt{2 π}} \exp (- \frac{(y_{i} - μ)^{2}}{2 σ^{2}})

\ln L (θ) = \sum_{i = 1}^{n} - \frac{1}{2} \ln (2 π) - \frac{1}{2} \ln (σ^{2}) - \frac{1}{2} \cdot \frac{{(y_{i} - μ)}^{2}}{σ^{2}}

在 Stata 中执行最大似然估计的第一步是将不包括线性回归函数的似然函数的基本设定编写为 ado 文件,需要注意的是，程序中只定义概似函数的基本设定，如何将线性回归函数加入模型由其他管道设定。

cap program drop mymean_lf
program define mymean_lf
    version 9.1        
    args  lnf mu sigma 
    quietly replace `lnf' = ln(normalden(($ML_y1 - `mu')/`sigma')) - ln(`sigma')
end

程序的第一行用于定义程序的名称，对应最后一行的end 语句，表明程序到此结束。

程序中 version 语句声明编写程序时 Stata 版本，以便更高版本能够兼容运行该程序。

程序的第三行用于声明输入项，关键命令是 args，其后的 lnf 是一个暂元，用于存放每次迭代得到的似然函数值，mu 和 sigma 表明我似然函数中有两个线性等式。

quietly replace 语句中 Stata 对每一个观察值求得似然函数值后会按观察值的顺序依次累加后存入暂元 lnf 中。

ML_y1 为默认写法，表示被解释变量，是 Stata 指定的一个全局暂元。如果模型中包含多个被解释变量， Stata 会自动将它们依次存放于 ML_y2、 ML_y3 … ML_yj中。

上述程序在保存时的文件名一定要为 mymean_lf（与我们所定义的程序名称保持一致），保存类型为 ado 文档。

上面的程序只是定义了最大似然估计的“基本架构”，而要实现对具体模型的估计还需要通过“输入项”的设定来定义最大似然估计的“完整架构”。进行程序完整架构定义的指令结构如下：

ml model lf [d0|d1|d2] progname (eq1: y = x1 x2 …) (eq2: x3 x4) (eq3:) ()

ml 方法的关键命令是 ml model 命令和 ml maximize 命令，ml model 命令用来定义需要拟合的模型， ml maximize 命令用来执行最大化。

Stata 对最大似然估计法的执行提供了四种不同的架构， lf 是最常用的架构，通常对样本满足独立同分布的条件时均可使用。在括号中，等号左边为被解释变量，等号右边为解释变量；冒号表示参数名称。

由于 Stata 进行最大似然估计使用的是数值算法，初始值的设定非常重要。我们可以使用如下命令让 Stata 找到进行最大似然估计的初始值。

我们调入汽车价格的数据，演示 Stata 中如何执行 MLE 估计。

sysuse auto, clear
ml model lf mymean_lf (price=mpg wei len) (sigma: )
ml search
ml maximize
est store mle_reg

MLE 线性回归模型与 OLS 结果对比。

reg price mpg wei len
est store ols_reg
esttab mle_reg ols_reg, mtitle(mle_reg ols_reg)

运算结果如下：

--------------------------------------------
                      (1)             (2)   
                  mle_reg         ols_reg   
--------------------------------------------
main                                        
mpg                -86.79          -86.79   
                  (-1.06)         (-1.03)   
weight              4.365***        4.365***
                   (3.84)          (3.74)   
length             -104.9**        -104.9*  
                  (-2.71)         (-2.64)   
_cons             14542.4*        14542.4*  
                   (2.54)          (2.47)   
--------------------------------------------
sigma                                       
_cons              2348.4***                
                  (12.17)                   
--------------------------------------------
N                      74              74   
--------------------------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

从结果可以看出，采用MLE和OLS估计得到的参数没有明显差别，但是 MLE 估计得到的参数的标准差要比 OLS 估计的小。

2.2 范例 2：面板随机边界模型的 MLE 估计

面板随机边界模型是一种测算各类效率的参数模型，完美假设下的产出函数为

Y_{i t} = f (X_{i t}, β)

实际生产中，生产受到各种随机因素的影响，往往还存在效率损失，产出函数为

Y_{i t} = f (X_{i t}, β) w_{i t} e^{v_{i t}}

$w_{i t}$ 表示效率水平，显然，0< $w_{i t}$ <=1。 $v_{i t}$ 为随机干扰的因素，采用指数形式 $e^{v_{i t}}$ 是为了保证产出大于 0。

对产出函数两边取对数，可得

\ln (Y_{i t}) = \ln f (X_{i t}, β) + \ln (w_{i t}) + v_{i t}

转化为线性计量模型：

y_{i t} = x_{i t}^{'} β + \ln (w_{i t}) + v_{i t}

令 $u_{i t} = - \ln (w_{i t})$ ， $u_{i t} >= 0$ ，

y_{i t} = x_{i t}^{'} β - u_{i t} + v_{i t}

$w_{i t}$ 表示效率水平，由 $u_{i t} = - \ln (w_{i t})$ 可知， $w_{i t} = \exp (- u_{i t})$ 。 $w_{i t}$ 被定义为在相同的投入下，实际产出与完全有效产出的比值。从效率的定义中可以看出，它测量的是一定投入、产出下决策单位之间的“相对效率”而不是“绝对效率”。

实证分析中，通常设定： $v_{i t}$ 服从 $N (0, σ_{v}^{2})$ 为一般意义上的随机干扰项 , $u_{i t}$ 服从单边分布，反映无效率的干扰项。

由于 $u_{i t}$ 无效率项的取值必须大于零，文献中多将其设定成单边分布。典型的分布形式包括: 指数分布、半正态分布、截断型半正态分布。截断型半正态分布可以很方便地研究影响非效率项的因素。

早期面板随机边界模型分为效率不随时间变化的模型和效率可随时间变化两大类。

无效率成分不随时间改变的模型如下:

y_{i t} = x_{i t}^{'} β + ε_{i t} = x_{i t}^{'} β - u_{i} + v_{i t}

其中： $v_{i t}$ 为一般意义上的随机干扰项，服从独立的正态分布 $N (0, σ_{v}^{2})$ 。 $u_{i}$ 服从半正态分布，即 $u_{i} \sim N^{+} (m u, σ_{u}^{2})$ 。

无效率成分随时间改变的模型如下:

y_{i t} = x_{i t}^{'} β - u_{i t} + v_{i t}

文献中常用的假设: $u_{i t} = γ_{t} u_{i}$ ， $u_{i}$ 服从半正态分布， $γ_{t} = \exp (- η (t - T_{i}))$ 。 $t$ 表示观察年份， $T_{i}$ 表示每个个体的最大时间长度。 $η$ 称为延迟参数，用以衡量非效率项随时间的下降程度。

若 $η$ 小于零则表示随着时间的推移，非效率程度在增强。

真正的固定效应模型，简称 TFE － SFA ( True Fixed Effect SFA)。

早期提出的所谓的面板随机边界模型有很明显的局限性 ——— 在产出边界函数的设定中没有考虑不可观测的个体效应，将公司个体效应归入了非效率项会导致非效率程度估计偏误。为了控制不可观测的个体效应，可以将产出边界函数设定为:

y_{i t} = α_{i} + x_{i t} β + v_{i t}

$α_{i}$ 表示不随时间变化且不可观测的个体效应。具体估计方法边际极大似然估计法（MMLE）和模拟边际极大似然估计法 (MMSLE)。具体命令参考 sftfe 及帮助文档。

在 Stata 执行随机边界模型估计主要使用 frontier 和 xtfrontier 命令，外部命令有: sfcross，sfpanel，sftfe，sfkk，SFA2tier 等。

我们以 Stata 官方自带数据集为例，演示 xtfrontier命令执行随机边界模型估计：

use http://www.stata-press.com/data/r13/xtfrontier1 ,clear
xtfrontier lnwidgets lnmachines lnworkers, ti
est store xtsfa_ti
xtfrontier lnwidgets lnmachines lnworkers, tvd
est store xtsfa_tvd
esttab xtsfa_ti xtsfa_tvd, mtitle(xtsfa_ti xtsfa_tvd) 
* 选项 ti 表示不随时间变化， tvd 表示随时间变化。

估计结果如下：

--------------------------------------------
                      (1)             (2)   
                 xtsfa_ti       xtsfa_tvd   
--------------------------------------------
lnwidgets                                   
lnmachines          0.290***        0.291***
                  (17.69)         (17.69)   

lnworkers           0.294***        0.294***
                  (19.07)         (19.06)   

_cons               3.031***        3.029***
                  (21.03)         (21.09)   
--------------------------------------------
lnsigma2                                    
_cons               1.422***        1.411***
                   (5.32)          (5.26)   
--------------------------------------------
lgtgamma                                    
_cons               1.139**         1.124** 
                   (3.20)          (3.14)   
--------------------------------------------
mu                                          
_cons               1.126           1.111   
                   (1.74)          (1.72)   
--------------------------------------------
eta                                         
_cons                             0.00168   
                                   (0.39)   
--------------------------------------------
N                     948             948   
--------------------------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

$η$ 不显著，表明时间效应不明显。

2.3 范例 3：非线性模型的 MLE 估计

指数分布是典型的非线性分布，常用于连续数据，特别是持续时间数据研究。指数分布的概率密度函数形式如下：

f (y) = λ e^{- λ y}, y > 0, λ > 0

均值为 1 $/ λ$ ，方差为 1 $/ λ^{2}$ ，进一步将解释变量加入模型：

λ = \exp (x^{'} β)

非线性指数模型可以表示如下：

y = \exp (- x^{'} β) + u

我们通过模拟数据方式演示如何用 MLE 估计非线性模型。

y | x \sim exponential [λ]

λ = \exp (β_{1} + β_{2} x)

其中： $x \sim N [1, 1^{2}]$ ， $(β_{1}, β_{2}) = (2, - 1)$ ，产生 10,000 笔数据，代码如下:

scalar a = 2
scalar b = -1
scalar mux = 1  
scalar sigx = 1 
set obs 10000
set seed 2003
gen x = mux + sigx*invnorm(uniform()) 
gen lamda = exp(a + b*x)
gen Ey = 1/lamda
* To generate exponential with mean mu=Ey use
 *   Integral 0 to a of (1/mu)exp(-x/mu) dx   by change of variables
 * = Integral 0 to a/mu of exp(-t)dt
 * = incomplete gamma function P(0,a/mu) in the terminology of Stata
gen y = Ey*invgammap(1,uniform())
outfile y x using mma05data.dta, replace

ado 文件基本设定编写如下：

cap program drop mleexp0
program define mleexp0
	version 8.0
	args lnf theta      /* Must use lnf while could use name other than theta */
	quietly replace `lnf' = `theta' - $ML_y1*exp(`theta')
end

MLE 估计命令如下：

sysuse mma05data.dta, clear
ml model lf mleexp0 (y = x)
ml search
ml maximize
estimates store rmle
*Obtain robust standard errors 
ml model lf mleexp0 (y = x), robust
ml search
ml maximize
estimates store rmlerobust
local mm "rmle rmlerobust"
esttab `mm',mtitle(`mm') compress nogap  b(%10.4f) se(%10.4f) t stats(N ll) keep(_cons x)