* File to produce figures 1-2 in Roodman and Morduch, JDS

use "Pitt data set 2011.dta", clear  // Pitt data set, formerly at http://www.pstc.brown.edu/~mp/RMdo/replication.zip

mata
mata set matalnum off
mata set mataoptimize on
// Routine to do lowess regressions with sampling weights. Variable [xbname] must exist before calling and will be overwritten.
// Ditto for [stdpname] unless it is "" in which case standard error of prediction is not returned
// Formulas from Stata manual entry for -lowess-.
// Data must be sorted by [xname]
mata drop wlowess()
void wlowess(string scalar yname, string scalar xname, string scalar wname, real scalar bw, string scalar xbname, string scalar stdpname, | string scalar selectvarname) {
	real scalar N, k, i, xi, imin, imax, Delta; real matrix data, XX; real colvector beta, Y, X, W, xb, stdp, E
	pragma unset X; pragma unset Y; pragma unset xb; pragma unset data; pragma unset stdp
	st_view(data, ., strlen(stdpname)? (xname, yname, wname, xbname, stdpname): (xname, yname, wname, xbname), selectvarname)
	N = rows(data)
	st_subview(xb, data, ., 4)
	if (strlen(stdpname))
		st_subview(stdp, data, ., 5)
	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 = (XX = cholinv(cross(X, 1, W, X, 1))) * cross(X, 1, W, Y, 0)
			xb[i] = data[i,1] * beta[1] + beta[2]
			if (strlen(stdpname)) {
				E = (Y - X * beta[1]) :- beta[2]
				stdp[i] = sqrt((cross(E, W, E)/sum(W)*N/(N-1)) * ((XX[1,1]*xi + 2*XX[1,2])*xi + XX[2,2])) // stdp = sig2*N/(N-1) * X_i * (X'X)^-2 * X_i' where X_i includes 1
			}
		}
	}
}
end

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

ren round wave

bysort vill: egen fprogvill = max(wchoice) // dummy for villages with female credit groups
by     vill: egen mprogvill = max(mchoice) // ditto for male
gen byte progid = mod(int((vill-10)/80)+1, 4) // 1 = BRAC village, 2 = BRDB, 3 = Grameen, 0 = none
gen byte fonly = fprogvill & !mprogvill
gen byte monly = mprogvill & !fprogvill
gen byte progvill = fprogvill | mprogvill
replace loghland = log(123) if nh==32111 & wave==3 // fix data error. loghland==log(10123) in round 3 and log(123) in rounds 1-2
replace loghland = loghland - ln(100) // convert from decimals to acres
gen hland = exp(loghland)
gen credit = ln(hcrprog)
recode credit (. = 0)

graph set window fontface "Calibri"

* Figure 1: Scatter plot of borrowing vs consumption with ln(1) and ln(1000) censoring
gen creditln1000 = cond(credit, credit, ln(1000))
scatter lexptot creditln1000, msize(vtiny) || scatter lexptot credit if credit==0, msize(vtiny) scheme(s1mono) ytitle("") legend(off) ///
	text(7 0 " Nonborrowers," " censored" " with log 1" " as in PK", place(se) just(left) size(medsmall)) ///
	text(7 `=ln(1000)' "Nonborowers, " "if censored " "with log 1,000 " "instead ", place(sw) just(right) size(medsmall)) ///
	text(7 9 "Borrowers", place(s) just(center) size(medsmall)) ///
	xtitle("Log cumulative borrowing (1992 taka)", size(medium)) plotregion(style(none)) yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	title("Log weekly household" "expenditure/capita" "(1992 taka)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) ylab(, angle(horizontal) labsize(medium)) ///
	xlab(, labsize(medium))


*** Figure 2: Lowess plots of borrowing on land
sort loghland
wlowess creditm loghland if mprogvill & wave==1 & hland<=0.5 [pw=weight], xb(low_creditm_loghland_l) stdp(low_creditm_loghland_l_stdp)
wlowess creditf loghland if fprogvill & wave==1 & hland<=0.5 [pw=weight], xb(low_creditf_loghland_l) stdp(low_creditf_loghland_l_stdp)
wlowess creditm loghland if mprogvill & wave==1 & hland> 0.5 [pw=weight], xb(low_creditm_loghland_u) stdp(low_creditm_loghland_u_stdp)
wlowess creditf loghland if fprogvill & wave==1 & hland> 0.5 [pw=weight], xb(low_creditf_loghland_u) stdp(low_creditf_loghland_u_stdp)

foreach var in low_creditm_loghland_l low_creditf_loghland_l low_creditm_loghland_u low_creditf_loghland_u {
	gen `var'_CIL = `var' - 1.96 * `var'_stdp
	gen `var'_CIU = `var' + 1.96 * `var'_stdp
}
	
* borrowing-land Lowess, women
twoway rarea low_creditf_loghland_u_CIU low_creditf_loghland_u_CIL loghland, fcolor(green) fintensity(inten30) astyle(ci) || ///
	rarea low_creditf_loghland_l_CIU low_creditf_loghland_l_CIL loghland, fcolor(green) fintensity(inten30) astyle(ci) || ///
	line low_creditf_loghland_? loghland, lcolor(green..) ||  ///
	scatter creditf loghland if progvill & wave==1, mcolor(green) msize(tiny..) msymbol(circle..) ///
	xline(`=ln(0.5)', lwidth(vvthin) lpattern(dash)) legend(off) ///
	title("Cumulative borrowing" "(1992 taka)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) scheme(s1color) ///
	plotregion(style(none) margin(zero)) ysize(4) xsize(6.5) ///
	yscale(lpattern(blank)) xscale(lpattern(blank) range(-6.214608 2.060514)) xlabel(none) xtitle("") ///
	ylabel(`=ln(1)' "1" `=ln(10)' "10" `=ln(100)' "100" `=ln(1000)' "1,000" `=ln(10000)' "10,000" , labsize(medium) angle(horizontal)) 
* ppt
twoway rarea low_creditf_loghland_u_CIU low_creditf_loghland_u_CIL loghland, fcolor(green) fintensity(inten30) astyle(ci) || ///
	rarea low_creditf_loghland_l_CIU low_creditf_loghland_l_CIL loghland, fcolor(green) fintensity(inten30) astyle(ci) || ///
	line low_creditf_loghland_? loghland, lcolor(green..) ||  ///
	scatter creditf loghland if progvill & wave==1, mcolor(green) msize(tiny..) msymbol(circle..) ///
	xline(`=ln(0.5)', lwidth(vvthin) lpattern(dash)) legend(off) ///
	title("Cumulative borrowing" "(1992 taka)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) scheme(s1rcolor) ///
	plotregion(style(none) margin(zero)) ysize(7.5) xsize(10) ///
	yscale(lpattern(blank)) xscale(lpattern(blank) range(-6.214608 2.060514)) xlabel(none) xtitle("") ///
	ylabel(`=ln(1)' "1" `=ln(10)' "10" `=ln(100)' "100" `=ln(1000)' "1,000" `=ln(10000)' "10,000" , labsize(medium) angle(horizontal)) ///
	xtitle("Household landholdings before borrowing (acres)") /// // keep last two lines for ppt
	xlabel(`=ln(.001)' "0.001" `=ln(.01)' "0.01" `=ln(.1)' "0.1" `=ln(.5)' "0.5" `=ln(1)' "1" `=ln(10)' "10", labsize(medium))


* borrowing-land Lowess, men
twoway rarea low_creditm_loghland_u_CIU low_creditm_loghland_u_CIL loghland, fcolor(navy) fintensity(inten30) astyle(ci) || ///
	rarea low_creditm_loghland_l_CIU low_creditm_loghland_l_CIL loghland, fcolor(navy) fintensity(inten30) astyle(ci) || ///
	line low_creditm_loghland_? loghland, lcolor(navy..) ||  ///
	scatter creditm loghland if progvill & wave==1, mcolor(navy) msize(tiny..) msymbol(circle..) ///
	xline(`=ln(0.5)', lwidth(vvthin) lpattern(dash)) legend(off) scheme(s1color) ///
	plotregion(style(none) margin(zero)) ysize(4) xsize(6.5) ///
	yscale(lpattern(blank)) xscale(lpattern(blank) range(-6.214608 2.060514)) xlabel(none) xtitle("") ///
	ylabel(`=ln(1)' "1" `=ln(10)' "10" `=ln(100)' "100" `=ln(1000)' "1,000" `=ln(10000)' "10,000" , labsize(medium) angle(horizontal)) ///
	xtitle("Household landholdings before borrowing (acres)", size(medium)) /// 
	xlabel(`=ln(.001)' "0.001" `=ln(.01)' "0.01" `=ln(.1)' "0.1" `=ln(.5)' "0.5" `=ln(1)' "1" `=ln(10)' "10", labsize(medium))
* ppt
twoway rarea low_creditm_loghland_u_CIU low_creditm_loghland_u_CIL loghland, fcolor(navy) fintensity(inten30) astyle(ci) || ///
	rarea low_creditm_loghland_l_CIU low_creditm_loghland_l_CIL loghland, fcolor(navy) fintensity(inten30) astyle(ci) || ///
	line low_creditm_loghland_? loghland, lcolor(navy..) ||  ///
	scatter creditm loghland if progvill & wave==1, mcolor(navy) msize(tiny..) msymbol(circle..) ///
	xline(`=ln(0.5)', lwidth(vvthin) lpattern(dash)) legend(off) scheme(s1rcolor) ///
	plotregion(style(none) margin(zero)) ysize(7.5) xsize(10) ///
	yscale(lpattern(blank)) xscale(lpattern(blank) range(-6.214608 2.060514)) xlabel(none) xtitle("") ///
	ylabel(`=ln(1)' "1" `=ln(10)' "10" `=ln(100)' "100" `=ln(1000)' "1,000" `=ln(10000)' "10,000" , labsize(medium) angle(horizontal)) ///
	xtitle("Household landholdings before borrowing (acres)", size(medium)) /// 
	xlabel(`=ln(.001)' "0.001" `=ln(.01)' "0.01" `=ln(.1)' "0.1" `=ln(.5)' "0.5" `=ln(1)' "1" `=ln(10)' "10", labsize(medium)) ///
	title("Cumulative borrowing" "(1992 taka)", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) 
	
* consumption-land Lowess (not in paper)
bysort nh: egen lexptotm = mean(lexptot)
cap drop low_lexptot_loghland*
wlowess lexptotm loghland if hland<=.5 & progvill          [pw=weight], xb(low_lexptot_loghland_l) stdp(low_lexptot_loghland_l_stdp)
wlowess lexptotm loghland if hland> .5 & progvill          [pw=weight], xb(low_lexptot_loghland_u) stdp(low_lexptot_loghland_u_stdp)
foreach var in low_lexptot_loghland_l low_lexptot_loghland_u {
	gen `var'_CIL = `var' - 1.96 * `var'_stdp
	gen `var'_CIU = `var' + 1.96 * `var'_stdp
}
sort loghland
scatter lexptot loghland if progvill, msize(tiny) mcolor(sienna) || /// 
	rarea low_lexptot_loghland_l_CIL low_lexptot_loghland_l_CIU loghland, fintensity(inten30) astyle(ci) || ///
	rarea low_lexptot_loghland_u_CIL low_lexptot_loghland_u_CIU loghland, fintensity(inten30) astyle(ci) || ///
	line low_lexptot_loghland_? loghland, lcolor(sienna..) ///
	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)) ysize(4) xsize(6.5) ///
	yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	ylabel(`=ln(20)' "20" `=ln(50)' "50" `=ln(100)' "100" `=ln(200)' "200" `=ln(500)' "500" `=ln(1000)' "1,000", labsize(medium) angle(horizontal)) ///
	xlabel(`=ln(.001)' "0.001" `=ln(.01)' "0.01" `=ln(.1)' "0.1" `=ln(.5)' "0.5" `=ln(1)' "1" `=ln(10)' "10", labsize(medium))
* ppt
scatter lexptot loghland if progvill, msize(tiny) mcolor(sienna) || /// 
	rarea low_lexptot_loghland_l_CIL low_lexptot_loghland_l_CIU loghland, fintensity(inten30) astyle(ci) || ///
	rarea low_lexptot_loghland_u_CIL low_lexptot_loghland_u_CIU loghland, fintensity(inten30) astyle(ci) || ///
	line low_lexptot_loghland_? loghland, lcolor(sienna..) ///
	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(s1rcolor) ///
	xtitle("Household landholdings before borrowing (acres)") ///
	plotregion(style(none) margin(zero)) ysize(7.5) xsize(10) ///
	yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	ylabel(`=ln(20)' "20" `=ln(50)' "50" `=ln(100)' "100" `=ln(200)' "200" `=ln(500)' "500" `=ln(1000)' "1,000", labsize(medium) angle(horizontal)) ///
	xlabel(`=ln(.001)' "0.001" `=ln(.01)' "0.01" `=ln(.1)' "0.1" `=ln(.5)' "0.5" `=ln(1)' "1" `=ln(10)' "10", labsize(medium))