Stata：内生变量与工具变量非线性关系处理-discretize
为此，我们邀请到了中山大学岭南学院郭凯明副教授，与大家一同学习理论模型的构建。郭老师一直专注于经济转型与中国经济方面的研究，发表论文近 40 篇，其中《经济研究》7 篇。
郭老师将从模型设定初衷、最基本的假设条件入手，通过讨论各种可能的建模思路和弯路，让学生不自觉中已经建立起理论分析的思维模式。最终的目标是：让学生不仅能「读懂模型」，还能「折腾模型」—— 可以自己修改甚至新设模型。
扫码直达课程主页：
作者：吕文渊 (华盛顿大学/韩国发展研究院公共政策学院 KDI School)
邮箱：lwy1022@uw.edu
目录
1. 理论背景
2. 命令介绍
3. 论文举例
4. 参考文献
5. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 理论背景
首先，我们先来复习一下工具变量 (Instrument Variable) 的基本知识。最简单的 OLS 回归模型：
其中自变量 与误差项 必须互相独立，即自变量必须是外生的，不能由模型本身所决定。然而，现实中的数据往往很难符合 OLS 模型中的严格外生假设。如果误差项 与 之间并不独立， 则对系数 的估计就会出现偏误。我们将这样的偏误便称之为内生性偏误 (Endogeneity Bias)。进一步，产生内生性偏误的主要原因有三个：
遗漏变量 (Omitted variable)； 测量误差 (Measurement error)； 互为因果 (Simultaneity)。
对于内生性问题最直接的解决方法之一是工具变量的使用 (下图所示)。因为工具变量可以通过两阶段最小二乘法将 的内生性去除掉，从而使得使得误差项能够独立于 。但是，如果深究一下，使用工具变量本身便存在着诸多的隐患——由于缺乏先验知识或者信息 (Prior Knowledge)，我们并不能彻底了解 与 之间的确切关系。也就是说，工具变量和自变量之间 (第一阶段，first-stage) 依然是一个黑箱中的关系，我们仍然需要对两者的关系进行处理。
对此，有两种可行的方案：(1) 将 与 之间的关系估计成为一个复杂的非线性模型。这样处理的优点是，第一阶段的假设较为强假设。同时，这也是一个纯数据驱动的过程 (data-driven procedure)，并不需要考虑其他因素。但缺点也是显而易见，这样的做法存在过拟合 (Overfitting) 问题——误差项的方差会过大，模型存在没有必要的复杂程度。(2) 将 与 之间的关系估计成为一个线性模型。优点是过拟合问题不再存在。但缺点是，如果两者关系是非线性的，则工具变量 效果便大打折扣，并且很可能会被认为是弱工具变量 (Weak In strument Variable)。
至此，可以说两种方案都不够令人满意，所以学者们考虑采取第 (3) 种方案：将连续型工具变量转换成为二值虚拟变量。这样处理的好处颇多：
首先，这是简约非参数模型 (Parsimonious non-parametric model， Angrist 和 Pischke，2009)； 其次，这样的处理便于解释方差的来源。但是风险也依然存在，即二值工具变量的边界如何确定是比较随意的，这取决于研究者和研究本身，以及第二阶段估计的敏感度会存在偏高的情况。
2. 命令介绍
Curci 等 (2019) 开发了将连续型工具变量转换为虚拟变量的命令 discretize。该命令能提供一个数据驱动过程来构建一个离散的工具变量，并且边界取的是第一阶段 F-statstics 最大的值，这解决了上文所说的边界决定因人而已异的问题。该命令的主要优势在于：
第一阶段的简约非参数模型 (Parsimonious non-parametric model)； 最小化了弱工具变量的问题； 不依靠研究者的自行决策，过程透明； 可以用图像来来检验第二阶段的稳定性。
*命令安装
lxhget discretize.ado, replace //将 discretize.ado 放到 Stata 对应的程序文件夹下
*第一阶段估计的语法discretize contvarname, endogenous(varname) range(min/max) interval(min(step)max)contvarname：连续型工具变量 (取整数，因为循环运算不能很好地处理小数)；Endogenous(varname)：内生变量；range(min/max)：取值范围；Interval(min(step)max)：检验虚拟变量 “门槛” 的最小和最大的区间，和检测步长 (step)。
*第二阶段估计的语法
discretize contrvarname, endogenous(varname)
range(min/max) interval(min(step)max)
interval(min(step)max) second depvar(varname)
在第二阶段估计中，仅加入了第二阶段因变量命令 second depvar，并使用 ivregress 进行两阶段最小二乘法的估计。在这个主命令下还有其他可使用的选项：
exogenous (varlist)：表示在一阶段和第二阶段使用到的外生变量 (可多个)；interact(varname)：与 IV 相关的交互变量；xt(estimator)：面板数据估计也可以使用xtreg和xtivreg；vce(vcetype)：增加聚类标准误；print：显示结果值；save：保存结果文件 (结果是在 95% 置信区间下)；graph：呈现系数估计或 F 统计量的图像。
3. 论文举例
该命令是 Curci 和 Masera (2018) 的产物。如果要举例说明命令 discretize 的使用，莫过于使用这篇论文作为例子。我们先来解一下作者论文的故事大概。该文主要是研究城市的 “郊区化” (suburbanization，即大量人口从城市中心转移到郊区的过程)，作者认为这种郊区化趋势是可以用犯罪率来解释的——由于城市中心的犯罪率高，导致富裕阶级会迁移至郊区。但是犯罪率与郊区化趋势存在内生性问题，为此作者选取土壤灰尘里的铅含量 (lead poisoning) 作为工具变量来克服上述内生性问题。之所以选取该工具变量，主要是因为：
重金属沉积在人体，会较高概率引致行为失当 (violent crime)； 重金属暴露的主要原因是汽车尾气排放 (car emission)； 尾气中的有害成分四乙铅 (tetraethyl lead) 主要混合在土壤和灰尘中，能够被仪器检测到，所以可以找到权威环保专业机构提供的数据； 如果四乙铅混合中性土壤 (neutral Ph soil) 则危害性会减弱一些。
综上所述，化学理论可以让我们做出一个合理预测，即城市的土壤性质如果越趋向中性，则犯罪率理论上应该会相应降低。所以该工具变量与高危犯罪率之间存在着明显的非线性关系，可能会在某一个区间内会呈现更明显的相关性。
我们来看一下 “郊区化” 的衡量指标和其 OLS 模型：
其中， 代表就是郊区化程度，即 市 时市区中心 (city center，缩写为 cc) 人口除以市中心和近郊人口 ()。在 OLS 模型中，自变量为平均犯罪率 (violent crimes per capita)，即为 。此外，模型中还加入了地区固定变量 ，时间固定变量 ，以及二者的交互项 。
然而，由于模型中的平均犯罪率和郊区化程度存在互为因果和遗漏变量偏误等问题，所以使用四乙铅做为工具变量。土壤中的铅含量与郊区化本身没有关系 (1960 年代的人也没有这个意识和相关科学知识，不会有自选择偏误)，但是含铅量与犯罪率密切相关，并且含铅量本身是一个连续型变量。所以，使用工具变量后的第一阶段估计模型为：
上述式子表示，城市犯罪率 与历史上的土壤铅含量 (t-19) 相关。并且，如果土壤 Ph 值在 6.8-7.7 之间，则取 1，其他取 0，即在该 Ph 值范围内，待估计系数应该为负数。因为按照前述理论，中性土壤能有效降低对人体的危害，进而减少犯罪率。又由于土壤的 Ph 值因城市位置和时间所决定，故在第一阶段估计中加入时间效应和城市固定效应。
对模型有了一定了解后，接下来进入核心代码编写时间。其中关键便是工具变量转换成虚拟变量的技术处理。作者使用 discretize 命令进行第一阶段估计的前期处理。Ph 值作为工具变量，取 6.5-8 的范围内进行测试 (最小为 5，最大为 10，间隔为 1，即 Ph 值 0.1 ，进行反复测试)，加入各种内外生控制变量后，结果如下所示。我们可以看到 r1 为 F-statistics 最大值的区域，即 Ph 值上限为 7.7，下限为 6.8 的范围内。由此，我们可以确定 IV 进行分界的标准，即 Ph 值 6.8-7.7 区间。
基于上述测试结果，建立好土壤 (good_soil) 变量 (该地区的生物可容性的代理变量) —— Ph 值在 6.8 到 7.7 区间。然后使用 xtivreg 命令 (针对面板数据的工具变量 2SLS) 进行模型估计。其中，自变量为标准化后的犯罪率 (standardized_vc)，而工具变量为好土壤*铅含量 (good soil * Lead)，于是可以看到系数估计为 -0.072 左右 (原文图表中系数经过一系列固定效应调整后，稳定在 -0.0018) 。
然后，我们再进行第二阶段分析。在重复第一阶段命令基础上，加入因变量市中心人口密度 (perc_cc)。由下图可以看出，在之前的界限内估计系数的绝对值最大。
作者还用 graph 命令将第一阶段 IV 转换的数据可视化。我们可以清晰地观察到原始 IV 进行转换的上下界，其泡泡图最大值 (即 F 统计量) 确实是在下界 (lower bound) 68，和上界（upper bound) 77 的位置。而第二阶段在这个范围内的系数估计也是最大的蓝色气泡之一 (其他较大蓝色气泡点存在更大的标准误)，且周围系数波动小，所以较为平均。从敏感性分析的角度而言也证实了系数的稳定性。由此, 原本对于第二阶段系数敏感性过高的担忧也可以烟消云散了。
4. 参考文献
Fontenay S. discretize: Command to convert a continuous instrument into a dummy variable for instrumental-variable estimation[C]//London Stata Conference 2019. Stata Users Group, 2019 (20). -PDF- Curci F, Masera F. Flight from urban blight: lead poisoning, crime and suburbanization[J]. 2018. -PDF-
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 工具变量, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Lasso一下：再多的控制变量和工具变量我也不怕-T217 Stata新命令-pdslasso：众多控制变量和工具变量如何挑选？ 专题：IV-GMM Stata：当工具变量小于内生变量时，该如何估计？-mmeiv 找不到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