Stata：最大似然估计入门教程(MLE)-ml

发布时间：2022-08-18 阅读 2199

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

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

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

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

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：曹昊煜 (兰州大学)
邮箱：caohy19@lzu.edu.cn

1. 最大似然估计与条件最大似然估计

1.1 M 估计量

计量经济学中一类重要的估计量称为 “极值估计量”，该估计量通过最大化某个目标函数获得，即：

\hat{θ} = \underset{θ \in Θ}{\arg max} Q_{n} (θ)

极值估计量又可以分为 M 估计量和 GMM 估计量。最大似然估计属于 M 估计量的一种，其目标函数的一般形式为样本均值：

Q_{n} (θ) = \frac{1}{n} \sum_{t = 1}^{n} m (w_{t}; θ)

1.2 最大似然估计

最大似然估计主要是将目标函数中的 $m (w_{t}; θ)$ 替换为密度函数的对数形式。令 ${w_{t}}$ 为独立同分布的序列，且概率密度为 $f (w_{t}; θ_{0})$ ，则数据 ${w_{1}, w_{2}, \dots, w_{n}}$ 的联合密度为：

f (w_{1}, w_{2}, \dots, w_{n}) = \prod_{t = 1}^{n} f (w_{t}; θ_{0})

当使用参数的假象值替代真实参数时，联合密度称为似然函数。真实参数 $θ_{0}$ 的 ML 估计是选择相应的 $θ$ ，使得似然函数最大化。由于对数变换是单调的，最大化似然函数等价于对数似然函数最大化：

\log f (w_{1}, w_{2}, \dots, w_{n}; θ) = \log [\prod_{t = 1}^{n} f (w_{t}; θ)] = \sum_{t = 1}^{n} \log f (w_{t}; θ)

因此 $θ_{0}$ 的 ML 估计量为 M 估计量：

Q_{n} (θ) = \frac{1}{n} \sum_{t = 1}^{n} m (w_{t}; θ) = \frac{1}{n} \sum_{t = 1}^{n} \log f (w_{t}; θ)

在一定的正则假定下，ML 估计量为一致渐进正态估计量。

1.3 条件最大似然估计

在一般的研究问题中，数据 $w_{t} = (y_{t}, x_{t}^{'})^{'}$ 可以区分为解释变量和被解释变量，而参数估计的主要目标是考察解释变量 $x_{t}$ 对被解释变量 $y_{t}$ 条件分布的影响。根据条件概率密度公式，数据的联合密度可以表示为：

f (y_{t}, x_{t}; θ_{0}, Ψ_{0}) = f (y_{t} | x_{t}; θ_{0}) f (x_{t}; Ψ_{0})

对等式两边取对数，并使用假想的参数替代真实参数得到对数似然函数：

\frac{1}{n} \sum_{t = 1}^{n} \log f (w_{t}; θ, Ψ) = \frac{1}{n} \sum_{t = 1}^{n} \log f (y_{t} | x_{t}; θ) + \frac{1}{n} \sum_{t = 1}^{n} \log f (x_{t}; Ψ)

若条件概率密度和解释变量的概率密度之间不存在函数相关，忽略等式右侧的第二项得到简化的对数自然函数。对应的 M 估计量的形式为：

Q_{n} (θ) = \frac{1}{n} \sum_{t = 1}^{n} \log f (y_{t} | x_{t}; θ)

因此该方法称为条件极大似然估计。在一般的应用中，条件极大似然估计更为普遍，其原因在于我们对联合分布中的参数建模后，被解释变量的分布会依赖于解释变量的变化。

2. 最大似然估计两种方法

2.1 最大似然估计：极简示例

首先使用一个简单的例子来说明 MLE，该示例摘自 B 站 Up 主「小崔说数」的讲解内容。假设城市居民只有两种出行方式，分别是步行 ( $y_{i} = 1$ ) 或开车 ( $y_{i} = 0$ )。二者发生的概率分别为 $θ$ 和 $1 - θ$ 。随机询问五个人，结果为 ${1, 0, 0, 1, 1}$ ，从而可以得到：

$y_{i}$	1	0	0	1	1
$p_{i}$	$θ$	$1 - θ$	$1 - θ$	$θ$	$θ$

在该模型中，相当于估计参数 $P (y_{i} = 1) = θ$ 。由于样本是随机抽取的，因此联合概率密度等于所有样本概率之积，定义似然函数：

L (θ) = θ \cdot (1 - θ) \cdot (1 - θ) \cdot θ \cdot θ = (1 - θ)^{2} θ^{3}

对等式两边取对数，则有：

\hat{θ} = \underset{θ}{\arg max} \log L (θ) = 3 \log θ + 2 \log (1 - θ)

求取一阶条件：

\frac{\partial \log L (θ)}{\partial θ} = \frac{3}{θ} + \frac{2}{1 - θ} = 0 \Rightarrow \hat{θ} = \frac{3}{5}

估计结果 $\hat{θ}$ 与 $y_{i}$ 的样本均值相同，因为在二值观测中，期望与取值为 1 的概率相同。通过上面这个简单的例子，可以看出 MLE 的基本步骤：

第一步：确定 $y_{i}$ 的概率分布；
第二步：写出样本中所有观测值的联合分布 $L (θ)$ ；
第三步：极大化 $\log L (θ)$ 从而得到 $θ$ 的估计。

2.2 最大似然估计：网格搜索

求解最大似然估计 (极大似然估计) 的结果可以使用一阶条件，或称为似然方程，但在某些情形下，并不存在参数的显式解，因此需要使用数值方法求解参数。一种简单的，便于考察 MLE 整个过程的方式是网格搜索，当真实参数范围已知时，可以通过搜索的方式估计参数。其过程如下：

以某个特定的步长在参数的取值范围内选择不同的参数值；
计算特定参数取值下的似然函数值；
选取最大化似然函数取值对应的参数作为真实参数的估计。

最大似然估计的一个简单例子是估计正态分布的均值和标准差。由正态分布的概率密度函数，平均对数似然函数为：

\frac{1}{n} \sum_{t = 1}^{n} \log f (w_{t}; μ, σ^{2}) = - \frac{1}{2} \log (2 π) - \frac{1}{2} \log (σ^{2}) - \frac{1}{n} \sum_{t = 1}^{n} \frac{(w_{t} - μ)^{2}}{2 σ^{2}}

假定数据生成过程服从均值为 1，标准差为 0.8 的正态分布，则网格搜索的基本思路是同时选择不同的 $μ$ 和 $σ$ ，并计算对数似然值并比较其大小。过程如下：

. clear
. set obs 10000
. set seed 1010110 // 设定种子值
. gen y = rnormal(1,0.8) in 1/1000 // 生成随机样本 y ~ N(1,0.8)
. gen mu = .
. gen sigma = .
. gen logL = . // 生成用于保存均值、标准差和对数似然值的
. local n = 1

. * 开始网格搜索
. forvalue mu = 0.01(0.01)1{
  2.   forvalues sigma = 0.01(0.01)1{    
  3.     gen term1 = (y-`mu')^2/(2*`sigma'^2)
  4.     egen term2 = sum(term1)     
  5.     qui sum term2
  6.     local last = r(max)
  7.     local N = 1000    
  8.     replace logL = -0.5*`N'*log(2*_pi) - 0.5*`N'*log(`sigma'^2) - `last' in `n'  // 保存对数似然函数值
  9.     replace mu = `mu' in `n'
 10.     replace sigma = `sigma' in `n'
 11.     drop term1 term2
 12.     local n = `n'+1
 13.   }  
 14. }

. * 搜索结果
. gsort - logL
. list mu sigma logL in 1/10

     +-------------------------+
     |  mu   sigma        logL |
     |-------------------------|
  1. | .99     .81   -1208.978 |
  2. | .99     .82   -1209.109 |
  3. | .99      .8   -1209.152 |
  4. | .99     .83   -1209.527 |
  5. | .99     .79   -1209.652 |
     |-------------------------|
  6. | .98     .81   -1209.714 |
  7. | .98     .82   -1209.827 |
  8. | .98      .8   -1209.907 |
  9. | .99     .84   -1210.215 |
 10. | .98     .83   -1210.228 |
     +-------------------------+

可以看出，均值和标准差的搜索结果与真实的参数非常接近。当然，在具体的 MLE 中，我们并不知道真实参数的取值范围，因此需要使用牛顿-拉夫森等其他数值计算方法，在以下的例子中 Stata 会自动实现数值计算。

3. Stata 范例

在 Stata 中，可以使用 ml 命令灵活地实现各类最大似然估计。关于 ml 命令的使用和设定方式可以参考帮助文档 help ml 或连享会推文「Stata：数值求解极大值及 MLE 示例」。本文将通过五个简单的范例，介绍如何在 Stata 中实现基本的 MLE 方法。

3.1 范例一：估计正态分布的参数

第一个例子仍然是对正态分布中均值和标准差的估计。此时我们不再使用网格搜索，而是使用其他数值方法，概率密度函数仍然服从正态分布：

\frac{1}{n} \sum_{t = 1}^{n} \log f (w_{t}; μ, σ^{2}) = - \frac{1}{2} \log (2 π) - \frac{1}{2} \log (σ^{2}) - \frac{1}{n} \sum_{t = 1}^{n} \frac{(w_{t} - μ)^{2}}{2 σ^{2}}

在 Stata 中使用 ml 命令，只需要定义特定观测值下的概率密度，而不是整体的对数似然函数，例如观测值 $t$ 的概率密度是均值为 $μ$ ，标准差为 $σ$ 的正态分布概率密度函数。

. clear
. set obs 1000
. set seed 1010110
. gen y = rnormal(1,0.8)       // 生成随机样本 y ~ N(1,0.8)
. cap program drop mynormal_lf  
. program define mynormal_lf   // 定义对数似然函数
  1.   version 16
  2.   args lnf mu sigma
  3.   qui replace `lnf' = ln(normalden($ML_y1,`mu',`sigma')) 
  4. end
. ml model lf mynormal_lf (y = ) (sigma: ) // 模型形式
. ml maximize                              // 参数求解

Iteration 0:   log likelihood = -1247.1045  
Iteration 1:   log likelihood = -1207.7554  
Iteration 2:   log likelihood = -1207.5499  
Iteration 3:   log likelihood = -1207.5496  
Iteration 4:   log likelihood = -1207.5496  
                                                         Number of obs = 1,000
                                                         Wald chi2(0)  =     .
Log likelihood = -1207.5496                              Prob > chi2   =     .
------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
eq1          |
       _cons |      1.033      0.026    40.37   0.000        0.983       1.083
-------------+----------------------------------------------------------------
sigma        |
       _cons |      0.809      0.018    44.72   0.000        0.774       0.845
------------------------------------------------------------------------------

关于似然函数的定义需要注意以下几点：

cap program drop：用于删去可能存在的重复程序；
progname_lf：是程序命名的一般形式，即自己确定的程序名称加上下划线和 lf 后缀；
version 16：声明 Stata 版本；
args：后面紧跟对数似然函数值 lnf 和参数名称；
$ML_y1：代表被解释变量。

估计结果中显示了 MLE 的迭代过程，最大化的对数似然值为 -1207.5496，估计出的均值为 1.033，标准差为 0.809，与设定的数据生成过程基本相同。

3.2 范例二：线性模型

对于均值和标准差的估计属于无条件最大似然估计，因为我们并没有对分布中的参数进行建模，也相当于对截距项回归。接下来考察普通的线性模型：

y_{t} = β_{0} + β_{1} x_{1 t} + β_{2} x_{2 t} + \dots + β_{k} x_{k t} + ε_{t}

假设随机干扰项 $ε_{t} \sim N (0, σ_{0}^{2})$ 。从而 $y_{t} \sim N (x_{t}^{'} β, σ_{0}^{2})$ ，即对正态分布的均值进行了线性建模。相应的条件对数似然函数为：

\log f (y_{t} | x_{t}; β) = \log {\frac{1}{\sqrt{2 π σ^{2}}} \exp [\frac{- (2 - x_{t}^{'} β)^{2}}{2 σ^{2}}]}

使用 Stata 进行线性模型的 ML 估计，并与 regress 命令做对比：

. sysuse auto.dta, clear
. cap program drop mylinear_lf
. program define mylinear_lf
  1.   version 16
  2.   args lnf xb sigma
  3.   qui replace `lnf' = ln(normalden($ML_y1,`xb',`sigma'))
  4. end
. ml model lf mylinear_lf (price = foreign trunk weight length) (sigma: )
. ml maximize, nolog

                                                        Number of obs =     74
                                                        Wald chi2(4)  =  90.07
Log likelihood = -666.25274                             Prob > chi2   = 0.0000
------------------------------------------------------------------------------
       price | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
eq1          |
     foreign |   3580.051    624.022     5.74   0.000     2356.990    4803.111
       trunk |    -10.282     78.824    -0.13   0.896     -164.774     144.210
      weight |      5.769      0.934     6.18   0.000        3.938       7.600
      length |    -89.653     34.534    -2.60   0.009     -157.339     -21.968
       _cons |   4672.931   3852.845     1.21   0.225    -2878.506   12224.368
-------------+----------------------------------------------------------------
sigma        |
       _cons |   1967.417    161.721    12.17   0.000     1650.451    2284.384
------------------------------------------------------------------------------

. reg price foreign trunk weight length

      Source |       SS           df       MS      Number of obs   =        74
-------------+----------------------------------   F(4, 69)        =     21.00
       Model |   348631330         4  87157832.5   Prob > F        =    0.0000
    Residual |   286434066        69  4151218.35   R-squared       =    0.5490
-------------+----------------------------------   Adj R-squared   =    0.5228
       Total |   635065396        73  8699525.97   Root MSE        =    2037.5
------------------------------------------------------------------------------
       price | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     foreign |   3580.051    646.236     5.54   0.000     2290.845    4869.256
       trunk |    -10.282     81.630    -0.13   0.900     -173.129     152.565
      weight |      5.769      0.967     5.96   0.000        3.839       7.699
      length |    -89.653     35.763    -2.51   0.015     -160.999     -18.307
       _cons |   4672.931   3989.999     1.17   0.246    -3286.899   12632.762
------------------------------------------------------------------------------

可见 MLE 与 OLS 的估计结果完全相同，仅仅在标准误上存在细微的差异，如果考察二者的一阶条件，则可以发现二者面临的优化问题是相同的。

3.3 范例三：定性反应模型

接下来考察几个非线性模型的例子。当被解释变量取值为 0-1 时，需要构建定性反应模型，其基本的概率分布为：

f (y_{t}) = π^{y_{t}} (1 - π)^{1 - y_{t}}

对参数 $π$ 进行建模，由于其代表的是 $y_{t}$ 取 1 的概率，因此其函数形式可以定义为正态分布或 Logistic 分布。当假定 $π$ 服从正态分布时，该模型称为 Probit 模型：

{\begin{aligned} f (y_{t} = 1 | x_{t}; θ_{0}) = Φ (x_{t}^{'} θ_{0}) \\ f (y_{t} = 0 | x_{t}; θ_{0}) = 1 - Φ (x_{t}^{'} θ_{0}) \end{aligned}

当假定 $π$ 服从 Logistic 分布时，该模型称为 Logit 模型：

{\begin{aligned} f (y_{t} = 1 | x_{t}; θ_{0}) = Λ (x_{t}^{'} θ_{0}) \\ f (y_{t} = 0 | x_{t}; θ_{0}) = 1 - Λ (x_{t}^{'} θ_{0}) \end{aligned}

其中，

Λ (x_{t}^{'} θ_{0}) = \frac{\exp (x_{t}^{'} θ_{0})}{1 + \exp (x_{t}^{'} θ_{0})}

在 Stata 中，normal 和 invlogit 分别表示正态分布累积分布函数和拟 Logit 变换，probit 和 logit 命令分别用于估计上述模型。

. * probit 模型
. sysuse auto.dta, clear
. cap program drop myprobit_lf
. program define myprobit_lf
  1.   version 16
  2.   args lnf xb
  3.   qui replace `lnf' = ln(normal(`xb')) if $ML_y1 == 1
  4.   qui replace `lnf' = ln(1 - normal(`xb')) if $ML_y1 == 0
  5. end
. ml model lf myprobit_lf (foreign = price trunk weight length)
. ml maximize, nolog

                                                        Number of obs =     74
                                                        Wald chi2(4)  =  15.11
Log likelihood = -17.72839                              Prob > chi2   = 0.0045
------------------------------------------------------------------------------
     foreign | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       price |      0.001      0.000     3.32   0.001        0.000       0.001
       trunk |      0.010      0.083     0.12   0.901       -0.153       0.174
      weight |     -0.004      0.002    -2.71   0.007       -0.007      -0.001
      length |      0.033      0.045     0.74   0.460       -0.055       0.121
       _cons |      1.201      5.070     0.24   0.813       -8.735      11.138
------------------------------------------------------------------------------

. probit foreign price trunk weight length, nolog

Probit regression                                       Number of obs =     74
                                                        LR chi2(4)    =  54.61
                                                        Prob > chi2   = 0.0000
Log likelihood = -17.72839                              Pseudo R2     = 0.6063
------------------------------------------------------------------------------
     foreign | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       price |      0.001      0.000     3.32   0.001        0.000       0.001
       trunk |      0.010      0.083     0.12   0.901       -0.153       0.174
      weight |     -0.004      0.002    -2.71   0.007       -0.007      -0.001
      length |      0.033      0.045     0.74   0.460       -0.055       0.121
       _cons |      1.201      5.070     0.24   0.813       -8.735      11.138
------------------------------------------------------------------------------

. * logit 模型 
. cap program drop mylogit_lf
. program define mylogit_lf
  1.   version 16
  2.   args lnf xb
  3.   qui replace `lnf' = ln(invlogit(`xb')) if $ML_y1 == 1
  4.   qui replace `lnf' = ln(1 - invlogit(`xb')) if $ML_y1 == 0
  5.     
.     * qui replace `lnf' = ln(exp(`xb')/(1+exp(`xb'))) if $ML_y1 == 1
.     * qui replace `lnf' = ln(1 - exp(`xb')/(1+exp(`xb'))) if $ML_y1 == 0
. end
. ml model lf mylogit_lf (foreign = price trunk weight length)
. ml maximize, nolog

                                                        Number of obs =     74
                                                        Wald chi2(4)  =  12.61
Log likelihood = -17.786494                             Prob > chi2   = 0.0134
------------------------------------------------------------------------------
     foreign | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       price |      0.001      0.000     3.14   0.002        0.000       0.002
       trunk |      0.009      0.142     0.06   0.949       -0.269       0.287
      weight |     -0.007      0.003    -2.55   0.011       -0.013      -0.002
      length |      0.050      0.081     0.61   0.539       -0.109       0.209
       _cons |      3.388      9.308     0.36   0.716      -14.855      21.631
------------------------------------------------------------------------------

. logit foreign price trunk weight length, nolog

Logistic regression                                     Number of obs =     74
                                                        LR chi2(4)    =  54.49
                                                        Prob > chi2   = 0.0000
Log likelihood = -17.786494                             Pseudo R2     = 0.6050
------------------------------------------------------------------------------
     foreign | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       price |      0.001      0.000     3.14   0.002        0.000       0.002
       trunk |      0.009      0.142     0.06   0.949       -0.269       0.287
      weight |     -0.007      0.003    -2.55   0.011       -0.013      -0.002
      length |      0.050      0.081     0.61   0.539       -0.109       0.209
       _cons |      3.388      9.308     0.36   0.716      -14.855      21.631
------------------------------------------------------------------------------

从估计结果中，无论是 Probit 还是 Logit 模型，直接编写的最大似然估计结果与官方命令提供的结果基本没有差异。

3.4 范例四：断尾回归

有时候，我们获取的样本并非从总体中随机抽样得到，而是从总体的子样本中获取的，例如规模以上的企业调查。假设随机变量 $y$ 的分布为 $f (y)$ ，则 $y > c$ 的断尾分布密度函数为：

f (y | y > c) = \frac{f (y)}{Prob (y > c)}

假设无断尾 $y_{t}$ 的分布为 $N (x_{t}^{'} β_{0}, σ_{0}^{2})$ ，则在 $c$ 处断尾的分布为：

f (y_{t} | y_{t} > c, x_{t}; β_{0}, σ_{0}^{2}) = \frac{\frac{1}{σ_{0}} ϕ (\frac{y_{t} - x_{t}^{'} β_{0}}{σ_{0}^{2}})}{1 - Φ (\frac{c - x_{t}^{'} β_{0}}{σ_{0}})}

从而对数似然函数为：

\log f (y_{t} | y_{t} > c, x_{t}; β_{0}, σ_{0}^{2}) = \log {\frac{\frac{1}{σ_{0}} ϕ (\frac{y_{t} - x_{t}^{'} β_{0}}{σ_{0}^{2}})}{1 - Φ (\frac{c - x_{t}^{'} β_{0}}{σ_{0}})}}

在 Stata 中，truncreg 命令可以用于估计断尾回归，因此可以使用该命令验证自己编写的断尾回归代码。

. lxhuse laborsub.dta, clear
. keep if whrs > 0      // 假设在 0 处断尾
. cap program drop mytrunc_lf
. program define mytrunc_lf
  1.   version 16
  2.   args lnf xb sigma
  3.   qui replace `lnf' = ln(((1/`sigma')*normalden(($ML_y1-`xb')/`sigma'))/(1-normal((0-`xb')/`sigma')))    
  4. end
. ml model lf mytrunc_lf (whrs = kl6 k618 wa we) (sigma: )
. ml maximize, nolog

                                                        Number of obs =    150
                                                        Wald chi2(4)  =  10.05
Log likelihood = -1200.9157                             Prob > chi2   = 0.0395
------------------------------------------------------------------------------
        whrs | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
eq1          |
         kl6 |   -803.004    321.361    -2.50   0.012    -1432.860    -173.148
        k618 |   -172.875     88.729    -1.95   0.051     -346.780       1.031
          wa |     -8.821     14.368    -0.61   0.539      -36.983      19.341
          we |     16.529     46.504     0.36   0.722      -74.617     107.674
       _cons |   1586.260    912.354     1.74   0.082     -201.922    3374.441
-------------+----------------------------------------------------------------
sigma        |
       _cons |    983.726     94.443    10.42   0.000      798.621    1168.830
------------------------------------------------------------------------------

. truncreg whrs kl6 k618 wa we, ll(0) nolog

Truncated regression
Limit: Lower =    0                             Number of obs     =        150
       Upper = +inf                             Wald chi2(4)      =      10.05
Log likelihood = -1200.9157                     Prob > chi2       =     0.0395
------------------------------------------------------------------------------
        whrs | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         kl6 |   -803.004    321.361    -2.50   0.012    -1432.861    -173.147
        k618 |   -172.875     88.729    -1.95   0.051     -346.781       1.031
          wa |     -8.821     14.368    -0.61   0.539      -36.983      19.341
          we |     16.529     46.504     0.36   0.722      -74.617     107.674
       _cons |   1586.260    912.355     1.74   0.082     -201.923    3374.442
-------------+----------------------------------------------------------------
      /sigma |    983.726     94.443    10.42   0.000      798.621    1168.831
------------------------------------------------------------------------------

3.5 范例五：截取回归

截取由 Tobin (1958) 首次引入经济学，与断尾分布不同，此时样本仍然是从总体中随机抽取的，而当被解释变量的取值小于某个特定的值时，全部被归并到该取值上。截取模型又称为归并模型或 Tobit 模型。其形式为潜变量模型：

y_{t}^{*} = x_{t}^{'} β_{0} + ε_{t}, ε_{t} \sim N (0, σ_{0}^{2})

y_{t} = {\begin{aligned} y_{t}^{*} & i f y_{t}^{*} > c \\ c & i f y_{t}^{*} < c \end{aligned}

可以将 $y_{t}$ 的分布视为两个部分， $y_{t} > c$ 的部分是正常的观测值，其分布与潜变量 $y_{t}^{*}$ 相同，而 $y_{t} = c$ 的概率与 $y_{t}^{*} \leq c$ 的概率相同。因此 Tobit 模型的概率密度由两个部分组成：

f (y_{t} | x_{t}; β_{0}, σ_{0}) = {[\frac{1}{σ_{0}} ϕ (\frac{y_{t} - x_{t}^{'} β_{0}}{σ_{0}^{2}})]}^{1 (y_{t} > c)} \times {[Φ (\frac{c - x_{t}^{'} β_{0}}{σ_{0}})]}^{1 (y_{t} = c)}

相应地，条件对数似然函数为：

\log f (y_{t} | x_{t}; β_{0}, σ_{0}) = \log {{[\frac{1}{σ_{0}} ϕ (\frac{y_{t} - x_{t}^{'} β_{0}}{σ_{0}^{2}})]}^{1 (y_{t} > c)} \times {[Φ (\frac{c - x_{t}^{'} β_{0}}{σ_{0}})]}^{1 (y_{t} = c)}}

使用 tobit 命令来验证自己编写的 MLE 程序：

. lxhuse womenwk.dta, clear
. cap program drop mytobit_lf
. program define mytobit_lf
  1.   version 16
  2.   args lnf xb sigma
  3.   qui replace `lnf' = ln((1/`sigma')*normalden(($ML_y1-`xb')/`sigma')) if $ML_y1 != 0
  4.   qui replace `lnf' = ln(normal((0-`xb')/`sigma')) if $ML_y1 == 0
  5. end
. ml model lf mytobit_lf (lwf = age married children education) (sigma: )
. ml maximize, nolog

                                                        Number of obs =  2,000
                                                        Wald chi2(4)  = 472.59
Log likelihood = -3349.9685                             Prob > chi2   = 0.0000
------------------------------------------------------------------------------
         lwf | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
eq1          |
         age |      0.052      0.006     9.08   0.000        0.041       0.063
     married |      0.484      0.104     4.68   0.000        0.281       0.687
    children |      0.486      0.032    15.33   0.000        0.424       0.548
   education |      0.115      0.015     7.62   0.000        0.085       0.145
       _cons |     -2.808      0.263   -10.67   0.000       -3.324      -2.292
-------------+----------------------------------------------------------------
sigma        |
       _cons |      1.873      0.040    46.80   0.000        1.794       1.951
------------------------------------------------------------------------------

. tobit lwf age married children education, ll(0) nolog

Tobit regression                                    Number of obs     =  2,000
                                                           Uncensored =  1,343
Limits: Lower =    0                                    Left-censored =    657
        Upper = +inf                                   Right-censored =      0
                                                    LR chi2(4)        = 461.85
                                                    Prob > chi2       = 0.0000
Log likelihood = -3349.9685                         Pseudo R2         = 0.0645
------------------------------------------------------------------------------
         lwf | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |      0.052      0.006     9.08   0.000        0.041       0.063
     married |      0.484      0.104     4.68   0.000        0.281       0.687
    children |      0.486      0.032    15.33   0.000        0.424       0.548
   education |      0.115      0.015     7.62   0.000        0.085       0.145
       _cons |     -2.808      0.263   -10.67   0.000       -3.324      -2.291
-------------+----------------------------------------------------------------
   var(e.lwf)|      3.507      0.150                         3.225       3.814
------------------------------------------------------------------------------

MLE 的估计结果与 Tobit 模型的估计结果完全相同。

4. 结束语

本文主要介绍了最大似然估计和条件最大似然估计的理论过程，并使用网格搜索对 MLE 中的数值计算进行了简单的示例，最后使用五个范例说明了在 Stata 中如何使用 ml 命令进行最大似然估计。

本文的范例主要涉及线性和非线性的单方程模型。实际上多方程模型同样可以使用 MLE，具体理论包括似无相关回归、完全信息极大似然估计和有限信息极大似然估计等。此外，ml 命令还可以进行更加复杂的设定，我们将在之后的推文中介绍。

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh mle probit tobit, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]