Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: mixed logit with halton sequence - program


From   "Arne Risa Hole" <arnehole@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: mixed logit with halton sequence - program
Date   Mon, 16 Jul 2007 13:41:50 +0100

Hi Alex,

I've got a paper forthcoming in the Stata Journal on estimating mixed
logit models using maximum simulated likelihood. I'll send you the
draft paper + code off-list.

Arne

On 16/07/07, Alexander Staus <Alexander.Staus@uni-hohenheim.de> wrote:
Dear Stata users,

I want to estimate a Mixed Logit Model for panel data with a dependent variable with 4
alternatives and different independent variables (some are constant over time, some are
varying over time).

First I want to estimate a model with just one term capturing the unobserved heterogeneity:

From utility-maximizing:U_{ijt} = a_i + b x_{ijt}  + e_{ijt}
i: individual
j: chosen alternative
t: time

I want to estimate the Logit function with Maximum Simulated Likelihood with Halton draws.

My program looks very similar to "Estimation of multinomial logit models with unobserved
heterogeneity using maximum simulated likelihood" by Peter Haan and Arne Uhlendorff,
Stata Journal 6 (2), 2006; rewritten for 4 alternatives and only the line with beginning *
changed, because of nonunderstanding.

my programme:

program define sim_mlogit_d0_one_random_cov
       version 9.2

       args todo b lnf
       tempvar etha1 etha2 etha3 random1 random2 random3 pi1 pi2 pi3 pi4 sum lnpi last lj
L1 L2

       tempname lnsig1 lnsig2 lnsig3 rho12 rho13 rho23 sigma1 sigma2 sigma3 cov12
cov13 cov23

       mleval `etha1' = `b', eq(1)
       mleval `etha2' = `b', eq(2)
       mleval `etha3' = `b', eq(3)
       mleval `lnsig1' = `b', eq(4) scalar
       mleval `lnsig2' = `b', eq(5) scalar
       mleval `lnsig3' = `b', eq(6) scalar
       mleval `rho12' = `b', eq(7) scalar
       mleval `rho13' = `b', eq(8) scalar
       mleval `rho23' = `b', eq(9) scalar


       qui     {

               scalar `sigma1' = (exp(`lnsig1'))^2
               scalar `sigma2' = (exp(`lnsig2'))^2
               scalar `sigma3' = (exp(`lnsig3'))^2
               scalar `cov12' = [exp(2*`rho12')-
1]/[exp(2*`rho12')+1]*(exp(`lnsig2'))*(exp(`lnsig1'))
               scalar `cov13' = [exp(2*`rho13')-
1]/[exp(2*`rho13')+1]*(exp(`lnsig3'))*(exp(`lnsig1'))
               scalar `cov23' = [exp(2*`rho23')-
1]/[exp(2*`rho23')+1]*(exp(`lnsig3'))*(exp(`lnsig2'))

               gen double `random1' = 0
               gen double `random2' = 0
               gen double `random3' = 0

               gen double `lnpi' = 0
               gen double `sum' = 0
               gen double `L1' = 0
               gen double `L2' = 0

               by hh: gen byte `last' = (_n==_N)

               gen double `pi1' = 0
               gen double `pi2' = 0
               gen double `pi3' = 0
               gen double `pi4' = 0

               }

       matrix W = (`sigma1' , `cov12' , `cov13' \ `cov12' , `sigma2' , `cov23' \ `cov13' ,
`cov23' , `sigma3')

       capture matrix L = cholesky(W)

       if _rc != 0 {
               di "Warning: cannot do Cholesky factorization of rho matrix" _rc
               }

       local l11=L[1,1]
       local l21=L[2,1]
       local l31=L[3,1]
       local l22=L[2,2]
       local l32=L[3,2]
       local l33=L[3,3]



       local repl=${draws}
       local r = 1

       forvalues r=1/$draws    {
               qui     {
                       replace `random1' = random_1`r' * `l11'
                       replace `random2' = random_2`r' * `l22'
                       replace `random3' = random_3`r' * `l33'

*                       replace `random2' = random_2`r'*`l22' + random_1`r'*`l21'

                       replace `pi1' = 1/(1+exp(`etha1'+`random1') + exp(`etha2' +
`random2') + exp(`etha3' + `random3'))
                       replace `pi2' = exp(`etha1' + `random1') * `pi1'
                       replace `pi3' = exp(`etha2' + `random2') * `pi1'
                       replace `pi4' = exp(`etha3' + `random3') * `pi1'

                       replace `lnpi' = ln(`pi1'*bet1 + `pi2'*bet2 + `pi3'*bet3 + `pi4'*bet4)

                       by hh: replace `sum' = sum(`lnpi')
                       by hh: replace `L1' = exp(`sum'[_N]) if _n==_N

                       by hh: replace `L2' = `L2' + `L1' if _n==_N
                       }
               }

       qui gen `lj' = cond(!`last',0,ln(`L2'/`repl'))
       qui mlsum `lnf' = `lj'
       if (`todo'==0 | `lnf'>=.) exit

end



I tried the same programme as in Haan/Uhlendorf for 3 alternatives, but received the same
output like for 4 alternatives:



.         ml model d0 sim_mlogit_d0_one_random_cov (dep= indep) (dep= indep) (dep= indep) > /lnsig1 /lnsig2 /lnsig3  /rho12 /rho13 /rho
> 23
r; t=0.10 11:51:27

.
.
.         ml search
initial:       log likelihood = -15670.836
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 504
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 506
Warning: cannot do Cholesky factorization of rho matrix 506
improve:       log likelihood =          0
Warning: cannot do Cholesky factorization of rho matrix 504
Warning: cannot do Cholesky factorization of rho matrix 504
rescale:       log likelihood =          0
Warning: cannot do Cholesky factorization of rho matrix 504
Warning: cannot do Cholesky factorization of rho matrix 504
--Break--




My second problem:

How to write the programme for an additional individual varying coefficient vector for some
variables:

From utility-maximizing:U_{ijt} = a_i + b_i x_{ijt} + b x_{ijt} + e_{ijt}
i: individual
j: chosen alternative
t: time


Some help?

Thanks in advance,

Alex




*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index