Hi All,

    Sorry to bother you. I'm trying to estimate a set of discrete choice
data in R with mixed logit models where one coefficient is random and
normally distributed. I've searched on the R help archive and don't see much
information very specific to what I'm doing, so I write the code myself,
which involves simulated maximum likelihood. But it doesn't work, as I
compare the results with Limdep. Could anyone please take a look at it? It's
a bit longer though......Comments are highly appreciated. Thank you very
much!

 /*1,000 people, Three alternatives, two variables, their coefficients are
B1 and B2 respectively, B1 is random*/

    library(maxLik)

    RP<-function(theta,y,X)
{
m1<-theta[1]                /* mean
s1<-theta[2]                 /*standard deviations
b2<-theta[3]                /*B2

P<-NULL

b1<-rnorm(500,mean=m1,sd=s1)       /*generate 500 random draws for B1

for(m in 0:999)
{
Dm<-X[(1+3*m):(3+3*m),]    /*Extract the data for one person

Pn<-NULL
for(n in 1:500)
{
b<-rbind(b1[n],b2)
an<-sum(exp(Dm%*%b))
Pmn<-exp(Dm%*%b)/an         /*Under each B1, compute the choice
probabilities
Pn<-cbind(Pn,Pmn)
}
Pm<-rowMeans(Pn)    /* The simulated probabilities for one person
P<-rbind(P,Pm)   /* Obtain the choice probabilities for all 1,000 people
}
sum(log(P)*(as.numeric(y)))   /* Log-likelihood function, where y is the
variable for choices
}

A<-matrix(c(0,1,0),1,3)
B<-0
rp<-maxBFGS(RP,start=c(-0.072,0.8,0.539),y=Y,X=EV,constraints=list(ineqA=A,ineqB=B))

 /* Constrained MLE, to ensure that the estimated standard deviation is
positive; the start values are taken from conditional logit estimation


-- 
Min Chen
Graduate Student
Department of Agricultural, Food, and Resource Economics
208 Cook Hall
Michigan State University
chenmin5@msu.edu / chenmin0905@gmail.com

	[[alternative HTML version deleted]]

