Hi everyone,
I am trying to decompose variance based on random normally
distributed data and want to replicate the simulations by
say 250 times. However, when I tried to run execute the
'simul' command, I only ended up with 0 and 1 observations
[shouldn't no of obs from this simulation equals to number
of replications?]. I am attaching what I have done here and
hopefully can get some programming insight from the Stata
listers.
Thanks in advance,
Susan
. // Generating 100 clusters //
.
. qui set obs 100
. gen id_psu = _n
.
. tempvar sdalpha1 sdalpha2
. gen `sdalpha1' = uniform()
. gen `sdalpha2' = uniform()
.
. // Expand to generate individual observations - assuming
each cluster has 10 obs, which implies total obse
> rvation is 1000 //
.
. expand 10
(900 observations created)
. sort id_psu, stable
. gen no = _n
.
.
. // Randomly generate latitude and longitude //
.
. gen xcoord = uniform()*100
. gen ycoord = uniform()*100
.
. nspatwmat, name(W)xcoord(xcoord) ycoord(ycoord) band(0
100) eigenval(E) standardize cluster(id_psu)
The following matrices have been created:
1. Inverse distance weights matrix W (row-standardized)
Dimension: 1000x1000
Euclidean distances are used to generate the weight
matrix
Distance band: 0 < d <= 100
Friction parameter: 1
Minimum distance: 1.052
1st quartile distance: 32.742
Median distance: 50.582
3rd quartile distance: 70.258
Maximum distance: 131.742
Largest minimum distance: 61.042
Smallest maximum distance: 36.750
2. Eigenvalues matrix E
Dimension: 1000x1
.
. * Generate the predictor --- Note that am assuming x
doesn't change in each simulation
.
. gen x = 5 + `sdalpha1'*invnorm(uniform()) + `sdalpha2'
. mkmat x
.
.
. gen ones = 1
. mkmat ones
. matrix X = x, ones
.
. scalar lambda = 0.2
. gen constant = 5
. mkmat constant
.
. gen a = 1/1000
. scalar N = _N
. matrix Idtot = I(N)
.
. matrix IlambdaW = Idtot-(lambda*W)
. matrix invIlambdaW = inv(IlambdaW)
.
.
.
. local id_psu = id_psu
. local a = a
.
.
. program mcsimul, rclass
1.
. *** Generating error terms which to be used to generate Y
***
. *** Here the generated Y is assumed to be a function of
the error term and imposed spatial autocorrelation
> using the Weighting Matrix ***
.
.
. gen e_psu = 0.15*invnorm(uniform())
2. mkmat e_psu
3. gen e_h = sqrt( (exp( .5*x - .01*x^2) )/ (1 + exp(
5*x - .01*x^2)) )*invnorm(uniform())
4. mkmat e_h
5. gen u = e_h + e_psu
6. mkmat u
7. matrix epsilon = invIlambdaW*u
8. matrix y = constant + x + epsilon
9. svmat y
10.
. matrix drop e_psu e_h IlambdaW invIlambdaW
11.
. // OLS approach //
.
. reg y1 x
12.
. *** Decompose the variance, getting the cluster
specific plus the contribution of cluster error to t
> otal variance ***
.
. predict uols, resid
13. egen ubarols = mean(uols), by (id_psu)
14. gen difols = uols - ubarols
15. gen dif2ols = (uols - ubarols)^2
16. egen sdif2ols = total(dif2ols)
17.
. gen taucols = (1/(100*99))*sdif2ols
18. gen ubar2ols = ubarols^2
19. egen subar2ols = total(ubar2ols)
20.
. gen nume1ols = ubar2ols * a
21. gen denomols = a*(1-a)
22. gen nume2ols = a*(1-a)*taucols
23.
. egen snume1ols = total(nume1ols)
24. egen snume2ols = total(nume2ols)
25. egen sdenomols = total(denomols)
26.
. gen varetaols = (snume1ols - snume2ols)/sdenomols
27.
.
. gen ua2ols = uols^2
28. qui summ ua2ols
29. gen ua2barols = r(mean)
30. gen rmseols = sqrt(ua2barols * (1000/98))
31. gen mseols = rmseols^2
32.
. gen mseaols = e(rmse)^2
33. gen ratiools = varetaols/mseaols
34.
. qui sum varetaols
35. return scalar varetols = r(mean)
36.
. qui sum ratiools
37. return scalar ratiools = r(mean)
38.
.
.
. end
. simulate vareta=r(vareta) ratiools=r(ratiools) , reps(250)
dots: mcsimul
command: mcsimul
vareta: r(vareta)
ratiools: r(ratiools)
Simulations (250)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 50
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 100
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 150
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 200
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 250
. desc
Contains data
obs: 250 simulate:
mcsimul
vars: 2 3 Jun 2009
09:09
size: 3,000 (99.9% of memory free)
-------------------------------------------------------------------------------
storage display value
variable name type format label variable label
-------------------------------------------------------------------------------
vareta float %9.0g r(vareta)
ratiools float %9.0g r(ratiools)
-------------------------------------------------------------------------------
Sorted by:
. sum
Variable | Obs Mean Std. Dev. Min
Max
-------------+--------------------------------------------------------
vareta | 0
ratiools | 1 .0958402 . .0958402
.0958402
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
