Stata 程序：数值求解极值的基本思想
温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。
作者： 徐云娇 (厦门大学)
邮箱： jilyo@stu.xmu.edu.cn
目录
1. while 语句基本介绍
极值求解背后的思想是当两次搜索得到的数值之差大于收敛判据时，则继续搜索；当差值小于收敛判据时，则结束循环，输出最终搜索结果。
以上思想的实现过程需要借助 Stata 中的 while 循环语句，while 后面会跟一个进行逻辑判断的条件，该条件为真时，该命令会一直执行括号里面的内容，直到条件不被满足时，该循环终止。
以下这个简单的 while 循环可以帮助大家理解其工作原理：
local j = 0 //定义暂元 j，初值为 0
while `j' < 5{ //当 j 小于 5 时，执行以下命令
dis `j' //在结果窗口列示 j
local j = `j'+1 //令新的 j 等于原始值加上 1
}
*-或
local j = 0
while `j' < 5{
dis `j++' //效果等同于 dis `j' + local j = `j'+1
}
在 Stata 中运行以上代码，结果窗口便会列示 5 行，分别显示 0、1、2、3、4，可见一旦 j = 5 时，循环便会停止。
值得注意的是，除了程序中的 j++ 外，还存在 ++j 这种语法，后者表示的是先把 j 加上 1 再列示 j 的当前值，所以当把以上程序中的 j++ 修改为 ++j 后，结果窗口将列示 1、2、3、4、5，应注意对二者的区分。
2. 数值求解极大值
在理解了 while 语句的使用规则之后，我们便可以利用 Stata 采用数值方法求取函数极大值：
*-示例: 采用数值法求取函数的极大值
twoway function y = -0.2*exp(x) + ln(x^2) + 7, ///
range(0 4) lw(*2)
假如我们想求解以上图示函数的极大值，那么在数值搜索的过程中涉及几个要点，分别是：初始值，步长，收敛判据。
求解的基本步骤是：
(1) 确定一个初始值
Stata 程序如下：
local delta = 0.05 // 步长
local x = 1 // x 的初始值
local j = 0 // 计数器：记录迭代次数
local e = 1 // y1-y0
local e0 = 0.01 // 收敛判据
while `e' > `e0'{
`trace'
local y0 = -0.2*exp(`x') + ln(`x'^2) + 7
local x = `x' + `delta'
local y1 = -0.2*exp(`x') + ln(`x'^2) + 7
local e = abs(`y1' - `y0')
dis in g "*" _c
local j = `j' + 1
}
dis "e = " `e'
dis "x = " `x' // x 的解
dis "y = " `y1' // y 的极大值
dis "j = " `j' // 迭代次数
*-图示：
twoway function y = -0.2*exp(x) + ln(x^2) + 7, ///
range(0 4) lw(thick) xline(`x') yline(`y1') ///
text(`=`y1'-0.5' `=`x'+0.8' "(`x'，`y1')")
执行以上程序，经过 14 次搜索，函数差值的绝对值低至 0.0063，此时我们便近似找到了函数的极大值点 1.7，以及其对应的函数极大值 6.966467，用图形表示如下：
在整个分析的过程中，有几个关键要素，分别是：
初始值：越靠近真实极值点的初始值无疑会减少搜索的次数从而让我们更快得到极值。
步长：假如设定一个大步长，那搜索的速度就会更快；但大步长也有风险，在某些情境下可能会直接跨过极值点导致求解出不够精确的极值。
收敛判据：越严格 (即越小) 的收敛判据会使收敛条件更难被满足，可能会出现循环无法停止从而找不到极值的情况。
更换上述程序中有关这三个关键变量的初始设定，进一步测试程序的稳定性：
(1) (更改初始值) 设定 x=2 为初始值，能否收敛？
执行程序之后未能收敛。观察函数图形易知，由于程序只能沿着坐标轴正向搜索，所以当把初始值设定为 x=2 时，两次搜索之间函数差值的绝对值只会越来越大，从而达不到收敛判据，最终未能收敛。
(2) (更改步长) 设定 delta=0.1，结果有何变化？
重新运行程序，此时经过 8 次搜索，函数差值的绝对值低至 0.00082，函数的极大值点为 1.8，以及其对应的函数极大值 y = 6.9656438，略小于步长为 0.05 时得到的极大值。
(3) (更改收敛判据) 设定 e0=0.0001，能否收敛?
执行程序之后未能收敛。在过小的收敛判据设定下，搜索过程中的收敛条件得不到满足，循环无法结束，所以极大值数值求解失败。
接下去进一步完善数值求解的过程，增加反向搜索功能：
*- 程序修改如下：
local h = 0.05 // 步长
local x = 1 // x 的初始值
local j = 0 // 计数器：记录迭代次数
local e = 1 // y1-y0
local e0 = `h'/10 // 收敛判据 (修改为动态数值)
while abs(`e')>`e0'{ // 修改：abs(`e')
local y0 = -0.2*exp(`x') + ln(`x'^2) + 7
local x = `x' + `h'
local y1 = -0.2*exp(`x') + ln(`x'^2) + 7
local e = `y1' - `y0' // 此前 e = abs(`y1'-`y0')
if (`e' < 0){
local h = -`h' // 新增：反向搜索
}
dis in g "*" _c
local j = `j' + 1
}
dis "e = " `e'
dis "x = " `x' // x 的解
dis "y = " `y1' // y 的极大值
dis "j = " `j' // 迭代次数
*-图示：
local x: dis %4.3f `x' // 新增：显示的更美观
local y: dis %4.3f `y1'
twoway function y = -0.2*exp(x) + ln(x^2) + 7, ///
range(0 4) lw(thick) xline(`x') yline(`y1') ///
text(`=`y'+0.5' `=`x'+0.4' "(`x', `y')")
执行以上程序，经过 15 次搜索，两次搜索的函数差值低至 0.00184，此时我们便近似找到了函数的极大值点 1.750，以及其对应的函数极大值 6.968，这比我们之前找到的极大值都略大一些，用图形表示如下：
修改后的程序不但帮助我们找到了更大的极大值，而且在其他方面优势也有很多，可以再次改变一些初始设定从而测试程序的稳定性：
(1) (更改初始值) 设定 x=3 为初始值，能否收敛？
由于新增了反向搜索的功能，所以此时执行程序之后可以收敛。经过 27 次搜索，再次锁定了极大值点 1.75，得到了与之前相同的极大值 6.968。
(2) (更改步长、收敛判据) 设定 delta=0.001，e0=0.0001，能否收敛？
可以收敛，执行以上程序，经过 691 次搜索，两次搜索的函数差值低至 0.0000986，小于收敛判据，此时我们确定了函数的极大值点 1.691，以及其对应的函数极大值 6.966。
需要指出的是，以上例子中只存在一个极值，但现实中还可能存在有多个极值点的情况，这个时候初始值、步长和收敛判据的选择就更加重要了。
3. MLE估计方法的实现
3.1 MLE 简介
用一句话来概括极大似然估计，那就是：利用已知的样本结果，反推最有可能 (最大概率) 导致这样结果的参数值。
记已知的独立同分布样本集为
与其对应的似然函数 (likelihood function) 为：
极大似然估计 (Maximum Likelihood Estimation, MLE) 就是在求解一个最大化问题：
为了更便于计算，通常会对似然函数进行对数转换：
3.2 Stata 实例
利用 Stata 实现 MLE 估计更加详细的介绍参见「极大似然估计 (MLE) 及 Stata 实现」，这里对 help ml 中有关最优化求解的选项进行介绍：
以简单的正态分布函数为例，并结合汽车价格数据，一个标准的 MLE 求解过程有两步：
第一步：将似然函数的基本设定编写为 ado 文件，并在 dofile 中选中相应代码，右键选择 "Tools" → "Execute selection quietly (run)"；
cap program drop mymean_lf
program define mymean_lf
version 9.1 //设定 Stata 9.0版本，以便更高版本能够兼容
args lnf mu sigma //输入项：似然函数，参数1，参数2，......
quietly replace `lnf' = ln(normalden($ML_y1, `mu', `sigma'))
end
/*$ML_y1 默认写法，表示被解释变量*/
第二步：调入汽车价格的数据，利用 ml 命令进行极大似然估计。
sysuse auto, clear
ml model lf mymean_lf (price = mpg wei len) (sigma: )
/* lf 是最常用的架构，通常对样本满足独立同分布的条件时均可使用；
假定汽车价格 (price) 服从均值为 mu，方差为 sigma 的正态分布，
并将 mu 设定为汽车每公里耗油量 (mpg)、重量 (wei)、长度 (len) 的线性函数*/
ml search
ml maximize
执行以上代码，便可得到估计结果：
Number of obs = 74
Wald chi2(3) = 41.15
Log likelihood = -679.35161 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
price | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
eq1 |
mpg | -86.789 81.643 -1.06 0.288 -246.807 73.228
weight | 4.365 1.135 3.84 0.000 2.139 6.590
length | -104.868 38.633 -2.71 0.007 -180.587 -29.149
_cons | 14542.430 5729.204 2.54 0.011 3313.396 25771.464
-------------+----------------------------------------------------------------
sigma |
_cons | 2348.394 193.036 12.17 0.000 1970.050 2726.738
------------------------------------------------------------------------------
改变最优化算法：
利用 technique(algorithm_spec) 选项可以指定似然函数最大化的方法，Stata 提供了四种不同的算法，分别是：
technique (nr): Stata's modified Newton-Raphson (NR) algorithm (默认算法) technique (bhhh): Berndt-Hall-Hall-Hausman (BHHH) algorithm technique (dfp): Davidon-Fletcher-Powell (DFP) algorithm technique (bfgs): Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm
sysuse auto, clear
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml maximize
est store nr
ml model lf mymean_lf (price = mpg wei len) (sigma:), technique(bfgs)
ml maximize
est store bfgs
ml model lf mymean_lf (price = mpg wei len) (sigma:), technique(nr bfgs)
ml maximize
est store nr_bfgs
local model nr bfgs nr_bfgs
esttab `model', mtitle(`model')
上面三种算法的估计结果完全一致，故不再展示。
但是在以上实例中，采用 bhhh 算法会在运行过程中提示出错，从而得不到结果；采用 dfp 算法则会出现无限循环的情况，最终无法收敛。可见：
(1) 采用不同的算法，迭代次数和结果会存在很大的差异； (2) 采用不同的数值算法便于我们找到全局最大值(尤其是对于较难收敛的问题)； (3) 若同时使用多种算法，Stata 默认情况下会将每种算法迭代 5 次，以便找出最佳的迭代方式。
改变初始值：
ml maximize：将初始值设定为； ml search：Stata 自身提供了多种快速搜索初始值的方法；ml init：可以根据理论或前期估计结果设定初始值。
sysuse auto, clear
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml max //maximize 简写为 max
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml search
ml max
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml init /sigma = 3 //设定初始值
ml max
不管我们如何改变初始值的设定方法，得到的估计结果都是一致的，说明此时估计方法的稳定性较高，利用数值搜索方法得到的极大值点很可能就是全局的极大值点，模型估计出的系数是可信的。
参考资料
连享会推文，郭李鹏，「极大似然估计 (MLE) 及 Stata 实现」
相关课程
连享会-直播课 上线了！
http://lianxh.duanshu.com
免费公开课：
直击面板数据模型 - 连玉君，时长：1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT，dofiles等)
课程一览
支持回看，所有课程可以随时购买观看。
| 专题 | 嘉宾 | 直播/回看视频 |
|---|---|---|
| ⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
| ⭕ Stata数据清洗 | 游万海 | 直播, 2 小时，已上线 |
| 研究设计 | 连玉君 | 我的特斯拉-实证研究设计，-幻灯片- |
| 面板模型 | 连玉君 | 动态面板模型，-幻灯片- |
| 面板模型 | 连玉君 | 直击面板数据模型 [免费公开课，2小时] |
Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。
关于我们
Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。 连享会-主页 和 知乎专栏，300+ 推文，实证分析不再抓狂。 公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词： 课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法等
连享会小程序：扫一扫，看推文，看视频……
扫码加入连享会微信群，提问交流更方便
✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD
热门资讯
