I came across a clear introduction to regression discontinuity in chapter 16 of Khandker, Shahidur R., Gayatri B. Koolwal, and Hussain A. Samad. 2010. Handbook on Impact Evaluation: Quantitative Methods and Practices. Washington, D.C: World Bank. Which is available for free at https://documents1.worldbank.org/cur...0Use0Only1.pdf

The book includes (in the appendix) code that I wanted to explore. However, after playing with it, my question isn't how to make this older code work, but how do I do something similar in Stata 15.1?

I'll include the code for fuzzy RD below, but I couldn't get it to work. I got the error statement "insufficient observations to compute bootstrap standard errors." The data file is "hh_98.dta" and it's available at: https://microdata.worldbank.org/index.php/catalog/436 . The code is in the Appendix for chapter 16 starting on page 227 of the book (page 248 of the PDF). You need the -locpoly- user-written command at -net sj 5-2 st0053_2- if you want to try this. Also, if you cut-and-paste from the PDF you'll need to find-and-replace the quotation marks (singles and doubles).

Again, my question isn't how to make this older code work, but how do I do something similar in Stata 15.1?

Thank you.

Code:
 use ... "hh_98.dta", clear

*****Program for Fuzzy Discontinuity;
#delimit ;
gen lexptot=ln(1+exptot);
gen lnland=ln(1+hhland/100);

capture prog drop rd_fuzzy;
prog rd_fuzzy, rclass;
version 8.2;
args treatment outcome;
confirm var `treatment';
confirm var `outcome';
tempname treatrd1 treatrd0 outrd1 outrd0 treat1 treat0 outcome1
outcome0;
locpoly `treatment' lnland if hhland<50, gen(`treatrd1')  at(lnland) nogr tri  w(3) d(1);
locpoly `treatment' lnland if hhland>=50, gen(`treatrd0') at(lnland) nogr tri  w(3) d(1);
locpoly `outcome' lnland if hhland<50, gen(`outrd1') at(lnland) nogr tri  w(3) d(1);
locpoly `outcome' lnland if hhland>=50, gen(`outrd0') at(lnland) nogr tri  w(3) d(1);
sum `treatrd1' if hhland>=45 & hhland<=55, meanonly;
scalar `treat1'=r(mean);
sum `treatrd0' if hhland>=45 & hhland<=55, meanonly;
scalar `treat0'=r(mean);
sum `outrd1' if hhland>=45 & hhland<=55, meanonly;
scalar `outcome1'=r(mean);
sum `outrd0' if hhland>=45 & hhland<=55, meanonly;
scalar `outcome0'=r(mean);
return scalar impact=(`outcome1'-`outcome0')/(`treat1'-`treat0');
end;

***Male participation;
set seed 12345
bootstrap "rd_fuzzy dmmfd lexptot" impact_fuzzy_m=r(impact), reps(100) nowarn
gen t_impact_fuzzy_m=_b[impact_fuzzy_m]/_se[impact_fuzzy_m]
sum t_impact_fuzzy_m
***Female participation;
set seed 123
bootstrap “rd_fuzzy dfmfd lexptot” impact_fuzzy_f=r(impact), reps(100) nowarn
gen t_impact_fuzzy_f=_b[impact_fuzzy_f]/_se[impact_fuzzy_f]
sum t_impact_fuzzy_f