/********************************************************************
********************   LOGIT ESTIMATION  ****************************
USER  MUST PROVIDE:
  _dataset  -  string, name of GAUSS dataset consisting of the
               dependent variable followed by the independent
               variable
USER MAY PROVIDE:
  _savedb  -  string, name of *.fmt file containing initial
              parameter values
  _saveb   -  string, name of *.fmt file for storing final
              parameter values
  _outfile -  string, name of output file
  _program -  string, name of the program file

Program written by: Alan G. Isaac
                    Associate Professor of Economics
                    The American University
                    Washington, DC 20016

This source code is provided free of any charge for public,
 non-commercial use. Please provide appropriate acknowledgment
 if it supports published research or if you incorporate
 portions of this code into your own.
This code has not been subject to formal test criteria and is
 provided without any performance guarantees.
A useful reference is
  William H. Greene, _Econometric Analysis_ (NY: Macmillan,1993).

--Alan G. Isaac, 20 July 1994


************************************************************************
***********************************************************************/

new;
_dataset = "";
_savedb = "";
_saveb = "c:\\gauss\\logitb";
_outfile = "c:\\gauss\\logit.out";
_program = "c:\\gauss\\logit.src";


output file = ^_outfile reset;
print "Results of " $_program " in " $_outfile;

if _dataset $== ""; "USER MUST SET _DATASET"; end; endif;

_names = getname(_dataset);
_dep = _names[1];
_ind = _names[2:rows(_names)];
"OUTPUT FROM THE LOGIT REGRESSION MODEL:";
_sv = SVlogit(_dataset,_dep,_ind);
{ _b,_vc,_llf } =
      newton(&llflogit,&grdlogit,&hsslogit,_dataset,_sv,.00001,100);
save ^_saveb = _b;
call results(_b,_vc,_llf,_dep,_ind,_dataset,_program);


output off;


/****************************************************************
***************  PROCEDURE DEFINITION AREA **********************
************************ BELOW **********************************
****************************************************************/


/*****************************************************************
******************************************************************
This procedure maximizes the log likelihood function using the
Newton- Raphson technique.
******************************************************************
*****************************************************************/
proc (3) = newton(&f,&grd,&hss,dataset,b0,tol,maxit);
  local f:proc,grd:proc,hss:proc;
  local b,i,inc,toltest,fhandle,nobs,vof,gradient,hessian,k1,readdata,covtheta;
    b=b0;
    i = 0;
    inc = 0;
    TOLTEST = 1;                              /* ARBITRARY */

        fhandle = -1;
        open fhandle = ^dataset;
        if fhandle == -1;
            errorlog dataset $+ " could not be opened";
            end;
        endif;
        fhandle = close(fhandle);

  DO WHILE TOLTEST>TOL AND i< MAXIT;
        open fhandle = ^dataset;
        call seekr(fhandle,1);
        nobs = 0; vof=0;gradient=0;hessian=0;

        k1 = getnr(6,colsf(fhandle));
        do until eof(fhandle);
            readdata = readr(fhandle,k1);
            vof = vof + sumc(f(b,readdata));
            gradient=gradient + sumc(grd(b,readdata))';
            hessian = hessian + hss(b,readdata);
            nobs = nobs + rows(readdata);
        endo;
        clear readdata;
        if fhandle > 0;
           fhandle = close(fhandle);
        endif;

        print "Value of function is " vof " after " i " iterations.";

        local minor;
        minor = 1;
        DO WHILE minor <= rows(b);
                if det(submat(-HESSIAN,seqa(1,1,minor),seqa(1,1,minor))) <= 0;
                     print "         implied covariance matrix not p.d.";
                     print "         adjusting covariance matrix...";
                     HESSIAN = diagrv(0*HESSIAN,-abs(diag(HESSIAN)));
                     minor = rows(b);
                endif;
                minor = minor +1;
        ENDO;

        COVTHETA = invpd(-HESSIAN);

        print "ESTIMATION: QUADRATIC APPROXIMATION.";

        inc = COVTHETA * GRADIENT';

        {GRADIENT,b,vof,TOLTEST}
         = IMPROVE(&f,&grd,b,inc,vof,dataset);
        i = i + 1;
  ENDO;

  if _saveb $/= "";
        save ^_saveb = b;
  endif;

  if TOLTEST > TOL;
     PRINT "DID NOT CONVERGE IN  " i-1 " ITERATIONS.";
        end;
  else;
     PRINT "CONVERGENCE ACHIEVED IN " i " ITERATIONS.";
     open fhandle = ^dataset;
     call seekr(fhandle,1);
     hessian=0;
     k1 = getnr(6,colsf(fhandle));
     do until eof(fhandle);
         readdata = readr(fhandle,k1);
         hessian = hessian + hss(b,readdata);
     endo;
     clear readdata;
     fhandle = close(fhandle);

     COVTHETA = invpd(-HESSIAN);
  endif;
retp(b,covtheta,vof);
endp;



/*****************************************************************
******************************************************************
This procedure is a crude search for an acceptable step size.
It also returns the relative gradient for a tolerance test.
******************************************************************
*****************************************************************/
PROC (4) = IMPROVE(&LLFn,&gradient,beta,inc,LLF,dataset);
    local LLFn:proc,gradient:proc,LLFtest,stepsize,stepnum,TOLTEST,GRAD;
     LLFtest = sum_it(&LLFn,dataset,beta + inc);
     DO WHILE NOT(FEQ(beta,beta+inc));
        stepsize = 1;
        DO WHILE (LLFtest <= LLF) AND NOT(FEQ(beta,beta+inc));
                /* "FAILURE TO IMPROVE: TRYING SMALLER STEP.";*/
                inc = inc/10;
                stepsize = stepsize +1;
                LLFtest = sum_it(&LLFn,dataset,beta + inc);
        ENDO;
        stepnum=1;
        if NOT FEQ(beta,beta+inc);
          DO WHILE LLFtest > LLF;
                stepnum = stepnum +1;
                print "IMPROVEMENT ACHIEVED IN LLFn.";
                "TRYING AN ADDITIONAL INCREMENT IN THE SAME DIRECTION...";
                LLF = LLFtest;
                LLFtest = sum_it(&LLFn,dataset,beta + inc*stepnum);
          ENDO;
         "TRYING A SMALLER STEP...";
        else;
          output on;
          print "TRYING A NEW DIRECTION...";
        endif;
        beta = beta + inc*(stepnum-1);
    ENDO;

    output on;
    print "Calculating new relative gradient for tolerance test.";
    GRAD = sum_it(&gradient,dataset,beta);
    TOLTEST = maxc(abs(GRAD' .* beta));
    /*
    print "tolerance test is now" TOLTEST;
    */
    print;
retp(GRAD,beta,LLF,TOLTEST);
endp;



/*****************************************************************
******************************************************************
This procedure returns a vector of logit log likelihoods.
All non-zero observations of the dependent varialbe are set to 1.
******************************************************************
*****************************************************************/
proc llflogit(b,dta);
    local n,c,y,x,xb;
    n = rows(dta);
    c = cols(dta);
    y = dta[.,1] .ne 0;
    x = ones(n,1)~dta[.,2:c];
    xb = x * b;
retp((y .* xb) - ln(1+exp(xb)));
endp;

/*****************************************************************
******************************************************************
This procedure returns an n by c matrix of stacked row gradients.
All non-zero observations of the dependent varialbe are set to 1.
******************************************************************
*****************************************************************/
proc grdlogit(b,dta);
    local n,c,y,x;
    n = rows(dta);
    c = cols(dta);
    y = dta[.,1] .ne 0;
    x = ones(n,1)~dta[.,2:c];
retp((y-(1./(1+exp(-x*b)))).*x);
endp;

/*****************************************************************
******************************************************************
This procedure returns a hessian for the likelihood function
given the parameters and data.
All non-zero observations of the dependent varialbe are set to 1.
******************************************************************
*****************************************************************/
proc hsslogit(b,dta);
    local n,c,y,x,xb;
    n = rows(dta);
    c = cols(dta);
    y = dta[.,1] .ne 0;
    x = ones(n,1)~dta[.,2:c];
    xb = x*b;
retp(-(x./(2+exp(-xb)+exp(xb)))'x);
endp;


/*****************************************************************
******************************************************************
This procedure returns the initial parameter values for the
optimization.
******************************************************************
*****************************************************************/
proc SVlogit(dataset,dep,ind);
    local b,b2;
    b = ls(dataset);
    if _savedb $/= "";
      load b2 = ^_savedb;
      if rows(b2) == rows(b);
        b=b2; print "Using saved parameters.";
      endif;
    endif;
retp(b);
endp;




/*****************************************************************
******************************************************************
This procedure returns the value of a function for the entire
dataset, given the parameters. The function may be vector valued
but must then return a row vector.
******************************************************************
*****************************************************************/
proc (1) = sum_it(&f,dataset,b);
   local f:proc,fhandle,sum,readdata,k1;
        open fhandle = ^dataset;
        call seekr(fhandle,1);
        sum=0;

        k1 = getnr(6,colsf(fhandle));
        do until eof(fhandle);
            readdata = readr(fhandle,k1);
            sum=sum + sumc(f(b,readdata));
        endo;
        clear readdata;
        fhandle = close(fhandle);
retp(sum);
endp;

/*****************************************************************
******************************************************************
This procedure returns the ordinary least squares estimates of
the parameters, regressing the first column of the dataset on
the remaining columns. A constant is _added_.
******************************************************************
*****************************************************************/
proc (1) = ls(dataset);
   local fhandle,x,y,mm,xy,k1,readdata,b;
        open fhandle = ^dataset;
        call seekr(fhandle,1);
        x=0;y=0;mm=0;xy=0;
        k1 = getnr(6,colsf(fhandle));
        do until eof(fhandle);
            readdata = readr(fhandle,k1);
            y = readdata[.,1];
            x = ones(rows(readdata),1)~readdata[.,2:cols(readdata)];
            mm = mm + x'x;
            xy = xy + x'y;
        endo;
        clear readdata;
        b = invpd(mm)*xy;
        fhandle = close(fhandle);
retp(b);
endp;

/*****************************************************************
******************************************************************
This procedure formats and prints the results of the optimization.
******************************************************************
*****************************************************************/
proc (0) = results(b,vc,llf,dep,ind,dataset,program);
  local fmt,mask,fh;
   call header(program,dataset,0);
   print "Dependent variable is: " $dep;
   print "There are" rows(ind) "independent variables: "; $ind';
   open fh = ^dataset;
      print "Number of observations: " rowsf(fh); fh = close(fh);
   print "Log of the likelihood function: " llf;

/*
  print "             Parameters w/ T-ratios:";
  let mask[1,3] = 0 1 1;
  let fmt[3,3] =
  "-*.*s " 8 8
  "*.*lf " 10 3
  "*.*lf " 10 3;
  call printfm(("Constant"|ind)~b~b./sqrt(diag(vc)),mask,fmt);
  print;

  print "             Parameters w/ standard errors & T-ratios:";
  let mask[1,4]= 0 1 1 1;
  call printfmt(("Constant"|ind)~b~sqrt(diag(vc))~b./sqrt(diag(vc)),mask);
*/
  print "       Parameters     T-ratios     P-value    exp(b):";
  let mask[1,5]= 0 1 1 1 1;
  let fmt[5,3] =
  "-*.*s " 8 8
  "*.*lf " 10 6
  "*.*lf " 10 3
  "*.*lf " 10 3
  "*.*lf " 10 3;
  call printfm(("Constant"|ind)~b
            ~b./sqrt(diag(vc))~2*cdfnc(abs(b)./sqrt(diag(vc)))~exp(b),mask,fmt);
  print "
=============================================================================";
  print "
=============================================================================";
endp;

end;
