* Routines used in other code files for Roodman & Morduch, JDS

* core linear LIML estimator
cap program drop EstimateLIML
program define EstimateLIML
	syntax [if]
	ivreg2 $depvar $covsy $nontarvar _Ivill* ($creditvars = $insts) [pw=weightpk] `if', cluster(nh) liml partial($covsy $nontarvar _Ivill*) ffirst
end

* core PK nonlinear LIML estimator
cap program drop EstimatePK
program define EstimatePK
	syntax [if], [bs femaleonly maleonly TECHnique(string) *]
	if "`technique'"=="" local technique dfp nr
	cmp  ($depvar = $creditvars $nontarvar $covsy $extracovsy _Ivill*) ///
	     `=cond("`maleonly'"==""  , "($creditw = $covsf)", "")' ///
		 `=cond("`femaleonly'"=="", "($creditm = $covsm)", "")' ///
		 [pw=weightpk] `if', tech(`technique') `=cond("`bs'"=="", "cluster(nh)", "")' ///
		 ind("$indy" `=cond("`maleonly'"=="","choicefpk*indf","")' `=cond("`femaleonly'"=="","choicempk*indm", "")') `options' ghkdraws(101, anti) nolr
end

* PK nonlinear LIML with credit pooled across gender
cap program drop EstimatePKpooled
program define EstimatePKpooled
	syntax [if], [bs *]
	cmp  ($depvar = $creditvars $nontarvar $covsy $extracovsy _Ivill*) ///
	     ($credit = $covsfm) ///
		 [pw=weightpk] `if', tech(dfp nr) `=cond("`bs'"=="", "cluster(nh)", "")' ///
		 ind("$indy" choicepk*indfm) `options' ghkdraws(101, anti)
end

* Perform PK-type estimate, then "flip" estimates and rerun search starting from that point. Save both results.
cap program drop EstimatePKFlip
program define EstimatePKFlip
	syntax [if], runname(string) [maleonly femaleonly pooled nolendersplit noOUTREG *]
	EstimatePK`pooled' `if', `options' `maleonly' `femaleonly'
	SaveEstimates, runname(`runname') addstat(`"`addstat'"') `outreg'
	tempname b
	mat `b' = e(b)
	if "`lendersplit'" == "" mat `b' = -`b'[1,1..6], `b'[1,7...]
		else mat `b' = -`b'[1,1..2], `b'[1,3...]
	if "`maleonly'`femaleonly'" == "" {
		mat `b'[1, colsof(`b') - 2] = -`b'[1, colsof(`b') - 2]
		mat `b'[1, colsof(`b') - 1] = -`b'[1, colsof(`b') - 1]
	}
	else {
		mat `b'[1, colsof(`b')] = -`b'[1, colsof(`b')]
	}
	EstimatePK`pooled' `if', init(`b') `options' `maleonly' `femaleonly'
	SaveEstimates, runname(`runname'_flip) addstat(`"`addstat'"') `outreg'
end

* wrapper for estimation routines that starts from estimates of a given run, estimates on bootstrap sample, and drops coefficients for village dummies to keep returned parameter vector stable
cap program drop EstimateForBootstrap
program EstimateForBootstrap, eclass
	syntax [if], estimator(string) runname(string) [*]
	tempname b V b_bs t
	tempvar sample

	if "`estimator'" == "PK" {
		est restore `runname'
		mat `b' = e(b)
		qui Estimate`estimator', bs `options' noest quietly
		mat `b_bs'=e(b)
		local bsparams: colfullnames `b_bs'
		forvalues i = 1/`=colsof(`b_bs')' {
			cap mat `t' = `b'[1,"`: word `i' of `bsparams''"]
			if !_rc mat `b_bs'[1, `i'] = `t'
		}
		Estimate`estimator' `if', bs `options' init(`b_bs') quietly
	}
	else Estimate`estimator' `if', bs `options'
	
	mat `b' = e(b)
	mat `V' = e(V)
	cap drop `sample'
	gen byte `sample' = e(sample)
	local bnames : colnames `b'
	forvalues i=`=colsof(`b')'(-1)1 { // drop entries for village dummies since the collection varies by iteration
		if strpos("`: word `i' of `bnames''", "_Ivill") {
			mat `b' = `b'[1,1..`=`i'-1'], `b'[1,`=`i'+1'...]
			mat `V' = (`V'[1..`=`i'-1',1..`=`i'-1'], `V'[1..`=`i'-1',`=`i'+1'...]) \ (`V'[`=`i'+1'...,1..`=`i'-1'], `V'[`=`i'+1'...,`=`i'+1'...])
		}
	}
	local ll = e(ll)
	eret post `b' `V', esample(`sample')
	ereturn scalar ll = `ll'
end

* run EstimatePKFlip and keep estimate with higher likelihood
cap program drop EstimateForBootstrapHigherMode
program define EstimateForBootstrapHigherMode, eclass
	syntax [if], runname(string) [*]
	qui EstimateForBootstrap `if', estimator(PK) runname(`runname') `options'
	mat b = e(b)
	mat V = e(V)
	scalar ll1 = e(ll)
	qui EstimateForBootstrap `if', estimator(PK) runname(`runname'_flip) `options'
	mat b2 = e(b)
	di b[1,1] " " b[1,2] " " b[1,3] " " ll1 " " b2[1,1] " " b2[1,2] " " b2[1,3] " " e(ll)
	if e(ll) < ll1 eret repost b=b V=V // if flipped estimate has lower LL, restore initial one
end

cap program drop SaveEstimates
program define SaveEstimates
	syntax, runname(string) [addstat(string) noOUTREG]
	est store `runname'
	if "`outreg'" != "" exit
	
	if `"`addstat'"' != "" local addstat `" `addstat' , "'
	local bnames : colnames e(b)

	if e(cmd)=="cmp" {
		tempname sd skew kurtosis skewp kurtosisp temp Utest vardfp vardmp skewdfp skewdmp
		cap drop ey
		qui predict ey if e(sample), resid
		sum ey [aw=weightpk], detail
		scalar `skew' = string(r(skewness), "%12.4f")
		scalar `kurtosis' = string(r(kurtosis), "%12.4f")
		scalar `sd' = string(r(sd), "%12.4f")
		sktest ey [aw=weightpk]
		mat `Utest' = r(Utest)
		local addstat `"`addstat' "rmse",`=`sd'',"skew",`=`skew'',"kurtosis",`=`kurtosis'',"skew p",`=`Utest'[1,1]',"kurtosis p",`=`Utest'[1,2]',"Log likelihood",`=e(ll)', "'

/*		est store `temp'
		preserve
		keep if e(sample)
		bootstrap vardf=r(vardf) vardm=r(vardm) skewdf=r(skewdf) skewdm=r(skewdm), seed(987654321) cluster(villagepk) reps(1000) saving(PKSkewDiff, replace): DistDiff
		restore
		est restore `temp'
		est drop `temp'
		preserve
		use PKSkewDiff, clear
		foreach stat in vardf vardm skewdf skewdm {
			count if `stat'<0
			scalar ``stat'p' = ((r(N)<=1000/2)*r(N) + (r(N)>1000/2)*(1000 - r(N)))/1000
		}
		restore
		DistDiff
		local addstat `" `addstat' "Var diff f", `r(vardf)', "Var diff f p", `=`vardfp'', "Var diff m", `r(vardm)', "Var diff m p", `=`vardmp'', "Skew diff f", `r(skewdf)', "Skew diff f p", `=`skewdfp'', "Skew diff m", `r(skewdm)', "Skew diff m p", `=`skewdmp'', "'
*/
		if e(k_dv) == 3 {
			local rho12 = string(tanh(_b[/atanhrho_12]), "%12.4f")
			local rho12se = string( _se[/atanhrho_12]/cosh(_b[/atanhrho_12])^2, "%12.4f")
			local rho13 = string(tanh(_b[/atanhrho_13]), "%12.4f")
			local rho13se = string( _se[/atanhrho_13]/cosh(_b[/atanhrho_13])^2, "%12.4f")
			test [atanhrho_12]_cons [atanhrho_13]_cons
			local rhobothp = string(r(p), "%12.4f")
			local addstat `"`addstat' "rho f",`=string(`rho12', "%12.4f")',"rho f se",`rho12se',"rho m",`=string(`rho13',"%12.4f")',"rho m se",`rho13se',"rho both p",`rhobothp', "'
		}
	}
	else if e(cmd) == "ivreg2" {
		local addstat `"`addstat' "Kleibergen-Paap underid p", e(idp), "Kleibergen-Paap weak ID F", e(widstat), "'
	}

	if "$extracovsy" != "" & (e(cmd)=="cmp" | e(cmd)=="bootstrap") {
		unab extracovsy : $extracovsy
		cap test `: list extracovsy & bnames'
		if !_rc {
			local addstat `"`addstat' "F test extra covs p",`=string(r(p), "%12.4f")', "'
		}
	}
	local addstat `"`addstat' "nothing", 0 "'
	unab t : _Ivill*
	local t `t' _cons
	local bnames: list uniq bnames
	outreg2 `:list bnames - t' using "$TempOutregFile", onecol noparen noaster ctitle("`runname'") dec(10) adec(10) addstat(`addstat')
end

cap program drop ResetGlobals
program define ResetGlobals
	global depvar lpcnsexppk
	global  covsy scoheadpk  afedhighpk  amedhighpk  afadultdpk  amadultdpk  sexheadpk  ageheadpk  llandbefpk  edheadpk  spsislnddpk  spbrolnddpk  spparlnddpk  hdsislnddpk  hdbrolnddpk  hdparlnddpk  _Iwave*
	global  covsf scoheadpk1 afedhighpk1 amedhighpk1             amadultdpk1 sexheadpk1 ageheadpk1 llandbefpk1 edheadpk1 spsislnddpk1 spbrolnddpk1 spparlnddpk1 hdsislnddpk1 hdbrolnddpk1 hdparlnddpk1 _Ivill*
	global  covsm scoheadpk1 afedhighpk1 amedhighpk1 afadultdpk1             sexheadpk1 ageheadpk1 llandbefpk1 edheadpk1 spsislnddpk1 spbrolnddpk1 spparlnddpk1 hdsislnddpk1 hdbrolnddpk1 hdparlnddpk1 _Ivill*
	global covsfm scoheadpk1 afedhighpk1 amedhighpk1 afadultdpk1 amadultdpk1 sexheadpk1 ageheadpk1 llandbefpk1 edheadpk1 spsislnddpk1 spbrolnddpk1 spparlnddpk1 hdsislnddpk1 hdbrolnddpk1 hdparlnddpk1 _Ivill*
	global creditvars lfbraclvpk1 lfbrdblvpk1 lfgramlvpk1 lmbraclvpk1 lmbrdblvpk1 lmgramlvpk1
	global creditw lfproglvpk1
	global creditm lmproglvpk1
	global credit  lproglvpk1
	global nontarvar nontarpk
	global extracovsy
	global insts zf* zm* choicefpk1 choicempk1
end

cap program drop DefineIndicators
program define DefineIndicators
	cap drop indf indm indfm
	gen byte indf = cond($creditw<=ln(1000)+.001, $cmp_left, $cmp_cont)
	gen byte indm = cond($creditm<=ln(1000)+.001, $cmp_left, $cmp_cont)
	cap gen byte indfm = cond($credit<=ln(1000)+.001, $cmp_left, $cmp_cont)
end

* Dichotomize borrowing and adjust cmp model specification
cap program drop Probitize
program define Probitize
	foreach suffix in "" 1 {
		foreach var in fbrac mbrac fbrdb mbrdb fgram mgram { 
			replace l`var'lvpk`suffix' = l`var'lvpk`suffix'>0
		}
	}
	bysort nh: replace lfproglvpk1 = fproglvpk[1] > 0
	by     nh: replace lmproglvpk1 = mproglvpk[1] > 0
	replace lfproglvpk1ln1 = lfproglvpk1
	replace lmproglvpk1ln1 = lmproglvpk1
	replace lfproglvpk = fproglvpk > 0
	replace lmproglvpk = mproglvpk > 0
	replace indf = $cmp_probit
	replace indm = $cmp_probit
end

* load and prepare household-level data set
cap program drop PrepHHData
program define PrepHHData
	qui {
*		odbc load, exec("select * from [Roodman & Morduch HH]") dsn(PK) clear
	 use "Roodman & Morduch HH 2.dta", clear
		keep if wave<4 & pksample
		sort nh wave
		ren pksamplempk choicempk
		ren pksamplefpk choicefpk
		gen byte choicepk = choicefpk | choicempk
		gen byte progidpk = mod(int((villagepk-10)/80)+1, 4) // 1 = BRAC village, 2 = BRDB, 3 = Grameen, 0 = none

		xi i.wave i.villagepk
		gen double fasset = fnlasset + flandvalb

		foreach var in fproglvpk mproglvpk {
			gen double l`var' = ln(`var')
			recode l`var' (. = `=cond(inlist("`var'", "fproglvpk", "mproglvpk"), `=C_t', `=C_v')' )
		}
		replace llandbefpk = log(123) if nh==321111 & wave==3 // fix data error. log land=log(10123) in wave 3 and log(123) in waves 1-2
		gen double landbefpk = exp(llandbefpk)
		gen double pcnsexppk = exp(lpcnsexppk)

		bysort villagepk: egen byte fprogvillpk = max(choicefpk) // dummy for villages with female credit groups
		by     villagepk: egen byte mprogvillpk = max(choicempk) // ditto for male
		gen byte progvillpk = fprogvillpk | mprogvillpk

		* For first-stage eqs, make variables based on first-wave data only
		sort nh wave
		foreach var in choicepk choicefpk choicempk lfproglvpk lmproglvpk lfbraclvpk lfbrdblvpk lfgramlvpk lmbraclvpk lmbrdblvpk lmgramlvpk scoheadpk afedhighpk amedhighpk afedhighnspk amedhighnspk afadultdpk amadultdpk sexheadpk ageheadpk llandbefpk edheadpk edheadnspk hdbrolnddpk hdsislnddpk hdparlnddpk spbrolnddpk spsislnddpk spparlnddpk monlypk fonlypk partpk {
			cap by nh: gen double `var'1 = `var'[1]
		}
		foreach prog in brac brdb gram {             // credit variables aggregated accross gender
			bysort nh (wave): gen double l`prog'lvpk1 = ln(`prog'lvpk[1])
			recode l`prog'lvpk1 (. = 0)
		}
		bysort nh: gen double lproglvpk1 = ln(proglvpk[1])
		recode lproglvpk1 (. = `=ln(1000)')

		recode lfproglvpk lmproglvpk (0 = `=ln(1000)')     // First-stage LHS variables take censoring *threshold* (not censoring *value*)
		bysort nh (wave): gen double lfproglvpk1ln1 = log(fproglvpk[1])
		bysort nh (wave): gen double lmproglvpk1ln1 = log(mproglvpk[1])
		bysort nh (wave): gen double lproglvpk1ln1 = log(proglvpk[1])
		recode l*proglvpk1ln1 (. = 0)
		
		bysort nh (wave): gen byte nontrgth_dejure = landbefpk[1] >= 50.01

		ResetGlobals
		DefineIndicators

		* interact controls with choice dummies to make instruments
		foreach var of varlist $covsf {
			gen double zf`var' = `var' * choicefpk
		}
		foreach var of varlist $covsm {
			gen double zm`var' = `var' * choicempk
		}
		foreach var of varlist $covsf afadultdpk1 _Ivill* {
			cap gen double zb`var' = `var' * choicepk
		}
	}
end

cap program drop DistDiff
program define DistDiff, rclass
	sum ey if choicefpk [aw=weightpk], detail
	scalar varf = r(Var)
	scalar skewf = r(skewness)
	sum ey if !choicefpk [aw=weightpk], detail
	return scalar vardf = cond(varf - r(Var)<., varf - r(Var), 0)
	return scalar skewdf = cond(skewf - r(skewness)<., skewf - r(skewness), 0)
	sum ey if choicempk [aw=weightpk], detail
	scalar varm = r(Var)
	scalar skewm = r(skewness)
	sum ey if !choicempk [aw=weightpk], detail
	return scalar vardm = cond(varm - r(Var)<., varm - r(Var), 0)
	return scalar skewdm = cond(skewm - r(skewness)<., skewm - r(skewness), 0)
end