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 <[email protected]> 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/
