Friday, July 13, 2012

Simulating Multinomial logit in Stata - Updated

* This post has been corrected based on suggestions form:
* Ali Hashemi Economics & Finance Ashland University and Mateusz (comment below)
* Thanks!

* This code simulates the multinomial logit model pioneered by McFadden (1973, 1983, 1984)
* (   - 1984 paper)

* The model allows for the estimation of a model that selects between multiple binary alternatives.

* Imagine that you are at ice cream shop and you are getting a rootbeer float.  The rootbeer is given.

* Now you need to selected between multiple mutually exclusive binary and complete outcomes.

* This model is used widely in such things as choosing travel destinations, insurance policies, etc.

* The structure is a little more abstract than we are used to:

* By definition sum{j=1 to J}[(P(Yj)]=1

* P(Yj=1)=exp(X'Bj)/{1+sum{i=1 to J}[exp(X'Bi)]}

* and

* P(Yj=0)=1/{1+sum{i=1 to J}[exp(X'Bi)]}

* Thus the more easily interpretable "Log-Odds Ratio is" P(1)/P(0) = exp(X'Bj)

* Note: obviously this is a highly non-linear setup

* So let's start the simulation:

set obs 1000

* First a few explanatory variables

gen fruit_loving = rnormal()

gen sweet_tooth = rnormal()

gen bean_loving = rnormal()

gen cocoa_loving = rnormal()

* Now let us calculate the XB for each of the dependent variables:
* Vanilla is the base outcome so for simplicity we will not have any coefficient values.
gen XBvanilla = 0

gen XBchocolate = 1 -2 * fruit_loving + -.1 * sweet_tooth + 1.75 * bean_loving + 2.2 * cocoa_loving

gen XBstrawberry = -1 + 2 * fruit_loving + .1 * sweet_tooth - 1.5 * bean_loving - 2.5 * cocoa_loving

* Now let's generate the Exp(XB)

gen expvan = exp(XBvanilla)

gen expchoc = exp(XBchocolate)

gen expstraw = exp(XBstrawberry)

* Now let's generate the denominator:

gen denom = expvan + expchoc + expstraw

* Finally let us generate the probabilities of seeing any particular outcome:
gen pvan = expvan/denom
gen pchoc = expchoc/denom
gen pstraw = expstra/denom

sum p*
* We can see chocolate is the most likely ice-cream flavor followed by strawberry, followed by vanilla.

* Now let's make the draws that we observe.  This is a little tricky.

gen runif = runiform()

gen flavor = 3
replace flavor = 2 if runif < pvan+pchoc
replace flavor = 1 if runif < pvan

label define flavor 1 "Vanilla" 2 "Chocolate" 3 "Strawberry"

label var flavor flavor
tab flavor
* Looking pretty good.

mlogit flavor fruit_loving sweet_tooth bean_loving cocoa_loving, baseoutcome(1)

predict pvan_hat ,equation(#1)
predict pchoc_hat ,equation(#2)
predict pstraw_hat ,equation(#3)

sum p*
* We can see the average predicted outcome probabilites are similar

pwcorr pchoc_hat pchoc, sig
* pchoc is statiscally significant

reg pvan  fruit_loving sweet_tooth bean_loving cocoa_loving
reg pvan_hat  fruit_loving sweet_tooth bean_loving cocoa_loving
* Looks very good
* We can see that how the independent variables affect the predicted variables is extremely similar.
* This is probably the best evidence yet that the mlogit command is doing what we want it to be doing.

* And
reg pchoc  fruit_loving sweet_tooth bean_loving cocoa_loving
reg pchoc_hat  fruit_loving sweet_tooth bean_loving cocoa_loving
* Looks also good

* And
reg pstraw  fruit_loving sweet_tooth bean_loving cocoa_loving
reg pstraw_hat  fruit_loving sweet_tooth bean_loving cocoa_loving


  1. You just saved my day. Thank you so much for this amazing blog, it's such an inspirational way to learn Stata (and to teach too).

    Greetings from Chile!

  2. Hi,
    drawing the unobserved component from the normal distribution while generating the XB's isn't really ok. I mean, you wouldn't see much difference, as the parameter are estimated up to scale anyway (that's the reason why, as you said: "the estimator does not work as well as expected in terms of absolute magnitude"), but the logit model implies that the unobserved part of the utility be distributed type I extreme value. So instead of normal distribution, I believe you should use: -ln(-ln(runiform())) which is the inverse of the CDF in the EV distribution.

    1. Thanks a lot Mateusz! This comment really helped me correct an error in my thinking about these types of problems.

    2. Dear Mateusz,
      Interesting point.. however, when I tried the Gumbel variate as you suggest, my conditional logit model does not recover the parameters well.. in contrast, when using the normal distribution as provided by Francis, it works well. I guess this is because this step..

      gen runif = runiform()
      gen flavor = 3
      replace flavor = 2 if runif < pvan+pchoc
      replace flavor = 1 if runif < pvan

      .. induces a Gumbel error distribution? Or am I missing something?

  3. This simulation (by drawing numbers from a normal (0,1) is equivalent to simulating Gumbel errors with what variance?