DID偏误问题:多时期DID的双重稳健估计量(下)-csdid

发布时间:2022-07-26 阅读 2525

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 系列推文:

PDF下载 - 推文合集

作者:张子楠 (浙江财经大学)
邮箱zinanzh@gmail.com

编者按:本文主要摘译自下文,特此致谢!
Source:Callaway B, Sant’Anna P H C. Difference-in-differences with multiple time periods[J]. Journal of Econometrics, 2021, 225(2): 200-230. -PDF-


目录


1. 绪论

在连享会推文「DID偏误问题:两时期DID的双重稳健估计量(上)-drdid」中,我们介绍了对于两时期 DID 模型,双重固定效应 (TWFE) 估计 DID 处理效应可能存在的偏误问题,以及如何通过双重稳健估计方法来缓解。此外,我们还介绍了如何通过 Sant 和 Zhao(2020) 提供的命令 drdid 来实现上述过程。

drdid 命令有两个不足:一是只适用两期 DID 情形,无法直接拓展到多期 DID 中;二是要求面板数据为强平衡数据,不能有缺失值。而这两点要求,在实践应用中都往往难以满足。特别是第一点,多时期 DID 要远比两时期 DID 更为常见。为了缓解上述问题,Callaway 和 Sant (2021) 提供了 drdid 命令的多时期版本——csdid

与其它多时期 DID 估计的纠偏估计命令类似,csdid 命令也是将样本分为不同的子组 (group),分别估计不同组别的 ATT(g),再通过特定策略将不同组别的 ATT(g) 加总算出样本期的 ATT。加总策略原则同样也是对那些可能存在偏误组的 ATT(g),降低它们的加总权重。

与其它纠偏估计量不同的是,csdid 是基于 drdid 的双重稳健估计思路,避免了 TWFE 估计量的问题。同时只需要 Outcome Regression 或 Inverse Propesnity Weight Regression 两种估计方法的依赖条件至少有一个满足,就可以显著缓解 DID 估计偏误问题。

csdid 命令还有一个优势在于,提供了多种加总不同子组 DID 回归结果的选项。除了常见的利用事件发生法去分析政策的动态效应以外,命令还允许我们计算不同子组处理效应的异质性,不同年份处理效应的异质性,以及处理效应的在一段时期的累计值。这不仅丰富了 DID 回归结果,还可以更为深入地了解 DID 总平均处理效应的底层来源。

注:Callaway 和 Sant Anna(2021) 同时分析了面板数据和混合截面数据两种数据结构下的双重稳健估计问题。本推文只介绍了面板数据部分,对混合截面数据相关内容感兴趣的读者可以自行阅读原文。

2. csdid 命令介绍

* 命令的安装
ssc install csdid, all replace
* 命令的语法
csdid depvar [indepvars] [if] [in] [weight], [ivar(varname)] time(varname) gvar(varname) [ options ]

对于面板数据,模型设定的选项包括以下四类:

  • depvar 为被解释变量,indepvars 为解释变量和控制变量;
  • ivar(varname)time (varname) 为面板设置选项。其中 ivar(varname) 设置个体标识变量,time (varname) 设置时间标识变量;
  • gvar(varname) 为多时期 DID 标识变量,和通常的多时期 DID 赋值不太一样。如果在样本期受到过政策冲击,赋值等于其受到冲击的期数序号。例如,在第 1 期受到冲击,赋值为 1。需要注意的是,如果某个个体在样本期开始前就已经受到冲击,即属于所谓的 always-treated 组,默认是不进入回归中;
  • notyet 要求对照组里不仅要有 not-yet 组里尚未受到冲击的样本,同时也要包括 always-treated 组样本。默认是不包括后者;
  • longlong2 对于冲击发生前的样本期的样本,使用 long gaps 而非 shortgaps;
  • 估计方法的选项包括:增进型双重稳健估计 (drimp)、基于 ipw 的双重稳健估计量 (dripw)、outcome regression 估计量 (reg)、标准化 ipw 估计量 (stdipw)、ipw 估计量 (ipw);
  • 标准误的选项默认是 robust and asymptotic standard errors,同时也可以选择 wbootcluster 等;
  • 不同组别政策效应加总 agg(aggtype)。加总方式 (aggtype) 包括:加总所有组所有时期 (simple)、每一组或者一群在所有时期加总 (group)、所有组在每一个时期分别加总 (calendar)、采用事件发生法在每一个时期分别加总 (event) 四大类。

3. csdid 基本用法

我们采用作者在帮助文档里提供的数据 mpdta.dta 来演示 csdid 的用法。

. use https://friosavila.github.io/playingwithstata/drdid/mpdta.dta, clear
. des

Contains data from https://friosavila.github.io/playingwithstata/drdid/mpdta.dta
 Observations:         2,500                  Written by R.              
    Variables:             6                  17 May 2021 11:45
-------------------------------------------------------------------------------------
Variable      Storage   Display    Value
    name         type    format    label      Variable label
-------------------------------------------------------------------------------------
year            int     %9.0g                 year
countyreal      long    %9.0g                 countyreal
lpop            double  %9.0g                 lpop
lemp            double  %9.0g                 lemp
first_treat     int     %9.0g                 first.treat
treat           byte    %9.0g                 treat
-------------------------------------------------------------------------------------

. tab  year

       year |      Freq.     Percent        Cum.
------------+-----------------------------------
       2003 |        500       20.00       20.00
       2004 |        500       20.00       40.00
       2005 |        500       20.00       60.00
       2006 |        500       20.00       80.00
       2007 |        500       20.00      100.00
------------+-----------------------------------
      Total |      2,500      100.00

. tab first_treat treat

first.trea |         treat
         t |         0          1 |     Total
-----------+----------------------+----------
         0 |     1,545          0 |     1,545 
      2004 |         0        100 |       100 
      2006 |         0        200 |       200 
      2007 |         0        655 |       655 
-----------+----------------------+----------
     Total |     1,545        955 |     2,500 

数据共包含 6 个变量,其中 year 表示年份、countyreal 表示个体、lpop 为被解释变量、lemp 为协变量。treat 为处理组指示变量,如果样本期内受到冲击,则为 1,否则为 0。first_treat 为个体第一次受到冲击的年份,取值包括 0,2004,2006 和 2007,0 表示样本期内没有受到冲击。观察上面结果可知,样本量为 2500,样本时期为 2003-2007,共五年。总共有 2500/5 个个体,其中分别有 100/5、200/5 和 655/5 数量的个体在 2004、2006 和 2007 年受到冲击。

接下来,我们用 csdid 命令来识别处理组平均处理效应 (ATT)。

. * 需要使用 agg(simple) 选项, 以估计总 ATT
. * method(dripw) 表示使用增进型双重稳健估计
. csdid lemp lpop , ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(simple)
                                                     
Outcome model  : least squares                          Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         ATT |     -0.042      0.012    -3.63   0.000       -0.064      -0.019
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details

此外,我们也用 TWFE 方法识别 ATT 来作为对比。

. * 生成 did 选项: 受到冲击的当期及以后为 1, 其它情况为 0
. gen did=0 
. replace did=1 if first_treat!=0 & year>= first_treat

. * TWFE 回归
. reghdfe lemp did lpop, absorb(countyreal year)  cluster(countyreal)

HDFE Linear regression                            Number of obs   =      2,500
Absorbing 2 HDFE groups                           F(   1,    499) =       7.59
Statistics robust to heteroskedasticity           Prob > F        =     0.0061
                                                  R-squared       =     0.9932
                                                  Adj R-squared   =     0.9915
                                                  Within R-sq.    =     0.0042
Number of clusters (countyreal) =        500      Root MSE        =     0.1391
                           (Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
             |               Robust
        lemp | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         did |     -0.037      0.013    -2.76   0.006       -0.063      -0.010
        lpop |      0.000  (omitted)
       _cons |      5.777      0.002  3741.28   0.000        5.774       5.780
------------------------------------------------------------------------------
Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
  countyreal |       500         500           0    *|
        year |         5           1           4     |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

对比可以发现,两种回归结果的回归系数和标准误都比较接近,置信区间基本重合。这说明对于这个例子,TWFE 估计偏误问题并不严重。注:有兴趣的读者可以用 Bacon 分解方法来理解背后原因。

4. csdid 进阶用法

除了更加稳健,以及能有更大概率获得无偏估计结果,csdid 命令的另一大优势在于,可以更为灵活地划分子组,并加总获得不同特征的 ATT,这在异质性分析等方面极大丰富 DID 估计结果。对这部分细节感兴趣的读者,可以阅读 Callaway 和 Sant Anna (2021) 的第 3 章,对此有详细的叙述。

接下来我们分别演示常见的几种加总获得样本总 ATT 的思路。第一、利用事件发生法的思路,以接受冲击时期的长短来来区分不同组别,估计事件冲击的动态效应。

. * 需使用 agg(event) 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(event)
                                                         
Outcome model  : least squares                          Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         T-3 |      0.027      0.014     1.90   0.057       -0.001       0.054
         T-2 |     -0.004      0.013    -0.28   0.780       -0.029       0.022
         T-1 |     -0.023      0.014    -1.60   0.109       -0.052       0.005
         T+0 |     -0.021      0.011    -1.83   0.067       -0.044       0.001
         T+1 |     -0.053      0.016    -3.24   0.001       -0.085      -0.021
         T+2 |     -0.140      0.035    -3.97   0.000       -0.210      -0.071
         T+3 |     -0.107      0.033    -3.25   0.001       -0.171      -0.042
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details

类似地,我们也用 TWFE 估计方法进行了对比,如下所示。由结果可知,对于我们这个例子,两种估计量的估计结果基本相一致。

. * 生成事件发生发所需要的变量
. gen didlength=year-first_treat
. replace didlength=0 if first_treat==0
. tab didlength,gen(event)
. drop event1  // 删除第一个 event 变量, 以和 csdid 命令处理策略一致

. * TWFE 回归
. reghdfe lemp event* lpop, absorb(countyreal year)  cluster(countyreal)

HDFE Linear regression                            Number of obs   =      2,500
Absorbing 2 HDFE groups                           F(   7,    499) =       3.60
Statistics robust to heteroskedasticity           Prob > F        =     0.0009
                                                  R-squared       =     0.9933
                                                  Adj R-squared   =     0.9915
                                                  Within R-sq.    =     0.0103
Number of clusters (countyreal) =        500      Root MSE        =     0.1388
                           (Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
             |               Robust
        lemp | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
      event2 |      0.021      0.015     1.39   0.166       -0.009       0.051
      event3 |      0.020      0.019     1.02   0.308       -0.018       0.058
      event4 |     -0.004      0.023    -0.16   0.877       -0.048       0.041
      event5 |     -0.022      0.025    -0.85   0.393       -0.072       0.028
      event6 |     -0.047      0.029    -1.63   0.103       -0.104       0.010
      event7 |     -0.135      0.039    -3.44   0.001       -0.213      -0.058
      event8 |     -0.096      0.042    -2.27   0.024       -0.179      -0.013
        lpop |      0.000  (omitted)
       _cons |      5.788      0.022   263.18   0.000        5.745       5.831
------------------------------------------------------------------------------
Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
  countyreal |       500         500           0    *|
        year |         5           1           4     |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation

第二、计算不同年份的所有组别平均效应。即对于某一年份的所有样本,其年平均处理效应是多少。需要注意的是,事件发生法是以受冲击的时期长度分组,这里是以某年份来区分不同的组别。

. * 需使用 agg(calendar) 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(calendar)
                                          
Outcome model  : least squares                          Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       T2004 |     -0.015      0.022    -0.66   0.511       -0.058       0.029
       T2005 |     -0.076      0.029    -2.67   0.008       -0.133      -0.020
       T2006 |     -0.046      0.021    -2.18   0.029       -0.088      -0.005
       T2007 |     -0.040      0.013    -3.06   0.002       -0.065      -0.014
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details

第三、计算不同组别的同一年份平均效应。在本示例里,事件冲击先后发生在 2004、2006 和 2007 这三年,因而有三个子组,分别用 G2004、G2006和 G2007 来标识。

. * 需使用 agg(group) 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(group)
                                                   
Outcome model  : least squares                           Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       G2004 |     -0.085      0.025    -3.44   0.001       -0.133      -0.036
       G2006 |     -0.020      0.017    -1.15   0.248       -0.054       0.014
       G2007 |     -0.029      0.016    -1.77   0.076       -0.061       0.003
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details

第四、不仅使用 Never Treated,也使用 Not yet Treated 作为控制组。需要注意的是,虽然回归结果表格下显示的是 "Control: Not yet Treated"。但根据帮助文件里选项设定的描述来看,这里面应该是同时使用了 Never Treated 和 Not yet Treated 来作为控制组的。

. * 需使用 notyet 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw)  agg(simple) notyet
                                                        
Outcome model  : least squares                           Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         ATT |     -0.041      0.011    -3.62   0.000       -0.064      -0.019
------------------------------------------------------------------------------
Control: Not yet Treated
See Callaway and Sant'Anna (2021) for details

第五、对于非平衡面板数据进行回归。在 drdid 命令中,是无法使用非平衡面板数据回归的。csdid 则在这方面拓展了双重稳健估计量在 DID 分析中的应用。

. * 生成随机数 sample, 赋值为 0 或者 1, 均值为 0.9
. set seed 1
. gen sample = runiform()<.9

. * 需使用 if sample==1 选项来挑选部分样本回归, 以实现类似存在缺失值的效果
. csdid lemp lpop if sample==1, cluster(countyreal) time(year) gvar(first_treat) method(dripw) agg(simple)
                                                         
Outcome model  : least squares                            Number of obs = 2,245
Treatment model: inverse probability
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         ATT |     -0.052      0.018    -2.86   0.004       -0.087      -0.016
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details

5. 进一步讨论

在帮助文件里,作者还给出了一些注意事项,主要有以下三点:

  • drdid 命令类似,csdid 假定协变量不随时间变化。当使用面板数据回归时,即使协变量是时变的,csdid 命令也只使用协变量基期值来估计;
  • 要使用 csdid 命令来获得 DID 无偏估计量,一般只能使用 never-treated 和 (或) not-yet-treated 的数据作为控制组,而不能使用 always-treated 组。否则由于异质性处理效应,平行趋势假设并不能满足,从而估计结果仍旧有偏;
  • 当使用面板数据时,与 drdid 命令不同,csdid 命令并没有要求数据为强平衡面板数据,即并没有要求没有缺失值。但是,在估计不同组的 ATT 时,存在缺失值的样本并不会进入回归。这和 TWFE 里对于缺失值的处理方式是一致的。

6. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh did, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh