Stata：各种DID估计量的比较分析

发布时间：2023-03-24 阅读 1506

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

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

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

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

☝ PDF下载 - 推文合集

作者：周嘉怡 (中山大学)
邮箱：sysu_zhoujy@126.com

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) 对双向固定效应估计的有效分解，即培根分解。该定理指出双向固定效应估计量是所有潜在因素的加权平均值， $2 \times 2$ DD 估计权重均基于组大小和处理方差。在方差加权共同趋势 (VWCT) 和时间不变处理效果的假设下，方差加权 ATT 是所有可能的 ATT 的加权平均值。

在更严格的假设下，该估计与 ATT 完全匹配。但当存在时变处理效果时，情况并非如此，因为双向固定效应估计的微分时序设计中的时变处理效果会产生偏差。因此，双向固定效应模型可能存在严重偏差。

假设在此设计中有三组：早期处理组 ( $k$ )，后来处理的一组 ( $l$ )，以及一个从未被处理过的组别 ( $U$ )，组 $k$ 和 $l$ 相似之处在于它们都得到了处理，但它们的不同之处在于 $k$ 处理早于 $l$ 。

以下三个公式分别代表上述三个组别相互比对所得出的 $2 \times 2$ ：

{\hat{δ}}_{k U}^{2 \times 2} = | {\bar{u}}_{k}^{post (k)} - {\bar{u}}_{k}^{pre (k)} | - | {\bar{u}}_{U}^{post (k)} - {\bar{u}}_{U}^{pre (k)} |

{\hat{δ}}_{k l}^{2 \times 2} = | {\bar{u}}_{k}^{m i d (k . l)} - {\bar{u}}_{k}^{pre (k)} | - | {\bar{u}}_{l}^{m i d (k . l)} - {\bar{u}}_{l}^{pre (k)} |

{\hat{δ}}_{l k}^{2 \times 2} = | {\bar{u}}_{l}^{post (l)} - {\bar{u}}_{l}^{mid (k . l)} | - | {\bar{u}}_{k}^{post (l)} - {\bar{u}}_{k}^{m i d (k . l)} |

考虑到这种表示法，DD 参数估计可以分解如下：

{\hat{δ}}^{D D} = \sum_{k \neq U} s_{k U} {\hat{δ}}_{k U}^{2 \times 2} + \sum_{k \neq U} \sum_{l > k} s_{k l} | μ_{k l} {\hat{δ}}_{k l}^{2 \times 2 . k} + (1 - μ_{k l}) {\hat{δ}}_{k l}^{2 \times 2 . l} |

而该式中的权重系数表示如下：

s_{k u} = \frac{n_{k} n_{u} {\bar{D}}_{k} (1 - {\bar{D}}_{k})}{\hat{Var} ({\tilde{D}}_{i t})}

s_{k l} = \frac{n_{k} n_{l} ({\bar{D}}_{k} - {\bar{D}}_{l}) (1 - ({\bar{D}}_{k} - {\bar{D}}_{l}))}{\hat{Var} ({\tilde{D}}_{i t})}

μ_{k l} = \frac{1 - {\bar{D}}_{k}}{1 - ({\bar{D}}_{k} - {\bar{D}}_{l})}

根据以上分析，可以得出：延长或缩短面板实际上可以纯粹通过改变组处理方差来改变点估计。那么，多种修正后的估计方法所解决的究竟是实验组处理时间不同所带来的差异，还是由于二者处理效用不同带来的影响，亦或是二者兼而有之。本文将对这一问题进行解答。

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 指定第一阶段公式，包括用于估计 $Y_{i} t (0)$ 的固定效应和协变量，不能放入处理变量。两阶段双重差分的原理是先识别组别效应和时期效应，然后在第二阶段将其剔除后，再对处理变量进行回归。主要用于解决的问题是处理组个体接受处理的时间是交错的，而且平均处理效应随着组别以及时间发生变化的情况。

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

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