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

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)", "")' ///
		 qui ind("$indy" `=cond("`maleonly'"=="","choicefpk*indf","")' `=cond("`femaleonly'"=="","choicempk*indm", "")') `options' ghktype(halton) ghkdraws(100) ghkanti
end

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' qui ghkdraws(100) ghkanti
end

cap program drop EstimatePKFlip
program define EstimatePKFlip
	syntax [if], runname(string) [maleonly femaleonly pooled nolendersplit *]
	EstimatePK`pooled' `if', `options' `maleonly' `femaleonly'
	SaveEstimates, runname(`runname')
	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)
end

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)
		Estimate`estimator', bs `options' noest
		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')
	}
	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 substr("`: word `i' of `bnames''", 1, 6) == "_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'...])
		}
	}
	eret post `b' `V', esample(`sample')
	ereturn scalar ll = e(ll)
end

cap program drop EstimateForBootstrapHigherMode
program define EstimateForBootstrapHigherMode, eclass
	syntax [if], runname(string)
	EstimateForBootstrap `if', estimator(PK) runname(`runname')
	mat b = e(b)
	mat V = e(V)
	scalar ll1 = e(ll)
	EstimateForBootstrap `if', estimator(PK) runname(`runname'_flip)
	if e(ll) < ll1 eret repost b=b V=V
	
noi mat list e(b)
end

cap program drop SaveEstimates
program define SaveEstimates
	syntax, runname(string)
	
	local bnames : colnames e(b)

	if e(cmd)=="cmp" {
		tempname skew kurtosis Utest
		cap drop ey
		qui predict ey, resid
		sum ey, detail
		scalar `skew' = string(r(skewness), "%9.4f")
		scalar `kurtosis' = string(r(kurtosis), "%9.4f")
		sktest ey
		mat `Utest' = r(Utest)
		qui sum ey
		local addstat `""rmse",`=string(r(sd),"%9.4f")',"skew",`=`skew'',"kurtosis",`=`kurtosis'',"skew p",`=`Utest'[1,1]',"kurtosis p",`=`Utest'[1,2]',"Log likelihood",`=e(ll)' "'

		if e(k_dv) == 3 {
			local rho12 = string(tanh([atanhrho_12]_coef[_cons]), "%9.4f")
			testnl tanh([atanhrho_12]_cons)=0
			local rho12se = string(`rho12' / sqrt(r(chi2)), "%9.4f")
			local rho13 = string(tanh([atanhrho_13]_coef[_cons]), "%9.4f")
			testnl tanh([atanhrho_13]_cons)=0
			local rho13se = string(`rho13' / sqrt(r(chi2)), "%9.4f")
			test [atanhrho_12]_cons [atanhrho_13]_cons
			local rhobothp = string(r(p), "%9.4f")
			local addstat `"`addstat',"rho f",`=string(`rho12', "%9.4f")',"rho f se",`rho12se',"rho m",`=string(`rho13',"%9.4f")',"rho m se",`rho13se',"rho both p",`rhobothp'"'
		}
	}
	else if e(cmd) == "ivreg2" {
		local 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 {
			if `"`addstat'"' != "" local addstat `" `addstat' , "'
			local addstat `"`addstat' "F test extra covs p",`=string(r(p), "%9.4f")'"'
		}
	}

	est store `runname'
	unab t : _Ivill*
	local t `t' _cons
	local bnames: list uniq bnames
	if `"`addstat'"' == ""  outreg2 `:list bnames - t' using "$resultsfile", onecol noparen noaster ctitle("`runname'") dec(10)
	  else                  outreg2 `:list bnames - t' using "$resultsfile", 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

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

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