st: How do I generate data for count models?

[Matthijs De Zwaan <m.dezwaan@gmail.com>]

Dear Stata-listers,

To get a feel for the behavior of different count model estimators, I
am trying to set up a small Monte Carlo experiment. I would like to
compare the Poisson pseudo ML estimate (-poisson, vce(robust)-) and
the NegBin2 model (-nbreg-) for variance specifications that do not
correspond to the Poisson or NegBin models, and see how they behave
for different sample sizes.

Cameron and Trivedi (2009, microeconometrics for Stata, ch. 17)
describe how to set a Poisson model using the -rpoisson()- command or
a NegBin1 model using -rgamma()- and then -rpoisson()- (see C&T,
section 17.2.2). I am not sure if I understand them correctly. Below
this text is some basic code that I use to make first a Poisson set up
and then a NegBin1 set up.

I have the following questions:
(1) is the set up below correct for those models? If not, what am I doing wrong?
(2) How could I get a set up that corresponds to the NegBin2 model?
(3) How do I make a set up where the mean is still exp(xb), but the
variance does not correspond nicely to either the Poisson or the
NegBin2 model? C&T (p. 567) mention that variance could for example be
lognormally distributed, but I have no preference for any
distribution, as long as is does not follow the NegBin2 model.

I am using Stata 10.1 for Mac OSX.

Any help is greatly appreciated,
Thanks,
Matthijs de Zwaan

**** SET UP ***
* Poisson model
clear
set obs 1000
gen x = .2*invnorm(uniform())
gen e  = .2*invnorm(uniform())
gen ypois = exp(1 + x + e) // E(y) = exp(xb)
replace ypois = rpoisson(ypois) // y~poisson(exp(xb))
summ ypois
poisson ypois x

*NegBin model
gen ynb = exp(1 + x + e) // E(y) = exp(xb)
replace ynb = rgamma(ynb,1) //
replace ynb = rpoisson(ynb) // y~NegBin(mu = exp(xb), alpha = 1)
summ ynb
nbreg ynb x, d(c)
exit
*** END ***
*
*   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/