Stata：控制函数法及其Stata应用
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会课程 · 2023 暑期班
作者：张弛 (中山大学)
邮箱：zhangch369@mail2.sysu.edu.cn
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 何时选择控制函数法？
2. 常系数线性模型
2.1 EEV 为连续型变量
2.2 EEV 为离散型变量
3. 相关随机系数模型
3.1 模型简介
3.2 模型设定
3.3 操作步骤和 Stata 实现
3.4 模型拓展
4. 非线性模型
4.1 EEV 是连续的
4.2 EEV 是非线性的
5. 控制函数法与工具变量法比较
6. 参考文献
7. 相关推文
1. 何时选择控制函数法？
在对政策效果进行评估时，选择参与到一个项目的主体可能在某些方面与选择不参与到这个项目的主体存在系统性差异，并且这样的差异很难被观察到，这会导致内生性问题。
例如，Sapelli and Vial (2002) 在研究智利上公立学校和私立学校对二年级学生成绩影响时，不同的父母教育程度、不同的家庭经济状况等因素，一方面会直接影响到学生的成绩，另一方面还会影响到学生是选择读公立学校还是私立学校。因此，“选择学校类型”这一变量本身就存在内生性。
控制函数法的核心思想就在于，通过对内生解释变量 (Endogenous Explanary Variables，下文简称 EEV) 的内生性来源进行提取，然后加入原始模型进行估计，来获得对政策效果的无偏估计。
2. 常系数线性模型
2.1 EEV 为连续型变量
2.1.1 原理简介
假设 为被解释变量， 为 EEV， 是 的外生变量 (其中， 是 的子向量)，我们估计的原始方程为：
外生性假设为：
对 的内生性假设：
其中， 是 的向量。此时，定义残差项之间的关系：
其中， 为总体回归系数。将 (4) 代入 (1)，得：
此时，该方程所估计的 即与残差项无关。因为我们将 中与残差项有关的，存在内生性的部分 已单独提取出，作为一个新变量加入我们的估计模型中。
2.1.2 检验变量是否为 EEV
原假设 相当于 Hausman 检验中
在这样的线性方法中，对 和 的估计和 2SLS 方法是一样的。注意 IV 方法和 CF 方法不一定等价。例如，原方程加入 的平方项
若我们使用IV方法，我们需要加入 作为 的工具变量，此时工具变量向量为 。但是若使用 CF 方法，则需要施加更强的假设：
在 CF 方法中，第一阶段回归为：
在前述条件下，我们有新的估计方程 (11)，此时，与 2SLS 的估计方法不一致。
用第一阶段回归所估计得到的 代入上式，用 OLS 方法估计。
这样的假设会使得 CF 方法更有效，但是稳健性不如 IV 方法，因为 IV 方法只需要满足 ，这个假设更弱。
2.1.3 估计方法和 Stata 实现
第一步：将 对 做回归，获得简化方程的残差项 。第二步：将 对 ，， 做回归，对 (6) 式的 OLS 估计就是 Control Function 的估计方法。在第二步的估计中，，， 都是一致的。
. bcuse wage2 //一个关于教育回报的数据集，作为演示的范例
. gen exper2=exper^2
. sum educ
. gen educ_average=r(mean)
. // 控制函数法stage1 获取存在内生性的残差项
. bootstrap, reps(1000): reg educ meduc feduc black
. predict v2, residual //残差项吸收了内生性
. // 教育回报率方程其他变量的设置
. gen black_educ=black*(educ-educ_average)
. // stage2 将内生的残差项加入原始模型回归
. bootstrap, reps(1000): reg lwage educ exper exper2 black tenure black_educ v2
Linear regression Number of obs = 722
Replications = 1,000
Wald chi2(7) = 165.62
Prob > chi2 = 0.0000
R-squared = 0.1824
Adj R-squared = 0.1744
Root MSE = 0.3811
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
lwage | coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
educ | 0.120 0.015 7.96 0.000 0.090 0.149
exper | 0.008 0.016 0.50 0.620 -0.024 0.040
exper2 | 0.000 0.001 0.66 0.508 -0.001 0.002
black | -0.149 0.057 -2.61 0.009 -0.261 -0.037
tenure | 0.009 0.003 3.18 0.001 0.004 0.015
black_educ | -0.030 0.026 -1.19 0.234 -0.080 0.020
v2 | -0.056 0.016 -3.59 0.000 -0.087 -0.026
_cons | 4.946 0.230 21.48 0.000 4.495 5.398
------------------------------------------------------------------------------
2.2 EEV 为离散型变量
2.2.1 原理方法
CF 方法相比于 IV 方法的优势在于，它可以针对非线性的 EEV 进行一致性估计，例如 EEV 可以是一个二值变量。当 为非线性形式时，我们无法使用 IV 估计，CF 方法是更好的选择。例如，待估模型为：
在该模型中，我们假设 为二值变量，满足以下方程：
运用 CF 方法，需施加以下假设：
此时， 取决于 的联合分布，且 独立于 ，，此时我们估计的方程为：
并且
其中， 是逆米尔斯比率 (inverse Mills ratio，IMR)。
2.2.2 估计步骤和 Stata 实现
第一步：估计 ，获得广义残差：
第二步：将 对 ，， 进行 OLS 回归，获得对 ，， 的一致估计。
. webuse drugexp, clear
. etregress lndrug chron age lninc, treat(ins=age married lninc work) poutcomes cfunction
Linear regression with endogenous treatment Number of obs = 6,000
Estimator: Control function
------------------------------------------------------------------------------
| Robust
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
lndrug |
chron | 0.467 0.032 14.61 0.000 0.405 0.530
age | 0.102 0.003 34.98 0.000 0.096 0.108
lninc | 0.055 0.023 2.45 0.014 0.011 0.099
1.ins | -0.860 0.348 -2.47 0.014 -1.543 -0.177
_cons | 1.666 0.253 6.59 0.000 1.170 2.161
-------------+----------------------------------------------------------------
ins |
age | 0.021 0.002 9.21 0.000 0.017 0.026
married | 0.085 0.036 2.35 0.019 0.014 0.155
lninc | 0.102 0.023 4.55 0.000 0.058 0.146
work | 0.288 0.037 7.75 0.000 0.215 0.361
_cons | -0.623 0.109 -5.73 0.000 -0.836 -0.410
-------------+----------------------------------------------------------------
/athrho0 | 0.404 0.172 2.34 0.019 0.066 0.742
/lnsigma0 | 0.316 0.050 6.31 0.000 0.218 0.414
/athrho1 | 0.793 0.299 2.66 0.008 0.208 1.378
/lnsigma1 | 0.187 0.061 3.04 0.002 0.066 0.307
-------------+----------------------------------------------------------------
rho0 | 0.383 0.147 0.065 0.630
sigma0 | 1.372 0.069 1.243 1.513
lambda0 | 0.525 0.226 0.082 0.969
rho1 | 0.660 0.169 0.205 0.881
sigma1 | 1.205 0.074 1.069 1.359
lambda1 | 0.795 0.251 0.303 1.288
------------------------------------------------------------------------------
Wald test of indep. (rho0 = rho1 = 0): chi2(2) = 8.88 Prob > chi2 = 0.0118
. // 在该估计中，内生解释变量为ins，它受到二值变量age，married，lninc，work的影响
. // poutcomes：表示用潜在结果模型估计第一阶段回归，如果不加该选项，即表示第一阶段回归为线性回归
. // cfunction：表示用控制函数法估计处理效应
2.2.3 EEV 是非线性的情况
假设 和其他外生变量是交互的：
如果 中， 是线性的，可以按照最原始的方法，如果是二值变量，可以使用离散型随机变量的方法，均为以上两种情况的简单变形。
3. 相关随机系数模型
3.1 模型简介
有时候， 的效果也会受到不可观测的因素影响。例如，就读于教会学校的教育回报率，可能会因为入学学生的个体差异而不同，且这些不同是无法观测的。
相关随机系数 (correlated random coefficient，CRC)，指允许一个随机的系数影响解释变量。这个概念最早由 Heckman and Vytlac 在1998年提出。在研究处理效应时， CRC 模型允许异质性的处理效应和自选择结合。
3.2 模型设定
假设估计模型：
其中， 与 ， 相关，并且 EEV 系数满足以下两个条件：
代入 (18) 得：
与此同时，针对 (21) 式，我们还需要施加如下的外生性假设：
注意，此处不能运用 2SLS 估计，因为此时估计出来的 是不一致的，残差为 ，仍然内生。
假设 的内生性设定如下：
存在不可观测因素 , ，他们和 满足以下线性关系：
并且，我们认为，所有不可观测的变量都与 独立，所以，我们实际上需要估计的式子是：
接着，我们将 对 做 OLS 回归，得到CF估计。此处 的系数，衡量的就是随机系数 对 的影响。但是这里必须注意，如果 的系数不为 0，对交互项 的 检验在第一阶段的回归是不可靠的，需要在第二步使用 bootstrap 去调整方差。
在这里，为了方便，我们是对交互项和残差项 () 做联合检验，检验是否存在外生性。
3.3 操作步骤和 Stata 实现
操作步骤：与线性常系数模型中EEV连续的情况基本类似，只是在第二阶段回归的时候加入 EEV 和残差的交互项。在这里，我们依然考虑一个讨论教育回报率的模型。
. bcuse wage2, clear
. gen exper2=exper^2
. sum educ
. gen educ_average=r(mean)
. // 控制函数法 stage1 获取存在内生性的残差项
. bootstrap, reps(1000): reg educ meduc feduc black
. predict v2, residual // 残差项吸收了内生性
. // 教育回报率方程其他变量的设置
. gen v2_educ=v2*educ // 构造残差和 eev 的交互项
. gen black_educ=black*(educ-educ_average)
. // stage2 将内生的残差项加入原始模型回归
. bootstrap,reps(1000): reg lwage educ exper exper2 black tenure black_educ v2 v2_educ
Linear regression Number of obs = 722
Replications = 1,000
Wald chi2(8) = 187.97
Prob > chi2 = 0.0000
R-squared = 0.1866
Adj R-squared = 0.1775
Root MSE = 0.3803
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
lwage | coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
educ | 0.124 0.014 8.56 0.000 0.095 0.152
exper | 0.007 0.016 0.42 0.672 -0.024 0.037
exper2 | 0.001 0.001 0.83 0.407 -0.001 0.002
black | -0.160 0.057 -2.82 0.005 -0.272 -0.049
tenure | 0.010 0.003 3.20 0.001 0.004 0.016
black_educ | -0.036 0.026 -1.37 0.170 -0.087 0.015
v2 | 0.040 0.050 0.80 0.423 -0.058 0.137
v2_educ | -0.007 0.003 -2.02 0.043 -0.013 -0.000
_cons | 4.922 0.215 22.86 0.000 4.500 5.344
------------------------------------------------------------------------------
3.4 模型拓展
(28) 式的模型也可以拓展到非线性的 (, )，包括多次方项和交互项。更一般化的模型如下：假设存在一个 的 ，存在这样一个方程 ， 的随机系数为 。此时，我们估计的方程为：
在第一阶段回归获得残差 后，我们将 对 ，，，， 做 OLS 回归，并使用非参数的 bootstrap 方法调整标准差，从而得到有效的估计。如果 是 的 维的向量，我们可以估计 。
类似的，对于 是二值变量的情况，该方法也一样成立。由前所示， 需满足前面的条件：
然后，获取标准化的残差：
最后，将 对 做回归，中心化的 帮助我们求出政策的平均效应。注意，针对 (26) 式的估计，我们不直接讨论 和 两种情况下分别的效果，ATE (Average Treatment Effect) 在这个模型下就是 。
4. 非线性模型
当 EEV 是连续的，CF 方法经常被运用于估计 Probit 和 Tobit 模型。
4.1 EEV 是连续的
4.1.1 Rivers-Vuong 估计方法
假设我们要估计的 为 Probit 模型：
EEV 是连续的：
并且，(，) 服从二元正态分布，期望为 0， 且与 独立。令 ，此时我们用 CF 方法的估计为：
在此模型之中每个系数都被乘以了 。
4.1.2 操作步骤和 Stata 实现
第一步：用 OLS 方法将 对 做回归，得到残差 。第二步：将 对 ,， 用 Probit 方法估计系数。并且，原假设 H0: ，对 做 检验。
若想求平均处理效应 (APE，average partial effects)，其实就是要求下面这个式子的变化：
从而可以求出原始模型 的估计值 (ASF，average structural function)：
. lxhuse mus17data.dta, clear
. // 控制函数法 stage1 获取存在内生性的残差项
. global xlist2 medicaid age age2 educyr actlim totchr
. regress private $xlist2 income ssiratio, vce(robust)
. predict lpuhat, residual
. // 控制函数法 stage2 将残差项加入原模型回归
. poisson docvis private $xlist2 lpuhat, vce(robust) nolog
Poisson regression Number of obs = 3,677
Wald chi2(8) = 718.87
Prob > chi2 = 0.0000
Log pseudolikelihood = -15010.614 Pseudo R2 = 0.1303
------------------------------------------------------------------------------
| Robust
docvis | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
private | 0.551 0.245 2.24 0.025 0.070 1.031
medicaid | 0.263 0.120 2.20 0.028 0.028 0.498
age | 0.335 0.070 4.81 0.000 0.199 0.471
age2 | -0.002 0.000 -4.79 0.000 -0.003 -0.001
educyr | 0.019 0.008 2.31 0.021 0.003 0.034
actlim | 0.205 0.041 4.96 0.000 0.124 0.287
totchr | 0.241 0.013 18.69 0.000 0.216 0.267
lpuhat | -0.417 0.249 -1.67 0.095 -0.905 0.072
_cons | -11.906 2.661 -4.47 0.000 -17.123 -6.690
------------------------------------------------------------------------------
4.1.3 拓展：Blundell and Powell 估计方法
BP 方法允许 是一个和 有关的函数，且不一定是线性的。BP 方法假设一个结构模型：
它的 ASF 定义为：
这个 ASF 随着 (，) 的变化而变化，从而求出 APE。另一个关键假设是关于 ：
并且，() 与 是独立的，且 的期望等于 0，从而 条件在 的期望为 。() 独立于 时，不可观测的 与 () 相关的部分，只体现在 上，即满足：
因此， 可以作为 的代理变量，由此，第一步，定义如下条件期望：
我们要求出的核心结果是：
因此，要先求出:
然后根据 ASF 估计量的定义，我们求出：
接下来，我们放宽假定，考虑一个二值响应模型。令 是所有内生和外生变量的一个线性组合， 是对 约简式的一个误差项，我们可以估计一个二值相应模型，这个模型的设定是灵活的。举个例子：
在这里，CF 方法的第一步要求我们从最初的简约式中估计出 ，第二步要求我们将 对 做 Probit 估计。检验外生性的方法也是类似的，对后三项做联合检验。之后，我们可以估计 ASF：
4.2 EEV 是非线性的
前面所述的 RV 方法无法适用于 是非连续的情况。在传统方法中，我们需要用 MLE 的方法去估计。接下来，我们介绍一种使用 CF 估计的方法。假设内生解释变量是一个二值变量：
在标准二值变量的 Probit 假设中，没有一致估计参数的 CF 方法，但是可以通过正常的 MLE 方法中做对标准化残差的t检验，来判断 是否外生。在第一步，我们如前所述，需要估计：
但是在这里，我们不需要通过估计 的系数来检验内生性，而是通过估计下面这个方程，原假设 ：
此时，用 去替代 。同样的，我们需要假设 的内生性全部被 捕捉，我们加入和前文类似的假设：
接着，再根据 BP 方法，估计 ASF：
和连续 EEV 的情况类似，我们一样也可以使用更灵活的方法去估计以下方程：
其中， 是 () 的函数。此时，在用 去替代 后，我们使用 Wald Test，检验 ，。此时，我们估计的 ASF 和前面的方法类似，将 替换成 即可。
控制函数法对模型估计的灵活设置，有助于解决更复杂的问题。例如，一个二值估计模型之中，同时存在一个连续的 EEV 和一个二值 EEV ，我们在第一阶段估计中估计出来的残差或者标准化残差，可以在第二阶段中以各种形式加入 的 probit 模型估计，例如二次项、三次项或者交互项。
5. 控制函数法与工具变量法比较
假设条件：
控制函数法：在两个阶段做回归时，残差都需要服从一个特定分布，但无需满足排斥性约束； IV：第一阶段对内生变量做回归时，残差无需服从一个特定分布，但需要满足排斥性约束。
估计过程：
控制函数法：模型设定更为灵活，它可以针对线性，非线性，常系数和随机系数进行估计； IV：对模型设定的要求不如 CF 灵活。例如，EEV 为二值变量时，CF 方法在保持一致估计的情况下，比传统 2SLS 方法更有效。但是，在多重非线性方程下的 EEV，IV 方法在一致性方面更加稳健，但是也有可能有效性降低。
估计结果：
控制函数法：可以求出 ATT ( average effect of the treatment on the treated )，ATE( average treatment effect )，LATE ( local average treatment effect)； IV：只能求 LATE。
6. 参考文献
Jeffrey M. Wooldridge, 2015. "Control Function Methods in Applied Econometrics", Journal of Human Resources, University of Wisconsin Press, vol. 50(2), pages 420-445. -PDF- Jeffrey M. Wooldridge, 2011. "Control Function And Related Methods". -PDF- Glewwe, P., P. Todd. Control function methods[C], Impact evaluation in international development: Theory, methods, and practice 267-279. -PDF- Heckman, J., S. Navarro-Lozano, 2004, Using matching, instrumental variables, and control functions to estimate economic choice models, Review of Economics and Statistics, 86 (1): 30-57. -PDF- Altonji, Joseph G. et al. “An Evaluation of Instrumental Variable Strategies for Estimating the Effects of Catholic Schooling.” The Journal of Human Resources XL (2002): 791 - 821. -PDF- Sapelli, Claudio, and Bernardita Vial. 2002. “The performance of private and public Schools in the Chilean Voucher System.” Cuardenos do Economía 39 (118): 423–54. -PDF- Jeffrey M. Wooldridge, 横截面与面板数据的经济计量分析[M]. 中国人民大学出版社, 2007. 计量经济圈推文：非线性模型及离散内生变量处理利器-应用计量经济学中的控制函数法！-Link-
7. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 内生性, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：IV-GMM IV专题- 内生性检验与过度识别检验 IV专题: 内生性检验与过度识别检验 专题：内生性-因果推断 内生性！内生性！解决方法大集合 第三种内生性：衡量偏误(测量误差)如何解决？-eivreg-sem 内生性：来源及处理方法-幻灯片下载 locmtest：非线性模型的内生性检验 第三种内生性：衡量偏误(测量误差)如何检验-dgmtest？ 专题：分位数回归 Stata：面板分位数模型估计及内生性初探
课程推荐：2023 暑期班
主讲老师：连玉君，王群勇
🍓 课程主页：https://www.lianxh.cn/news/fdc69c3695aec.html
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。