/*  This program compares the ML estimates and 
**  standard errors for probit estimation using
**  three techniques:  Newton-Raphson, OPG, and 
**  Gauss Newton regression as suggested by 
**  Davidson and MacKinnon, 2003 p. 461  */

/*  The Tuce data  */


let data[32,5]=  
  1      2.66      20      0        0  
  2      2.89      22      0        0  
  3      3.28      24      0        0  
  4      2.92      12      0        0  
  5      4.00      21      0        1  
  6      2.86      17      0        0  
  7      2.76      17      0        0  
  8      2.87      21      0        0  
  9      3.03      25      0        0  
 10      3.92      29      0        1  
 11      2.63      20      0        0  
 12      3.32      23      0        0  
 13      3.57      23      0        0  
 14      3.26      25      0        1  
 15      3.53      26      0        0  
 16      2.74      19      0        0  
 17      2.75      25      0        0  
 18      2.83      19      0        0  
 19      3.12      23      1        0  
 20      3.16      25      1        1  
 21      2.06      22      1        0  
 22      3.62      28      1        1  
 23      2.89      14      1        0  
 24      3.51      26      1        0  
 25      3.54      24      1        1  
 26      2.83      27      1        1  
 27      3.39      17      1        1  
 28      2.67      24      1        0  
 29      3.65      21      1        1  
 30      4.00      23      1        1  
 31      3.10      21      1        0  
 32      2.39      19      1        1  
;

y=data[.,5];
const=ones(32,1);
x=const~data[.,2 3 4];


/*  Define ln likelihood function for probit */

   proc (1) = lnl(beta);
   local zhat, cdf, lnl; 
        zhat = x*beta;
        cdf = cdfn(zhat);
        lnl = (y.*ln(cdf) + (1-y).*ln(1-cdf));
        retp(sumc(lnl));
   endp;

/*  This computes the observation by observation derivatives. 
**  Column summation yields the function below. 
*/

proc (1) = Zgrad(beta);
local zhat, pdf, cdf, del;
zhat = x*beta;
        pdf = pdfn(zhat);
        cdf = cdfn(zhat);
        del=(y-cdf).*pdf./(cdf.*(1-cdf));
        del = (del.*.ones(1,cols(x)).*x);
        retp(del);
 endp;

/*  Calculate gradient */

fn grad(beta)= sumc(zgrad(beta));

/*  Calculate Hessian     */

  proc (1) = hessi(beta);
  local k, zhat, pdf, cdf, d, H;
        k = cols(x);
        zhat = x*beta;
        pdf = pdfn(zhat);
        cdf = cdfn(zhat);
        d = y.*((pdf+zhat.*cdf)/cdf.^2) + (1-y).*((pdf - zhat.*(1-cdf))/(1-cdf).^2);
        d = pdf.*d;
        H = -x'*((ones(1,k).*.d).*x);
        retp(H);
  endp;

/*  Calculate Outer Product of the Gradients matrix (BHHH)   */

  proc (1) = opg(beta);
  local k, zhat, pdf, cdf, d, H;
        h=zgrad(beta)'*zgrad(beta);
        retp(H);
  endp;

/*  Simple routine to compute the MLE iteratively using OPG */ 

ii=1;
b=invpd(x'x)*x'y;
do while abs(maxc(grad(b))) > .000001 ;
b=b+invpd(opg(b))*grad(b);
ii=ii+1;
endo;

stderr=sqrt(diag(invpd(opg(b))));
"The OPG MLE :"  b~stderr;
"Iterations:"  ii ;

/*  Similar routine to compute MLE iteratively using NR */

ii=1;
b=invpd(x'x)*x'y;
do while abs(maxc(grad(b))) > .000001 ;
b=b+invpd(-hessi(b))*grad(b);
ii=ii+1;
endo;

stderr=sqrt(diag(invpd(-hessi(b))));
"The NR MLE :"  b~stderr;
"Iterations:"  ii ;

/* Routine that computes the GNR */

gnr=1;
b=invpd(x'x)*x'y;               /* Starting value         */
ii=1;
do while abs(maxc(gnr)) > .000001; 
zhat=x*b;                       /* Index function         */
cdf=cdfn(zhat);                 /* std normal cdf         */
pdf=pdfn(zhat);                 /*  std normal pdf        */
vt=sqrt(cdf.*(1-cdf));          /* std dev of GNR errors  */ 
ystar=(y-cdf)./vt;              /*  transform y           */
xstar=((pdf./vt).*.ones(1,cols(x))).*x;  /*  transform X  */
gnr=invpd(xstar'xstar)*xstar'*ystar;   /*  The GNR        */
b=b+gnr;               /* updating the estimates          */
ii=ii+1;               /* counter to obtain # iterations  */
endo;


/*  Computation of the Covariance matrix  */

ehat=ystar-xstar*b;               /*  residuals            */
s2=ehat'ehat/(rows(x));           /*  estimated variance--
                                  **  should converge to 1 and is optional, theoretically */
cov=s2.*invpd(xstar'xstar);       /*  covariance estimate  */
se=sqrt(diag(cov));               /*  standard errors      */
"GNR Estimates " b~SE;
"iterations " ii;

{b,covb,retcode} = probit(x,y,&probitli);
beta0 = b;