# st: Monte Carlo Simulations

 From "Susan Olivia" To statalist@hsphsun2.harvard.edu Subject st: Monte Carlo Simulations Date Wed, 03 Jun 2009 09:17:54 -0700

```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.

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/
```