? Newton iterations
proc MLNewt maxit tol;
local grad vcov delt;

const maxsqz,30;
supres @smpl;

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

const ifhess,1;
mform(nrow=novar,type=sym) vcov;

dologl loglold;  ? LogL at initial values
const iter,0;
write iter,loglold;
const nfeval,1;

do iter=1,maxit;

?  derivatives of LogL
  derivs;

  ?smpl smpall;
  smpl smpw;
  msd(noprint) dnames;
  copy @sum grad;
  msd(noprint) dl2;
  set js = 1;
  do j=1,novar;
    do i=j,novar;
      set vcov(i,j) = @sum(js);
      set js = js+1;
    enddo;
  enddo;
  mat vcov = vcov";
  mat delt = vcov*grad;
  mat crit = grad'delt;

  ? test for convergence and squeeze to determine stepsize.

  mat relchg = abs(delt)/(abs(b) + 10*tol);
  mat worst = max(relchg);
  if worst<tol; then; goto 20;

  set lambda = 1;
  set lamfac = .5;
  do isqz=1,maxsqz;
    mat testb = b + delt*lambda;
    unmake testb bnames;

    ?----------------------------
    ? evaluate log likelihood here
    ?----------------------------

    dologl loglnew;
    if (.not.miss(loglnew)); then; set nfeval = nfeval+1;

    ?----------------------------
    ?write loglold loglnew isqz;
    ?write loglold loglnew lambda lamfac isqz;

    ? plain cea method
    ?if loglnew>loglold; then; goto 15;

    ? implement ceab below
    if (miss(loglnew) .or. (loglnew<=loglold)) .and. (lamfac=2); then; do;
      set lambda = lambda*.5;
      mat testb = b + delt*lambda;
      unmake testb bnames;
      dologl loglnew;
      set nfeval = nfeval+1;
      goto 15;
    enddo;
    if loglnew>loglold; then; do;
      if (lamfac=.5) .and. (isqz>1); then; goto 15;
      if (lamfac=.5) .and. (isqz=1); then; set lamfac=2;  ? ceab kicks in
      set loglold = loglnew;
    enddo;

    set lambda = lambda*lamfac;
  enddo;
  if lamfac=2; then; goto 15;
  write(form="(' Failure to improve after ',F4.0,' squeezes.')") isqz;
  write relchg b delt;
  goto 18;
 
?15 write iter,loglnew,isqz;
15 write iter,loglnew,lambda crit;
  set loglold = loglnew;
  copy testb b;
enddo;

18 write(form="(' Failure to converge after ',F4.0,' iterations.')") iter;
write relchg b delt;
goto 30;

20 write(form="(' Convergence after ',F4.0,' iterations.')") iter;

30 unmake b bnames;
write nfeval loglold bnames crit;
title 'Newton SEs';
tstats(names=bnames) b vcov;

mmake gradi dnames;
mat bhhh = gradi'gradi;
title 'White SEs';
mat vcovw = vcov'bhhh*vcov;
tstats(names=bnames) b vcovw;

title 'BHHH SEs';
mat vcov = bhhh";
tstats(names=bnames) b vcov;

endproc;