/* gauss\src\ghmnl.src  
   Procedures for maxlik estimation of
      logit model
      multinomial logit with 3 choices

   For logit it provides ....
     loglli(theta,f1)                 Log-likelihood function for logit MLE
     logrd(theta,f1)                  Gradient for logit MLE
     loghess(theta,indy,indx,f1,nobs) Hessian for logit MLE

    For multinomial logit with 3 choices it provides
     mnl3lli(theta,f1)                 Log-likelihood function for logit MLE
     mnl3rd(theta,f1)                  Gradient for logit MLE
     mnl3hess(theta,indy,indx,f1,nobs) Hessian for logit MLE

    Note lli is an nobs x 1 vector  (not a scalar)
         grd is an nobs x k matrix  (not a k x 1 vector) 
         hess is a k x k matrix     (as usual) 
 
    The program assumes that data in a nobs x p matrix z are ordered as
      The first column of z is the dependent variable y
      (for logit this takes value 0 or 1 
       and for 3-choice MNL this takes value 0, 1 or 2)
      The remaining columns of z are regressors.

   Note that theta is a column vector, so x*theta is nobs x 1
*/



/********** (1) Logit model (two choices)        ***********
*/

/* (1a) Log-likelihood function for logit */
proc(1) = loglli(theta,z);
local li,x,y,p,prob;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
prob = exp(x*theta)./(1+exp(x*theta));     
li = y.*prob + (1-y).*(1-prob);
@ print y~prob~li; @
retp(li);
endp;

/* (1b) Gradient for logit */
proc(1) = loggrd(theta,z);
local gr,p,x,y,prob;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
prob = exp(x*theta)./(1+exp(x*theta));       
gr = x.*(y-prob);
retp(gr);
endp;

/* (1c) Hessian for logit */
proc(1) = loghess(theta,z);
local H,p,x,y,prob;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
prob = exp(x*theta)./(1+exp(x*theta));       
H = -(x.*prob.*(1-prob))'x;
retp(H);
endp;



/********** (2) Multinomial logit model (three choices)     ***********
*/

/* (2a) Log-likelihood function for mnl 3-choice */
proc(1) = mnl3lli(theta,z);
local li,x,y,p,d1,d2,theta1,theta2,denom,prob1,prob2;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
/* create indicator variables for y = 0, 1 or 2 */
d1 = 1*(y.==1);
d2 = 1*(y.==2);
/* partition theta, which has 2 times (p-1) entries */
theta1 = theta[1:p-1];
theta2 = theta[p:2*p-2];
denom = 1+exp(x*theta1)+exp(x*theta2);
prob1 = exp(x*theta1)./denom;       
prob2 = exp(x*theta2)./denom;
li = (1-d1-d2).*(1-prob1-prob2) + d1.*prob1 + d2.*prob2;
retp(li);
endp;

/* (2b) Gradient for mnl 3-choice */
proc(1) = mnl3grd(theta,z);
local gr,x,y,p,d1,d2,theta1,theta2,denom,prob1,prob2,gr1,gr2;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
/* create indicator variables for y = 0, 1 or 2 */
d1 = 1*(y.==1);
d2 = 1*(y.==2);
/* partition theta, which has 2 times (p-1) entries */
theta1 = theta[1:p-1];
theta2 = theta[p:2*p-2];
denom = 1+exp(x*theta1)+exp(x*theta2);
prob1 = exp(x*theta1)./denom;       
prob2 = exp(x*theta2)./denom;
gr1 = x.*(d1-prob1);
gr2 = x.*(d2-prob2);
gr = gr1~gr2;
retp(gr);
endp;

/* (2c) Hessian for mnl 3-choice */
proc(1) = mnl3hess(theta,z);
local H,x,y,p,d1,d2,theta1,theta2,denom,prob1,prob2,
      H11,H12,H22,H1,H2;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
/* create indicator variables for y = 0, 1 or 2 */
d1 = 1*(y.==1);
d2 = 1*(y.==2);
/* partition theta, which has 2 times (p-1) entries */
theta1 = theta[1:p-1];
theta2 = theta[p:2*p-2];
denom = 1+exp(x*theta1)+exp(x*theta2);
prob1 = exp(x*theta1)./denom;       
prob2 = exp(x*theta2)./denom;
H11 = -(x.*(prob1.*(1-prob1)))'x;
H12 = (x.*(prob1.*prob2))'x;
H22 = -(x.*(prob2.*(1-prob2)))'x;
H1 = H11~H12;
H2 = H12~H22;
H = H1|H2;
retp(H);
endp;

