Statalist


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

Re: st: Monte Carlo Simulations


From   Stas Kolenikov <skolenik@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Monte Carlo Simulations
Date   Wed, 3 Jun 2009 13:25:20 -0500

My guess is Stata cannot overwrite the variables you created when you
first ran your program. (Sorry, it's just too damn long to go line by
line... although I am sure you can avoid -mkmat- commands these days
with Mata.) You could

1. generate those error terms as missing before simulations, and code
-replace- within your program, or

2. -capture drop- the variables you need to create.

Type -simulate, noisily- to see what's going on within your program.

On Wed, Jun 3, 2009 at 11:17 AM, Susan Olivia <olivia@primal.ucdavis.edu> wrote:
> 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/
>



-- 
Stas Kolenikov, also found at http://stas.kolenikov.name
Small print: I use this email account for mailing lists only.

*
*   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