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;