目录
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 提供了如下的选项：
controls(strings)：控制变量组，用于 partialling out 获得残差，可设定多组。如 controls('var1' 'var1 var2' 'var1 var2 var3')；absorb(strings)：固定效应变量组，用于吸收固定效应，可设定多组。如 absorb('f1' 'fe1 fe2')；
如果产业权重是 'long' 格式，那么必须设定如下四个选项：
sfilename：'long' 数据格式的产业权重数据文件；
同时，我们可以使用 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-
