* this file produces Figures 3-7

// Routine to do lowess regressions with sampling weights. Formulas from Stata manual entry for -lowess-.
mata
mata set matalnum off
mata set mataoptimize on
mata drop wlowess()
// Data must be sorted by [xname]. Variable [yhatname] must exist before calling and will be overwritten.
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

* odbc load, table("Roodman & Morduch HH") dsn(PK) clear
use "Roodman & Morduch HH.dta", clear
keep if wave<4 & pksample

replace pzwflour=0 if pzwflour==. // needed to match sample sizes. These prices are actually stored as 0 in original data.
replace pzmilk=0 if pzmilk==.
xi i.wave i.village
gen double fasset = fnlasset + flandvalb
gen byte ineligible = !eligible
scalar C_0 = ln(1000)
quietly {
	foreach var in prog fprog mprog fbrac mbrac fbrdb mbrdb fgram mgram {
		gen double l`var'lv = ln(`var'lv)
		recode l`var'lv (. = `=C_0')
		gen byte `var'd = `var'lv > 0
	}
	replace landaf = .1 if landaf==0 // according to Pitt e-mail, 1/4/08
	replace landbef = .1 if landbef==0 
	gen double llandaf = ln(landaf)
	gen double llandbef = ln(landbef)
	gen double lpcnsexp = ln(pcnsexp)
	gen byte crcensored = progid<4 & proglv==0
	gen byte nontar = progid==5
}

recode hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd (. = 0)

gen lproglvln1 = cond(lproglv<=ln(1000)+.0001, 0, lproglv)
* Figure 3
scatter lpcnsexp lproglv, msize(vtiny) || scatter lpcnsexp lproglvln1 if lproglvln1==0, msize(vtiny) scheme(s1color) ytitle("") legend(off) ///
	text(5 0 " Censored" " values" " at log 1", place(e) just(left) size(medium)) text(5 `=ln(1000)' "All data, with " "censored values " "at log 1,000 ", place(w) just(right) size(medium)) ///
	xtitle("Log weekly household expenditure/capita (1992 taka)", size(medium)) plotregion(style(none)) yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	title("Log cumulative borrowing" "(1992 taka)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) ylab(, angle(horizontal) labsize(medium)) ///
	xlab(, labsize(medium))

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 llandbef edhead hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd _Iwave* crcensored
local covsm scohead afedhigh amedhigh afadultd         `sexhead' agehead llandbef edhead hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd _Iwave*
local covsf scohead afedhigh amedhigh          amadultd sexhead  agehead llandbef edhead hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd _Iwave* 

* Simple bivariates
reg lfproglv lpcnsexp [pw=weight]
predict OLS_lfproglv_lpcnsexp
reg lmproglv lpcnsexp [pw=weight]
predict OLS_lmproglv_lpcnsexp
reg lfproglv lpcnsexp [pw=weight] if progid!=5
predict OLS_lfproglv_lpcnsexp_tar if progid!=5
reg lmproglv lpcnsexp [pw=weight] if progid!=5
predict OLS_lmproglv_lpcnsexp_tar if progid!=5

reg fprogd lpcnsexp [pw=weight]
predict OLS_fprogd_lpcnsexp
reg mprogd lpcnsexp [pw=weight]
predict OLS_mprogd_lpcnsexp
reg fprogd lpcnsexp [pw=weight] if progid!=5
predict OLS_fprogd_lpcnsexp_tar if progid!=5
reg mprogd lpcnsexp [pw=weight] if progid!=5
predict OLS_mprogd_lpcnsexp_tar if progid!=5

* First-stage projection of credit vars
reg lfproglv `covsf' `villagechars' if pksamplef [pw=weight]
predict lfproglvhat if pksamplef
replace lfproglvhat = C_0 if !pksamplef
reg lmproglv `covsm' `villagechars' if pksamplem [pw=weight]
predict lmproglvhat if pksamplem
replace lmproglvhat = C_0 if !pksamplem
reg lfproglv `covsf' `villagechars' if pksamplef & progid!=5 [pw=weight]
predict lfproglvhat_tar if pksamplef & progid!=5 
replace lfproglvhat_tar = C_0 if !pksamplef
reg lmproglv `covsm' `villagechars' if pksamplem & progid!=5 [pw=weight]
predict lmproglvhat_tar if pksamplem
replace lmproglvhat_tar = C_0 if !pksamplem

probit fprogd `covsf' `villagechars' if pksamplef [pw=weight]
predict fprogdhat if pksamplef
replace fprogdhat = 0 if !pksamplef
probit mprogd `covsm' `villagechars' if pksamplem [pw=weight]
predict mprogdhat if pksamplem
replace mprogdhat = 0 if !pksamplem
probit fprogd `covsf' `villagechars' if pksamplef & progid!=5 [pw=weight]
predict fprogdhat_tar if pksamplef & progid!=5
replace fprogdhat_tar = 0 if !pksamplef
probit mprogd `covsm' `villagechars' if pksamplem & progid!=5 [pw=weight]
predict mprogdhat_tar if pksamplem
replace mprogdhat_tar = 0 if !pksamplef

* Partialling out from projected credit vars
reg lfproglvhat `covsy' `villagechars' [pw=weight]
predict lfproglvhatr, resid
reg lmproglvhat `covsy' `villagechars' [pw=weight]
predict lmproglvhatr, resid
reg lfproglvhat `covsy' `villagechars' if progid!=5 [pw=weight]
predict lfproglvhatr_tar if progid!=5 , resid
reg lmproglvhat `covsy' `villagechars' if progid!=5  [pw=weight]
predict lmproglvhatr_tar if progid!=5 , resid

reg fprogdhat `covsy' `villagechars' [pw=weight]
predict fprogdhatr, resid
reg mprogdhat `covsy' `villagechars' [pw=weight]
predict mprogdhatr, resid
reg fprogdhat `covsy' `villagechars' if progid!=5 [pw=weight]
predict fprogdhatr_tar if progid!=5 , resid
reg mprogdhat `covsy' `villagechars' if progid!=5  [pw=weight]
predict mprogdhatr_tar if progid!=5 , resid

* Partialling out from consumption
reg lpcnsexp `covsy' `villagechars' [pw=weight]
predict lpcnsexpr, resid
reg lpcnsexp `covsy' `villagechars' if progid!=5 [pw=weight]
predict lpcnsexpr_tar if progid!=5 , resid

* OLS of partialled, projected credit on partialled consumption
reg lfproglvhatr lpcnsexpr [pw=weight]
predict OLS_lfproglvhatr_lpcnsexpr
reg lmproglvhatr lpcnsexpr [pw=weight]
predict OLS_lmproglvhatr_lpcnsexpr
reg lfproglvhatr_tar lpcnsexpr_tar [pw=weight] if progid!=5
predict OLS_lfproglvhatr_lpcnsexpr_tar if progid!=5
reg lmproglvhatr_tar lpcnsexpr_tar [pw=weight] if progid!=5
predict OLS_lmproglvhatr_lpcnsexpr_tar if progid!=5

reg fprogdhatr lpcnsexpr [pw=weight]
predict OLS_fprogdhatr_lpcnsexpr
reg mprogdhatr lpcnsexpr [pw=weight]
predict OLS_mprogdhatr_lpcnsexpr
reg fprogdhatr_tar lpcnsexpr_tar [pw=weight] if progid!=5
predict OLS_fprogdhatr_lpcnsexpr_tar if progid!=5
reg mprogdhatr_tar lpcnsexpr_tar [pw=weight] if progid!=5
predict OLS_mprogdhatr_lpcnsexpr_tar if progid!=5

wlowess fprogd lpcnsexp [pw=weight], gen(low_fprogd_lpcnsexp)
wlowess mprogd lpcnsexp [pw=weight], gen(low_mprogd_lpcnsexp)
wlowess fprogd lpcnsexp if progid!=5 [pw=weight], gen(low_fprogd_lpcnsexp_tar)
wlowess mprogd lpcnsexp if progid!=5 [pw=weight], gen(low_mprogd_lpcnsexp_tar)

wlowess lfproglv lpcnsexp [pw=weight], gen(low_lfproglv_lpcnsexp)
wlowess lmproglv lpcnsexp [pw=weight], gen(low_lmproglv_lpcnsexp)
wlowess lfproglv lpcnsexp if progid!=5 [pw=weight], gen(low_lfproglv_lpcnsexp_tar)
wlowess lmproglv lpcnsexp if progid!=5 [pw=weight], gen(low_lmproglv_lpcnsexp_tar)

wlowess fprogdhatr lpcnsexpr [pw=weight], gen(low_fprogdhatr_lpcnsexpr)
wlowess mprogdhatr lpcnsexpr [pw=weight], gen(low_mprogdhatr_lpcnsexpr)
wlowess fprogdhatr_tar lpcnsexpr_tar if progid!=5 [pw=weight], gen(low_fprogdhatr_lpcnsexpr_tar)
wlowess mprogdhatr_tar lpcnsexpr_tar if progid!=5 [pw=weight], gen(low_mprogdhatr_lpcnsexpr_tar)

wlowess lfproglvhatr lpcnsexpr [pw=weight], gen(low_lfproglvhatr_lpcnsexpr)
wlowess lmproglvhatr lpcnsexpr [pw=weight], gen(low_lmproglvhatr_lpcnsexpr)
wlowess lfproglvhatr_tar lpcnsexpr_tar if progid!=5 [pw=weight], gen(low_lfproglvhatr_lpcnsexpr_tar)
wlowess lmproglvhatr_tar lpcnsexpr_tar if progid!=5 [pw=weight], gen(low_lmproglvhatr_lpcnsexpr_tar)

scalar povline = 1.08 * 83.88985/81.48471 * 10.2 * 7 // convert $1.08 from 1993 to 1992 $, then use average 1991-92 taka PPP exchange rate

* Figure 4
scatter low_lfproglv_lpcnsexp_tar lpcnsexp, mcolor(green) msize(vtiny) || scatter OLS_lfproglv_lpcnsexp_tar lpcnsexp, mcolor(green) msize(vtiny) || ///
	scatter low_lmproglv_lpcnsexp_tar lpcnsexp, mcolor(navy) msize(vtiny) || scatter OLS_lmproglv_lpcnsexp_tar lpcnsexp, mcolor(navy) msize(vtiny) ///
	title("Cumulative borrowing" "(1992 taka, fitted values)", orient(horizontal) just(left) place(w) size(medium) span) ///
	text(7.42 4 "Women", place(e) size(medium) color(green)) text(7.17 4.2 "Men", size(medium) color(blue)) sort(lpcnsexp) ///
	scheme(s1color) plotregion(style(none) margin(zero)) legend(off) yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	xtitle("Average weekly household consumption/capita (1992 taka, observed values)", size(medium)) ///
	xlab(`=ln(povline)' "`=string(povline,"%2.0f")'" `=ln(20)' "20" `=ln(30)' "30" `=ln(50)' "50" `=ln(100)' "100" `=ln(200)' "200" `=ln(300)' "300" `=ln(500)' "500" `=ln(1000)' "1,000", labsize(medium))  ///
	ylab(`=ln(1000)' "1,000" `=ln(1200)' "1,200" `=ln(1400)' "1,400" `=ln(1600)' "1,600" `=ln(1800)' "1,800", labsize(medium) angle(horizontal) format(%2.1f)) ///
	text(6.85 `=ln(povline)' "($1/day" "PPP)", color(black) place(n) size(medium) just(center))

* Figure 5
scatter low_lfproglv_lpcnsexp lpcnsexp, mcolor(green) msize(vtiny) || scatter OLS_lfproglv_lpcnsexp lpcnsexp, mcolor(green) msize(vtiny) || ///
	scatter low_lmproglv_lpcnsexp lpcnsexp, mcolor(navy) msize(vtiny) || scatter OLS_lmproglv_lpcnsexp lpcnsexp, mcolor(navy) msize(vtiny) ///
	title("Cumulative borrowing" "(1992 taka, fitted values)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) ///
	text(7.22 4.65 "Women", place(e) size(medium) color(green)) text(7.09 4.2 "Men", size(medium) color(blue)) sort(lpcnsexp) ///
	scheme(s1color) plotregion(style(none) margin(zero)) legend(off) yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	xtitle("Average weekly household consumption/capita (1992 taka, observed values)", size(medium)) ///
	xlab(`=ln(povline)' "`=string(povline,"%2.0f")'" `=ln(20)' "20" `=ln(30)' "30" `=ln(50)' "50" `=ln(100)' "100" `=ln(200)' "200" `=ln(300)' "300" `=ln(500)' "500" `=ln(1000)' "1,000", labsize(medium))  ///
	ylab(`=ln(1000)' "1,000" `=ln(1100)' "1,100" `=ln(1200)' "1,200" `=ln(1300)' "1,300" `=ln(1400)' "1,400" `=ln(1500)' "1,500", labsize(medium) angle(horizontal) format(%2.1f)) ///
	text(`=ln(1000)' `=ln(povline)' "($1/day" "PPP)", color(black) place(n) size(medium) just(center))

* Figure 6
scatter low_lfproglvhatr_lpcnsexpr_tar lpcnsexpr_tar, mcolor(green) msize(vtiny) || scatter OLS_lfproglvhatr_lpcnsexpr_tar lpcnsexpr_tar, mcolor(green) msize(vtiny) || ///
	scatter low_lmproglvhatr_lpcnsexpr_tar lpcnsexpr_tar, mcolor(navy) msize(vtiny) || scatter OLS_lmproglvhatr_lpcnsexpr_tar lpcnsexpr_tar, mcolor(navy) msize(vtiny) ///
	title("Log cumulative borrowing" "(1992 taka, fitted values)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) ///
	text(-.0007 1 "Women", place(e) size(medium) color(green)) text(.002 0 "Men", place(e) size(medium) color(blue)) sort(lpcnsexpr_tar) ///
	scheme(s1color) plotregion(style(none) margin(zero)) legend(off) yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	xtitle("Log weekly household consumption/capita (1992 taka, observed values)", size(medium)) ///
	xlab(-1 `"`=char(150)+"1"'"' 0 0 1 1 2 2, labsize(medium)) ylab(-.01 `"`=char(150)+"0.010"'"' -.005 `"`=char(150)+"0.005"'"' 0 0 .005 "0.005" .01 "0.010", labsize(medium) angle(horizontal))

* Figure 7
scatter low_lfproglvhatr_lpcnsexpr lpcnsexpr, mcolor(green) msize(vtiny) || scatter OLS_lfproglvhatr_lpcnsexpr lpcnsexpr, mcolor(green) msize(vtiny) || ///
	scatter low_lmproglvhatr_lpcnsexpr lpcnsexpr, mcolor(navy) msize(vtiny) || scatter OLS_lmproglvhatr_lpcnsexpr lpcnsexpr, mcolor(navy) msize(vtiny) ///
	title("Log cumulative borrowing" "(1992 taka, fitted values)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) ///
	text(-.06 1 "Women", place(e) size(medium) color(green)) text(.02 1 "Men", place(e) size(medium) color(blue)) sort(lpcnsexpr) ///
	scheme(s1color) plotregion(style(none) margin(zero)) legend(off) yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	xtitle("Log weekly household consumption/capita (1992 taka, observed values)", size(medium)) ///
	xlab(-1 `"`=char(150)+"1"'"' 0 0 1 1 2 2, labsize(medium)) ylab(-.05 `"`=char(150)+"0.05"'"' 0 0 .05 "0.05" .1 "0.10", labsize(medium) angle(horizontal))