* This file runs the simulation reported in the R&M appendix, Table 11.
* It uses -cmp- to perform WESML-LIML-FE on Monte-Carlo-simulated data sets constructed according to Pitt's sim7.do.
* See appendix of www.pstc.brown.edu/~mp/reply.pdf.

program define sim, rclass
	clear
	set obs `=Nvill'
	gen vill = _n
	xi i.vill
	* whether a village has a program for a gender depends on a gender-specific
	* component and a common component
	gen common = invnorm(uniform())
	gen mu = sqrt(0.8)*invnorm(uniform()) + sqrt(0.2)*common
	gen mupoorf = sqrt(0.5)*invnorm(uniform()) + sqrt(0.5)*common
	gen mupoorm = sqrt(0.5)*invnorm(uniform()) - sqrt(0.5)*common
	* programs are allocated according to the mu of the poor by gender (mupoor)
	expand NHH
	gen x1 = uniform() - 0.5
	gen x2 = uniform() - 0.5
	gen land = uniform() * .7
	gen choice = land < 0.5
	* more women's programs than men's
	gen choicef = choice & (mupoorf < 1)
	gen choicem = choice & (mupoorm < 0)
	gen e = invnorm(uniform())
	gen em = sqrt(2)*e + sqrt(2)*invnorm(uniform())
	gen ef = sqrt(2)*e + sqrt(2)*invnorm(uniform())
	
	* first-stage coefficients to be estimated
	gen creditf = 2*x1 + 3*x2 + 1*land - 1 + ef + mupoorf
	gen creditm = 1*x1 + 2*x2 + 2*land - 1 + em + mupoorm
	
	replace creditf = C_0 if creditf < C_0 | choicef==0
	replace creditm = C_0 if creditm < C_0 | choicem==0

	* second-stage coefficients to be estimated
	gen y = 2*x1 + 2*x2 + 1*land+ 1.5*creditf + 0.5*creditm + 2 + e + mu
	
	bysort vill: egen byte t = max(creditf) if choicef // see if anyone with choice in the village borrowed--needed for first-stage convergence
	gen byte indf = cond(choicef & t, cond(creditf <= C_0, $cmp_left, $cmp_cont), $cmp_out)
	drop t
	by vill: egen byte t = max(creditm) if choicem
	gen byte indm = cond(choicem & t, cond(creditm <= C_0, $cmp_left, $cmp_cont), $cmp_out)
	drop t

	cmp (y = x1 x2 creditf creditm land _Ivill*) (creditf = x1 x2 land _Ivill*) (creditm = x1 x2 land _Ivill*), ind($cmp_cont indf indm) nolr qui
	return scalar creditf_liml = _b[y:creditf]
	return scalar creditm_liml = _b[y:creditm]
	
	foreach var of varlist _Ivill* x1 x2 land {
		gen byte zm`var' = choicem * `var'
		gen byte zf`var' = choicef * `var'
	}
	ivreg y x1 x2 land _Ivill* (creditf creditm = z* choice?)
	return scalar creditf_2sls = _b[creditf]
	return scalar creditm_2sls = _b[creditm]

	ivreg y x1 x2 land _Ivill* (creditf creditm = choice?)
	return scalar creditf_2sls_min = _b[creditf]
	return scalar creditm_2sls_min = _b[creditm]
end

set seed 0
cmp setup
scalar C_0 = 0    // credit censoring level
scalar Nvill = 90 // # of villages
scalar NHH = 20   // # of HH/village

simulate creditf_liml=r(creditf_liml) creditm_liml=r(creditm_liml) creditf_2sls=r(creditf_2sls) creditm_2sls=r(creditm_2sls) creditf_2sls_min=r(creditf_2sls_min) creditm_2sls_min=r(creditm_2sls_min), reps(100): sim
sum, sep(0)