Stata：最牛IV-Shift Share IV 实操-ssaggregate
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会寒假班
作者：甘徐沁 (复旦大学)
邮箱：jesse_gan@foxmail.com
目录
1. 简介
2. Shift Share IV 等价理论
3. ssaggregate 命令
4. Shift Share IV 应用
5. 结论
6. 参考文献
7. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
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 等价理论。
考虑研究主题是劳动力的供给水平和工资之间的关系。主要被解释变量 是某地区 工资增长率，主要解释变量 是地区 的劳动力就业增长率。在下面的估计中，感兴趣的因果效应系数是 ，即逆劳动供给弹性：
其中 和控制变量 正交。考虑到同一地区内影响劳动力工资水平变化的因素和影响劳动力就业的因素，可能由多种因素共同决定：比如移民流入情况，宏观经济走势，人力资本变化等，从而导致 OLS 的估计不一致，因此我们需要考虑使用工具变量。Bartik 工具变量被定义为：
其中 是一系列产业冲击， 是地区 在初始时期产业 的就业份额，或称冲击份额。产业 的冲击对于不同地区的作用取决于此地区产业 的初始份额。假设每个地区的冲击份额之和为常数，即 。同时假设冲击是随机变量，并且数据非 iid 生成，因此不同地区所受的冲击可能是相关的，并且不同地区的不可观测的冲击也可能相关 (误差项自相关)。
当我们假设数据并非 iid 生成时，Bartik 工具变量回归的矩条件可以写成：
这个条件描述了 Bartik 工具和第二阶段误差项之间的正交关系，其中 为地区 的劳动力在全国市场中所占份额。根据 Frisch-Waugh-Lovell 定理，我们将 Bartik 工具变量回归系数写为残差化的 和 的二元关系：
其中 和 分别表示将 和 对控制变量 做回归 (以 作为权重) 所得到的残差。
在产业 (冲击) 层面的正交等价条件可以表达为：
其中， 可以解读为某产业在总经济中的所占比例， 则是冲击-权重的误差项平均。这个等式的左边描述了 Bartik 工具 和误差项 的正交关系，以 为权重，右边则描述了冲击项 和 的正交关系，以 为权重。
由此可以看出，右边的条件成立时，Bartik 工具的外生性成立。上式中的等价关系为我们引出 Shift Share IV 估计等价理论提供了灵感。接着，使用 Bartik 工具变量对原始方程进行回归所得的系数 、以 为权重的冲击层面、和以 为工具对以下方程进行工具变量回归等价：
其中 表示冲击-权重的变量 的均值。
证明：
上述证明中最后边的等式即以 为权重的冲击层面，以 为工具的工具变量回归。
由此我们可以概括性地将 Bartik 工具变量的估计等价关系记忆为：地区层面的工具变量回归，等价于将变量转换为冲击 (或称产业) 层面后的工具变量回归。前者的估计单位，或称每个观测值为地区，而后者中，每个观测值为产业。当然需要注意，在这里两个工具变量回归使用的权重并不相同，前者使用的权重是地区的劳动力份额，而后者使用的权重为 ，可以理解为各冲击在不同地区的 (劳动力份额) 加权平均。
当然 Shift Share IV 的理论建立并不仅仅包括了此估计等价理论，更为重要的是我们可以通过将地区层面的 Bartik 工具变量估计视为等同于冲击层面的 Shift Share IV 估计，于是 “随机分配” 的冲击便可以使得 Shift Share IV ——以及等价的原 Bartik 工具变量，得到真实 的一致估计。
3. ssaggregate 命令
既然 Shift Share IV 的实施需要先将数据从产业层面转换为冲击层面， Shift Share IV 的几位作者便编写了 ssaggregate 命令来满足这一前提需求。这个变量实际上是先将地区层面的被解释变量和解释变量对控制变量残差化，再根据冲击比例在冲击层面上进行加权平均，即 。
*命令安装
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 的文章中，作者关注的变量如下：
被解释变量 ，是 时期地区 的制造业就业人数占总体工作年龄人口的比例的增加值； 解释变量 ，是每名美国工人水平的地区 的进口冲击大小。值得注意的是，这里的解释变量本身也是一个内积的形式，表达为 。其中 为劳动力雇佣数量， 和 分别表示地区与产业， 表示两期内 (1990-2000，2000-2007) 美国对中国的产品进口额对应增加值； 工具变量为 。其中 是地区 产业 的就业份额，使用两期分别十年前的数据来衡量。 是同期世界上另外八个高收入国家在产业 对中国的产品进口增加值 (每名美国工人水平) 之和，用于衡量因为中国的产品竞争力提高而对美国的各产业造成的冲击。
此处使用的数据来源于 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 实际上是 ，即各地区劳动力比例。
. 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：排他性约束及经典文献解读
连享会：因果推断专题（报名持续进行中）
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下：连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。