clear
set more off 
set matsize 11000
cap log close
estimates clear

#delimit ;

*****************************************************************************************************************;
*Paths;
*****************************************************************************************************************;

*access files in subdirectory generated by data creating dofiles;
local arsenicdataset "dtafiles/HouseholdArsenic_2004.dta";

*generate subdirectory to store csv files;
cap mkdir csvfiles;

*****************************************************************************************************************;
*LOCALS and PROGRAMS;
*****************************************************************************************************************;

local figure2 "csvfiles/figure2access";
local figure3 "csvfiles/figure3accessbycontamination";

capture program drop bootlpoly;
program define bootlpoly, rclass;
	syntax varlist [if] [in], iter(string) gridpoints(string) min_xvar(string) max_xvar(string);
	
	if "`iter'"!="101"  {;preserve; };
	if "`if'"!="" {;
		keep `if';
	};
	bsample, cluster(dhsid);
	
	gen grid=`min_xvar'+(_n-1)/(`gridpoints'-1)*(`max_xvar'-`min_xvar') if _n<=`gridpoints';
	
	lpoly `varlist' `if', generate(y_`iter') nograph at(grid);
	local gridsize=r(ngrid);
	
	assert y_`iter'==. if _n>`gridpoints';
if "`iter'"=="101"  {; stop; };
	mkmat y_`iter' if _n<=`gridpoints', matrix(y_`iter');
	
	return matrix y_`iter'=y_`iter';
	
	if "`iter'"!="101"  {;restore;};
end
};

cap program drop makegraph;
program define makegraph;
	syntax varlist [if], varindex(string) gridpoints(string) numiters(string);
	
	local xvar=word("`varlist'",2);

	sum `xvar' `if';
	local min_xvar=r(min);
	local max_xvar=r(max);

	gen x_`varindex'=`min_xvar'+(_n-1)/(`gridpoints'-1)*(`max_xvar'-`min_xvar') if _n<=`gridpoints';

	lpoly `varlist' `if', generate(y_`varindex') se(se_`varindex') nograph at(x_`varindex');
	gen b_`varindex'=r(bwidth) if _n<=`gridpoints';
	set seed 12345;
	
	bootlpoly `varlist' `if', iter(1) gridpoints(`gridpoints') min_xvar(`min_xvar')
		max_xvar(`max_xvar');
	
	mat y=r(y_1);

	#delimit ;
	forvalues i=2(1)`numiters' {;
		di c(current_time);
		di `i';
		bootlpoly `varlist' `if', iter(`i') gridpoints(`gridpoints') min_xvar(`min_xvar')
			max_xvar(`max_xvar');
		mat y=(y, r(y_`i'));
	};

	svmat y;
	egen CI_low_`varindex'=rowpctile(y1-y`numiters'), p(2.5);
	egen CI_high_`varindex'=rowpctile(y1-y`numiters'), p(97.5);
	egen CI_median_`varindex'=rowpctile(y1-y`numiters'), p(50);
	egen CI_mean_`varindex'=rowmean(y1-y`numiters');
	
	drop y1-y`numiters';
	
end;

** lpolys with bootstrapped standard errors **;
use "`arsenicdataset'", clear;

/*foreach varname of varlist   as_closest- maxdepthwalkable10 
	pctcontaminated5 pctcontaminated10 mean_mindistanceU_G- wmean_mindistanceU_inC2radius_R {;
	
	replace `varname'=. if inlist(shdistr,3,46,84);
	
};

gen surfacewaterD=inlist(sourcetestedwater,31,32) if sourcetestedwater!=.;
gen surfacewaterDW=inlist(dishwashingsource,31,32) if dishwashingsource!=.;
gen contaminatedorsurface=contaminatedwell==1 | surfacewaterD==1 if contaminatedwell!=. & surfacewaterD!=.;

order dhsid shdistr 
contaminatedorsurface numwalkable5
wfraction_mindistanceC_1m_G	fraction_mindistanceC_1m_G  numcontaminated5  // contamination
heardofarsenic numarsenictest contaminatedwell redwell paintedwell // household arsenic
fraction_mindistanceU_1m_G wfraction_mindistanceU_1m_G mean_mindistanceU_G wmean_mindistanceU_G // clean access 
year 
;
gen last=1;
drop year-last;
*/;
*****************************************************************************************************************;
** figure 2, access to a clean well, no contamination measure;
*****************************************************************************************************************;

local gridpoints=50;
local numiters=100;
local switch="contaminatedorsurface";


makegraph `switch' fraction_mindistanceU_1m_G, varindex(fractionclean) gridpoints(`gridpoints') numiters(`numiters');
makegraph `switch' wfraction_mindistanceU_1m_G, varindex(wfractionclean) gridpoints(`gridpoints') numiters(`numiters');
outsheet b_fractionclean x_fractionclean y_fractionclean CI_low_fractionclean CI_high_fractionclean b_wfractionclean 
	x_wfractionclean y_wfractionclean CI_low_wfractionclean CI_high_wfractionclean using "`figure2'.csv", comma replace;


*****************************************************************************************************************;
** figure 3, access to a clean well, conditioning on contamination;
*****************************************************************************************************************;

local contamination = "wfraction_mindistanceC_1m_G";

sum `contamination',d;
gen med`contamination'=`contamination'>0;
*shortname;
local measure="fractionclean";

makegraph `switch' fraction_mindistanceU_1m_G if med`contamination'==1 & numwalkable5>0 & heardofarsenic==1, varindex(`measure'_Ch) gridpoints(`gridpoints') numiters(`numiters');
makegraph `switch' fraction_mindistanceU_1m_G if med`contamination'==0 & numwalkable5>0 & heardofarsenic==1, varindex(`measure'_Uh) gridpoints(`gridpoints') numiters(`numiters');
makegraph `switch' fraction_mindistanceU_1m_G if med`contamination'==1 & numwalkable5>0 & heardofarsenic==0, varindex(`measure'_Cnh) gridpoints(`gridpoints') numiters(`numiters');
makegraph `switch' fraction_mindistanceU_1m_G if med`contamination'==0 & numwalkable5>0 & heardofarsenic==0, varindex(`measure'_Unh) gridpoints(`gridpoints') numiters(`numiters');

drop se*;
order b_*;
outsheet b_`measure'_Ch x_`measure'_Ch-CI_high_`measure'_Ch 
	b_`measure'_Uh x_`measure'_Uh-CI_high_`measure'_Uh
	b_`measure'_Cnh x_`measure'_Cnh-CI_high_`measure'_Cnh
	b_`measure'_Unh x_`measure'_Unh-CI_high_`measure'_Unh
	using "`figure3'.csv", comma replace;
		



