【应用计量系列131】不平行的DD，平行的DDD与stata应用

过年期间，我写过“未来几年，三重差分的理论进展或许是面板数据领域最重要的发展方向之一” 【应用计量系列130】春节假期读的DID新论文列表 。虽然DDD与DD都是在1994年AER上出现的，但是，DD已经经历了三波革命性的发展（传统的DD，事件研究DD和异质性处理效应DD），而DDD似乎并没有这么热闹，在应用文献中，甚至成了DD的扩展。

DDD的应用实际上非常广泛，包括我自己的文章（彭飞,许文立,吕鹏等.未预期的非税负担冲击:基于“营改增”的研究[J].经济研究,2020,55(11):67-83.）。最近，NBER上的两篇研究中国经济社会的文章（Li et al.(2024,w32162)研究高考中引入英语听力对城乡学生的入学不平等的影响；Beraja et al.(2024,w27723)研究了中国安全技术政府采购对创新的影响）都采用了DDD的研究设计。

我从去年下半年开始，就向国内介绍最新的DDD理论进展和新的认识。总结起来：经验应用文献对DDD的传统用发主要是基于安慰剂检验和基准回归。但是大家并没有真正认识到，如果DD的平行趋势满足，DD得到无偏的估计量，DDD估计量并不能为经验研究提供任何额外的有用的信息。更为重要的是，DDD作为基准回归时，并不需要满足DD和安慰剂DD同时满足平行趋势假设，而且即使DD不满足平行趋势时，我们使用DDD可以更大的发挥其优势，同样可以识别出处理效应。

基于这些理论与认知的进展，“DDD研究设计的应用应该如何展开”成为对应用经济学者来说非常关心的问题。

我自己有一些自己的想法，大家可以去B站（博主： 宏观经济研学会 ： 【CIMERS公开课】不平行的DD，平行的DDD ）看看相关的讲座视频。下面，给出视频配套的do文件：

****************************************************

* DDD.do
* 作者：许文立 [email protected] 安徽大学经济学院
* 描述：用潜在结果框架来模拟ddd——无偏的DD，不需要DDD：最低工资调整
* 参考：Scott Cunningham，Series of “Triple Differences”
* 日期：2024-3-18
****************************************************


clear
capture log close
set seed 20240318

* 省，县和时间
set obs 30
gen province = _n

* 设定处理组,30个省中有15个省提高了最低工资
gen treated = 0
replace treated = 1 in 1/15


* 每个省有100个县
expand 100
bysort province: gen county_no = _n
egen county =group(county_no province)
drop county_no

* 每个县有三组人群：低收入（1），高收入（2），退休（3）
expand 3
bysort county province: gen worker = _n
egen id = group(worker county province)

* 时间：10年
expand 10
sort province
bysort province county worker: gen year = _n

* 设定年份
foreach y of numlist 1/10 {
	local year 2014 +`y'-1
	replace year = `year' if year == `y'
}

* 定义处理后时期（2019年及以后）
gen post = year >= 2019

* 基期收入（2014）：不同的人群收入不同
gen baseline = 36000 if worker == 1     // 低收入人群
replace baseline = 60000 if worker == 2   // 高收入人群
replace baseline = 48000 if worker == 3   // 退休人群

* 处理地区的基期收入要更高
replace baseline = 2 * baseline if treated == 1

* Y（0）的年增长率
gen year_diff = year - 2014
gen y0 = baseline +2000 * year_diff

* Y（0）的随机误差项
gen error = rnormal(0,2500)
replace y0 = y0 + error

* 定义处理效应为10000(处理地区的低收入人群)，Y(1)
gen y1 = y0
replace y1 = y0 + 10000 if treated ==1 & worker == 1 & post == 1
gen delta = y1 - y0
su delta if treated == 1 & worker ==1 & post == 1

* 处理变量
gen d = 0
replace d = 1 if treated == 1 & worker ==1 & post == 1

* 最终的收入
gen income = d*y1 + (1-d)*y0


*****************************************
* Case 1：基于无偏DD和空的安慰剂DD的DDD
*****************************************

* 1. 处理后，低收入者，处理地区
egen avg_wage_plt = mean(income) if post == 1 & treated == 1 & worker ==1

* 2. 处理前，低收入者，处理地区
egen avg_wage_blt = mean(income) if post == 0 & treated == 1 & worker ==1

* 3. 处理后，高收入者和退休，处理地区
egen avg_wage_pht = mean(income) if post == 1 & treated == 1 & worker !=1

* 4. 处理前，高收入者和退休，处理地区
egen avg_wage_bht = mean(income) if post == 0 & treated == 1 & worker !=1

* 5. 处理后，低收入者，控制组地区
egen avg_wage_plc = mean(income) if post == 1 & treated == 0 & worker ==1

* 6. 处理前，低收入，控制组地区
egen avg_wage_blc = mean(income) if post == 0 & treated == 0 & worker ==1

* 7. 处理后，高收入者和退休，控制组地区
egen avg_wage_phc = mean(income) if post == 1 & treated == 0 & worker != 1

* 8. 处理前，高收入者和退休，控制组地区
egen avg_wage_bhc = mean(income) if post == 0 & treated == 0 & worker != 1



** 无偏的DID效应为 10000
* 概述和储存各个均值
su avg_wage_plt, meanonly
local post_lowwage_tre = r(mean)

su avg_wage_blt, meanonly
local befor_lowwage_tre = r(mean)

su avg_wage_plc, meanonly
local post_lowwage_con = r(mean)

su avg_wage_blc, meanonly
local before_lowwage_con = r(mean)

* 计算DID
local DID = (`post_lowwage_tre' - `befor_lowwage_tre') - (`post_lowwage_con' - `before_lowwage_con')

* 展示DID估计量
dis "DID估计量：" `DID'        // 有限样本量下，均值计算有一定的误差，10031.195


* 安慰剂DID
su avg_wage_pht,meanonly
local post_control_tre = r(mean)

su avg_wage_bht,meanonly
local before_control_tre = r(mean)

su avg_wage_phc,meanonly
local post_control_con = r(mean)

su avg_wage_bhc,meanonly
local before_control_con = r(mean)

* 计算处理地区的低收入人群和控制组人群（高收入、退休）的DID
local DID_lowwage_tre = (`post_lowwage_tre' - `befor_lowwage_tre')
local DID_control_tre = (`post_control_tre' - `before_control_tre')

* 计算控制组地区的低收入人群和控制组人群（高收入、退休）的DID
local DID_lowwage_con = (`post_lowwage_con' - `before_lowwage_con')
local DID_control_con = (`post_control_con' - `before_control_con')

* 计算DDD
local DDD = (`DID_lowwage_tre' - `DID_control_tre') - (`DID_lowwage_con' - `DID_control_con')

* 展示DDD结果
dis "DDD估计量："`DDD'       // 有限样本量下，均值计算有一定的误差，10036.598


***********************************
* DDD回归
***********************************

gen lowwage = 0
replace lowwage = 1 if worker == 1

reg income post##treated##lowwage,cluster(province)

su delta if treated == 1 & worker == 1 & post == 1

dis `r(mean)'


* DDD事件研究图

reg income i.post##treated##lowwage,cluster(province)

* 事件研究图

gen T = 0 
replace T = 1 if treated == 1 & lowwage == 1

reg income treated## lowwage year##treated year##lowwage T##ib2018.year,cluster(province)


coefplot, keep(1.T#*) omitted baselevels cirecast(rcap) ///
    rename(1.T#([0-9]+).year = \1, regex) at(_coef) ///
    yline(0, lp(solid)) xline(2018.5, lpattern(dash)) ///
    xlab(2014(1)2023)

	
*****************************************************************
* 结论：有上述结果可以看出，DD无偏的（10031），结果非常接近真实的处理效应10000，
* 而DDD的结果是10036，与DD非常接近，此时，由于DD是无偏的，满足平行趋势，
* 所以DDD是在做无用功。
*****************************************************************

capture log close



******************************************************
* 下面，我们来看看基于有偏的DD（不满足平行趋势），DDD的作用



clear
capture log close
set seed 20240319

* 省，县和时间
set obs 30
gen province = _n

* 设定处理组,30个省中有15个省提高了最低工资
gen treated = 0
replace treated = 1 in 1/15


* 每个省有100个县
expand 100
bysort province: gen county_no = _n
egen county =group(county_no province)
drop county_no

* 每个县有三组人群：低收入（1），高收入（2），退休（3）
expand 3
bysort county province: gen worker = _n
egen id = group(worker county province)

* 时间：10年
expand 10
sort province
bysort province county worker: gen year = _n

* 设定年份
foreach y of numlist 1/10 {
	local year 2014 +`y'-1
	replace year = `year' if year == `y'
}

* 定义处理后时期（2019年及以后）
gen post = year >= 2019

* 基期收入（2014）：不同的人群收入不同
gen baseline = 36000 if worker == 1     // 低收入人群
replace baseline = 60000 if worker == 2   // 高收入人群
replace baseline = 48000 if worker == 3   // 退休人群

* 处理地区的基期收入要更高
replace baseline = 2 * baseline if treated == 1

* Y（0）的年增长率
gen year_diff = year - 2014

* 设定省份和人群的差异化趋势
gen prov_trend = 10000 if treated == 1
replace prov_trend = 15000 if treated == 0

gen group_trend = 5000 if worker == 1
replace group_trend = 8000 if worker == 2 | worker == 3

* 收入水平Y0，合并趋势
gen y0 = baseline + prov_trend*year_diff + group_trend*year_diff

* Y（0）的随机误差项
gen error = rnormal(0,2500)
replace y0 = y0 + error

* 定义处理效应为10000(处理地区的低收入人群)，Y(1)
gen y1 = y0
replace y1 = y0 + 10000 if treated ==1 & worker == 1 & post == 1
gen delta = y1 - y0
su delta if treated == 1 & worker ==1 & post == 1

gen ATT = `r(mean)'
su ATT


* 处理变量
gen d = 0
replace d = 1 if treated == 1 & worker ==1 & post == 1

* 最终的收入
gen income = d*y1 + (1-d)*y0


*****************************************
* Case 2：比较处理地区和控制地区的低收入人群：DD
*****************************************

* 1. 处理后，低收入者，处理地区
egen avg_wage_plt = mean(income) if post == 1 & treated == 1 & worker ==1

* 2. 处理前，低收入者，处理地区
egen avg_wage_blt = mean(income) if post == 0 & treated == 1 & worker ==1

* 3. 处理后，高收入者和退休，处理地区
egen avg_wage_pht = mean(income) if post == 1 & treated == 1 & worker !=1

* 4. 处理前，高收入者和退休，处理地区
egen avg_wage_bht = mean(income) if post == 0 & treated == 1 & worker !=1

* 5. 处理后，低收入者，控制组地区
egen avg_wage_plc = mean(income) if post == 1 & treated == 0 & worker ==1

* 6. 处理前，低收入，控制组地区
egen avg_wage_blc = mean(income) if post == 0 & treated == 0 & worker ==1

* 7. 处理后，高收入者和退休，控制组地区
egen avg_wage_phc = mean(income) if post == 1 & treated == 0 & worker != 1

* 8. 处理前，高收入者和退休，控制组地区
egen avg_wage_bhc = mean(income) if post == 0 & treated == 0 & worker != 1



** 有偏的DID，效应为 10000
* 概述和储存各个均值
su avg_wage_plt, meanonly
local post_lowwage_tre = r(mean)

su avg_wage_blt, meanonly
local befor_lowwage_tre = r(mean)

su avg_wage_plc, meanonly
local post_lowwage_con = r(mean)

su avg_wage_blc, meanonly
local before_lowwage_con = r(mean)

* 计算DID
local DID = (`post_lowwage_tre' - `befor_lowwage_tre') - (`post_lowwage_con' - `before_lowwage_con')

* 展示DID估计量
dis "DID估计量：" `DID'        // DD不满足平行趋势，-14969.547，真实处理效应10000

* 回归结果

gen T = 0
replace T = 1 if treated ==1 & worker == 1
replace T =. if worker ~= 1

reg income T##post,cluster(province)       // 相同的结果-14969.55

* 有偏的事件研究
reg income T##ib2018.year, cluster(province)

coefplot, keep(1.T#*) omitted baselevels cirecast(rcap) ///
    rename(1.T#([0-9]+).year = \1, regex) at(_coef) ///
    yline(0, lp(solid)) xline(2018.5, lpattern(dash)) ///
    xlab(2014(1)2023)


	
*****************************************
* Case 3：比较处理地区的低收入人群和高收入人群：DD
*****************************************	
	
* 安慰剂DID
su avg_wage_plt,meanonly
local post_lowwage_tre = r(mean)

su avg_wage_blt,meanonly
local before_lowwage_tre = r(mean)

su avg_wage_pht,meanonly
local post_control_tre = r(mean)

su avg_wage_bht,meanonly
local before_control_tre = r(mean)

* 计算处理地区的低收入人群和控制组人群（高收入、退休）的DID
local DID_tre = (`post_lowwage_tre' - `before_lowwage_tre') - (`post_control_tre' - `before_control_tre')

* 展示DD结果
dis "DID估计量：" `DID_tre'       // -4942.5156

*


* 说明平行趋势偏误
* 计算处理地区的低收入人群的趋势差异
su y0 if treated ==1 & worker == 1 & post ==0
scalar avg_y0_before_lowwage_tre = r(mean)
su y0 if treated ==1 & worker == 1 & post ==1
scalar avg_y0_post_lowwage_tre = r(mean)

scalar l_t = avg_y0_post_lowwage_tre - avg_y0_before_lowwage_tre

* 计算处理地区的高收入人群和退休人群的趋势差异
su y0 if treated ==1 & (worker == 2 |worker ==3) & post ==0
scalar avg_y0_before_lowwage_con = r(mean)
su y0 if treated ==1 & (worker == 2 |worker ==3)  & post ==1
scalar avg_y0_post_lowwage_con = r(mean)

scalar s_t = avg_y0_post_lowwage_con - avg_y0_before_lowwage_con

dis "l_t=" l_t
dis "s_t=" s_t

scalar non_pt =l_t - s_t
dis "非平行趋势偏离=" non_pt        // -14942.519 =-4942.519 - 10000(ATT)


* 回归结果

gen T1 = 0
replace T1 = 1 if treated ==1 & worker == 1
replace T1 =. if treated == 0

reg income T1##post,cluster(province)       // 相同的结果-4942.519

* 有偏的事件研究
reg income T1##ib2018.year, cluster(province)

coefplot, keep(1.T1#*) omitted baselevels cirecast(rcap) ///
    rename(1.T1#([0-9]+).year = \1, regex) at(_coef) ///
    yline(0, lp(solid)) xline(2018.5, lpattern(dash)) ///
    xlab(2014(1)2023)



******************************************************
* 无偏的DDD
******************************************************

su avg_wage_pht,meanonly
local post_control_tre = r(mean)

su avg_wage_bht,meanonly
local before_control_tre = r(mean)

su avg_wage_phc,meanonly
local post_control_con = r(mean)

su avg_wage_bhc,meanonly
local before_control_con = r(mean)

* 计算处理地区的低收入与控制组人群的dd
local DID_lowwage_tre = (`post_lowwage_tre' - `befor_lowwage_tre')
local DID_control_tre = (`post_control_tre' - `before_control_tre')


* 计算控制组地区的低收入与控制组人群的DD
local DID_lowwage_con = (`post_lowwage_con' - `before_lowwage_con')
local DID_control_con = (`post_control_con' - `before_control_con')

* 计算DDD
local DDD = (`DID_lowwage_tre' - `DID_control_tre') - (`DID_lowwage_con' - `DID_control_con')

*展示ddd结果
dis "DDD估计量=" `DDD'



****** DDD回归

gen 	T2 = 0
replace T2 = 1 if worker==1

reg income post##treated##T2, cluster(province)

su delta if treated==1 & worker==1 & post==1
di `r(mean)'

****** DDD事件研究图

reg income i.year##treated##T2, cluster(province)


gen 	T3=0
replace T3=1 if treated==1 & T2==1

reg income treated##T2 year##treated year##T2 T3##ib2018.year, cluster(province)


coefplot, keep(1.T3#*) omitted baselevels cirecast(rcap) ///
    rename(1.T3#([0-9]+).year = \1, regex) at(_coef) ///
    yline(0, lp(solid)) xline(2018.5, lpattern(dash)) ///
    xlab(2014(1)2023)

	
***********************************************************
* 从上述结果可以看出，DD不满足平行趋势
* DDD的结果仍然是无偏的
***********************************************************


***********************************************************
* 将有偏的DD、有偏的安慰剂DD和DDD放在一张事件研究图中
***********************************************************

* 定义临时性文件
tempfile event1 event2 event3


********************************
* 有偏的DD
********************************

* 运行回归
reg income T##ib2018.year,cluster(province)

* 储存回归结果
matrix b=r(table)

* 提取系数和置信区间

local coef00 = b[1, 23] 
local coef01 = b[1, 24]
local coef02 = b[1, 25]
local coef03 = b[1, 26]
local coef04 = b[1, 27]
local coef05 = b[1, 28]
local coef06 = b[1, 29]
local coef07 = b[1, 30]
local coef08 = b[1, 31]
local coef09 = b[1, 32]

local ll00 = b[5, 23]  // Lower limit of confidence interval for coef00
local ll01 = b[5, 24]  // Lower limit of confidence interval for coef01
local ll02 = b[5, 25]  // Lower limit of confidence interval for coef02
local ll03 = b[5, 26]  // Lower limit of confidence interval for coef03
local ll04 = b[5, 27]  // Lower limit of confidence interval for coef04
local ll05 = b[5, 28]  // Lower limit of confidence interval for coef05
local ll06 = b[5, 29]  // Lower limit of confidence interval for coef06
local ll07 = b[5, 30]  // Lower limit of confidence interval for coef07
local ll08 = b[5, 31]  // Lower limit of confidence interval for coef08
local ll09 = b[5, 32]  // Lower limit of confidence interval for coef09

local ul00 = b[6, 23]  // Upper limit of confidence interval for coef00
local ul01 = b[6, 24]  // Upper limit of confidence interval for coef01
local ul02 = b[6, 25]  // Upper limit of confidence interval for coef02
local ul03 = b[6, 26]  // Upper limit of confidence interval for coef03
local ul04 = b[6, 27]  // Upper limit of confidence interval for coef04
local ul05 = b[6, 28]  // Upper limit of confidence interval for coef05
local ul06 = b[6, 29]  // Upper limit of confidence interval for coef06
local ul07 = b[6, 30]  // Upper limit of confidence interval for coef07
local ul08 = b[6, 31]  // Upper limit of confidence interval for coef08
local ul09 = b[6, 32]  // Upper limit of confidence interval for coef09


preserve

clear
set obs 10
gen year = _n + 2013
gen coef = .
gen ll = .
gen ul = .

replace coef = `coef00' in 1
replace coef = `coef01' in 2
replace coef = `coef02' in 3
replace coef = `coef03' in 4
replace coef = `coef04' in 5
replace coef = `coef05' in 6
replace coef = `coef06' in 7
replace coef = `coef07' in 8
replace coef = `coef08' in 9
replace coef = `coef09' in 10


replace ul = `ul00' in 1
replace ul = `ul01' in 2
replace ul = `ul02' in 3
replace ul = `ul03' in 4
replace ul = `ul04' in 5
replace ul = `ul05' in 6
replace ul = `ul06' in 7
replace ul = `ul07' in 8
replace ul = `ul08' in 9
replace ul = `ul09' in 10

replace ll = `ll00' in 1
replace ll = `ll01' in 2
replace ll = `ll02' in 3
replace ll = `ll03' in 4
replace ll = `ll04' in 5
replace ll = `ll05' in 6
replace ll = `ll06' in 7
replace ll = `ll07' in 8
replace ll = `ll08' in 9
replace ll = `ll09' in 10

* Save this as a tempfile
tempfile event1
save `event1'

* Restore the original dataset
restore



********************************
* 有偏的安慰剂DD
********************************

* 运行回归
reg income T1##ib2018.year,cluster(province)

* 储存回归结果
matrix b=r(table)

* 提取系数和置信区间

local coef00 = b[1, 23] 
local coef01 = b[1, 24]
local coef02 = b[1, 25]
local coef03 = b[1, 26]
local coef04 = b[1, 27]
local coef05 = b[1, 28]
local coef06 = b[1, 29]
local coef07 = b[1, 30]
local coef08 = b[1, 31]
local coef09 = b[1, 32]

local ll00 = b[5, 23]  // Lower limit of confidence interval for coef00
local ll01 = b[5, 24]  // Lower limit of confidence interval for coef01
local ll02 = b[5, 25]  // Lower limit of confidence interval for coef02
local ll03 = b[5, 26]  // Lower limit of confidence interval for coef03
local ll04 = b[5, 27]  // Lower limit of confidence interval for coef04
local ll05 = b[5, 28]  // Lower limit of confidence interval for coef05
local ll06 = b[5, 29]  // Lower limit of confidence interval for coef06
local ll07 = b[5, 30]  // Lower limit of confidence interval for coef07
local ll08 = b[5, 31]  // Lower limit of confidence interval for coef08
local ll09 = b[5, 32]  // Lower limit of confidence interval for coef09

local ul00 = b[6, 23]  // Upper limit of confidence interval for coef00
local ul01 = b[6, 24]  // Upper limit of confidence interval for coef01
local ul02 = b[6, 25]  // Upper limit of confidence interval for coef02
local ul03 = b[6, 26]  // Upper limit of confidence interval for coef03
local ul04 = b[6, 27]  // Upper limit of confidence interval for coef04
local ul05 = b[6, 28]  // Upper limit of confidence interval for coef05
local ul06 = b[6, 29]  // Upper limit of confidence interval for coef06
local ul07 = b[6, 30]  // Upper limit of confidence interval for coef07
local ul08 = b[6, 31]  // Upper limit of confidence interval for coef08
local ul09 = b[6, 32]  // Upper limit of confidence interval for coef09


preserve

clear
set obs 10
gen year = _n + 2013
gen coef2 = .
gen ll2 = .
gen ul2 = .

replace coef2 = `coef00' in 1
replace coef2 = `coef01' in 2
replace coef2 = `coef02' in 3
replace coef2 = `coef03' in 4
replace coef2 = `coef04' in 5
replace coef2 = `coef05' in 6
replace coef2 = `coef06' in 7
replace coef2 = `coef07' in 8
replace coef2 = `coef08' in 9
replace coef2 = `coef09' in 10


replace ul2 = `ul00' in 1
replace ul2 = `ul01' in 2
replace ul2 = `ul02' in 3
replace ul2 = `ul03' in 4
replace ul2 = `ul04' in 5
replace ul2 = `ul05' in 6
replace ul2 = `ul06' in 7
replace ul2 = `ul07' in 8
replace ul2 = `ul08' in 9
replace ul2 = `ul09' in 10

replace ll2 = `ll00' in 1
replace ll2 = `ll01' in 2
replace ll2 = `ll02' in 3
replace ll2 = `ll03' in 4
replace ll2 = `ll04' in 5
replace ll2 = `ll05' in 6
replace ll2 = `ll06' in 7
replace ll2 = `ll07' in 8
replace ll2 = `ll08' in 9
replace ll2 = `ll09' in 10

* Save this as a tempfile
tempfile event2
save `event2'

* Restore the original dataset
restore

********************************
* 无偏的DDD
********************************

* 运行回归
reg income treated##T2 year##treated year##T2 T3##ib2018.year,cluster(province)

* 储存回归结果
matrix b=r(table)

* 提取系数和置信区间

* Extract coefficients and their confidence intervals
local coef00 = b[1, 71] 
local coef01 = b[1, 72]
local coef02 = b[1, 73]
local coef03 = b[1, 74]
local coef04 = b[1, 75]
local coef05 = b[1, 76]
local coef06 = b[1, 77]
local coef07 = b[1, 78]
local coef08 = b[1, 79]
local coef09 = b[1, 80]

local ll00 = b[5, 71]  // Lower limit of confidence interval for coef00
local ll01 = b[5, 72]  // Lower limit of confidence interval for coef01
local ll02 = b[5, 73]  // Lower limit of confidence interval for coef02
local ll03 = b[5, 74]  // Lower limit of confidence interval for coef03
local ll04 = b[5, 75]  // Lower limit of confidence interval for coef04
local ll05 = b[5, 76]  // Lower limit of confidence interval for coef05
local ll06 = b[5, 77]  // Lower limit of confidence interval for coef06
local ll07 = b[5, 78]  // Lower limit of confidence interval for coef07
local ll08 = b[5, 79]  // Lower limit of confidence interval for coef08
local ll09 = b[5, 80]  // Lower limit of confidence interval for coef09

local ul00 = b[6, 71]  // Upper limit of confidence interval for coef00
local ul01 = b[6, 72]  // Upper limit of confidence interval for coef01
local ul02 = b[6, 73]  // Upper limit of confidence interval for coef02
local ul03 = b[6, 74]  // Upper limit of confidence interval for coef03
local ul04 = b[6, 75]  // Upper limit of confidence interval for coef04
local ul05 = b[6, 76]  // Upper limit of confidence interval for coef05
local ul06 = b[6, 77]  // Upper limit of confidence interval for coef06
local ul07 = b[6, 78]  // Upper limit of confidence interval for coef07
local ul08 = b[6, 79]  // Upper limit of confidence interval for coef08
local ul09 = b[6, 80]  // Upper limit of confidence interval for coef09

* Display some of the coefficients and confidence intervals to verify
display "Coefficient for coef00: " `coef00'
display "Lower limit for coef00: " `ll00'
display "Upper limit for coef00: " `ul00'



preserve

clear
set obs 10
gen year = _n + 2013
gen coef3 = .
gen ll3 = .
gen ul3 = .

replace coef3 = `coef00' in 1
replace coef3 = `coef01' in 2
replace coef3 = `coef02' in 3
replace coef3 = `coef03' in 4
replace coef3 = `coef04' in 5
replace coef3 = `coef05' in 6
replace coef3 = `coef06' in 7
replace coef3 = `coef07' in 8
replace coef3 = `coef08' in 9
replace coef3 = `coef09' in 10


replace ul3 = `ul00' in 1
replace ul3 = `ul01' in 2
replace ul3 = `ul02' in 3
replace ul3 = `ul03' in 4
replace ul3 = `ul04' in 5
replace ul3 = `ul05' in 6
replace ul3 = `ul06' in 7
replace ul3 = `ul07' in 8
replace ul3 = `ul08' in 9
replace ul3 = `ul09' in 10

replace ll3 = `ll00' in 1
replace ll3 = `ll01' in 2
replace ll3 = `ll02' in 3
replace ll3 = `ll03' in 4
replace ll3 = `ll04' in 5
replace ll3 = `ll05' in 6
replace ll3 = `ll06' in 7
replace ll3 = `ll07' in 8
replace ll3 = `ll08' in 9
replace ll3 = `ll09' in 10

* Save this as a tempfile
tempfile event3
save `event3'

* Restore the original dataset
restore

clear
use `event1'
merge using `event2'
drop _merge
merge using `event3'
drop _merge

* Event study
twoway (rcap ul ll year, sort lcolor(black) lwidth(thin) lpattern(solid)) ///
       (scatter coef year, sort mcolor(black) msize(small) msymbol(circle)) ///
       (rcap ul2 ll2 year, sort lcolor(red) lwidth(thin) lpattern(solid)) ///
       (scatter coef2 year, sort mcolor(red) msize(small) msymbol(triangle)) ///
       (rcap ul3 ll3 year, sort lcolor(green) lwidth(thin) lpattern(solid)) ///
       (scatter coef3 year, sort mcolor(green) msize(small) msymbol(square)), ///
       yline(0, lwidth(medium) lpattern(solid) lcolor(blue)) ///
       yline(10000, lcolor(gs12) lpattern(dash) lwidth(thin)) ///
       text(10000 2021 "ATT=10000", color(black) size(medium)) ///
       xtitle("Year") ///
       xline(2018.5, lwidth(medium) lpattern(dash) lcolor(blue)) ///
       xlabel(2014(1)2023) ///
       title("两个有偏DD vs. 无偏的DDD") ///
       subtitle("平行偏误假设") ///
       legend(order(2 "Biased DD 1" 4 "Biased DD 2" 6 "Unbiased DDD"))