/*  This program implements the Manski IV bound and SV bound estimators (including inference) in BSV (2012)  */
/*  Date:  October 9, 2010  */
/*  Authors:  Jay Bhattacharya, Azeem Shaikh, and Ed Vytlacil  */
 

/*  This version replaces the max verison of the Q(theta) definition with a sum of squares version */
#delimit ;
capture clear;
capture log close;
set mem 50m;
set seed 210391;
set matsize 500;
set more off;
capture log close;
log using "manskibound_subsample95 v5.log", replace;

global alpha 0.05;
global crit1 = invnorm(1 - $alpha/2);
global crit2 = invnorm(1 - $alpha);
global instrument weekday2;

global numsubs = 1000;
global subsampsz = 50;

*;
*   	Manski bounds;
*;

capture program drop manski;
program define manski;

	/*  This program calculates the 90th and 95th percentiles of the Q(theta) distribution  */;
	/*  for Manski bounds using the subsample draws in sub_sample.dta.   			*/;
	/*  These are the critical values in the test in BSV paper.     			*/;
	
	args t;  /*  The single argument is a number reflecting mortality horizon  */;
	         /*  e.g. if t=7, then the program uses dth7                       */;

	preserve;  /*  The program assumes that the data set consists of subsample draws from rhc.dta */;

	/*  Calculate pieces of Qb(theta) for each subsample */;
	gen dy_`t' = swang*dth`t';			/*  E[DY]          */;
	gen notdy_`t' = (1-swang)*dth`t';		/*  E[(1-D)Y]      */;
	gen notdypd_`t' = (1-swang)*dth`t' + swang; 	/*  E[(1-D)*Y + D] */;
	gen dypnotd_`t' = swang*dth`t' + (1-swang);	/*  E[DY + (1-D)]  */;
	gen mlb_`t' = dy_`t' - notdypd_`t';		/*  Manski LB	   */;
	gen mub_`t' = dypnotd_`t' - notdy_`t';		/*  Manski UB      */; 

	/* We only need the means (conditional on the instrument) from the subsample draws */;
	collapse (count) n_`t' = dth`t' (sd) sd_dy = dy_`t' sd_notdy = notdy_`t' sd_notdypd = notdypd_`t' sd_dypnotd = dypnotd_`t' 
				sd_mlb = mlb_`t' sd_mub = mub_`t' (mean) dy_`t' notdy_`t' notdypd_`t' dypnotd_`t', by(subsampnum $instrument); 
 
	/* Put the conditional means from each subsample draws on a single row */;		
	quietly reshape wide n_`t' sd_dy sd_notdy sd_notdypd sd_dypnotd sd_mlb sd_mub dy_`t' notdy_`t' notdypd_`t' dypnotd_`t', i(subsampnum) j($instrument);
	
	/* Evaluate Q(theta) at different values of theta for each subsample. 	*/;
	/* The grid over theta goes from -1 to 1 by 0.02.			*/;
	quietly expand 101;
	sort subsamp;
	gen theta = -1.0 + 0.02*(_n-1) - 2.02*(subsampnum-1);
	quietly replace theta = 0 if abs(theta) < 1.0e-12;
	
	/* Calculate Q(theta) for each subsample using theta and the conditional means. */;
	manski_qtheta `t';	
	quietly replace q_theta = q_theta* sqrt($subsampsz);

	/* Calculate the 90th and 95th percentile of the Q(theta) distribution for each theta */;
	collapse (p90) manski`t'_90 = q_theta (p95) manski`t'_95 = q_theta, by(theta);

	/* Save the result and restore the original subsample draws dataset */;
	sort theta;	
	quietly save manski`t'_temp.dta, replace;
	restore;
	
end;

capture program drop svbound;
program define svbound;

	/*  This program calculates the 90th and 95th percentiles of the Q(theta) distribution  */;
	/*  for SV bounds using the subsample draws in sub_sample.dta.   			*/;
	/*  These are the critical values in the test in BSV paper.     			*/;
	
	args t;  /*  The single argument is a number reflecting mortality horizon  */;
	         /*  e.g. if t=7, then the program uses dth7                       */;

	preserve;  /*  The program assumes that the data set consists of subsample draws from rhc.dta */;

	/*  Calculate pieces of Qb(theta) for each subsample */;
	gen d_`t' = swang;				/*  E[D]           */;
	gen dy_`t' = swang*dth`t';			/*  E[DY]          */;
	gen notdy_`t' = (1-swang)*dth`t';		/*  E[(1-D)Y]      */;
	gen notdypd_`t' = (1-swang)*dth`t' + swang; 	/*  E[(1-D)*Y + D] */;
	gen dypnotd_`t' = swang*dth`t' + (1-swang);	/*  E[DY + (1-D)]  */;
	gen mlb_`t' = dy_`t' - notdypd_`t';		/*  Manski LB	   */;
	gen mub_`t' = dypnotd_`t' - notdy_`t';		/*  Manski UB      */; 
	
	/* We only need the means (conditional on the instrument) from the subsample draws */;
	collapse (count) n_`t' = dth`t' (sd) sd_d`t' = d_`t' sd_dypnotd = dypnotd_`t' sd_notdy = notdy_`t' sd_dy = dy_`t' sd_mlb = mlb_`t' sd_mub = mub_`t'
			sd_notdypd = notdypd_`t' sd_dth`t' = dth`t' (mean) d_`t' dth`t' dy_`t' notdy_`t' notdypd_`t' dypnotd_`t', by(subsampnum $instrument); 

	/* Put the conditional means from each subsample draws on a single row */;		
	quietly reshape wide n_`t' sd_d`t' sd_dypnotd sd_notdy sd_dy sd_notdypd sd_mlb sd_mub sd_dth`t' d_`t' dth`t' notdy_`t' dypnotd_`t' dy_`t' notdypd_`t', i(subsampnum) j($instrument);

	/* Evaluate Q(theta) at different values of theta for each subsample. 	*/;
	/* The grid over theta goes from -1 to 1 by 0.02.			*/;
	quietly expand 101;
	sort subsamp;
	gen theta = -1.0 + 0.02*(_n-1) - 2.02*(subsampnum-1);
	quietly replace theta = 0 if abs(theta) < 1.0e-12;

	
	/* Calculate Q(theta) for each subsample using theta and the conditional means. */;
	sv_qtheta `t';	
	quietly replace q_theta1 = q_theta1* sqrt($subsampsz);
	quietly replace q_theta2 = q_theta2* sqrt($subsampsz);

	/* Calculate the 90th and 95th percentile of the Q(theta) distribution for each theta */;
	collapse (p90) sv1bound`t'_90 = q_theta1 sv2bound`t'_90 = q_theta2 (p95) sv1bound`t'_95 = q_theta1 sv2bound`t'_95 = q_theta2, by(theta);

	gen svbound`t'_90 = sv1bound`t'_90;
	gen svbound`t'_95 = sv1bound`t'_95;

	drop sv1bound`t'_95 sv2bound`t'_95 sv1bound`t'_90 sv2bound`t'_90;
	
	/* Save the result and restore the original subsample draws dataset */;
	sort theta;	
	quietly save svbound`t'_temp.dta, replace;
	restore;
	
end;


capture program drop manski2;
program define manski2;

	/*  This program calculates Q(theta) for Manski bounds using the the original dataset rhc.dta	*/;
	
	args t numobs;  /*  The argument t is a number reflecting mortality horizon  	  */;
	         	/*  e.g. if t=7, then the program uses dth7                       */;
	         	/*  The argument numobs is the number of observations in the      */;
	         	/*  original data.						  */;

	preserve;  
	
	gen dy_`t' = swang*dth`t';			/*  E[DY]          */;
	gen notdy_`t' = (1-swang)*dth`t';		/*  E[(1-D)Y]      */;
	gen notdypd_`t' = (1-swang)*dth`t' + swang; 	/*  E[(1-D)*Y + D] */;
	gen dypnotd_`t' = swang*dth`t' + (1-swang);	/*  E[DY + (1-D)]  */;
	gen mlb_`t' = dy_`t' - notdypd_`t';		/*  Manski LB	   */;
	gen mub_`t' = dypnotd_`t' - notdy_`t';		/*  Manski UB      */; 

	/* We only need the means (conditional on the instrument) from the subsample draws */;
	collapse (count) n_`t' = dth`t' (sd) sd_dy = dy_`t' sd_notdy = notdy_`t' sd_notdypd = notdypd_`t' sd_dypnotd = dypnotd_`t' 
				sd_mlb = mlb_`t' sd_mub = mub_`t' (mean) dy_`t' notdy_`t' notdypd_`t' dypnotd_`t', by(subsampnum $instrument); 
 
	/* Put the conditional means from each subsample draws on a single row */;		
	gen id = 1;
	quietly reshape wide n_`t' sd_dy sd_notdy sd_notdypd sd_dypnotd sd_mlb sd_mub dy_`t' notdy_`t' notdypd_`t' dypnotd_`t', i(subsampnum) j($instrument);


	/* Evaluate Q(theta) at different values of theta 	*/;
	/* The grid over theta goes from -1 to 1 by 0.02.	*/;
	quietly expand 101;
	gen theta = -1.0 + 0.02*(_n-1);
	quietly replace theta = 0 if abs(theta) < 1.0e-12;
	
	/* Calculate Q(theta) using theta and the conditional means. */;
	manski_qtheta `t';	
	quietly replace q_theta = q_theta * sqrt(`numobs');
	
	/* Save the result and restore the original dataset */;
	sort theta;	
	rename q_theta q_theta`t';
	manski_point `t';
	summ bl_manski`t' bu_manski`t';
	drop dy_`t'0 notdy_`t'0 notdypd_`t'0 dypnotd_`t'0 dy_`t'1 notdy_`t'1 notdypd_`t'1 dypnotd_`t'1 id;
	quietly save manski`t'_rhc.dta, replace;
	restore;
	
end;

capture program drop svbound2;
program define svbound2;

	/*  This program calculates Q(theta) for Manski bounds using the the original dataset rhc.dta	*/;
	
	args t numobs;  /*  The argument t is a number reflecting mortality horizon  	  */;
	         	/*  e.g. if t=7, then the program uses dth7                       */;
	         	/*  The argument numobs is the number of observations in the      */;
	         	/*  original data.						  */;

	preserve;  
	
	/*  Calculate pieces of Qb(theta) for each subsample */;
	gen d_`t' = swang;				/*  E[D]           */;
	gen dy_`t' = swang*dth`t';			/*  E[DY]          */;
	gen notdy_`t' = (1-swang)*dth`t';		/*  E[(1-D)Y]      */;
	gen notdypd_`t' = (1-swang)*dth`t' + swang; 	/*  E[(1-D)*Y + D] */;
	gen dypnotd_`t' = swang*dth`t' + (1-swang);	/*  E[DY + (1-D)]  */;
	gen mlb_`t' = dy_`t' - notdypd_`t';		/*  Manski LB	   */;
	gen mub_`t' = dypnotd_`t' - notdy_`t';		/*  Manski UB      */; 


	/* We only need the means (conditional on the instrument) from the subsample draws */;
	collapse (count) n_`t' = dth`t' (sd) sd_d`t' = d_`t' sd_dypnotd = dypnotd_`t' sd_notdy = notdy_`t' sd_dy = dy_`t' sd_mlb = mlb_`t' sd_mub = mub_`t'
			sd_notdypd = notdypd_`t' sd_dth`t' = dth`t' (mean) d_`t' dth`t' dy_`t' notdy_`t' notdypd_`t' dypnotd_`t', by(subsampnum $instrument); 

	/* Put the conditional means from each subsample draws on a single row */;		
	gen id = 1;
	quietly reshape wide n_`t' sd_d`t' sd_dypnotd sd_notdy sd_dy sd_notdypd sd_mlb sd_mub sd_dth`t' d_`t' dth`t' notdy_`t' dypnotd_`t' dy_`t' notdypd_`t', i(subsampnum) j($instrument);

	/* Evaluate Q(theta) at different values of theta 	*/;
	/* The grid over theta goes from -1 to 1 by 0.02.	*/;
	quietly expand 101;
	gen theta = -1.0 + 0.02*(_n-1);
	quietly replace theta = 0 if abs(theta) < 1.0e-12;
	
	/* Calculate Q(theta) using theta and the conditional means. */;
	sv_qtheta `t';	
	quietly replace q_theta1 = q_theta1 * sqrt(`numobs');
	quietly replace q_theta2 = q_theta2 * sqrt(`numobs');
	quietly gen q_theta = q_theta1;
	
	/* Save the result and restore the original dataset */;
	sort theta;	
	rename q_theta q_theta`t';
	sv_point `t';
	summ bl_svbound`t' bu_svbound`t';
	drop d_`t'0 d_`t'1 notdy_`t'0 dypnotd_`t'0 notdy_`t'1 dypnotd_`t'1 dy_`t'1 notdypd_`t'1 dy_`t'0 notdypd_`t'0 id deltay_`t' dth`t'0 dth`t'1;
	quietly save svbound`t'_rhc.dta, replace;
	restore;
	
end;



capture program drop manski_qtheta;
program define manski_qtheta;

	/* This program calculates q_theta for the Manski IV bounds estimator using the formulas in BSV (2012) */;
	/*  t is the number of days since hospital admission */;
	
	args t;

	gen se1_00 = sd_mlb0/sqrt(n_`t'0);
	gen se1_11 = sd_mlb1/sqrt(n_`t'1);
	gen se1_01 = sqrt((sd_dy0^2)/n_`t'0 + (sd_notdypd1^2)/n_`t'1);
	gen se1_10 = sqrt((sd_dy1^2)/n_`t'1 + (sd_notdypd0^2)/n_`t'0);
	
	gen se2_00 = sd_mub0/sqrt(n_`t'0);
	gen se2_11 = sd_mub1/sqrt(n_`t'1);
	gen se2_01 = sqrt((sd_dypnotd0^2)/n_`t'0 + (sd_notdy1^2)/n_`t'1);
	gen se2_10 = sqrt((sd_dypnotd1^2)/n_`t'1 + (sd_notdy0^2)/n_`t'0);
	
	quietly gen q_theta = 	max((dy_`t'0 - notdypd_`t'0 - theta)/se1_00, 0)^2 +
				max((dy_`t'1 - notdypd_`t'0 - theta)/se1_10, 0)^2 +
				max((dy_`t'0 - notdypd_`t'1 - theta)/se1_01, 0)^2 +
				max((dy_`t'1 - notdypd_`t'1 - theta)/se1_11, 0)^2 +
				max((theta - dypnotd_`t'0 + notdy_`t'0)/se2_00, 0)^2 +
				max((theta - dypnotd_`t'1 + notdy_`t'0)/se2_10, 0)^2 +
				max((theta - dypnotd_`t'0 + notdy_`t'1)/se2_01, 0)^2 +
				max((theta - dypnotd_`t'1 + notdy_`t'1)/se2_11, 0)^2
				;
		
end;

capture program drop sv_qtheta;
program define sv_qtheta;

	/* This program calculates q_theta for the SV bounds estimator using the formulas in BSV (2012) */;
	/*  t is the number of days since hospital admission */;
	
	args t;
	
	quietly gen deltad`t' = d_`t'1 - d_`t'0;
	quietly gen deltay_`t' = dth`t'1 - dth`t'0;
	quietly gen se_deltad = sqrt( (((sd_d`t'1)^2) / n_`t'1) + (((sd_d`t'0)^2) / n_`t'0) ) ; 
	quietly gen se_delta = sqrt( (((sd_dth`t'1)^2) / n_`t'1) + (((sd_dth`t'0)^2) / n_`t'0) ) ;
*	summ deltay_`t' deltad`t';
	
	quietly gen se_positive10 = sqrt((sd_dypnotd1^2 / n_`t'1) + (sd_notdy0^2 / n_`t'0));
	quietly gen se_positive11 = sd_mub1/sqrt(n_`t'1);
	quietly gen se_positive01 = sqrt((sd_dypnotd0^2 / n_`t'0) + (sd_notdy1^2 / n_`t'1));
	quietly gen se_positive00 = sd_mub0/sqrt(n_`t'0);
	
	quietly gen se_negative10 = sqrt((sd_dy1^2 / n_`t'1) + (sd_notdypd0^2 / n_`t'0));
	quietly gen se_negative11 = sd_mlb1/sqrt(n_`t'1);
	quietly gen se_negative01 = sqrt((sd_dy0^2 / n_`t'0) + (sd_notdypd1^2 / n_`t'1));
	quietly gen se_negative00 = sd_mlb0/sqrt(n_`t'0);

	quietly gen q_theta1 = 	max(-deltay_`t'/se_delta, 0)^2 +
				max((deltay_`t' - theta)/se_delta, 0)^2 +
				max((theta - dypnotd_`t'1 + notdy_`t'0)/se_positive10, 0)^2  if theta > 0;

	quietly replace q_theta1 = 	max(deltay_`t'/se_delta, 0)^2 +
					max((-deltay_`t' + theta)/se_delta, 0)^2 +
					max((- theta + dy_`t'1 - notdypd_`t'0)/se_negative10, 0)^2 if theta < 0;


	quietly replace q_theta1 =  abs(deltay_`t')/se_delta   if theta == 0;


	quietly gen q_theta2 = max(	abs(deltad`t')/se_deltad, 
						(dy_`t'0 - notdypd_`t'0 - theta)/se_negative00,
						(dy_`t'1 - notdypd_`t'0 - theta)/se_negative10,
						(dy_`t'0 - notdypd_`t'1 - theta)/se_negative01,
						(dy_`t'1 - notdypd_`t'1 - theta)/se_negative11,
						(theta - dypnotd_`t'0 + notdy_`t'0)/se_positive00,	
						(theta - dypnotd_`t'1 + notdy_`t'0)/se_positive10,	
						(theta - dypnotd_`t'0 + notdy_`t'1)/se_positive01,	
						(theta - dypnotd_`t'1 + notdy_`t'1)/se_positive11	
					);

	
end;


capture program drop manski_point;
program define manski_point;

	/*  This program calculates the Manski IV bounds "point" estimator (no inference) */;
	/*  t is the number of days since hospital admission */;
	
	args t;
	
	gen bl_manski`t' = max(dy_`t'0, dy_`t'1) - min(notdypd_`t'0, notdypd_`t'1);
	gen bu_manski`t' = min(dypnotd_`t'0, dypnotd_`t'1) - max(notdy_`t'0, notdy_`t'1);
	

end;


capture program drop sv_point;
program define sv_point;

	/*  This program calculates the SV bounds "point" estimator (no inference) */;
	/*  t is the number of days since hospital admission */;
	
	args t;
	
	gen bl_svbound`t' = deltay_`t' if deltay_`t' > 0;
	gen bu_svbound`t' = dypnotd_`t'1 - notdy_`t'0 if deltay_`t' > 0;
	replace bl_svbound`t' = dy_`t'1 - notdypd_`t'0 if deltay_`t' < 0;
	replace bu_svbound`t' = deltay_`t' if deltay_`t' < 0;

end;

capture program drop calc_bounds;
program define calc_bounds;

	/*  This program runs the svbound or manskibound programs for a given group of patients */;

	args grp estim;  /*  grp is the subset of data to be analyzed (all, cat11, cat13, cat18, cat19) */;
			 /*  estim is the estimator to be used (manski, svbound, pqd) 			*/;

	use "sub_sample_`grp'.dta";

	/*  Calculate Q(theta) thresholds for the Manski bound for mortality at varying horizons  */; 
	`estim' 7;    /* For dth7   */;
	`estim' 30;   /* For dth30  */;
	`estim' 60;   /* For dth60  */;
	`estim' 90;   /* For dth90  */;
	`estim' 180;  /* For dth180 */;

	/*  Merge all the horizons in to a single dataset 	   */;
	/*  Each row has thresholds for a different value of theta */;
	use `estim'7_temp, clear;
	foreach i of numlist 30 60 90 180 {;
		quietly merge theta using `estim'`i'_temp;
		quietly summ _merge;
		if ((r(mean) != 3)|(r(sd) !=0)) {;
			display "Merge with `estim'`i'_temp.dta for `grp' failed" as error;
			exit;
		};		
		drop _merge;
		sort theta;
	};

	quietly save `estim'_thresholds_`grp', replace;

	!del `estim'7_temp.dta
	!del `estim'30_temp.dta
	!del `estim'60_temp.dta
	!del `estim'90_temp.dta
	!del `estim'180_temp.dta

	/*  Calculate Q(theta) for the original data  */;
	use "limited_data_`grp'.dta";

	`estim'2 7 _N;
	`estim'2 30 _N;
	`estim'2 60 _N;
	`estim'2 90 _N;
	`estim'2 180 _N;

	/*  Merge all the horizons in to a single dataset 	   */;
	/*  Each row has q_theta for a different value of theta    */;
	use `estim'7_rhc, clear;
	foreach i of numlist 30 60 90 180 {;
		quietly merge theta using `estim'`i'_rhc;
		quietly summ _merge;
		if ((r(mean) != 3)|(r(sd) !=0)) {;
			display "Merge with `estim'`i'_rhc.dta for `grp' failed" as error;
			exit;
		};		
		drop _merge;
		sort theta;
	};

	quietly save `estim'_qtheta_`grp'.dta, replace;

	!del `estim'7_rhc.dta
	!del `estim'30_rhc.dta
	!del `estim'60_rhc.dta
	!del `estim'90_rhc.dta
	!del `estim'180_rhc.dta

	/* Merge with threshold file */;
	quietly merge theta using `estim'_thresholds_`grp'.dta;
	quietly summ _merge;
	if ((r(mean) != 3)|(r(sd) !=0)) {;
		display "Merge with `estim'_thresholds_`grp'.dta failed" as error;
		exit;
	};
	drop _merge;

*	set more on;
*	list theta q_theta* *95;
*	set more off;
	
	/*  Characterize the confidence interval */;
	foreach i of numlist 7 30 60 90 180 {;
		quietly gen in`i'_90 = q_theta`i' <= `estim'`i'_90;
		quietly gen in`i'_95 = q_theta`i' <= `estim'`i'_95;
*		list theta q_theta`i' in`i'_90 in`i'_95 `estim'`i'_90 `estim'`i'_95;

		gen in_`estim'`i' = (bu_`estim'`i' >= theta) * (bl_`estim'`i' <= theta);
		replace in`i'_90 = 1 if in_`estim'`i' == 1;
		replace in`i'_95 = 1 if in_`estim'`i' == 1;

		gen t_in`i'_90 = theta*in`i'_90;
		replace t_in`i'_90 = . if t_in`i'_90==0;

		gen t_in`i'_95 = theta*in`i'_95;
		replace t_in`i'_95 = . if t_in`i'_95==0;

		
		list theta q_theta`i' in`i'_90 in`i'_95 in_`estim'`i'; 

		};

	
	foreach i of numlist 7 30 60 90 180 {;
		egen lb`i'_90 = min(t_in`i'_90);
		egen ub`i'_90 = max(t_in`i'_90);

		egen lb`i'_95 = min(t_in`i'_95);
		egen ub`i'_95 = max(t_in`i'_95);
	};

	save `estim'bound_raw_conf_`grp'.dta, replace;
	
	collapse 
		(max) lb7_90 lb30_90 lb60_90 lb90_90 lb180_90 lb7_95 lb30_95 lb60_95 lb90_95 lb180_95
		(min) ub7_90 ub30_90 ub60_90 ub90_90 ub180_90 ub7_95 ub30_95 ub60_95 ub90_95 ub180_95;

	/*  Put 90% and 95% confidence bounds into different rows  */;
	quietly gen id = 1;
	quietly reshape long lb7_ lb30_ lb60_ lb90_ lb180_ ub7_ ub30_ ub60_ ub90_ ub180_, i(id) j(size);
	foreach i of numlist 7 30 60 90 180 {;
		rename lb`i'_ `grp'_lb_conf`i';
		rename ub`i'_ `grp'_ub_conf`i';
	};


	drop id;

	quietly reshape long `grp'_lb_conf `grp'_ub_conf, i(size) j(period);
	quietly gen estimator = "`estim'";

	quietly replace `grp'_ub_conf = 1 if `grp'_ub_conf == .;
	quietly replace `grp'_lb_conf = -1 if `grp'_lb_conf == .;

	!del `estim'_qtheta_`grp'.dta
	!del `estim'_thresholds_`grp'.dta

	
	/* Save the final file */;
	sort size period;
	save `estim'_conf_`grp'.dta, replace;
	list;

end;


/* ******************************************************************************** */;


/*  Calculate Manski bound confidence intervals */;
calc_bounds all manski;
forval i = 11/19 {;
	calc_bounds cat`i' manski;
};

use manski_conf_all.dta;
foreach i of numlist 11/19 {;
	sort size period;
	merge size period using manski_conf_cat`i'.dta;
	quietly summ _merge;
	if ((r(mean) != 3)|(r(sd) !=0)) {;
		display "Merge with manski_conf_cat`i'.dta failed" as error;
		exit;
	};		
	drop _merge;
};

save manski_conf.dta, replace;

!del manski_conf_all.dta
!del manski_conf_cat11.dta
!del manski_conf_cat12.dta
!del manski_conf_cat13.dta
!del manski_conf_cat14.dta
!del manski_conf_cat15.dta
!del manski_conf_cat16.dta
!del manski_conf_cat17.dta
!del manski_conf_cat18.dta
!del manski_conf_cat19.dta

clear;

use manskibound_raw_conf_all.dta;
keep theta in*;
gen str group = "all";
save manskibound_raw_conf_all.dta, replace;

foreach i of numlist 11/15 17/19 {;
	use manskibound_raw_conf_cat`i';
	keep theta in*;
	gen str group = "cat`i'";
	save manskibound_raw_conf_cat`i', replace;
};

clear;

use manskibound_raw_conf_all.dta;
foreach i of numlist 11/15 17/19 {;
	append using manskibound_raw_conf_cat`i';
};

foreach i of numlist 7 30 60 90 180 {;
	foreach j of numlist 90 95 {;
		replace in`i'_`j' = . if in`i'_`j' == 0;
		if (`i' != `j') {;
			rename in`i'_`j' in`j'_`i';
		};
	};
};

reshape long in95_ in90_ in_manski, i(theta group) j(days);
rename in95_ in95;
rename in90_ in90;

replace in95 = theta if in95 == 1;
replace in90 = theta if in90 == 1;
replace in_manski = theta if in_manski == 1;
replace in_manski = . if in_manski == 0;
rename in95 in_manski95;
rename in90 in_manski90;

save manskibound_raw_conf.dta, replace;

!del manskibound_raw_conf_all.dta
!del manskibound_raw_conf_cat11.dta
!del manskibound_raw_conf_cat12.dta
!del manskibound_raw_conf_cat13.dta
!del manskibound_raw_conf_cat14.dta
!del manskibound_raw_conf_cat15.dta
!del manskibound_raw_conf_cat16.dta
!del manskibound_raw_conf_cat17.dta
!del manskibound_raw_conf_cat18.dta
!del manskibound_raw_conf_cat19.dta


/*  Calculate SV bound confidence intervals */;
calc_bounds all svbound;
forval i=11/19 {;
	calc_bounds cat`i' svbound;
};

use svbound_conf_all.dta;
foreach i of numlist 11/19 {;
	sort size period;
	merge size period using svbound_conf_cat`i'.dta;
	quietly summ _merge;
	if ((r(mean) != 3)|(r(sd) !=0)) {;
		display "Merge with svbound_conf_cat`i'.dta failed" as error;
		exit;
	};		
	drop _merge;
};

save svbound_conf.dta, replace;

!del svbound_conf_all.dta
!del svbound_conf_cat11.dta
!del svbound_conf_cat12.dta
!del svbound_conf_cat13.dta
!del svbound_conf_cat14.dta
!del svbound_conf_cat15.dta
!del svbound_conf_cat16.dta
!del svbound_conf_cat17.dta
!del svbound_conf_cat18.dta
!del svbound_conf_cat19.dta


clear;

use svboundbound_raw_conf_all.dta;
keep theta in*;
gen str group = "all";
save svboundbound_raw_conf_all.dta, replace;

foreach i of numlist 11/15 17/19 {;
	use svboundbound_raw_conf_cat`i';
	keep theta in*;
	gen str group = "cat`i'";
	save svboundbound_raw_conf_cat`i', replace;
};

clear;

use svboundbound_raw_conf_all.dta;
foreach i of numlist 11/15 17/19 {;
	append using svboundbound_raw_conf_cat`i';
};

foreach i of numlist 7 30 60 90 180 {;
	foreach j of numlist 90 95 {;
		replace in`i'_`j' = . if in`i'_`j' == 0;
		if (`i' != `j') {;
			rename in`i'_`j' in`j'_`i';
		};
	};
};

reshape long in95_ in90_ in_svbound, i(theta group) j(days);
rename in95_ in95;
rename in90_ in90;

replace in95 = theta if in95 == 1;
replace in90 = theta if in90 == 1;
replace in_svbound = theta if in_svbound == 1;
replace in_svbound = . if in_svbound == 0;
rename in95 in_sv95;
rename in90 in_sv90;
rename in_svbound in_sv;

save svbound_raw_conf.dta, replace;

!del svboundbound_raw_conf_all.dta
!del svboundbound_raw_conf_cat11.dta
!del svboundbound_raw_conf_cat12.dta
!del svboundbound_raw_conf_cat13.dta
!del svboundbound_raw_conf_cat14.dta
!del svboundbound_raw_conf_cat15.dta
!del svboundbound_raw_conf_cat16.dta
!del svboundbound_raw_conf_cat17.dta
!del svboundbound_raw_conf_cat18.dta
!del svboundbound_raw_conf_cat19.dta

/* Done! */;
set more on;