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

Stata连享会   主页 || 视频 || 推文 || 知乎

温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装命令如下：
  . ssc install lianxh
详情参见帮助文件 (有惊喜)：
  . help lianxh

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

• 全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
• 计量专题 | 论文写作 | 数据分享 | 专题课程
• 结果输出 | Stata绘图 | 数据处理 | Stata程序
• 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
• 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
• Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
• 文本分析-爬虫 | Python-R-Matlab | 机器学习
• Markdown | 工具软件 | 其它

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

目录

本文来自 连享会专题 - 倍分法 DID

这是连享会「倍分法(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，并且每年增加3
bysort 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 <= 20
bysort id: replace y1 = 10 + 5*x1 + 3*x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40
bysort 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 r
predict e, ue
binscatter e time, line(connect) by(G)

///输出生成的图片，令格式为800*600
graph 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-多维固定效应。四个命令的比较放在了下方的表格中。

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

HDFE Linear regression                              Number of obs   =	   600
Absorbing 2 HDFE groups                             F(   3,    528) = 46468.92
Statistics 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, r
eststo reg
xtreg y c.D x1 x2 i.time, r fe
eststo xtreg_fe
areg y c.D x1 x2 i.time, absorb(id) robust
eststo areg
reghdfe 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	           600   
Groups|                          60		                 
+------------------------------------------------------------------------+
* p<0.05, ** p<0.01, *** p<0.001
+------------------------------------------------------------------------+
			

$D_{it}$ 的四个估计系数都为 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，并且每年增加3
bysort 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 <= 20
bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40
bysort 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 <= -4
tab 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   =	   600
Absorbing 2 HDFE groups                        F(  11,    520) =   76749.93
Statistics 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	.5082091
eventt3 | -.0993291   .2633063    -0.38   0.706    -.6166038	.4179457
eventt4 | -.1192429   .2381015    -0.50   0.617    -.5870021	.3485162
eventt5 |  2.555801   .2667126     9.58   0.000     2.031834	3.079768
eventt6 |  6.008552    .258689    23.23   0.000     5.500348	6.516756
eventt7 |  8.697876   .2586339    33.63   0.000     8.189781	9.205972
eventt8 |  11.72037   .3015466    38.87   0.000     11.12797	12.31277
eventt9 |  14.94957   .3702594    40.38   0.000     14.22218	15.67696
eventt10|  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   =	   600
Absorbing 2 HDFE groups	                          F(  11,    520) = 75233.93
Statistics 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.184053
c.D#c.year3 |   .6093975   .3984095     1.53	0.127    -.1732924	1.392087
c.D#c.year4 |   .4808495   .3948783     1.22	0.224    -.2949033	1.256602
c.D#c.year5 |   .1168801   .4088713     0.29	0.775    -.6863626	.9201227
c.D#c.year6 |   23.81018   .3870237    61.52	0.000     23.04986	24.5705
c.D#c.year7 |   28.48194   .3664986    77.71	0.000     27.76194	29.20194
c.D#c.year8 |   31.9992   .3978656     80.43	0.000     31.21758	32.78082
c.D#c.year9 |   36.2474   .4087051     88.69	0.000     35.44448	37.05031
c.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 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

///生成协变量x1, x2

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

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

gen D = 0 
replace D = 1 if  id > 29
gen post = 0
replace 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 == 1
tab 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   =	   600
Absorbing 2 HDFE groups                           F(  11,    520) =   75233.93
Statistics 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.184053
eventt3 |  .6093975   .3984095     1.53   0.127    -.1732924	1.392087
eventt4 |  .4808495   .3948783     1.22   0.224    -.2949033	1.256602
eventt5 |  .1168801   .4088713     0.29   0.775    -.6863626	.9201227
eventt6 |  23.81018   .3870237    61.52   0.000     23.04986	24.5705
eventt7 |  28.48194   .3664986    77.71   0.000     27.76194	29.20194
eventt8 |   31.9992   .3978656    80.43   0.000     31.21758	32.78082
eventt9 |   36.2474   .4087051    88.69   0.000     35.44448	37.05031
eventt10|  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. 总结和拓展

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

这也就构成本系列文章的第一阶段。第一篇文章主要展示 Standard DID 模型和 ESA 方法结合的思路，第二篇和第三篇文章分别对两种情况（是否存在一直处于控制组的个体）下 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. 参考资料

1. Stata: reghdfe-多维固定效应
2. 开学礼包：如何使用双重差分法的交叉项（迄今最全攻略）
3. 连享会专题 - 倍分法 DID
4. 连享会-倍分法DID详解 (一)：传统 DID
5. 倍分法DID详解 (二)：多时点 DID (渐进DID)
6. 倍分法DID详解 (三)：多时点 DID (渐进DID) 的进一步分析

附：文中相关代码

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，并且每年增加3
bysort 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 <= 20
bysort id: replace y1 = 10 + 5*x1 + 3*x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40
bysort 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 r
predict e, ue
binscatter e time, line(connect) by(G)

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

reg y c.D x1 x2 i.time i.id, r
eststo reg
xtreg y c.D x1 x2 i.time, r fe
eststo xtreg_fe
areg y c.D x1 x2 i.time, absorb(id) robust
eststo areg
reghdfe 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，并且每年增加3
bysort 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 <= 20
bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth_date + 1) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40
bysort 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 <= -4
tab 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 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

///生成协变量x1, x2

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

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


gen D = 0 
replace D = 1 if  id > 29
gen post = 0
replace 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 == 1
tab 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)

相关课程

连享会-直播课 上线了！
http://lianxh.duanshu.com

免费公开课：

• 直击面板数据模型 - 连玉君，时长：1小时40分钟
• Stata 33 讲 - 连玉君, 每讲 15 分钟.
• 部分直播课 课程资料下载 (PPT，dofiles等)

课程一览

支持回看，所有课程可以随时购买观看。

Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。
• 连享会-主页 和 知乎专栏，300+ 推文，实证分析不再抓狂。
• 公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

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

✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD