DID会固定年份吗_倍分法DID详解 (三)：多时点 DID (渐进DID) 的进一步分析

作者：王昆仑 (天津大学)
E-mail: shawn0513@163.com

连享会专题课程：DSGE 模型及应用

这是连享会「倍分法(DID)专题推文」系列的第三篇文章。前两篇文章分别对多种 DID 模型的估计，平行趋势的检验以及政策的动态效果的展示，通过模拟的方式给出了较为详尽的解答。

具体而言，在第一篇推文「连享会-倍分法DID详解 (一)：传统 DID」中，我们对具有相同政策时点的 标准倍分法 (Standard DID) 模型进行了详细的介绍，第二篇推文「倍分法DID详解 (二)：多时点 DID (渐进DID)」又对政策时点存在差异，实验组个体又不断变动的 时变倍分法 (Time-varying DID) 模型进行了系统的介绍。

本文作为本系列的第三篇文章，我们继续使用同一套数据结构，利用模拟的方法将 Standard DID 和 Time-varying DID 方法统一到一起，继续分析关于 Time-varying DID 方法剩下的一些待解决的问题。

本系列的推文

连享会-倍分法DID详解 (一)：传统 DID

倍分法DID详解 (二)：多时点 DID (渐进DID)

倍分法DID详解 (三)：多时点 DID (渐进DID) 的进一步分析 (本文)

1. 引言

本文讨论的背景依然是 Time-varying DID 模型，为了解决当政策干预的时点和群体不断发生变化时，Standard DID 的交互项设置会导致估计系数有偏的问题。在「倍分法DID详解 (二)：多时点 DID (渐进DID)」的末尾，我们提出了如果样本中存在一直处于控制组的个体，那么 Time-varying DID 模型的设置又是否会发生的问题。由于政策效果是否随时间变化不会改变我们的分析，出于行文简洁的考虑，本文中我们就以政策效果随时间发生变化的情况来对上述的问题进行解答。

2. Time-varyig DID Simulation: 政策效果随时间发生改变

2.1 模拟数据的生成

再次仿照「连享会-倍分法DID详解 (一)：传统 DID」生成基础的数据结构，依然为60个体*10年=600个观察值的平衡面板数据。本文的 Time-varying DID 设置为, id 编号为 1-20 的个体在 2004 年接收政策干预，编号 21-40 的个体在 2006 年接受干预，编号为 41-60 一直不接受干预。因此，三组个体接受政策干预的时长分别为 6 年，4 年和 0 年。

///设定60个观测值，设定随机数种子. clear all. set obs 60. set seed 10101. gen id =_n

/// 每一个数值的数量扩大11倍，再减去前六十个观测值，即60*11-60 = 600，为60个体10年的面板数据. expand 11. drop in 1/60. count

///以id分组生成时间标识. bys id: gen time = _n+1999. xtset id time

///生成协变量以及个体和时间效应. gen x1 = rnormal(1,7). gen x2 = rnormal(2,5)

. sort time id. by time: gen ind = _n. sort id time. by id: gen T = _n. gen y = 0

///生成处理变量,此时D为Dit,设定1-20在2004年接受冲击，21-40为2006年，41-60个体一直不接受干预

. gen D = 0. gen birth_date = 0. gen G = 0

forvalues i = 1/20{ replace D = 1 if id  == `i' & time >= 2004   replace birth_date = 2004 if id == `i'    replace  G = 1 if id == `i'}

forvalues i = 21/40{ replace D = 1 if id  == `i' & time >= 2006   replace birth_date = 2006 if id == `i'    replace  G = 2 if id == `i'}

///将基础数据结构保存成dta文件，命名为DID_Basic_Simu.dta,默认保存在当前的 working directory 路径下

save "DID_Basic_Simu_2.dta", replace

生成结果变量，设定接受干预个体的政策效果当年为 3，且每年增加 3。利用固定效应模型消除个体效应和 x1、x2 两个变量对结果变量的影响，得到残差，画出更加干净分组的结果变量的时间趋势图。这个图形可以作为具有平行趋势的一个旁证。

///调用生成的基础数据文件use "DID_Basic_Simu_2.dta", clear

///Y的生成，设定真实的政策效果为当年为3，并且每年增加3bysort id: gen     y0 = 10 + 5*x1 + 3*x2 + T + ind + rnormal()bysort id: gen     y1 = 10 + 5*x1 + 3*x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5*x1 + 3*x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5*x1 + 3*x2 + T + ind + rnormal() if y1 == .

replace y = y0 + D * (y1 - y0)

///去除个体效应和协变量对Y的影响，得到残差并画图xtreg y x1 x2 , fe rpredict e, uebinscatter e time, line(connect) by(G)

///输出生成的图片，令格式为800*600graph export "article_3.png", as(png) replace width(800) height(600)

2.2 多种估计结果的展示

依然使用双向固定效应(Two-Way Fixed Effects)的方法进行估计，使用reg,xtreg,areg,reghdfe等四个Stata命令进行估计。在本文中，主要展示命令 'reghdfe'的输出结果，该命令的具体介绍可以参考 Stata: reghdfe-多维固定效应。四个命令的比较放在了下方的表格中。

此处依然使用双向固定效应 (two-way fixed effects) 方法估计交互项的系数，其中用到 reg, xtreg, areg, reghdfe 等多个 Stata 命令。在本文中，主要展示命令 reghdfe 的输出结果，该命令的具体介绍可以参考 Stata: reghdfe-多维固定效应。四个命令的比较放在了下方的表格中。

	reg	xtreg	areg	reghdfe
个体固定效应	i.id	xtreg,fe	areg,absorb(id)	absorb(id time)
时间固定效应	i.time	i.time	i.time	absorb(id time)
估计方法	OLS	组内去平均后OLS	OLS	OLS
优点	命令熟悉，逻辑清晰	固定效应模型的官方命令	官方命令，可以提高组别不随样本规模增加的估计效率	高维固定效应模型，可以极大提到估计效率，且选项多样，如支持多维聚类
缺点	运行速度慢，结果呈现过多不太需要的固定效应的结果	需要手动添加时间固定效应	需要手动添加时间固定效应	无

. reghdfe y c.D x1 x2, absorb(id time) vce(robust)+-----------------------------------------------------------------------------+   (converged in 3 iterations)

HDFE Linear regression                              Number of obs   =    600Absorbing 2 HDFE groups                             F(   3,    528) = 46468.92Statistics robust to heteroskedasticity             Prob > F        =  0.0000                                                    R-squared       =  0.9972                                                    Adj R-squared   =  0.9969                                                    Within R-sq.    =  0.9963                                                    Root MSE        =  2.3658+-----------------------------------------------------------------------------+                 |              Robust        y |     Coef.   Std. Err.      t     P>t     [95% Conf.   Interval]+-----------------------------------------------------------------------------+                D |  6.933194   .3673665    18.87   0.000     6.211515    7.654873       x1 |  4.986501    .014677   339.75   0.000     4.957669    5.015334       x2 |  3.025778   .0201695   150.02   0.000     2.986155    3.0654+-----------------------------------------------------------------------------+   

Absorbed degrees of freedom:+--------------------------------------------------------+  Absorbed FE|   Num. Coefs.  =   Categories  -   Redundant+--------------------------------------------------------+          id |          60              60              0      time |           9              10              1+--------------------------------------------------------+

reg y c.D x1 x2 i.time i.id, reststo regxtreg y c.D x1 x2 i.time, r feeststo xtreg_feareg y c.D x1 x2 i.time, absorb(id) robusteststo aregreghdfe y c.D x1 x2, absorb(id time) vce(robust)eststo reghdfe

estout *, title("The Comparison of Actual Parameter Values") ///         cells(b(star fmt(%9.3f)) se(par)) ///         stats(N N_g, fmt(%9.0f %9.0g) label(N Groups)) ///        legend collabels(none) varlabels(_cons Constant) keep(x1 x2 D)

+------------------------------------------------------------------------+                      The Comparison of Actual Parameter Values+------------------------------------------------------------------------+            |         reg           xtreg_fe       areg         reghdfe+------------------------------------------------------------------------+           D |       6.933***        6.933***      6.933***      6.933***      |      (0.367)         (0.450)     (0.367)     (0.367)    x1|       4.987***        4.987***    4.987***      4.987***      |      (0.015)         (0.015)   (0.015)           (0.015)    x2|       3.026***        3.026***    3.026***      3.026***      |      (0.020)         (0.020)   (0.020)        (0.020)+------------------------------------------------------------------------+       N |         600             600                 600              600Groups|                          60                      +------------------------------------------------------------------------+* p<0.05, ** p<0.01, *** p<0.001+------------------------------------------------------------------------+

的四个估计系数都为 6.993，每年的平均处理效应见第三章的图形。从表格展示的结果可以知道，四个命令估计的系数大小是一致的，唯一的区别在于固定效应模型的估计系数其标准误略有变化。再次验证了这四个命令在估计 Two-Way Fixed Effects 上是等价的。

3. 一直存在控制组个体的 Time-varyig DID 和 ESA 方法的结合

根据本文第二节的结果可以知道，是否一直存在控制组个体对于 Time-varying DID 模型的单方程设定没有任何的影响。因此，在本节，我们将进入本文的关键问题————是否一直存在控制组个体会对 ESA 方法的方程设定产生影响吗？如果两种方程设定存在差异，那么其来源在何处呢？

关于 Time-Varying DID 和 ESA 结合的思路与 Standard DID 和 ESA 方法结合的思路的异同，我相信已经在本系列的第二篇文章中讲述清楚。具体的内容可以参考「倍分法DID详解 (二)：多时点 DID (渐进DID)」的 2.3 节。

在设定相对时间值的时候，我们会发现如果存在一直在控制组的个体，那么对于这些个体的相对时间进行赋值就是一个难题——因为它一直没有接受干预，那么自然就不存在有政策前 N 期和政策后 N 期。因此无论将其赋值成什么都会造成问题，那么将此类个体的相对时间设定成缺失值是不是可以呢？答案是不可以。因为将这些个体的相对时间值设定为缺失值，会导致这些个体的相对时间值的虚拟变量也会是缺失值。这进一步会使得 ESA 方程可用的个体仅剩下了样本中的实验组个体，那么自然也无法实现我们想要的估计出每年度的政策效果，进而对平行趋势进行检验的目的。

那么如何来解决这个问题呢？最简单的办法就是在生成每一个个体的相对时间值的虚拟变量后，我们可以手动将虚拟变量中的缺失值赋值为0。这样就可以在回归中利用上这些个体。另外，将这些个体的缺失值赋值为0还可以表示它们是控制组个体。

(Note: 如果从双重差分的角度去理解，那么我们所生成的相对时间值的虚拟变量中的0的含义是只有表示为控制组个体这一个含义吗？这个问题也留给大家进行思考)

///调用生成的基础数据文件use "DID_Basic_Simu_2.dta", clear

///Y的生成，设定真实的政策效果为当年为3，并且每年增加3bysort id:     gen y0 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal()bysort id:     gen y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal() if y1 == .

replace y = y0 + D * (y1 - y0)

///Time-varying DID 和 Event Study Approach 的结合///用当前年份减去个体接受处理的年份，得到相对的时间值 event，将 -4 期之前的时间归并到第 -4 期，由于部分个体没有多于 -4 期的时间///由于存在一直处于控制组的个体，因此 event 变量对于这些个体为缺失值///然后生成相对时间值的虚拟变量 eventt，对于一直存在控制组的个体 eventt 虚拟变量的值全取为 0 ，并将首期设定为基准对照组

gen event = time - birth_date if id <= 40

replace event = -4 if event <= -4tab event, gen(eventt)

forvalues i = 1/10{    replace eventt`i' = 0 if eventt`i' == .}

drop eventt1

. reghdfe y eventt* x1 x2, absorb(id time) vce(robust)+--------------------------------------------------------------------------+             (converged in 3 iterations)

HDFE Linear regression                         Number of obs   =     600Absorbing 2 HDFE groups                        F(  11,    520) =   76749.93Statistics robust to heteroskedasticity        Prob > F        =  0.0000                                               R-squared       =   0.9995                                               Adj R-squared   =   0.9995                                               Within R-sq.    =   0.9994                                               Root MSE        =   0.9782+--------------------------------------------------------------------------+                          |               Robust      y |   Coef.      Std. Err.      t    P>t     [95% Conf.  Interval]+--------------------------------------------------------------------------+                   eventt2 |  .0406161   .2380169     0.17   0.865    -.4269769  .5082091eventt3 | -.0993291   .2633063    -0.38   0.706    -.6166038  .4179457eventt4 | -.1192429   .2381015    -0.50   0.617    -.5870021  .3485162eventt5 |  2.555801   .2667126     9.58   0.000     2.031834  3.079768eventt6 |  6.008552    .258689    23.23   0.000     5.500348  6.516756eventt7 |  8.697876   .2586339    33.63   0.000     8.189781  9.205972eventt8 |  11.72037   .3015466    38.87   0.000     11.12797  12.31277eventt9 |  14.94957   .3702594    40.38   0.000     14.22218  15.67696eventt10|  17.99397   .3880296    46.37   0.000     17.23167  18.75626     x1 |  4.988479   .0057955   860.75   0.000     4.977094  4.999865     x2 |  3.006429   .0085394   352.07   0.000     2.989653  3.023205+--------------------------------------------------------------------------+                

Absorbed degrees of freedom:+-------------------------------------------------------------+         

Absorbed FE|  Num. Coefs.  =   Categories  -   Redundant+-------------------------------------------------------------+                    id |          60              60              0      time |           9              10              1+-------------------------------------------------------------+           

///图示法coefplot, ///   keep(eventt*)  ///   coeflabels(eventt2 = "-3"   ///   eventt3 = "-2"              ///   eventt4 = "-1"              ///   eventt5 = "0"               ///   eventt6  = "1"              ///   eventt7  = "2"              ///   eventt8  = "3"              ///   eventt9  = "4"              ///   eventt10 = "5")             ///   vertical                    ///   yline(0)                    ///   ytitle("Coef")              ///   xtitle("Time passage relative to year of adoption of implied contract exception") ///   addplot(line @b @at)        ///   ciopts(recast(rcap))        ///   ///rescale(100)             ///   scheme(s1mono)

///输出生成的图片，令格式为800*600>graph export "article3_3.png",as(png) replace width(800) height(600)

4. 稳健性检验：以 Standard DID 为例

尽管第三节的处理十分的直观，而且画出的估计系数的动态变化也与预期相符合，但是我们可能依然对上述的操作存在疑虑，所以本文增加了一个稳健性检验，以保证本文第三节的处理方法是合理的。

之所以选取 Standard DID 模型的模拟作为稳健性检验的原因有二：

其一，Standard DID 从另一个角度可以认为是一种一直有控制组个体的 DID 模型，尽管此时不是 Time-varying DID 模型，但是从相对时间值这个角度来说二者相通的；

其二，Standard DID 的动态政策效果我们已经可以通过成熟的办法来得到，那么就等于我们预先知道了答案，从而有了明确的标准对第三节的处理合适与否进行判断。

本节首先呈现第一篇文章中 5.2 节的估计结果和图 4，作为估计结果比较的基础。

///第一篇文章中 5.2 节的估计结果reghdfe y c.D#(c.year2-year10) x1 x2 , absorb(id time) vce(robust)+--------------------------------------------------------------------------+       (converged in 3 iteration

HDFE Linear regression                             Number of obs   =       600Absorbing 2 HDFE groups                           F(  11,    520) = 75233.93Statistics robust to heteroskedasticity            Prob > F        = 0.0000                                                  R-squared       =    0.9996                                                  Adj R-squared   =    0.9996                                                  Within R-sq.    =    0.9994                                                  Root MSE        =    1.0147+--------------------------------------------------------------------------+                          |              Robust          y |     Coef.   Std. Err.      t     P>t     [95% Conf.   Interval]       +--------------------------------------------------------------------------+        c.D#c.year2 |   .4381916   .3796627     1.15  0.249    -.3076697  1.184053c.D#c.year3 |   .6093975   .3984095     1.53  0.127    -.1732924  1.392087c.D#c.year4 |   .4808495   .3948783     1.22  0.224    -.2949033  1.256602c.D#c.year5 |   .1168801   .4088713     0.29  0.775    -.6863626  .9201227c.D#c.year6 |   23.81018   .3870237    61.52  0.000     23.04986  24.5705c.D#c.year7 |   28.48194   .3664986    77.71   0.000     27.76194  29.20194c.D#c.year8 |   31.9992   .3978656     80.43  0.000     31.21758  32.78082c.D#c.year9 |   36.2474   .4087051     88.69  0.000     35.44448  37.05031c.D#c.year10|   40.45248   .3979999   101.64  0.000     39.67059  41.23436          x1|   4.996797   .0061877   807.54  0.000     4.984641  5.008953          x2|   3.004127   .0087679   342.63  0.000     2.986902  3.021352+--------------------------------------------------------------------------+        

Absorbed degrees of freedom:+-------------------------------------------------------------+           Absorbed FE|   Num. Coefs.  =   Categories   Redundant+-------------------------------------------------------------+                id |          60              60          0      time |           9              10         1+-------------------------------------------------------------+

首先，仿造第一篇文章的 5.2 节，生成结构为60个体*10年共600个观测值的平衡面板数据。我们在数据中设定 2005 年为政策发生当年，id 30-60 的个体为处理组，剩余的个体形成控制组。

///设定60个观测值，设定随机数种子

clear allset obs 60set seed 10101gen id =_n

/// 每一个数值的数量扩大11倍，再减去前六十个观测值，即60*11-60 = 600，为60个体10年的面板数据

expand 11drop in 1/60count

///以id分组生成时间标识

bys id: gen time = _n+1999xtset id time

///生成协变量x1, x2

gen x1 = rnormal(1,7)gen x2 = rnormal(2,5)

///生成个体固定效应和时间固定效应

sort time idby time: gen ind = _nsort id timeby id: gen T = _n

///生成treat和post变量，以2005年为接受政策干预的时点，id为30-60的个体为处理组，其余为控制组

gen D = 0replace D = 1 if  id > 29gen post = 0replace post = 1 if time >= 2005

///将基础数据结构保存成dta文件，命名为DID_Basic_Simu.dta,默认保存在当前的 working directory 路径下

save "DID_Basic_Simu.dta"，replace

///调用刚生成的基础数据结构use "DID_Basic_Simu.dta", clear

///生成两种潜在结果，并且合成最终的结果变量，令政策的真实效果随时间发生变化，即(5*T-T)，由于从2005年开始接受干预，因此，每年的政策效果应为24，28，32，36，40.

bysort id: gen     y0 = 10 + 5*x1 + 3*x2 + T + ind  + rnormal()bysort id: gen     y1 = 10 + 5*x1 + 3*x2 + T + ind  + rnormal() if time < 2005

bysort id: replace y1 = 10 + 5*x1 + 3*x2 + T * 5 + ind + rnormal() if time >= 2005

gen y = y0 + D * (y1 - y0)

/// Standard DID 仿造本文第三节生成相对时间值/// 对于一直存在控制组的个体 eventt 虚拟变量的值全取为 0 ，并将首期设定为基准对照gen event = time - 2005 if D == 1tab event, gen(eventt)

forvalues i = 1/10{  replace eventt`i' = 0 if eventt`i' == .}

drop eventt1

. reghdfe y eventt* x1 x2, absorb(id time) vce(robust)+-----------------------------------------------------------------------------+(converged in 3 iterations)+-----------------------------------------------------------------------------+HDFE Linear regression                            Number of obs   =       600Absorbing 2 HDFE groups                           F(  11,    520) =   75233.93Statistics robust to heteroskedasticity             Prob > F        =  0.0000                                                    R-squared       =  0.9996                                                    Adj R-squared   =  0.9996                                                    Within R-sq.    =  0.9994                                                    Root MSE        =  1.0147+-----------------------------------------------------------------------------+        |               Robust      y |     Coef.    Std. Err.      t    P>t     [95% Conf. Interval]+-----------------------------------------------------------------------------+eventt2 |  .4381916   .3796627     1.15   0.249    -.3076697  1.184053eventt3 |  .6093975   .3984095     1.53   0.127    -.1732924  1.392087eventt4 |  .4808495   .3948783     1.22   0.224    -.2949033  1.256602eventt5 |  .1168801   .4088713     0.29   0.775    -.6863626  .9201227eventt6 |  23.81018   .3870237    61.52   0.000     23.04986  24.5705eventt7 |  28.48194   .3664986    77.71   0.000     27.76194   29.20194eventt8 |   31.9992   .3978656    80.43   0.000     31.21758  32.78082eventt9 |   36.2474   .4087051    88.69   0.000     35.44448  37.05031eventt10|  40.45248   .3979999   101.64   0.000     39.67059  41.23436     x1 |  4.996797   .0061877   807.54   0.000     4.984641  5.008953     x2 |  3.004127   .0087679   342.63   0.000     2.986902  3.021352+-----------------------------------------------------------------------------+

Absorbed degrees of freedom:+--------------------------------------------------------+   Absorbed FE|  Num. Coefs.  =   Categories  -   Redundant+--------------------------------------------------------+           id |          60              60              0      time |           9              10              1+--------------------------------------------------------+    

///图示法coefplot, ///   keep(eventt*)  ///   coeflabels(eventt2 = "-4"   ///   eventt3 = "-3"              ///   eventt4 = "-2"              ///   eventt5 = "-1"              ///   eventt6  = "0"              ///   eventt7  = "1"              ///   eventt8  = "2"              ///   eventt9  = "3"              ///   eventt10 = "4")             ///   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)

///输出生成的图片，令格式为800*600>graph export "article3_3png",as(png) replace width(800) height(600)

观察本节通过生成相对时间值的办法对 Standard DID 的 ESA 方程进行改写的结果可以知道，本文的结果与第一篇文章的图 4 的结果完全一致。因此，可以让我们更加确信本文处理方法的正确性。

5. 总结和拓展

这是本系列的第三篇文章，因此是时候对这三篇文章做一个小结了。

这也就构成本系列文章的第一阶段。第一篇文章「连享会-倍分法DID详解 (一)：传统 DID」主要展示 Standard DID 模型和 ESA 方法结合的思路，第二篇「倍分法DID详解 (二)：多时点 DID (渐进DID)」和第三篇文章分别对两种情况(是否存在一直处于控制组的个体)下 Time-Varying DID 模型和 ESA 方法的结合做了介绍。这三种情况下 ESA 所构建的回归方程的处理略有差别，这也是文中着重在讲的问题。

另外，Standard DID 模型和 ESA 方法的结合思路看似与 Time-varying DID 方法格格不入，但是上篇文章已经从获得相对时间值这个角度对二者的联系进行了说明，本文想再从第三篇文章讨论的问题对二者的联系展开一定的阐释：从一直存在控制组个体这一点上看，本文所讨论的的对象和 Standard DID 模型是一致的，也就是说，恰好成熟的 Standard DID 模型和 ESA 方法的结合给本文提供了一个可参照的可靠结果，让我们对 Time-Varying DID 的处理更加具有信心。因此，关于以上这两点，三篇文章是有机统一的整体，作为本系列推文的第一阶段恰到好处。

在陈强老师的推文中，几乎将常见的 DID 模型都进行了介绍，包括而不限于本系列的前三篇文章所讨论的类型。本系列后续的文章依然会在 DID 模型的大背景下，通过模拟的方法介绍相关的内容，比如介绍其他形式的 DID 模型，判断遗漏变量问题，DID 模型的系数与加权回归的关系以及如何处理不满足平行趋势的 DID 模型等等问题。

6. 参考资料

Stata: reghdfe-多维固定效应
开学礼包：如何使用双重差分法的交叉项(迄今最全攻略)
连享会-倍分法DID详解 (一)：传统 DID
倍分法DID详解 (二)：多时点 DID (渐进DID)

Stata连享会计量专题 || 公众号合集

附：文中相关代码

A.1 基础数据结构的生成和基础回归结果

///设定60个观测值，设定随机数种子. clear all. set obs 60. set seed 10101. gen id =_n

/// 每一个数值的数量扩大11倍，再减去前六十个观测值，即60*11-60 = 600，为60个体10年的面板数据. expand 11. drop in 1/60. count

///以id分组生成时间标识. bys id: gen time = _n+1999. xtset id time

///生成协变量以及个体和时间效应. gen x1 = rnormal(1,7). gen x2 = rnormal(2,5)

. sort time id. by time: gen ind = _n. sort id time

. by id: gen T = _n. gen y = 0

///生成处理变量,此时D为Dit,设定1-20在2004年接受冲击，21-40为2006年，41-60个体一直不接受干预

. gen D = 0. gen birth_date = 0. gen G = 0

forvalues i = 1/20{    replace D = 1 if id  == `i' & time >= 2004   replace birth_date = 2004 if id == `i'   replace  G = 1 if id == `i'}

forvalues i = 21/40{  replace D = 1 if id  == `i' & time >= 2006   replace birth_date = 2006 if id == `i'   replace  G = 2 if id == `i'}

///将基础数据结构保存成dta文件，命名为DID_Basic_Simu.dta,默认保存在当前的 working directory 路径下

save "DID_Basic_Simu_2.dta", replace  // 保存模拟数据

///调用生成的基础数据文件use "DID_Basic_Simu_2.dta", clear

///Y的生成，设定真实的政策效果为当年为3，并且每年增加3bysort id: gen     y0 = 10 + 5*x1 + 3*x2 + T + ind + rnormal()bysort id: gen     y1 = 10 + 5*x1 + 3*x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5*x1 + 3*x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5*x1 + 3*x2 + T + ind + rnormal() if y1 == .

replace y = y0 + D * (y1 - y0)

///去除个体效应和协变量对Y的影响，得到残差并画图xtreg y x1 x2 , fe rpredict e, uebinscatter e time, line(connect) by(G)

///输出生成的图片，令格式为800*600graph export "article_3.png",as(png) replace width(800) height(600)

reg y c.D x1 x2 i.time i.id, reststo regxtreg y c.D x1 x2 i.time, r feeststo xtreg_feareg y c.D x1 x2 i.time, absorb(id) robusteststo aregreghdfe y c.D x1 x2, absorb(id time) vce(robust)eststo reghdfe

estout *, title("The Comparison of Actual Parameter Values") ///        cells(b(star fmt(%9.3f)) se(par)) ///         stats(N N_g, fmt(%9.0f %9.0g) label(N Groups)) ///        legend collabels(none) varlabels(_cons Constant) keep(x1 x2 D)

A.2 一直存在控制组个体的 Time-varying DID 和 ESA 方法的结合

///调用生成的基础数据文件use "DID_Basic_Simu_2.dta", clear

///Y的生成，设定真实的政策效果为当年为3，并且每年增加3bysort id:     gen y0 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal()bysort id:     gen y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal() if y1 == .

replace y = y0 + D * (y1 - y0)

///Time-varying DID 和 Event Study Approach 的结合///用当前年份减去个体接受处理的年份，得到相对的时间值 event，将 -4 期之前的时间归并到第 -4 期，由于部分个体没有多于 -4 期的时间///由于存在一直处于控制组的个体，因此 event 变量对于这些个体为缺失值///然后生成相对时间值的虚拟变量 eventt，对于一直存在控制组的个体 eventt 虚拟变量的值全取为 0 ，并将首期设定为基准对照组

gen event = time - birth_date if id <= 40

replace event = -4 if event <= -4tab event, gen(eventt)

forvalues i = 1/10{ replace eventt`i' = 0 if eventt`i' == .}

drop eventt1

reghdfe y eventt* x1 x2, absorb(id time) vce(robust)

///图示法coefplot, ///   keep(eventt*)  ///   coeflabels(eventt2 = "-3"   ///   eventt3 = "-2"              ///   eventt4 = "-1"              ///   eventt5 = "0"               ///   eventt6  = "1"              ///   eventt7  = "2"              ///   eventt8  = "3"              ///   eventt9  = "4"              ///   eventt10 = "5")             ///   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)

///输出生成的图片，令格式为800*600>graph export "article3_3.png",as(png) replace width(800) height(600)

A3. 稳健性检验部分的代码

///设定60个观测值，设定随机数种子clear allset obs 60set seed 10101gen id =_n

/// 每一个数值的数量扩大11倍，再减去前六十个观测值，即60*11-60 = 600，为60个体10年的面板数据

expand 11drop in 1/60count

///以id分组生成时间标识

bys id: gen time = _n+1999xtset id time

///生成协变量x1, x2

gen x1 = rnormal(1,7)gen x2 = rnormal(2,5)

///生成个体固定效应和时间固定效应

sort time idby time: gen ind = _nsort id timeby id: gen T = _n

///生成treat和post变量，以2005年为接受政策干预的时点，id为30-60的个体为处理组，其余为控制组

gen D = 0replace D = 1 if  id > 29gen post = 0replace post = 1 if time >= 2005

///将基础数据结构保存成dta文件，命名为DID_Basic_Simu.dta,默认保存在当前的 working directory 路径下

save "DID_Basic_Simu.dta", replace

///调用刚生成的基础数据结构use "DID_Basic_Simu.dta", clear

///生成两种潜在结果，并且合成最终的结果变量，令政策的真实效果随时间发生变化，即(5*T-T)，由于从2005年开始接受干预，因此，每年的政策效果应为24，28，32，36，40.

bysort id: gen     y0 = 10 + 5 * x1 + 3 * x2 + T + ind  + rnormal()bysort id: gen     y1 = 10 + 5 * x1 + 3 * x2 + T + ind  + rnormal() if time < 2005

bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T * 5 + ind   + rnormal() if time >= 2005

gen y = y0 + D * (y1 - y0)

/// Standard DID 仿造本文第三节生成相对时间值///对于一直存在控制组的个体 eventt 虚拟变量的值全取为 0 ，并将首期设定为基准对照组gen event = time - 2005 if D == 1tab event, gen(eventt)

forvalues i = 1/10{    replace eventt`i' = 0 if eventt`i' == .}

drop eventt1

. reghdfe y eventt* x1 x2, absorb(id time) vce(robust)

///图示法coefplot, ///   keep(eventt*)  ///   coeflabels(eventt2 = "-4"   ///   eventt3 = "-3"              ///   eventt4 = "-2"              ///   eventt5 = "-1"              ///   eventt6  = "0"              ///   eventt7  = "1"              ///   eventt8  = "2"              ///   eventt9  = "3"              ///   eventt10 = "4")             ///   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)

///输出生成的图片，令格式为800*600>graph export "article3_3png", as(png) replace width(800) height(600)

连享会专题课程：DSGE 模型及应用

DID会固定年份吗_倍分法DID详解 (三)：多时点 DID (渐进DID) 的进一步分析

连享会专题课程：DSGE 模型及应用

1. 引言

2. Time-varyig DID Simulation: 政策效果随时间发生改变

2.1 模拟数据的生成

2.2 多种估计结果的展示

3. 一直存在控制组个体的 Time-varyig DID 和 ESA 方法的结合

4. 稳健性检验：以 Standard DID 为例

5. 总结和拓展

6. 参考资料

附：文中相关代码

A.1 基础数据结构的生成和基础回归结果

A.2 一直存在控制组个体的 Time-varying DID 和 ESA 方法的结合

A3. 稳健性检验部分的代码

连享会专题课程：DSGE 模型及应用

DID会固定年份吗_倍分法DID详解 (三)：多时点 DID (渐进DID) 的进一步分析相关推荐

最新文章

热门文章