/* This SAS IML program conducts a duration analysis of the lengths of strikes as a function of
   the deviation of output from its trend level, an indicator of the business cycle 
   position of the economy.  The data was downloaded from the CD provided in the Greene
   textbook, Econometric Analysis, 4th Ed. Table A20.1 in his data appendix.  The data was
   originally analyzed by J. Kennan (1985) in his paper "The Duration of Contract Strikes
   in U.S. Manufacturing," Journal of Econometrics, 28, 5 - 28. */ 


data strike;
  input dur eco;
datalines; 
  7.00000      .0113800
  9.00000      .0113800
  13.0000      .0113800
  14.0000      .0113800
  26.0000      .0113800
  29.0000      .0113800
  52.0000      .0113800
  130.000      .0113800
  9.00000      .0229900
  37.0000      .0229900
  41.0000      .0229900
  49.0000      .0229900
  52.0000      .0229900
  119.000      .0229900
  3.00000     -.0395700
  17.0000     -.0395700
  19.0000     -.0395700
  28.0000     -.0395700
  72.0000     -.0395700
  99.0000     -.0395700
  104.000     -.0395700
  114.000     -.0395700
  152.000     -.0395700
  153.000     -.0395700
  216.000     -.0395700
  15.0000     -.0546700
  61.0000     -.0546700
  98.0000     -.0546700
  2.00000      .00535000
  25.0000      .00535000
  85.0000      .00535000
  3.00000      .0742700
  10.0000      .0742700
  1.00000      .0645000
  2.00000      .0645000
  3.00000      .0645000
  3.00000      .0645000
  3.00000      .0645000
  4.00000      .0645000
  8.00000      .0645000
  11.0000      .0645000
  22.0000      .0645000
  23.0000      .0645000
  27.0000      .0645000
  32.0000      .0645000
  33.0000      .0645000
  35.0000      .0645000
  43.0000      .0645000
  43.0000      .0645000
  44.0000      .0645000
  100.000      .0645000
  5.00000     -.104430    
  49.0000     -.104430    
  2.00000     -.00700000
  12.0000     -.00700000
  12.0000     -.00700000
  21.0000     -.00700000
  21.0000     -.00700000
  27.0000     -.00700000
  38.0000     -.00700000
  42.0000     -.00700000
  117.000     -.00700000
;

data strike;
  set strike;
  dum = 1;

proc iml;
start mle;
use strike;
read all into t var{dur};
read all into x var{eco};
read all into d var{dum};

/* Calculation of Unrestricted MLE estimates using Newton-Raphson Method */

theta= {1,4,-9};
crit =1;
n=nrow(t);
ones = j(n,1,1);
result=j(10,9,0);

do iter=1 to 10 while (crit>1.0e-10);
 sigma=theta[1,1];
 beta1=theta[2,1];
 beta2=theta[3,1];
 w = (ones/sigma)#(log(t) - ones#beta1 - x#beta2);
 lnL = d#(w - log(sigma))- exp(w);
 lnL = sum(lnL);
 g1 = (ones/sigma)#(w#exp(w) - d#(w + ones));
 g2=(ones/sigma)#(exp(w) - d);
 g3 = (ones/sigma)#x#((exp(w) - d));
 g1=sum(g1);
 g2=sum(g2);
 g3=sum(g3);
 g=g1//g2//g3;
 h11= -(ones/sigma**2)#((w##2)#exp(w) + 2#w#exp(w) - 2#w#d - d);
 h11= sum(h11);
 h12= -(ones/sigma**2)#(exp(w) - d + w#exp(w));
 h12 = sum(h12);
 h13= -(ones/sigma**2)#x#(exp(w) - d + w#exp(w));
 h13 = sum(h13);
 h21 = h12;
 h31 = h13;
 h22 = -(ones/sigma**2)#exp(w);
 h22 = sum(h22);
 h23 = -(ones/sigma**2)#x#exp(w);
 h23 = sum(h23);
 h32 = h23;
 h33 = -(ones/sigma**2)#(x##2)#exp(w);
 h33 = sum(h33);
 h=(h11||h12||h13)//(h21||h22||h23)//(h31||h32||h33);
 db=-inv(h)*g;
 thetanew = theta + db;
 crit = sqrt(ssq(thetanew-theta));
 theta=thetanew;
 result[iter,] = iter||(theta`)||g1||g2||g3||crit||lnL;
 end;

 lnLu = lnL;

cnames = {iter,sigma,beta1,beta2,g1,g2,g3,crit,lnLu};

print "Calculation of Unrestricted MLE estimates using Hessian-Based Newton-Raphson Method";
print "Iteration steps ", result[colname=cnames];

print , "Unrestricted Log-likelihood = ", lnLu;

 /*Covariance matrix from Hessian*/

cov = -inv(h);
se_sigma_h = sqrt(cov[1,1]);
se_beta1_h = sqrt(cov[2,2]);
se_beta2_h = sqrt(cov[3,3]);
z_sigma_h = sigma/se_sigma_h;
z_beta1_h = beta1/se_beta1_h;
z_beta2_h = beta2/se_beta2_h;

/*Covariance matrix from BHHH*/

g1 = (ones/sigma)#(w#exp(w) - d#(w + ones));
g2=(ones/sigma)#(exp(w) - d);
g3 = (ones/sigma)#x#((exp(w) - d));
gmat = g1||g2||g3;
bhhh = gmat`*gmat;
covbh3 = inv(bhhh);
se_sigma_b=sqrt(covbh3[1,1]);
se_beta1_b=sqrt(covbh3[2,2]);
se_beta2_b=sqrt(covbh3[3,3]);
z_sigma_b = sigma/se_sigma_b;
z_beta1_b = beta1/se_beta1_b;
z_beta2_b = beta2/se_beta2_b;

pnames = {sigma,beta1, beta2};
print , "The Maximum Likelihood Estimates: Hessian-Based Newton-Raphson Iteration",
      theta [rowname=pnames];
print , "Asymptotic Covariance Matrix-From Hessian", cov
      [rowname=pnames colname=pnames];
print "Standard errors: ",se_sigma_h,se_beta1_h,se_beta2_h;
print ,"Asymptotic Covariance Matrix-From bhhh", covbh3
      [rowname=pnames colname=pnames];
print "Standard errors: ",se_sigma_b,se_beta1_b, se_beta2_b;

print "Wald test of hypothesis of constant hazard (sigma=1)";
   Wald = (sigma-1)*inv(cov[2,2])*(sigma-1);        * Wald test;
   critval = cinv(.95,1);      * calculates the 95th percentile of chi-square 1;
   pval = 1 - probchi(wald,1); * calculates the probability value of Wald;
   print "Results of Wald test Using Hessian" Wald critval pval;

   Wald = (sigma-1)*inv(covbh3[2,2])*(sigma-1);        * Wald test;
   critval = cinv(.95,1);      * calculates the 95th percentile of chi-square 1;
   pval = 1 - probchi(wald,1); * calculates the probability value of Wald;
   print "Results of Wald test Using BHHH" Wald critval pval;
   
/* ML Estimation of Restricted Model*/

   print , "Maximum Likelihood Estimation of Restricted Model";
     print "*************************************************";
   theta = {4,-9};
   crit = 1;
   n = nrow(t);
   result = j(10,7,0);
   
   do iter = 1 to 10 while (crit > 1.0e-10);

   beta1=theta[1,1];
   beta2=theta[2,1];
   w = (log(t) - ones#beta1 - x#beta2);
   lnLr = d#w - exp(w);
   lnLr = sum(lnLr);
   g1 = -(d - exp(w));
   g1 = sum(g1);
   g2 = -x#(d - exp(w));
   g2 = sum(g2);
   g = g1//g2;
   h11 = -exp(w);
   h12 = -x#exp(w);
   h22 = -(x##2)#exp(w);
   h11 = sum(h11);
   h12 = sum(h12);
   h21 = h12;
   h22 = sum(h22);
   h = (h11||h12)//(h21||h22);
   db = -inv(h)*g;
   thetanew = theta + db;
   crit = sqrt(ssq(thetanew - theta));
   result[iter,] = iter||(theta`)||g1||g2||crit||lnLr;
   theta = thetanew;
   end;

   cov = -inv(h);
   cnames = {iter,beta1,beta2,g1,g2,crit,lnLr};
   print "Iteration steps",result [colname=cnames];
   pnames = {beta1,beta2};
   print , "The Maximum Likelihood Estimates-Restricted Model", (theta`)
                      [colname=pnames];
   print , "Asymptotic Covariance Matrix-From Hessian of Restricted Model", cov
           [rowname=pnames colname=pnames];

   /* Gradient evaluated at restricted MLE estimates  */

   sigma = 1;
   w = (ones/sigma)#(log(t) - ones#beta1 - x#beta2);
   g1 = (ones/sigma)#(w#exp(w) - d#(w + ones));
   g2=(ones/sigma)#(exp(w) - d);
   g3 = (ones/sigma)#x#((exp(w) - d));
   gmat = g1||g2||g3;
   g1=sum(g1);
   g2=sum(g2);
   g3=sum(g3);
   g=g1//g2//g3;

/* Hessian evaluated at restricted MLE estimates */

 h11= -(ones/sigma**2)#((w##2)#exp(w) + 2#w#exp(w) - 2#w#d - d);
 h11= sum(h11);
 h12= -(ones/sigma**2)#(exp(w) - d + w#exp(w));
 h12 = sum(h12);
 h13= -(ones/sigma**2)#x#(exp(w) - d + w#exp(w));
 h13 = sum(h13);
 h21 = h12;
 h31 = h13;
 h22 = -(ones/sigma**2)#exp(w);
 h22 = sum(h22);
 h23 = -(ones/sigma**2)#x#exp(w);
 h23 = sum(h23);
 h32 = h23;
 h33 = -(ones/sigma**2)#(x##2)#exp(w);
 h33 = sum(h33);
 h=(h11||h12||h13)//(h21||h22||h23)//(h31||h32||h33);

   
   LM = g`*(-inv(h))*g;                         * LM test;
   critval = cinv(.95,1);
   pval = 1 - probchi(LM,1);
   print "Results of LM test Using Hessian" LM critval pval;

/* BHHH evaluated at Restricted MLE*/

   bhhh = gmat`*gmat;
   covbh3r = inv(bhhh);

   LM = g`*covbh3r*g;                         * LM test;
   critval = cinv(.95,1);
   pval = 1 - probchi(LM,1);
   print "Results of LM test Using BHHH" LM critval pval;



   LR = -2*(lnLr-lnLu);                       * Likelihood Ratio test;
   pval = 1 - probchi(LR,1);
   print "Results of LR test" LR critval pval;

/*  Let's see if we get essentially the same maximum likelihood estimates if we use
    a BHHH-based Newton-Raphson iteration.  */ 

   theta= {1,3.77,-9.35};
   crit =1;
   n=nrow(t);
   ones = j(n,1,1);
   result=j(60,9,0);

 do iter= 1 to 60 while (crit>1.0e-10);
 sigma=theta[1,1];
 beta1=theta[2,1];
 beta2=theta[3,1];
 w = (ones/sigma)#(log(t) - ones#beta1 - x#beta2);
 lnL = d#(w - log(sigma))- exp(w);
 lnL = sum(lnL);
 g1 = (ones/sigma)#(w#exp(w) - d#(w + ones));
 g2=(ones/sigma)#(exp(w) - d);
 g3 = (ones/sigma)#x#((exp(w) - d));
 gmat = g1||g2||g3;
 g1 = sum(g1);
 g2 = sum(g2);
 g3 = sum(g3);
 g = g1//g2//g3;
 bhhh = gmat`*gmat;
 db= inv(bhhh)*g;
 thetanew = theta + db;
 crit = sqrt(ssq(thetanew-theta));
 theta = thetanew;
 result[iter,] = iter||(theta`)||g1||g2||g3||crit||lnL;
 end;

cnames = {iter,sigma,beta1,beta2,g1,g2,g3,crit,lnL};

print "Calculation of Unrestricted MLE estimates using BHHH-Based Newton-Raphson Method";
print "Iteration steps ", result[colname=cnames]; 

finish;
run mle;

quit;

/*  Here we compare the results produced by the SAS Proc Lifereg.  */
    
proc lifereg data=strike;
  model dur = eco / dist = weibull;
    
 run;