-------------------------------------------------------------------------------
log: /Applications/Stata/771simul2.smcl
log type: smcl
opened on: 23 Jan 2004, 13:13:18
. 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 reg
> ions
.
.
. * define the simulation experiment
. capture program drop mcsimul2
. program define mcsimul2, rclass
1. version 8.0
2. syntax [, c(real 1)]
3. tempvar e1 e2
4.
. gen `e1' = invnorm(uniform())*`c'*zmu
5. gen `e2' = invnorm(uniform())*`c'*z_factor
6. replace y1 = true_y + `e1'
7. replace y2 = true_y + `e2'
8. ttest y1 = 0
9. return scalar p1 = r(p)
10. ttest y2 = 0
11. return scalar p2 = r(p)
12. return scalar c = `c'
13. end
.
. forv i=1/10 {
2. * 50 states' data from US Census
. webuse census2, clear
3. * true_y indicates the variable being used as base for the noisy measureme
> nt
. gen true_y = $true_mu
4. * z_factor is the variable generating heteroskedasticity in e2
. gen z_factor = region
5. sum z_factor, meanonly
6. scalar zmu = r(mean)
7. qui gen y1 = .
8. qui gen y2 = .
9. local c = `i'*10
10. simulate "mcsimul2, c(`c')" c= r(c) p1=r(p1) p2=r(p2) ///
> , saving(cc`i') reps(`reps') replace
11. }
(1980 Census data by state)
command: mcsimul2 , c(10)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(20)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(30)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(40)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(50)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(60)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(70)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(80)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(90)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
(1980 Census data by state)
command: mcsimul2 , c(100)
statistics: c = r(c)
p1 = r(p1)
p2 = r(p2)
.
. * combine the results files into a single file for tabulation
. use cc1
(simulate: mcsimul2 , c(10))
. forv i=2/10 {
2. append using cc`i'
3. }
. save cc_1_10,replace
file cc_1_10.dta saved
. 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)
Summary statistics: mean
by categories of: c (r(c))
c | p1 p2 RfNull_1 RfNull_2
---------+----------------------------------------
10 | 1.90e-14 1.86e-12 100 100
20 | .0000184 .0000911 99.99816 99.99089
30 | .0020743 .0039576 99.79257 99.60424
40 | .0217047 .0280593 97.82953 97.19407
50 | .0635537 .0737673 93.64463 92.62327
60 | .1142537 .1291062 88.57463 87.08938
70 | .1711062 .1843764 82.88938 81.56236
80 | .2225545 .239061 77.74455 76.0939
90 | .2479417 .2686859 75.20583 73.13141
100 | .295781 .2993188 70.4219 70.06812
---------+----------------------------------------
Total | .1138988 .1226424 88.61012 87.73576
--------------------------------------------------
. tabstat R5pc_1 R5pc_2,stat(sum) by(c) nototal
Summary statistics: sum
by categories of: c (r(c))
c | R5pc_1 R5pc_2
---------+--------------------
10 | 100 100
20 | 100 100
30 | 99.3 98
40 | 91.3 86.9
50 | 73.8 71.3
60 | 57.9 52.9
70 | 45.5 42.1
80 | 34.5 32.8
90 | 29.6 27.1
100 | 23.6 23.4
------------------------------
.
. set more on
. log close
log: /Applications/Stata/771simul2.smcl
log type: smcl
closed on: 23 Jan 2004, 13:28:28
-------------------------------------------------------------------------------