/*****************************************************************
******************************************************************
Proc: CORC
Author: Alan G. Isaac   aisaac@american.edu
Last Update: August 1996
Format: bt = corc(y,x,rho);
Purpose: computes iterated Cochrane-Orcutt regression estimates
         and t-ratios.
Remarks: A constant is _added_.
         There is _no_ checking for missing values.
         User may wish to change tol or max_iter.
Inputs:
  y - the LHS variable as a nobs x 1 vector,
      (first element must be the first observation)
  x - the RHS variables as a nobs x K matrix
      (no constant: a constant will be *added*)
  rho - an initial value of rho (usu. zero)
Output:
  bt - (K+2)x2 matrix, b~t, the parameter estimates b with
       associated t-ratios t
       (constant first,rho last, the rest ordered as is x)
References:
  Hamilton, J., Time Series Analysis, 1994, p.223
  Davidson R. & J. MacKinnon, 1993, Estimation and Inference
    in Econometrics, p.339
******************************************************************
*****************************************************************/
proc (1) = corc(y,x,rho);
  local tol,max_iter,iter,x_,x_1,y_,y_1,oldrho,
        nobs,b,xx,xy,xr,yr,u,u_1,e,df,ssr,sst,ee,t,invxx,cov;
  tol=10^-8; @convergence criterion@
  max_iter=25; @maximum allowable iterations@
  x = ones(rows(x),1)~x; @add constant term@
  x_=trimr(x,1,0);y_=trimr(y,1,0); @transform data for use in do loop@
  x_1=trimr(lag(x),1,0);y_1=trimr(lag(y),1,0);
  oldrho=rho+1; @so test criterion satisfied on first iteration@
  iter=0;
  do while abs(rho-oldrho)>tol;
    iter=iter+1;
    if iter>max_iter;
     "     Did not converge in "$+ftocv(max_iter,1,0)$+" iterations.
     Exit without convergence. Consider increasing max_iter.";stop;
    endif;
    oldrho=rho;  @store value of rho at beginning of interation@
    xr=x_-rho*x_1;
    yr=y_-rho*y_1;
    xx=xr'xr;xy=xr'yr;
    invxx=invpd(xx);
    b=invxx*xy;  @OLS regression of yr on xr@
    u=y-x*b;
    u_1=trimr(lag1(u),1,0);u=trimr(u,1,0);
    rho=(u_1'u)/(u_1'u_1);
  endo;
  nobs=rows(yr); @rows(y)-1@
  e=yr-xr*b;
  @Davidson & MacKinnon covariance matrix calculation:@
  xx=(xr~u_1)'(xr~u_1);
  invxx=invpd(xx);
  cov = (e'e/(nobs-cols(invxx))) * invxx;
  b=b|rho;
  t = b./sqrt(diag(cov));
  retp(b~t);
endp; @corc()@