Two-Equation Maximum Likelihood Problem

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: IMNLPEX7                                            */
/*   TITLE: Two-Equation Maximum Likelihood Problem             */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswmh/sasghg               UPDATE: frwick 5/2015   */
/*     REF: Technical Report Examples                           */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

proc iml;
z={ 1.33135 0.64629 0.4026 -20 0.24447,
    1.39235 0.66302 0.4084 -19 0.23454,
    1.41640 0.65272 0.4223 -18 0.23206,
    1.48773 0.67318 0.4389 -17 0.22291,
    1.51015 0.67720 0.4605 -16 0.22487,
    1.43385 0.65175 0.4445 -15 0.21879,
    1.48188 0.65570 0.4387 -14 0.23203,
    1.67115 0.71417 0.4999 -13 0.23828,
    1.71327 0.77524 0.5264 -12 0.26571,
    1.76412 0.79465 0.5793 -11 0.23410,
    1.76869 0.71607 0.5492 -10 0.22181,
    1.80776 0.70068 0.5052  -9 0.18157,
    1.54947 0.60764 0.4679  -8 0.22931,
    1.66933 0.67041 0.5283  -7 0.20595,
    1.93377 0.74091 0.5994  -6 0.19472,
    1.95460 0.71336 0.5964  -5 0.17981,
    2.11198 0.75159 0.6554  -4 0.18010,
    2.26266 0.78838 0.6851  -3 0.16933,
    2.33228 0.79600 0.6933  -2 0.16279,
    2.43980 0.80788 0.7061  -1 0.16906,
    2.58714 0.84547 0.7567   0 0.16239,
    2.54865 0.77232 0.6796   1 0.16103,
    2.26042 0.67880 0.6136   2 0.14456,
    1.91974 0.58529 0.5145   3 0.20079,
    1.80000 0.58065 0.5046   4 0.18307,
    1.86020 0.62007 0.5711   5 0.18352,
    1.88201 0.65575 0.6184   6 0.18847,
    1.97018 0.72433 0.7113   7 0.20415,
    2.08232 0.76838 0.7461   8 0.18847,
    1.94062 0.69806 0.6981   9 0.17800,
    1.98646 0.74679 0.7722  10 0.19979,
    2.07987 0.79083 0.8557  11 0.21115,
    2.28232 0.88462 0.9925  12 0.23453,
    2.52779 0.95750 1.0877  13 0.20937,
    2.62747 1.00285 1.1834  14 0.19843,
    2.61235 0.99329 1.2565  15 0.18898,
    2.52320 0.94857 1.2293  16 0.17203,
    2.44632 0.97853 1.1889  17 0.18140,
    2.56478 1.02591 1.2249  18 0.19431,
    2.64588 1.03760 1.2669  19 0.19492,
    2.69105 0.99669 1.2708  20 0.17912 };

start fiml(pr) global(z);
   c1 = pr[1]; c2 = pr[2]; c3 = pr[3]; c4 = pr[4]; c5 = pr[5];
   /* 1. Compute Jacobian */
   lndet = 0 ;
   do t= 1 to 41;
      j11 = (-c3/c4) * c1 * 10 ##(c2 * z[t,4]) * (-c4) * c5 *
            z[t,1]##(-c4-1) * (c5 * z[t,1]##(-c4) + (1-c5) *
            z[t,2]##(-c4))##(-c3/c4 -1);
      j12 = (-c3/c4) * (-c4) * c1 * 10 ##(c2 * z[t,4]) * (1-c5) *
            z[t,2]##(-c4-1) * (c5 * z[t,1]##(-c4) + (1-c5) *
            z[t,2]##(-c4))##(-c3/c4 -1);
      j21 = (-1-c4)*(c5/(1-c5))*z[t,1]##( -2-c4)/ (z[t,2]##(-1-c4));
      j22 = (1+c4)*(c5/(1-c5))*z[t,1]##( -1-c4)/ (z[t,2]##(-c4));

      j = (j11 || j12 ) // (j21 || j22) ;
      if any(j = .) then detj = 0.;
         else detj = det(j);
      if abs(detj) < 1.e-30 then do;
         return(.);
      end;
      lndet = lndet + log(abs(detj));
   end;

   /* 2. Compute Sigma */
   sb = j(2,2,0.);
   do t= 1 to 41;
      eq_g1 = c1 * 10##(c2 * z[t,4]) * (c5*z[t,1]##(-c4)
              + (1-c5)*z[t,2]##(-c4))##(-c3/c4) - z[t,3];
      eq_g2 = (c5/(1-c5)) * (z[t,1] / z[t,2])##(-1-c4) - z[t,5];
      resid = eq_g1 // eq_g2;
      sb = sb + resid * resid`;
   end;
   sb = sb / 41;
   /* 3. Compute log L */
   const = 41. * (log(2 * 3.1415) + 1.);
   lnds = 0.5 * 41 * log(det(sb));
   logl = const - lndet + lnds;
   return(logl);
finish fiml;

pr = j(1,5,0.001);
optn = {0 2};
tc = {. . . 0};
call nlpdd(rc, xr,"fiml", pr, optn,,tc);
quit;