1simnhlogit=function(theta,lnprices,Xexpend) {
2#   function to simulate non-homothetic logit model
3#       creates y  a n x 1 vector with indicator of choice (1,...,m)
4#	lnprices is n x m array of log-prices faced
5#       Xexpend is n x d array of variables predicting expenditure
6#
7#       non-homothetic model specifies ln(psi_i(u))= alpha_i - exp(k_i)u
8#
9#       structure of theta vector:
10#           alpha  (m x 1)
11#           k      (m x1 )
12#           gamma  (k x 1)   expenditure function coefficients
13#           tau   -- scaling of v
14#
15   m=ncol(lnprices)
16   n=nrow(lnprices)
17   d=ncol(Xexpend)
18   alpha=theta[1:m]
19   k=theta[(m+1):(2*m)]
20   gamma=theta[(2*m+1):(2*m+d)]
21   tau=theta[length(theta)]
22   iotam=c(rep(1,m))
23   c1=as.vector(Xexpend%*%gamma)%x%iotam-as.vector(t(lnprices))+alpha
24   c2=c(rep(exp(k),n))
25   u=callroot(c1,c2,.0000001,20)
26   v=alpha - u*exp(k)-as.vector(t(lnprices))
27   vmat=matrix(v,ncol=m,byrow=TRUE)
28   vmat=tau*vmat
29   Prob=exp(vmat)
30   denom=Prob%*%iotam
31   Prob=Prob/as.vector(denom)
32# draw y
33   y=vector("double",n)
34   ind=1:m
35   for (i in 1:n)
36        {
37        yvec=rmultinom(1,1,Prob[i,])
38        y[i]=ind%*%yvec
39        }
40
41return(list(y=y,Xexpend=Xexpend,lnprices=lnprices,theta=theta,prob=Prob))
42
43}