Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.

# Re: st: Simulating Multinomial Logit in Stata

 From Maarten Buis To statalist@hsphsun2.harvard.edu Subject Re: st: Simulating Multinomial Logit in Stata Date Tue, 12 Mar 2013 11:56:38 +0100

```On Mon, Mar 11, 2013 at 6:06 PM, Ali Hashemi wrote:
> I'm trying to simulate multinomial logit data in Stata. In this
> particular case, each individual has to pick from a set of 10
> alternatives defined by a combination of two attributes (x1 and x2).
> My question is regarding the way I should create the choice variable.
> Here is what I currently do:
>
> 1) Calculate the probability of each alternative as exp(XB)/(1+exp(XB))
> 1) Calculate the cumulative probability of each alternative (CDF) for
> each individual
> 2) Generate a random number (r)  for each individual (which is fixed
> across all the alternatives that he faces)
> 3) The first alternative with the cumulative probability greater than
> this random number would the choice. (CDF>r)
>
>
> I would appreciate your comment on the above procedure. Also, I'm not
> clear why don't we simply calculate the probability of each
> alternative as exp(XB)/(1+exp(XB)) and then pick the alternative with
> highest probability? Your help is much appreciated.

Below is an example of how to simulate a multinomial process in Stata.
You can obviously do what you proposed, but -mlogit- will not recover
those parameters, as it estimates a different model.

In the example I have also added some code on how I would analyse the
results of such a simulation. That part requires -simsum- (SSC)
-qplot- (SJ) and simpplot (SSC). You can use -findit- to find those
programs and install them.

*------------------ begin example ------------------
clear all
set more off

program define sim
drop _all
set obs 500
gen x1 = rnormal()
gen x2 = rnormal()

gen xb1 = ln(.5) + ln(1.5)*x1 + ln(.5)*x2
gen xb2 = ln(2) + ln(.33)*x1 + ln(2)*x2

gen p1 = 1        / ( 1 + exp(xb1) + exp(xb2) )
gen p2 = exp(xb1) / ( 1 + exp(xb1) + exp(xb2) )
gen p3 = exp(xb2) / ( 1 + exp(xb1) + exp(xb2) )

gen u = runiform()

gen y = cond(u < p1, 1, ///
cond(u < p1 + p2, 2, 3 ))

mlogit y x1 x2, base(1) rrr
end

simulate  b1_2= _b[2:x1]  b2_2= _b[2:x2]  b0_2= _b[2:_cons] ///
se1_2=_se[2:x1] se2_2=_se[2:x2] se0_2=_se[2:_cons] ///
b1_3= _b[3:x1]  b2_3= _b[3:x2]  b0_3= _b[3:_cons] ///
se1_3=_se[3:x1] se2_3=_se[3:x2] se0_3=_se[3:_cons] ///
, reps(20000) : sim

// get an overview of how well -mlogit- performed
simsum b1_2, se(se1_2) true(ln(1.5) ) mcse
simsum b2_2, se(se2_2) true(ln(0.5) ) mcse
simsum b0_2, se(se0_2) true(ln(0.5) ) mcse
simsum b1_3, se(se1_3) true(ln(0.33)) mcse
simsum b2_3, se(se2_3) true(ln(2)   ) mcse
simsum b0_3, se(se0_3) true(ln(2)   ) mcse

// take a closer look at the sampling distributions of the coefficients
qplot b1_2 b2_2 b0_2,   trscale(invnormal(@)) ///
yline(`=ln(1.5)' `=ln(.5)'  ) xline(0)
qplot b1_3 b2_3 b0_3,   trscale(invnormal(@)) ///
yline(`=ln(2)'    `=ln(.33)') xline(0)

// take a closer look at the p-values
gen p1_2 = 2*normal(- abs(( b1_2 - ln(1.5 ) ) / se1_2) )
gen p2_2 = 2*normal(- abs(( b2_2 - ln(0.5 ) ) / se2_2) )
gen p0_2 = 2*normal(- abs(( b0_2 - ln(0.5 ) ) / se0_2) )

gen p1_3 = 2*normal(- abs(( b1_3 - ln(0.33) ) / se1_3) )
gen p2_3 = 2*normal(- abs(( b2_3 - ln(2   ) ) / se2_3) )
gen p0_3 = 2*normal(- abs(( b0_3 - ln(2   ) ) / se0_3) )

gen id = _n
reshape long p1_ p2_ p0_, i(id) j(eq)

label define eq  2 "2 versus 1" ///
3 "3 versus 1"
label value eq eq
label var p0_ "constant"
label var p1_ "x1"
label var p2_ "x2"

simpplot p0_ p1_ p2_, by(eq)       ///
scheme(s2color) legend(cols(3)) ///
ylabel(,angle(horizontal))
*------------------- end example -------------------
(For more on examples I sent to the Statalist see:
http://www.maartenbuis.nl/example_faq )

---------------------------------
Maarten L. Buis
WZB
Reichpietschufer 50
10785 Berlin
Germany

http://www.maartenbuis.nl
---------------------------------
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/
```