最新Sun和Abraham(2020)和TWFE估计多期或交错DID并绘图展示结果！详细解读code！

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

邮箱：[e-mail]

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

关于DID的最新进展，参看：1.120篇DID双重差分方法的文章合集, 包括代码,程序及解读, 建议收藏！2.诚实双重差分法DID, 面板事件研究法和Bacon分解的经典应用文！3.前沿: 多期或渐进或交叠DID, 如何进行平行趋势检验呢？4.多期DID或渐进DID或交叠DID, 最新Stata执行命令整理如下供大家学习，5.DID前沿: 5种方法估计事件研究的因果效应, 并使用绘制系数和置信区间, 详细代码和数据，6.事件研究法开展政策评估和因果识别, 分享8篇提供数据和代码的文章，7.推荐用渐进(多期)DID和事件研究法开展政策评估的论文及其实现数据和代码！8.机器学习已经与政策评估方法, 例如事件研究法结合起来识别政策因果效应了！9.多期DID前沿方法大讨论, e.g., 进入-退出型DID, 异质性和动态性处理效应DID, 基期选择问题等,10.面板数据政策评估方法(DID等)的最新进展与相关应用, 包括机器学习和因子估计法，11.DID相关前沿问题“政策交错执行+堆叠DID+事件研究”, 附完整slides，12.前沿, 合成双重差分法SDID方法介绍和示例, 附code和数据！13.前沿, 模糊双重差分法FDID方法介绍和示例, 附code和数据！

正文

接着“合成控制法SCM新进展, 直接通过代码演示, 附数据和code！”，今天，通过详细示例和代码注释等方式给大家分享“如何用TWFE和最新 Sun和Abraham (2020)估计多期或交错DID并绘图展示结果”。下面的内容简单易懂，无论你是有无DID相关基础，根据下面的数据和code都可以看懂这一方法的运行过程。

*完整版code放在计量社群里的，可以自行下载运行。

① TWFE估计多期DID

*1.使用TWFE估计多期DID或者说交错DID，代码如下，以下结果经过调试都可以直接运行*

*导入数据

use bacon_example.dta, clear

*为受到政策影响的州创建滞后/领先项，所谓的lag/lead项，对于未受到政策影响的州用0填充，这允许在每个州发生 `treat` 和 `time_to_treat` 之间的交互，否则，可能会有一些NA值出现。

g time_to_treat = year - _nfd  //_nfd为政策发生的年份，这里就是时间上的动态项

replace time_to_treat = 0 if missing(_nfd)

g treat = !missing(_nfd)  //这个州是否受到政策影响

如下，这个州在1971年受到政策影响：

如下，这个州1984年受到政策影响：

如下，这个州未受到政策影响：

*通过构建shifted_tt变量，把time_to_treat为负数的变成非负数

summ time_to_treat

g shifted_ttt = time_to_treat - r(min)

*使用reghdfe进行回归，同时控制年份和州固定效应，在州级别进行聚类

reghdfe asmrs treat#shifted_ttt pcinc asmrh cases, a(stfips year) vce(cluster stfips)

*获得回归系数和标准误

g se = .

g coef = .

levelsof shifted_ttt, l(times)

foreach t in `times' {

replace coef = _b[1.treat#`t'.shifted_ttt] if shifted_ttt == `t'

replace se = _se[1.treat#`t'.shifted_ttt] if shifted_ttt == `t'

}

*构建置信区间

g ci_top = coef+1.96*se

g ci_bottom = coef - 1.96*se

*根据time_to_treat把重复了的剔除掉

keep time_to_treat coef se ci_*

duplicates drop

sort time_to_treat

* 创建连接的系数散点图，使用 rcap 中包含的 CI和水平和垂直方向为0的线

summ ci_top

local top_range = r(max)

summ ci_bottom

local bottom_range = r(min)

twoway (sc coef time_to_treat, connect(line)) ///

(rcap ci_top ci_bottom time_to_treat)///

(function y = 0, range(time_to_treat)) ///

(function y = 0, range(`bottom_range' `top_range') horiz), ///

xtitle("Time to Treatment") caption("95% Confidence Intervals Shown")

②  Sun 和 Abraham (2020)估计多期DID

当然，如前所述，这种TWFE分析受到 Sun 和 Abraham (2020) 的批评。我们使用 Stata 中的 eventstudyinteract 命令估计 Sun 和 Abraham建议的方法。先自行安装 eventstudyinteract命令。

*Sun and Abraham (2020)的方法进行估计的代码如下，都已经经过调试，可以直接使用，并给予了注释。

*导入数据

use bacon_example.dta, clear

*为受到政策影响的州创建滞后/领先项，所谓的lag/lead项

* 对于未收到政策影响的州用0填充

* 这允许在每个州发生 `treat` 和 `time_to_treat` 之间的交互

* 否则，可能会有一些NA出现，并且估计值将被关闭

g time_to_treat = year - _nfd

replace time_to_treat = 0 if missing(_nfd)

* 这将确定差异

* btw 控制和处理状态

g treat = !missing(_nfd)

g never_treat = missing(_nfd)

* 手动为处理组创建相对时间指标，这里只用了-9到16区间，你可以将全部区间-21到27用起来

*由于缺乏观测值而忽略遥远的领先和滞后项

*（请注意，假设这些领先/滞后项之外的任何影响为 0）

tab time_to_treat

forvalues t = -9(1)16 {

if `t' < -1 {

local tname = abs(`t')  //令tname为t的绝对值，免得生成变量名称为g_m-1

g g_m`tname' = time_to_treat == `t'  //当time_to_treat等于1，则生成g_m1 = 1

}

else if `t' >= 0 {

g g_`t' = time_to_treat == `t'

}

}

eventstudyinteract asmrs g_*, cohort(_nfd) control_cohort(never_treat) covariates(pcinc asmrh cases) absorb(i.stfips i.year) vce(cluster stfips)   //cohort()分类变量，对应于每个单位的初始政策处理时间。如果有单位接受多次政策处理，Sun and Abraham (2020)将初始政策处理时间定义为基于第一次处理。对于从未受政策处理过的单位，应将此分类变量设置为缺失。

//control_cohort()二元变量，可以是从未接受政策处理的单位或最后接受处理的单位。如果使用最后一个接受政策处理的单位作为对照队列，则排除最后一个队列接受政策处理的时间段。

* 获取政策效应和绘图

* 在撰写本文时，系数矩阵未标记，因此我们不能使用 _b[] 和 _se[]

* 相反，我们将使用结果表

matrix T = r(table)

g coef = 0 if time_to_treat == -1

g se = 0 if time_to_treat == -1

forvalues t = -9(1)16 {

if `t' < -1 {

local tname = abs(`t')

replace coef = T[1,colnumb(T,"g_m`tname'")] if time_to_treat == `t'

replace se = T[2,colnumb(T,"g_m`tname'")] if time_to_treat == `t'

}

else if `t' >= 0 {

replace coef =  T[1,colnumb(T,"g_`t'")] if time_to_treat == `t'

replace se = T[2,colnumb(T,"g_`t'")] if time_to_treat == `t'

}

}

* 获取置信区间

g ci_top = coef+1.96*se

g ci_bottom = coef - 1.96*se

*根据time_to_treat把重复了的剔除掉

keep time_to_treat coef se ci_*

duplicates drop

sort time_to_treat

keep if inrange(time_to_treat, -9, 16)

* 创建连接的系数散点图，使用 rcap 中包含的 CI和水平和垂直方向为 0 的线

summ ci_top

local top_range = r(max)

summ ci_bottom

local bottom_range = r(min)

twoway (sc coef time_to_treat, connect(line)) ///

(rcap ci_top ci_bottom time_to_treat)///

(function y = 0, range(time_to_treat)) ///

(function y = 0, range(`bottom_range' `top_range') horiz), ///

xtitle("Time to Treatment with Sun and Abraham (2020) Estimation") caption("95% Confidence Intervals Shown")

*完整版code放在计量社群里的，可以自行下载运行。

**R软件操作代码

library(data.table)

wranglinglibrary(fixest)    

dat = fread("https://raw.githubusercontent.com/LOST-STATS/LOST-STATS.github.io/master/Model_Estimation/Data/Event_Study_DiD/bacon_example.csv")

dat[, treat := ifelse(is.na(_nfd), 0, 1)]

dat[, time_to_treat := ifelse(treat==1, year - _nfd, 0)]

mod_twfe = feols(asmrs ~ i(time_to_treat, treat, ref = -1) + pcinc + asmrh + cases | stfips + year,  cluster = ~stfips,  data = dat)

iplot(mod_twfe,  xlab = 'Time to treatment',  main = 'Event study: Staggered treatment (TWFE)')

---

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

2.5年，计量经济圈近1000篇不重类计量文章，

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

Econometrics Circle

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