|
斷點(diǎn)回歸設(shè)計(jì)RDD是當(dāng)前最熱門的因果推斷計(jì)量方法,，最主要的原因在于它的透明性和強(qiáng)因果識(shí)別性,，里面的每一步都可以成功運(yùn)行出來(lái),，若需要do文件和數(shù)據(jù)dta的請(qǐng)進(jìn)入計(jì)量經(jīng)濟(jì)圈社群直接提取(文末),。 
gen y = outcome // 結(jié)果變量
gen d = running>0 // 處理變量（0/1種類） gen v = running // 分配變量或參考變量
gen vd = v*d // 交互項(xiàng) local i=1
forvalues i=2/4 {
gen v`i'=v^`i'
gen v`i'd=v`i'*d
} // 產(chǎn)生分配變量的三次方,、四次方和他們與處理變量的交互項(xiàng)
qui tab year, gen(dyear) // 如果在面板數(shù)據(jù)中,，想要控制年份可以產(chǎn)生虛擬變量
gen pop2 = pop^2 // 將來(lái)用在回歸中作為協(xié)變量,，pop的平方項(xiàng) ————————————————
*************************************
*圖形識(shí)別，提供三種方式
************************************* **1.結(jié)果變量是不是在斷點(diǎn)處跳躍---------
global sizebin 0.2 //根據(jù)你的那個(gè)running variable選擇箱體,，這個(gè)你自己設(shè)定參數(shù)
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 // 產(chǎn)生置信區(qū)間的上邊界
gen downfit=fit-1.645*fitsd // 產(chǎn)生置信區(qū)間的下邊界 preserve // 第一種方式繪制斷點(diǎn)回歸圖
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 // 第二種方式繪制斷點(diǎn)回歸圖形 
rdplot y v, cut(0) nbins(10) // 第三種方式繪制斷點(diǎn)回歸圖 
/**通過(guò)圖形識(shí)別,，我們發(fā)現(xiàn)在斷點(diǎn)處結(jié)果變量y發(fā)生了跳躍**/ ——————————————— ********************************
*估計(jì)結(jié)果，使用三種方式
*********************************** **1. 非參數(shù)估計(jì)-------------- rdrobust y v,c(0) kernel(uni) bwselect(mserd) all // 使用rdrobust進(jìn)行的非參數(shù)估計(jì)
rdrobust y v, c(0) kernel(tri) bwselect(mserd) all // 這里使用的是triangular密度估計(jì)
rdrobust y v, c(0) kernel(epa) bwselect(mserd) all // 這里使用的是epanechnikov密度估計(jì)
**2. 非參數(shù)估計(jì)---------------------- rd y v, mbw(50 100 200) gr z0(0) kernel(tri) // 這個(gè)根據(jù)最優(yōu)帶寬計(jì)算了三個(gè)相應(yīng)帶寬,，感覺(jué)比較方便
rd y v, mbw(50 100 200) gr z0(0) kernel(rec) // 這里使用的是rectangle密度估計(jì)
**3. 參數(shù)估計(jì):局部線性回歸------ rdbwselect y v, c(0) kernel(uni) bwselect(mserd) // 選擇最優(yōu)帶寬 preserve
keep if v>= -0.216 & v<= 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階多項(xiàng)式
eststo x4:reg y d##c.(v v2 v3), robust // 局部線性回歸法,，選擇3階多項(xiàng)式
eststo x5:reg y d##c.(v v2 v3 v4), robust // 局部線性回歸法,，選擇4階多項(xiàng)式
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) //輸出結(jié)果到rtf格式
restore
********************************
*穩(wěn)健性檢驗(yàn)
********************************* **1. 加入?yún)f(xié)變量后看看回歸結(jié)果是不是依然顯著----- *1.1 非參估計(jì)加入?yún)f(xié)變量
rd y v, cov(pop pop2) mbw(50 100 200) z0(0) kernel(tri) // 加入?yún)f(xié)變量pop和pop2
*1.2 參數(shù)估計(jì)加入?yún)f(xié)變量 preserve
eststo x11: reg y d pop pop2, robust // 加入?yún)f(xié)變量pop和它的平方項(xiàng)
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) //輸出加入?yún)f(xié)變量后的結(jié)果到rtf格式
restore ———————————————— **2.檢驗(yàn)其中的協(xié)變量是不是在斷點(diǎn)處連續(xù)-------
**2.1 繪制圖形檢驗(yàn)一下協(xié)變量pop是不是連續(xù)的
cmogram pop v,cut(0) scatter lineat(0) qfitci // 第二種方式繪制斷點(diǎn)回歸圖形 rdplot pop v, cut(0) nbins(10) // 第三種方式繪制斷點(diǎn)回歸圖
**2.2 使用估計(jì)方法估計(jì)出來(lái)具體系數(shù)看顯著不 ** 非參數(shù)估計(jì)-----------------
rdrobust pop v,c(0) kernel(uni) bwselect(mserd) all // 使用rdrobust進(jìn)行的非參數(shù)估計(jì)
** 參數(shù)估計(jì):局部線性回歸------
rdbwselect pop v, c(0) kernel(uni) bwselect(mserd) // 最優(yōu)帶寬的選擇 preserve
keep if v>= -0.175 & v<= 0.175 ="">
eststo xa:reg pop d, robust
eststo xb:reg pop d##c.v, robust // 用協(xié)變量作為偽結(jié)果變量，進(jìn)行斷點(diǎn)回歸,，選擇1階多項(xiàng)式
eststo xb:reg pop d##c.(v v2), robust // 用協(xié)變量作為偽結(jié)果變量,，進(jìn)行斷點(diǎn)回歸，選擇2階多項(xiàng)式
eststo xc:reg pop d##c.(v v2 v3), robust // 用協(xié)變量作為偽結(jié)果變量,，進(jìn)行斷點(diǎn)回歸,，選擇3階多項(xiàng)式
eststo xd:reg pop d##c.(v v2 v3 v4), robust // 用協(xié)變量作為偽結(jié)果變量，進(jìn)行斷點(diǎn)回歸,，選擇4階多項(xiàng)式
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) //輸出加入?yún)f(xié)變量后的結(jié)果到rtf格式 /**結(jié)果顯示pop回歸方程不是顯著的,，所以rdd是適用于此的**/ ————————————————
**3.Mccracy檢驗(yàn):操縱running variable檢驗(yàn)--- net install DCdensity, from('http://www./DCdensity') // 安裝McCrary檢驗(yàn)命令
*注意：以下這個(gè)關(guān)于分配變量在斷點(diǎn)處跳躍的操縱檢驗(yàn)會(huì)隨著下面的binsize和bandwidth設(shè)置而不同的
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 // 產(chǎn)生t值，這個(gè)需要你根據(jù)系數(shù)提取出來(lái)
display 2*ttail(2651, t) // 得到p值,，2651是自由度
/**可以看出在5%顯著性水平下實(shí)際上Mccrary檢驗(yàn)是通不過(guò)的,，證明沒(méi)有操縱**/
** 把鄰近斷點(diǎn)處的那些密度分布放大一些看，這樣可以更能清楚地看見(jiàn)是不是有操縱—-
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.安慰劑檢驗(yàn)-----------------------
**4.1 改變斷點(diǎn)的位置----------------- rdplot y v if v<0, c(-0.25) ="">
rdplot y v if v>0, c(0.25) // 將原來(lái)的斷點(diǎn)0改變?yōu)樾碌臄帱c(diǎn)0.25 
rdrobust y v,c(-0.25) kernel(uni) bwselect(mserd) all // 新斷點(diǎn)處使用rdrobust進(jìn)行的非參數(shù)估計(jì)
rdrobust y v,c(0.25) kernel(uni) bwselect(mserd) all // 新斷點(diǎn)處使用rdrobust進(jìn)行的非參數(shù)估計(jì) /** 通過(guò)以上發(fā)現(xiàn)改變斷點(diǎn)后不顯著了,，所以我們的斷點(diǎn)選擇是有道理的**/
**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 /** 通過(guò)以上發(fā)現(xiàn)改變帶寬并沒(méi)有影響其顯著性,，因此我們識(shí)別的因果效應(yīng)很穩(wěn)健**/
|
