* Code for Roodman & Morduch JDS appendix

set seed 987654321
cmp setup
cap program drop sim
program sim, rclass
	syntax, [skew]
	drop _all
	set obs 10000
	gen cf = rnormal()>0                // exogenous female choice
	gen cm = rnormal()>0                // exogenous male choice
	drawnorm ef em                      // draw uncorrelated normal errors for female and male credit equations
	gen e2 = (ef + em + rnormal())/sqrt(3)        // second-stage error correlated with first-stage errors
	if "`skew'" != "" {                 // if specified, skew the errors to the right
		replace e2 = (invchi2(15,normal(e2))-15)/sqrt(2*15)  // skewed mean-0, variance-1 distribution
	}
	gen xstarf = cf*ef                  // latent female borrowing propensity, where credit available
	gen xstarm = cm*em                  // latent male borrowing propensity, where credit available
	gen xf = cond(xstarf>0, xstarf, 0)  // observed borrowing, censored from below at 0
	gen xm = cond(xstarm>0, xstarm, 0)  // observed borrowing, censored from below at 0
	gen y = 1*xf + 1*xm + e2            // second-stage equation. coefficients on xf and xm are of interest
	ivregress liml y (xf xm = cf cm)    // classical linear LIML/2SLS
	return scalar bf_LIML = _b[xf]
	return scalar bm_LIML = _b[xm]
	cmp (y=xf xm) (xf= ) (xm= ), nolr ind($cmp_cont "cf*cond(xf>0,$cmp_cont,$cmp_left)" "cm*cond(xm>0,$cmp_cont,$cmp_left)") tech(dfp) qui // PK estimation
	return scalar bf_PK = _b[xf]
	return scalar bm_PK = _b[xm]
	ereturn clear
end
simulate, reps(100): sim
di "Correct coefficients are 1"
sum
simulate, reps(100): sim, skew
di "Correct coefficients are 1"
sum

/* Results:
. simulate, reps(100): sim

      command:  sim
        bm_PK:  r(bm_PK)
        bf_PK:  r(bf_PK)
      bm_LIML:  r(bm_LIML)
      bf_LIML:  r(bf_LIML)

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. di "Correct coefficients are 1"
Correct coefficients are 1

. sum

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       bm_PK |       100    .9984311    .0248468    .927655    1.05477
       bf_PK |       100    .9964832    .0248407   .9439897   1.068589
     bm_LIML |       100    .9900967    .0543169   .8136391   1.094452
     bf_LIML |       100    .9957779    .0455746   .9109817   1.099662

. simulate, reps(100): sim, skew

      command:  sim, skew
        bm_PK:  r(bm_PK)
        bf_PK:  r(bf_PK)
      bm_LIML:  r(bm_LIML)
      bf_LIML:  r(bf_LIML)

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. di "Correct coefficients are 1"
Correct coefficients are 1

. sum

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       bm_PK |       100    1.123902    .0310911   1.058255   1.200891
       bf_PK |       100    1.111628    .0301147   1.023171   1.184058
     bm_LIML |       100    1.007845    .0538146   .8949308   1.158872
     bf_LIML |       100    .9899866     .049715    .867712   1.116753
*/