RDD断点回归, 实证完整程序百科全书式的宝典

凡是搞计量经济的，都关注这个号了

邮箱：[e-mail]

所有计量经济圈方法论丛的do文件, 宏微数据库和各种软件都放在社群里.欢迎到因果推断研究小组交流访问.感谢北京大学新结构经济学研究中心的青年AP作为社群群友对计量经济圈的大力支持.

通告：各位计量社群群友，咱们的社群群公告对软件目录进行了大幅度的增添，几乎涵盖了所有你能够想得到的社科软件。若你对某些有需要可以直接提取使用，咱们共同的计量经济圈社群讲究互帮互助互相分享交流。

1.实证研究中用到的135篇文章, 社科学者toolkit

断点回归设计RDD是当前最热门的因果推断计量方法，最主要的原因在于它的透明性和强因果识别性，里面的每一步都可以成功运行出来，完整版do文件和数据的放在计量社群里可以提取使用.

导入数据(在社群里)

gen y = outcome   // 结果变量
gen d = running>0     // 处理变量（0/1种类）

gen v = running   // 分配变量或参考变量
gen vd = v*d    // 交互项

local i=1
forvalues i=2/4 {
gen v`i'=v^`i'
gen v`i'd=v`i'*d
}                // 产生分配变量的三次方、四次方和他们与处理变量的交互项
qui tab year, gen(dyear)   // 如果在面板数据中，想要控制年份可以产生虚拟变量
gen pop2 = pop^2   // 将来用在回归中作为协变量，pop的平方项

————————————————
*************************************
*图形识别，提供三种方式
*************************************

**1.结果变量是不是在断点处跳跃---------

global sizebin 0.2   //根据你的那个running variable选择箱体，这个你自己设定参数  
gen bin=floor(v/$sizebin)
gen midbin=bin*$sizebin+0.5*$sizebin

bys bin: egen mean=mean(y)
reg y d v v2 vd v2d, robust
predict fit    
predict fitsd, stdp
gen upfit=fit+1.645*fitsd  // 产生置信区间的上边界
gen downfit=fit-1.645*fitsd  // 产生置信区间的下边界

preserve            // 第一种方式绘制断点回归图
twoway (rarea upfit downfit v, sort fcolor(gs12) lcolor(gs12)) ///
(line fit v if v<0, sort lcolor(green) lwidth(thick)) ///
(line fit v if v>0, sort lcolor(red) lwidth(thick)) ///
(scatter mean midbin, msize(large) mcolor(black) msymbol(circle_hollow)), ///
ytitle("") xtitle("treatment, X (cutoff: X=0)") xline(0, lcolor(black)) ///
legend(off) xlabel(-1(0.2)1) title("policy implementation")
graph copy all, replace 
restore

cmogram y v,cut(0) scatter lineat(0) qfitci // 第二种方式绘制断点回归图形

rdplot y v, cut(0) nbins(10)  // 第三种方式绘制断点回归图

/**通过图形识别，我们发现在断点处结果变量y发生了跳跃**/

———————————————

********************************
*估计结果，使用三种方式
***********************************

**1. 非参数估计--------------

rdrobust y v,c(0) kernel(uni) bwselect(mserd) all  // 使用rdrobust进行的非参数估计
rdrobust y v, c(0) kernel(tri) bwselect(mserd) all // 这里使用的是triangular密度估计
rdrobust y v, c(0) kernel(epa) bwselect(mserd) all // 这里使用的是epanechnikov密度估计

**2. 非参数估计----------------------

rd y v, mbw(50 100 200) gr z0(0) kernel(tri) // 这个根据最优带宽计算了三个相应带宽，感觉比较方便
rd y v, mbw(50 100 200) gr z0(0) kernel(rec) // 这里使用的是rectangle密度估计

**3. 参数估计:局部线性回归------

rdbwselect y v, c(0) kernel(uni) bwselect(mserd) // 选择最优带宽

preserve
keep if v>= -0.216 & v<= 0.216  // 根据上式的结果显示带宽为0.216
eststo x1: reg y d, robust   // 面板的话选择xtreg，如果是2sls选择xtivregre
eststo x2:reg y d##c.v, robust
eststo x3:reg y d##c.(v v2), robust // 局部线性回归法，选择2阶多项式
eststo x4:reg y d##c.(v v2 v3), robust  // 局部线性回归法，选择3阶多项式
eststo x5:reg y d##c.(v v2 v3 v4), robust  // 局部线性回归法，选择4阶多项式
esttab x1 x2 x3 x4 x5 using y.rtf, star(* .1 ** .05  * .01) nogap nonumber replace ///
se(%5.4f) ar2 aic(%10.4f) bic(%10.4f)  //输出结果到rtf格式
restore

********************************
*稳健性检验
*********************************

**1. 加入协变量后看看回归结果是不是依然显著-----

*1.1 非参估计加入协变量
rd y v, cov(pop pop2) mbw(50 100 200) z0(0) kernel(tri) // 加入协变量pop和pop2

*1.2 参数估计加入协变量

preserve
eststo x11: reg y d pop pop2, robust  // 加入协变量pop和它的平方项
eststo x21:reg y d##c.v pop pop2, robust
eststo x31:reg y d##c.(v v2) pop pop2, robust 
eststo x41:reg y d##c.(v v2 v3)pop pop2, robust  
eststo x51:reg y d##c.(v v2 v3 v4) pop pop2, robust 
esttab x11 x21 x31 x41 x51 using y1.rtf, star(* .1 ** .05  * .01) nogap nonumber replace ///
se(%5.4f) ar2 aic(%10.4f) bic(%10.4f)  //输出加入协变量后的结果到rtf格式
restore

————————————————

**2.检验其中的协变量是不是在断点处连续-------

**2.1 绘制图形检验一下协变量pop是不是连续的
cmogram pop v,cut(0) scatter lineat(0) qfitci // 第二种方式绘制断点回归图形

rdplot pop v, cut(0) nbins(10)  // 第三种方式绘制断点回归图

**2.2 使用估计方法估计出来具体系数看显著不

** 非参数估计-----------------
rdrobust pop v,c(0) kernel(uni) bwselect(mserd) all  // 使用rdrobust进行的非参数估计

** 参数估计:局部线性回归------
rdbwselect pop v, c(0) kernel(uni) bwselect(mserd) // 最优带宽的选择

preserve
keep if v>= -0.175 & v<= 0.175   // 根据上式选择的最优带宽为0.175
eststo xa:reg pop d, robust
eststo xb:reg pop d##c.v, robust  // 用协变量作为伪结果变量，进行断点回归，选择1阶多项式
eststo xb:reg pop d##c.(v v2), robust  // 用协变量作为伪结果变量，进行断点回归，选择2阶多项式
eststo xc:reg pop d##c.(v v2 v3), robust  // 用协变量作为伪结果变量，进行断点回归，选择3阶多项式
eststo xd:reg pop d##c.(v v2 v3 v4), robust  // 用协变量作为伪结果变量，进行断点回归，选择4阶多项式
restore

esttab x11 x21 x31 x41 x51 using m.rtf, star(* .1 ** .05  * .01) nogap nonumber replace ///
se(%5.4f) ar2 aic(%10.4f) bic(%10.4f)  //输出加入协变量后的结果到rtf格式

/**结果显示pop回归方程不是显著的，所以rdd是适用于此的**/

————————————————
**3.Mccracy检验:操纵running variable检验---

net install DCdensity, from("http://www.czxa.top/DCdensity") // 安装McCrary检验命令

*注意：以下这个关于分配变量在断点处跳跃的操纵检验会随着下面的binsize和bandwidth设置而不同的

preserve 
DCdensity v, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) b(0.2) h(0.216) // McCracy test
gen upfhat=fhat+1.645*se_fhat
gen lowfhat=fhat-1.645*se_fhat

twoway (rarea upfhat lowfhat r0 if r0<0, sort fcolor(gs12) lcolor(gs12)) ///
 (rarea upfhat lowfhat r0 if r0>0, sort fcolor(gs12) lcolor(gs12)) ///
 (line fhat r0 if r0<0, lcolor(red)) (line fhat r0 if r0>0, lcolor(blue)) ///
 (scatter Yj Xj if Yj>0, mcolor(gs4) msymbol(circle_hollow)),  ///
 ytitle("Density") xtitle("") xline(0) legend(off)

restore

gen t= .079111002/.143889525  // 产生t值，这个需要你根据系数提取出来

display 2*ttail(2651, t)  // 得到p值，2651是自由度

/**可以看出在5%显著性水平下实际上Mccrary检验是通不过的，证明没有操纵**/

** 把邻近断点处的那些密度分布放大一些看，这样可以更能清楚地看见是不是有操纵—-
preserve

DCdensity v, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) b(0.2) h(0.216) // McCracy test
local breakpoint 0
local cellmpname Xj
local cellvalname Yj
local evalname r0
local cellsmname fhat
local cellsmsename se_fhat
drop if `cellmpname' < -1 | `cellmpname' > 0.5  // 把小于-1和大于0.5的部分都去掉
drop if `evalname' < -1 | `evalname' > 0.5
tempvar hi
quietly gen `hi' = `cellsmname' + 1.96*`cellsmsename'
tempvar lo
quietly gen `lo' = `cellsmname' - 1.96*`cellsmsename'
gr twoway (scatter `cellvalname' `cellmpname', msymbol(circle_hollow) mcolor(gray))             ///
(line `cellsmname' `evalname' if `evalname' < `breakpoint', lcolor(black) lwidth(medthick))   ///
(line `cellsmname' `evalname' if `evalname' > `breakpoint', lcolor(black) lwidth(medthick))   ///
(line `hi' `evalname' if `evalname' < `breakpoint', lcolor(black) lwidth(vthin))              ///
(line `lo' `evalname' if `evalname' < `breakpoint', lcolor(black) lwidth(vthin))              ///
(line `hi' `evalname' if `evalname' > `breakpoint', lcolor(black) lwidth(vthin))              ///
(line `lo' `evalname' if `evalname' > `breakpoint', lcolor(black) lwidth(vthin)),             ///
 xline(`breakpoint', lcolor(black)) legend(off)

restore

——————————————

** 4.安慰剂检验-----------------------

**4.1 改变断点的位置-----------------

rdplot y v if v<0, c(-0.25)  // 将原来的断点0改变为新的断点-0.25

rdplot y v if v>0, c(0.25)  // 将原来的断点0改变为新的断点0.25

rdrobust y v,c(-0.25) kernel(uni) bwselect(mserd) all  // 新断点处使用rdrobust进行的非参数估计
rdrobust y v,c(0.25) kernel(uni) bwselect(mserd) all  // 新断点处使用rdrobust进行的非参数估计

/** 通过以上发现改变断点后不显著了，所以我们的断点选择是有道理的**/

**4.2 改变带宽-----------------

rdrobust y v,c(0) kernel(uni) h(0.1) all  // 改变带宽为0.1
rdrobust y v,c(0) kernel(uni) h(0.4) all  // 改变带宽为0.4

/** 通过以上发现改变带宽并没有影响其显著性，因此我们识别的因果效应很稳健**/

完整版do文件和数据的放在计量社群里可以提取使用.

所有计量经济圈方法论丛的do文件, 社科数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.

可以到计量经济圈社群进一步访问交流各种学术问题，这年头，我们不能强调一个人的英雄主义，需要多多汲取他人的经验教训来让自己少走弯路。

计量经济圈当前有几个阵地，他们分别是如下4个matrix：

①小鹅社群：数据软件书籍等所有资料(最多且更新频繁)，

②微信群：服务于计量经济圈社群群友(最活跃)，

③研究小组：因果推断, 空间计量, 面板数据(最专业)，

④QQ群：2000人大群服务于社群群友(最大)。

计量经济圈是中国计量第一大社区，我们致力于推动中国计量理论和实证技能的提升，圈子以海内外高校研究生和教师为主。计量经济圈绝对六多精神：社科资料最多、社科数据最多、科研牛人最多、海外名校最多、热情互助最多、前沿趋势最多。如果你热爱计量并希望长见识，那欢迎你加入到咱们这个大家庭(戳这里)，要不然你只能去其他那些Open access圈子了。注意：进去之后一定要看小鹅社群“群公告”，不然接收不了群息，也不知道怎么进入咱们独一无二的微信群和QQ群。在规则框架下社群交流讨论无时间限制。

只有进去之后才能够看见这个群公告