* This file replicates Morduch 1998

* odbc load, exec("select * from [Roodman & Morduch HH]") dsn(PK) clear
use "Roodman & Morduch HH.dta", clear
keep if wave<4

replace landaft = 264 if nh==2620151  // according to Morduch's land13.do
replace landbef = 264 if nh==2620151

global HHcovs scohead afadultd amadultd afedhigh amedhigh sexhead agehead llandaf edhead hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd
global villcovs primcoed hpruralh hpfampln hpmidwif bamtr pzrice pzwflour pzmoil pzhegg pzmilk pzpotato wagem wagef nowagef

replace pzwflour=0 if pzwflour==. // needed to match PK sample sizes. These prices are actually stored as 0 in original data.
replace pzmilk=0 if pzmilk==.
scalar C_0 = ln(1)
quietly {
	gen pcworkf = fworkhh/fworkers
	gen pcworkm = fworkhh/mworkers
	gen pcwork = (fworkhh+mworkhh)/(mworkers+fworkers)
	gen lpcworkf = ln(pcworkf)  // Morduch appears to leave ln(0)=.
	gen lpcworkm = ln(pcworkm)  // From DIFFS1.do: /* note lsfv lsmv = . if lsma, lsfa = 0 */
	gen lpcwork = ln(pcwork)
* recode l*work* (. = 0)

	foreach var in fproglv mproglv braclv fbraclv mbraclv brdblv fbrdblv mbrdblv gramlv fgramlv mgramlv {
		gen double l`var' = ln(`var')
		recode l`var' (. = `=C_0')
	}
	replace landaf = .0001 if landaf==0 // according to Morduch's DIFFS1.DO. Pitt used 0.1.
	replace landbef = .0001 if landbef==0 
	gen double llandaf = ln(landaf)
	gen double llandbef = ln(landbef)
	gen double lpcnsexp = ln(pcnsexp)
	gen byte nontar = progid==5
*	recode hdbrolndd hdsislndd hdparlndd spbrolndd spsislndd spparlndd (. = 0)  // Pitt does this, but Morduch doesn't
}

tsset nh wave // make sure obs are ordered by wave for "first" below
collapse lpcnsexp lpcwork $HHcovs $villcovs (first) villprog nontar weight q r landaf landbef village lfproglv lmproglv lfbraclv lmbraclv ///
		lfbrdblv lmbrdblv lfgramlv lmgramlv pksample pksamplef pksamplem fworkhh mworkhh (sd) vlpcnsexp=lpcnsexp ///
		vlpcwork=lpcwork (count) Nlpcnsexp=lpcnsexp Nlpcwork=lpcwork, by(nh)

qui foreach var in lpcnsexp lpcwork {
	replace v`var' = . if N`var'<3
	drop N`var'
	replace v`var' = v`var'^2 * (3-1)/3 // -collapse- bias-corrects s.d. Undo that for comparability
	gen lv`var' = ln(v`var')
}

* Number of de jure non-target, non-borrowing HH, of which how many in control villages:
count if landbef>50 & !q
count if landbef>50 & !q & villprog==0

local eligrule de facto
* diff-diffs
* This makes Table 7 when classification=2 (de jure) and depvar=vlpcnsexp
forvalues classification=1/2 {
	xi i.nontar*i.villprog, noomit
	foreach depvar in vlpcnsexp lpcwork vlpcwork lpcnsexp {
		di _n as res "Villprog coding: 0=control, 1=BRAC, 2=BRDB, 3=Grameen"
		di _n as res "Averages (upper left)"
		reg  `depvar' _I*X* [pw=weight] if `if' 1, nocons
		
		di _n as res "Differences by program (bottom left; look at first 4 coefs)"
		reg `depvar' _InonXvil_0_* _Ivillprog_* [pw=weight] if `if' 1, nocons
		
		di _n as res "Differences by eligibility (upper right; look at first 6 coefs)"
		reg `depvar' _InonXvil_*_1 _InonXvil_*_2 _InonXvil_*_3 _Inontar_* [pw=weight] if `if' 1, nocons
		
		di _n as res "Diffs in diffs (bottom right; look at first coef in each regression"
		forvalues p=1/3 {
			reg `depvar' _InonXvil_0_`p' _Inontar_0 _Ivillprog_`p' [pw=weight] if `if' (villprog==0 | villprog==`p')
		}
	}
	local if !(q & !r) &  // exclude mistargeted for de jure
	local eligrule de jure
	replace nontar = landbef>50
}

* Table 6 top half. Look at first 3 coefs of each regression.
gen nontaraft = landaft>=50
xi i.nontar*i.villprog i.village, noomit
foreach depvar in lpcnsexp vlpcnsexp lpcwork vlpcwork {
	reg `depvar' _Ivillprog_1 _Ivillprog_2 _Ivillprog_3 $HHcovs [pw=weight] if !(q & !r) & landbef<=50, cluster(village)
	reg `depvar' _Ivillprog_1 _Ivillprog_2 _Ivillprog_3 $HHcovs $villcovs [pw=weight] if !(q & !r) & landbef<=50, cluster(village)
	reg `depvar' _InonXvil_0_1 _InonXvil_0_2 _InonXvil_0_3 nontaraft $HHcovs _Ivillage* [pw=weight] if !(q & !r), cluster(village)
}


*** Individual-level regressions

* odbc load, exec("select * from [Roodman & Morduch ind]") dsn(PK) clear
use "Roodman & Morduch ind.dta", clear
keep if wave<4

scalar C_0 = ln(1000)
replace pzwflour=0 if pzwflour==.
replace pzmilk=0 if pzmilk==.
quietly {
	foreach var in fproglv mproglv fbraclv mbraclv fbrdblv mbrdblv fgramlv mgramlv {
		gen double l`var' = ln(`var')
		recode l`var' (. = `=C_0')
	}
	replace landaf = .1 if landaf==0 // according to Pitt e-mail, 1/4/08
	replace landbef = .1 if landbef==0 
	gen double llandaf = ln(landaf)
	gen double llandbef = ln(landbef)
	gen byte nontar = progid==5
}

scalar Y_0 = ln(1)      // labor supply censoring level? (1 is lowest observed value.)

egen long id = group (nh pid)
tsset id wave // make sure obs are ordered by wave for "first" below
collapse villprog nontar weight q r medec517 fedec517 mwork fwork $HHcovs $villcovs (first) landaf landbef village lfproglv lmproglv lfbraclv lmbraclv ///
		lfbrdblv lmbrdblv lfgramlv lmgramlv pksample pksamplef pksamplem, by(nh pid)

* diff-diffs
replace nontar = landbef>50
xi i.nontar*i.villprog, noomit
foreach depvar in mwork fwork medec517 fedec517 {
	di _n as res "Villprog coding: 0=control, 1=BRAC, 2=BRDB, 3=Grameen"
	di _n as res "Averages"
	reg  `depvar' _I*X* [pw=weight] if !(q & !r), nocons
	
	di _n as res "Differences by program"
	reg `depvar' _Ivillprog_* _InonXvil_0_* [pw=weight] if !(q & !r), nocons
	
	di _n as res "Differences by eligibility."
	reg `depvar' _Inontar_* _InonXvil_*_1 _InonXvil_*_2 _InonXvil_*_3 [pw=weight] if !(q & !r), nocons
	
	di _n as res "Diffs in diffs. Look at first coef in each regression."
	forvalues p=1/3 {
		reg `depvar' _InonXvil_0_`p' _Inontar_0 _Ivillprog_`p' [pw=weight] if !(q & !r) & (villprog==0 | villprog==`p')
	}
}


* Table 6 bottom half. Look at first 3 coefs of each regression.
gen nontaraft = landaft>=50
xi i.nontar*i.villprog i.village, noomit
local regcmd reg
foreach depvar in medec517 fedec517 mwork fwork {
	if "`depvar'" == "mwork" {
		local regcmd tobit
		local ll ll
	}
	`regcmd' `depvar' _Ivillprog_1 _Ivillprog_2 _Ivillprog_3 $HHcovs [pw=weight] if !(q & !r) & landbef<=50, cluster(village) `ll'
	`regcmd' `depvar' _Ivillprog_1 _Ivillprog_2 _Ivillprog_3 $HHcovs $villcovs [pw=weight] if !(q & !r) & landbef<=50, cluster(village) `ll'
	`regcmd' `depvar' _InonXvil_0_1 _InonXvil_0_2 _InonXvil_0_3 nontaraft $HHcovs _Ivillage* [pw=weight] if !(q & !r), cluster(village) `ll'
}