/* EVALMNL.G: evaluation of multinomial logit likelihood and derivatives

              By John Rust, Yale University, March, 1997.    */

    proc (6)=evalmnl(dt,dloc,iloc,cloc,q,ll,dll,ls,lle,pnp,pcp,h,hll);

    local y,x,p,n,i,lnsum,d,cp;     /* number of choices in MNL model is n  */
    local dcp,hcp,ex;

    n=3;                            /* Specify special case of trinomial
                                       logit model */

    y=dt[.,dloc];                   /* dependent variable y in {1,...,n}    */

    x=dt[.,iloc'];                  /* independent variables in x: n copies
                                       for each alternative 1,...,n, i.e.
                                       x=(x(1),...,x(n)) */
    d=rows(iloc)/n;

    if (d /= rows(q));
       "Error: number of independent variables specified does
               not match the dimension of q vector";
        goto skip;
    endif;

    p=zeros(rows(x),n);
    lnsum=zeros(rows(x),1);

    i=1; do until i > n;           /* construct choice probabilities       */

       lnsum=lnsum+exp(x[.,1+d*(i-1):d*i]*q);

    i=i+1; endo;

    i=1; do until i > n;

     if (i == 1);

       p[.,1]=ones(rows(x),1)./lnsum;
       cp=(y.==1).*p[.,1];
       dcp=((y.==1)-p[.,1]).*x[.,1+d*(i-1):d*i];
       hcp=moment(x[.,1+d*(i-1):d*i].*sqrt(p[.,1]),1);
       ex=x[.,1+d*(i-1):d*i].*p[.,1];

     else;

       p[.,i]=exp(x[.,1+(i-1)*d:d*i]*q)./lnsum;
       cp=cp+(y.==i).*p[.,i];
       dcp=dcp+((y.==i)-p[.,i]).*x[.,1+d*(i-1):d*i];
       hcp=hcp+moment(x[.,1+d*(i-1):d*i].*sqrt(p[.,i]),1);
       ex=ex+x[.,1+d*(i-1):d*i].*p[.,i];

     endif;

    i=i+1; endo;

    ll=ll+sumc(ln(cp));            /* cumulate log-likeihood        */

    dll=dll+sumc(dcp);             /* derivative of log-likelihood          */

    if ls > 0; goto skip; endif;   /* don't compute hessian in linesearch   */

    if lle == -1;                  /* last likelihood evaluation iteration  */

                                   /* percent correctly predicted           */
       pcp=pcp+sumc((maxindc(p').==y));

                                   /* percent 'naively' predicted           */
       pnp=sumc(y.== 1)+pnp;

                                   /* cumulate outer product of 1st derivs  */
       h=h+moment(dcp,1);

    endif;
                                   /* cumulate hessian of likelihood        */
    hll=hll-hcp+moment(ex,1);
/*    hll=hll-moment(dcp,1); */

    skip:

    retp(ll,dll,pnp,pcp,h,hll);

    endp;