?  Hessian Check code
?  Requires Procs  dologl  and  derivs  .

Proc HESSCHEC;

supres @smpl;

mmake b bnames;
mat novar = nrow(b);

? compute numeric second derivatives, by varying first derivatives
const ifhess,0;
NUMHESS numeric;

?  analytic derivatives of LogL
dologl loglold;
const ifhess,1;
derivs;

smpl smpw;
msd(noprint) dl2;
copy numeric analytic;
set js = 1;
do j=1,novar;
  do i=j,novar;
    set analytic(i,j) = @sum(js);
    set js = js+1;
  enddo;
enddo;

? print comparison of analytic and numeric
mat error = (analytic-numeric)/analytic;
print analytic numeric error;

endproc;

?--------------------------------------
? Proc for numeric derivatives
Proc NUMHESS hess;
local b novar bold bi eps1 eps2 eps gradp gradm ii jj loglp loglm;

mmake b bnames;
mat novar = nrow(b);

copy b bold;
copy b gradp;
copy b gradm;
mform(nrow=novar,type=sym) hess;
do ii=1,novar;
  copy bold b;
  set bi = b(ii);
  ?set eps1 = 1.e-4;
  set eps1 = 1.e-6;
  ?set eps2 = abs(bi*1.e-3);
  set eps2 = abs(bi*1.e-5);
  set eps = eps1 + pos(eps2-eps1);  ? eps = max(eps1,eps2)
  set b(ii) = bi+eps;
  unmake b bnames;
  dologl loglp;  ? evaluate h0 and u at new param values
  derivs;
    msd(noprint) dnames;
    copy @sum gradp;
  set b(ii) = bi-eps;
  unmake b bnames;
  dologl loglm;
  derivs;
    msd(noprint) dnames;
    copy @sum gradm;
  do jj=1,ii;
    set hess(ii,jj) = (gradm(jj)-gradp(jj))/(2*eps);
  enddo;
  set b(ii) = bi;
  unmake b bnames;
enddo;

endproc;