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}