version 9.2
#delimit;
set more off;
  quietly log;
  local logon = r(status);
  if "`logon'" == "on" {; log close; };
log using boehmke2009jop-montecarlo, replace text;


/*	****************************************************************	*/
/*     	File Name:	boehmke2009jop-analyze.do				*/
/*     	Date:   	December 01, 2011					*/
/*      Author: 	Frederick J. Boehmke					*/
/*      Purpose:	This file combines the Monte Carlo results to create	*/
/*			Table 1 and Figure 1 comparing the unconditional and 	*/
/*			conditional dyadic EHA models. It assumes the Monte	*/
/*			Carlo trials have been run in the folder real/ and	*/
/*			simulated/.						*/
/*      Input File:	real/montecarlo-dyadsXX-sims.dta,			*/
/*      		simulated/montecarlo-dyadsXX-sims.dta			*/
/*      Output File:	boehmke2009jop-montecarlo.log, 				*/
/*			boehmke2009jop-montecarlo.dta,				*/
/*			boehmke2009jop-montecarlo-collapsed.dta,		*/
/*			boehmke2009jop-montecarlo-fig01.gph,			*/
/*			boehmke2009jop-montecarlo-table01.csv,			*/
/*	Version:	Stata 9 or above.					*/
/*	****************************************************************	*/


	/************************************************************************/
	/* Combine the data from the MCs. Each model is run by a separate 	*/
	/* file, so we have to combine the results from each.			*/
	/************************************************************************/

	
		/* First, the data when the independent variables are 		*/
		/* simulated. Each contains the results for different values of */
		/* rho, the correlation between the independent variables. 	*/
	
clear;
save boehmke2009jop-montecarlo, replace emptyok;

forvalues num=1/12 {;

  if `num' < 10 {;

	use simulated/montecarlo-dyads0`num'-sims.dta, clear;

	};

  else if `num' >= 10 {;

	use simulated/montecarlo-dyads`num'-sims.dta, clear;

	};

  generat sim = `num';

  append using boehmke2009jop-montecarlo;

  save boehmke2009jop-montecarlo, replace;

  };
  
		/* Second, with the real world values of the independent	*/
		/* variables. This has correlation of about 0.99.		*/

forvalues num=1/12 {;

  if `num' < 10 {;

	use real/montecarlo-dyads0`num'-sims.dta, clear;
	
	};

  else if `num' >= 10 {;

	use real/montecarlo-dyads`num'-sims.dta, clear;

	};

  generat sim = `num';
  
  generat corrparam = 0;

  append using boehmke2009jop-montecarlo;

  save boehmke2009jop-montecarlo, replace;

  };
  

		/* Recode significance to get percentage Type I errors. 	*/
  
recode sig_* (-1=1);

foreach var of varlist b_* {;

  generat sd_`var' = `var';

  };

foreach var of varlist se_* {;

  replace `var' = `var'^2;

  };

  save boehmke2009jop-montecarlo, replace;

  
	/************************************************/
	/* Arrange the data to generate Table 1. 	*/
	/************************************************/
	
  
collapse (mean) b_* se_* sig_* numobs (sd) sd_*, by(sim corrparam);

  foreach var of varlist se_* {;

	replace `var' = sqrt(`var');

	};

  order b* sd* se* sig*;
  order sim corrparam *rpcpinc01* *rpcpinc02* *rpcpinc_jnt*;

  keep if sim == 1 | sim == 2 | sim == 3 | sim == 5 | sim == 10 | sim == 11 | sim == 8;
  
		/* Put them in the same order as in the table. */
  
	replace sim = 13 if sim == 8;
  
  gsort +sim -corrparam;
  
  order sim corrparam numobs;
  
  save boehmke2009jop-montecarlo-collapsed, replace;

	outsheet using boehmke2009jop-montecarlo-table01.csv, comma replace;

  
	/****************************************/
	/* Generate plots to create Figure 1. 	*/
	/****************************************/
	
	
use boehmke2009jop-montecarlo.dta, clear;

  keep sim corrparam b_rpcpinc01 b_rpcpinc02 b_rpcpinc_jnt;
  
	bysort corrparam sim: gen i=_n;
  
  reshape wide b_rpcpinc01 b_rpcpinc02 b_rpcpinc_jnt, i(i sim) j(corrparam);

  twoway kdensity b_rpcpinc011 if sim==3, scheme(s1mono) nodraw
	lcolor(gs0) lwidth(medthick) lpattern(shortdash)
    || 	kdensity b_rpcpinc013 if sim==3,
	lcolor(gs3) lwidth(medthick) lpattern(dash)
    || 	kdensity b_rpcpinc015 if sim==3,
	lcolor(gs6) lwidth(medthick) lpattern(longdash)
    || 	kdensity b_rpcpinc0110 if sim==3,
	lcolor(gs9) lwidth(medthick) lpattern(solid)
	ylabel(0(5)25, grid)
	xlabel(-.1(.1).3) xline(.1) xtitle("")
	title("Coefficients for X for State 1") l1title("Unconditional Dyadic EHA", size(medlarge))
	xsize(4) ysize(3)
	legend(label(1 "rho=0.85") label(2 "rho=0.65") label(3 "rho=0.50") label(4 "rho=0.05")
  	  cols(4) rows(1))
	name(b1sim3, replace);

  twoway kdensity b_rpcpinc021 if sim==3, scheme(s1mono)  nodraw
	lcolor(gs0) lwidth(medthick) lpattern(shortdash)
    || 	kdensity b_rpcpinc023 if sim==3,
	lcolor(gs3) lwidth(medthick) lpattern(dash)
    || 	kdensity b_rpcpinc025 if sim==3,
	lcolor(gs6) lwidth(medthick) lpattern(longdash)
    || 	kdensity b_rpcpinc0210 if sim==3,
	lcolor(gs9) lwidth(medthick) lpattern(solid)
	ylabel(0(20)80, grid)
	xlabel(-.1(.1).2) xline(0) xtitle("")
	title("Coefficients for X for State 2") l1title(" ", size(medlarge))
	xsize(4) ysize(3)
	legend(label(1 "rho=0.85") label(2 "rho=0.65") label(3 "rho=0.50") label(4 "rho=0.05")
  	  cols(4) rows(1))
	name(b2sim3, replace);

  twoway kdensity b_rpcpinc011 if sim==11, scheme(s1mono) nodraw 
	lcolor(gs0) lwidth(medthick) lpattern(shortdash)
    || 	kdensity b_rpcpinc013 if sim==11,
	lcolor(gs3) lwidth(medthick) lpattern(dash)
    || 	kdensity b_rpcpinc015 if sim==11,
	lcolor(gs6) lwidth(medthick) lpattern(longdash)
    || 	kdensity b_rpcpinc0110 if sim==11,
	lcolor(gs9) lwidth(medthick) lpattern(solid)
	ylabel(0(5)25, grid)
	xlabel(-.1(.1).3) xline(.1) xtitle("")
	title(" ") l2title("Conditional Dyadic EHA", size(medlarge))
	xsize(4) ysize(3)
	legend(label(1 "rho=0.85") label(2 "rho=0.65") label(3 "rho=0.50") label(4 "rho=0.05")
  	  cols(4) rows(1))
	name(b1sim11, replace);

  twoway kdensity b_rpcpinc021 if sim==11, scheme(s1mono) nodraw 
	lcolor(gs0) lwidth(medthick) lpattern(shortdash)
    || 	kdensity b_rpcpinc023 if sim==11,
	lcolor(gs3) lwidth(medthick) lpattern(dash)
    || 	kdensity b_rpcpinc025 if sim==11,
	lcolor(gs6) lwidth(medthick) lpattern(longdash)
    || 	kdensity b_rpcpinc0210 if sim==11,
	lcolor(gs9) lwidth(medthick) lpattern(solid)
	ylabel(0(20)80, grid)
	xlabel(-.1(.1).2) xline(0) xtitle("")
	title(" ") l1title(" ", size(medlarge))
	xsize(4) ysize(3)
	legend(label(1 "rho=0.85") label(2 "rho=0.65") label(3 "rho=0.50") label(4 "rho=0.05")
  	  cols(4) rows(1))
	name(b2sim11, replace);

grc1leg b1sim3 b2sim3 b1sim11 b2sim11, scheme(s1mono)
	cols(2) rows(2)
	imargin(zero)
	xsize(6) ysize(4.5);

	graph save boehmke2009jop-montecarlo-fig1, replace asis;

log close;
clear;
exit, STATA;
