# Stata：断点回归RDD简明教程

Stata 连享会   主页 || 视频 || 推文

E-mail: StataChina@163.com

## 2. 图形观察

### 2.1 生成模拟数据

``````	clear all
global dir d:/RDDStata
capture mkdir \$dir
cd \$dir

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  //保存一份数据以备后用

``````

### 2.2 断点效应的图形观察

``````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)
``````

``````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 局部线性回归

``````//由于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)
``````

### 3.2 局部多项式回归

``````use "RDD_simu_data0.dta", clear

rdrobust y1 xc  //自动选择阶数
rdrobust y1 xc, p(2) //二阶拟合
rdrobust y1 xc, p(3) //三阶拟合

``````

``````*-------------myic------------------
program define myic
version 13
qui estat ic
mat a = r(S)
end
*-------------myic------------------
*-Note: 调用自定义程序myic的方法：
* (1) 选中上述代码，
* (2) 按快捷键 Ctrl+R, 将程序读入内存

use "RDD_simu_data0.dta", clear
set matsize 2000
set seed 135
sample 10    //rdcv回归较为耗时，仅随机抽取10%的观察值来演示。

#d ;
rdcv y1 xc, thr(0) deg(1);	myic;   est store m1;
rdcv y1 xc, thr(0) deg(2);	myic;   est store m2;
rdcv y1 xc, thr(0) deg(3);	myic;   est store m3;
#d cr    // #d 表示 #delimit

*-对比回归结果
local m "m1 m2 m3"
esttab `m', mtitle(`m') b(%6.3f) t(%6.3f)  ///
s(N r2 r2_a AIC BIC) nogap compress
``````

AIC 761.160 446.327 819.628
BIC 775.494 464.740 849.097

### 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 局部平滑性的检验

``````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
``````

### 4.2 驱动变量不受人为控制的检验

``````use "RDD_simu_data0.dta", clear

rdrobust y1 xc
local h = e(h_l)   //获取最优带宽
rddensity xc, p(1) hl(`h') hr(`h')
``````

## 5. 稳健性检验

### 5.1 断点的安慰剂检验

``````
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 样本选择的敏感性检验

``````
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 带宽选择的敏感性检验

``````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
``````

## 相关课程

http://lianxh.duanshu.com

### 课程一览

Stata数据清洗 游万海 直播, 2 小时，已上线

Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

#### 关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。
• 连享会-主页知乎专栏，300+ 推文，实证分析不再抓狂。
• 公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：`课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法`

✏ 连享会学习群-常见问题解答汇总：
https://gitee.com/arlionn/WD