多期DID文献解读：含铅汽油与死亡率和社会成本-L113

发布时间：2021-05-05 阅读 7975

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

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

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

✌ 课程详情： 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 | 工具软件 | 其它

作者： 邹恬华 (对外经济贸易大学)
邮箱： 201802060@uibe.edu.cn

Source: Hollingsworth, Alex., and Ivan Rudik, The effect of leaded gasoline on elderly mortality: Evidence from regulatory exemptions. 2021, Forthcoming at American Economic Journal: Economic Policy. -Link-, code, Non-technical summary

注：官方 code 链接在文章于 2021 年 8 月正式发表于 AEJ 前不可用。邮件询问作者后，作者提供了部分数据集（基于 R）与代码（基于 Stata）。

小插曲： 我是在 Github 上闲逛时看到这篇文章的，觉得很有意思，研究思路和方法都有很多值得借鉴的地方。于是就在连享会的 BAT 项目组发布了这个选题。邹恬华同学认领了这个任务。然而，兴冲冲地读完论文后，才发现作者提供的数据和 codes 下载链接无法打开，导致论文的复现工作搁浅。于是，便有了下面的小插曲。其实，在给作者发邮件之前，我们访问了作者的主页，也到连享会 - 论文重现网址页面提供的各个平台搜索了这篇文章的相关信息，一无所获。幸运的是，作者很 Nice，很快回复了我们的邮件，也附上了文中使用的备份数据和代码。所以，任何时候，只要不轻言放弃，总会有好结果的。(by 连玉君, 2021/5/2 23:30)

本文第一作者 Hollingsworth, Alex，https://hollina.github.io/

文章计量亮点：

政策动态效应分析：政策动态影响图（本文 4.2.2 节）
政策异质性效应分析：时间与空间滞后效应图（本文 4.1.3 节）

文章复现：

模型构建： 使用美国《清洁空气法》中的赛车免税条款，构建准自然实验。2007 年，美国全国汽车赛车协会 (NASCAR) 和美国汽车赛车俱乐部 (ARCA) ，从含铅汽油转换为无铅汽油。
结论： 在比赛中使用含铅汽油会增加环境铅，升高的血铅率和老年人死亡率；死亡率估计表明，每克汽油中添加的铅造成的损失超过 1,100 美元。

1. 文章简评与亮点速览

Hollingsworth et al. (2021) 利用 2007 年含铅汽油转换为无铅汽油的准自然实验，通过评估该事件对环境铅含量、血铅含量、老年人死亡率的政策动态影响，从而评估含铅汽油的环境成本。

从计量方法来看，文章使用多期 DID 评估了政策的动态影响，并进行以下分析：

政策动态效应分析：政策动态影响图（本文 4.2.2 节）
政策异质性效应分析：时间与空间滞后效应图（本文 4.1.3 节）

上述图形的经济含义与绘图代码，详见以下文章复现部分。

2. 背景介绍

铅作为一种有毒重金属，被认为与多种健康隐患有关。世界各国从 20 世纪 70 年代开始逐步采取行动减少铅排放。尽管几乎所有国家都禁止公路含铅汽油的使用，但含铅汽油仍在航空业与赛车业拥有广泛应用。

赛车业的铅排放问题较为严重。下图比较了赛车业、航空业的铅排放量与有毒物质释放清单 (Toxic Release Inventory , TRI) 中各设施的年空气铅排放量。结果显示，一场 3 小时的赛车比赛所排放的铅含量超过了 TRI 中 70%的含铅设施的年排放量，相当于一个机场全年的平均含铅量。若将所有赛车比赛铅排放加总，其每年预估铅排放将达到 200 万克，在有毒物质释放清单 (Toxic Release Inventory , TRI) 所有 5000 多个设施中排名第 30。

3. 数据获取与处理

文章主要进行了三个多期 DID 检验（即政策分别对环境铅含量、血铅含量、老年人死亡率的动态影响）。以下列出所有数据来源与处理，处理后的数据也可以直接在作者主页获取。官方 code 链接在文章今年八月于 AEJ 正式发表前不可用，作者提供了基于 $R$ 的部分数据集。

3.1 比赛数据

指标： 驾驶员数量，完成比赛所需的圈数，每圈长度，每位驾驶员在比赛中实际行驶的英里数，赛道位置和比赛日期。

样本：

处理组：NASCAR 三个 national series 和 ARCA Menards Series ，这四个系列于 2007 年改用无铅汽油。
对照组：NASCAR 和 ARCA 其它系列。
其它对照组 (稳健性检验使用) ：NASCAR K&N Pro Series，未知含铅燃料状态。NTT IndyCar Series，从不使用含铅汽油，且与 NASCAR、ARCA 使用相同赛道。

来源： 主要比赛数据来自 Racing Reference 网站，NTT IndyCar Series 比赛数据来自 ChampCarStats ，赛车场地理位置信息来自 Google Maps 手动匹配。

3.2 污染数据

污染数据来自环保署 (EPA) 的 AirData 数据集，包含 1980 年到 2018 年，每个监测日每个监测点的 CO，NO2，PM10，PM2.5，Pb，VOC 和 O3 数据。1980 年之前的每日铅水平数据由作者向环保署 (EPA) 索要得到。

3.3 气象数据

气象数据来自国家海洋大气管理局 (National Oceanic and Atmospheric Administration , NOAA) 的 NCEP-NCAR Reanalysis 1 模型，该模型提供了 2.5° 经纬度网格下全球每日气象数据。本文使用了 1948 年-2018 年的地表层面 (数据库中 0.995 sigma 水平) 数据，包括平均空气温度 (C°) ，压力 (kPA) ，相对湿度 (%) 和风的数据。我们使用数据库中风的矢量数据计算了风速与风向。

由于降水量无地面层面上的数据，本文使用 2.5° 经纬度网格下大气柱每日降水量作为替代。

最后根据 2.5° 经纬度网格位置，将获得与每个污染物监测点匹配的气象数据。

3.4 血铅含量数据

血铅含量数据来自疾病控制和预防中心 (Centers for Disease Control and Prevention , CDC) ，包含 2005 年-2015 年中县级-年份水平的数据。数据报告了在给定的县-年度接受测试的儿童被观察到血铅水平升高的百分比。

关注的指标为血铅升高确诊病例数，即静脉抽血检查中血铅浓度超过 10µg/dl，或两周内两次毛细血管抽血检查中血铅浓度超过 10µg/dl。

3.5 死亡率数据

死亡率数据来自疾病控制和预防中心(CDC)国家生命统计系统。死亡率数据详细至个人层面。我们使用数据库中的死因、死亡年龄、死亡县、死亡月份、死亡年份与居住的县的数据。

人口数据来自国家癌症研究所的监测、流行病学和最终结果 (Surveillance, Epidemiology, and End Result , SEER) 项目。SEER 数据跨度为 1999 年-2016 年。本文使用的指标是每个人群中每 10 万人的年龄标准化死亡率，使用 2000 年每个人群的美国人口进行标准化。

本文将重点放在老年人口 (65 岁以上) ，但同时也在附录部分报告了婴儿(1 岁以下)的结果。

3.6 其它数据

社会经济变量：收入中位数的县级数据来自美国人口普查，失业率数据来自劳工统计局，非白人百分比的数据来自 SEER 项目。

其它行业的铅排放数据：来自有毒物质释放清单 (Toxic Release Inventory , TRI) 。由于 TRI 数据样本须达到一定规模门槛，我们还使用县商业模式数据集中制造业的县级薪资数据来补充铅排放数据，解决了对较小设施的审查问题。同时我们还解决了 TRI 数据库中的已知测量误差问题。

4. 实证检验

文章主要进行了三个多期 DID 检验（即政策分别对环境铅含量、血铅含量、老年人死亡率的动态影响），分别汇报于 4.1-4.3 节。

4.1 环境铅浓度

4.1.1 描述性统计

下图展示了不同时间和空间的环境铅浓度描述性统计结果。左图 (A)、右图 (B) 分别标明了 1996-2006 年、2007-2017 年比赛开始前 4 周和比赛结束后 4 周的平均环境铅浓度。不同的是，在 2007-2017 年间，比赛使用无铅汽油。X 轴上数值为对于比赛举办的相对周数。

左图 (A) 表明在含铅的比赛年，赛前的铅浓度在赛道 50 英里内（图中蓝线）呈轻微下降趋势，而在更远的地方持平或略有上升（图中红线）。在比赛结束后的一周内，赛道 50 英里内监视器的浓度上升了 20%以上，之后缓慢减弱，直到四周后达到赛前水平。在 50 英里以外的地方，浓度基本上没有变化。

右图 (B) 表明在 NASCAR 和 ARCA 开始不适用含铅燃料（2007 年）后，比赛后周围环境铅浓度的跃升就消失了。

4.1.2 计量模型

建立计量模型如下：

\begin{aligned} asinh (P b)_{m d w y} = & β_{l} {leaded race miles}_{m d w y} \\ + β_{u} {unleaded race miles}_{m d w y} \\ + γ X_{m d w y} + Θ_{m y} + Ω_{w y} + ε_{m d w y} . \end{aligned} (1)

其中，

$asinh (P b)_{m d w y}$ ：第 $y$ 年第 $w$ 周第 $d$ 天在 $m$ 点监测到的铅浓度的反双曲正弦变换。使用反双曲正弦变化的原因是由于铅浓度分布右偏的，并且有大量的零值观测值。
${leaded race miles}_{m d w y}$ ：含铅里程，即 $y$ 年 $w$ 周 $d$ 天 $m$ 监测点，前一周内、50 英里内的使用含铅汽油比赛的赛车总里程。
${unleaded race miles}_{m d w y}$ ：不含铅里程，即 $y$ 年 $w$ 周 $d$ 天 $m$ 监测点，前一周内、50 英里内的不使用含铅汽油比赛的赛车总里程。
$X_{m d w y}$ ：控制变量，包括平均温度、可降水量、相对湿度、表面压力和风速。
$Θ_{m y}$ ：监测点-年度固定效应，用来控制监测点铅浓度的时间趋势。
$Ω_{w y}$ ：星期-年份固定效应，用来控制不同监测点间的周内冲击。
$ε_{m d w y}$ ：干扰项。

感兴趣的系数包括：

$β_{l}$ ：过去的一周里驾驶的含铅里程数对 50 英里内环境铅读数的影响

$β_{u}$ ：过去的一周里驾驶的不含铅里程数对 50 英里内环境铅读数的影响

采用多重固定效应命令(reghdfe)，回归结果如下图所示：

(1) 列仅控制监测点和年份固定效应。结果显示每 10 万公里含铅里程对 $asinh (P b)_{m d w y}$ 的影响为 0.16，每 10 万公里无铅里程对 $asinh (P b)_{m d w y}$ 的影响为 0.07。

(2) - (3) 列还控制了天气与星期固定效益，两个系数都没有发生显著变化。

(4) 列还控制了州-年份固定效应，系数分别变为 0.13 和 0.02。

进一步，(5) - (7) 列将空间固定效应由州细化为县、监测点 (列 6) ，并添加年份-星期固定效应。结果显示无铅比赛里程和环境铅水平之间的估计关系逐渐削弱，而让含铅比赛里程对环境铅浓度的估计影响相对不变。

4.1.3 事件研究法与空间滞后模型

文章进一步估计了含铅里程与不含铅里程对环境铅浓度的时间动态影响与距离递减效应。

针对时间动态影响建立了事件研究法模型，在基准回归 (1) 式中加入 14 个新变量：比赛前 2~4 周、比赛后 1~4 周中每周分别对应的 50 英里内含铅里程数与不含铅里程数。（比赛前 1 周对应的两个变量已在基准回归方程中）

针对距离递减效应建立了空间滞后模型，在基准回归 (1) 式中加入 4 个新变量：50-75 公里、75-100 公里分别对应的比赛前一周含铅里程数与不含铅里程数。

以空间滞后模型（右图）为例，其估计代码与解释如下：

// 为了展示主要步骤，此处省略数据处理代码
// 设置全局变量
global weather_controls air_temp pr_wtr pres rhum wind_speed
global cluster_variable fips
global fixed_effects i.site_id##i.year i.week##i.year

// 进行回归
// inv_sine指铅浓度的反双曲正弦变换，m_l开头变量指含铅历程，m_u开头变量指不含铅里程
// 变量结尾的数字指在哪个范围内
// eg. m_l_past_7_top_arca_0_50指0-50公里内比赛前一周含铅里程数
reghdfe inv_sine ///
	m_l_past_7_top_arca_0_50  m_u_past_7_top_arca_0_50 ///
	m_l_past_7_top_arca_50_75  m_u_past_7_top_arca_50_75 ///
	m_l_past_7_top_arca_75_100  m_u_past_7_top_arca_75_100 ///
	$weather_controls ///
	if year > 1995 & regional_miles == 0 ///
	, absorb($fixed_effects) vce(cluster $cluster_variable)

// 储存回归结果
// 若存在order,b,high,low变量，则删除
capture drop order
capture drop b
capture drop high
capture drop low
// 设定order,b,high,low变量，默认值为空
// b指回归系数，high和low分别指置信区间上下界，order为横坐标顺序
gen order = .
gen b =.
gen high =.
gen low =.
// 对order,b,high,low变量赋值
replace order = 1 in 1
replace b    = _b[m_l_past_7_top_arca_0_50] in 1
replace high = _b[m_l_past_7_top_arca_0_50] + 1.96*_se[m_l_past_7_top_arca_0_50] in 1
replace low  = _b[m_l_past_7_top_arca_0_50] - 1.96*_se[m_l_past_7_top_arca_0_50] in 1

replace order = 2 in 2
replace b    = _b[m_l_past_7_top_arca_50_75] in 2
replace high = _b[m_l_past_7_top_arca_50_75] + 1.96*_se[m_l_past_7_top_arca_50_75] in 2
replace low  = _b[m_l_past_7_top_arca_50_75] - 1.96*_se[m_l_past_7_top_arca_50_75] in 2

replace order = 3 in 3
replace b    = _b[m_l_past_7_top_arca_75_100] in 3
replace high = _b[m_l_past_7_top_arca_75_100] + 1.96*_se[m_l_past_7_top_arca_75_100] in 3
replace low  = _b[m_l_past_7_top_arca_75_100] - 1.96*_se[m_l_past_7_top_arca_75_100] in 3

replace order = 5 in 4
replace b    = _b[m_u_past_7_top_arca_0_50] in 4
replace high = _b[m_u_past_7_top_arca_0_50] + 1.96*_se[m_u_past_7_top_arca_0_50] in 4
replace low  = _b[m_u_past_7_top_arca_0_50] - 1.96*_se[m_u_past_7_top_arca_0_50] in 4

replace order = 6 in 5
replace b    = _b[m_u_past_7_top_arca_50_75] in 5
replace high = _b[m_u_past_7_top_arca_50_75] + 1.96*_se[m_u_past_7_top_arca_50_75] in 5
replace low  = _b[m_u_past_7_top_arca_50_75] - 1.96*_se[m_u_past_7_top_arca_50_75] in 5

replace order = 7 in 6
replace b    = _b[m_u_past_7_top_arca_75_100] in 6
replace high = _b[m_u_past_7_top_arca_75_100] + 1.96*_se[m_u_past_7_top_arca_75_100] in 6
replace low  = _b[m_u_past_7_top_arca_75_100] - 1.96*_se[m_u_past_7_top_arca_75_100] in 6

// 绘图
twoway 	///
 scatteri 0 0 0 7.5, recast(line) lcol(gs8) lp(longdash) lwidth(0.5) /// y=0基准线
	|| rcap low high order if order<8, col(white) msize(5) lwidth(.5) ///  设置置信区间
	|| rcap low high order if order<8, col(black) msize(4) lwidth(.25) /// 设置置信区间
	|| scatter b order if order<4, col(white) msize(5) msymbol(s) /// 设置含铅里程组
	|| scatter b order if order<4, col(black) msize(4) msymbol(s) /// 设置含铅里程组
	|| scatter b order if order<8 & order>3, col(white) msize(5) msymbol(T) /// 设置无铅里程组
	|| scatter b order if order<8 & order>3, col(black) msize(4) msymbol(T) /// 设置无铅里程组
	text(.25 1.6 "Leaded", place(e) si(5)) ///
	text(.25 5.3 "Unleaded", place(e) si(5)) ///
		xlab(0 " " 1 "0-50" 2 "50-75  " 3 "    75-100" 4 " " 5 "0-50" 6 " 50-75  " 7 "    75-100" ///横坐标设置
				, nogrid labsize(5) notick) ///
		ylab(-.1(.1).25 , nogrid labs(5)) ///纵坐标设置
		legend(off) ///
		xtitle("Distance from racetrack to monitor (mi)", size(5)) ///横坐标标题
	subtitle("Effect of 100,000 race miles on asinh(Pb)", size(5) pos(11)) ///纵坐标标题
	graphregion(margin(r+10))

// 导出与保存
graph export "$results_path/figure_3b_distance_gradient.png", replace
graph export "$results_path/figure_3b_distance_gradient.pdf", replace

估计结果如下图所示：

左图 (A) 为事件研究法结果。含铅里程（左半边），在比赛前几周没有引起显著的铅浓度变化，而在比赛后第一周， $asinh (P b)_{m d w y}$ 增加了 0.15 ，平均浓度增加了 16% 。在之后的几周内，这种效应逐渐趋于零。所有无铅里程的估计（右半边）都以零为中心，没有明显的趋势，这表明无铅比赛不影响环境铅水平。

右图 (B) 为空间滞后模型结果。含铅里程（左半边），对 $asinh (P b)_{m d w y}$ 的影响在距离上有一个直观的递减梯度。我们估计在 0-50 英里的距离范围内影响为 14% ，在 50-75 英里范围内下降为 3% ，在 75-100 英里范围内大约为零。所有无铅里程的估计（右半边）都以零为中心，没有明显的趋势。

4.2 血铅浓度

4.2.1 描述性统计

下图绘制了从 2005 年到 2009 年县级平均高血铅患病率（加权的权重是参与测试儿童数的平方根）。在无铅改革前（2007 年前，图中虚线左侧），比赛县（蓝线）、邻县有比赛的县（绿线）、控制组（红线）所有这三种类型的县的高血铅患病率都有轻微的下降趋势，而且几乎是平行的，尽管邻县有比赛的县和比赛县的高血铅患病率分别高于控制组 1 和 2 个百分点。在 2007 年，当 NASCAR 和 ARCA 开始使用无铅汽油时，控制县的高血铅患病率继续沿相同的趋势变化，但比赛县、邻县有比赛的县的高血铅患病率下降至几乎与控制组相同的水平。在随后的几年里，不同组别的高血铅患病率有相似的水平并具有平行趋势。

4.2.2 计量模型

由于血铅浓度数据水平为县级，故建立事件研究法模型如下：

\begin{aligned} asinh (E B L L)_{s c y} = & \sum_{t \in Y} β_{t}^{r} 1 (race in county)_{s c t} \\ + \sum_{t \in Y} β_{t}^{b} 1 (race in border county)_{s c t} \\ + γ X_{s c y} + Θ_{c} + Ω_{s y} + ε_{s c y} \end{aligned} (2)

其中，

$asinh (E B L L)_{s c y}$ ： $s$ 州 $c$ 县 $y$ 年高血铅患病率的反双曲正弦变换。我们再次使用反双曲正弦变换，因为高血铅率的分布是右偏的，并且包含大量的零值观测值。
$\sum_{t \in Y} 1 (race in county)_{s c t}$ ：比赛县虚拟变量，如果 $c$ 县在 $t$ 年有比赛，则该虚拟变量等于 1。 $Y$ 代表 2005-2015 年的年份集合，但不包括 2007 年。
$β_{t}^{r}$ ：反映了 $t$ 年本县赛车比赛对 $y$ 年高血铅患病率的影响。第 $t$ 年的控制县是本县与邻县都没有比赛的县。这允许一些县在一年中作为对照组，在另一年中作为一个治疗组。然而，由于血铅的半衰期很短，而且在一年的前 40 天或后 40 天一般没有 NASCAR 或 ARCA 比赛，所以含铅比赛对次年血铅升高发生率的影响不大。
$\sum_{t \in Y} 1 (race in border county)_{s c t}$ ：邻县有比赛的县虚拟变量，如果 $c$ 县的邻县在 $t$ 年有比赛，且本县没有比赛，则该虚拟变量等于 1。这组变量是用来研究高血铅患病率是否和前文研究的环境铅一样具有距离递减效应。
$β_{t}^{b}$ ：反映了 $t$ 年邻县赛车比赛（但本县无赛车比赛）对 $y$ 年高血铅患病率的影响。
$X_{s c y}$ ：控制变量，包括失业率、中位数收入、非白人百分比、制造业工资和 TRI 铅排放量。前三项是控制社会经济状态和潜在的跨地区差异，而后两项是控制铅排放的其它可能来源。
$Θ_{c}$ ：县固定效应
$Ω_{s y}$ ：州-年份固定效应
$ε_{s c y}$ ：干扰项

回归采用异方差稳健标准误，并汇报县级聚类标准误。回归结果绘图代码：


***********************************
**两期DID，用来估计前两年平均效应**
***********************************

// DID回归，被解释变量ihs_bll，交叉项i.pre_2007##i.race（是否先于2007年、是否是比赛县）
reghdfe ihs_bll        ///
	i.pre_2007##i.race ///
	$control_list      ///
	[aw=weight]        ///
	, absorb(county_id i.StateFIPS##i.year) ///
      vce(cluster county_id)

// 回归系数与p统计量
local nascar_estimate = round(_b[1.pre_2007#2.race],.01)
local border_estimate = round(_b[1.pre_2007#1.race],.01)
local nascar_p = 2*ttail(e(df_r),_b[1.pre_2007#2.race]/_se[1.pre_2007#2.race])
local border_p = 2*ttail(e(df_r),_b[1.pre_2007#1.race]/_se[1.pre_2007#1.race])

// 设置显著性水平星星个数
if `nascar_p' <= .01 {
	    local nascar_stars = "***"
	}
	else if `nascar_p' <= .05 {
	    local nascar_stars = "**"
	}
	else if `nascar_p' <= .1 {
	    local nascar_stars = "*"
	}
	else {
	    local nascar_stars = ""
	}

if `border_p' <= .01 {
	    local border_stars = "***"
	}
	else if `border_p' <= .05 {
	    local border_stars = "**"
	}
	else if `border_p' <= .1 {
	    local border_stars = "*"
	}
	else {
	    local border_stars = ""
	}
	di "`border_stars'"

***********
**多期DID**
***********

// 多期DID回归
// ib2007.year：因子变量，以2007年为基年
reghdfe ihs_bll ///
	ib2007.year##i.race ///
	$control_list ///
	[aw=weight]   ///
	, absorb(county_id i.StateFIPS##i.year) vce(cluster county_id)

// 若存在order,b,high,low变量，则删除
capture drop order
capture drop border*
capture drop b high low

// 设定order,b,high,low变量，默认值为空
// b指回归系数，high和low分别指置信区间上下界，order为年份绘图顺序
// 无border前缀为比赛县数据，有border前缀为邻县比赛的县数据
gen order = .
gen border_b =.
gen border_high =.
gen border_low =.
gen b =.
gen high =.
gen low =.
local i = 1

// 对2005-2006年order,b,high,low变量赋值
forvalues yr = 2005(1)2006{
	replace order = `i' in `i'

	replace border_b    = _b[`yr'.year#1.race] in `i'
	replace border_high = _b[`yr'.year#1.race] + 1.96*_se[`yr'.year#1.race] in `i'
	replace border_low 	= _b[`yr'.year#1.race] - 1.96*_se[`yr'.year#1.race] in `i'

	replace b    = _b[`yr'.year#2.race] in `i'
	replace high = _b[`yr'.year#2.race] + 1.96*_se[`yr'.year#2.race] in `i'
	replace low  = _b[`yr'.year#2.race] - 1.96*_se[`yr'.year#2.race] in `i'

	local i = `i' + 1
}

// 对2007年order,b,high,low变量赋值（此时`i'=3）
replace order = `i' in `i'

replace border_b    = 0  in `i'
replace border_high = 0  in `i'
replace border_low  = 0  in `i'

replace b    = 0  in `i'
replace high = 0  in `i'
replace low  = 0  in `i'
local i = `i' + 1

// 对2008-2015年order,b,high,low变量赋值（此时`i'=4）
forvalues yr = 2008(1)2015{
	replace order = `i' in `i'

	replace border_b 	= _b[`yr'.year#1.race] in `i'
	replace border_high = _b[`yr'.year#1.race] + 1.96*_se[`yr'.year#1.race] in `i'
	replace border_low 	= _b[`yr'.year#1.race] - 1.96*_se[`yr'.year#1.race] in `i'

	replace b    = _b[`yr'.year#2.race] in `i'
	replace high = _b[`yr'.year#2.race] + 1.96*_se[`yr'.year#2.race] in `i'
	replace low  = _b[`yr'.year#2.race] - 1.96*_se[`yr'.year#2.race] in `i'

	local i = `i' + 1
}


**********************
**绘制上图（比赛县）**
**********************

// 设置配色方案
colorpalette hsv, dark n(3) nograph
return list
// 设置文字位置
local year_pointer_label = _b[2008.year#2.race] - 1.96*_se[2008.year#2.race]
local county_label = _b[2008.year#2.race] - 1.96*_se[2008.year#2.race] -.07
local year_pointer_end = `year_pointer_label' + .03
local pointer_lower = `nascar_estimate' + .02
local pointer_upper = `nascar_estimate' + .07
local text_location = `pointer_upper' + 0.025

// 绘图
twoway rarea low high order if order<18, fcol(gs14) lcol(white) msize(3) /// 估计系数
	|| connected b order if order<18, lw(1.1) col(white) msize(7) msymbol(s) lp(solid) /// 强调
		text(`county_label' 1 "Race counties", place(e) si(7)) ///
	|| connected b order if order<18, lw(0.6) col("71 71 179") msize(5) msymbol(s) lp(solid) /// 连线
	|| scatteri 0 1 0 11, recast(line) lcol(gs8) lp(longdash) lwidth(0.5) /// y=0基准线
	|| pcarrowi `year_pointer_end' 3 -.03 3, mlcolor(black) mfcolor(black) lcolor(black) barbsize(1.5) /// 指向第一年使用无铅汽油的箭头
		text(`year_pointer_label' 2.8 "First majority unleaded year", place(e) si(large)) /// 添加文字"First majority unleaded year"
	|| scatteri `nascar_estimate' 1 `nascar_estimate' 2.55, recast(line) lcol(black) lp(---)  /// 2005-2006年的平均效应，数值由两期DID给出
	|| pcarrowi `pointer_upper' 2.2 `pointer_lower' 2.2, mlcolor(black) mfcolor(black) lcolor(black) barbsize(1.5) /// 指向平均效应的箭头
		text(`text_location' 1.8 "Average effect relative to 2007-2015: `nascar_estimate'{sup:`nascar_stars'}", place(e) si(large)) /// 添加文字"Average effect relative to 2007-2015:.16***"
		xlab(1 "2005" 2 "2006" 3 "2007" 4 "2008" 5 "2009" 6 "2010" 7 "2011" ///
			 8 "2012" 9 "2013" 10 "2014"  11 "2015" /// 横坐标设置
				, nogrid labsize(7) angle(0)) ///
		ylab(, nogrid labs(7)) ///
		legend(off) ///
		xtitle("", size(7))
// 文件保存
graph save "$results_path/bll.gph", replace

****************************
**绘制下图（邻县有比赛县）**
****************************

// 设置文字位置
local year_pointer_label = _b[2015.year#1.race] - 1.96*_se[2015.year#1.race]
local year_pointer_end = `year_pointer_label' + .03
local pointer_lower = `border_estimate' + .02
local pointer_upper = `border_estimate' + .07
local text_location = `pointer_upper' + 0.025

// 绘图
twoway rarea border_low border_high order if order<18, fcol(gs14) lcol(white) msize(3) /// 估计系数
	|| connected border_b order if order<18, lw(1.1) col(white) msize(7) msymbol(d) lp(solid) /// 强调
		text(`county_label' 1 "Border counties", place(e) si(7)) ///
	|| connected border_b order if order<18, lw(0.6) col("71 179 71") msize(5) msymbol(d) lp(solid) /// 连线
	|| scatteri 0 1 0 11, recast(line) lcol(gs8) lp(longdash) lwidth(0.5) /// y=0基准线
	|| scatteri `border_estimate' 1 `border_estimate' 2.55, recast(line) lcol(black) lp(---) /// 2005-2006年的平均效应，数值由两期DID给出
	|| pcarrowi `pointer_upper' 2.2 `pointer_lower' 2.2, mlcolor(black) mfcolor(black) lcolor(black) barbsize(1.5) /// 指向平均效应的箭头
	text(`text_location' 1.8 "Average effect relative to 2007-2015: `border_estimate'{sup:`border_stars'}", place(e) si(large)) /// 添加文字"Average effect relative to 2007-2015:.05"
		xlab(1 "2005" 2 "2006" 3 "2007" 4 "2008" 5 "2009" 6 "2010" 7 "2011" ///
			 8 "2012" 9 "2013" 10 "2014"  11 "2015" /// 横坐标设置
				, nogrid labsize(7) angle(0)) ///
		ylab(, nogrid labs(7)) ///
		legend(off) ///
		xtitle("", size(7))

graph save "$results_path/bll_border.gph", replace

**************
**上下图合并**
**************

graph combine ///
	"$results_path/bll" /// 上图文件名（比赛县）
	"$results_path/bll_border" /// 下图文件名（邻县比赛的县）
	, altshrink  ycommon col(1) ysize(4.5) ///
	subtitle("     Effect of a race on asinh(EBLL)", size(4) pos(11))
graph export "$results_path/figure_5_bll_event_study.png", replace
graph export "$results_path/figure_5_bll_event_study.pdf", replace

回归结果如下图所示：

上半部分汇报了比赛县结果。2005 年和 2006 年比赛县的高血铅患病率较 2007 年分别高 18% 和 13% 。2005 年与 2006 年的差异原因是 2006 年比赛里程中有 14% 为无铅里程。2007 年后的所有年份，点估计量都接近于零。图中虚线是 DID 的点估计量，DID 回归中的比赛县、邻县有比赛的县的虚拟变量替换为另外两个指标，一个是该县在 2005 年或 2006 年是否有赛车比赛的虚拟变量，另一个是邻县在 2005 年或 2006 年是否有赛车比赛且本县无赛车比赛的虚拟变量。回归结果显示，高血铅患病率在含铅比赛年比无铅比赛年平均高 17% 。

下半部分汇报了邻县有比赛的县的结果。2005 年和 2006 年邻县有比赛的县的高血铅患病率较 2007 年分别高 6% 和 5% 。2007 年后的所有年份，点估计量都接近于零。DID 结果显示高血铅患病率在含铅比赛年比无铅比赛年平均高 6%。此结果相比比赛县 17% 的结果更小，这与存在距离递减效应的结论一致。

事件分析法的结论阐明了比赛县具有更高高血铅患病率的机制。高血铅症可能由赛车比赛中增量铅与之前的存量铅引起。比赛改用无铅汽油减少了铅的空气传播，同时避免了增加铅的存量。如果高血铅发病率与当地铅存量具有函数关系，我们应当观察到比赛县相比控制组拥有更低的高血铅发病率下降速度，或高血铅发病率根本不下降。我们在描述性统计与回归结果中实际观察到的 2007 年及之后年份高血铅发病率的迅速下降，印证了这一假设。

4.3 老年人死亡率

4.3.1 描述性统计

下图绘制了经过年龄标准化的平均全因素死亡率，加权权重为各县按老年人口数的平方根。

在抽样期间，比赛县（蓝线）、邻县有比赛的县（绿线）、控制组（红线）三种类型的县都呈现出急剧下降的趋势。在更换无铅燃料前，比赛县和控制组的全因素死亡率几乎相同，而邻县有比赛的县的死亡率稍低，但呈平行趋势。 2007 年，比赛县的死亡率降到了控制组以下，而邻县有比赛的县的死亡率也相对于控制组下降了。 2007 年后的趋势也大致平行。

4.3.2 计量模型

使用与研究高血铅发病率相同的计量模型：

\begin{aligned} {mortality rate}_{s c y} = & \sum_{t \in Y} β_{t} 1 (race in county)_{s c t} \\ + \sum_{t \in Y} β_{t}^{b} 1 (race in border county)_{s c t} \\ + γ X_{s c y} + Θ_{c} + Ω_{s y} + ε_{s c y} \end{aligned} (3)

其中，

${mortality rate}_{s c y}$ ： $s$ 州 $c$ 县 $y$ 年的某特定原因「老年人死亡率」，数据经过年龄标准化。老年人定义为 65 岁以上人群，死亡率以 10 万人中死亡人数计。由于数据大致呈现正态分布且零值较少，不进行反双曲正弦变换。
其余变量与研究高血铅发病率相同的计量模型设定完全相同。

文章探讨了六种类型的老年人死亡率 ${mortality rate}_{s c y}$ ：

全因素死亡率。以避免潜在的选择偏误。
心血管疾病死亡率。
缺血性心脏病死亡率 (IHD) 。
呼吸道疾病死亡率。
绝望死亡率 (Case and Deaton , 2015)。如抑郁、惊恐发作、自杀等。
糖尿病死亡率。由于糖尿病是慢性病，不太可能受环境铅含量的突然变化的影响，故作为本文的安慰剂。

回归结果如下图所示：

Panel A 汇报了全因素死亡率结果。2007 年前，比赛县（蓝线）、邻县有比赛的县（绿线）与控制组无明显差距，2007 年比赛县、邻县有比赛的县全因素死亡率相比控制组下降，2007 年后持续下降。与前文类似，虚线指 DID 结果。

Panel B-E 汇报了其它四种死亡率结果。除呼吸道疾病死亡率 (Panel D) 2007 年前有一些差距外，结论与全因素死亡率结果相似。

使用无铅燃料对降低死亡率的效果如下表所示，邻县有比赛的县的死亡率降低幅度小于比赛县，这与存在距离递减效应的结论一致。

降低死亡率(人/10 万人)	比赛县	邻县有比赛的县
全因素死亡率	91	38
心血管疾病死亡率	37	12
缺血性心脏病死亡率	53	20
呼吸道疾病死亡率	20	5
绝望死亡率	4	2

5.含铅汽油的社会成本

文章分两部分考虑含铅汽油的社会成本。

1.含铅汽油导致的血铅含量上升，及其引致的 IQ 下降。

2.含铅汽油对老年人全因素死亡率的影响，用来估算含铅汽油社会成本的下限。

文章采用结合现有文献与自己计算数据逐步计算的方式估算含铅汽油的社会成本，我们列出作者计算的流程图如下：

5.1 血铅、IQ 与未来收入

5.2 老年人死亡率

6.结论

6.1 研究结果与意义

对含铅汽油对老年人死亡率的影响进行了因果分析
额外的铅暴露会导致更高的死亡风险
通过分析含铅汽油对老年人死亡率的影响，得出每克铅的社会成本下限为 1100 美元
含铅汽油的社会成本远大于市场价格中反映的私人成本。即使航空汽油的外部影响仅为我们赛车燃料的 1％（每克 10 美元），其社会有效价格也将是当前市场价格的三倍。这表明含铅航空汽油会对全球造成重大损害。
尽管目前的环境政策多关注儿童，但老年人实际承受了高昂的环境成本

6.2 注意事项与可能的改进

社会成本对于环境浓度的响应是非线性的
每克铅社会成本实际上取决于边际效应，这意味着了解成本函数及其边际效应至关重要
含铅汽油社会成本存在地点与时间上的异质性
缺乏含铅汽油销售数量与销售位置的信息
缺乏无铅汽油铅含量数据，无铅汽油并非完全无铅

7. 参考文献

Hollingsworth, Alex., and Ivan Rudik, The effect of leaded gasoline on elderly mortality: Evidence from regulatory exemptions. Forthcoming at American Economic Journal: Economic Policy. -Link-, code, Non-technical summary

8. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 论文重现重现经典文献 bad 多期DID
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]