"TNLEST.GPR: program to compute maximum likelihood estimates of
            trinomial logit model as a hint for Econ 551 midterm exam

            By John Rust, Yale University, April, 1997. ";

 /* STEP 1: Load in ASCII datasets and convert to GAUSS data sets
            with alternative-specific covariates for use in general
            form of multinomial logit model coded in EVAL.G.
	    Note the 3 coefficients (const, slope_a and slope_b) are
	    normalized to 0 for alternative 1. So estimate 3
	    coefficients for alternative 2.*/

   load dat[]=data1.asc;
   dat=reshape(dat,2000,3);
   x=ones(2000,1)~dat[.,2:3];
   let vnames=const slope_a slope_b;
   dt=saved(dat[1:1500,1]~x[1:1500,.],"data1","vnames");

   /* Step 2: set up input variables to estimate MNL model in
              DATA1.ASC, performing a derivative check first
              to make sure likelihood and derivatives are correctly
              coded in EVAL.G */

   iloc=seqa(2,1,3);
   cloc=0;
   dloc=1;
   q=zeros(6,1);
   infile="data1";
   outfile="data1.est";
   qsav="q_data1";
   dvar="y";
   let ivar=c2 sa2 sb2 c3 sa3 sb3;


/* Uncomment this line if you would like to run a derivative check,
   i.e. a comparison of numerical and analytical derivatives of the
   log-likelihood function  

   call drvchk(dloc,iloc,cloc,q,infile,2000,outfile,qsav,&evalmnl1);  */


   /* Below is the main procedure that maximizes the likelihood function */


   call max(dvar,dloc,ivar,iloc,cloc,q,infile,500,outfile,qsav,&evalmnl1); 
   

/* See how well model predicts both in and out of sample */

   /* construct 2000x3 matrix of predicted choice probabilities */

   load q=q_data1;
   x1=x[.,1:3]*q[1:3];
   x2=x[.,1:3]*q[4:6];
   cv=(dat[.,1].==1)~(dat[.,1].==2)~(dat[.,1].==3);
   cp=ones(2000,1)~exp(x1)~exp(x2);
   cp=cp./sumc(cp');

   "Mean prediction errors in sample, by alternative";
   meanc(cv[1:1500,.]-cp[1:1500,.]);

   "Mean prediction errors out of sample, by alternative";
   meanc(cv[1501:2000,.]-cp[1501:2000,.]);