DID新进展：异质性多期DID估计的新方法-csdid
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2023 Stata寒假班
作者：彭晴 (暨南大学)
邮箱：natasha_pong@163.com
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 背景简介
2. 前提假设
3. 多期 DID 估计量
4. Stata 实操
4.1 csdid 命令介绍
4.2 csdid 命令实操
5. R 语言实操
5.1 R 包的安装和导入
5.2 R 语言的实操展示
6. 参考资料
7. 相关推文
1. 背景简介
双重倍差法 (Difference-in-Differences，DID)，是目前实证分析中用于识别因果关系的流行方法之一。标准的 DID 模型将样本分为两组：实验组和对照组；将时间分为两个阶段：政策发生前和政策发生后。所有的实验组样本都在同一时间点受到政策冲击。
随着 DID 方法的拓展，许多实证研究将其拓展为多期 DID，即实验组并非在同一时点遭受政策冲击。但是，自 2019 年来，不少学者纷纷指出这种多期 DID 有可能会产生有偏估计 (Athey and Imbens，2022；Baker et al.，2022；Goodman-Bacon，2021)。
其主要原因在于，多期 DID 估计的本质是多个不同处理效应的加权平均，权重可能存在为负的情形。在权重为负的情下，不同处理效应加权平均后得到的平均处理效应，可能会与真实的平均处理效应方向相反。Baker et al. (2021) 通过数据模拟发现，多期 DID 估计出来有偏误的处理效应甚至会与真实处理效应的符号相反。
为此，Callaway and Sant'Anna (2021) 提出了一种用于识别异质性多期 DID 的新方法，该新方法适用于以下三种情形：
时间分为多期； 实验组受到政策冲击的时间并非同一； 实验组和对照组只有在控制了协变量之后才满足平行趋势假定。
下面本文将详细介绍 Callaway and Sant'Anna (2021) 提出的多期 DID 估计量，以及对应的 Stata 实操和 R 语言实操。
2. 前提假设
假设时间用 来表示，个体样本用 来表示。并且，
为虚拟变量， 表示个体 在 期为 “被处理” 的状态，反之为 “未被处理” 的状态； 表示样本 在时间 首次变成 “被处理” 的状态； 表示个体样本 从 1 到 期都处于 “未被处理” 的状态，这类个体样本被称为 “从未被处理” 的对照组 (Nevered-treated)； 表示在给定协变量 时，样本在时间 受到处置的概率：
Callaway and Sant'Anna (2021) 提出的多期 DID 估计量必须满足以下几个假设：
处理效应不可逆转：一旦样本受到政策冲击，那么其在下一期会继续保留被处理的状态，即 ，意味着 也必然等于1。 独立同分布假设。 从未被处理的对照组满足有条件平行趋势假设：
还未被处理的对照组满足有条件平行趋势假设：
其中， 和 为样本在时间 被处理和未被处理的潜在结果。
值得注意的是，假设 3 和 4 是两个不同的有条件平行趋势假设。假设 3 意味着，在控制协变量 之后，在时间 首次 “被处理” 的样本与 “从未被处理” 的样本 (Nevered-treated) 需满足平行趋势假设。假设 4 意味着，在控制协变量 之后，在时间 首次 “被处理” 的样本与在时间 “还未被处理” 的样本 (Not-yet-treated) 需满足平行趋势假设。
3. 多期 DID 估计量
假设一个双向固定效应模型，具体形式如下：
其中， 为被解释变量， 为时间固定效应， 为个体固定效应， 为 0-1 虚拟变量，等于 1 表示受到政策影响，0 则反之， 是个体特征的向量， 为随机扰动项。
首先，Callaway and Sant’Anna (2021) 按照首次被处理的时间将处理组分为不同组别 (用字母 标识)，估计每个组别 在时间 的处理效应 。采用逆概率加权计算的 非参点估计量为：
需要注意的是，作者提到可以采用三种方法进行估计：回归 (OR)、逆概率加权法 (IPW) 以及双重稳健法 (DR) 。此处只介绍了逆概率加权法的估计量，其他两种方法详见 Callaway and Sant’Anna (2021) 的论文。
其次，Callaway and Sant’Anna (2021) 提供了不同的加总方法，用于计算平均处理效应：
选择不同的权重 ，可以计算出不同类型的平均处理效应，具体包括如下四种：
简单加权平均处理效应 (Simple ATT)：将 等权重简单加权求和； 分组平均处理效应 (Group ATT)：按照首次被处理的时间 () 分组加权求和的平均处理效应； 动态平均处理效应 (Dynamic ATT)：按照距离首次被处理时间分组加权求和的平均处理效应，类似于 Event Study； 日历平均处理效应 (Calendar Time ATT)：按照正常年份分组加权求和的平均处理效应。
4. Stata 实操
4.1 csdid 命令介绍
csdid 命令安装：
ssc install csdid, replace
csdid 命令语法：
csdid depvar [indepvars] [if] [in] [weight], [ivar(varname)] time(varname) gvar(varname) [options]
其中，
depvar：指定回归的被解释变量；indepvars：指定回归的解释变量；ivar：指定面板回归中的个体标识，如国家 ID、企业 ID 等；time：指定面板回归中的时间标识；gvar：分组标识，按首次被处理的时间分组；notyet: 定义 “从未被处理” 的样本 (Nevered-treated) 和 “还未被处理” 的样本 (Not-yet-treated) 为对照组。当不添加notyet时 (默认情况)，只选择 “从未被处理” 的样本 (Nevered-treated) 作为对照组。
此外，添加 method(method) 选项，用于选择估计方法。可选择的估计方法包括：drimp 为基于逆概率加权最小二乘法得到的双重稳健 DID 估计量，为默认估计方法；dripw 为基于逆概率的普通最小二乘法得到的双重稳健 DID 估计量；reg 为普通最小二乘法；stdipw 为标准化的逆概率加权法；ipw 为逆概率加权法。
添加 agg(aggtype) 选项，用于选择计算平均处理效应的加权方法。可选择的加权方法包括：simple 对应上述的 Simple ATT；group 对应上述的 Group ATT；calendar 对应上述的 Calendar Time ATT；event 对应上述的 Dynamic ATT。
关于标准误，默认情形是计算稳健标准误。添加 wboot 选项，可用 WildBootstrap 方法计算标准误；添加 cluster 选项，可计算聚类标准误。
4.2 csdid 命令实操
本文选择作者提供的样本数据 mpdta.dta 来展示 csdid 命令的基本操作。
mpdta.dta 数据集包括 500 个 country，时间范围是 2004年—2007年，一共有 2500 个观测值。其包括如下 6 个变量：year 为时间标识，countryreal 为个体标识，lpop 为协变量，lemp 为被解释变量，first_treat 为首次被处理的时间。
当 first_treat 取值为 0 时表示，该样本为 “从未被处理” 的样本 (Nevered-treated)。treat 为是否为处理组的标识，取值为 1 时表示该样本为处理组。
. use "https://friosavila.github.io/playingwithstata/drdid/mpdta.dta", clear
. des
Contains data from https://friosavila.github.io/playingwithstata/drdid/mpdta.dta
Observations: 2,500 Written by R.
Variables: 6 17 May 2021 11:45
---------------------------------------------------------------------------------
Variable Storage Display Value
name type format label Variable label
---------------------------------------------------------------------------------
year int %9.0g year
countyreal long %9.0g countyreal
lpop double %9.0g lpop
lemp double %9.0g lemp
first_treat int %9.0g first.treat
treat byte %9.0g treat
---------------------------------------------------------------------------------
首先，我们采用 dripw 方法估计出每个组别 在时间 的处理效应 。
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw)
Difference-in-difference with Multiple Time Periods
Number of obs = 2,500
Outcome model : least squares
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
g2004 |
t_2003_2004 | -0.015 0.022 -0.66 0.511 -0.058 0.029
t_2003_2005 | -0.076 0.029 -2.67 0.008 -0.133 -0.020
t_2003_2006 | -0.140 0.035 -3.97 0.000 -0.210 -0.071
t_2003_2007 | -0.107 0.033 -3.25 0.001 -0.171 -0.042
-------------+----------------------------------------------------------------
g2006 |
t_2003_2004 | -0.000 0.022 -0.02 0.983 -0.044 0.043
t_2004_2005 | -0.006 0.018 -0.34 0.737 -0.042 0.030
t_2005_2006 | 0.001 0.019 0.05 0.961 -0.037 0.039
t_2005_2007 | -0.041 0.020 -2.09 0.036 -0.080 -0.003
-------------+----------------------------------------------------------------
g2007 |
t_2003_2004 | 0.027 0.014 1.90 0.057 -0.001 0.054
t_2004_2005 | -0.005 0.016 -0.29 0.771 -0.035 0.026
t_2005_2006 | -0.028 0.018 -1.56 0.118 -0.064 0.007
t_2006_2007 | -0.029 0.016 -1.77 0.076 -0.061 0.003
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details
从估计结果可知，组别 按照首次被处理的时间进行划分，分为三组：在 2004 年首次受到冲击 (标识 )；在 2006 年首次受到冲击 (标识 )；在 2007 年首次受到冲击 (标识 )。以 组的结果为例，在 2004 年之后，冲击对该组样本存在显著的负向影响。
进一步地，我们还可以选择不同的加总方式，计算平均处理效应。此处，我们以计算动态平均处理效应为例进行展示。在上述命令的基础上，添加选项 agg(event) 可计算动态平均处理效应。
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(event)
Difference-in-difference with Multiple Time Periods
Number of obs = 2,500
Outcome model : least squares
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
T-3 | 0.027 0.014 1.90 0.057 -0.001 0.054
T-2 | -0.004 0.013 -0.28 0.780 -0.029 0.022
T-1 | -0.023 0.014 -1.60 0.109 -0.052 0.005
T+0 | -0.021 0.011 -1.83 0.067 -0.044 0.001
T+1 | -0.053 0.016 -3.24 0.001 -0.085 -0.021
T+2 | -0.140 0.035 -3.97 0.000 -0.210 -0.071
T+3 | -0.107 0.033 -3.25 0.001 -0.171 -0.042
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details
5. R 语言实操
5.1 R 包的安装和导入
首先，安装作者提供的 did-R-Package。
install.packages("did")
安装好之后，可调用 did 包。
library("did")
5.2 R 语言的实操展示
调用 did 包之后，导入作者提供的样本数据 mpdta.rda。
load("mpdta.rda")
接着，att_gt 命令可计算每个组别 在时间 的处理效应 。其中，参数 yname 为被解释变量；tname 为时间标识；idname 为个体标识；gname 分组标识，按首次被处理的时间分组；xformla 为协变量；data 为所使用的数据集。
除此之外，我们还可以增加一些 option。例如，可增加 control_group = "notyettreated" 选项，表示定义 “从未被处理” 的样本 (Nevered-treated) 和 “还未被处理” 的样本 (Not-yet-treated) 为对照组。与 Stata 一致，默认情形为只选择 “从未被处理” 的样本 (Nevered-treated) 为对照组。
下面为使用 R 包计算出来的结果，与 Stata 计算出来的结果一致。
out <- att_gt(yname="lemp",
tname="year",
idname="countyreal",
gname="first.treat",
xformla=~lpop,
data=mpdta)
Group-Time Average Treatment Effects:
Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band]
2004 2004 -0.0145 0.0229 -0.0760 0.0470
2004 2005 -0.0764 0.0308 -0.1590 0.0061
2004 2006 -0.1404 0.0396 -0.2467 -0.0342 *
2004 2007 -0.1069 0.0331 -0.1958 -0.0180 *
2006 2004 -0.0005 0.0223 -0.0603 0.0594
2006 2005 -0.0062 0.0194 -0.0582 0.0458
2006 2006 0.0010 0.0204 -0.0538 0.0557
2006 2007 -0.0413 0.0198 -0.0945 0.0119
2007 2004 0.0267 0.0135 -0.0094 0.0628
2007 2005 -0.0046 0.0169 -0.0498 0.0407
2007 2006 -0.0284 0.0181 -0.0769 0.0200
2007 2007 -0.0288 0.0150 -0.0689 0.0113
---
Signif. codes: `*' confidence band does not cover 0
P-value for pre-test of parallel trends assumption: 0.23267
Control Group: Never Treated, Anticipation Periods: 0
Estimation Method: Doubly Robust
最后，我们使用 aggte 命令计算动态平均处理效应，其结果也与 Stata 计算的一致。其中，out 为上面 att_gt 命令运行的结果，type 为选择的加总方式。
aggte(out, type="dynamic")
Overall summary of ATT's based on event-study/dynamic aggregation:
ATT Std. Error [ 95% Conf. Int.]
-0.0804 0.0206 -0.1207 -0.04 *
Dynamic Effects:
Event time Estimate Std. Error [95% Simult. Conf. Band]
-3 0.0267 0.0146 -0.0109 0.0644
-2 -0.0036 0.0129 -0.0369 0.0297
-1 -0.0232 0.0146 -0.0611 0.0146
0 -0.0211 0.0113 -0.0502 0.0081
1 -0.0530 0.0170 -0.0969 -0.0091 *
2 -0.1404 0.0401 -0.2440 -0.0369 *
3 -0.1069 0.0345 -0.1960 -0.0178 *
---
Signif. codes: `*' confidence band does not cover 0
Control Group: Never Treated, Anticipation Periods: 0
Estimation Method: Doubly Robust
需要注意的是，在 R 语言中，可选择的加总方式也分为四种 (与上述介绍的一致)。但与 Stata 命令不同的是，计算动态平均处理效应的 option 在 R 中为 dynamic (在 Stata 中为 event)。其余三种加总方式，在 Stata 和 R 语言中的表示方式一致，分别为 simple、group 和 calendar。
6. 参考资料
Callaway B, Sant’Anna P H C. Difference-in-differences with multiple time periods[J]. Journal of Econometrics, 2021, 225(2): 200-230. -PDF- Sant’Anna P H C, Zhao J. Doubly robust difference-in-differences estimators[J]. Journal of Econometrics, 2020, 219(1): 101-122. -PDF- Athey S, Imbens G W. Design-based analysis in difference-in-differences settings with staggered adoption[J]. Journal of Econometrics, 2022, 226(1): 62-79. -PDF- Baker A C, Larcker D F, Wang C C Y. How much should we trust staggered difference-in-differences estimates?[J]. Journal of Financial Economics, 2022, 144(2): 370-395. -PDF- Goodman-Bacon A. Difference-in-differences with variation in treatment timing[J]. Journal of Econometrics, 2021, 225(2): 254-277. -PDF-
7. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh did, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：最新课程 ⏫NEW-面板数据因果推断：从入门到精通 专题：倍分法DID Stata：异质性稳健DID估计量方法汇总 DID偏误问题：多时期DID的双重稳健估计量(下)-csdid DID偏误问题：两时期DID的双重稳健估计量(上)-drdid DID的陷阱和注意事项 Stata：事件研究法的稳健有效估计量-did_imputation DID最新进展：异质性处理条件下的双向固定效应DID估计量 (TWFEDD) 专题：公开课 ⏩公开课：双重差分方法的新进展——交错型DID
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。