Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:曹昊煜 (兰州大学)
邮箱:caohy19@lzu.edu.cn
目录
计量经济学中一类重要的估计量称为 “极值估计量”,该估计量通过最大化某个目标函数获得,即:
极值估计量又可以分为 M 估计量和 GMM 估计量。最大似然估计属于 M 估计量的一种,其目标函数的一般形式为样本均值:
最大似然估计主要是将目标函数中的
当使用参数的假象值替代真实参数时,联合密度称为似然函数。真实参数
因此
在一定的正则假定下,ML 估计量为一致渐进正态估计量。
在一般的研究问题中,数据
对等式两边取对数,并使用假想的参数替代真实参数得到对数似然函数:
若条件概率密度和解释变量的概率密度之间不存在函数相关,忽略等式右侧的第二项得到简化的对数自然函数。对应的 M 估计量的形式为:
因此该方法称为条件极大似然估计。在一般的应用中,条件极大似然估计更为普遍,其原因在于我们对联合分布中的参数建模后,被解释变量的分布会依赖于解释变量的变化。
首先使用一个简单的例子来说明 MLE,该示例摘自 B 站 Up 主「小崔说数」的讲解内容。假设城市居民只有两种出行方式,分别是步行 (
1 | 0 | 0 | 1 | 1 | |
---|---|---|---|---|---|
在该模型中,相当于估计参数
对等式两边取对数,则有:
求取一阶条件:
估计结果
求解最大似然估计 (极大似然估计) 的结果可以使用一阶条件,或称为似然方程,但在某些情形下,并不存在参数的显式解,因此需要使用数值方法求解参数。一种简单的,便于考察 MLE 整个过程的方式是网格搜索,当真实参数范围已知时,可以通过搜索的方式估计参数。其过程如下:
最大似然估计的一个简单例子是估计正态分布的均值和标准差。由正态分布的概率密度函数,平均对数似然函数为:
假定数据生成过程服从均值为 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 会自动实现数值计算。
在 Stata 中,可以使用 ml
命令灵活地实现各类最大似然估计。关于 ml
命令的使用和设定方式可以参考帮助文档 help ml
或连享会推文「Stata:数值求解极大值及 MLE 示例」。本文将通过五个简单的范例,介绍如何在 Stata 中实现基本的 MLE 方法。
第一个例子仍然是对正态分布中均值和标准差的估计。此时我们不再使用网格搜索,而是使用其他数值方法,概率密度函数仍然服从正态分布:
在 Stata 中使用 ml
命令,只需要定义特定观测值下的概率密度,而不是整体的对数似然函数,例如观测值
. 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,与设定的数据生成过程基本相同。
对于均值和标准差的估计属于无条件最大似然估计,因为我们并没有对分布中的参数进行建模,也相当于对截距项回归。接下来考察普通的线性模型:
假设随机干扰项
使用 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 的估计结果完全相同,仅仅在标准误上存在细微的差异,如果考察二者的一阶条件,则可以发现二者面临的优化问题是相同的。
接下来考察几个非线性模型的例子。当被解释变量取值为 0-1 时,需要构建定性反应模型,其基本的概率分布为:
对参数
当假定
其中,
在 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 模型,直接编写的最大似然估计结果与官方命令提供的结果基本没有差异。
有时候,我们获取的样本并非从总体中随机抽样得到,而是从总体的子样本中获取的,例如规模以上的企业调查。假设随机变量
假设无断尾
从而对数似然函数为:
在 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
------------------------------------------------------------------------------
截取由 Tobin (1958) 首次引入经济学,与断尾分布不同,此时样本仍然是从总体中随机抽取的,而当被解释变量的取值小于某个特定的值时,全部被归并到该取值上。截取模型又称为归并模型或 Tobit 模型。其形式为潜变量模型:
可以将
相应地,条件对数似然函数为:
使用 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
模型的估计结果完全相同。
本文主要介绍了最大似然估计和条件最大似然估计的理论过程,并使用网格搜索对 MLE 中的数值计算进行了简单的示例,最后使用五个范例说明了在 Stata 中如何使用 ml
命令进行最大似然估计。
本文的范例主要涉及线性和非线性的单方程模型。实际上多方程模型同样可以使用 MLE,具体理论包括似无相关回归、完全信息极大似然估计和有限信息极大似然估计等。此外,ml
命令还可以进行更加复杂的设定,我们将在之后的推文中介绍。
Note:产生如下推文列表的 Stata 命令为:
lianxh mle probit tobit, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh