DID从经典到前沿方法的保姆级教程, 释放最完整数据和代码！

凡是搞计量经济的，都关注这个号了

邮箱：[e-mail]

所有计量经济圈方法论丛的code程序, 宏微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.

关于多期DID或交叠DID: 1.DID相关前沿问题“政策交错执行+堆叠DID+事件研究”, 附完整slides，2.交错(渐进)DID中, 用TWFE估计处理效应的问题, 及Bacon分解识别估计偏误，3.典范! 这篇AER在一图表里用了所有DID最新进展方法, 审稿人直接服了！4.最新Sun和Abraham(2020)和TWFE估计多期或交错DID并绘图展示结果！详细解读code！5.多期DID或渐进DID或交叠DID, 最新Stata执行命令整理如下供大家学习，6.多期DID前沿方法大讨论, e.g., 进入-退出型DID, 异质性和动态性处理效应DID, 基期选择问题等，7.交叠DID中平行趋势检验, 事件研究图绘制, 安慰剂检验的保姆级程序指南！8.欣慰! 营养午餐计划终于登上TOP5! 交叠DID+异质性稳健DID！9.用事件研究法开展政策评估的过程, 手把手教学文章！10.从双重差分法到事件研究法, 双重差分滥用与需要注意的问题，11.系统梳理DID最新进展: 从多期DID的潜在问题到当前主流解决方法和代码! 12.标准DID中的平行趋势检验,动态效应, 安慰剂检验, 预期效应教程

今天，我们根据一份《DID handbook》再从头到尾整理一下，从经典DID到Bacon分解到交叠DID最新进展的各种估计方法适用条件及其实现代码，并跟随这份Handbook的例子复现各种DID估计方法的结果并对其差异进行对比。下方所有数据和代码都放在网盘里，可以直接放在Stata中运行出来结果。这真的是保姆级教程，如果看完这个你还不懂DID，那就真的需要再回到课堂学学了😄。

正文

《DID保姆级手册》

*感谢计量社群群友@regmoneky 整理，增添了部分内容的注释，原文来自下方Source。

Source：DID_Handbook，https://github.com/IanHo2019/DID_Handbook；欢迎对更好的编码或错误纠正提出意见，作者联系方式：[e-mail]。

双重差分法（也可缩写为DID、DiD或DD，我更倾向于使用DID）是当今社会科学中最流行的统计技术之一，广泛应用于定量研究中。其受欢迎的主要原因在于使用者对理解和应用DID开展实证研究方面相对"容易上手"。然而，在阅读了一系列最近（2017年至今）发表的高质量计量经济学论文后，我意识到双重差分并没有我之前所想的那么简单。构建这个资源库的主要目标是分享我对双重差分的更深入理解，以及我在运行双重差分时使用的Stata代码。请注意，我在这里只对每种双重差分格式的一些细节进行了概述，更详细的内容请阅读那些论文。

经典、传统和教科书双重差分

双重差分，顾名思义，涉及在两个时期对两个群体进行比较。第一个差分是在群体之间进行的，而第二个差分总是在时间之间进行的。在处理时间之后成为接受处理的群体，被称为接受处理的群体。其他单位被称为对照组。双重差分通常侧重于识别和估计接受治疗群体上的平均处理效应（ATT），即测量治疗对那些从未接受治疗转变为接受治疗的人的平均效应。在实证研究中，实施双重差分规范的主导方法是运行双向固定效应（TWFE）回归：

双重差分（Difference in Difference，简称DID）假设平行趋势、无预期和处理效应同质性，TWFE回归中的估计结果等于参与政策处理的因果效应。然而，不幸的是，这种TWFE回归通常对处理效应的异质性不具有稳健性。这是当前DID文献中一个热门的研究课题。

需要注意的是，这里我们只考虑吸收型政策处理（absorbing treatment）：一旦一个个体接受政策处理，它在未来的任何时间段内都不能退出政策处理。一些研究人员尝试将DID推广到非吸收型政策处理（也称为“转换型政策处理（switching treatment）”）设计；这是可行的，但需要额外的假设。

在Stata中，可以使用命令xtreg、areg、reghdfe或xtdidregress轻松运行TWFE回归。需要注意的是，xtdidregress仅适用于Stata 17或更高版本，而reghdfe只能在安装了reghdfe和ftools包之后使用。我喜欢使用reghdfe，因为它具有灵活的选项并且计算速度较快。reghdfe的基本语法如下：

reghdfe Y D, absorb(id t) cluster(id)

absorb选项指定固定效应，cluster选项指定标准误差聚类的级别。

这种双重差分的格式仅涉及两个时间段和两个群体，这就是为什么我称之为经典、传统和教科书格式。经典意味着“标准”：它是最早、最简单的一种，展示了DID背后的关键理论；传统意味着“传统”：它已经建立了很长时间，但在现代应用中逐渐过时；教科书意味着：它已经在许多教科书中被介绍过，比如Jeffrey Wooldrige的《横截面和面板数据的计量经济分析》（2010）和Bruce Hansen的《计量经济学》（2022）。如果那些接受处理的个体在同一时期进行处理，TWFE DID可以推广到多时期多个体的情况。这种处理有时被称为分块政策处理（block treatment）。

在2017年之前，许多研究人员天真地认为TWFE DID也可以轻松推广到个体在不同时期接受政策处理的情况（即交错处理staggered treatment）。然而，事实并非如此！正如下文所描述的那样，这绝对不容易。

静态双重差分下的Bacon分解

首先，什么是静态双重差分（Static DID）？静态双重差分的估计的是，一个不随时间变化的单一处理效应。也就是说，通过运行静态双重差分模型，我们只得到一个β值，并且使用一个系数总结自政策实施以来的处理效应。经典的双重差分正是一种静态双重差分。

其次，什么是Bacon分解（Bacon Decomposition，交错(渐进)DID中, 用TWFE估计处理效应的问题, 及Bacon分解识别估计偏误）？这个概念来自Goodman-Bacon（2021）。该作者提出并证明了双重差分分解定理，指出TWFE双重差分估计量（适用于交错处理情况）等于所有可能的两组/两时期双重差分估计量的加权平均值。为了强调这一发现，作者将双重差分估计量写为。

双重差分分解定理非常重要，因为它告诉我们，在经典的双重差分中，如果包括在多个时间段接受政策处理的个体，就存在一种“不好”的比较情况，即在晚期处理之前和之后将晚期处理组与早期处理组进行比较。建议在运行静态双重差分模型时进行Bacon分解，可以看出我们的估计结果更多来自何处的贡献。例如，一个负的双重差分估计结果可能只是因为在一个权重较大的不好的比较中出现了负的结果。

在Stata中进行Bacon分解，Andrew Goodman-Bacon（明尼阿波利斯联邦储备银行）、Thomas Goldring（乔治亚州立大学）和Austin Nichols（亚马逊）编写了bacondecomp包。基本语法如下：

bacondecomp Y D, ddtail

Y为因变量，D为处理虚拟变量，ddtail选项用于更详细的分解。

Stata 18（发布于2023年4月25日）引入了一个新的估计命令estat bdecomp，用于执行Bacon分解。它可用于didregress或xtdidregress命令之后，并且通过添加graph选项可以轻松创建图表。

遗憾的是，Stata中的Bacon分解仅在具有强平衡面板数据的情况下有效。对于具有不平衡面板数据的情况，为了进行Bacon分解并找到附加每个系数的权重，可以尝试由Clement de Chaisemartin、Xavier D'Haultfoeuille和Antoine Deeb（世界银行）编写的twowayfeweights。

平衡面板数据的合成双重差分（Synthetic DID）

Arkhangelsky等人（2022）提出了一种方法，即合成双重差分（SDID，前沿, 合成双重差分法SDID方法介绍和示例, 附code和数据！），它结合了双重差分（DID）和合成控制（SC）方法吸引人的特点。

与DID类似，SDID允许在所有处理前期期间处理组和对照组之间存在恒定的差异。

与SC类似，SDID通过重新加权和匹配处理前趋势来放宽传统的“平行趋势”假设。

SDID的关键步骤是找到个体权重（wi），以确保比较是在政策处理前近似遵循平行趋势的个体和对照组之间进行的，并找到时间权重（λt），给那些与政策处理后期间更相似的“政策处理前期间”更高的权重。通过将这些权重应用于标准双重差分（简单地对所有组和期间分配相等权重），我们可以得到政策处理的因果效应的一致估计。幸运的是，该方法可以应用于有交错政策处理的情况。不幸的是，在不平衡面板数据的情况下无法使用该方法，并且它无法估计动态效应。

在Stata中，可以使用sdid包（由Damian Clarke（Exeter大学）和Daniel Pailañir（智利大学）编写）来实施SDID。基本语法如下：

sdid depvar groupvar timevar treatment, vce(vcetype)

其中：depvar是因变量；groupvar是指示个体的变量，timevar是指示时间期间的变量；treatment指示在特定时间期间受到政策处理的个体。vce(vcetype)是一个必选选项，用于指定用于估计方差的方法。可用的推断方法包括自助法（bootstrap）、刀切法（jackknife）、安慰剂法（placebo）和无推断法（noinference）。如果只有一个个体接受处理，则无法使用自助法和刀切法。对于使用安慰剂法，我们需要至少一个比受处理个体更多的对照个体。

值得注意的是，可以使用graph和g1on选项创建图表（如Arkhangelsky等人，2022中的图1），显示个体权重（散点图）、时间权重（面积图）和结果趋势（线图）。

动态双重差分（Dynamic DID）

通常，研究人员对仅估计静态效应并不满意；他们可能还想看到政策的长期效应。例如，一旦Sibal Yang在Canvas上发布了新的问题集，那么在截止日期之前的每一天，问题集对学生幸福感的影响是什么？无论如何，经典的动态双重差分模型允许处理政策效应随时间而变化的情况。下面是动态双重差分模型的一个示例：

在动态模型中，研究人员需要解决多重共线性的问题。避免多重共线性的最常见方法是如上所示排除-1期间的处理虚拟变量（即政策处理前的最后一个期间）。此外，上面我对较远的相对期间进行了分组，这也是解决不平衡问题的常见方法。

在这种双重差分的格式中，要想取得无偏估计就要比经典双重差分复杂得多。许多计量经济学家试图解决这个问题，而且他们现在仍在努力。接下来，我只介绍几个全新的（动态的）双重差分估计值以及在Stata中使用它们的相应命令。请注意，以下某些Stata包还在开发中，因此在将其应用于研究之前，最好阅读它们的最新文档。

基于交互权重的双重差分估计（Interaction-Weighted Estimator for DID）

Sun＆Abraham（2021）为静态/动态固定效应TWFE模型中的系数识别提出了假设。具体而言，如果处理效应随时间变化（即βk随着k的变化），那么静态TWFE固定效应模型中的估计值将存在偏差，但在动态TWFE固定效应模型中，估计值在同质性假设下仍然有效（即βk在受处理组之间不会发生变化）。遗憾的是，当允许处理效应的异质性时，静态和动态TWFE固定效应模型中的估计值都存在偏差。

为了应对TWFE模型中的系数被污染问题，Sun＆Abraham（2021）提出了一种基于交互权重（IW）的估计方法。他们的估计值通过估计队列特定处理组的平均处理效应”（cohort-specific average treatment effect on the treated, CATT）的可解释加权平均值。在这里，队列被定义为同时首次接受政策处理的所有个体组成的群组。

其中一位作者，Liyang Sun（麻省理工学院），编写了一个Stata软件包eventstudyinteract，用于实施他们的IW估计方法并构建估计的置信区间。要使用eventstudyinteract命令，还必须安装另一个软件包：avar。基本语法如下：

eventstudyinteract y rel_time_list, absorb(id t) cohort(variable) control_cohort(variable) vce(vcetype)

请注意，我们必须包含一个相对时间指标列表，就像在经典的动态双重差分回归中一样。

不幸的是，这个命令与estout软件包不太兼容；因此，要在图表中报告结果，可能需要先将结果存储在一个矩阵中，然后再处理该矩阵。

DID的双重稳健估计值（Doubly Robust Estimator for DID）

Callaway和Sant'Anna（2021）特别关注分解的因果参数，即组别g在时间t的平均处理效应，其中“组别”是根据个体首次接受政策处理的时间段定义的。他们将该参数称为“组别-时间平均处理效应”，用ATT(g, t)表示，并提出了三种不同类型的DID估计方法来估计该参数：

结果回归（outcome regression，OR）；

逆概率加权（inverse probability weighting，IPW）；

双重稳健（doubly robust，DR）。

在估计中，Callaway和Sant'Anna建议使用DR方法，因为该方法只要求我们正确地指定比较组的结果或倾向得分模型中的一个（但不一定是两者）。

这两位作者与Fernando Rios-Avila合作编写了一个Stata软件包csdid，用于实施Callaway和Sant'Anna（2021）中提出的DID估计方法。在内部，所有2✖2的DID估计都是使用drdid命令获得的，因此要运行csdid，我们必须安装两个软件包：csdid和drdid。

基本语法如下：

csdid Y covar, ivar(id) time(t) gvar(group)

要运行这个模型，我们需要gvar变量：它为接受处理的个体首次接受政策处理的时间，对于未接受政策处理的个体等于0。请注意，该命令允许我们将协变量包括在回归中；在某些情况下，平行趋势假设可能仅在对观察到的协变量进行调整后成立。

该命令具有几种内置的估计系数方法；默认方法是dripw，即基于逆概率加权和普通最小二乘法的双重稳健DID估计值，该方法来自于Sant'Anna和Zhao（2020）的研究。可以使用method( )选项将其更改为其他可用的方法。此外，默认情况下，估计的是稳健和渐近标准误。然而，还有其他选项可供选择，例如使用wboot选项进行乘积Wild Bootstrap估计。在Stata命令窗口中输入help csdid以了解更多详细信息。

Callaway和Sant'Anna（2021）还提供了一些聚合方案，用于形成更加总体的因果参数。在Stata中，我们可以使用post-estimation estat或csdid_estat命令生成聚合估计结果。如果使用自助法程序（bootstrap procedures）估计标准误，推荐使用第二个命令。

双重差分的多组-时间估计值

自2020年以来，Clément de Chaisemartin（SciencePo）和Xavier D'Haultfœuille（CREST）撰写了一系列论文，提出了不同的双重差分估计技术。他们的主要贡献包括：

当处理效应随时间或跨群体呈现异质时，他们的估计值是有效的。

他们的估计值允许政策处理的转换（即允许个体退出）。当然，需要额外的假设；de Chaisemartin和D'Haultfœuille（2020）明确提出了关于Y(1)的强外生性、Y(1)的共同趋势以及始终在特定时间窗口中接受政策处理的稳定群组的存在。

他们的估计值考虑了面板数据中较晚出现的政策处理的折扣（受政策处理时间短，当然受到政策的影响相对就小一些）。这可能是经济学和统计学能够很好结合的证据。

他们提出了几种安慰剂估计值（通过模仿实际估计值构造），以检验“无预期”和“平行趋势”假设。

请注意，当使用尚未接受处理的群组作为对照组且回归中没有控制变量时，他们的估计值在数值上等价于上一节中介绍的Callaway和Sant'Anna（2021）的估计值。此外，请注意，de Chaisemar和D'Haultfœuille使用“staggered”一词来称呼不允许政策转换的处理设计。这与文献非常不一致；通常我们称这种政策处理设计为“吸收型处理”，而将“交错型处理”用于描述不同组别受到政策影响的时间存在差异的设计。

de Chaisemar和D'Haultfœuille编写了一个Stata软件包did_multiplegt，用于应用他们的估计器。基本语法如下：

did_multiplegt Y G T D

其中Y是结果变量，G是群组变量，T是时间变量，D是处理虚拟变量。该命令有很多选项，这里只介绍几个重要的选项：

如果指定了robust_dynamic选项，则将使用de Chaisemartin和D'Haultfoeuille（2022）中提出的估计值；否则将使用de Chaisemar和D'Haultfœuille（2020）中的估计值。如果要估计动态处理效应，则必须指定robust_dynamic。我们可以使用dynamic(#)选项指定要估计的动态效应的期数。

当指定了robust_dynamic时，Stata使用de Chaisemartin和D'Haultfoeuille（2022）中提出的长差分安慰剂（long difference placeboes）。我们可以使用firstdiff_placebo选项使Stata使用de Chaisemar和D'Haultfœuille（2020）中提出的一阶差分安慰剂。通过placebo(#)选项，可以指定要估计的安慰剂估计值的数量。该数量最多可以等于数据中的时间期数。

该命令通过使用自助法来估计标准误。我们可以使用cluster（）选项来要求Stata在特定级别上使用块自助法（block bootstrap）。不能直接在cluster（）中使用交互项，但可以在运行回归之前使用group（）函数生成交互项。

discount（#）选项允许我们对面板中后期估计的处理效应进行折扣。

该命令默认在每个回归之后生成图形。通过graphoptions（）选项，我们可以修改图形的外观。

这两位作者为他们的软件包编写了一个信息丰富的文档（包含一个长的常见问题解答部分）。可以在命令窗口中输入help did_multiplegt来查看文档。

基于插补的双重差分估计值（Imputation Estimator for DID）

Borusyak、Jaravel和Spiess（2023）提出了一种使用插补过程的有限样本高效稳健双重差分估计值。插补方法受到欢迎的原因包括：

它在计算上高效（只需要估计一个简单的TWFE模型）；

插补方法可以轻松地将平行趋势和无预期假设与估计值联系起来。

其中一位作者，Kirill Borusyak（伦敦大学学院），编写了一个Stata软件包did_imputation，用于实施他们的插补方法，以估计动态处理效应并进行事件研究中的前期趋势检验。基本语法如下：

did_imputation Y id t Ei, fe(id t) horizons(#) pretrends(#)

其中，horizons选项告诉Stata我们希望估计多少个前瞻期的处理效应，而pretrends选项告诉Stata要对某些时期进行前期趋势检验。前期趋势系数报告为pre1、pre2等。与前面提到的方法不同，这里的前期趋势系数的数量不会影响处理效应估计后的其他估计（他们始终在平行趋势和无预期假设下计算）。

此外，Borusyak、Jaravel和Spiess（2022）是一篇出色的论文，指出了经典双重差分中臭名昭著的“负权重”问题。这个问题的出现是因为OLS估计对处理效应的同质性施加了非常强的限制。这就是为什么有些计量经济学家称经典的动态双重差分为一个受污染的估计值。

待续...

潜在候选论文：Dube等人（2023）。

实例

在本节中，我将展示如何在实证分析中使用上述估计方法，特别是使用Stata中的特定命令/软件包。

时间固定效应（TWFE）与合成双重差分（SDID）

在这里，我将使用三个真实世界的数据集来展示如何运行TWFE回归。如果可行，我还将展示运行SDID的代码，然后进行比较。

我将使用的数据来自三篇论文：

"OW05_prop99.dta"来自Orzechowski和Walker（2005），你可以从这里获取最新版本。Abadie等人（2010）和Arkhangelsky等人（2021）使用这些数据来估计提案99（增加香烟税）对加利福尼亚州人均销售的香烟的影响。请注意，这是一个分块政策处理（block treatment）的案例。

"BCGV22_gender_quota.dta"来自Bhalotra等人（2022）。他们使用这些数据来估计议会性别配额（在议会中为妇女保留席位）对议会中妇女比例的影响。请注意，这是一个交错处理的案例。

"SZ18_state_taxes.dta"来自Serrato和Zidar（2018）。他们使用这些数据来估计州企业税收减税/增税对税收和经济活动的影响。请注意，这是一个交错处理的案例；然而，Serrato和Zidar（2018）使用动态标准差分模型（即，未解决负权重问题），因此他们的结果可能存在偏差。此外，请注意，不幸的是sdid命令无法运行动态模型，因此我们无法使用SDID来更新Serrato和Zidar（2018）的结果。

我将使用的回归命令包括：

xtdidregress，Stata内置命令，用于在面板数据上运行双重差分回归。使用此命令后，我们可以使用estat创建趋势图并进行一些基本检验。

xtreg、areg和reghdfe最初用于运行固定效应模型，但也可以轻松应用于运行双重差分回归。

sdid，用于运行SDID的外部命令。通过method()选项，我们还可以使用它来运行标准双重差分和合成控制模型。

xthdidregress，是Stata 18引入的一种用于估计异质性ATT的命令。请注意，xthdidregress命令允许多种加权方式，我选择使用aipw（增强逆概率加权，也称为“双重稳健”）。尽管Stata 18的文档没有明确说明此命令是由哪位作者提出，但我猜测这个命令的思想源于Callaway和Sant'Anna（2021）。请注意，回归后，estat aggregation允许我们在同一组或时间内汇总ATT，分析指定时间窗口内的动态效应，并创建图表。

在此处可以找到在这三个数据集上运行回归的完整代码。

下面这些代码经过测试，可以直接在Stata运行，数据附在文后网盘里（网盘里的代码1）。

cd "C:\Users\xiwan\Desktop\DataforDID"  

**# Block Treatment: Orzechowski & Walker (2005)

use "OW05_prop99.dta", clear

encode state, gen(state_code)

xtset state_code year

* TWFE DID

xtdidregress (packspercapita) (treated), group(state_code) time(year) vce(cluster state_code)

estat trendplots // visualization

estat ptrends // parallel-trends test

estat granger // anticipation test

*estat grangerplot // time-specific treatment effects

eststo twfe1: qui xtreg packspercapita treated i.year, fe cluster(state_code)

eststo twfe2: qui areg packspercapita treated i.year, absorb(state_code) cluster(state_code)

eststo twfe3: qui reghdfe packspercapita treated, absorb(state_code year) cluster(state_code)

estout twfe*, keep(treated) ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N r2_a, nostar labels("Observations" "Adjusted R2") fmt("%9.0fc"3))

* Sythetic DID

eststo syn_did: sdid packspercapita state year treated, vce(placebo) seed(1) ///

graph g1on ///

g1_opt(xtitle("") ///

plotregion(fcolor(white) lcolor(white)) ///

graphregion(fcolor(white) lcolor(white)) ///

) ///

g2_opt(ylabel(0(25)150) ytitle("Packs per capita") ///

plotregion(fcolor(white) lcolor(white)) ///

graphregion(fcolor(white) lcolor(white)) ///

) ///

graph_export("prop99_did_", .pdf)

ereturn list

matrix list e(omega) // unit-specific weights

matrix list e(lambda) // time-specific weights

estout twfe* syn_did, keep(treated) ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N r2_a, nostar labels("Observations" "Adjusted R2") fmt("%9.0fc"3))

* Compare SDID, DID, and SC

foreach m in sdid did sc {

if "`m'"=="did" {

local g_opt = "msize(small)"

}

else {

local g_opt = ""

}

sdid packspercapita state year treated, ///

method(`m') vce(noinference) graph g1on `g_opt' ///

g1_opt(xtitle("") ytitle("") ///

xlabel(, labsize(tiny)) ///

ylabel(-100(25)50, labsize(small)) ///

plotregion(fcolor(white) lcolor(white)) ///

graphregion(fcolor(white) lcolor(white)) ///

) ///

g2_opt(title("`m'") xtitle("") ytitle("") ///

xlabel(, labsize(small)) ///

ylabel(0(25)150, labsize(small)) ///

plotregion(fcolor(white) lcolor(white)) ///

graphregion(fcolor(white) lcolor(white)) ///

)

graph save g1_1989 "`m'_1.gph", replace

graph save g2_1989 "`m'_2.gph", replace

}

graph combine "sdid_2.gph" "did_2.gph" "$figdir\sc_2.gph" "sdid_1.gph" "did_1.gph" "sc_1.gph", ///

cols(3) xsize(3.5) ysize(2) ///

graphregion(fcolor(white) lcolor(white))

graph export "compare_sdid_did_sc.pdf", replace

*****************************************************

**# Staggered Treatment: Bhalotra et al. (2022)

use "BCGV22_gender_quota.dta", clear

drop if lngdp==.

tab year

isid country year

* Balanced panel: 115 countries, years from 1990 to 2015.

* Treated group: 9 countries.

* Control group: 106 other states.

eststo stagg1: sdid womparl country year quota, vce(bootstrap) seed(3) ///

graph g1on ///

g1_opt(xtitle("") ytitle("") xlabel(, labsize(tiny)) ///

plotregion(fcolor(white) lcolor(white)) ///

graphregion(fcolor(white) lcolor(white)) ///

) ///

g2_opt(ytitle("") ///

plotregion(fcolor(white) lcolor(white)) ///

graphregion(fcolor(white) lcolor(white)) ///

)

ereturn list

matrix list e(tau) // adoption-period specific estimate

matrix list e(lambda)

matrix list e(omega)

matrix list e(adoption) // different adoption years

eststo stagg2: sdid womparl country year quota, covariates(lngdp, optimized) vce(bootstrap) seed(3)

eststo stagg3: sdid womparl country year quota, covariates(lngdp, projected) vce(bootstrap) seed(3)

estout stagg*, ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N, nostar labels("Observations") fmt("%9.0fc"))

******************************************************************

**# Staggered Treatment: Serrato & Zidar (2018)

use "SZ18_state_taxes.dta", replace

keep if year >= 1980 & year <= 2010

drop if fips_state == 11 | fips_state == 0 | fips_state > 56

xtset fips_state year

* Balanced panel: 50 states, years from 1980 to 2010.

* Treated group: 15 states.

* Control group: 35 other states.

* Construct dependent variables

gen log_rev_corptax = 100*ln(rev_corptax+1)

gen log_gdp = 100*ln(GDP)

g r_g = 100*rev_corptax/GDP

* Screen out the tax change with a pre-determined threshold

local threshold = 0.5

gen ch_corporate_rate = corporate_rate - L1.corporate_rate

replace ch_corporate_rate = 0 if abs(ch_corporate_rate) <= `threshold'

gen ch_corporate_rate_inc = (ch_corporate_rate > 0 & !missing(ch_corporate_rate))

gen ch_corporate_rate_dec = (ch_corporate_rate < 0 & !missing(ch_corporate_rate))

* Contruct treatment dummies

** Static

gen change_year = year if ch_corporate_rate_dec==1

bysort fips_state: egen tchange_year = min(change_year)

gen treated = (year >= tchange_year)

** Dynamic

gen period = year - tchange_year

gen Dn5 = (period < -4)

forvalues i = 4(-1)1 {

gen Dn`i' = (period == -`i')

}

forvalues i = 0(1)5 {

gen D`i' = (period == `i')

}

gen D6 = (period >= 6 & period != .)

* Static DID

local ylist = "log_rev_corptax log_gdp r_g"

local i = 1

foreach yvar in `ylist' {

quietly{

eststo areg`i': areg `yvar' treated i.fips_state, absorb(year) cluster(fips_state)

eststo hdreg`i': reghdfe `yvar' treated, absorb(fips_state year) cluster(fips_state)

eststo sdid`i': sdid `yvar' fips_state year treated, vce(bootstrap) seed(2018)

local ++i

}

}

estout areg1 hdreg1 sdid1 areg2 hdreg2 sdid2 areg3 hdreg3 sdid3, keep(treated) ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N r2_a, nostar labels("Observations" "Adj. R2") fmt("%9.0fc" 3))

* Dynamic DID (sdid cannot do this)

local ylist = "log_rev_corptax log_gdp r_g"

local i = 1

foreach yvar in `ylist' {

quietly{

eststo areg`i': areg `yvar' Dn5-Dn2 D0-D6 i.fips_state, absorb(year) cluster(fips_state)

eststo hdreg`i': reghdfe `yvar' Dn5-Dn2 D0-D6, absorb(fips_state year) cluster(fips_state)

local ++i

}

}

estout areg1 hdreg1 areg2 hdreg2 areg3 hdreg3, keep(D*) ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N r2_a, nostar labels("Observations" "Adj. R2") fmt("%9.0fc" 3))

**# Stata 18 new command: xthdidregress

local controls = "FederalIncomeasStateTaxBase sales_wgt throwback FedIncomeTaxDeductible Losscarryforward FranchiseTax"

local ps_var = "FederalIncomeasStateTaxBase sales_wgt throwback FedIncomeTaxDeductible FranchiseTax"

xthdidregress aipw (log_gdp `controls') (treated `ps_var'), group(fips_state) vce(cluster fips_state)

* Visualizing ATT for each cohort

estat atetplot

graph export "SZ18_cohort_ATT.pdf", replace

* Visualizing ATT over cohort

estat aggregation, cohort ///

graph(xlab(, angle(45) labsize(small)) legend(rows(1) position(6)))

graph export "SZ18_cohort_agg_ATT.pdf", replace

* Visualizing ATT over time

estat aggregation, time ///

graph(xlab(, angle(45) labsize(small)) legend(rows(1) position(6)))

graph export "SZ18_time_agg_ATT.pdf", replace

* Visualizing dynamic effects

estat aggregation, dynamic(-5/6) ///

graph( ///

title("Dynamic Effects on Log State GDP", size(medlarge)) ///

xlab(, labsize(small) nogrid) ///

legend(rows(1) position(6)) ///

xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///

)

graph export "SZ18_dynamic_ATT.pdf", replace

---

Bacon分解、事件研究图和安慰剂检验

可以通过运行以下代码将要使用的数据集加载到Stata中：

use "http://pped.org/bacon_example.dta", clear

面板数据包含1964年至1996年美国49个州（包括华盛顿特区，但不包括阿拉斯加和夏威夷）的州级信息（特别是无过错离婚开始年份和自杀死亡率）。这些数据最初由Stevenson＆Wolfers（2006）用于估计无过错（或单方面）离婚对女性自杀率的影响。

这里，首先做一个无过错离婚改革（一种分散处理）对女性自杀进行静态的双重差分估计：

处理组包括采用单方面离婚法的州，而对照组包括其余州。

所有以下命令的估计系数应相同（但由于不同的算法，标准误差和R平方不同）。

xtdidregress (asmrs pcinc asmrh cases) (post), group(stfips) time(year) vce(cluster stfips)

xtreg asmrs post pcinc asmrh cases i.year, fe vce(cluster stfips)

areg asmrs post pcinc asmrh cases i.year, absorb(stfips) vce(cluster stfips)

reghdfe asmrs post pcinc asmrh cases, absorb(stfips year) cluster(stfips)

其中，asmrs是自杀死亡率，post是政策处理虚拟变量。其他变量都是控制变量。Stata报告的双重差分系数为-2.516（标准误差为2.283），与零无显著差异。

然后，我们可以对TWFE DID模型应用Bacon分解定理。

bacondecomp asmrs post pcinc asmrh cases, ddetail

下面这个Bacon分解结果为新版本bacondecomp得出的。其中，Timing-groups表示不同政策处理时间带来的差异，它占据了总效应的 37.76%，对应的是旧版本bacondecomp中的Earlier T vs. Later C 和Later T vs. Earlier C两部分的权重，Always_v_timing的权重为37.8%，对应的是旧版本bacondecomp中的T vs. Already treated，即新处理组与已处理组的对比差异。Never_v_timing的权重为23.9%，对应的是旧版本bacondecomp中的T vs. Never treated，即处理组对应着从未接受处理过的对照组的差异。而同组内的差异占比为 0.509%，始终处理与从未接受处理组的效应仅占 0.0018%。

要看偏差，主要是看T vs. Already treated或新版本中的Always_v_timing的权重。

它报告数据集中有14个时间组（在不同时间接受处理可以作为彼此的对照组），包括一个从未处理的组和一个始终处理的组。最大的权重分配给始终处理组和时间组之间的比较。

我们还可以使用以下编码（在xtdidregress之后）进行分解。散点图可以在这里找到。

estat bdecomp, graph

请记住，Bacon分解作为一种诊断工具，而不是一种修正方法。分解告诉我们DID模型中“不好的组别的比较”的严重程度，但它不能解决此问题。

接下来，我进行的是相应的动态DID回归。我使用eventdd命令，因为它可以同时运行模型并生成图表。该命令允许进行一些基本的回归（例如xtreg和reghdfe），对于从高级DID回归中绘制结果，我推荐使用event_plot包（将在下一个示例中详细介绍）。要使用eventdd，必须安装两个软件包，即eventdd和matsort。

最后，我通过随机选择安慰剂处理时间并重复进行1000次TWFE回归来进行时间安慰剂检验。这项工作是使用Stata内置命令permute完成的。检验结果显示，我上述的估计可能不是来自于一个不可观察的时间趋势。安慰剂检验的图表（由附加命令dpplot创建）可以在这里找到。

完整的编码可以在这里找到。

下面这些代码经过测试，可以直接在Stata运行（网盘里的代码2）。

*** Bacon Decomposition *************************

use "http://pped.org/bacon_example.dta", clear

* We see multiple treatment years.

tab year

tab _nfd // the year of the passage of no-fault divorce law

* Regression of female suicide on no-fault divorce reforms

eststo reg1: quietly reghdfe asmrs post pcinc asmrh cases, absorb(stfips year) cluster(stfips)

estout reg1, keep(post pcinc asmrh cases _cons) ///

varlabels(_cons "Constant") ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) ///

stats(N r2_a, nostar labels("Observations" "R-Square") fmt("%9.0fc" 3))

bacondecomp asmrs post pcinc asmrh cases, ddetail

graph export "DID_Decomposition_bacondecomp.pdf", replace

xtdidregress (asmrs pcinc asmrh cases) (post), group(stfips) time(year) vce(cluster stfips)

estat bdecomp, graph

graph export "DID_Decomposition_estat.pdf", replace

*** Event Study Plots *****************************************************

gen rel_time = year - _nfd

* The following uses "xtreg".

eventdd asmrs pcinc asmrh cases i.year, ///

timevar(rel_time) method(fe, cluster(stfips)) ///

noline graph_op( ///

xlabel(-20(5)25, nogrid) ///

xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///

legend(order(1 "Point Estimate" 2 "95% CI") size(*0.8) position(6) rows(1) region(lc(black))) ///

)

* The following uses "reghdfe".

eventdd asmrs pcinc asmrh cases, ///

timevar(rel_time) method(hdfe, cluster(stfips) absorb(stfips year)) ///

noline graph_op( ///

xlabel(-20(5)25, nogrid) ///

xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///

legend(order(1 "Point Estimate" 2 "95% CI") size(*0.8) position(6) rows(1) region(lc(black))) ///

)

* Only balanced periods in which all units have data are shown in the plot.

eventdd asmrs pcinc asmrh cases i.year, ///

timevar(rel_time) method(hdfe, cluster(stfips) absorb(stfips year)) balanced ///

noline graph_op( ///

xlabel(, nogrid) ///

xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///

legend(order(1 "Point Estimate" 2 "95% CI") size(*0.8) position(6) rows(1) region(lc(black))) ///

)

* Only specified periods are shown in the plot; periods beyond the window are accumulated.

eventdd asmrs pcinc asmrh cases i.year, ///

timevar(rel_time) method(hdfe, cluster(stfips) absorb(stfips year)) ///

accum leads(5) lags(10) ///

noline graph_op( ///

xlabel(, nogrid) ///

xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///

legend(order(1 "Point Estimate" 2 "95% CI") size(*0.8) position(6) rows(1) region(lc(black))) ///

)

*** Placebo Test *****************************************

* Randomly select a placebo treatment time and run TWFE regression for 1000 times

permute post coefficient=_b[post], reps(1000) seed(1) saving("placebo_test.dta", replace): reghdfe asmrs post pcinc asmrh cases, absorb(stfips year) cluster(stfips)

use "placebo_test.dta", clear

dpplot coefficient, ///

xline(-2.516, lc(red) lp(dash)) xline(0, lc(gs12) lp(solid)) ///

xtitle("Effect estimate") ytitle("Density of distribution") ///

xlabel(-3(1)2 -2.516, nogrid labsize(small)) ///

ylabel(, labsize(small)) caption("")

graph export "placebo_test_plot.pdf", replace

交错处理的动态 DID

如果一家公司以低于其在国内市场上通常销售的价格出口产品，我们称该公司在倾销产品。这种不公平的外国定价行为可能对进口市场的企业和经济产生不利的扭曲作用。考虑到这种不利影响，WTO反倾销协议允许政府采取一些反倾销行动来应对外国的倾销行为。最典型的反倾销行动是对特定产品从特定出口国征收较高的进口关税，希望将不公平的低价提高到正常价格，从而减轻对进口国的损害。

在这里，我将使用动态DID模型来估计美国对中国出口商实施的反倾销关税在2000年到2009年的动态效应。Bown和Crowley（2007）称这种效应为“贸易破坏”，并通过使用美国和日本的数据运行IV、FE和GMM模型来估计静态政策效应。有关该论文的文献综述可以在这里找到。

我将使用一个“产品-年”层面的数据集，该数据集是从全球反倾销数据库和中国海关数据合并而来（感谢清华大学中国数据中心）。我将运行的动态DID模型如下：

除非另有说明，否则标准误差在“产品-年”的层面进行聚类。

处理组是来自中国的那些受到美国反倾销关税的产品，而对照组是来自中国的那些经历了反倾销调查但最终未受到反倾销关税的产品。请注意，我没有将那些从未接受调查的产品纳入对照组。原因是反倾销调查是非随机的；受调查的产品的出口价格始终低于没有接受调查的产品。如果我将受到反倾销关税的产品与未经调查的产品进行比较，那么我的估计很可能存在偏误。

关于对特定产品（在6位数级别上进行编码）实施反倾销关税的年份的信息存储在变量year_des_duty中。我使用这个变量和年份来构建一系列相对时间虚拟变量：

gen period_duty = year - year_des_duty

gen Dn3 = (period_duty < -2)

forvalues i = 2(-1)1 {

gen Dn`i' = (period_duty == -`i')

}

forvalues i = 0(1)3 {

gen D`i' = (period_duty == `i')

}

gen D4 = (period_duty >= 4) & (period_duty != .)

传统上，研究人员使用政策处理前的系数来检验预趋势（pre-trends）：如果处理前的系数与0没有显著差异，那么他们会得出平行趋势假设成立的结论。然而，Sun和Abraham（2021）证明了这种做法存在严重缺陷，并需要进行修正。

使用双重差分的交互加权估计的（interaction-weighted estimation of DID）编码如下：

gen first_union = year_des_duty

gen never_union = (first_union == .)

local dep = "value quantity company_num m_quantity"

foreach y in `dep'{

eventstudyinteract ln_`y' Dn3 Dn2 D0-D4, \\\

cohort(first_union) control_cohort(never_union) \\\

absorb(product year) vce(cluster product#year)

}

我们需要告诉Stata每个个体的初始处理时间对应的变量是哪个，将其命名为first_union。对于从未接受政策处理的个体，该变量应设置为缺失。此外，我们还需要提供一个虚拟变量，对应于对照组，可以是从未接受处理的个体或最后接受处理的个体。在这里，我将从未接受处理的个体作为对照组，并构建一个名为never_union的变量来表示它。

值得注意的是，Kirill Borusyak编写了一个名为event_plot的包，用于方便地绘制双重差分估计结果，包括后处理系数和（如果有的话）预趋势系数，以及置信区间。我使用这个命令创建了四张图，他们展示了四个结果变量的动态效应。对于绘图，我通常会自定义绘图类型，但实际上，如果你对可视化的要求不像我这么高，可以节省很多时间，使用默认类型（使用default_look选项）。

双重稳健估计（doubly robust estimation of DID）的编码如下：

gen gvar = year_des_duty

recode gvar (. = 0)

local dep = "value quantity company_num m_quantity"

foreach y in `dep'{

quietly csdid ln_`y', ivar(product) time(year) gvar(gvar) \\\

method(dripw) wboot(reps(10000)) rseed(1)

csdid_estat event, window(-3 4) estore(cs_`y') wboot(reps(10000)) rseed(1)

}

cs_did命令在Stata结果窗口中可能显示一个非常长的输出表格，因此我在csdid之前加上了quietly命令。此外，在许多应用中，我更关心不同时间点的异质性处理效应，而不是跨不同组的效应（而不是由Callaway＆Sant'Anna, 2021定义的组-时间平均处理效应）；因此，我使用csdid_estat仅在-3到4期间产生汇总估计结果。现在结果窗口中的输出表格较短。还值得注意的是，我使用wboot选项来估计Wild Bootstrap标准误差，重复10,000次。

与之前一样，我使用event_plot命令创建四面张图来展示动态效应。这次，我使用default_look选项节省时间；此外，我使用together选项将前后滞后项显示为一个连续的曲线。

de Chaisemartin和D'Haultfoeuille的DID估计的编码如下：

egen clst = group(product year)

local ylist = "value quantity m_quantity company_num"

foreach y in `ylist'{

did_multiplegt ln_`y' year_des_duty year treated, ///

robust_dynamic dynamic(4) placebo(2) jointtestplacebo ///

seed(1) breps(100) cluster(clst)

}

这里使用de Chaisemartin和D'Haultfoeuille（2022）的估计方法，因为我想估计动态效应。值得注意的是，我使用长差分安慰剂进行安慰剂检验；动态效应是使用长差分的双重差分估计值来估计的，因此使用长差分安慰剂是正确。相比之下，一阶差分安慰剂估计值是连续时间段之间的双重差分；如果在这里添加firstdiff_placebo，用于说明动态处理效应的图形将是无意义的（即不可比较的）。

相关讨论可以在这里找到（https://www.statalist.org/forums/forum/general-stata-discussion/general/1599964-graph-for-the-dynamic-treatment-effect-using-did_multiplegt-package）。

我个人不喜欢did_multiplegt自动生成的图形，而且不幸的是，graphoptions( )选项不如我期望的灵活。因此，我选择撤回并将估计结果存储在矩阵中。然后，可以使用event_plot命令创建图形。创建图形的过程与我在应用Sun＆Abraham（2021）的估计值时所做的非常相似。我的四张图形展示如下：

imputation DID估计的编码如下：

gen id = product

gen Ei = year_des_duty

local ylist = "value quantity company_num m_quantity"

foreach y in `ylist'{

did_imputation ln_`y' id year Ei, \\\

fe(product year) cluster(clst) horizons(0/4) pretrends(2) minn(0) autosample

}

我们需要为个体的处理日期提供一个变量，其中缺失值表示从未接受处理的个体。我按照软件文档将其命名为Ei。

可以在这里找到使用插补方法估计的动态效应的四张图形，具体如下。这次，我使用了default_look选项，但不使用together选项，所以前瞻和滞后期显示为两条分别用不同颜色表示的曲线。

总结一下，无论采用哪种估计方法，结果都显示出美国反倾销税征收对四个结果变量产生了持续且负面的影响。

下面这些代码经过测试，可以直接在Stata运行（网盘里的代码3）。

cd "C:\Users\xiwan\Desktop\DataforDID"

clear all

*************************************************

use "CCD_GAD_USA_affirm.dta", clear

**# Create a series of dummies for duty impositions

gen period_duty = year - year_des_duty

gen treated = (period_duty < . & period_duty >= 0)

gen Dn3 = (period_duty < -2)

forvalues i = 2(-1)1 {

gen Dn`i' = (period_duty == -`i')

}

forvalues i = 0(1)3 {

gen D`i' = (period_duty == `i')

}

gen D4 = (period_duty >= 4) & (period_duty != .)

**# Classical Dynamic DID

* install reghdfe ftools

encode hs06, generate(product) // string variables not allowed for xtset and regression

local dep = "value quantity company_num m_quantity"

foreach y in `dep'{

quietly{

eststo reg_`y': reghdfe ln_`y' Dn3 Dn2 D0-D4, absorb(product year) vce(cluster product#year)

}

}

estout reg*, keep(D*) ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N r2_a, nostar labels("Observations" "Adjusted R-Square") fmt("%9.0fc" 3))

******************************************************

**# Sun and Abraham (2021)

* install eventstudyinteract avar

gen first_union = year_des_duty

gen never_union = (first_union == .)

local dep = "value quantity company_num m_quantity"

foreach y in `dep'{

* Regression

eventstudyinteract ln_`y' Dn3 Dn2 D0-D4, cohort(first_union) control_cohort(never_union) absorb(product year) vce(cluster product#year)

* Visualization

if "`y'"=="value"{

local panel = "A)"

local title = "ln(Value)"

}

if "`y'"=="quantity"{

local panel = "B)"

local title = "ln(Total Quantity)"

}

if "`y'"=="company_num"{

local panel = "C)"

local title = "ln(Number of Exporters)"

}

if "`y'"=="m_quantity"{

local panel = "D)"

local title = "ln(Mean Quantity)"

}

forvalue i=1/7 {

local m_`i' = e(b_iw)[1,`i']

local v_`i' = e(V_iw)[`i',`i']

}

matrix input matb_sa= (`m_1',`m_2',0,`m_3',`m_4',`m_5',`m_6',`m_7')

mat colnames matb_sa= ld3 ld2 ld1 lg0 lg1 lg2 lg3 lg4

matrix input mats_sa= (`v_1',`v_2',0,`v_3',`v_4',`v_5',`v_6',`v_7')

mat colnames mats_sa= ld3 ld2 ld1 lg0 lg1 lg2 lg3 lg4

event_plot matb_sa#mats_sa, ///

stub_lag(lg#) stub_lead(ld#) ///

ciplottype(rcap) plottype(scatter) ///

lag_opt(msymbol(D) mcolor(black) msize(small)) ///

lead_opt(msymbol(D) mcolor(black) msize(small)) ///

lag_ci_opt(lcolor(black) lwidth(medthin)) ///

lead_ci_opt(lcolor(black) lwidth(medthin)) ///

graph_opt( ///

title("`panel' `title'", size(medlarge) position(11)) ///

xtitle("Period", height(5)) xsize(5) ysize(4) ///

ytitle("Coefficient", height(5)) ///

xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///

yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) ///

xlabel(-3 "< -2" -2(1)3 4 "> 3", labsize(small)) ///

ylabel(, labsize(small)) ///

legend(order(1 "Coefficient" 2 "95% CI") size(*0.8)) ///

name(sa_`y', replace) ///

graphregion(color(white)) ///

)

}

* install grc1leg: net install grc1leg.pkg, replace

grc1leg sa_value sa_quantity sa_company_num sa_m_quantity, ///

legendfrom(sa_value) cols(2) ///

graphregion(fcolor(white) lcolor(white)) ///

name(sa_fig, replace)

gr draw sa_fig, ysize(5) xsize(6.5)

graph export "SA_DID_Trade_Destruction.pdf", replace

***************************************************************

**# Callaway & Sant'Anna (2021)

* install csdid drdid

gen gvar = year_des_duty

recode gvar (. = 0)

local dep = "value quantity company_num m_quantity"

foreach y in `dep'{

quietly csdid ln_`y', ivar(product) time(year) gvar(gvar) method(dripw) wboot(reps(10000)) rseed(1)

csdid_estat event, window(-3 4) estore(cs_`y') wboot(reps(10000)) rseed(1)

}

estout cs_*, keep(T*) ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N r2_a, nostar labels("Observations" "Adjusted R-Square") fmt("%9.0fc" 3))

* Visualization

local ylist = "value quantity company_num m_quantity"

foreach y in `ylist'{

if "`y'"=="value"{

local panel = "A)"

local title = "ln(Value)"

}

if "`y'"=="quantity"{

local panel = "B)"

local title = "ln(Total Quantity)"

}

if "`y'"=="company_num"{

local panel = "C)"

local title = "ln(Number of Exporters)"

}

if "`y'"=="m_quantity"{

local panel = "D)"

local title = "ln(Mean Quantity)"

}

event_plot cs_`y', default_look ///

stub_lag(Tp#) stub_lead(Tm#) together ///

graph_opt( ///

xtitle("Period") ytitle("ATT") ///

title("`panel' `title'", size(medlarge) position(11)) ///

xlab(-3(1)4, labsize(small)) ///

ylab(, angle(90) nogrid labsize(small)) ///

legend(lab(1 "Coefficient") lab(2 "95% CI") size(*0.8)) ///

name(cs_`y', replace) ///

)

}

grc1leg cs_value cs_quantity cs_company_num cs_m_quantity, ///

legendfrom(cs_value) cols(2) ///

graphregion(fcolor(white) lcolor(white)) ///

name(cs_fig, replace)

gr draw cs_fig, ysize(5) xsize(6.5)

graph export "CS_DID_Trade_Destruction.pdf", replace

*********************************************************

**# de Chaisemartin & D'Haultfoeuille (2020, 2022)

* install did_multiplegt

egen clst = group(product year)

local dep = "value quantity company_num m_quantity"

foreach y in `dep'{

* Regression

did_multiplegt ln_`y' year_des_duty year treated, ///

robust_dynamic dynamic(4) placebo(2) jointtestplacebo ///

seed(1) breps(100) cluster(clst)

* Visualization

if "`y'"=="value"{

local panel = "A)"

local title = "ln(Value)"

}

if "`y'"=="quantity"{

local panel = "B)"

local title = "ln(Total Quantity)"

}

if "`y'"=="company_num"{

local panel = "C)"

local title = "ln(Number of Exporters)"

}

if "`y'"=="m_quantity"{

local panel = "D)"

local title = "ln(Mean Quantity)"

}

forvalue i=1/8 {

local m_`i' = e(estimates)[`i',1]

local v_`i' = e(variances)[`i',1]

}

matrix input matb_DIDl= (`m_1',`m_2',`m_3',`m_4',`m_5',0,`m_7',`m_8')

mat colnames matb_DIDl= lg0 lg1 lg2 lg3 lg4 ld1 ld2 ld3

matrix input mats_DIDl= (`v_1',`v_2',`v_3',`v_4',`v_5',0,`v_7',`v_8')

mat colnames mats_DIDl= lg0 lg1 lg2 lg3 lg4 ld1 ld2 ld3

event_plot matb_DIDl#mats_DIDl, ///

stub_lag(lg#) stub_lead(ld#) ///

ciplottype(rcap) plottype(scatter) ///

lag_opt(msymbol(D) mcolor(black) msize(small)) ///

lead_opt(msymbol(D) mcolor(black) msize(small)) ///

lag_ci_opt(lcolor(black) lwidth(medthin)) ///

lead_ci_opt(lcolor(black) lwidth(medthin)) ///

graph_opt( ///

title("`panel' `title'", size(medlarge) position(11)) ///

xtitle("Period", height(5)) xsize(5) ysize(4) ///

ytitle("Average Effect", height(5)) ///

xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) ///

yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) ///

xlabel(-3/4, nogrid labsize(small)) ///

ylabel(, labsize(small)) ///

legend(order(1 "Coefficient" 2 "95% CI") size(*0.8) position(6) rows(1) region(lc(black))) ///

name(DIDl_`y', replace) ///

graphregion(color(white)) ///

)

}

grc1leg DIDl_value DIDl_quantity DIDl_company_num DIDl_m_quantity, ///

legendfrom(DIDl_value) cols(2) ///

graphregion(fcolor(white) lcolor(white)) ///

name(DIDl_fig, replace)

gr draw DIDl_fig, ysize(5) xsize(6.5)

graph export "CD_DIDl_Trade_Destruction.pdf", replace

*****************************************************

**# Borusyak, Jaravel & Spiess (2022)

* install did_imputation

gen id = product

gen Ei = year_des_duty

local ylist = "value quantity company_num m_quantity"

foreach y in `ylist'{

quietly{

eststo imp_`y': did_imputation ln_`y' id year Ei, fe(product year) cluster(clst) horizons(0/4) pretrends(2) minn(0) autosample

}

}

estout imp*, ///

coll(none) cells(b(star fmt(3)) se(par fmt(3))) ///

starlevels(* .1 ** .05 *** .01) legend ///

stats(N r2_a, nostar labels("Observations" "Adjusted R-Square") fmt("%9.0fc" 3))

* Visualization

* install event_plot

local ylist = "value quantity company_num m_quantity"

foreach y in `ylist'{

if "`y'"=="value"{

local panel = "A)"

local title = "ln(Value)"

}

if "`y'"=="quantity"{

local panel = "B)"

local title = "ln(Total Quantity)"

}

if "`y'"=="company_num"{

local panel = "C)"

local title = "ln(Number of Exporters)"

}

if "`y'"=="m_quantity"{

local panel = "D)"

local title = "ln(Mean Quantity)"

}

event_plot imp_`y', default_look ///

graph_opt( ///

xtitle("Period") ytitle("Coefficient estimate") ///

title("`panel' `title'", size(medlarge) position(11)) ///

xlab(-2(1)3 4 "> 3", labsize(small)) ///

ylab(, angle(90) nogrid labsize(small)) ///

yline(0, lcolor(gs8) lpattern(dash)) ///

legend(size(*0.8)) ///

name(imp_`y', replace) ///

)

}

grc1leg imp_value imp_quantity imp_company_num imp_m_quantity, ///

legendfrom(imp_value) cols(2) ///

graphregion(fcolor(white) lcolor(white)) ///

name(imp_fig, replace)

gr draw imp_fig, ysize(5) xsize(6.5)

graph export "Imputation_DID_Trade_Destruction.pdf", replace

在Stata上就是这样，可以直接使用。

数据和代码链接: 

https://pan.baidu.com/s/1uA5b2xxJ9gADU6kPtQZ9-A

提取码: qy2q 

---

关于交叠或多期DID，参看：1.120篇DID双重差分方法的文章合集, 包括代码,程序及解读, 建议收藏！2.诚实双重差分法DID, 面板事件研究法和Bacon分解的经典应用文！3.前沿: 多期或渐进或交叠DID, 如何进行平行趋势检验呢？4.多期DID或渐进DID或交叠DID, 最新Stata执行命令整理如下供大家学习，5.DID前沿: 5种方法估计事件研究的因果效应, 并使用绘制系数和置信区间, 详细代码和数据，6.事件研究法开展政策评估和因果识别, 分享8篇提供数据和代码的文章，7.推荐用渐进(多期)DID和事件研究法开展政策评估的论文及其实现数据和代码！8.机器学习已经与政策评估方法, 例如事件研究法结合起来识别政策因果效应了！9.前沿, 模糊双重差分法FDID方法介绍和示例, 附code和数据！10.双重差分法和事件研究法的区别主要在哪里？11.前沿, 合成双重差分法SDID方法介绍和示例, 附code和数据！12.具有空间溢出效应的双重差分法估计最全综述, 理论和操作尽有！13.最新Sun和Abraham(2021)和TWFE估计多期或交错DID并绘图展示结果！详细解读code！

---

下面这些短链接文章属于合集，可以收藏起来阅读，不然以后都找不到了。

5年，计量经济圈近1500篇不重类计量文章，

可直接在公众号菜单栏搜索任何计量相关问题,

Econometrics Circle

数据系列：空间矩阵 | 工企数据 | PM2.5 | 市场化指数 | CO2数据 |  夜间灯光 | 官员方言  | 微观数据 | 内部数据计量系列：匹配方法 | 内生性 | 工具变量 | DID | 面板数据 | 常用TOOL | 中介调节 | 时间序列 | RDD断点 | 合成控制 | 200篇合辑 | 因果识别 | 社会网络 | 空间DID数据处理：Stata | R | Python | 缺失值 | CHIP/ CHNS/CHARLS/CFPS/CGSS等 |干货系列：能源环境 | 效率研究 | 空间计量 | 国际经贸 | 计量软件 | 商科研究 | 机器学习 | SSCI | CSSCI | SSCI查询 | 名家经验计量经济圈组织了一个计量社群，有如下特征：热情互助最多、前沿趋势最多、社科资料最多、社科数据最多、科研牛人最多、海外名校最多。因此，建议积极进取和有强烈研习激情的中青年学者到社群交流探讨，始终坚信优秀是通过感染优秀而互相成就彼此的。