new; cls; /* BHHH method */
x={4.2,4.8,6.0,3.7,4.5,3.2,5.3,4.9,4.8,6.8};
y={5.4,5.9,5.9,6.7,3.7,2.3,4.2,3.2,7.3,6.6};
x1=ones(rows(x),1)~x;
b0={1,1,1};
print "b*=" bhhh(&ll,b0)';


proc ll(b);
   local beta,s,e;
   beta=b[1:2];
   s=b[3];
   e=y-x1*beta;
   retp( -( -1/2*ln(2*pi)-1/2*ln(s^2)-1/2*(e/s)^2 ) );
endp;

proc bhhh(&ll,beta);
   local step,maxiter,tol,i,g,g_1,h,ll:proc;
   step=0.5; maxiter=1e+7; tol=1e-5; g=1;
   i=1; 
   do while i<=maxiter;
      g_1=g;
      g=gradp(&ll,beta);
      H=g'g;
      if abs(g-g_1)<tol;
         break;
      endif;
      beta=beta-step*sumc(g)/H;
      i=i+1;
   endo;
   retp(beta);
endp;