Stata：最牛IV-Shift Share IV 实操-ssaggregate

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

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

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
  . ssc install lianxh
详情参见帮助文件 (有惊喜)：
  . help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg
 

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

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

⛳ Stata 系列推文：

• 全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
• 计量专题 | 论文写作 | 数据分享 | 专题课程
• 结果输出 | Stata绘图 | 数据处理 | Stata程序
• 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
• 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
• Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
• 文本分析-爬虫 | Python-R-Matlab | 机器学习
• Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：甘徐沁 (复旦大学)
邮箱：jesse_gan@foxmail.com

目录

1. 简介

在之前推文中，我们简单介绍了 Bartik 工具变量的缘起，以及基于其最新发展出的 Shift Share IV ( Bartik 工具变量和 Shift Share IV 仅是两个不同的名字，指代同一事物，但理解角度不同)。在本文中，我们来介绍一下基于 Shift Share IV 等价理论编写的 Stata 命令：ssaggregate。

2. Shift Share IV 等价理论

为了更加清晰地理解 Shift Share IV 背后的原理，我们先来回顾一下 Borusyak 等 (2020) 提出的 Bartik 工具的 Shift Share IV 等价理论。

考虑研究主题是劳动力的供给水平和工资之间的关系。主要被解释变量 $y_{\ell }$ 是某地区 $\ell$ 工资增长率，主要解释变量 $x_{\ell }$ 是地区 $\ell$ 的劳动力就业增长率。在下面的估计中，感兴趣的因果效应系数是 $\beta$ ，即逆劳动供给弹性：

其中 $\epsilon$ 和控制变量 $w_{\ell }$ 正交。考虑到同一地区内影响劳动力工资水平变化的因素和影响劳动力就业的因素，可能由多种因素共同决定：比如移民流入情况，宏观经济走势，人力资本变化等，从而导致 OLS 的估计不一致，因此我们需要考虑使用工具变量。Bartik 工具变量被定义为：

其中 $g_n$ 是一系列产业冲击，$s_{\ell n}$ 是地区 $\ell$ 在初始时期产业 $n$ 的就业份额，或称冲击份额。产业 $n$ 的冲击对于不同地区的作用取决于此地区产业 $n$ 的初始份额。假设每个地区的冲击份额之和为常数，即 $\sum _n{s}_{\ell n}=1$。同时假设冲击是随机变量，并且数据非 iid 生成，因此不同地区所受的冲击可能是相关的，并且不同地区的不可观测的冲击也可能相关 (误差项自相关)。

当我们假设数据并非 iid 生成时，Bartik 工具变量回归的矩条件可以写成：

这个条件描述了 Bartik 工具和第二阶段误差项之间的正交关系，其中 $e_{\ell }$ 为地区 $\ell$ 的劳动力在全国市场中所占份额。根据 Frisch-Waugh-Lovell 定理，我们将 Bartik 工具变量回归系数写为残差化的 $y$ 和 $x$ 的二元关系：

其中 $y_{\ell }^{\perp }$ 和 $x_{\ell }^{\perp }$ 分别表示将 $y$ 和 $x$ 对控制变量 $w$ 做回归 (以 $e_{\ell }$ 作为权重) 所得到的残差。

在产业 (冲击) 层面的正交等价条件可以表达为：

其中，$s_n=\sum _{\ell }{e}_{\ell }{s}_{\ell n}$ 可以解读为某产业在总经济中的所占比例，${\overline{\epsilon }}_n=\frac{\sum _{\ell }{e}_{\ell }{s}_{\ell n}{\epsilon }_{\ell }}{\sum _{\ell }{e}_{\ell }{s}_{\ell n}}$ 则是冲击-权重的误差项平均。这个等式的左边描述了 Bartik 工具 $z_{\ell }$ 和误差项 ${\epsilon }_{\ell }$ 的正交关系，以 $e_{\ell }$ 为权重，右边则描述了冲击项 $g_n$ 和 ${\overline{\epsilon }}_n$ 的正交关系，以 $s_n$ 为权重。

由此可以看出，右边的条件成立时，Bartik 工具的外生性成立。上式中的等价关系为我们引出 Shift Share IV 估计等价理论提供了灵感。接着，使用 Bartik 工具变量对原始方程进行回归所得的系数 $\hat{\beta }$、以 $s_n$ 为权重的冲击层面、和以 $g_n$ 为工具对以下方程进行工具变量回归等价：

其中 ${\overline{v}}_n=\frac{\sum _{\ell }{e}_{\ell }{s}_{\ell n}{v}_{\ell }}{\sum _{\ell }{e}_{\ell }{s}_{\ell n}}$ 表示冲击-权重的变量 $v_{\ell }$ 的均值。

证明：

上述证明中最后边的等式即以 $s_n$ 为权重的冲击层面，以 $g_n$ 为工具的工具变量回归。

由此我们可以概括性地将 Bartik 工具变量的估计等价关系记忆为：地区层面的工具变量回归，等价于将变量转换为冲击 (或称产业) 层面后的工具变量回归。前者的估计单位，或称每个观测值为地区，而后者中，每个观测值为产业。当然需要注意，在这里两个工具变量回归使用的权重并不相同，前者使用的权重是地区的劳动力份额，而后者使用的权重为 $s_n=\sum _{\ell }{e}_{\ell }{s}_{\ell n}$，可以理解为各冲击在不同地区的 (劳动力份额) 加权平均。

当然 Shift Share IV 的理论建立并不仅仅包括了此估计等价理论，更为重要的是我们可以通过将地区层面的 Bartik 工具变量估计视为等同于冲击层面的 Shift Share IV 估计，于是 “随机分配” 的冲击便可以使得 Shift Share IV ——以及等价的原 Bartik 工具变量，得到真实 $\beta$ 的一致估计。

3. ssaggregate 命令

既然 Shift Share IV 的实施需要先将数据从产业层面转换为冲击层面， Shift Share IV 的几位作者便编写了 ssaggregate 命令来满足这一前提需求。这个变量实际上是先将地区层面的被解释变量和解释变量对控制变量残差化，再根据冲击比例在冲击层面上进行加权平均，即 ${\overline{v}}_n=\frac{\sum _{\ell }{e}_{\ell }{s}_{\ell n}{v}_{\ell }}{\sum _{\ell }{e}_{\ell }{s}_{\ell n}}$。

*命令安装
ssc install ssaggregate, replace

*命令语法
Using "long" exposure weights, saved separately (executes faster):

	ssaggregate varlist [if] [in] [weight], n(varlist) s(varname) l(varlist) sfilename(filename) [other options]
                 
Using "wide" exposure weights, saved in memory:

        ssaggregate varlist [if] [in] [weight], n(varname) s(stubname) [other options]

ssaggregate 提供了如下的选项：

• t(varlist)：时期识别变量；
• controls(strings)：控制变量组，用于 partialling out 获得残差，可设定多组。如 controls("var1" "var1 var2" "var1 var2 var3")；
• absorb(strings)：固定效应变量组，用于吸收固定效应，可设定多组。如 absorb("f1" "fe1 fe2")；
• addmissing：生成 “缺失产业” 变量。

如果产业权重是 "long" 格式，那么必须设定如下四个选项：

• n(varlist)：产业识别变量；
• sfilename："long" 数据格式的产业权重数据文件；
• s(varname)：产业权重变量；
• l(varlist)：地区识别符。

同时，我们可以使用 Stata 常规的加权选项，将地区劳动力份额作为权重包括进来。由于作者提示若将产业权重数据以 "long" 格式储存在额外的一份数据集中，会使得命令执行速度更快，故我们在本文中仅演示 "long" 格式的命令。

4. Shift Share IV 应用

在 Autor 等 (2013) (以下简称 ADH) 的文章中，三位作者研究了美国 “日益增长的中国进口产品” 这一事实如何影响了美国劳动力市场。根据贸易理论，在 1990-2010 年间，中国因为转型市场经济体带来的生产率上升，以及加入世界贸易组织带来的贸易成本下降，大大增强了中国制造业产品的国际竞争力。

在这 20 年间，美国对于中国产品的进口大大增加，这可能会使得美国对于本土制造业的需求下降，从而对美国本土劳动力市场产生一系列负面效果，如制造业劳动力就业率下降，劳动力从制造业转移到其他产业等。

ADH 的文章中，作者关注的变量如下：

• 被解释变量 $y_{\ell t}$，是 $t$ 时期地区 $\ell$ 的制造业就业人数占总体工作年龄人口的比例的增加值；
• 解释变量 $x_{\ell t}$，是每名美国工人水平的地区 $\ell$ 的进口冲击大小。值得注意的是，这里的解释变量本身也是一个内积的形式，表达为 $x_{\ell t}=\sum _n\frac{L_{\ell nt}}{L_{nt}}\frac{\mathrm{\Delta }{M}_{ucnt}}{L_{\ell t}}$。其中 $L$ 为劳动力雇佣数量，$\ell$ 和 $n$ 分别表示地区与产业，$\mathrm{\Delta }{M}_{uc}$ 表示两期内 (1990-2000，2000-2007) 美国对中国的产品进口额对应增加值；
• 工具变量为 $z_{\ell t}=\sum _n{s}_{\ell nt}{g}_{nt}$。其中 $s_{\ell nt}$ 是地区 $\ell$ 产业 $n$ 的就业份额，使用两期分别十年前的数据来衡量。$g_{nt}$ 是同期世界上另外八个高收入国家在产业 $n$ 对中国的产品进口增加值 (每名美国工人水平) 之和，用于衡量因为中国的产品竞争力提高而对美国的各产业造成的冲击。

此处使用的数据来源于 Shift Share Share IV 作者 Borusyak 的 Github 仓库，链接为：https://github.com/borusyak/shift-share。我们首先使用地区层面的回归来复制 ADH 原文表 3 的第 6 列，回归为两期混合截面。在下面的回归中：t2 表示第二期，l_shind_manuf_cbp 为期初的地区制造业就业份额，reg* 为一系列普查区域类别。

l_sh_popedu_c、l_sh_popfborn、l_sh_empl_f、l_sh_routine33、l_task_outsource 为5个控制变量，分别是大学教育人群比例、外国出生人群比例、女性就业比例、程序化工作比例、工作离岸外包指数。 wei 实际上是 $e_l$，即各地区劳动力比例。

. lxhuse location_level.dta, clear
. ivreg2 y (x=z) t2 l_shind_manuf_cbp reg* l_sh_popedu_c l_sh_popfborn ///
>   l_sh_empl_f l_sh_routine33 l_task_outsource  [aw=wei], cluster(clus)
(sum of wgt is     2.0000e+00)

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

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

Number of clusters (clus) =         48                Number of obs =     1444
                                                      F( 16,    47) =    42.53
                                                      Prob > F      =   0.0000
Total (centered) SS     =  4396.587068                Centered R2   =   0.3429
Total (uncentered) SS   =   12720.4953                Uncentered R2 =   0.7729
Residual SS             =  2889.068933                Root MSE      =    1.414

-----------------------------------------------------------------------------------
                  |               Robust
                y | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
                x |     -0.596      0.099    -6.04   0.000       -0.790      -0.403
               t2 |     -0.242      0.401    -0.61   0.545       -1.028       0.543
l_shind_manuf_cbp |     -4.022      1.314    -3.06   0.002       -6.598      -1.446
       reg_midatl |      0.313      0.281     1.11   0.266       -0.238       0.864
        reg_encen |      1.261      0.337     3.74   0.000        0.600       1.921
        reg_wncen |      1.624      0.372     4.36   0.000        0.894       2.353
         reg_satl |     -0.288      0.234    -1.23   0.217       -0.746       0.170
        reg_escen |      1.076      0.335     3.22   0.001        0.420       1.732
        reg_wscen |      0.732      0.231     3.16   0.002        0.278       1.185
        reg_mount |      0.402      0.257     1.56   0.118       -0.102       0.906
        reg_pacif |      0.027      0.191     0.14   0.888       -0.348       0.401
    l_sh_popedu_c |      0.013      0.012     1.08   0.281       -0.011       0.037
    l_sh_popfborn |      0.030      0.011     2.80   0.005        0.009       0.052
      l_sh_empl_f |     -0.006      0.024    -0.24   0.811       -0.054       0.042
   l_sh_routine33 |     -0.245      0.064    -3.84   0.000       -0.370      -0.120
 l_task_outsource |     -0.059      0.237    -0.25   0.803       -0.523       0.405
            _cons |      6.279      1.937     3.24   0.001        2.482      10.075
-----------------------------------------------------------------------------------
Underidentification test (Kleibergen-Paap rk LM statistic):             16.880
                                                   Chi-sq(1) P-val =    0.0000
------------------------------------------------------------------------------
Weak identification test (Cragg-Donald Wald F statistic):              533.322
                         (Kleibergen-Paap rk Wald F statistic):         47.643
Stock-Yogo weak ID test critical values: 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.
------------------------------------------------------------------------------
Hansen J statistic (overidentification test of all instruments):         0.000
                                                 (equation exactly identified)
------------------------------------------------------------------------------
Instrumented:         x
Included instruments: t2 l_shind_manuf_cbp reg_midatl reg_encen reg_wncen
                      reg_satl reg_escen reg_wscen reg_mount reg_pacif
                      l_sh_popedu_c l_sh_popfborn l_sh_empl_f l_sh_routine33
                      l_task_outsource
Excluded instruments: z
------------------------------------------------------------------------------

下面，我们使用 ssaggregate 将原始数据转换为冲击，即产业层面。其中，sic87dd 是四位数的产业代码，Lshares 是外部分年的地区各产业比例数据集，czone 是地区指示变量，ind_share 是产业比例变量。

. lxhget adh.zip, replace
. unzipfile adh.zip, replace
. use ./adh/location_level.dta, clear

. ssaggregate y z x  [aw=wei], n(sic87dd) t(year) s(ind_share) sfilename(./adh/Lshares) l(czone) addmissing ///
>   controls("t2 l_shind_manuf_cbp reg* l_sh_popedu_c l_sh_popfborn l_sh_empl_f l_sh_routine33 l_task_outsource")
. replace sic87dd = 0 if mi(sic87dd)
. merge 1:1 sic87dd year using ./adh/shocks, assert(1 3) nogen  // sic87dd==0 is _m==1
. gen sic3 = int(sic87dd/10)
. foreach v of varlist g* sic3{
  2.   replace `v'= 0 if sic87dd==0
  3. }
. label var g "Industry China Shock (binned)"

. *冲击层面的 Shift Share IV 回归
. ivreg2 y (x=g) [aw=s_n], cluster(sic3)

IV (2SLS) estimation
--------------------
Estimates efficient for homoskedasticity only
Statistics robust to heteroskedasticity and clustering on sic3
Number of clusters (sic3) =        137                Number of obs =      796
                                                      F(  1,   136) =    27.12
                                                      Prob > F      =   0.0000
Total (centered) SS     =  44.66656614                Centered R2   =   0.2111
Total (uncentered) SS   =  44.66656614                Uncentered R2 =   0.2111
Residual SS             =  35.23809809                Root MSE      =    .2104
------------------------------------------------------------------------------
             |               Robust
           y | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
           x |     -0.596      0.114    -5.23   0.000       -0.820      -0.373
       _cons |      0.000      0.011     0.00   1.000       -0.021       0.021
------------------------------------------------------------------------------
Underidentification test (Kleibergen-Paap rk LM statistic):              8.285
                                                   Chi-sq(1) P-val =    0.0040
------------------------------------------------------------------------------
Weak identification test (Cragg-Donald Wald F statistic):              122.362
                         (Kleibergen-Paap rk Wald F statistic):         39.649
Stock-Yogo weak ID test critical values: 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.
------------------------------------------------------------------------------
Hansen J statistic (overidentification test of all instruments):         0.000
                                                 (equation exactly identified)
------------------------------------------------------------------------------
Instrumented:         x
Excluded instruments: g
------------------------------------------------------------------------------

可以看出，同估计等价的理论相一致，在地区层面的工具变量回归和转换到冲击层面的工具变量回归的结果完全相同，都是 -0.596。但是标准误并不相同，这是因为数据转换后样本量产生了变化 (同时聚类的类定义也发生了改变)。感兴趣并想要深入探究的读者可以阅读 Borusyak 等 (2020) 文章的第五节 "Shock-Level Inference and Testing"。

5. 结论

在本文中，我们回顾了 Bartik 工具变量的估计等价定理，即地区层面的工具变量回归，等价于将变量转换为冲击 (或称产业) 层面后的工具变量回归。同时我们还介绍了将数据层级进行转换的 ssaggregate 命令。我们可以通过此命令将数据从地区层面转换到冲击层面，再使用诸如外部命令 ivreg2 或 Stata 提供的内置命令 ivregress 进行工具变量估计。

Bartik 工具变量的理论研究在目前还是一个关于因果推断的新兴话题，如果有感兴趣的读者，可以关注更多的此方面的文献。

6. 参考文献

• Borusyak, Kirill, Peter Hull, and Xavier Jaravel. “Quasi-Experimental Shift-Share Research Designs.” The Review of Economic Studies, June 8, 2021, rdab030. -Link-
• Autor, David H, David Dorn, and Gordon H Hanson. “The China Syndrome: Local Labor Market Effects of Import Competition in the United States.” American Economic Review 103, no. 6 (October 1, 2013): 2121–68. -Link-

7. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh 工具变量, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：Stata命令

  • Stata新命令-pdslasso：众多控制变量和工具变量如何挑选？
• 专题：IV-GMM

  • Stata：无需工具变量的IV估计-kinkyreg-
  • Stata：当工具变量小于内生变量时，该如何估计？-mmeiv
  • Lasso一下：再多的控制变量和工具变量我也不怕-T217
  • IV在哪里？奇思妙想的工具变量
  • twostepweakiv：弱工具变量有多弱？
  • 多个(弱)工具变量如何应对-IV-mivreg？
  • IV：工具变量不满足外生性怎么办？
  • IV-工具变量法：第一阶段系数符号确定时的小样本无偏估计
  • IV：可以用内生变量的滞后项做工具变量吗？
  • Stata: 工具变量法 (IV) 也不难呀！
  • IV-估计：工具变量不外生时也可以用！
• 专题：内生性-因果推断

  • Stata：内生变量与工具变量非线性关系处理-discretize
  • 工具变量-IV：排他性约束及经典文献解读

相关课程

免费公开课

• 直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页，Bilibili 站版
• Stata 33 讲 - 连玉君, 每讲 15 分钟. 课程主页，课件，Bilibili 站版
• Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页
• 部分直播课 课程资料下载 (PPT，dofiles等)

最新课程-直播课

• Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

⛳ 课程主页

⛳ 课程主页

关于我们

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

连享会小程序：扫一扫，看推文，看视频……

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

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

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