/* This functions calculates autocorrelation estimates for lag k */
proc autocor(series, k);
  local rowz,y,x,rho;
  rowz = rows(series);
  y = series[k+1:rowz];
  x = series[1:rowz-k];
  rho = inv(x'x)*x'y;            /* compute autocorrelation by OLS */
  retp(rho);
endp;

BHHH usually refers to the procedure explained in Berndt, E., Hall, B.,
Hall, R., & Hausman, J. (1974), "Estimation and Inference in Nonlinear
Structural Models," Annals of Economic and Social Measurement, 3/4: 653-665.

BHHH provides a method of estimating the asymptotic covariance matrix of a
Maximum Likelihood Estimator. In particular, the covariance matrix for a MLE
depends on the second derivatives of the log-likelihood function. However,
the second derivatives tend to be complicated nonlinear functions. BHHH
estimates the asymptotic covariance matrix using first derivatives instead
of analytic second derivatives. Thus, BHHH is usually easier to compute than
other methods.

In addition to the original BHHH article referenced above, BHHH is also
discussed in Greene, W.H., Econometric Analysis, 3rd Edition, Prentice-Hall,
1997. Greene's econometric software program, LIMDEP, uses BHHH for some of
the estimation routines.

Someone (perhaps BHHH themselves?) wrote a Fortran subroutine in the 1970's
to do BHHH. I do not have a copy of this subroutine at the present time. You
may want to check out Green's econometric software, LIMDEP, to see if it
will do what you require, rather than writing your own program to use an
existing BHHH subroutine. The Web address for LIMDEP is:
  http://www.limdep.com/index.htm

Cheers,
Paul.

--
  Paul L. Schumann, Ph.D., Professor of Management
  Minnesota State University, Mankato (formerly Mankato State University)
  Mankato, MN  56002
  mailto:paul.schumann@mankato.msus.edu
  http://krypton.mankato.msus.edu/~schumann/www/welcome.html

/* A series of residuals eps_hat[] is generated from a regression, e.g.: */

eps_hat = y - X*betaols;

/* Then the Box-Pierce statistic for each k can be calculated this way: */

print "Box-Pierce statistic for k=1 is" BP(eps_hat,1);
print "Box-Pierce statistic for k=2 is" BP(eps_hat,2);
print "Box-Pierce statistic for k=3 is" BP(eps_hat,3);

proc BP(series, k);
  local beep, rho;
  beep = 0;
  do until k < 1;
    rho = autocor(series, k);
    beep = beep + rho * rho;
    k = k - 1;
  endo;
  beep = beep * rows(series);     /* BP = T* (the sum) */
  retp(beep);
endp;

/* This functions calculates autocorrelation estimates for lag k */
proc autocor(series, k);
  local rowz,y,x,rho;
  rowz = rows(series);
  y = series[k+1:rowz];
  x = series[1:rowz-k];
  rho = inv(x'x)*x'y;            /* compute autocorrelation by OLS */
  retp(rho);
endp;