【应用计量系列20】DID设计与双向固定效应（TWFE）估计量
All is Study, All is Research
许文立，CIMERS，cimers_dsge@econmod.cn
全文约3045字
预计需12分钟阅读全文
今天推送的内容是我在安徽大学经济学院开设的课程《科研写作创新训练》中有关“双重差分(DID)”的节选。更详细的内容，可以免费下载我的《应用计量经济学讲稿》。
应用计量系列请参见：经验分析方法及Stata命令汇总
正
文
1、经典的2×2DID设计
我们首先来看看我的家乡——武汉的高等教育。大家可能有所耳闻，英雄的城市武汉是中国的高等教育重镇，坐拥众多全球知名、全国一流的高等学府，其中全国最美的大学——武汉大学就在珞珈山下，东湖畔，在武大学习可以漫步在樱花大道，欣赏东湖和磨山的风景，聆听古编钟奏响的乐曲，经过张之洞、周恩来、张培刚等伟人和学术巨擎们的生活旧址与铜像前，时刻会想起那声振聋发聩地宣言”为中华之崛起而读书“！
21世纪之初，武汉很多高校由于扩招而在武汉郊区假设新校区。假设2004年，武汉大学在江夏区建设了新校区，我们想看看高校新区的经济（旅游）效应如何。这个时候，我们就选择一个与江夏区比较类似的武汉郊区，例如，我的故乡黄陂区（花木兰的故乡）——没有建设高校新区——来进行对照。这个时候就形成了经典的2×2DID：两个时期——建设新校区前（2003年）和后（2004年）、两个组群——江夏区和黄陂区。
（图1 武汉赏花地图。注：图片来源于网络，如侵权请联系删除）
我们定义两个虚拟变量，是否建设新校区NU（=1，是；=0，否），建设时间前后Period（=1，后；=0，前），这个时候，我们就可以写出经典的双向固定效应（TWFE）模型：
其中，y表示结果变量——经济发展，β表示回归系数。上述模型也可以用下列表格来表示：
| NU = 0 | NU = 1 | 差分 | |
|---|---|---|---|
| Period = 0 | β0 | β0+β2 | β2 |
| Period = 1 | β0+β1 | β0+β1+β2+β3 | β2+β3 |
| 差分 | β1 | β1+β3 | β3 |
上述表格右下角的β3就是两次差分后的结果，也是我们关注的TWFE估计量。下面，我们用模拟数据来看看stata的估计结果。
首先，我们用stata生产一套2×2DID的数据集。也就是数据结构为两期、两个组群的面板数据。
clear // 清除stata已存在的数据local units = 2local start = 1local end = 2
local time = `end' - `start' + 1local obsv = `units' * `time'set obs `obsv'
egen NU = seq(), b(`time') egen Period = seq(), f(`start') t(`end')
sort NU Periodxtset NU Period // 声明面板数据类型
lab var NU 'Panel variable'lab var Period 'Time variable'第二步，定义处理变量T和结果变量y：
* 创建处理变量T和结果变量ygen T = NU==2 & Period==2 //双等号表示恒等于，即处理发生在第二个id，第二期，这时T=1，其它所有情况T=0gen btrue = cond(T==1, 3, 0) //cond表示当T==1为真时，btru这个变量才赋值为3，否则赋值为0gen Y = NU + 3*Period + btrue*T //利用这个函数来生成结果变量y的数据
从y的数据生成函数中，我们可以很容易的看出，在建新校区后，也就是2004年黄陂区和江夏区的结果Y的差异为4。
下面，我们用图形来看看这个数据：
lab de prepost1 1 '前' 2 '后'lab val Period prepost1
twoway /// (connected Y Period if NU==1) /// (connected Y Period if NU==2) /// , /// legend(order(1 'NU=黄陂区' 2 'NU=江夏区')) /// xlabel(1 2, valuelabel) ylabel(4(1)11)
从上图我们可以看到，蓝色的线条（黄陂区）和红色的线条（江夏区）在新建校区后的差异为4，而在新校区建设前的差异为1。因此，建设新校区前后的净变化量应该为4-1=3，即新建校区对当地经济的影响效应为3。
下面，我们可以用面板数据回归来看看结果：
* 面板数据回归xtreg Y T Period,fe
我们可以很清晰地看到，T的回归的系数就是3。
也可以使用另一种双向固定效应命令：
* 常用的双向固定效应命令reghdfe* 如果没有安装这个程序包，请先安装：* ssc install reghdfe,replace
reghdfe Y T, absorb(NU Period)从双向固定效应模型的回归结果也可以看出，建设新校区对经济的效应为3。
2、推广到多期的Periods×2DID设计
我们将经典的2×2推广到多个时期，例如，2000-2019年的20年。多期DID的数据：
*********************************************** 二、多期的Periods×2DID**********************************************clear // 清除stata已存在的数据local units = 2local start = 1local end = 20local time = `end' - `start' + 1local obsv = `units' * `time'set obs `obsv'egen NU = seq(), b(`time')egen Period = seq(), f(`start') t(`end')sort NU Periodxtset NU Period // 声明面板数据类型lab var NU 'Panel variable'lab var Period 'Time variable'
假设第二个个体（江夏)在第5期新建校：
* 假设在第5期之后，有新校区假设
* 创建处理变量T和结果变量ygen T = NU==2 & Period>=5 //双等号表示恒等于，即处理发生在第二个id，第二期，这时T=1，其它所有情况T=0
lab var T 'Treated'
gen btrue = cond(T==1, 4, 0) //cond表示当T==1为真时，btru这个变量才赋值为3，否则赋值为0 gen Y = NU + Period + btrue*T //利用这个函数来生成结果变量y的数据
lab var Y 'Outcome variable'我们画出多期DID的图：
从这个图可以看出，在第5期前，黄陂区和江夏区的差异是1，而第5期之后，这个差异扩大到了5，因此，新建校区的效应是4。
我们用回归也可以得到同样的结果：
* 面板数据回归xtreg Y T Period,fe* 常用的双向固定效应命令reghdfe* 如果没有安装这个程序包，请先安装：* ssc install reghdfe,replacereghdfe Y T, absorb(NU Period)
3、多期-多组群，相同处理期的Periods×NU DID设计
下面，我们来看看只有一个对照组，另个处理组，两个处理组受到处理的时间相同，但是，可能处理强度有所差异：
*********************************************** 三、多个组群-多期的Periods×NU DID*********************************************** 同时处理，但处理强度不同
clear // 清除stata已存在的数据local units = 3local start = 1local end = 20
local time = `end' - `start' + 1local obsv = `units' * `time'set obs `obsv'
egen NU = seq(), b(`time') egen Period = seq(), f(`start') t(`end')
sort NU Periodxtset NU Period // 声明面板数据类型
lab var NU 'Panel variable'lab var Period 'Time variable'
* 假设在第5期之后，有新校区假设
* 创建处理变量T和结果变量ygen T = 0replace T = 1 if NU>=2 & Period>=5 //处理发生在第二和三个id，第5期，这时T=1，其它所有情况T=0
lab var T 'Treated'
cap drop Y
gen Y = 0
replace Y = cond(T==1, 2, 0) if NU==2
replace Y = cond(T==1, 4, 0) if NU==3
lab var Y 'Outcome variable'
twoway /// (connected Y Period if NU==1) /// (connected Y Period if NU==2) /// (connected Y Period if NU==3) /// , /// xline(4.5) /// xlabel(1(1)20) /// legend(order(1 'NU=1' 2 'NU=2' 3 'NU=3'))
* 面板数据回归
xtreg Y T Period,fe * 常用的双向固定效应命令reghdfe* 如果没有安装这个程序包，请先安装：* ssc install reghdfe,replace
reghdfe Y T, absorb(NU Period)
图形如下所示。从图中可以卡到，对于个体2（红色线）来说，处理后的效应时2，而对于个体3（绿色线）来说，处理后的效应为4。这就意味着，平均处理效应(ATT)为3，即2和3的处理效应的平均值，(2+4)/2=3，在最终的ATT里，两个处理组的权重都是0.5。
我们用回归来看看ATT：
* 面板数据回归xtreg Y T Period,fe* 常用的双向固定效应命令reghdfe* 如果没有安装这个程序包，请先安装：* ssc install reghdfe,replacereghdfe Y T, absorb(NU Period)
上面的例子，很容易就得到ATT，这是因为它没有个体固定效应和时间固定效应，我们从图里就可以看出，这ATT等于3。如果我们加入一些控制变量，例如时间趋势和个体固定效应，我们在来看看结果。
* 首先，舍弃上例中生成的数据，重新来生成Y的数据cap drop Y
gen Y = 0
replace Y = NU + Period + cond(T==1, 0, 0) if NU==1replace Y = NU + Period + cond(T==1, 2, 0) if NU==2replace Y = NU + Period + cond(T==1, 4, 0) if NU==3
lab var Y 'Outcome variable'
* 可视化数据
twoway /// (connected Y Period if NU==1) /// (connected Y Period if NU==2) /// (connected Y Period if NU==3) /// , /// xline(4.5) /// xlabel(1(1)10) /// legend(order(1 'NU=1' 2 'NU=2' 3 'NU=3'))结果如下图所示，同样的三个组群，两个在第5期接受处理（2和3），一个对照组。但是，这个时候的平均处理效应为多少呢？好像从图里不能直观的看出来，下面我们就用回归来得到结果。
* 面板数据回归xtreg Y T Period,fe* 常用的双向固定效应命令reghdfe* 如果没有安装这个程序包，请先安装：* ssc install reghdfe,replacereghdfe Y T, absorb(NU Period)
上述两个回归的结果都得到ATT=3。我们来比较一下在回归中控制平行趋势假设——控制固定效应和不控制固定效应的差异：
* 我们在上述回归中控制平行趋势假设
reg Y T // 不控制任何因素reg Y T i.Period // 仅仅控制时间固定效应reg Y T i.NU // 仅仅控制个体固定效应reg Y T i.Period i.NU // 双向固定效应(正确!)大家可以去比较一下上面的结果差异。
4、多期-多组群，不相同处理期的Periods×NU DID设计
其实，在现实经济中，很多政策都是在不同时间开始实施的。下面我们就来看看这种情形——不同处理时点的DID设计。下面，我们还是用模拟数据来演示：
*********************************************** 四、多个组群-多期、不同处理时点的Periods×NU DID**********************************************clear // 清除stata已存在的数据local units = 3local start = 1local end = 20local time = `end' - `start' + 1local obsv = `units' * `time'set obs `obsv'egen NU = seq(), b(`time')egen Period = seq(), f(`start') t(`end')sort NU Periodxtset NU Period // 声明面板数据类型lab var NU 'Panel variable'lab var Period 'Time variable'* 假设在第5期之后，有新校区假设* 创建处理变量T和结果变量ygen T = 0replace T = 1 if NU==2 & Period>=5 //处理发生在第二个id，第5期，这时T=1，其它所有情况T=0replace T = 1 if NU==3 & Period>=10 //处理发生在第三个id，第10期，这时T=1，其它所有情况T=0lab var T 'Treated'cap drop Ygen Y = 0replace Y = cond(T==1, 2, 0) if NU==2replace Y = cond(T==1, 4, 0) if NU==3lab var Y 'Outcome variable'twoway ///(connected Y Period if NU==1) ///(connected Y Period if NU==2) ///(connected Y Period if NU==3) ///, ///xline(4.5) ///xlabel(1(1)20) ///legend(order(1 'NU=1' 2 'NU=2' 3 'NU=3'))
图形结果如下：
上述图形结果显示，第2个个体处理发生在第5期，处理后的效应为2，而第3个个体的处理发生在第10期，处理后的效应为4。那这个时候的平均处理效应ATT还是3吗？现在可以肯定的是，ATT不是3，如果是3的话，那不就与上面的第三节的效应相同吗，可是此时的处理时点并不是如同第三节都在第5期发生处理呀，所以两个ATT肯定有所差异。如果不是，那应该是多少呢？
我们用回归来看看，ATT到底是多少？
* 面板数据回归
xtreg Y T Period,fe * 常用的双向固定效应命令reghdfe* 如果没有安装这个程序包，请先安装：* ssc install reghdfe,replace
reghdfe Y T, absorb(NU Period)我们从上述回归结果可以看出，此时的ATT=3.21。
下面，我考虑一下个体因素和时间趋势。
* 加入一些时间趋势和个体效应cap drop Ygen Y = 0replace Y = NU + Period + cond(T==1, 0, 0) if NU==1replace Y = NU + Period + cond(T==1, 2, 0) if NU==2replace Y = NU + Period + cond(T==1, 4, 0) if NU==3lab var Y 'Outcome variable'* 可视化数据twoway ///(connected Y Period if NU==1) ///(connected Y Period if NU==2) ///(connected Y Period if NU==3) ///, ///xline(4.5) ///xlabel(1(1)10) ///legend(order(1 'NU=1' 2 'NU=2' 3 'NU=3'))
这个时候的ATT=？
我们在用回归来看看结果：
* 面板数据回归
xtreg Y T Period,fe * 常用的双向固定效应命令reghdfe* 如果没有安装这个程序包，请先安装：* ssc install reghdfe,replace
reghdfe Y T, absorb(NU Period)回归结果显示，ATT=3.29。
从上面两个结果可以看出，它们都不同于第三节的ATT=3。因此，肯定有什么不一样的地方。我们再稍微仔细想想。哦！原来此时处理时点不同，个体2是在第5期受到的处理，而个体3是在第10期受到的处理。那么，我们现在就不能清晰地区分出“处理前”和“处理后”，也不能很好地划分“处理组”和“对照组”。例如，在第5-10期，只有NU=2这个组群发生了变化，其它两个组群都保持不变，而在第10期之后，只NU=3组群发生了变化，其它两个组群也没有变。
那么，我们应该如何理解不同处理时点的ATT=3.29这个估计值呢？
Goodman-Bacon（2020，JoE）给了我们答案，且听下回分解！
注：（1）更多内容，请阅读《应用计量经济学讲稿》，在主页下载；
（2）上述案例的dofile请在主页下载。
参考文献：
A. Naqvi(2021): The classic 2x2 DiD or the Twoway Fixed Effects Model (TWFE).
●图片贴纸来源 | 秀米IUMI