Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: Monte Carlo Simulations


From   "Susan Olivia" <olivia@primal.ucdavis.edu>
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. 

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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index