“众里寻他千百度”——工具变量法（IV）的Stata操作

工具变量法是经济学家的神奇创造，是一味解决内生性问题的良药，但事实上，想要找到一个合适的工具变量是相当困难的，工具变量法不需要你有多么高超的计量分析技术，它需要的是丰富的知识结构和突发奇想的艺术创造力。经历王国维先生所说“悬思”和“苦索”两境之后，你才有可能达到“顿悟”之境（众里寻他千百度，蓦然回首，那人却在，灯火阑珊处），找到“她”的身影。之前跟大家介绍过一些经典的、有趣的工具变量以及找寻工具变量的思路，今天咱们来瞅瞅工具变量法的Stata操作，让我们看一个经典的案例！

数据说明

陈婷、龚启圣和马驰骋三位老师（2020）发表在《The Economic Journal》杂志上的论文《Long Live Keju ! The Persistent Effects of China’s Civil Examination System》是一篇妙趣横生的工具变量论文，EJ杂志官网上公布了这篇论文使用的数据和代码，接下来我就使用公布的数据和代码跟大家分享一下工具变量法的Stata操作。

Replication Data for: Ting Chen, James Kai-sing Kung, Chicheng Ma.  Long Live Keju ! The Persistent Effects of China’s Civil Examination System[J]. The Economic Journal,  2020,  130 (631) : 2030–2064.

建议大家在看下面的内容之前，最好先看一下“科举万岁！妙哉工具变量也！”这篇推文，这样可能理解起来更加顺畅。

识别策略

在这篇论文中，作者研究的问题是明清时期各府在科举考试中的成功率（以明清时期的进士人数衡量）与当代人力资本（以今天的受教育年限衡量）之间的关系。但是囿于遗漏变量问题带来的内生性问题，我们需要找到一个科举教育的工具变量。

作者将目光投向了书籍印刷：四书五经是中国科举考试的关键，然而光有教材是不够的，想要成为科举考试的赢家，你还需要一大堆参考书（类似于今天的《五年高考三年模拟》）。在某种程度上，各府在科举考试中的成功与印刷和获取书籍的便捷性密切相关，而主要印刷中心位于松竹产地附近（中国的印刷技术主要依靠松木和竹子来生产油墨和纸张），印刷所需的原料主要通过水路运输。因此，各府到最近的松木和竹子产地的河流距离是科举教育的一个可行的工具变量。

被解释变量：各府今天平均受教育年限对数lneduyear

核心解释变量（内生变量）：各府每万人的进士数量（进士密度）对数lnjinshipop

工具变量：各府到最近的松木和竹子产地的河流距离bprvdist

控制变量：各府夜间灯光亮度lnnightlight、到海岸线的距离lncoastdist、地形崎岖程度tri、农业产量suitability 、人口密度lnpopdensity、城市化率urbanrates和省级固定效应

Stata操作

工具变量法的难点在于找到一个合适的工具变量并说明其合理性，Stata操作其实相当简单，只需一行命令就可以搞定，我们通常使用的工具变量法的Stata命令主要就是ivregress命令和ivreg2命令。

ivregress命令

ivregress命令是Stata自带的命令，支持两阶段最小二乘（2SLS）、广义矩估计（GMM）和有限信息最大似然估计（LIML）三种工具变量估计方法，我们最常使用的是两阶段最小二乘法（2SLS），因为2SLS最能体现工具变量的实质，并且在球形扰动项的情况下，2SLS是最有效率的工具变量法。

顾名思义，两阶段最小二乘法（2SLS）需要做两个回归：

（1）第一阶段回归：用内生解释变量对工具变量和控制变量回归，得到拟合值。

（2）第二阶段回归：用被解释变量对第一阶段回归的拟合值和控制变量进行回归。

如果要使用2SLS方法，我们只需在ivregress后面加上2sls即可，然后将内生解释变量lnjinshipop和工具变量bprvdist放在一个小括号中，用=号连接。选项first表示报告第一阶段回归结果，选项cluster()表示使用聚类稳健的标准误。

ivregress 2sls lneduyear (lnjinshipop=bprvdist) lnnightlight lncoastdist tri suitability lnpopdensity urbanrates i.provid , first cluster(provid)

第一阶段回归结果

First-stage regressions
-----------------------

                                                Number of obs     =        274
                                                No. of clusters   =         28
                                                F(   7,    239)   =      85.27
                                                Prob > F          =     0.0000
                                                R-squared         =     0.6487
                                                Adj R-squared     =     0.5988
                                                Root MSE          =     0.4442

------------------------------------------------------------------------------
             |               Robust
 lnjinshipop |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
lnnightlight |    .183385   .0682506     2.69   0.008     .0489354    .3178346
 lncoastdist |   .0350333    .077158     0.45   0.650    -.1169634    .1870299
         tri |    1.06676   .5637082     1.89   0.060    -.0437105    2.177231
 suitability |  -.0769726   .0549697    -1.40   0.163    -.1852596    .0313144
lnpopdensity |    .196144   .0843727     2.32   0.021     .0299349    .3623532
  urbanrates |   3.352916   1.687109     1.99   0.048      .029414    6.676419
             |
      provid |
         12  |   .2051006   .0551604     3.72   0.000      .096438    .3137632
         13  |  -1.890425   .0951146   -19.88   0.000    -2.077795   -1.703055
......
         64  |  -1.301895   .1581021    -8.23   0.000    -1.613346   -.9904433
             |
    bprvdist |  -.0846917   .0107859    -7.85   0.000    -.1059393   -.0634441
       _cons |   2.126233   .9791046     2.17   0.031     .1974567     4.05501
------------------------------------------------------------------------------

从表中可以看出，工具变量bprvdist的系数为-0.085，标准误为0.011，在1%的水平上显著。

第二阶段回归结果

Instrumental variables (2SLS) regression          Number of obs   =        274
                                                  Wald chi2(34)   =      13.62
                                                  Prob > chi2     =     0.9993
                                                  R-squared       =     0.7874
                                                  Root MSE        =     .05434

                                (Std. Err. adjusted for 28 clusters in provid)
------------------------------------------------------------------------------
             |               Robust
   lneduyear |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 lnjinshipop |   .0834221   .0116773     7.14   0.000      .060535    .1063092
lnnightlight |   .0592777   .0102138     5.80   0.000     .0392589    .0792964
 lncoastdist |   .0093399   .0118791     0.79   0.432    -.0139426    .0326224
         tri |  -.1044346    .063917    -1.63   0.102    -.2297097    .0208404
 suitability |  -.0008944   .0125824    -0.07   0.943    -.0255555    .0237666
lnpopdensity |  -.0472648   .0115102    -4.11   0.000    -.0698245   -.0247051
  urbanrates |  -.0353689   .1784687    -0.20   0.843    -.3851612    .3144233
             |
      provid |
         12  |    -.10732   .0088755   -12.09   0.000    -.1247157   -.0899243
         13  |   .0244675   .0276005     0.89   0.375    -.0296284    .0785634
......
         64  |  -.0740796   .0350354    -2.11   0.034    -.1427477   -.0054116
             |
       _cons |   2.014146   .1490819    13.51   0.000     1.721951    2.306341
------------------------------------------------------------------------------

第二阶段回归结果就是使用工具变量法对内生性问题进行“治疗”之后得到的估计结果，很多论文一般也只报告第二阶段回归结果。从表中可以看出，进士密度lnjinshipop的系数为0.083，标准误为0.012，在1%的水平上显著。

弱工具变量检验

检验工具变量的有效性，首先要检验工具变量的相关性。存在弱工具变量问题时，2SLS 估计不仅难以矫正OLS 估计的偏差，反因有更大的标准差而有更低的效率， 导致“治疗比疾病本身更坏”。我们可以使用estat firststage命令对第一阶段的结果进行分析，以判断是否存在弱工具变量问题。

. estat firststage,forcenonrobust

  First-stage regression summary statistics
  --------------------------------------------------------------------------
               |            Adjusted      Partial       Robust
      Variable |   R-sq.       R-sq.        R-sq.       F(1,27)   Prob > F
  -------------+------------------------------------------------------------
   lnjinshipop |  0.6487      0.5988       0.3276        61.914    0.0000
  --------------------------------------------------------------------------
  (F statistic adjusted for 28 clusters in provid)


  Minimum eigenvalue statistic = 116.419     

  Critical Values                      # of endogenous regressors:    1
  Ho: Instruments are weak             # of excluded instruments:     1
  ---------------------------------------------------------------------
                                     |    5%     10%     20%     30%
  2SLS relative bias                 |         (not available)
  -----------------------------------+---------------------------------
                                     |   10%     15%     20%     25%
  2SLS Size of nominal 5% Wald test  |  16.38    8.96    6.66    5.53
  LIML Size of nominal 5% Wald test  |  16.38    8.96    6.66    5.53
  ---------------------------------------------------------------------

结果显示，偏=0.3276（偏是扣除了其他外生变量后工具变量的解释力度），说明工具变量bprvdist对内生变量lnjinshipop有很强的解释力度。F统计量=61.914>10，根据经验准则可以判断，我们的工具变量不是一个弱工具变量。除此之外，estat firststage命令还会报告一个最小特征值统计量（Minimum eigenvalue statistic），一般大于2SLS Size of nominal 5% Wald test中10%对应的临界值就表明不存在弱工具变量问题。

特别说明：豪斯曼检验（即内生性检验）和过度识别检验（即外生性检验）在此就不予以介绍了，因为这两个检验都是很“鸡肋”的，并且这里只有一个工具变量，也没法做过度识别检验，后面我将专门写一篇推文对这两个检验进行详细说明。

ivreg2命令

ivreg2命令是对ivregress命令的改进和优化，功能更加强大，支持的估计方法更多（默认使用2SLS），并且会直接报告工具变量的几个统计检验结果，陈婷、龚启圣和马驰骋三位老师使用的就是ivreg2命令。

ivreg2命令是一个外部命令，所以使用之前需要安装（ssc install ivreg2）。它的语法格式和ivregress基本一致：

ivreg2 lneduyear (lnjinshipop=bprvdist) lnnightlight lncoastdist tri suitability lnpopdensity urbanrates i.provid , first cluster(provid)

第一阶段回归结果

First-stage regression of lnjinshipop:

Statistics robust to heteroskedasticity and clustering on provid
Number of obs =                    274
Number of clusters (provid) =       28
------------------------------------------------------------------------------
             |               Robust
 lnjinshipop |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    bprvdist |  -.0846917   .0107633    -7.87   0.000    -.1058948   -.0634886
lnnightlight |    .183385   .0681077     2.69   0.008     .0492169    .3175531
 lncoastdist |   .0350333   .0769964     0.45   0.650     -.116645    .1867116
         tri |    1.06676   .5625276     1.90   0.059    -.0413849    2.174906
 suitability |  -.0769726   .0548546    -1.40   0.162    -.1850328    .0310876
lnpopdensity |    .196144    .084196     2.33   0.021      .030283    .3620051
  urbanrates |   3.352916   1.683576     1.99   0.048     .0363742    6.669459
             |
      provid |
         12  |   .2051006   .0550449     3.73   0.000     .0966656    .3135357
         13  |  -1.890425   .0949154   -19.92   0.000    -2.077402   -1.703447
......
         64  |  -1.301895    .157771    -8.25   0.000    -1.612694   -.9910956
             |
       _cons |   2.126233   .9770541     2.18   0.031      .201496    4.050971
------------------------------------------------------------------------------

第二阶段回归结果

IV (2SLS) estimation
--------------------

Estimates efficient for homoskedasticity only
Statistics robust to heteroskedasticity and clustering on provid

Number of clusters (provid) =       28                Number of obs =      274
                                                      F( 34,    27) =     0.34
                                                      Prob > F      =   0.9984
Total (centered) SS     =  3.805186838                Centered R2   =   0.7874
Total (uncentered) SS   =  1282.736705                Uncentered R2 =   0.9994
Residual SS             =  .8089841945                Root MSE      =   .05434

------------------------------------------------------------------------------
             |               Robust
   lneduyear |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 lnjinshipop |   .0834221   .0116773     7.14   0.000      .060535    .1063092
lnnightlight |   .0592777   .0102138     5.80   0.000     .0392589    .0792964
 lncoastdist |   .0093399   .0118791     0.79   0.432    -.0139426    .0326224
         tri |  -.1044346    .063917    -1.63   0.102    -.2297097    .0208404
 suitability |  -.0008944   .0125824    -0.07   0.943    -.0255555    .0237666
lnpopdensity |  -.0472648   .0115102    -4.11   0.000    -.0698245   -.0247051
  urbanrates |  -.0353689   .1784687    -0.20   0.843    -.3851612    .3144233
             |
      provid |
         12  |    -.10732   .0088755   -12.09   0.000    -.1247157   -.0899243
         13  |   .0244675   .0276005     0.89   0.375    -.0296284    .0785634
......
         64  |  -.0740796   .0350354    -2.11   0.034    -.1427477   -.0054116
             |
       _cons |   2.014146   .1490819    13.51   0.000     1.721951    2.306341
------------------------------------------------------------------------------

弱工具变量检验结果

ivreg2命令会直接报告Cragg-Donald Wald F 统计量和Kleibergen-Paap Wald rk F统计量两个用于弱工具变量检验的统计量，其中Cragg-Donald Wald F 统计量假设扰动项iid（独立同分布），而Kleibergen-Paap Wald rk F统计量没有做扰动项iid的假设。需要注意的是，这里的原假设是所选工具变量是弱工具变量，当两个F统计量大于Stock-Yogo weak ID test critical values中10%偏误的临界值时，可以拒绝原假设，认为不存在弱工具变量问题。

Weak identification test
Ho: equation is weakly identified
Cragg-Donald Wald F statistic                                     116.42
Kleibergen-Paap Wald rk F statistic                                61.91

Stock-Yogo weak ID test critical values for K1=1 and L1=1:
                                   10% maximal IV size             16.38
                                   15% maximal IV size              8.96
                                   20% maximal IV size              6.66
                                   25% maximal IV size              5.53
Source: Stock-Yogo (2005).  Reproduced by permission.
NB: Critical values are for Cragg-Donald F statistic and i.i.d. errors.

需要本篇推文使用的数据和代码的朋友，请在后台对话框内回复关键词“IV”。

补充：对于面板数据的工具变量法操作，Stata还有专门的两个命令——xtivreg和xtivreg2两个命令，语法格式与ivreg2基本类似，大家有兴趣的话可以去看一下这两个命令的help文件。