参考文献：
1. Cunningham，Causal Inference:The Mixtape, Yale Press, 2021
2. Baker, Andrew. Difference-in-Differences Methodology,2019
3. Coady Wing, Seth M. Freedman, and Alex Hollingsworth. Stacked Difference-in-Differences,NBER Working Paper No. 32054, 2024
一、堆叠DID概述
交叠处理下，异质性处理效应使得TWFE有偏，很多新的估计量已经出现了（参考最新的推文）。
堆叠DID（stacked did）是应对交叠处理偏误的一种方法。近些年，在应用研究中广泛采用，例如，Cengiz et al.,2019; Deshpande and Li, 2019; Butters et al., 2022; Callison and Kaestner, 2014。
例如Cengiz et al.（2019）估计了最低工资变化对低收入劳动的影响，基于美国1979-2016年的138次州水平的最低工资变化，采用DID。在附录中，Cengiz et al.（2019）指出，加总OLS估计的DID估计量会存在一些问题。作为稳健性检验，作者分别估计了每个事件的单个处理效应，并做出分布图。
作者创建了138个最低工资变化事件的处理地区数据集（将总的数据集分成不同的子集），其中，包括对应事件地区和所有“干净的”控制地区（在估计事件时间窗口还没有发生最低工资变化的州）的结果变量和控制变量。对于每个事件，作者都跑一个单一处理个体的DID：
然后，Cengiz et al.（2019）画出所有处理效应的分布，其中，置信区间使用Ferman and Pinto（2022）的异方差稳健聚类余值bootstrap方法（因为在每个DID回归中，只有单一处理个体）。Cengiz et al.（2019）也将每个事件数据集堆叠在一起来计算平均处理效应。正如作者指出，堆叠DID用更加严格的标准来筛选控制组，且对交叠处理的异质性处理效应偏误问题更加稳健。依靠堆叠事件，该方法等价于事件同时发生的情形，因此，防止了交叠处理的负权重问题（Baker，2019）。
也就是说，堆叠DID对每个有效的处理事件（子样本，排除“后处理 vs 先处理”的比较）都构造一个子样本数据集，然后对子样本数据集跑DID回归。
Joshua Bleiberg基于Cengiz et al.（2019）的文章写了一个stata命令stackedev，我们用 Asjad Naqvi 的模拟数据来看看这个命令的用法：
* 模拟数据
clear
local units = 30
local start = 1
local end = 60
local time = `end' - `start' + 1
local obsv = `units' * `time'
set obs `obsv'
egen id = seq(), b(`time')
egen t = seq(), f(`start') t(`end')
sort id t
xtset id t
set seed 20211222
gen Y = 0 // outcome variable
gen D = 0 // intervention variable
gen cohort = . // treatment cohort
gen effect = . // treatment effect size
gen first_treat = . // when the treatment happens for each cohort
gen rel_time = . // time - first_treat
levelsof id, local(lvls)
foreach x of local lvls {
local chrt = runiformint(0,5)
replace cohort = `chrt' if id==`x'
}
levelsof cohort , local(lvls)
foreach x of local lvls {
local eff = runiformint(2,10)
replace effect = `eff' if cohort==`x'
local timing = runiformint(`start',`end' + 20) //
replace first_treat = `timing' if cohort==`x'
replace first_treat = . if first_treat > `end'
replace D = 1 if cohort==`x' & t>= `timing'
}
replace rel_time = t - first_treat
replace Y = id + t + cond(D==1, effect * rel_time, 0) + rnormal()
// generate leads and lags (used in some commands)
summ rel_time
local relmin = abs(r(min))
local relmax = abs(r(max))
// leads
cap drop F_*
forval x = 2/`relmin' { // drop the first lead
gen F_`x' = rel_time == -`x'
}
//lags
cap drop L_*
forval x = 0/`relmax' {
gen L_`x' = rel_time == `x'
}
// generate the control_cohort variables (used in some commands)
gen never_treat = first_treat==.
sum first_treat
gen last_cohort = first_treat==r(max) // dummy for the latest- or never-treated cohort
// generate the gvar variabls (used in some commands)
gen gvar = first_treat
recode gvar (. = 0)
xtline Y, overlay legend(off)
stata命令的使用：
* 安装命令
ssc install stackedev, replace
help stackedev
*运行堆叠DID
stackedev Y F_* L_*, cohort(first_treat) time(t) never_treat(never_treat) unit_fe(id) clust_unit(id)
event_plot, default_look graph_opt(xtitle("相对事件时间") ytitle("平均处理效应") xlabel(-10(1)10) ///
title("stackedev")) stub_lag(L_#) stub_lead(F_#) trimlag(10) trimlead(10) together
Cengiz et al. (2019)和Deshpande and Li (2019)用“干净的”控制组来构建堆叠数据集，但是会存在事件研究估计系数的样本不平衡的问题，这可能会产生样本结构变化，进而不可比，使得事件研究系数没有因果含义。Butters et al. (2022)用了干净的控制组，也施加了结构平衡的限制。Callison and Kaestner (2014)用两个事件时间的堆叠数据集，并用基准期的平均结果来匹配处理组和控制组。
此外，Cengiz et al. (2019)报告的标准误实在group×sub-event层面聚类，Callison and Kaestner (2014)、Deshpande and Li (2019)和Butters et al. (2022)则聚类在group层面。
这些都只是应用研究文献，并没有给出堆叠DID估计量的识别公式。Wing et al.（2024）正式给出了堆叠DID估计量的正式推导，并提出了一种加权方法来计算平均处理效应。而且应用型文献的回归方程仅仅包括group×sub-event、time×sub-event固定效应，并没有使用饱和式堆叠回归方程。
Wing et al.（2024）提出了一类构建堆叠数据集的包容性条件来删节事件研究窗口，以确保每个子事件（sub-event）数据集处理前后的时期数是平衡的。然后，作者提出了一种新方法来估计总的平均处理效应（ATT），这种单一堆叠回归法允许使用传统的方法来进行统计推断。
Wing et al.（2024）认为，他们的删节总ATT（trimmed aggregate ATT）有三个优势：（1）所有的估计量都有一致的因果理解，因为它是一系列平衡因果效应的凸结合；（2）由于平均处理效应来自于删节的处理效应，因此，事件时间上的系数裱花反映的就是处理效应动态，而不是每个事件时间上的结构变化；（3）在DD假设下，处理前的系数应该等于0。因为堆叠数据的结构在事件时间上是非常稳定的，处理前处理效应系数的值就是差异化的处理前趋势，这可以理解成处理前共同趋势或者预期。
二、 stata例子
我们用Wing et al.（2024）提供的例子来看看他们的删节堆叠DID估计量：美国ACA医疗补助计划扩围对19-60岁的成年人未保险率的影响。数据集跨度2008-2021年，美国51个州的面板数据，共714个样本。
***********************************************************
*Getting Started With Stacked DID
*A side-by-side comparison of R and Stata code
*AUTHOR:Coady Wing, Alex Hollingsworth, and Seth Freedman
***********************************************************
// Load the data
clear
cd "/Users/xuwenli/Library/CloudStorage/OneDrive-个人/DSGE建模及软件编程/教学大纲与讲稿/应用计量经济学讲稿/应用计量经济学讲稿与code/DID与SC/"
import delimited "acs1860_unins_2008_2021.csv",clear
sum
tab adopt_year
其中，unins是结果变量——19-60岁的成年人没有参加健康保险的比例；adopt_year是初次实施ACA计划的年份；还未实施ACA计划的州state用NA表示。如下图所示，ACA扩围计划分别在2014、2015、2016、2019、2020、2021年实施。
Wing et al.（2024）指出，堆叠DID最关键的是要重组数据集，从原始面板数据中重新形成“干净”控制组的子数据集。因此，作者提供了一个stata命令create_sub_exp来重塑堆叠DID所需的子数据集。如果对这个命令比较感兴趣，请去网站下载（https://rawcdn.githack.com/hollina/stacked-did-weights/18a5e1155506cbd754b78f9cef549ac96aef888b/stacked-example-r-and-stata.html）。需要注意的是，作者并没有以.ado的形式呈现上述命令。如果有需要，可以给我发邮件。我跟Wing邮件交流的时候，发现的他们的手册写得并不是很清楚，初次使用可能会报错。我自己的步骤是：
① 首先，将Wing et al.（2024）的“create_sub_exp()”函数命令copy到stata的dofile中，然后保存成create_sub_exp.ado；
② 将create_sub_exp.ado放置在stata/ado/personal文件夹中
命令格式：
create_sub_exp,
timeID() groupID() adoptionTime() focalAdoptionTime() kappa_pre() kappa_post()
- timeID:时间变量
- groupID：个体变量
- adoptionTime：每个个体的处理时点，对于从未处理的个体是NA；
- focalAdoptionTime：用户声明的子事件的处理时点；
- kappa_pre：用户声明的处理前时期；
- kappa_post：用户声明的处理后时期；
我们知道，堆叠DID就是按照每个事件时间重新构建数据子集。下面，来看看这个命令重塑的ACA医疗扩围计划的数据子集。例如，我们先看看2014年实施的ACA计划的子集。
/* Create sub-experiment data for stack */
//Making sub-experimental data sets
* Save dataset
preserve
* Run this function with focal year 2014
create_sub_exp, ///
timeID(year) ///
groupID( statefip) ///
adoptionTime(adopt_year) ///
focalAdoptionTime(2014) ///
kappa_pre(3) ///
kappa_post(2)
* Open temp dataset created with function
use temp/subexp2014.dta, clear
* Summarize
sum statefip year adopt_year unins treat post event_time feasible sub_exp
* Restore dataset
restore
从上述结果，我们可以看出，2014年处理时点的子样本集有276个样本。这个时候，我们就可以跑DID回归了:
* twfe
xtset statefip year
gen d= treat*post
reghdfe unins d,ab(statefip year) cluster(statefip)
从上述结果可以看出，2014年处理时点的数据子集的DD估计系数为-0.019，在95%的置信水平上显著。意味着ACA扩围计划降低了未保险比例，统计上显著。
实践中，我们可能对单一数据子集并不感兴趣。我们想要所有的单一处理时点的数据子集堆叠在一起。下面，用一个stata循环来运行上述命令得到每一个实施年份的数据子集。然后，将单一数据子集堆叠在一起，形成一个“垂直”结合的数据集stacked_dtc。注意，在这个数据集中并不包含2020和2021两年，因为Wing et al.（2024）在删节的时候，选取了处理前3期（-3，-2，-1），处理后3期（0,1,2），而2020和2021处理后并没有3期数据可用，因此没有构造它们的数据子集。这就是Wing et al.（2024）提出的“包容性标准”——让事件研究设计不存在结构性变化，即是结构性平衡的。
// Build the stack of sub-experiments
//create the sub-experimental data sets
levelsof adopt_year, local(alist)
di "`alist'"
qui{
// Loop over the events and make a data set for each one
foreach j of numlist `alist' {
// Preserve dataset
preserve
// run function
create_sub_exp, ///
timeID(year) ///
groupID( statefip) ///
adoptionTime(adopt_year) ///
focalAdoptionTime(`j') ///
kappa_pre(3) ///
kappa_post(2)
// restore dataset
restore
}
// Append the stacks together, but only from feasible stacks
* Determine earliest and latest time in the data.
* Used for feasibility check later
sum year
local minTime = r(min)
local maxTime = r(max)
local kappa_pre = 3
local kappa_post= 2
gen feasible_year = adopt_year
replace feasible_year = . if adopt_year < `minTime' + `kappa_pre'
replace feasible_year = . if adopt_year > `maxTime' - `kappa_post'
sum feasible_year
local minadopt = r(min)
levelsof feasible_year, local(alist)
clear
foreach j of numlist `alist' {
display `j'
if `j' == `minadopt' use temp/subexp`j', clear
else append using temp/subexp`j'
}
// Clean up
* erase temp/subexp`j'.dta
}
* Summarize
sum statefip year adopt_year unins treat post event_time feasible sub_exp
* Treated, control, and total count by stack
preserve
keep if event_time == 0
gen N_treated = treat
gen N_control = 1 - treat
gen N_total = 1
collapse (sum) N_treated N_control N_total, by(sub_exp)
list sub_exp N_treated N_control N_total in 1/4
/*
sumup treat if event_time == 0, s(N)
stacked_dtc[event_time==0,
.(N_treated = fsum(treat),
N_control = fsum(1-treat),
N_total = .N
),
by = sub_exp][order(sub_exp)]
*/
restore
上表展示了堆叠数据中包含四个处理时点：2014、2015、2016、2019。且2014年有28个处理组，2015年有3个，2016年有2个，2019年也有2个。与此对比，2014、2015、2016年包含18个干净的控制组，2019年包含11个。
其实，我们此时可以分别利用上述4个数据子集来跑DID回归，得到TWFE DID估计量，这些估计量都是对应处理年份ACA的平均处理效应，大家可以自己去试试。
但是，我们可能更关心所有的处理时点的总的平均处理效应，也就是上述4个处理事件的ACA总的平均处理效应。这个时候，可以使用Wing et al.（2024）给出的另一个stata命令compute_weights。这个命令可以构造恰当的权重来加权堆叠DID回归，也就是说，这个命令利用上述构建的堆叠数据集来计算处理组和控制组的样本权重。命令格式：
compute_weights,
treatedVar(string) ///
eventTimeVar(string) ///
groupID(string) ///
subexpVar(string)
- treatedVar：处理变量，在给定数据子集中，个体是否作为处理组；
- eventTimeVar：处理时点
- groupID：个体变量
- subexpVar：每个事件时间变量
// Use compute_weights() to compute the stacked weights
compute_weights, ///
treatedVar(treat) ///
eventTimeVar(event_time) ///
groupID(statefip) ///
subexpVar(sub_exp)
* Summarize
by sub_exp:sum stack_weight if treat == 0 & event_time == 0
得到权重后，我们可以利用堆叠数据集来跑堆叠回归，得到总的ACA处理效应：
// Estimate the stacked regression
// Create dummy variables for event-time
char event_time[omit] -1
xi i.event_time
// Run regression
qui reghdfe unins i.treat##i._I* [aw = stack_weight], cluster(statefip) absorb(treat event_time)
est sto weight_stack
// Show results
esttab weight_stack, keep(1.treat#1*) se
// event study plot
coefplot,yline(0) vertical drop(_cons)
