其他
倍差法DID详解 (二)：多时点 DID (渐进DID)
///设定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年，36-60为2008年. gen D = 0. gen birth_date = 0
forvalues i = 1/20{ replace D = 1 if id == `i' & time >= 2004 replace birth_date = 2004 if id == `i'}
forvalues i = 21/40{ replace D = 1 if id == `i' & time >= 2006 replace birth_date = 2006 if id == `i'}
forvalues i = 41/60{ replace D = 1 if id == `i' & time >= 2008 replace birth_date = 2008 if id == `i'}
///将基础数据结构保存成dta文件，命名为DID_Basic_Simu.dta,默认保存在当前的 working directory 路径下
save "DID_Basic_Simu_1.dta", replace///设定60个观测值，设定随机数种子///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///Y的生成，使得接受冲击的个体的政策真实效果为10bysort id: gen y0 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal()bysort id: gen y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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(D)
///输出生成的图片，令格式为800*600graph export "article2_1.png",as(png) replace width(800) height(600)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 = 600Absorbing 2 HDFE groups F( 3, 528) = 282958.86Statistics robust to heteroskedasticity Prob > F = 0.0000 R-squared = 0.9996 Adj R-squared = 0.9995 Within R-sq. = 0.9994 Root MSE = 0.9630+-----------------------------------------------------------------------------+ | Robust y| Coef. Std. Err. t P>t [95% Conf. Interval]+-----------------------------------------------------------------------------+ D| 9.884974 .1571734 62.89 0.000 9.576212 10.19374 x1| 4.995274 .0056514 883.90 0.000 4.984172 5.006376 x2| 2.998722 .0083092 360.89 0.000 2.982399 3.015045+-----------------------------------------------------------------------------+ 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 | 9.885*** 9.885*** 9.885*** 9.885*** | (0.157) (0.139) (0.157) (0.157) x1 | 4.995*** 4.995*** 4.995*** 4.995*** | (0.006) (0.006) (0.006) (0.006) x2 | 2.999*** 2.999*** 2.999*** 2.999*** | (0.008) (0.009) (0.008) (0.008)+--------------------------------------------------------------------------------+ N | 600 600 600 600Groups | 60 +--------------------------------------------------------------------------------+ * p<0.05, ** p<0.01, *** p<0.001+--------------------------------------------------------------------------------+///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///Y的生成，使得接受冲击的个体的政策真实效果为10bysort id: gen y0 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal()bysort id: gen y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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 期的时间///然后生成相对时间值的虚拟变量，eventt，并将首期设定为基准对照组gen event = time - birth_datereplace event = -4 if event <= -4tab event, gen(eventt)drop eventt1
xtreg y eventt* x1 x2 i.time, r fe
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*600graph export "article2_2.png",as(png) replace width(800) height(600)///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///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 + 1 ) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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(D)
///输出生成的图片，令格式为800*600graph export "article2_3.png",as(png) replace width(800) height(600)3.2 多种估计结果的呈现
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 = 600Absorbing 2 HDFE groups F( 3, 528) = 49885.71Statistics robust to heteroskedasticity Prob > F = 0.0000 R-squared = 0.9973 Adj R-squared = 0.9970 Within R-sq. = 0.9964 Root MSE = 2.3386
+-----------------------------------------------------------------------------+ | Robust y | Coef. Std. Err. t P>t [95% Conf. Interval]+-----------------------------------------------------------------------------+ D | 2.69066 .3612682 7.45 0.000 1.980961 3.40036 x1| 4.987876 .0138214 360.88 0.000 4.960725 5.015028 x2| 3.026438 .0201725 150.03 0.000 2.98681 3.066066+-----------------------------------------------------------------------------+
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 | 2.691*** 2.691*** 2.691*** 2.691*** | (0.361) (0.180) (0.361) (0.361) x1| 4.988*** 4.988*** 4.988*** 4.988*** | (0.014) (0.013) (0.014) (0.014) 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+--------------------------------------------------------------------------------+///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///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 + 1 ) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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 期的时间///然后生成相对时间值的虚拟变量，eventt，并将首期设定为基准对照组
gen event = time - birth_datereplace event = -4 if event <= -4tab event, gen(eventt)drop eventt1
xtreg y eventt* x1 x2 i.time, r fe
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*600graph export "article2_4.png",as(png) replace width(800) height(600)xtset statefip wrkyrgen bb = 0gen aa = 11gen event = wrkyr - branch_reform
replace event = -10 if event <= -10replace event = 15 if event >= 15
///tab event, gen(eventt)gen event1 = event + 10gen y = log(gini)
qui xtreg y ib10.event1 i.wrkyr, fe r
///生成前十期系数均值xforvalues i = 0/9{ gen b_`i' = _b[`i'.event1]
}
gen avg_coef = (b_0+b_9+b_8+b_7+b_6+b_5+b_4+b_3+b_2+b_1)/10su avg_coefreturn list
coefplot, baselevels /// drop(_cons event1 *.wrkyr ) /// if(@at!=11) ///去除当年的系数点 coeflabels(0.event1 = "-10" /// 1.event1 = "-9" /// 2.event1 = "-8" /// 3.event1 = "-7" /// 4.event1 = "-6" /// 5.event1 = "-5" /// 6.event1 = "-4" /// 7.event1 = "-3" /// 8.event1 = "-2" /// 9.event1 = "-1" /// 10.event1 = "0" /// 11.event1 = "1" /// 12.event1 = "2" /// 13.event1 = "3" /// 14.event1 = "4" /// 15.event1 = "5" /// 16.event1 = "6" /// 17.event1 = "7" /// 18.event1 = "8" /// 19.event1 = "9" /// 20.event1 = "10" /// 21.event1 = "11" /// 22.event1 = "12" /// 23.event1 = "13" /// 24.event1 = "14" /// 25.event1 = "15") ///更改系数的label vertical ///转置图形 yline(0, lwidth(vthin) lpattern(dash) lcolor(teal)) ///加入y=0这条虚线 ylabel(-0.06(0.02)0.06) /// xline(11, lwidth(vthin) lpattern(dash) lcolor(teal)) /// ytitle("Percentage Changes", size(small)) ///加入Y轴标题,大小small xtitle("Years relative to branch deregulation", size(small)) ///加入X轴标题，大小small transform(*=@-r(mean)) ///去除前十期的系数均值，更好看 addplot(line @b @at if @at < 11, lcolor(gs0) || line @b @at if @at>11, lcolor(gs0) lwidth(thin)) ///增加点之间的连线 ciopts(lpattern(dash) recast(rcap) msize(medium)) ///CI为虚线上下封口 msymbol(circle_hollow) ///plot空心格式 scheme(s1mono)
///输出生成的图片，令格式为800*600graph export "article2_6.png",as(png) replace width(800) height(600)
参考资料
Beck, T., Levine, R., & Levkov, A. (2010). Big bad banks? The winners and losers from bank deregulation in the United States. The Journal of Finance, 65(5), 1637-1667. 黄河泉老师的帖子 多期DID：平行趋势检验图示 开学礼包：如何使用双重差分法的交叉项（迄今最全攻略） Stata: reghdfe-多维固定效应
. 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年，36-60为2008年. gen D = 0. gen birth_date = 0
forvalues i = 1/20{ replace D = 1 if id == `i' & time >= 2004 replace birth_date = 2004 if id == `i'}
forvalues i = 21/40{ replace D = 1 if id == `i' & time >= 2006 replace birth_date = 2006 if id == `i'}
forvalues i = 41/60{ replace D = 1 if id == `i' & time >= 2008 replace birth_date = 2008 if id == `i'}
///将基础数据结构保存成dta文件，命名为DID_Basic_Simu.dta,默认保存在当前的 working directory 路径下
save "DID_Basic_Simu_1.dta"，replace
///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///Y的生成，使得接受冲击的个体的政策真实效果为10bysort id: gen y0 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal()bysort id: gen y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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(D)
///输出生成的图片，令格式为800*600graph export "article2_1.png",as(png) replace width(800) height(600)
///保存并输出多个命令的结果///输出生成的图片，令格式为800*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)
///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///Y的生成，使得接受冲击的个体的政策真实效果为10bysort id: gen y0 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal()bysort id: gen y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + 10 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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 期的时间///然后生成相对时间值的虚拟变量，eventt，并将首期设定为基准对照组
gen event = time - birth_datereplace event = -4 if event <= -4tab event, gen(eventt)drop eventt1
xtreg y eventt* x1 x2 i.time, r fe
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*600graph export "article2_2.png",as(png) replace width(800) height(600)
. 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年，36-60为2008年. gen D = 0. gen birth_date = 0
forvalues i = 1/20{ replace D = 1 if id == `i' & time >= 2004 replace birth_date = 2004 if id == `i'}
forvalues i = 21/40{ replace D = 1 if id == `i' & time >= 2006 replace birth_date = 2006 if id == `i'}
forvalues i = 41/60{ replace D = 1 if id == `i' & time >= 2008 replace birth_date = 2008 if id == `i'}
///将基础数据结构保存成dta文件，命名为DID_Basic_Simu.dta,默认保存在当前的 working directory 路径下
save "DID_Basic_Simu_1.dta"，replace
///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///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 + 1 ) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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(D)
///输出生成的图片，令格式为800*600graph export "article2_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)
///调用生成的基础数据文件clearuse "DID_Basic_Simu_1.dta"
///Y的生成，使得接受冲击的个体的政策真实效果为10bysort id: gen y0 = 10 + 5 * x1 + 3 * x2 + T + ind + rnormal()bysort id: gen y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2004 & id >= 1 & id <= 20bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2006 & id >= 21 & id <= 40bysort id: replace y1 = 10 + 5 * x1 + 3 * x2 + T + ind + (time - birth + 1 ) * 3 + rnormal() if time >= 2008 & id >= 41 & id <= 60bysort 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 期的时间///然后生成相对时间值的虚拟变量，eventt，并将首期设定为基准对照组
gen event = time - birth_datereplace event = -4 if event <= -4tab event, gen(eventt)drop eventt1
xtreg y eventt* x1 x2 i.time, r fe
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*600graph export "article2_4.png",as(png) replace width(800) height(600)///本文并没有提供 Beck, Levine & Levkov(2010) 原文的代码和数据集，本部分的运行结果仅做生成图形对比。感兴趣的同学可以去引言中提到的经管之家黄河泉老师的帖子下载相关内容进行操作。
xtset statefip wrkyrgen bb = 0gen aa = 11gen event = wrkyr - branch_reform
replace event = -10 if event <= -10replace event = 15 if event >= 15
///tab event, gen(eventt)gen event1 = event + 10gen y = log(gini)
qui xtreg y ib10.event1 i.wrkyr, fe r
///生成前十期系数均值xforvalues i = 0/9{ gen b_`i' = _b[`i'.event1]
}
gen avg_coef = (b_0+b_9+b_8+b_7+b_6+b_5+b_4+b_3+b_2+b_1)/10su avg_coefreturn list
coefplot, baselevels /// drop(_cons event1 *.wrkyr ) /// if(@at!=11) ///去除当年的系数点 coeflabels(0.event1 = "-10" /// 1.event1 = "-9" /// 2.event1 = "-8" /// 3.event1 = "-7" /// 4.event1 = "-6" /// 5.event1 = "-5" /// 6.event1 = "-4" /// 7.event1 = "-3" /// 8.event1 = "-2" /// 9.event1 = "-1" /// 10.event1 = "0" /// 11.event1 = "1" /// 12.event1 = "2" /// 13.event1 = "3" /// 14.event1 = "4" /// 15.event1 = "5" /// 16.event1 = "6" /// 17.event1 = "7" /// 18.event1 = "8" /// 19.event1 = "9" /// 20.event1 = "10" /// 21.event1 = "11" /// 22.event1 = "12" /// 23.event1 = "13" /// 24.event1 = "14" /// 25.event1 = "15") ///更改系数的label vertical ///转置图形 yline(0, lwidth(vthin) lpattern(dash) lcolor(teal)) ///加入y=0这条虚线 ylabel(-0.06(0.02)0.06) /// xline(11, lwidth(vthin) lpattern(dash) lcolor(teal)) /// ytitle("Percentage Changes", size(small)) ///加入Y轴标题,大小small xtitle("Years relative to branch deregulation", size(small)) ///加入X轴标题，大小small transform(*=@-r(mean)) ///去除前十期的系数均值，更好看 addplot(line @b @at if @at < 11, lcolor(gs0) || line @b @at if @at>11, lcolor(gs0) lwidth(thin)) ///增加点之间的连线 ciopts(lpattern(dash) recast(rcap) msize(medium)) ///CI为虚线上下封口 msymbol(circle_hollow) ///plot空心格式 scheme(s1mono)
///输出生成的图片，令格式为800*600graph export "article2_6.png",as(png) replace width(800) height(600)►一周热文统计计量丨倍差法DID详解 (一)：传统 DID
数据呈现丨Python 手绘风格可视化神包：cutecharts
老姚专栏丨极大似然估计并不难理解
软件应用丨STATA经典资源快递已送达，千万别错过！
数据呈现 | R制图：日历热图，杭州空气质量知多少？统计计量丨部分线性函数型系数面板数据模型的估计与Stata应用数据呈现丨R与Office珠联璧合：如何实现R制图便捷导入PPT？
数据Seminar
这里是大数据、分析技术与学术研究的三叉路口
作者：王昆仑出处：Stata连享会推荐：简华（何年华）编辑：Rye
欢迎扫描👇二维码添加关注