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/ArsenicContamination_2004_BAMWSP_IDW.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(dhsidyear);
	
	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!=.;
glu;
gen nonswitcher=contaminated==1 | surfacewaterD==1 if contaminated!=. & surfacewaterD!=.;

** BGS measures of access;
ren fraction_mindistanceU_1m_G fractionclean;
ren wfraction_mindistanceU_1m_G wfractionclean;
ren mean_mindistanceU_G distclean;
ren wmean_mindistanceU_G wdistclean;

** BAMWSP measures of access: clustersafedepth lnsafedepth percentcontAs100_nw percentcontAs100_w ;
drop meanrasterdistclean;

** IDW measures of access: pctrasterover100 lnmeanrasterdistclean;
replace meanrasterdistclean200=meanrasterdistclean200/1000;
gen lnmeanrasterdistclean200=ln(meanrasterdistclean200+1);

gen wsim_meandistance=weightedsim_meandistance_simC20;
gen lnweightedsim=ln(weightedsim_meandistance_simC20+1);

*****************************************************************************************************************;
** figure 3, access to a clean well, no contamination measure;
*****************************************************************************************************************;

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

**BGS;
if 1 {;
	makegraph `switch' fractionclean, varindex(fractionclean) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' wfractionclean, 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 "`output3BGS'_fractionclean_`switch'.csv", comma replace;

	makegraph `switch' distclean if numwalkable5>=3, varindex(distclean) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' wdistclean if numwalkable5>=3, varindex(wdistclean) gridpoints(`gridpoints') numiters(`numiters');
	outsheet b_distclean x_distclean y_distclean CI_low_distclean CI_high_distclean b_wdistclean x_wdistclean y_wdistclean CI_low_wdistclean CI_high_wdistclean using "`output3BGS'_distclean_`switch'.csv", comma replace;
};

**BAMWSP;
if 1 {;
	*makegraph `switch' percentcontAs100_nw, varindex(percentcontAs100_nw) gridpoints(`gridpoints') numiters(`numiters');
	*makegraph `switch' percentcontAs100_w, varindex(percentcontAs100_w) gridpoints(`gridpoints') numiters(`numiters');
	*outsheet b_percentcontAs100_nw x_percentcontAs100_nw y_percentcontAs100_nw CI_low_percentcontAs100_nw CI_high_percentcontAs100_nw b_percentcontAs100_w x_percentcontAs100_w y_percentcontAs100_w CI_low_percentcontAs100_w CI_high_percentcontAs100_w using "`output3BAM'_percentcontAs100_`switch'.csv", nol comma replace;

	makegraph `switch' clustersafedepth, varindex(clustersafedepth) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' lnsafedepth, varindex(lnsafedepth) gridpoints(`gridpoints') numiters(`numiters');
	outsheet b_clustersafedepth x_clustersafedepth y_clustersafedepth CI_low_clustersafedepth CI_high_clustersafedepth b_lnsafedepth x_lnsafedepth y_lnsafedepth CI_low_lnsafedepth CI_high_lnsafedepth using "`output3BAM'_safedepth_`switch'.csv", comma replace;
};

*IDW;
if 1 {;
	local measure="meanrasterdistclean200";
	local measure2="lnmeanrasterdistclean";
	makegraph `switch' `measure', varindex(`measure') gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure2', varindex(`measure2') gridpoints(`gridpoints') numiters(`numiters');
	outsheet b_`measure' x_`measure' y_`measure' CI_low_`measure' CI_high_`measure' 
		b_`measure2' x_`measure2' y_`measure2' CI_low_`measure2' CI_high_`measure2'
		using "`output3IDW'_rasterdistclean_`switch'.csv", nol comma replace;
		
	*local measure="pctrasterover100";
	*makegraph `switch' `measure', varindex(`measure') gridpoints(`gridpoints') numiters(`numiters');
	*outsheet b_`measure' x_`measure' y_`measure' CI_low_`measure' CI_high_`measure' 
		using "`output3IDW'_`measure'_`switch'.csv", nol comma replace;
};

*IDW new;
if 0 {;
	local measure="distcleanpoint500";
	local measure2="distcleanpoint1000";
	makegraph `switch' `measure', varindex(`measure') gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure2', varindex(`measure2') gridpoints(`gridpoints') numiters(`numiters');
	outsheet b_`measure' x_`measure' y_`measure' CI_low_`measure' CI_high_`measure' 
		b_`measure2' x_`measure2' y_`measure2' CI_low_`measure2' CI_high_`measure2'
		using "`output3IDW'_distcleanpoint_`switch'.csv", nol comma replace;
};

*simulation new;
if 1 {;
	local measure="wsim_meandistance";
	local measure2="lnweightedsim";
	makegraph `switch' `measure' if `measure'<.5, varindex(`measure') gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure2' if `measure'<.5, varindex(`measure2') gridpoints(`gridpoints') numiters(`numiters');
	outsheet b_`measure' x_`measure' y_`measure' CI_low_`measure' CI_high_`measure' 
		b_`measure2' x_`measure2' y_`measure2' CI_low_`measure2' CI_high_`measure2'
		using "`output3IDW'_simulated_`switch'.csv", nol comma replace;
};
*****************************************************************************************************************;
** figure 4, access to a clean well, conditioning on contamination;
*****************************************************************************************************************;

** BGS: contamination options: wfraction_mindistanceC_1m_G;
local contamination1 = "wfraction_mindistanceC_1m_G";

** BAMWSP: contamination options: percentcont_w fractionknown50pctplus percentcont_nw;
local contamination2 = "percentcont_w";

** IDW: contamination options: pctrastercontaminated;
local contamination3 = "pctrastercontaminated";


*rename these just by heardofarsenic and just by medcontamination;
**BGS;
if 1 {;
sum `contamination1',d;
gen med`contamination1'=`contamination1'>0;
foreach measure in "fractionclean" "wfractionclean" {;

	makegraph `switch' `measure' if med`contamination1'==1 & numwalkable5>0, varindex(`measure'_C) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination1'==0 & numwalkable5>0, varindex(`measure'_U) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_C x_`measure'_C-CI_high_`measure'_C 
		b_`measure'_U x_`measure'_U-CI_high_`measure'_U
		using "`output3_bBGS'_`measure'_`switch'.csv", comma replace;
	
	makegraph `switch' `measure' if heardofarsenic==1, varindex(`measure'_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if heardofarsenic==0, varindex(`measure'_nh) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_h x_`measure'_h-CI_high_`measure'_h 
		b_`measure'_nh x_`measure'_nh-CI_high_`measure'_nh
		using "`output3_cBGS'_`measure'_`switch'.csv", comma replace;

	makegraph `switch' `measure' if med`contamination1'==1 & numwalkable5>0 & heardofarsenic==1, varindex(`measure'_Ch) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination1'==0 & numwalkable5>0 & heardofarsenic==1, varindex(`measure'_Uh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination1'==1 & numwalkable5>0 & heardofarsenic==0, varindex(`measure'_Cnh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination1'==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 "`output4BGS'_`measure'_`switch'.csv", comma replace;
		

};
drop med`contamination1';	
};

if 0 {;
**BAMWSP;
gen percentcont100_nw=percentcontAs100_nw ;
sum `contamination2',d;
*lnsafedepth2 not actually used;
gen lnsafedepth2=lnsafedepth if `contamination2'>r(mean) & clustersafedepth>7;
replace lnsafedepth2=lnsafedepth if `contamination2'<=r(mean);

foreach measure in "lnsafedepth" {;
*"percentcontAs100_w" "percentcont100_nw" 

	sum `contamination2',d;
	gen high`contamination2'=`contamination2'>r(mean);

	makegraph `switch' `measure' if high`contamination2'==1, varindex(`measure'_C) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination2'==0, varindex(`measure'_U) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_C x_`measure'_C-CI_high_`measure'_C 
		b_`measure'_U x_`measure'_U-CI_high_`measure'_U
		using "`output3_bBAM'_`measure'_`switch'.csv", comma replace;
	
	makegraph `switch' `measure' if heardofarsenic==1, varindex(`measure'_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if heardofarsenic==0, varindex(`measure'_nh) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_h x_`measure'_h-CI_high_`measure'_h 
		b_`measure'_nh x_`measure'_nh-CI_high_`measure'_nh
		using "`output3_cBAM'_`measure'_`switch'.csv", comma replace;
	
	sum `contamination2',d;
	gen med`contamination2'=`contamination2'>0;

	makegraph `switch' `measure' if med`contamination2'==1 & heardofarsenic==1, varindex(`measure'_Ch) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination2'==0 & heardofarsenic==1, varindex(`measure'_Uh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination2'==1 & heardofarsenic==0, varindex(`measure'_Cnh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination2'==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 "`output4BAM'_`measure'_`switch'.csv", comma replace;
	drop med`contamination2';
	drop high`contamination2';
};
drop percentcont100_nw;
};

if 0 {;
**IDW;
gen lnrstrdistclean=lnmeanrasterdistclean;
gen distclean500=distcleanpoint500;
foreach measure in "lnrstrdistclean" {;
*"lnrstrdistclean" "pctrasterover100" {;

	sum `contamination3',d;
	gen high`contamination3'=`contamination3'>r(mean);

	makegraph `switch' `measure' if high`contamination3'==1, varindex(`measure'_C) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination3'==0, varindex(`measure'_U) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_C x_`measure'_C-CI_high_`measure'_C 
		b_`measure'_U x_`measure'_U-CI_high_`measure'_U
		using "`output3_bIDW'_`measure'_`switch'.csv", comma replace;
	
	makegraph `switch' `measure' if heardofarsenic==1, varindex(`measure'_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if heardofarsenic==0, varindex(`measure'_nh) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_h x_`measure'_h-CI_high_`measure'_h 
		b_`measure'_nh x_`measure'_nh-CI_high_`measure'_nh
		using "`output3_cIDW'_`measure'_`switch'.csv", comma replace;
	
	makegraph `switch' `measure' if high`contamination3'==1 & heardofarsenic==1, varindex(`measure'_Ch_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination3'==0 & heardofarsenic==1, varindex(`measure'_Uh_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination3'==1 & heardofarsenic==0, varindex(`measure'_Cnh_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination3'==0 & heardofarsenic==0, varindex(`measure'_Unh_h) 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 "`output4IDW'_`measure'_high_`switch'.csv", comma replace;

	sum `contamination3',d;
	gen med`contamination3'=`contamination3'>0;

	makegraph `switch' `measure' if med`contamination3'==1 & heardofarsenic==1, varindex(`measure'_Ch) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination3'==0 & heardofarsenic==1, varindex(`measure'_Uh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination3'==1 & heardofarsenic==0, varindex(`measure'_Cnh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination3'==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 "`output4IDW'_`measure'_`switch'.csv", comma replace;
	drop med`contamination3';
	drop high`contamination3';
};
drop lnrstrdistclean distclean500;
};

if 0 {;
**Simulations;
**could do this with other contaminatino measures;
gen wsim_meand=weightedsim_meandistance_simC20;
foreach measure in "wsim_meand" {;

	sum `contamination2',d;
	gen high`contamination2'=`contamination2'>r(mean);

	makegraph `switch' `measure' if high`contamination2'==1& `measure'<.5, varindex(`measure'_C) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination2'==0& `measure'<.5, varindex(`measure'_U) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_C x_`measure'_C-CI_high_`measure'_C 
		b_`measure'_U x_`measure'_U-CI_high_`measure'_U
		using "`output3_bIDW'_`measure'_`switch'.csv", comma replace;
	
	makegraph `switch' `measure' if heardofarsenic==1& `measure'<.5, varindex(`measure'_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if heardofarsenic==0& `measure'<.5, varindex(`measure'_nh) gridpoints(`gridpoints') numiters(`numiters');

	drop se*;
	order b_*;
	outsheet b_`measure'_h x_`measure'_h-CI_high_`measure'_h 
		b_`measure'_nh x_`measure'_nh-CI_high_`measure'_nh
		using "`output3_cIDW'_`measure'_`switch'.csv", comma replace;
	
	makegraph `switch' `measure' if high`contamination2'==1 & heardofarsenic==1& `measure'<.5, varindex(`measure'_Ch_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination2'==0 & heardofarsenic==1& `measure'<.5, varindex(`measure'_Uh_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination2'==1 & heardofarsenic==0& `measure'<.5, varindex(`measure'_Cnh_h) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if high`contamination2'==0 & heardofarsenic==0& `measure'<.5, varindex(`measure'_Unh_h) 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 "`output4IDW'_`measure'_high_`switch'.csv", comma replace;

	sum `contamination2',d;
	gen med`contamination2'=`contamination2'>0;

	makegraph `switch' `measure' if med`contamination2'==1 & heardofarsenic==1& `measure'<.5, varindex(`measure'_Ch) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination2'==0 & heardofarsenic==1& `measure'<.5, varindex(`measure'_Uh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination2'==1 & heardofarsenic==0& `measure'<.5, varindex(`measure'_Cnh) gridpoints(`gridpoints') numiters(`numiters');
	makegraph `switch' `measure' if med`contamination2'==0 & heardofarsenic==0& `measure'<.5, 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 "`output4IDW'_`measure'_`switch'.csv", comma replace;
	drop med`contamination2';
	drop high`contamination2';
};
drop wsim_meand;
};
