/*--------------------------------------------------------
!! AUTOLAG.SET : COMPUTE SPECTRAL DENSITY MATRIX AT FREQUENCY
                 ZERO OF A VECTOR OF ORTHOGONALITY CONDITIONS
                 BY NEWEY AND WEST (1987 ECONOMETRICA) WITH THEIR
                 (1994 R.E. STUD) AUTOMATIC LAG SELECTION PROCEDURE

   USAGE:    {sw,nwlag,n}=autolag(h,prewite);

   OUTPUT:   sw:   spectral density matrix at freq. 0 of vector of
                   orthogonality conditions. (inv(sw) is the optimal
                   weighting matrix in gmm).
          nwlag:   bandwidth parameter.
              n:   lag length to compute s0 and s1.

   INPUT:     h:   T X K matrix of orthogonality conditions.
        prewite:   1 asks for prewhitening. 0--no prewhitening.
-------------------------------------------------------------------*/

proc(3) = autolag(h,prewite);

local wvec,h,hstar,htil,A,v,T1,k,n,sto,gam,m,
      sigma0,sigma,s0,s1,sbar,vcv,i,j,t,w,s,ht,sto_1,sto_2;

 k=cols(h);                    T1=rows(h);

 wvec=ones(k,1);               wvec[1]=0;

 hstar=h[1:T1-1,.];            htil=h[2:T1,.];


 sto_1=zeros(k,k); sto_2=zeros(k,k);     @ FITTING THE VAR(1)    @
 if prewite==1;
  t=2; do while t<=T1;
   sto_1=sto_1+h[t-1,.]'*h[t-1,.];
   sto_2=sto_2+h[t,.]'*h[t-1,.];
  t=t+1; endo;
  A=(sto_2)*inv(sto_1);
  v=htil-hstar*A';
 else;
  v=htil;
 endif;

 v=zeros(1,k)|v;

    n=int(4*(T1*.01)^(2/9));  
/*  n=int(12*(T1*.01)^(2/9)); */


wee:
 sigma=zeros(n,1); sigma0=0;
 
 sto=v*wvec;
 sigma0=sto[2:T1]'*sto[2:T1];
 sigma0=sigma0/(T1-1);


 j=1; do while j<=n;
 sigma[j]=sto[j+2:T1]'*sto[2:T1-j]/(T1-1);  
 j=j+1; endo;

 clear sto,s0,s1;
 j=1; do while j<=n;
  s1=s1+2*j*sigma[j];
  sto=sto+2*sigma[j];
 j=j+1; endo;

 s0=sigma0+sto;

 gam=1.1447*(((s1/s0)^2)^(1/3));

 m=int(gam*(T1^(1/3)));   @ Bandwidth parameter @

woe:
 w=zeros(k,k) ;
  t=2;
  do while t<=T1;
   w = v[t,.]'*v[t,.] + w ;
  t=t+1;  endo;
 sbar=w/(T1-1);

  j=1 ;
  do while j<=m;
   w=zeros(k,k) ;
   t=j+2;
   do while t<=T1;
    w = v[t,.]'*v[t-j,.] + w ;
   t=t+1;  endo;
   w=w/(T1-1);

   sbar = sbar + (w + w')*(1-(j/(m+1)));

  j=j+1;  endo;

  if prewite==1;  
   sbar=inv(eye(rows(A))-A)*sbar*(inv(eye(rows(A))-A))';
  endif;

  retp(sbar,m,n) ;
  endp ;

@=============================================================@