capt log close
log using 771simul2_9, replace

* EC771: univariate simulation experiment 2   cfb 4113 rev 4123
* rev for 9.1 syntax 6317

set more off
local reps 1000
global true_mu 50

* set up model: y1 contains a homoskedastic error
*               y2 contains an error term which is heteroskedastic across regions


* define the simulation experiment
capture program drop mcsimul2
program define mcsimul2, rclass
     version 9.1
     syntax [, c(real 1)]
     tempvar e1 e2

     gen `e1' = invnorm(uniform())*`c'*zmu
     gen `e2' = invnorm(uniform())*`c'*z_factor
     replace y1 = true_y + `e1'
     replace y2 = true_y + `e2'
     ttest y1 = 0
     return scalar p1 = r(p)
     ttest y2 = 0
     return scalar p2 = r(p)
     return scalar c = `c'
end
	
	forv i=1/10 {
* 50 states' data from US Census
		qui webuse census2, clear
* true_y indicates the variable being used as base for the noisy measurement
		gen true_y = $true_mu
* z_factor is the variable generating heteroskedasticity in e2
		gen z_factor = region
		sum z_factor, meanonly
		scalar zmu = r(mean)
		qui {
			gen y1 = .
			gen y2 = .
			local c = `i'*10
			simulate c=r(c)  p1=r(p1) p2=r(p2), /// 
			saving(ccc`i',replace) nodots reps(`reps'): mcsimul2, c(`c')
		}
	}
	
* combine the results files into a single file for tabulation
use ccc1
forv i=2/10 {
	append using ccc`i'
	}
gen RfNull_1 = (1-p1)*100
gen RfNull_2 = (1-p2)*100
gen R5pc_1 = (p1<0.05)/10
gen R5pc_2 = (p2<0.05)/10
format %6.2f R5pc_1
format %6.2f R5pc_2
save ccc_1_10,replace
* get descriptives on results file
tabstat p1 p2 RfNull_1 RfNull_2,stat(mean) by(c)
tabstat R5pc_1 R5pc_2,stat(sum) by(c) nototal

set more on
log close

log2html 771simul2_9, replace