Stata：当工具变量小于内生变量时，该如何估计？-mmeiv
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
理论 + 实证：从「读懂模型」到「折腾模型」
🎦 理论模型构建专题
📅 2021 年 10 月 2-3 日 (周六-周日)
🔑 郭凯明副教授 (中山大学)
🍓 课程主页：https://gitee.com/lianxh/emodel
🌴报名链接： http://junquan18903405450.mikecrm.com/QdtTXkm
理论模型可以简洁、凝练地抽离出经济现象的本质，使我们能够进行更深层次的思考和分析。然而，建立理论模型并非易事，若能将 理论和实证有机结合，那更加难能可贵了。
为此，我们邀请到了中山大学岭南学院郭凯明副教授，与大家一同学习理论模型的构建。郭老师一直专注于经济转型与中国经济方面的研究，发表论文近 40 篇，其中《经济研究》7 篇。
郭老师将从模型设定初衷、最基本的假设条件入手，通过讨论各种可能的建模思路和弯路，让学生不自觉中已经建立起理论分析的思维模式。最终的目标是：让学生不仅能「读懂模型」，还能「折腾模型」—— 可以自己修改甚至新设模型。
扫码直达课程主页：
作者：张子楠 (浙江财经大学)
邮箱：zinanzh@gmail.com
目录
1. 背景简介
2. mmeiv 命令
3. Stata 示例
3.1 模拟数据
3.2 具体示例
3.3 方法比较
3.4 补充说明
4 命令有效性检验
5. 参考资料
6. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 背景简介
工具变量法的可识别条件要求工具变量 (IV) 个数大于等于内生变量的个数，即所谓的工具变量秩条件。然而，在实际研究中，我们也可能碰到工具变量个数小于内生变量个数的情形。针对这种情况，Caetano and Escanciano (2021) 提出了一种思路，即通过使用工具变量与控制变量交叉项来构建新的工具变量，以满足秩条件要求。
下面以两个内生变量 ( 和 ) 和一个二元工具变量 情形来进行说明。回归方程如下所示：
结合模型和内生性情况的设定，我们设定数据生成过程 (DGP) 如下所示：
其中 为控制变量， 为遗漏变量，、、、、 为随机扰动项，且均服从 的标准正态分布。观察 DGP 可知，由于 、 和 都包含 和 ，存在由于遗漏变量导致的内生性问题。与此同时，变量 和 相关，且 和 无关， 只能通过 和 来影响 ，从而 满足工具变量的有效性和外生性要求。不妨设定 服从均值为 0.5 的二项分布。
由于这里有两个内生变量，只有一个工具变量，所以传统方法上无法进行识别。为理解这一点，我们用二阶段最小二乘法 (TSLS) 进行说明。TSLS 的第一阶段是将每个内生变量单独对所有工具变量和所有外生变量进行回归：
显然这种方法下，工具变量矩阵的秩为 1，难以识别出系数 和 。Caetano and Escanciano (2021) 思路则在于如果对于不同的 ， 与内生变量关系出现显著波动的话。那么就可以生成交叉项，并在第一阶段进行如下回归：
如果 ，即随着不同的 ， 对 的回归系数有显著不同，从而可以满足秩条件满足，我们就可以识别出回归系数。Caetano and Escanciano (2021) 还给出了命令 mmeiv 来进行上述操作。
2. mmeiv 命令
*命令安装
ssc install mmeiv, replace
*命令语法
mmeiv depvar [varlist_W2] (varlist_X = varname_T { varlist_W1 }) [if] [in] [weight] [, options]
其中，
depvar：为被解释变量；varlist_X：为内生变量；varname_T：为工具变量；varlist_W1：为用来和工具变量交乘的控制变量。注意varlist_W1要用大括号括起来；varlist_W2：为没有用来和工具变量进行交乘的控制变量，可以不写；opitions：包括五个选项，其中一个是异方差控制选项vce(vcetype)，可选robust、cluster等方式，默认是不控制。其它四个是有关回归结果的边际效应内容，和本文主题关系不大，在此不一一介绍。
3. Stata 示例
3.1 模拟数据
本文使用模拟数据来展示 mmeiv 命令的用法，模拟数据的数据生成过程和第 1 节中设定相一致。首先，设定方程真实回归系数为：
**设定环境
clear all
set seed 10101
global numobs 1000
set obs $numobs
********************************************
********* 第一部分：赋值真实回归系数值 *******
********************************************
// 方程(*)的系数设定
scalar alpha0=0.5
scalar beta01=0.8
scalar beta02=0.7
scalar beta03=3
// 方程(**)中X1的系数设定
scalar alpha01=1
scalar alpha11=1
scalar alpha21=0.5
scalar alpha31=0.2
scalar alpha41=0.1
// 方程(**)中X2的系数设定
scalar alpha02=0.9
scalar alpha12=0.1
scalar alpha22=0.1
scalar alpha32=1
scalar alpha42=0.7
// 方程(***)中的系数设定
scalar gama0=0.5
scalar gamaW=3
scalar gamaT=1.8
其次，生成模拟数据。具体如下所示：
********************************************
********** 第二部分：生成模拟数据 ************
********************************************
gen U1=rnormal()
gen U2=rnormal()
gen Uw=rnormal()
gen u=rnormal()
gen Z= rbinomial(1,0.5) //生成工具变量Z
gen T=gamaT*rnormal() //生成遗漏变量T
gen W=3*Uw+0 //生成控制变量W
gen U=gama0*W+T+u //回归方程随机扰动项U
// 最后，生成回归方程的被解释变量Y和内生变量X1、X2
gen X1=alpha01+alpha11*Z+alpha21*W+alpha31*Z*W+alpha41*T+U1
gen X2=alpha02+alpha12*Z+alpha22*W+alpha32*Z*W+alpha42*T+U2
gen Y=alpha0+beta01*X1+beta02*X2+beta03*W+U
// 删除后续演示中无需使用的变量
drop u Uw U2 U1 T
这样我们就生成了回归所需要的数据。其中被解释变量为 ，内生变量为 和 ，控制变量为 ，工具变量为 ，遗漏变量为 ，随机扰动项为 。易知有 、 和 相关，满足内生性要求。此外， 和 、 相关，满足工具变量有效性要求；且 和 不相关，满足外生性要求。
. pwcorr_a X* Z U
| X1 X2 Z U
-------------+-----------------------------------------------
X1 | 1.000
X2 | 0.669*** 1.000
Z | 0.259*** 0.061* 1.000
U | 0.567*** 0.706*** 0.002 1.000
3.2 具体示例
mmeiv 回归结果如下所示，可知 的回归系数为 0.68，95% 的置信区间为 [0.42 0.94]。虽然和真实回归系数 有所偏离，但置信区间包括了真实回归系数值。 的回归结果也相似，说明 mmeiv 回归结果还是可信的。
. ********************************************
. ********** 第三部分：mmeiv 回归 ***********
. ********************************************
. mmeiv Y (X1 X2=Z {W})
Multiple Marginal Effects IV Estimation
Number of obs = 1000
Wald chi2(3) = 45485.84
Prob > chi2 = 0.0000
Root MSE = 2.0072
------------------------------------------------------------------------------
Y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X1 | .6782149 .1332119 5.09 0.000 .4171248 .939305
X2 | .7398204 .0492009 15.04 0.000 .6433885 .8362522
------------------------------------------------------------------------------
Instrumented (X): X1 X2
Instrument (T): Z
Separable covariates (W1): W
Exogenous covariates (W2):
3.3 方法比较
作为对比，这一部分我们还做其它三种回归，看哪一种回归方法的结果更有效。第一个为不考虑内生性问题，直接 OLS 回归；第二个强行认定只有一个内生变量，用一个工具变量 进行 TSLS 回归；第三个认为是有两个内生变量，先生成 和 的交乘项 ，然而用 和 两个一起作为工具变量，进行 TSLS 回归。
第一种回归，也就是 OLS 回归的代码和结果如下所示。对于 回归系数，OLS 回归值为 0.778，95% 置信区间为 [0.682 0.874]。相对于 mmeiv 回归结果，其置信区间同样包含了真实回归系数 0.8，同时区间宽度更窄，回归系数偏误也更小。这个结果似乎表示 OLS 可能更为有效，但事实上这很大程度上是偶然因素导致。本文将在第 4 节通过 Monte Carlo 模拟来证明这点。
此外， 回归结果为 1.172，偏离真实值 0.7，同时置信区间也没有包含真实值 0.7。OLS 回归对 的识别有效性要显著弱于 mmeiv 回归。
. ********************************************
. ******** 第四部分：与其它回归方法的比较 ******
. ********************************************
. // 回归对比1：OLS 回归结果
. reg Y X1 X2 W
Source | SS df MS Number of obs = 1,000
-------------+---------------------------------- F(3, 996) = 20599.54
Model | 190074.927 3 63358.309 Prob > F = 0.0000
Residual | 3063.41247 996 3.07571533 R-squared = 0.9841
-------------+---------------------------------- Adj R-squared = 0.9841
Total | 193138.34 999 193.331671 Root MSE = 1.7538
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X1 | 0.778 0.049 15.85 0.000 0.682 0.874
X2 | 1.172 0.026 44.63 0.000 1.121 1.224
W | 3.224 0.036 90.40 0.000 3.154 3.294
_cons | 0.084 0.091 0.92 0.357 -0.095 0.264
------------------------------------------------------------------------------
第二种回归是强行认定只有一个内生变量，并且只用一个工具变量 进行 TSLS 回归。回归代码和结果如下所示，其中 的回归系数为 0.58，置信区间包括了真实值，但 的回归系数则存在较大偏误，置信区间也没有包含真实回归系数值。
. ivregress 2sls Y X2 W (X1 =Z)
Instrumental variables (2SLS) regression Number of obs = 1,000
Wald chi2(3) = 60859.31
Prob > chi2 = 0.0000
R-squared = 0.9839
Root MSE = 1.764
------------------------------------------------------------------------------
Y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X1 | 0.583 0.116 5.03 0.000 0.356 0.811
X2 | 1.199 0.030 39.71 0.000 1.140 1.258
W | 3.326 0.066 50.66 0.000 3.197 3.455
_cons | 0.350 0.170 2.06 0.040 0.016 0.683
------------------------------------------------------------------------------
Instrumented: X1
Instruments: X2 W Z
第三，用 和 交乘项 ，以及 共同作为工具变量来进行 IV 回归。回归代码和结果如下所示，观察回归结果可以发现，无论是回归系数、标准误、置信区间以及 等，均与 mmeiv 回归结果完全一致。这并非偶然，事实上如第 1 节介绍的一样，mmeiv 命令的解决思路，就在于构建交乘项 ，然后再将 和 同时作为工具变量，进行两阶段最小二乘法工具变量回归的。
. // 生成Z和W交乘项，命名为ZW
. gen ZW=Z *W
. // 利用工具变量Z ZW进行两阶段二乘法IV回归
. ivregress 2sls Y W (X1 X2=Z ZW)
Instrumental variables (2SLS) regression Number of obs = 1,000
Wald chi2(3) = 45485.84
Prob > chi2 = 0.0000
R-squared = 0.9791
Root MSE = 2.0072
------------------------------------------------------------------------------
Y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X1 | 0.678 0.133 5.09 0.000 0.417 0.939
X2 | 0.740 0.049 15.04 0.000 0.643 0.836
W | 3.567 0.074 48.11 0.000 3.422 3.712
_cons | 0.719 0.191 3.75 0.000 0.344 1.094
------------------------------------------------------------------------------
Instrumented: X1 X2
Instruments: W Z ZW
3.4 补充说明
补充说明 1：在利用协变量 生成交乘项时， 要满足协变量完备性 (convaiance completeness)，即 包含元素个数至少要比内生变量多一个。相关证明见 Caetano and Escanciano (2021)。
补充说明 2：在这里，我们并没有要求 是外生的。事实上从回归方程随机扰动项 的生成过程 可知， 并非外生的，但我们仍旧获得了 和 的一致回归结果。相对于再寻找一个好的 IV，mmeiv 命令只需要一个满足完备性的协变量，甚至允许该协变量可以不外生，从而大大降低了识别难度。
补充说明 3：上述例子可以很容易扩展 N 个内生变量但只有 1 个工具变量情形。只需要找到 N-1 个满足完备性的协变量，然后放入 { varlist_W1 } 地方再回归即可。
补充说明 4：Caetano and Escanciano (2021) 的思路也为一个常见的 IV 处理方法提供了理论支持。在计量应用文献里，有时候工具变量是不随时间而变动的变量，如小河流长度 (River) 等。如果这时再同时控制时间固定效应进行面板回归，工具变量会在组内去芯过程中被消除掉。对于这种问题，一个解决方法是将小河流长度交乘年份 (Year) 生成新的变量 ，然后用 作为新的工具变量进行工具变量回归。这种处理方法和 Caetano and Escanciano (2021) 的处理思路是完全一致的。进一步理解，我们不仅可以用 来交乘，我们也可以用其他协变量来交乘。(注：读者可思考一下，这里是否还要求 这类协变量满足协变量完备性？)
补充说明 5：mmeiv 回归还可以通过绘制图形，来表示工具变量回归的边际效应，特别是工具变量维度小于内生变量维度的情形。比如工具变量是政策这种 0-1 变量，而内生变量是连续变量的情形。下面来展示相应用法，代码和边际效应图分别如下所示：
. mmeiv Y (X1 X2=Z {W}), plot
4 命令有效性检验
在第 3 节，mmeiv 回归结果中 回归系数为 0.68，与真实系数 0.8 存在较大偏误，也低于 OLS 回归结果 (0.778)。为了看这种偏误是否是偶然因素，我们接下来通过 Monto Carlo 模拟来进行检验。设定模拟次数为 1000。模拟代码如下所示：
**********************************************
***********第五部分 1：mmeiv 的 MC 模拟********
**********************************************
** 注意1：回归前，可先将《第一部分：赋值真实回归系数值》的代码运行一遍。
** 注意2：本部分代码需要全选中，然后运行。
global numsims "1000" // number of simulations
capture program drop mscmmeiv
program mscmmeiv, rclass
version 11
drop _all
global numobs 1000 // sample size N
set obs $numobs
gen U1=rnormal()
gen U2=rnormal()
gen Uw=rnormal()
gen u=rnormal()
gen Z= rbinomial(1,0.5) //生成工具变量Z
gen T=gamaT*rnormal() //生成遗漏变量T
gen W=gamaW*Uw+0 //生成控制变量W
gen U=gama0*W+T+u //回归方程的随机扰动项U。
gen X1=alpha01+alpha11*Z+alpha21*W+alpha31*Z*W+alpha41*T+U1
gen X2=alpha02+alpha12*Z+alpha22*W+alpha32*Z*W+alpha42*T+U2
gen Y=alpha0+beta01*X1+beta02*X2+beta03*W+U
mmeiv Y (X1 X2=Z {W})
return scalar b2 =_b[X1] //提取X1的MSC情况
return scalar se2 = _se[X1]
return scalar t2 = (_b[X1]-beta01)/_se[X1]
return scalar r2 = abs(return(t2))>invttail($numobs-2,.025)
return scalar p2 = 2*ttail($numobs-2,abs(return(t2)))
return scalar b22=_b[X2] //提取X2的MSC情况
return scalar se22= _se[X2]
return scalar t22= (_b[X2]-beta02)/_se[X2]
return scalar r22= abs(return(t2))>invttail($numobs-2,.025)
return scalar p22= 2*ttail($numobs-2,abs(return(t2)))
end
simulate b2r=r(b2) se2r=r(se2) reject2r=r(r2) p2r=r(p2) ///
b2r2=r(b22) se2r2=r(se22) p22r=r(p22), ///
reps($numsims) nolegend nodots: mscmmeiv
. // 1000次模拟后获得的X1的回归系数均值(b2r)、以及拒绝回归系数等于0.8的p值(p2r)
. mean b2r p2r
Mean estimation Number of obs = 1,000
--------------------------------------------------------------
| Mean Std. err. [95% conf. interval]
-------------+------------------------------------------------
b2r | 0.801 0.004 0.793 0.809
p2r | 0.508 0.009 0.490 0.526
--------------------------------------------------------------
.
. // 1000次模拟后获得的X2的回归系数均值、以及拒绝回归系数等于0.7的p值
. mean b2r2 p22r
Mean estimation Number of obs = 1,000
--------------------------------------------------------------
| Mean Std. err. [95% conf. interval]
-------------+------------------------------------------------
b2r2 | 0.699 0.002 0.696 0.702
p22r | 0.508 0.009 0.490 0.526
--------------------------------------------------------------
1000 次模拟结果显示， 回归系数均值为 0.8，和真实系数值 0.80 一致。同时 95% 置信区间为 [0.792 0.809]，包含了真实系数。另外，拒绝回归系数等于 0.8 的 值为 0.505，无法拒绝。类似地，对于 ，回归系数均值为 0.702，同时 95% 置信区间为 [0.699 0.705]，包含了真实系数 0.7。另外，拒绝回归系数等于 0.7 的 值为 0.505，无法拒绝。可见，第 3.2 节中的偏误很大程度是偶然因素所导致。从 MC 模拟结果可以看出，mmeiv 命令的估计量具有较好的有效性。
作为对比，我们也对 OLS 回归进行 MC 模拟，代码如下所示。MC 模拟结果显示， 回归系数均值为 0.797， 置信区间为 [0.793 0.800]， 回归系数均值为 1.17，置信区间为 [1.168 1.171]。显然，OLS 估计有效性要低于 mmeiv 命令。
*********************************************
********第五部分 2：OLS 的 MC 模拟 ***********
*********************************************
** 注意1：回归前，可先将《第一部分：赋值真实回归系数值》的代码运行一遍。
** 注意2：本部分代码需要全选中，然后运行。
global numsims "1000" // number of simulations
capture program drop mscols
program mscols, rclass
version 11
drop _all
global numobs 1000
set obs $numobs
gen U1=rnormal()
gen U2=rnormal()
gen Uw=rnormal()
gen u=rnormal()
gen Z= rbinomial(1,0.5)
gen T=gamaT*rnormal()
gen W=gamaW*Uw+0
gen U=gama0*W+T+u
gen X1=alpha01+alpha11*Z+alpha21*W+alpha31*Z*W+alpha41*T+U1
gen X2=alpha02+alpha12*Z+alpha22*W+alpha32*Z*W+alpha42*T+U2
gen Y=alpha0+beta01*X1+beta02*X2+beta03*W+U
reg Y X1 X2 W
return scalar b2 =_b[X1] //提取X1的MSC情况
return scalar se2 = _se[X1]
return scalar t2 = (_b[X1]-beta01)/_se[X1]
return scalar r2 = abs(return(t2))>invttail($numobs-2,.025)
return scalar p2 = 2*ttail($numobs-2,abs(return(t2)))
return scalar b22=_b[X2] //提取X2的MSC情况
return scalar se22= _se[X2]
return scalar t22= (_b[X2]-beta02)/_se[X2]
return scalar r22= abs(return(t2))>invttail($numobs-2,.025)
return scalar p22= 2*ttail($numobs-2,abs(return(t2)))
end
simulate b2r=r(b2) se2r=r(se2) reject2r=r(r2) p2r=r(p2) ///
b2r2=r(b22) se2r2=r(se22) p22r=r(p22), ///
reps($numsims) nolegend nodots: mscols
. // 1000次模拟后获得的X1的回归系数均值、以及拒绝回归系数等于0.8的p值。
. mean b2r p2r
Mean estimation Number of obs = 1,000
--------------------------------------------------------------
| Mean Std. err. [95% conf. interval]
-------------+------------------------------------------------
b2r | 0.797 0.002 0.793 0.800
p2r | 0.487 0.009 0.468 0.505
--------------------------------------------------------------
. // 1000次模拟后获得的X2的回归系数均值、以及拒绝回归系数等于0.7的p值。
. mean b2r2 p22r
Mean estimation Number of obs = 1,000
--------------------------------------------------------------
| Mean Std. err. [95% conf. interval]
-------------+------------------------------------------------
b2r2 | 1.169 0.001 1.168 1.171
p22r | 0.487 0.009 0.468 0.505
--------------------------------------------------------------
5. 参考资料
Caetano C, Escanciano J C. Identifying multiple marginal effects with a single instrument[J]. Econometric Theory, 2021, 37(3): 464-494. -PDF-
Caetano C. & Escanciano J. C.(2020),"IDENTIFYING MULTIPLE MARGINAL EFFECTS WITH A SINGLE INSTRUMENT", Econometric Theory 37(3):464-494.PDF
陈勇吏，连享会推文，Stata：蒙特卡洛模拟分析 (Monte Carlo Simulation)
6. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 工具变量 蒙特卡洛模拟, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Lasso一下：再多的控制变量和工具变量我也不怕-T217 Stata新命令-pdslasso：众多控制变量和工具变量如何挑选？ 专题：Stata程序 Stata：蒙特卡洛模拟A-(Monte-Carlo-Simulation)没那么神秘 Stata：蒙特卡洛模拟分析 (Monte Carlo Simulation) 专题：IV-GMM 找不到IV？基于异方差构造工具变量 IV在哪里？奇思妙想的工具变量 twostepweakiv：弱工具变量有多弱？ 多个(弱)工具变量如何应对-IV-mivreg？ IV：工具变量不满足外生性怎么办？ IV-工具变量法：第一阶段系数符号确定时的小样本无偏估计 IV：可以用内生变量的滞后项做工具变量吗？ Stata: 工具变量法 (IV) 也不难呀！ IV-估计：工具变量不外生时也可以用！ 专题：内生性-因果推断 工具变量-IV：排他性约束及经典文献解读
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🎦 🎦 🎦
🎧 项目申报书专题：自科、社科和教育部基金项目
📅 2021 年 9 月 19-21 日
🔑 迟老师；孙老师；张老师
🎤 讲授方式：网络直播，每天 7 小时（8:00-12:00；14:00-17:00）
8:00-12:00 ：申报书写作要领 14:00-15:30：经典申报书解读 15:30-17:00：答疑解惑
🍓 课程主页：https://gitee.com/lianxh/kt
🍏 关于我们
连享会 ( www.lianxh.cn ) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下：连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。