* This file produces Figures 1, 2, and 8 and performs the fuzzy RDD regressions reported in text, in which the bandwidth is varied and the range for 
* key t stats is observed

mata
mata set matalnum off
mata set mataoptimize on
mata drop wlowess()
// Routine to do lowess regressions with sampling weights. Variable [yhatname] must exist before calling and will be overwritten.
// Formulas from Stata manual entry for -lowess-.
// Data must be sorted by [xname]
void wlowess(string scalar yname, string scalar xname, string scalar wname, real scalar bw, string scalar yhatname, | string scalar selectvarname) {
	real scalar N, k, i, xi, imin, imax, Delta; real matrix data, Y, X, W, Yhat; real colvector beta
	pragma unset X; pragma unset Y; pragma unset Yhat
	st_view(data, ., (xname, yname, wname, yhatname), selectvarname)
	N = rows(data)
	st_subview(Yhat, data, ., 4)
	k = trunc((N * bw - .5) / 2)
	for (i=N; i; i--) {
		xi = data[i,1]
		imin = max((1, i-k)); imax = min((i+k, N))
		if (Delta = 1.0001 * max((data[imax,1] - xi, xi - data[imin,1]))) {
			st_subview(X, data, (imin, imax), 1)
			st_subview(Y, data, (imin, imax), 2)
			W = data[|imin, 3 \ imax, 3|] :* (1 :- (abs(X :- xi)/Delta):^3):^3
			beta = cholinv(cross(X, 1, W, X, 1)) * cross(X, 1, W, Y, 0)
			Yhat[i] = data[i,1] * beta[1] + beta[2]
		}
	}
}
end

cap program drop wlowess
program define wlowess, sortpreserve
	syntax varlist(ts) [if] [in] [pw fw aw iw], gen(string) [bw(real 0.8)]
	marksample touse
	if "`weight'" != "" {
		tempvar wvar
		qui gen double `wvar' `exp' if `touse'
	}
	else {
		noi di as err "No weights."
		error 1
	}
	qui gen double `gen' = .
	tsrevar `varlist'
	tokenize `r(varlist)'
	sort `2'
	mata wlowess("`1'", "`2'", "`wvar'", `bw', "`gen'", "`touse'")
end

scalar C_0 = ln(1000)   // credit censoring level

* odbc load, exec("select * from [Roodman & Morduch HH] where wave in (1,2,3)") dsn(PK) clear
use "Roodman & Morduch HH.dta", clear
keep if wave<4 & pksample

replace landbef = landbef/100
quietly {
	foreach var in fproglv mproglv fbraclv mbraclv fbrdblv mbrdblv fgramlv mgramlv proglv {
		gen double l`var' = log(`var')
		recode l`var' (. = `=C_0')
	}
	gen double llandbef = log(landbef)
	gen double llandvalb = log(landvalb)
	recode llandvalb (. = `=ln(400)')
	gen double lpcnsexp = log(pcnsexp)
}

preserve
keep if wave==1

* Lowess plots of borrowing dummy on land
gen byte fprogd = (fproglv>0)*100
gen byte mprogd = (mproglv>0)*100
wlowess mprogd llandbef if pksamplem [pw=weight], gen(low_mprogd_llandbef)
wlowess fprogd llandbef if pksamplef [pw=weight], gen(low_fprogd_llandbef)
wlowess lmproglv llandbef if pksamplem [pw=weight], gen(low_lmproglv_llandbef)
wlowess lfproglv llandbef if pksamplef [pw=weight], gen(low_lfproglv_llandbef)
wlowess mprogd llandvalb if pksamplem [pw=weight], gen(low_mprogd_llandvalb)
wlowess fprogd llandvalb if pksamplef [pw=weight], gen(low_fprogd_llandvalb)
wlowess lmproglv llandvalb if pksamplem [pw=weight], gen(low_lmproglv_llandvalb)
wlowess lfproglv llandvalb if pksamplef [pw=weight], gen(low_lfproglv_llandvalb)

* Figure 1
twoway scatter low_mprogd_llandbef llandbef if llandbef>=-4 & llandbef<=2, mcolor(navy) msize(vtiny) || scatter low_fprogd_llandbef llandbef if llandbef>=-4 & llandbef<=2, mcolor(green) msize(vtiny) ///
	xline(`=ln(0.5)', lwidth(vvthin) lpattern(dash)) legend(off) scheme(s1color) sort(llandbef) plotregion(style(none) margin(zero)) ///
	yscale(lpattern(blank)) xscale(lpattern(blank)) ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%", angle(horizontal) labsize(medium)) ///
	xlabel(`=ln(.02)' "0.02" `=ln(.05)' "0.05" `=ln(.1)' "0.1" `=ln(.2)' "0.2" `=ln(.5)' "0.5" `=ln(1)' "1" `=ln(2)' "2" `=ln(5)' "5", labsize(medium)) text(15 -3 "Men", color(navy) size(medium)) text(32 -3 "Women", color(green) size(medium)) ///
	title("Probability of" "borrowing", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) xtitle("Landholdings before borrowing (acres)", size(medium)) ///
	text(0 `=ln(.5)' "Households with" "half acre or less" "formally eligible", place(nw) just(right) size(medium) margin(medsmall))

reg llandvalb llandbef
scalar line = _b[_cons]+ln(0.5)*_b[llandbef]
* Figure 2
twoway scatter low_mprogd_llandvalb llandvalb, mcolor(navy) msize(vtiny) || scatter low_fprogd_llandvalb llandvalb, mcolor(green) msize(vtiny) ///
	legend(off) scheme(s1color) sort(llandvalb) plotregion(style(none) margin(zero)) ///
	yscale(lpattern(blank)) xscale(lpattern(blank)) ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%", angle(horizontal) labsize(medium)) xline(`=line', lwidth(vvthin) lpattern(dash)) ///
	xlabel(`=line' "10.5" `=ln(500)' "0.5" `=ln(1000)' "1" `=ln(2000)' "2" `=ln(5000)' "5" `=ln(10000)' "10" `=ln(20000)' "20" `=ln(100000)' "100" `=ln(200000)' "200" `=ln(500000)' "500", labsize(medium)) ///
	text(15 7 "Men", color(navy) size(medium)) text(32 7 "Women", color(green) size(medium)) ///
	title("Probability of" "borrowing", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) xtitle("Land value before borrowing (1,000 taka)", size(medium)) ///
	text(0 `=line' "Value corresponding" "to half acre from" "full-sample regression" "of log value on log area", place(nw) just(right) size(medium) margin(medsmall))

restore


* Lowess plots of outcome vs. forcing variable.

wlowess lpcnsexp llandbef if pksamplef & landbef<=0.5 & fonly [pw=weight], gen(low_lpcnsexp_llandbef_fl)
wlowess lpcnsexp llandbef if pksamplef & landbef> 0.5 & fonly [pw=weight], gen(low_lpcnsexp_llandbef_fu)
* Figure 8
scatter lpcnsexp llandbef, msize(.5) || line low_lpcnsexp_llandbef_fl llandbef, sort(llandbef) lwidth(vvthin) || line low_lpcnsexp_llandbef_fu llandbef, sort(llandbef) lwidth(vvthin) ///
	|| if llandbef<4  & lpcnsexp<6.608 & pksamplef & fonly, ///
	xline(`=ln(0.5)', lwidth(vvthin) lpattern(dash)) legend(off) ///
	title("Weekly household" "consumption/capita" "(1992 taka)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) scheme(s1color) ///
	xtitle("Household landholdings before borrowing (acres)") ///
	plotregion(style(none) margin(zero)) ///
	yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	ylabel(`=ln(20)' "20" `=ln(30)' "30" `=ln(50)' "50" `=ln(100)' "100" `=ln(200)' "200" `=ln(300)' "300", labsize(medium) angle(horizontal)) ///
	xlabel(`=ln(.005)' "0.005" `=ln(.01)' "0.01" `=ln(.02)' "0.02" `=ln(.05)' "0.05" `=ln(.1)' "0.1" `=ln(.2)' "0.2" `=ln(.5)' "0.5" `=ln(1)' "1" `=ln(2)' "2" `=ln(5)' "5", labsize(medium))


* Fuzzy RDD per Imbens and Lemieux 2008, p. 627, citing Hahn, Todd and Van der Klaauw 2001
* note that 2SLS version is very poor fit 
gen double llandbef_l = (llandbef - ln(0.5)) * (llandbef<=ln(0.5))
gen double llandbef_u = (llandbef - ln(0.5)) * (llandbef>ln(0.5))

scalar bw = 1
ivreg lpcnsexp llandbef_l llandbef_u (lfproglv = r) if pksamplef & !pksamplem & abs(llandbef - ln(0.5)) <= bw [pw=weight], cluster(nh)
ivreg lpcnsexp llandbef_l llandbef_u (lmproglv = r) if pksamplem & !pksamplef & abs(llandbef - ln(0.5)) <= bw [pw=weight], cluster(nh)
ivreg lpcnsexp llandbef_l llandbef_u (lproglv = r) if (pksamplef | pksamplem) & abs(llandbef - ln(0.5)) <= 1.61 [pw=weight], cluster(nh)

xi i.village i.wave
recode hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd (. = 0)
local villagechars primcoed hpruralh hpfampln hpmidwif pzrice pzwflour pzmoil pzhegg pzmilk pzpotato wagef nowagef wagem bamtr
local villagechars _Ivill*
local covsy scohead afedhigh amedhigh afadultd amadultd sexhead agehead edhead hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd _Iwave* crcensored

cap mat drop Rm
cap mat drop Rm_covs
qui forvalues bw=.05(.02)2.001 {
	ivreg lpcnsexp llandbef_l llandbef_u (lproglv = r) if pksamplem & !pksamplef & abs(llandbef - ln(0.5)) <= `bw' [pw=weight], cluster(nh)
	mat Rm = nullmat(Rm) \ `bw', _b[lproglv], (_b[lproglv] - 1.96*_se[lproglv]), (_b[lproglv] + 1.96*_se[lproglv]), _b[lproglv]/_se[lproglv], e(N)
	ivreg lpcnsexp llandbef_l llandbef_u `covsy' `villagechars' (lproglv = r) if pksamplem & !pksamplef & abs(llandbef - ln(0.5)) <= `bw' [pw=weight], cluster(nh)
	mat Rm_covs = nullmat(Rm_covs) \ _b[lproglv], (_b[lproglv] - 1.96*_se[lproglv]), (_b[lproglv] + 1.96*_se[lproglv]), _b[lproglv]/_se[lproglv]
}
mat Rm = Rm, Rm_covs
mat colnames Rm = bandwidth coef "CI bottom" "CI top" t N "with controls:coef" "with controls:CI bottom" "with controls:CI top" "with controls:t"

* Results for male-only sample
mat list Rm
* full male-only sample size:
count if pksamplem & !pksamplef
* check the ranges on the t stats:
mata colmin(st_matrix("Rm")) \ colmax(st_matrix("Rm"))

cap mat drop Rf
cap mat drop Rf_covs
qui forvalues bw=.05(.02)2.001 {
	ivreg lpcnsexp llandbef_l llandbef_u (lproglv = r) if pksamplef & !pksamplem & abs(llandbef - ln(0.5)) <= `bw' [pw=weight], cluster(nh)
	mat Rf = nullmat(Rf) \ `bw', _b[lproglv], (_b[lproglv] - 1.96*_se[lproglv]), (_b[lproglv] + 1.96*_se[lproglv]), _b[lproglv]/_se[lproglv], e(N)
	ivreg lpcnsexp llandbef_l llandbef_u `covsy' `villagechars' (lproglv = r) if pksamplef & !pksamplem & abs(llandbef - ln(0.5)) <= `bw' [pw=weight], cluster(nh)
	mat Rf_covs = nullmat(Rf_covs) \ _b[lproglv], (_b[lproglv] - 1.96*_se[lproglv]), (_b[lproglv] + 1.96*_se[lproglv]), _b[lproglv]/_se[lproglv]
}
mat Rf = Rf, Rf_covs
mat colnames Rf = bandwidth coef "CI bottom" "CI top" t N "with controls:coef" "with controls:CI bottom" "with controls:CI top" "with controls:t"

* Results for female-only sample
mat list Rf
* full female-only sample size:
count if pksamplef & !pksamplem 
* check the ranges on the t stats:
mata colmin(st_matrix("Rf")) \ colmax(st_matrix("Rf"))