/* EVALMNL1.G: evaluation of multinomial logit likelihood and derivatives
              Specialized version for trinomial logit model with
              alternative-specific coefficients
              By John Rust, Yale University, March, 1997.    */

    proc (6)=evalmnl1(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;

    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 1 copy
                                       for all alternatives 1,...,n,
                                       with q's for alt 1 normalized to zero */
    d=cols(x);

    hcp=zeros(2*d,2*d);
    p=zeros(rows(x),n);
    lnsum=ones(rows(x),1);

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

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

    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];

     else;

       p[.,i]=exp(x*q[1+d*(i-2):(i-1)*d])./lnsum;
       cp=cp+(y.==i).*p[.,i];

     endif;

    i=i+1; endo;

    dcp=-(y.==1).*((x.*p[.,2])~(x.*p[.,3]))+
        (y.==2).*((x.*(ones(rows(x),1)-p[.,2]))~(-x.*p[.,3]))+
        (y.==3).*((-x.*p[.,2])~(x.*(ones(rows(x),1)-p[.,3])));


    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        */

    hcp[1:d,1:d]=moment(x.*sqrt(p[.,2].*(ones(rows(x),1)-p[.,2])),1);
    hcp[1:d,d+1:2*d]=-moment(x.*sqrt(p[.,2].*p[.,3]),1);
    hcp[d+1:2*d,1:d]=hcp[1:d,d+1:2*d];
    hcp[d+1:2*d,d+1:2*d]=moment(x.*sqrt(p[.,3].*(ones(rows(x),1)-p[.,3])),1);

    hll=hll-hcp;

    skip:

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

    endp;