Stata:当工具变量小于内生变量时,该如何估计?-mmeiv

发布时间:2021-09-18 阅读 3757

Stata连享会   主页 || 视频 || 推文 || 知乎

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装命令如下:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

作者:张子楠 (浙江财经大学)
邮箱zinanzh@gmail.com


目录


1. 背景简介

工具变量法的可识别条件要求工具变量 (IV) 个数大于等于内生变量的个数,即所谓的工具变量秩条件。然而,在实际研究中,我们也可能碰到工具变量个数小于内生变量个数的情形。针对这种情况,Caetano and Escanciano (2021) 提出了一种思路,即通过使用工具变量与控制变量交叉项来构建新的工具变量,以满足秩条件要求。

下面以两个内生变量 (X1 和 X2) 和一个二元工具变量 Z 情形来进行说明。回归方程如下所示:

结合模型和内生性情况的设定,我们设定数据生成过程 (DGP) 如下所示:

其中 W 为控制变量,T 为遗漏变量,U1U2UWuTu 为随机扰动项,且均服从 N(0,1) 的标准正态分布。观察 DGP 可知,由于 X1X2 和 U 都包含 W 和 T,存在由于遗漏变量导致的内生性问题。与此同时,变量 Z 和 X 相关,且 Z 和 U 无关,Z 只能通过 X1 和 X2 来影响 Y,从而 Z 满足工具变量的有效性和外生性要求。不妨设定 Z 服从均值为 0.5 的二项分布。

由于这里有两个内生变量,只有一个工具变量,所以传统方法上无法进行识别。为理解这一点,我们用二阶段最小二乘法 (TSLS) 进行说明。TSLS 的第一阶段是将每个内生变量单独对所有工具变量和所有外生变量进行回归:

显然这种方法下,工具变量矩阵的秩为 1,难以识别出系数 δ11 和 δ12。Caetano and Escanciano (2021) 思路则在于如果对于不同的 WZ 与内生变量关系出现显著波动的话。那么就可以生成交叉项ZW,并在第一阶段进行如下回归:

如果 δ11δ32δ12δ31,即随着不同的 WZ 对 X 的回归系数有显著不同,从而可以满足秩条件满足,我们就可以识别出回归系数。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),可选 robustcluster 等方式,默认是不控制。其它四个是有关回归结果的边际效应内容,和本文主题关系不大,在此不一一介绍。

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 

这样我们就生成了回归所需要的数据。其中被解释变量为 Y,内生变量为 X1 和 X2,控制变量为 W,工具变量为 Z,遗漏变量为 T,随机扰动项为 U。易知有 X1X2 和 U 相关,满足内生性要求。此外,Z 和 X1X2 相关,满足工具变量有效性要求;且 Z 和 U 不相关,满足外生性要求。

. 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 回归结果如下所示,可知 X1 的回归系数为 0.68,95% 的置信区间为 [0.42 0.94]。虽然和真实回归系数 β01=0.8 有所偏离,但置信区间包括了真实回归系数值。X2 的回归结果也相似,说明 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 回归;第二个强行认定只有一个内生变量,用一个工具变量 Z 进行 TSLS 回归;第三个认为是有两个内生变量,先生成 Z 和 W 的交乘项 ZW,然而用 ZW 和 Z 两个一起作为工具变量,进行 TSLS 回归。

第一种回归,也就是 OLS 回归的代码和结果如下所示。对于 X1 回归系数,OLS 回归值为 0.778,95% 置信区间为 [0.682 0.874]。相对于 mmeiv 回归结果,其置信区间同样包含了真实回归系数 0.8,同时区间宽度更窄,回归系数偏误也更小。这个结果似乎表示 OLS 可能更为有效,但事实上这很大程度上是偶然因素导致。本文将在第 4 节通过 Monte Carlo 模拟来证明这点。

此外,X2 回归结果为 1.172,偏离真实值 0.7,同时置信区间也没有包含真实值 0.7。OLS 回归对 X2 的识别有效性要显著弱于 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
------------------------------------------------------------------------------

第二种回归是强行认定只有一个内生变量,并且只用一个工具变量 Z 进行 TSLS 回归。回归代码和结果如下所示,其中 X1 的回归系数为 0.58,置信区间包括了真实值,但 X2 的回归系数则存在较大偏误,置信区间也没有包含真实回归系数值。

. 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

第三,用 Z 和 W 交乘项 ZW,以及 Z 共同作为工具变量来进行 IV 回归。回归代码和结果如下所示,观察回归结果可以发现,无论是回归系数、标准误、置信区间以及 R2 等,均与 mmeiv 回归结果完全一致。这并非偶然,事实上如第 1 节介绍的一样,mmeiv 命令的解决思路,就在于构建交乘项 ZW,然后再将 ZW 和 Z 同时作为工具变量,进行两阶段最小二乘法工具变量回归的。

. // 生成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:在利用协变量 W 生成交乘项时,W 要满足协变量完备性 (convaiance completeness),即 W 包含元素个数至少要比内生变量多一个。相关证明见 Caetano and Escanciano (2021)。

补充说明 2:在这里,我们并没有要求 W 是外生的。事实上从回归方程随机扰动项 U 的生成过程 U=γ0W+T+u 可知,W 并非外生的,但我们仍旧获得了 X1 和 X2 的一致回归结果。相对于再寻找一个好的 IV,mmeiv 命令只需要一个满足完备性的协变量,甚至允许该协变量可以不外生,从而大大降低了识别难度。

补充说明 3:上述例子可以很容易扩展 N 个内生变量但只有 1 个工具变量情形。只需要找到 N-1 个满足完备性的协变量,然后放入 { varlist_W1 } 地方再回归即可。

补充说明 4:Caetano and Escanciano (2021) 的思路也为一个常见的 IV 处理方法提供了理论支持。在计量应用文献里,有时候工具变量是不随时间而变动的变量,如小河流长度 (River) 等。如果这时再同时控制时间固定效应进行面板回归,工具变量会在组内去芯过程中被消除掉。对于这种问题,一个解决方法是将小河流长度交乘年份 (Year) 生成新的变量 RiverYear,然后用 RiverYear 作为新的工具变量进行工具变量回归。这种处理方法和 Caetano and Escanciano (2021) 的处理思路是完全一致的。进一步理解,我们不仅可以用 Year 来交乘,我们也可以用其他协变量来交乘。(注:读者可思考一下,这里是否还要求 Year 这类协变量满足协变量完备性?)

补充说明 5mmeiv 回归还可以通过绘制图形,来表示工具变量回归的边际效应,特别是工具变量维度小于内生变量维度的情形。比如工具变量是政策这种 0-1 变量,而内生变量是连续变量的情形。下面来展示相应用法,代码和边际效应图分别如下所示:

. mmeiv Y  (X1 X2=Z {W}), plot 

4 命令有效性检验

在第 3 节,mmeiv 回归结果中 X1 回归系数为 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 次模拟结果显示,X1 回归系数均值为 0.8,和真实系数值 0.80 一致。同时 95% 置信区间为 [0.792 0.809],包含了真实系数。另外,拒绝回归系数等于 0.8 的 p 值为 0.505,无法拒绝。类似地,对于 X2,回归系数均值为 0.702,同时 95% 置信区间为 [0.699 0.705],包含了真实系数 0.7。另外,拒绝回归系数等于 0.7 的 p 值为 0.505,无法拒绝。可见,第 3.2 节中的偏误很大程度是偶然因素所导致。从 MC 模拟结果可以看出,mmeiv 命令的估计量具有较好的有效性。

作为对比,我们也对 OLS 回归进行 MC 模拟,代码如下所示。MC 模拟结果显示,X1 回归系数均值为 0.797, 置信区间为 [0.793 0.800],X2 回归系数均值为 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

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,400+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxh 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh