* Code for Roodman & Morduch JDS figures 3-5

* global workpath "h:\microfinance\pk\temp"
* cd "$workpath"

set more off
set matsize 800
cmp setup
global indy $cmp_cont
scalar C_t = ln(1000)
scalar C_v = ln(1)

do "PK reproduction 2 routines.do"

PrepHHData
EstimatePKFlip, runname(replication)

// Figure 3: graph of bimodal likelihood, labelled with coefficient on female credit
di in red "Hit Ctrl-break after full-model fit begins and restart code at next line"
EstimatePK, inter
// manually restart code here:
cap drop lnl
gen double lnl = .
postfile bimodal lambda lfgramlvpk ll using bimodal, replace
est restore replication
mat b=e(b)
local ll = e(ll)
est restore replication_flip
mat b_flip=e(b)
local ll_flip=e(ll)
forvalues lambda=-1(.05)2.001 {
	mat b2 = (1-`lambda')*b + `lambda'*b_flip
	mata _first_call = 1
	cmp_lf1 0 b2 lnl
	sum lnl [aw=weightpk], meanonly
	post bimodal (`lambda') (b2[1,3]) (r(sum))
}
postclose bimodal
preserve
use bimodal, clear
est restore replication
mat b=e(b)
local ll = e(ll)
est restore replication_flip
mat b_flip=e(b)
local ll_flip=e(ll)
line ll lfgramlvpk if ll>-21000, scheme(s1color) /*xlabel(-.08(.02).1, labsize(medsmall) format("%3.2f"))*/ ylabel(, labsize(medsmall) angle(horizontal) format("%4.0f")) ///
	title("Log likelihood", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) lwidth(thin) ///
	xtitle("Estimated impact of log cumulative female microcredit from Grameen Bank", margin(medsmall) size(medium)) || ///
	scatteri `ll' `=b[1,3]' `ll_flip' `=b_flip[1,3]', msymbol(Oh) legend(off)  plotregion(style(none) margin(zero)) xsize(6.5) ///
	text(`ll' `=b[1,3]' `"(`=string(`=b[1,3]',"%4.3f")', `=string(`ll', "%4.0f")')"', place(s) margin(medium)) ///
	text(`ll_flip' `=b_flip[1,3]' `"(`=string(`=b_flip[1,3]',"%4.3f")', `=string(`ll_flip', "%4.0f")')"', place(s) margin(medium)) ///
	yscale(lpattern(blank))
restore

* histogram of log household per-capita consumption (dropped from paper)
preserve
keep weightpk lpcnsexppk
gen w100 =weightpk*100
expand w100
histogram lpcnsexppk, fraction scheme(s1color) xtitle("Log per-capita household consumption") ytitle("") plotregion(style(none) margin(zero)) xsize(6.5) ///
	yscale(lpattern(blank)) xscale(lpattern(blank)) title("Fraction", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) ///
	ylabel(, labsize(medsmall) angle(horizontal) format("%3.2f")) xlabel(, labsize(medsmall))
restore

* Figure 5: sensitivity to removing high-consumption obs
preserve
EstimatePKFlip, runname(replication) nooutreg
est restore replication
sort lpcnsexppk
cap postclose dropoutliers
postfile dropoutliers max b1 se1 b2 se2 skew kurt using dropoutliers, replace
tempname b Utest
local b1  =  _b[lfgramlvpk1]
local se1 = _se[lfgramlvpk1]
est restore A_flip	
local b2  =  _b[lfgramlvpk1]
local se2 = _se[lfgramlvpk1]
sum lpcnsexppk, detail
post dropoutliers (lpcnsexppk[_N]) (`b1') (`se1') (`b2') (`se2') (`r(skewness)') (`r(kurtosis)')
forvalues i=1/50 {
	drop if _n==_N
	est restore A
	mat `b' = e(b)
	EstimatePK, init(`b')
	est store A
	local b1  =  _b[lfgramlvpk1]
	local se1 = _se[lfgramlvpk1]
	est restore A_flip
	mat `b' = e(b)
	EstimatePK, init(`b')
	est store A_flip
	if _b[lfgramlvpk1] > `b1' {
		local b2 = `b1'
		local se2 = `se1'
		local b1  = _b[lfgramlvpk1]
		local se1 = _se[lfgramlvpk1]
	}
	else {
		local b2  =  _b[lfgramlvpk1]
		local se2 = _se[lfgramlvpk1]
	}
	sum lpcnsexppk, detail
	post dropoutliers (lpcnsexppk[_N]) (`b1') (`se1') (`b2') (`se2') (`r(skewness)') (`r(kurtosis)')
}
postclose dropoutliers
use "dropoutliers", clear
gen b1u=b1+1.96*se1
gen b1l=b1-1.96*se1
gen b2u=b2+1.96*se2
gen b2l=b2-1.96*se2
twoway rspike b2u b2l max || rspike b1u b1l max || scatter b1 max if abs(b1-b2)<.0001, mcolor(purple) || ///
	scatter b1 b2 max if abs(b1-b2)>=.0001, mcolor(blue*.5 red*.5) ///
	legend(off) xtitle("Maximum log per-capita household consumption allowed in sample") ///
	ytitle("Impact of log female microcredit from Grameen Bank") ysize(4) xsize(6.5) graphregion(fcolor(white)) ///
	xlabel(,format(%2.1f)) ylabel(,format(%2.1f))
restore
ResetGlobals

* Figure 4: bimodal distribution of estimator: bootstrap the estimator that finds both modes and picks the higher one
est restore replication
bootstrap, cluster(villagepk) reps(1000) seed(0) saving(PKreply_highermode, replace every(1)): EstimateForBootstrapHigherMode, runname(replication)
preserve
use "PKreply_highermode", clear
histogram  lpcnsexppk_b_lfgramlvpk1, fraction kdens scheme(s1color) xtitle("Estimated impact of log cumulative female microcredit from Grameen Bank", margin(medsmall) size(medium)) ///
	title("Fraction", margin(medsmall) orient(horizontal) just(left) place(w) size(medium) span) ytitle("") plotregion(style(none) margin(zero)) xsize(6.5) ///
	yscale(lpattern(blank)) xscale(lpattern(blank)) ///
	ylabel(0(.02).1, labsize(medsmall) angle(horizontal) format("%3.2f")) xlabel(, labsize(medsmall) format("%3.2f"))
di "Number converged:"
count if lpcnsexppk_b_lfgramlvpk1<.
di "Number negative:"
count if lpcnsexppk_b_lfgramlvpk1<0
restore

* display consumption details for 16 highest per-capita HH consumption obs (dropped from paper)
* requires a DSN called PK pointing to the Roodman & Morduch SQL Server database
preserve
#delimit ;
odbc load, clear dsn(PK) exec("
	SELECT     floor([Roodman & Morduch HH].nh/10) as [Household ID],  [Roodman & Morduch HH].wave as [Seasonal survey round], 
				[HH consumption details].Division, [HH consumption details].District, [Roodman & Morduch HH].fproglvpk as [Cumulative borrowing female], 
				[Roodman & Morduch HH].mproglvpk as [Cumulative borrowing male], [HH consumption details].hhitd as Purpose,
				[HH consumption details].Spent as Amount, [Roodman & Morduch HH].pcnsexp * [Roodman & Morduch HH].hhnmembers * 16 AS [Total spending], 
				[Roodman & Morduch HH].hhnmembers as [Household members], [Roodman & Morduch HH].lpcnsexppk as [Log consumption_capita_week]
		FROM         [HH consumption details] INNER JOIN
								  (SELECT     TOP (16) wave, nh, pcnsexp, hhnmembers, lpcnsexppk, pksample, 
										case when creditwpk=0 then 0 else EXP(creditwpk) end AS fproglvpk, 
										case when creditmpk=0 then 0 else EXP(creditmpk) end AS mproglvpk
									FROM          [Roodman & Morduch HH] AS [Roodman & Morduch HH_1]
									WHERE      (pksample = 1)
									ORDER BY lpcnsexppk DESC) AS [Roodman & Morduch HH] ON [HH consumption details].nh = [Roodman & Morduch HH].nh AND 
							  [HH consumption details].wave = [Roodman & Morduch HH].wave
		WHERE Spent>=9000
	ORDER BY [Roodman & Morduch HH].lpcnsexppk, [HH consumption details].hhitid
");
#delimit cr
di "HH consumption details for top 16:"
list
restore