Thanks Rick for the quick reply. I’m unfamiliar with Proc nlmixed. I used IML for a similar model (without the simulation part) because allows me to create my model in matrix language with all the details that I needed to include in the model. It might not be the only way to do it though. I put the main parts of the program here, without certain details, to briefly show what I’m trying to do. I hope it is clear. I have a multinomial logistic structural model. I estimate the model parameters and then use these parameters for predictions. In this model, each individual faces certain choices (say 10 choices). Data is cross sectional. * Program within Proc IML, to set up the likelihood ; start DM(b) global ( X1, X2, e1, e2, e3, D, … ); * rows are number of individuals; * D is the number of draws for the random component; do d=1 to D; * Using Cholesky decomposition; E1[,d]=t#e1[,d] ; E2[,d]=r#e1[,d] + m#e2[,d]; * e1 and e2 are standard normal; * Random coefficients are; b1[,1]=Z1*c + E1[,d] ; b2[,1]=Z2*d + E2[,d]; * where E1 and E2~MVN(0,G) and are the random component; * The main model is: Yk [,1]= b1#log(X1k) + b2#log(X2k) + b3log(X1k)^2 +…+ epsk; /* k=1,…,10, one equation per choice */ *epsk~type I extreme value therefore changes in eps across choices ~ logistic; * The logistic probability for each individual and each choice is ; prob(choice 1)=exp(Y1) / sum[exp(Yk)]; . . . prob(choice 10)=exp(Y7) / sum[exp(Yk)]; * then the likelihood for each individual in the matrix is, L[,1] = s # prob(choice k) ; * s is the observed choice for each individual; * L is conditional on a given random draw; /* L from each round (total of D rounds) are all added together (sum_L) to eventually get the approximate integral love the random component (SL) ; */ sum_L [ ,1]=sum_L [ , 1]+ L [ ,1]; * end the D loop; end; SL=sum_L / D; * The log-simulated likelihood is; SLL=sum(log(SL)); return (SLL); finish DM; Parameters to estimate are c, d, b3, t, r, m (G). I used nlpnrr for optimisation. As I mentioned earlier instead of using Cholesky, I would like to use a different method that directly estimates the covariance matrix for E1 and E2 along with the other parameters. Is there a way to do this within this program (IML) or it should be transferred to Proc nlmixed? About your second comment on the Cholesky estimates: as you mentioned the resulting covariance matrix by definition is positive semidefinite, regardless of the signs of t,r,m. But what I don't quit understand is how I should justify changing the signs of the estimates.
... View more