* EC771: univariate simulation experiment 2   cfb 4113 rev 4123

capt log close
log using 771simul2, replace
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 8.0
     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
		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 = .
		qui gen y2 = .
		local c = `i'*10
		simulate "mcsimul2, c(`c')" c= r(c) p1=r(p1) p2=r(p2) ///
		, saving(cc`i') reps(`reps') replace
		}
	
* combine the results files into a single file for tabulation
use cc1
forv i=2/10 {
	append using cc`i'
	}
save cc_1_10,replace
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
* 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