Stata：各种DID估计量的比较分析
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会课程 · 2023 暑期班
作者：周嘉怡 (中山大学)
邮箱：sysu_zhoujy@126.com
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 简介
2. 案例介绍
3. 对双向固定效应偏误的回顾
4. 不同稳健估计量分析
4.1 did_imputation
4.2 did_multiplegt
4.3 csdid
4.4 eventstudyinteract
4.5 reghdfe
4.6 bacondecomp
4.7 did2s
4.8 stackedev
5. 对比与总结
5.1 结果图片集合
5.2 与事实对比结果
6. 相关推文
1. 简介
双重差分模型 (Difference-in-Differences) 是用于政策评估 (Policy Evaluation) 的常用模型。最早的双重差分模型设计来源于 1855 年斯诺建立的一项自然实验，即验证霍乱是由于水而非空气传播的事实。在这一自然实验中，我们了解到了“反事实”的选取以及“平行趋势假设”对 DID 模型设计的重要性。
基于 DID 模型在现实生活中的广泛应用，不同学者基于不同的假设条件建立了多种稳健估计量。本文将基于给定的两个案例背景，介绍 Albert Alex Zevelev (2022) 在 GitHub 仓库 Compare-DiD-Estimators 中给出的各种 DID 估计方法的适用场景和性质。
该仓库中提供了两个案例，不同之处在于案例 1 中设定了 homogenous (constant) TEs，而案例 2 中设定了 heterogenous (constant) TEs。因此本文会将两个案例代码共同介绍，结果部分单独分析。
2. 案例介绍
本文使用的基础数据结构为 100 个体 × 15 时期 = 1500 个观察值的平衡面板数据。后文中的估计量均使用这份模拟数据。
clear all
timer clear
* 设定 1500 个观察值及其他暂元
set seed 10
global T = 15
global I = 100
global pre 5
global post 8
global ep event_plot
global g0 "default_look"
global g1 xla(-$pre (1) $post) // global g1 xla(-5(1)5)
global g2 xt("Periods since the event")
global g3 yt("Average causal effect")
global g $g1 $g2 $g3
global t "together"
set obs `=$I*$T'
* 生成 id 与时间
gen i = int((_n-1)/$T )+1 // unit id
gen t = mod((_n-1),$T )+1 // calendar period
tsset i t
* 设定 35 个单位处理时间为 7, 35 个单位处理时间为 11
gen Ei = 7 if t==1 & i>=1 & i <=35
replace Ei = 11 if t==1 & i>=36 & i <=70
bys i (t): replace Ei = Ei[1]
* 生成处理变量, K 为相对处理时间, D 为处理时间哑变量
gen K = t-Ei
gen D = K>=0 & Ei!=.
两个案例设定了不同的处理效应，其代码分别如下：
案例1
* 各组别接受处理后的效应相同
gen tau = cond(D==1, 49, 0)
* 生成误差项及结果变量 Y
gen eps = rnormal()
gen Y = i/100 + tau*D + 3*eps
* 对未处理组 Ei 部分设置新变量进行填充
gen gvar = cond(Ei==., 0, Ei)
gen never_tr = Ei==.
* 建立滞后期和领先期的时间序列变量
forval x = 1/$pre {
gen F_`x' = K == -`x'
replace F_`x' = 0 if never_tr==1
}
forval x = 0/$post {
gen L_`x' = K == `x'
replace L_`x' = 0 if never_tr==1
}
rename F_1 ref
* 将每个面板叠加到同一图形上
xtline Y, overlay legend(off) name(gY, replace)
案例2
* 不同组别接受处理后效应不同
gen tau = 0
replace tau = 20 if D==1 & i>=1 & i <=35
replace tau = 60 if D==1 & i>=36 & i <=70
* 生成误差项及结果变量 Y
gen eps = rnormal()
gen Y = i/100 + tau*D + 3*eps
* 对未处理组 Ei 部分设置新变量进行填充
gen gvar = cond(Ei==., 0, Ei)
gen never_tr = Ei==.
* 建立滞后期和领先期的时间序列变量
forval x = 1/$pre {
gen F_`x' = K == -`x'
replace F_`x' = 0 if never_tr==1
}
forval x = 0/$post {
gen L_`x' = K == `x'
replace L_`x' = 0 if never_tr==1
}
rename F_1 ref
* 将每个面板叠加到同一图形上
xtline Y, overlay legend(off) name(gY, replace)
以上两个图片可以使我们更直观地体会到两个案例所设定的不同背景：案例 1 中的两个处理组差异在于处理时间不同；而案例 2 中的两个处理组不仅设置了处理时间上的差异，处理效果也明显不同。
3. 对双向固定效应偏误的回顾
Causal Inference: The Mixtape 一书第九章节提到了 Goodman-Bacon (2019) 对双向固定效应估计的有效分解，即培根分解。该定理指出双向固定效应估计量是所有潜在因素的加权平均值， DD 估计权重均基于组大小和处理方差。在方差加权共同趋势 (VWCT) 和时间不变处理效果的假设下，方差加权 ATT 是所有可能的 ATT 的加权平均值。
在更严格的假设下，该估计与 ATT 完全匹配。但当存在时变处理效果时，情况并非如此，因为双向固定效应估计的微分时序设计中的时变处理效果会产生偏差。因此，双向固定效应模型可能存在严重偏差。
假设在此设计中有三组：早期处理组 ()，后来处理的一组 ()，以及一个从未被处理过的组别 ()，组 和 相似之处在于它们都得到了处理，但它们的不同之处在于 处理早于 。
以下三个公式分别代表上述三个组别相互比对所得出的 ：
考虑到这种表示法，DD 参数估计可以分解如下：
而该式中的权重系数表示如下：
根据以上分析，可以得出：延长或缩短面板实际上可以纯粹通过改变组处理方差来改变点估计。那么，多种修正后的估计方法所解决的究竟是实验组处理时间不同所带来的差异，还是由于二者处理效用不同带来的影响，亦或是二者兼而有之。本文将对这一问题进行解答。
4. 不同稳健估计量分析
本文所涉及的大部分估计方法所解决的主要问题是 TWFE 在交叠 DID 估计中的偏误分解，可以分为以下三种类型：
(一) 组别-时期平均处理效应
DeChaisemartin 和 d'Haultfœuille (2020) 提出的估计量 ( did_multiplegt)Sun 和 Abraham (2021) 提出的估计量 ( eventstudyinteract)Callaway 和 Sant’Anna (2021) 提出的估计量 ( csdid)
(二) 插补估计量
Borusyak et al.(2021) 提出的估计量 ( did_imputation)Gardner(2021) 提出的估计量 ( did2s)
(三) 堆叠回归估计量
Cengiz et al. (2019) 提出的估计量 ( stackedev)
注：部分原仓库代码在 stata17.0 版本执行下会报错，以下代码已针对相应问题进行更改处理，因此与原仓库中代码有一定差异，下文中已加入注释标识。
4.1 did_imputation
did_imputation Y i t Ei, horizons(0/$post) pretrend($pre) minn(0)
estimates store bjs // 存储估计量
$ep bjs, $t $g0 graph_opt($g ti("BJS 21") name(gBJS, replace)) // 绘图
其中，
horizons指定每个视界的处理效果的加权平均值/总和；pretrend指定对平行趋势进行预测；minn指定经过处理的观测值的最小有效数，设置为minn(0)指报告所有系数。
Borusyak et al. (2021) 提供了一种基于插补的反事实方法解决 TWFE 的估计偏误问题。通过改良原有假设使其最终估计量满足有效性，渐近正态性，并检验了平行趋势。
4.2 did_multiplegt
did_multiplegt Y i t D, robust_dynamic dynamic($post) placebo($pre) ///
breps(20) cluster(i) graphoptions(name(gCD, replace))
// 计算并绘图 (注：仓库中原 code 执行会出错，这里加入 graphoptions 使其符合 did_multiplegt 的语法)
matrix dcdh_b = e(estimates) // 存储估计量
matrix dcdh_v = e(variances)
其中，
robust_dynamic指定估计量为 Chaisemartin and D'Haultfoeuille (2020b) 版本，若不指定则为 2020a 版本；dynamic指定了要估计的动态处理效果的数量。该选项只能在指定了robust_dynamic选项时使用；placebo指定安慰剂效应的数量；breps指定输入重抽样的次数，用以计算标准误；cluster指定引导生成相应变量标准误；graphoptions绘图选项。
De Chaisemartin 和 D‘Haultfoeuille (2020) 提出通过加权计算两种处理效应的值得到平均处理效应的无偏估计，多期多个体 DIDM 模型使得放松了对样本处理的限制，满足了大部分情况下政策力度不同的现实要求，也放开了传统 DID 模型中一次性处理的条件，与现实的贴合度更高。
4.3 csdid
csdid Y, ivar(i) time(t) gvar(gvar) notyet
estat event, estore(cs) // 生成并存储估计量
$ep cs, stub_lag(Tp#) stub_lead(Tm#) $t $g0 graph_opt($g ti("CS 20") ///
name(gCS, replace)) // 绘图
其中，
ivar指定面板回归中的个体标识；time指定面板回归中的时间标识；gvar分组标识，按首次被处理的时间分组；notyet定义“从未被处理”的样本 (Nevered-treated) 和“还未被处理”的样本 (Not-yet-treated) 为对照组；event指定估计动态 ATT。
4.4 eventstudyinteract
eventstudyinteract Y L_* F_*, vce(cluster i) absorb(i t) ///
cohort(Ei) control_cohort(never_tr)
$ep e(b_iw)#e(V_iw), stub_lag(L_#) stub_lead(F_#) $t $g0 ///
graph_opt($g ti("SA 20") name(gSA, replace)) // 绘图
matrix sa_b = e(b_iw) // 存储估计量
matrix sa_v = e(V_iw)
其中，
vce指定方差-协方差分量估计 (variance-covariance component estimation，vce)；absorb指定个体和时间固定效应。
4.5 reghdfe
reghdfe Y L_* F_*, a(i t) cluster(i)
estimates store ols // 存储估计量
$ep ols, stub_lag(L_#) stub_lead(F_#) $t $g0 graph_opt($g ti("OLS") ///
name(gOLS, replace)) // 绘图
4.6 bacondecomp
bacondecomp Y D, ddetail gropt(legend(off) name(gGB, replace))
// 注：仓库中原 code 执行会出错，这里加入 gropt 使其符合 bacondecomp 的语法
4.7 did2s
did2s Y, first_stage(i t) second_stage(F_* L_*) treatment(D) cluster(i)
$ep, stub_lag(L_#) stub_lead(F_#) $t $g0 graph_opt($g ti("Gardner 21") name(gG, replace)) // 绘图
matrix did2s_b = e(b) // 存储估计量
matrix did2s_v = e(V)
其中，first_stage 指定第一阶段公式，包括用于估计 的固定效应和协变量，不能放入处理变量。两阶段双重差分的原理是先识别组别效应和时期效应，然后在第二阶段将其剔除后，再对处理变量进行回归。主要用于解决的问题是处理组个体接受处理的时间是交错的，而且平均处理效应随着组别以及时间发生变化的情况。
4.8 stackedev
stackedev Y F_* L_* ref, cohort(Ei) time(t) never_treat(never_tr) ///
unit_fe(i) clust_unit(i)
$ep, stub_lag(L_#) stub_lead(F_#) $t $g0 graph_opt($g ti("CDLZ 19") ///
name(gCDLZ, replace)) // 绘图
matrix stackedev_b = e(b) // 存储估计量
matrix stackedev_v = e(V)
Cengiz et al. (2019) 所设计的基本思路是将数据集重建为相对事件时间的平衡面板，然后控制组群效应和时间固定效应，以得到处理效应的加权平均值。
5. 对比与总结
5.1 结果图片集合
* 绘图集合
graph combine gY gOLS gGB gBJS gCD gCS gSA gG gCDLZ, ycommon name(combined, replace)
案例1
案例2
5.2 与事实对比结果
* 构建真实结果向量
matrix btrue = J(1,9,.)
matrix colnames btrue = tau0 tau1 tau2 tau3 tau4 tau5 tau6 tau7 tau8
qui forvalues h = 0/8 {
sum tau if K==`h'
matrix btrue[1,`h'+1]=r(mean)
}
* 利用存储估计量与真实结果做比对
event_plot ///
btrue# bjs dcdh_b#dcdh_v cs sa_b#sa_v did2s_b#did2s_v stackedev_b#stackedev_v ols, ///
stub_lag( tau# tau# Effect_# Tp# L_# L_# L_# L_#) ///
stub_lead(pre# pre# Placebo_# Tm# F_# F_# F_# F_#) ///
plottype(scatter) ciplottype(rcap) ///
together perturb(-0.325(0.1)0.325) trimlead(5) noautolegend ///
graph_opt( ///
title("Event study estimators in a simulated panel", size(med)) ///
xtitle("Periods since the event", size(small)) ///
ytitle("Average causal effect", size(small)) xlabel(-$pre(1)$post) ///
legend(order(1 "Truth" 2 "BJS" 4 "dCdH" ///
6 "CS" 8 "SA" 10 "G" 12 "CDLZ" 14 "TWFE") rows(2) position(6) region(style(none))) ///
/// the following lines replace default_look with something more elaborate
xline(-0.5, lcolor(gs8) lpattern(dash)) yline(0, lcolor(gs8)) graphregion(color(white)) ///
bgcolor(white) ylabel(, angle(horizontal)) ) ///
lag_opt1(msymbol(+) color(black)) lag_ci_opt1(color(black)) ///
lag_opt2(msymbol(O) color(cranberry)) lag_ci_opt2(color(cranberry)) ///
lag_opt3(msymbol(Dh) color(navy)) lag_ci_opt3(color(navy)) ///
lag_opt4(msymbol(Th) color(forest_green)) lag_ci_opt4(color(forest_green)) ///
lag_opt5(msymbol(Sh) color(dkorange)) lag_ci_opt5(color(dkorange)) ///
lag_opt6(msymbol(Th) color(blue)) lag_ci_opt6(color(blue)) ///
lag_opt7(msymbol(Dh) color(red)) lag_ci_opt7(color(red)) ///
lag_opt8(msymbol(Oh) color(purple)) lag_ci_opt8(color(purple))
在案例 1 中所设定的情况下，仓库中所提供的估计量方法最终得出的结果是无偏的。
在案例 2 中所设定的情况下，TWFE-OLS 和 stackedev 两种方法计算的估计量为有偏的，而其余方法得出的估计量是无偏的。
文中提出的改进后的 DID 估计量均通过不同方法改进了 TWFE 估计偏误的问题，同时也说明了基于案例2提供的条件，stackedev 所使用的堆叠法依然不能很好地解决有偏的问题。由于本文设定案例设置组别数较少，不能很好地体现各种估计方法在其他应用场景中的差异，还需基于各种估计方法下的基本假设灵活应用。
6. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh did, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：专题课程 DID新进展：异质性多期DID估计的新方法-csdid ⏩ 因果推断专题-RDD-DID-IV-合成控制 专题：倍分法DID Stata：异质性稳健DID估计量方法汇总 DID偏误问题：多时期DID的双重稳健估计量(下)-csdid DID偏误问题：两时期DID的双重稳健估计量(上)-drdid DID的陷阱和注意事项 Stata：事件研究法的稳健有效估计量-did_imputation DID最新进展：异质性处理条件下的双向固定效应DID估计量 (TWFEDD) Stata倍分法新趋势：did2s-两阶段双重差分模型 Stata-DID：不同处理时点不同持久期的倍分法（flexpaneldid） DIDM：多期多个体倍分法-did_multiplegt 面板PSM+DID如何做匹配？ 专题：公开课 ⏩公开课：双重差分方法的新进展——交错型DID
课程推荐：2023 暑期班
主讲老师：连玉君，王群勇
🍓 课程主页：https://www.lianxh.cn/news/fdc69c3695aec.html
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。