Monte Carlo simulation of a whole matrix using simulate and eclass

[FredPH]

I would like to compute a summary statistic of a set of groups within my data (such as age,sex) with confidence intervals. For that purpose I use monte carlo simulation drawing values from a Poisson distribution for every row in my data and then collapsing the rows to have the summary. The whole procedure works fine if the result of the simulation is just one value (using return scalar in rclass) but as soon as I try to simulate more than one result (using ereturn matrix in eclass) it is not working (see Stata code below). I get the error message: "type mismatch error in expression: e(A)". How could I simulate a whole vector or even matrix of results without more complex loops etc.?



program bootPGW, eclass
gen id=_n
sort id
gen N=_N
by id: gen DL2=floor(rpoisson(calung))
by id:gen D02=floor(rpoisson(D0))
by id:gen Dsmoking=floor(rpoisson(smoking))
by id:gen ML2=(DL2/numpyr)*1000
by id:gen AL2=(ML2-CPSIIrate)/ML2
by id:replace AL2=0 if AL2<0
by id:gen A02=1-exp(-PWGcoef*(ML2-CPSIIrate))
by id:gen A2=(AL2*DL2+A02*D02)/(DL2+D02)
gen Adeaths=totdeath*A2
collapse (sum) Adeaths=Adeaths totdeath=totdeath Dsmoking=Dsmoking, by(edu_3cat sex country year)
gen AF_PWG=Adeaths/totdeath
gen AF_simple=Dsmoking/totdeath
mkmat AF_PWG, matrix(A)
ereturn matrix A=A
end

simulate a=e(A), reps(1000) nodots seed(123): bootPGW

[Roberto Ferrer]

Cross-posted:

http://stackoverflow.com/questions/2...atrix-in-stata

Please see advice on cross-posting in the FAQ:

People posting on Statalist may also post the same question on other listservers or in web forums. There is absolutely no rule against doing that; it is not our business to constrain what you do elsewhere. But if you do post elsewhere, we ask that you provide cross-references in URL form to searchable archives. That way, people interested in your question can quickly check what has been said elsewhere and avoid posting similar comments. Being open about cross-posting saves everyone time.
Cross-posting does not affect the request elsewhere in this guide that you close threads on Statalist. If your question was answered well elsewhere, you are asked to post a cross-reference to that in a closure on Statalist.

[Brendan Halpin]

simulate can deal with multiple values being returned. See the example in the help:

simulate mean=r(mean) var=r(Var), reps(10000): lnsim, obs(100)

[FredPH]

@Brendan: Thx! This works for few variables. However, if I have hundreds of values and several variables, storing the results in a matrix isinevitable. Also, the code should be flexible towards changes in the level of aggregation (variables in collapse).

@Roberto: Sorry not reading FAQ. I included the URL in the other post.

[FredPH]

The solution has been postet at stackoverflow (see link). However, the simulation proram works for rowvectors only. But that is already much better than using single scalars.