双重差分法（DID）平行趋势检验绘图的相关操作

DID是近年来的“学术明星”，常用于政策评估领域。DID有一个重要的前提假设——平行趋势假定，规范地对平行趋势假定进行检验是一篇DID论文中至关重要的部分，很多人都不重视这一检验，或者对平行趋势假定有些错误的认识（诸如PSM之后就不用做平行趋势检验）。平行趋势检验的一般做法借鉴的是事件研究法的思想，我们首先需要生成年份虚拟变量与处理组虚拟变量的交互项，将这些交互项作为解释变量进行回归（特别注意要丢掉一期，作为基准组）。交互项的系数反映的就是特定年份处理组和控制组之间的差异（与基准组相比），我们特别希望看到的就是政策时点前的虚拟变量与处理组虚拟变量的交互项的系数不显著，这样才说明在政策时点前处理组和控制组不存在异质性的时间趋势。关于这一部分的更详细内容大家可以去看“中国式DID：不做平行趋势检验的DID不是好DID”和“DID大法：如何用Stata做平行趋势检验”这两篇推文。当然，直接看回归系数可能不够直观，所以我们一般都会根据回归结果绘制出回归系数的取值和置信区间，对政策在不同年份的动态经济效应进行呈现，今天想跟大家分享的就是DID平行趋势检验绘图的Stata操作。

数据来源

Nathan Nunn和Nancy Qian两位大神（2011）关于土豆对旧大陆人口增长和城市化进程贡献的QJE论文是一篇相当经典的DID论文，作者在其个人主页上(https://scholar.harvard.edu/nunn/home)公布了这篇论文的数据和代码，接下来我就使用作者提供的部分数据和代码给大家分享一下DID平行趋势检验绘图的Stata操作。

Replication Data for: Nathan Nunn, Nancy Qian. The potato's contribution to population and urbanization: evidence from a historical experiment[J]. The Quarterly Journal of Economics, 2011, 126(2):593-650.

建议大家在看下面这部分内容之前，最好先看一下“土豆的魔力：土豆对旧大陆人口增长和城市化进程的贡献”这篇推文，这样可能理解起来更加顺畅。

准备工作

首先，我们需要生成年份虚拟变量和处理组变量（在土豆这篇文章中是各国适合种植土豆的土地总面积的自然对数ln_wpot）的交互项，因为一共包含1000、1100、1200、1300、1400、1500、1600、1700、1750、1800、1850和1900年12期数据，所以一共会生成12个交互项。

foreach x in 1000 1100 1200 1300 1400 1500 1600 1700 1750 1800 1850 1900{
 gen ln_wpot_`x'=ln_wpot*[year==`x']
 }
drop ln_wpot_1000

接下来，我们需要使用被解释变量（城市化率city_pop_share）对这些交互项进行回归，我们可以使用reghdfe命令估计出各个交互项的系数（实际上就是一个动态DID模型的估计），特别注意的是，我们在回归时需要丢掉一期作为基准组，否则就会有多重共线性的问题。在这里，我选择的是第1期（1000年那一期）。关于基准组的选择，我个人比较推荐的是选择第1期或者-1期（政策时点前1期）。

. reghdfe city_pop_share ln_wpot_1*, absorb(year isocode c.ln_oworld#year c.ln_tropical#year c.ln_rugged#year c.ln_elevation#year) cluster(isocode)
(MWFE estimator converged in 14 iterations)

HDFE Linear regression                            Number of obs   =      1,552
Absorbing 6 HDFE groups                           F(  11,    129) =       3.68
Statistics robust to heteroskedasticity           Prob > F        =     0.0001
                                                  R-squared       =     0.4555
                                                  Adj R-squared   =     0.3749
                                                  Within R-sq.    =     0.0326
Number of clusters (isocode) =        130         Root MSE        =     0.0385

                              (Std. Err. adjusted for 130 clusters in isocode)
------------------------------------------------------------------------------
             |               Robust
city_pop_s~e |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
ln_wpot_1100 |  -.0006226   .0011987    -0.52   0.604    -.0029942    .0017489
ln_wpot_1200 |  -.0011701    .001165    -1.00   0.317    -.0034752    .0011349
ln_wpot_1300 |   .0013089   .0013653     0.96   0.340    -.0013924    .0040103
ln_wpot_1400 |   .0009813   .0013815     0.71   0.479     -.001752    .0037146
ln_wpot_1500 |   .0007654    .001226     0.62   0.534    -.0016602     .003191
ln_wpot_1600 |  -.0000298   .0027418    -0.01   0.991    -.0054546     .005395
ln_wpot_1700 |   .0022237   .0014821     1.50   0.136    -.0007087     .005156
ln_wpot_1750 |    .001334   .0016905     0.79   0.431    -.0020106    .0046787
ln_wpot_1800 |   .0017523   .0016173     1.08   0.281    -.0014476    .0049521
ln_wpot_1850 |   .0031259    .001659     1.88   0.062    -.0001565    .0064084
ln_wpot_1900 |   .0100279   .0030556     3.28   0.001     .0039824    .0160735
       _cons |   .0126719   .0040811     3.11   0.002     .0045974    .0207464
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-------------------------------------------------------------+
         Absorbed FE | Categories  - Redundant  = Num. Coefs |
---------------------+---------------------------------------|
                year |        12           0          12     |
             isocode |       130         130           0    *|
    year#c.ln_oworld |        12           0          12    ?|
  year#c.ln_tropical |        12           0          12    ?|
    year#c.ln_rugged |        12           0          12    ?|
 year#c.ln_elevation |        12           0          12    ?|
-------------------------------------------------------------+
? = number of redundant parameters may be higher
* = FE nested within cluster; treated as redundant for DoF computation

“自行车”：connect和line命令绘图

估计出回归系数后，我们需要做的就是在图中绘制出回归系数的取值情况和置信区间。第一种绘图方法就是使用connect和line等绘图命令进行绘图。这种方法说难不难，就是较为繁琐。首先，我们需要导出回归结果（主要就是需要回归系数和标准误），然后对回归结果进行简单整理，去除如Observations等没有用的信息。

outreg2 using "urbanization_figure.txt", replace sideway noparen se nonotes nocons noaster nolabel text keep(ln_wpot_1*)
insheet using "urbanization_figure.txt", clear
keep if inrange(_n,5,15)
gen year = substr(v1,9,4)
rename (v2 v3)(coef se)
destring, force replace
keep year coef se

接下来，我们需要计算出各个交互项系数的置信区间（t分布95%置信度的临界值可以利用invttail命令计算出来，e(df_r)是存储的自由度）。最后，我们就可以在图中绘制出回归系数的取值和置信区间了，在这里，回归系数使用connect命令进行绘制，95%置信区间使用line命令进行绘制，最后再将回归系数的图和置信区间的图拼在一张图上即可。

gen lb = coef - 1.96*se
gen ub = coef + 1.96*se

twoway (connect coef year,color(gs1) msize(small)) ///
(line lb year,lwidth(thin) lpattern(dash) lcolor(gs2)) ///
(line ub year,lwidth(thin) lpattern(dash) lcolor(gs2)), ///
yline(0,lwidth(vthin) lpattern(dash) lcolor(teal)) ///
xline(1700,lwidth(vthin) lpattern(dash) lcolor(teal)) ///
ylabel(,labsize(*0.85) angle(0)) xlabel(1100(100)1900,labsize(*0.75)) ///
ytitle("Coefficients") ///
xtitle("Year") ///
legend(off) ///图例
graphregion(color(white)) //白底

从图中可以看出，在1750年之前，交互项的系数基本在0左右（95%的置信区间包含了0值），这表明在土豆传入之前，土豆种植适宜性不同的地区其城市化进程并没有出现异质性的时间趋势（土豆种植适宜性与城市化率之间的关系随时间的推移是恒定的），这一点支持了我们的平行趋势假定。在1750年之后，土豆的传入对旧大陆城市化进程的影响在数量级上开始逐渐增加，并且在1850年后快速增加。这一点表明土豆对城市化的影响具有延迟效应，作者在文中给出了解释：在马尔萨斯时代（19世纪前），土豆的好处更多地表现为人口的增加，而不是人均收入和城市化水平的提高。
备注：这一部分代码参考了陈祎、范子英、顾晓敏和周黎安（2020，AER）《Arrival of Young Talent: The Send- Down Movement and Rural Education in China》一文公布的部分代码。

connect和line命令绘图（矩阵）

第一种方法使用的是outreg2命令导出回归结果，其实我们还可以使用矩阵导出回归结果。在回归后，我们可以通过ereturn list命令查看Stata储存的回归结果，包括标量、宏和矩阵等，进而导出回归系数和标准误信息。

matrix coef = e(b) //系数
matrix list coef
matrix cov = e(V) //协方差矩阵
matrix list cov 
gen coef = .
gen se = .
forvalues i = 1(1)11 {
 replace coef = coef[1,`i'] if _n==`i'
 replace se = sqrt(cov[`i',`i']) if _n==`i'
}
gen lb=coef-invttail(e(df_r),0.025)*se //置信区间下界
gen ub=coef+invttail(e(df_r),0.025)*se //置信区间上界
keep coef se lb ub
drop if coef == .

后面的绘图命令与第一种方法是一样的，回归系数使用connect命令进行绘制，95%置信区间使用line命令进行绘制。

gen year = _n
forvalues i =1(1)7 {
 replace year = 1100+(`i'-1)*100 if year==`i'
}
forvalues i =8(1)11 {
 replace year = 1750+(`i'-8)*50 if year==`i'
}
twoway (connect coef year,color(gs1) msize(small)) ///
(line lb year,lwidth(thin) lpattern(dash) lcolor(gs2)) ///
(line ub year,lwidth(thin) lpattern(dash) lcolor(gs2)), ///
yline(0,lwidth(vthin) lpattern(dash) lcolor(teal)) ///
xline(1700,lwidth(vthin) lpattern(dash) lcolor(teal)) ///
ylabel(,labsize(*0.85) angle(0)) xlabel(1100(100)1900,labsize(*0.75)) ///
ytitle("Coefficients") ///
xtitle("Year") ///
legend(off) ///图例
graphregion(color(white)) //白底

“电动车”：coefplot命令绘图
第三种绘图方法是使用coefplot命令绘图，coefplot命令可以便捷地根据回归结果帮助我们绘制回归系数的取值和置信区间，常用于DID平行趋势检验制图。我们不需要手动导出导入回归结果，直接在回归后使用coefplot命令就能进行绘图操作。

coefplot, baselevels ///
keep(ln_wpot*) ///
vertical ///转置图形
coeflabels(ln_wpot_1100=1100 ln_wpot_1200=1200 ///
ln_wpot_1300=1300 ln_wpot_1400=1400 ln_wpot_1500=1500 ///
ln_wpot_1600=1600 ln_wpot_1700=1700 ln_wpot_1750=1750 ///
ln_wpot_1800=1800 ln_wpot_1850=1850 ln_wpot_1900=1900) /// 
yline(0,lwidth(vthin) lpattern(dash) lcolor(teal)) ///
xline(7,lwidth(vthin) lpattern(dash) lcolor(teal)) ///
ylabel(,labsize(*0.85) angle(0)) xlabel(,labsize(*0.75)) ///
ytitle("Coefficients") ///
xtitle("Year") ///
msymbol(O) msize(small) mcolor(gs1) ///plot样式
addplot(line @b @at,lcolor(gs1) lwidth(medthick)) ///增加点之间的连线
ciopts(recast(rline) lwidth(thin) lpattern(dash) lcolor(gs2)) ///置信区间样式
graphregion(color(white)) //白底

coefplot命令有着诸多绘图选项，部分绘图选项的解释如下：

• keep：保留指定的系数，ln_wpot*指代“ln_wpot”开头的变量，我们只需要在图中绘制“ln_wpot”开头的这些交互项的系数和置信区间。
• coeflabels：为系数指定自定义标签，在这里用来修改横坐标。
• msymbol、msize和mcolor：设置点的样式、大小和颜色。
• addplot(line @b @at)：增加点之间的连线。
• ciopts：设置置信区间样式。

对比一下，可以发现三种方法绘制出来的图差别很小，我个人更加推荐coefplot命令，毕竟有了“电动车”，还要什么“自行车”，不过“自行车”也有“自行车”的好处，毕竟更加灵活和自由。
需要本文使用的数据和代码的朋友，请在公众号后台对话框内回复关键词“平行趋势2”。喜欢还请帮忙点个在看哦！