1. 断点回归基本原理
断点回归分析被认为是最接近随机实验的检验方法，能够缓解参数估计的内生性问题，近来在越来越多的研究中得到使用。现有资料已经对断点回归方法的基本原理和效应识别进行了较为广泛的介绍，但对阶数选择和稳健性检验等问题的仍相对较少涉及。本文将基于Stata软件来系统介绍断点回归方法的图形观测、效应识别和有效性和稳健性检验。限于篇幅，本文将内容限定于清晰断点回归方法（Sharp Regression Discontinuity Design ），且只考虑只有一个断点和一个分配变量的问题。
2.特征观察
我们先生成一份模拟数据，并保存为 RDD_simu_data0 。生成的数据中， z1 和 z2 为控制变量。 y1 为结果变量（outcome variable）。x 为分配变量（assignment vaiable）。分配点（cutoff point）设定为 0.5 ，从而x大于0.5 的为实验组，小于0.5的为对照组。
此外，在RDD检验中，我们通常还会对分配变量进行去中心化处理，即用分配变量减去分配点值。如本文中，令 xc=x-0.5 。进而 xc 大于 0 的位实验组，反之为对照组。
set obs 4000
set seed 123 //设立种子可重复结果
gen x = runiform() //分配变量
gen xc = x-0.5 //分配变量去中心化
gen e = rnormal()/5 // noise
gen z1 = rnormal()*0.5 //控制变量
gen z2=1+3*invnormal(uniform())+sin(x*5)/3+e //另一个控制变量
gen T=0
replace T=1 if x>0.5 //treatment
gen g0 = 0 + 3*log(x+1) + sin(x*6)/3
gen g1 = T + 3*log(x+1) + sin(x*6)/3
gen y1 = g1 + 0.5*z1 +0.3*z2+ e // outcome vaiable，with cutoff effect
gen y0 = g0 + 0.5*z1 +0.3*z2+ e // outcome variable, without cutoff effect
label var y1 "Outcome variable (y)"
label var y0 "Outcome variable (y)"
label var x "Assignment variable (x)"
label var xc "Centered Assignment variable (x-c)"
label var T "T=1 for x>0.5, T=0 otherwise"
drop e g*
save "RDD_simu_data0.dta", replace //保存一份数据以备后用使用 RDD 方法检验时，首先要确定结果变量在分配点存在跳跃现象，也即存在断点效应。可以用散点图来观察。下图中给出了不存在断点效应和存在断点效应两种情况。如下图所示，右侧图的结果变量y在分配点0.5处一个相对较为明显的跳跃，说明可能存在断点效应。
use "RDD_simu_data0.dta", clear
twoway (scatter y0 xc, msymbol(+) msize(*0.4) mcolor(black*0.3)) , title("无断点")
graph save y0, replace
twoway (scatter y1 xc, msymbol(+) msize(*0.4) mcolor(black*0.3)) , title("有断点")
graph save y1, replace
graph combine y0.gph y1.gph, row(1)
但用散点图来观察存在两个问题：一是样本太多时不够直观，二是实际分析时中跳跃现象可能不那么清晰。为此，我们可以利用拟合方法，对分配点左右分别拟合，通过观察两侧拟合线的的差异来更容易推测跳跃现象是否发生 RDD分析里提供了rdplot命令处理这项工作。下图中分别列出了利用散点图、 rdplot 命令 + 线性拟合、 rdplot命令 + 二阶多项式拟合图和rdplot命令 + 三阶多项式拟合图的结果。本部分相应代码如下。其中 rdplot 命令中， c() 选项表示断点位置，不设定则默认为 0 。 p() 选项表示拟合的阶数。
use "RDD_simu_data0.dta", clear
twoway (scatter y1 xc, msymbol(+) msize(*0.4) mcolor(black*0.3)), title("散点图")
graph save scatter.gph, replace
rdplot y1 xc, c(0) p(1) graph_options(title(线性拟合)) // 线性拟合图
graph save rd1, replace
rdplot y1 xc, c(0) p(2) graph_options(title(二次型拟合))//二次型拟合图
graph save rd2, replace
graph combine scatter.gph rd1.gph rd2.gph
3. 政策效应估计
3.1 局部线性回归
使用局部线性回归法，是假定在断点邻域中的处理效应为线性，通过在左右两侧邻域分别进行线性回归并比较两侧回归系数差异来进行识别。局部回归检验的一个重要环节在于断点邻域大小的选择，也即 RDD 分析里带宽选择 （bandwidth selection） 的权衡问题。这是因为带宽越大，则意味着有越多的样本被纳入检验中，参数估计更准确，但也意味着样本随机性要求越难满足，内生性问题可能更严重。
本文中断点xc的邻域为 ([xc-h1,xc+h2]) ， h1 和 h2 分别为左右两侧带宽。h1和h2可以相等，也可以不等。在断点分析中，可进行局部线性断点回归的命令有 rd、rdrobust 和 rdcv 三个命令。这三个都会自动给出该命令下最优带宽。本部分相应代码如下。
//由于rdc命令回归较为耗时，本文仅随机抽取模拟数据中10%的观察值来演示。
use "RDD_simu_data0.dta", clear
set matsize 2000
set seed 135
sample 10 //随机抽取10%的观察值
rdplot y1 xc, c(0) //检测一下，看看数据特征是否发生明显变化
// 不同局部线性断点回归命令
rd y1 xc, c(0)
rdrobust y1 xc, c(0) p(1)
rdcv y1 xc, thr(0) deg(1)操作结果较多，此处省略。比较三个命令的回归结果，可以发现回归系数分别为 0.982,0.978和0.978 ，不同命令的系数结果基本一致。对于最优带宽选择，三个命令下分别为 0.208，0.187和0.2 ，十分接近。
此外， rd 命令不仅给出了最优带宽，还同时给出了带宽取最优带宽50%和200%的回归结果。观察 rd 命令结果可发现，50%、100%和200%带宽下回归系数也较为接近，分别为 0.954,0.982和0.954 。可知模拟数据的清晰断点回归结果对带宽选择并不敏感，这也表明回归结果是稳健的。
3.2 局部多项式回归
线性假设可能会错误估计了断点左右的回归系数，我们可以采取非线性拟合的办法进行弥补，即使用局部多项式断点回归方法。上文介绍的 rd、rdrobust和rdcv 三个命令，同样可以用于局部多项式断点回归分析。本部分相应代码如下：
在这里插入use "RDD_simu_data0.dta", clear
rdrobust y1 xc //自动选择阶数
rdrobust y1 xc, p(2) //二阶拟合
rdrobust y1 xc, p(3) //三阶拟合代码片结果同样较多，此处省略
对于局部多项式断点回归，关键问题之一在于阶数的选择。我们利用赤池信息准则 （Akaike Information Criterion，AIC） 和贝叶斯信息准则 （Bayesian Information Criterion，BIC） ，选择不同阶数回归中AIC或BIC信息准则小的值。
3.3 全局多项式回归
全局多项式回归是使用了样本里所有数据来进行多项式回归。这在方法上，等价于局部回归分析里将左右带宽设置为分配变量的最小值和最大值。从而可以同样用上述命令来分析。需要注意的是，由于使用了全部样本，全局断点回归分析可能存在较为严重的内生性问题。本部分相应代码如下：
use "RDD_simu_data0.dta", clear
sum xc
local hvalueR=r(max)
local hvalueL= abs(r(min))
rdrobust y1 xc, h(`hvalueL' `hvalueR') //自动选择阶数
rdrobust y1 xc, h(`hvalueL' `hvalueR') p(2) //二阶拟合
rdrobust y1 xc, h(`hvalueL' `hvalueR') p(3) //三阶拟合结果略
4. RDD有效性检验
4.1 局部平滑性的检验
对于局部平滑假设，是指除了结果变量，所有所有其它变量在断点附近都不应该存在处理效应，也即没有出现跳跃现象。在检验方法上，我们可以利用图形直接观察，也可以将每一个协变量作为安慰剂结果变量 (placebo outcomes) ，使用断点回归方法进行检验。
use "RDD_simu_data0.dta", clear
rdplot y1 xc graph_options(title(z1平滑性检验))
graph save rdz1_smooth, replace
rdplot z2 xc graph_options(title(z2平滑性检验))/
graph save rdz2_smooth, replace
graph combine rdz1_smooth.gph rdz2_smooth.gph, title("变量z1 & z2的平滑性检验")
// 从图形，似乎是存在跳跃的，但这并不严格，要看回归结果
rdrobust z1 xc
rdrobust z2 xc
两个回归结果的p值分别 0.399 和 0.741 ，说明不能拒绝不存在断点的假设，可知局部平滑假设满足。
4.2 驱动变量不受人为控制的检验
检验的思路在于，如果不存在人为操控，那么在断点附近样本的数量应该相近，才符合随机性。我们可以用 rddensity 命令来检验断点两侧样本数量是否相近。本部分相应代码如下所示，回归结的 p 值为 0.195 ，不能拒绝断点附近两测样本量大致相等的假设，可知驱动变量不受人为控制的假设满足。
use "RDD_simu_data0.dta", clear
rdrobust y1 xc
local h = e(h_l) //获取最优带宽
rddensity xc, p(1) hl(`h') hr(`h')主要要安装包括rddensity命令的包，否则会报错，可在help里面搜关键词rddensity，选择安装第一个
5. 稳健性检验
5.1 断点的安慰剂检验
稳健性检验的一个自然而然的思路在于选择一个不同于断点的值作为安慰剂断点 (placcebo cutoff points) 。如果断点回归结果变得不显著，则表明断点的真实性。相应代码分别取真实断点两侧 20%、40%、60% 和 80% 样本分位数处作为断点。作为对比，我们也放入了真实断点在图形里。如下图所示，五 个placebo cutoffs 的回归系数都不显著异于0，从而在这些点处不存在处理效应。
use "RDD_simu_data0.dta", clear
sum xc
local xcmax=r(max)
local xcmin= r(min)
forvalues i=1(1)4{
local jr=`xcmax'/(`i'+1)
local jl=`xcmin'/(`i'+1)
rdrobust y1 xc if xc>0,c(`jr')
estimates store jl`i'
rdrobust y1 xc if xc<0,c(`jl')
estimates store jr`i'
}
//加上真实断点的回归结果，作为benchmark结果
rdrobust y1 xc ,c(0)
estimates store jbaseline
//输出图形
local vlist "jl1 jl2 jl3 jl4 jbaseline jr4 jr3 jr2 jr1 "
coefplot `vlist' , yline(0) drop(_cons) vertical
5.2 样本选择的敏感性检验
由于越接近断点的样本，越有动机去人为操控，我们删除最接近断点的样本，来观察回归是否显著（甜甜圈效应, donut hole approach ）。如果仍旧存在，说明即使存在人为操控，断点效应仍旧存在。下面代码里，我们分别删除了断点附近 5%，10%，15%，25% 和 30% 的样本，进行了 6 组稳健性检验。图形给出了回归系数和 95% 的置信区间。可知，在删除 20% 及以下时，回归结果都保持显著。
use "RDD_simu_data0.dta", clear
sum xc
local xcmax=r(max)
forvalues i=1(1)6{
local j=`xcmax'*0.05*`i'
rdrobust y1 xc if abs(xc)>`j'
estimates store obrob`i'
}
//输出图形
local vlist "obrob1 obrob2 obrob3 obrob4 obrob5 obrob6 "
coefplot `vlist' , yline(0) drop(_cons) vertical
5.3 带宽选择的敏感性检验
带宽长度会显著影响回归结果，一个稳健的结果要求对带宽长度不那么敏感。下面代码里，我们先通过rdrobust命令提取最优带宽h，然后分别手动设置带宽为 h 的 25%-400% 倍，看回归结果是否仍旧显著。图形给出了回归系数和95%的置信区间。可知，在最优带宽 25%-400% 范围内，回归结果保持显著，说明结论较为可靠。
use "RDD_simu_data0.dta", clear
rdrobust y1 xc //自动选择最优带宽
local h = e(h_l) //获取最优带宽
forvalues i=1(1)8{
local hrobust=`h'*0.25*`i'
rdrobust y1 xc ,h(`hrobust')
estimates store hrob`i'
}
//输出图形
local vlist "hrob1 hrob2 hrob3 hrob4 hrob5 hrob6 hrob7 hrob8 "
coefplot `vlist' , yline(0) drop(_cons) vertical
