作者：李闯 (新疆大学)
邮箱：chuangli0808@163.com
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 引言
2. 模型介绍
3. 数据准备
4. 基准回归
5. 平行趋势检验
5.1 时间趋势图法
5.2 事件分析法
6. 相关推文
双重差分法 (DID) 作为目前政策评估方法中非常流行的一种方法，广泛应用于社科研究的各个领域。现有文献中关于该方法的应用多而零碎，给初学者的学习与选择带来了一定困难。因此，为了提升文献阅读与方法应用效率，系统学习 DID 是十分有必要的。
本推文将从基本模型介绍、前期数据准备、基准回归、平行趋势检验、安慰剂检验以及 PSM-DID 等整体流程入手，采用模型分析与实例结合的形式对多时点 DID 模型进行讲解。
我们以经典 DID 模型为例：
为政策虚拟变量，本质是 。 为个体虚拟变量，政策实施个体取值为 1，其余为 0； 为时间虚拟变量，政策实施当年及以后取值为 1，其他为 0。( 与 在数值上相同)
我们以白俊红等 (2022) 的论文数据为例，以国家创新型城市试点政策实施为准自然实验，展示如何用代码构造变量 。
需要注意的是该文中 的城市政策实施时间为 2010 年，但是 DID 却从 2012 开始计算。后续分析中本推文将其调整至 2010 年，所得结果可能与原文有细微差异，但不影响方法学习。
准备面板数据，即不含 相关变量在内的论文实证数据，需要包括个体 id，便于下文匹配。
准备政策实施样本，并构造政策实施时间变量： (便于区分原文变量，用 表示)，示例如下：
. lxhuse pt, clear
. list in 1/8
-------------
| city pt |
|-------------|
1. | 3 2010 |
2. | 4 2010 |
3. | 5 2011 |
4. | 14 2010 |
5. | 25 2011 |
|-------------|
6. | 26 2010 |
7. | 34 2010 |
8. | 35 2010 |
-------------
. lxhuse inno_policy, clear
. merge m:1 city using 'https://file.lianxh.cn/data/i/inno_pt.dta'
. drop _merge
. replace pt=0 if pt==. // 对照组样本policy=0
. g Treat=0
. replace Treat=1 if pt~=.
. g Post=0
. replace Post=1 if pt==2008&year>=2008
. replace Post=1 if pt==2010&year>=2010
. replace Post=1 if pt==2011&year>=2011
. replace Post=1 if pt==2012&year>=2012
. replace Post=1 if pt==2013&year>=2013
. replace Post=1 if pt==2018&year>=2018 // 对创新政策实施不同年份分别进行替换
. gen did=Treat*Post // Post与did在数值上相等
需要注意的是在进行数据匹配时最好采用 id，直接用个体汉字名称可能会出错。
. * 验证结果是否一致
. xtset city year
. global xlist 'lnagdp indust_stru finance ainternet market '
. xtreg entre_activation did $xlist i.year, fe robust
Fixed-effects (within) regression Number of obs = 4,200
Group variable: city Number of groups = 280
R-squared: Obs per group:
Within = 0.5500 min = 15
Between = 0.6964 avg = 15.0
Overall = 0.5535 max = 15
F(20,279) = 116.98
corr(u_i, Xb) = 0.3706 Prob > F = 0.0000
(Std. err. adjusted for 280 clusters in city)
------------------------------------------------------------------------------
| Robust
entre_acti~n | Coefficient std. err. t P>|t| [95% conf. interval]
------------- ----------------------------------------------------------------
did | 0.297 0.039 7.72 0.000 0.222 0.373
lnagdp | 0.264 0.143 1.84 0.067 -0.019 0.546
indust_stru | -0.569 0.186 -3.06 0.002 -0.935 -0.203
finance | 0.110 0.040 2.78 0.006 0.032 0.188
ainternet | 0.010 0.003 2.92 0.004 0.003 0.016
market | -0.019 0.009 -2.02 0.044 -0.037 -0.001
------------- ----------------------------------------------------------------
基于构造的变量进行回归后，得到与原数据相同结论，表明本文对数据的处理具有合理性。
DID 模型中的基本回归命令与面板回归命令基本一致，应用比较普遍的是 xtreg 与 reghdfe。对于 xtreg 命令的使用已较为成熟，此处不再赘述。reghdfe 命令可以解决多维固定效应问题，而且运行速度较快，因而备受青睐。以下为两种命令代码示例：
. xtreg entre_activation did $xlist i.year, fe robust
. xtreg gdpr did invest consume export gov second agg innov i.year, fe cluster(id)
. reghdfe gdpr did invest consume export gov second agg innov, ///
ab(province year) vce(cluster id)
. reghdfe gdpr did invest consume export gov second agg innov, ///
ab(i.province#i.year province year id) vce(cluster id)
时间趋势图可以粗略地对政策实施前后处理组与对照组样本之间的变化趋势做出判断，是进行 DID 研究的前提所在。
. graph set window fontface 'Times New Roman'
. graph set window fontfacesans '宋体' // 设置图形输出的字体
. egen mean_y=mean(entre_activation), by(year group)
. graph twoway (connect mean_y year if group==1,sort) ///
> (connect mean_y year if group==0,sort lpattern(dash)), ///
> xline(2010,lpattern(dash) lcolor(gray)) ///
> ytitle('{stSans:创}''{stSans:业}''{stSans:活}''{stSans:跃}''{stSans:度}', ///
> orientation(h)) xtitle('{stSans:年度}') ///
> ylabel(0(0.5)3,labsize(*0.75) format(%02.1f)) xlabel(,labsize(*0.75)) ///
> legend(label(1 '{stSans:处理组}') label( 2 '{stSans:控制组}')) ///图例
> xlabel(2004 (1) 2018) graphregion(color(white)) //白底
. graph export '平行趋势事前检验.png'
如上图所示，以 2010 年国家创新型城市试点为例，2010 年之后，试点城市 (处理组) 的创业活跃度相对于非试点城市 (对照组) 经历了更高程度的增长，二者的创业活跃度差距逐渐扩大，而在 2010 年之前，二者具有较为一致的变化趋势。
事件分析法检验平行趋势最常见的一种方法，一般模型构造如下：
时间虚拟变量 表示各城市确立为试点城市前 年、当年和后 年的观测值，非试点城市虚拟变量均为 0。白俊红等 (2022) 考虑到政策实施时间与试点数量等因素，取 。本推文提供以下三种检验方法，详细经济学含义请参考原文。
手动计算各时期虚拟变量，coefplot 辅助作图。
. * 代码一
. xtset city year
. gen event = year - branch_reform if group==1
. replace event = -4 if event <= -4
. tab event, gen(eventt)
. forvalues i = 1/15{
2. replace eventt`i' = 0 if eventt`i' == .
3. }
. drop eventt1
. xtreg entre_activation eventt* $xlist i.year, r fe
. coefplot, keep(eventt2 eventt3 eventt4 eventt5 eventt6 eventt7 eventt8 eventt9 ///
> eventt10 eventt11 eventt12 eventt13) coeflabels(eventt2 = '-3' eventt3 = '-2' ///
> eventt4 = '-1' eventt5 = '0' eventt6 = '1' eventt7 = '2' eventt8 = '3' ///
> eventt9 = '4' eventt10 = '5' eventt11 = '6' eventt12 = '7' eventt13 = '8') ///
> vertical yline(0) ytitle('coef') ///
> xtitle('time passage relative to year of adoption of implied contract exception') ///
> addplot(line @b @at) ciopts(recast(rcap)) scheme(s1mono)
. graph export 'article3_3.png',as(png) replace width(800) height(600)
. * 代码二
. gen distance=year-pt
. replace distance = -4 if distance <= -4
. tab distance,missing
. forvalues i=1/4{
2. gen d_`i' = 0
3. replace d_`i' = 1 if Treat== 1 & distance== -`i'
4. }
. forvalues i=1/10 {
2. gen d`i' = 0
3. replace d`i' = 1 if Treat== 1 & distance== `i'
4. }
. gen d0 = 0
. replace d0 = 1 if Treat== 1 & distance== 0
. xtset city year
. xtreg entre_activation d_3 d_2 d_1 d0 d1-d10 $xlist i.year,fe r
. graph set window fontface 'Times New Roman'
. graph set window fontfacesans '宋体' // 设置图形输出的字体
. coefplot, keep( d_3 d_2 d_1 d0 d1 d2 d3 d4 d5 d6 d7 d8) ///
> coeflabels( d_3='-3' d_2='-2' d_1='-1' d0='0' d1='1' ///
> d2='2' d3='3' d4='4' d5='5' d6='6' d7='7' d8='8' ) ///
> yline(0,lcolor(edkblue*0.8)) /// 加入 y=0 这条虚线
> ylabel(0(0.5)1.5,labsize(medsmall) format(%02.1f)) ///
> vertical addplot(line @b @at) ytitle('{stSans:系}' ///
> '{stSans:数}', size(medsmall) orientation(h)) ///
> xtitle('{stSans:国家创新型城市试点政策的平行趋势检验}', ///
> size(medsmall)) yline(0) levels(95) scheme(s2color) ///
> ciopts(recast(rcap) lpattern(dash))
. graph save '平行趋势.gph', replace
上述两种方法对于时间虚拟变量的计算结果完全相同，不同之处在于代码一在运用 coefplot 绘图时需要将 eventt* 转换为基于政策前后的数值，代码二则可以直接绘图，平行趋势图如下：
据上图所知，政策发生前的相对时间虚拟变量系数均不显著且数值较小，表明政策发生之前，处理组与对照组在创业活跃度上无显著差异，满足平行趋势假设。
进一步从动态效应来看，在政策发生后当年，试点政策短暂产生创业效应但还未稳定，试点政策推行 3 年后，创新型城市试点政策的影响系数显著为正并不断提升，表明创新型城市试点政策能够产生促进创业活跃度提升的政策效应，但具有一定的滞后性。
tvdiff 能够实现仅运行一条命令就可以同时实现图示处理期前后处理效应，以及呈现平行趋势检验统计结果。更多详细介绍，请参考连享会推文 Stata：多期倍分法 (DID) 详解及其图示。
* 命令安装
ssc install eventdd, replace
* 命令语法
eventdd depvar [indepvars] [ if ] [ in ] [ weight ],
timevar(timevar) ci(type, ...) [ method baseline(#)
level(#) accum lags(#) leads(#) noend noline keepbal(varname)
absorb(varname) wboot wboot_op(string) balanced inrange
graph_op(string) ci_op(string) coef_op(string) endpoints_op(string)]
其中，
depvar：被解释变量；indepvars：控制变量；timevar：必须，设定时间政策变量；ci(type， ...)：指定置信区间的图形类型，rarea (带区域阴影)、rcap (带封顶尖钉)、rline (带线条)，rcap 是默认的图表类型；method(type，[absorb(absvars)] * ...)：指定估计方法，默认为 ols，xtreg 对应 fe，hdfe 对应absorb(absvars)；baseline(#)：指定政策时间相对基准期，默认为 - 1；level(#)：设置置信度；默认值为 95%；lags(#)：指明事件研究中要考虑的最大事件后时段，在 accum、keepbal、inrange 指定时，必须确定；leads(#)：指明事件研究中要考虑的最大事件前期数，同上；accum：指定超出某些指定值的所有时期都应累积到最终点数中；noend：当指定 accum 选项时，要求累积终点不显示在图形输出上；keepbal(varname)：指定只保留面板中平衡的单位以供估计。这里varname 表示指示单位的面板变量 (如 state)；balanced：要求只绘制“平衡的”提前和滞后期；inrange：要求只绘制指定期数的提前和滞后期图形；noline：删除 -1 这条基准线；wboot：要求用 wild bootstrap 估计置信区间；wboot_op(string)：指定任意选项用于 wild bootstrap 估计 (例如 seed()、bootcluster())；graph_op(string)：允许包含 twoway 中允许的任何其他通用图形选项，包括标题、线型、轴标签等；coef_op(string)：允许包含 scatter 中的图形选项；endpoints_op：允许包含端点系数的 scatter 图形选项，须在 accum 下使用。下面我们使用作者提供的数据来检验此命令。作者使用了 Stevenson 和 Wolfers (2006) 关于美国无过错离婚改革和女性自杀的数据。数据包括 49 个州从 1964 年至 1996 年的面板数据，各州单方面离婚改革的时间不同。
. webuse bacon_example.dta, clear
. xtset stfips year
. xtreg asmrs post pcinc asmrh cases i.year, fe cluster(stfips)
. gen timeToTreat = year - _nfd // 设定时间政策变量
. eventdd asmrs pcinc asmrh cases i.year, timevar(timeToTreat) ///
> method(fe, cluster(stfips)) ///
> graph_op(ytitle('Suicidesper 1m Women') ///
> xlabel(-20(5)25))
. graph save 'a.gph', replace
. graph export 'Event Study Example Based on No-fault Divorce Reforms.png', replace
. eventdd asmrs pcinc asmrh cases i.year, fe timevar(timeToTreat) ///
> ci(rcap) inrange leads(10) lags(10) cluster(stfips) ///
> graph_op(ytitle('Suicides per 1m Women'))
. graph save 'b.gph', replace
. graph export 'Show periods between -10 y 10.png', replace
. eventdd asmrs pcinc asmrh cases i.year, fe timevar(timeToTreat) ///
> ci(rcap) balanced cluster(stfips) ///
> graph_op(ytitle('Suicides per 1m Women'))
. graph save 'c.gph', replace
. graph export 'Show balanced periods.png', replace
. eventdd asmrs pcinc asmrh cases, hdfe timevar(timeToTreat) ci(rcap) ///
> cluster(stfips) absorb(i.stfips i.year) keepbal(stfips) leads(10) ///
> lags(10) graph_op(ytitle('Suicidesper 1m Women'))
. graph save 'd.gph', replace
. graph export 'Only balanced units between -10 y 10.png', replace
. eventdd asmrs pcinc asmrh cases i.year, fe timevar(timeToTreat) ci(rcap) ///
> cluster(stfips) accum leads(15) lags(10) ///
> graph_op(ytitle('Suicidesper 1m Women'))
. graph save 'e.gph', replace
. graph export 'Accumulating into end points -15 y 10.png', replace
. eventdd asmrs pcinc asmrh cases i.year, timevar(timeToTreat) ///
> method(fe, cluster(stfips)) accum noend leads(15) lags(10) ///
> graph_op(ytitle('Suicidesper 1m Women') )
. graph save 'f.gph', replace
. graph export 'Not showing end points.png', replace
. graph combine a.gph b.gph c.gph d.gph e.gph f.gph, ///
> col(3) row(2) iscale(1) xsize(45) ysize(20)
上述对应的 6 组代码分别对应下图 a-f 的平行趋势检验情况，表明 eventdd 命令在运用事件分析法检验平行趋势时可以灵活处理不同时期，功能十分强大。
. eventdd asmrs pcinc asmrh cases i.year, fe timevar(timeToTreat) ci(rarea) ///
> cluster(stfips) accum leads(15) lags(20) graph_op(xlabel(-15(5)20) ///
> scheme(s1mono) ytitle('Suicides per 1m Women')) ci_op(fcolor(ltblue%45)) ///
> coef_op(msymbol(Oh)) endpoints_op(msymbol(O))
. graph export 'Plot with alternative appearance options.png', replace
下图展示了调整外观之后的平行趋势检验图，可以很直观看出在政策实施前，散点在红线附近波动，表明政策实施前无过错离婚改革和女性自杀具有一致变化趋势。
Note：产生如下推文列表的 Stata 命令为：
lianxh 绘图 did, m
安装最新版lianxh命令：
ssc install lianxh, replace
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
联系客服
