/* --------------      Newey and West 1994 (Wong, Jun95)   --------------- **
**  Newey and West (1994): "Automatic Lag Selection in Covariance Matrix   **
**        Estimation," Review of Economic Studies, 64(4), pp.631-653.      **
** ----------------------------------------------------------------------- **

      Author: Ka-fu Wong                   email: kfwong@macc.wisc.edu
              Department of Economics,
              University of Wisconsin,
              Madison, WI 53706

    This code is written and submitted for public, non-commerical use.
    You are welcome to use and distribute them as long as you include proper
    attribution to the author.  Although this code has been tested
    extensively, we do not guarantee its performance.  Comments and
    suggestions are welcome.  It is last revised on June 12, 1995.
** ----------------------------------------------------------------------- **

The following procedures compute the covariance matrix nonparametrically.
Choices of kernels include Bartlett, Truncated and Quadratic Spectral.
See Newey and West(1994) for details.  A separate file EXAMPLE.PGM
demonstrates the use of these procedures.

Format: {var,gam}=BTBAR(m,vstar,ahat,ww)      for Bartlett
        {var,gam}=QSBAR(vstar,ahat,ww)        for Quadratic Spectral
        {var,gam}=TRU(m,vstar,ahat)           for Truncated

*
* INPUT:
*       m:= Bandwidth for Truncated and Bartlett.
*   vstar:= residuals used to compute the covariance matrix.
*    ahat:= the (rxr) matrix of VAR(1) regression coefficients obtained
*           regressing cross products of instruments and residuals on
*           their first lags.
*     ww := the weighting vector to reduce vstar to a scalar.
*           In general set ww=(0,1,1,...,1).  I.e. equal
*           weight to all regressors except the constant term (assuming that
*           the first regressor is a constant term).
*
* OUTPUT:
*   var:= covariance matrix.
*   gam:= same as the gamma for Bartlett and QS kernels in Newey and West(1994);
*         indicator for non-positive definiteness of the covariance matrix
*         for Truncated kernel for the specified bandwidth.  (When
*         this happens, the bandwidth is set to 1.)
*
-------------------------------------------------------------------------- */
/*
** Parameters for Bartlett and QS kernels
*/
   btq=1; btcgam=1.1447; qsq=2; qscgam=1.3221; @ c.f. Newey & West(1994) @



/*
** Kernel functions
*/
   @ ---- Truncated Kernel ---- @
     fn truker(x)=1;
   @ ---- Quadratic Spectral Kernel ---- @
     fn qsker(x)= ( 25/( 12*(pi^2)*((x[1,.]/x[2,.])^2) ) )
          * ( sin(6*pi*(x[1,.]/x[2,.])/5) / (6*pi*(x[1,.]/x[2,.])/5)
            - cos(6*PI*(x[1,.]/x[2,.])/5)                            ) ;
   @ ---- Bartlett Kernel ---- @
     fn btker(x)=1-(x[1,.]/(x[2,.]+1)) ;




/*
** Compute small s (for BT only)
** See Newey and West(1994) for details.
*/
proc (1) = ssmall(Vstar,p,rv,ww); @ only for Bartlett kernel @
local hh,sig,j,s,s1,s2,tempsig;

   hh=Vstar*ww; @ ##### h=(xx*w).*e;@

   /* generates sig ( p+1 x 1 ) */
   sig=(hh'hh)/rv;
   j=1; do until j>p;
            /* it is important not to adjust for degrees of freedom */
      tempsig= ( hh[1:(rv-j),.]'hh[(j+1):rv,.] ) / rv ;
      sig= sig | tempsig ;
   j=j+1; endo;

   s = 2 * sumc(sig) - sig[1,.];
   s1= 2 * sumc(   (0|seqa(1,1,p))        .* sig ) ;
   s2= 2 * sumc( ( (0|seqa(1,1,p)) .^ 2 ) .* sig ) ; @ used only by Parzen @

retp(s|s1|s2);
endp;




/*
** Compute small s for QS to be used to compute m
** See Newey and West(1994) for details.
*/
proc (1) = qsssmall(hh,rv,ww); @ only for QS kernel @
local j,s,s1,s2,pres2,pres0,be,uhat,usig2,rv_1,xx,yy;

   pres2={}; pres0={}; rv_1=rv-1;

   j=1; do until j>cols(hh);
       xx=hh[1:(rv_1),j]; yy=hh[2:rv,j];
       be=invpd(xx'xx)*xx'yy;
       uhat=yy-xx*be;
       usig2=(uhat'uhat)/(rv_1);  @ 1 obs lost due to AR(1) @
       pres2=pres2|((2*be*usig2/((1-be)^4))^2);
       pres0=pres0|((usig2/(1-be)^2)^2);
    j=j+1; endo;

    s2=ww'pres2;  s=ww'pres0;  s1=0;  @ s1, a dummy, never referenced @

retp(s|s1|s2);
endp;




/*
** Compute covariance matrix (for all kernels)
*/
proc (1) = cov(hh,&ker,m,mm,rv);
local ker:fn,scap,j,omej;

   scap=(hh'hh)/rv;
   j=1; do until j>mm;
      omej = ( ( hh[1:(rv-j),.]'hh[(j+1):rv,.]) ./ rv );
      scap = scap + ker(j|m) .* ( omej + omej') ;
   j=j+1; endo;

retp(scap);
endp;




/*
** Compute Jhat using the prewhiten h, then recolor
** See Newey and West(1994) and Andrews and Monahan(1992) for details.
** set ahat=0 if prewhitening has not been done.
*/
proc (1) = jpw(vstar,ahat,&ker,m,mm,rv);
local ker:fn,jhatstar,jhatpw,dhat;

   dhat=inv(eye(cols(ahat))-ahat);
   jhatstar=cov(Vstar,&ker,m,mm,rv);
   jhatpw=dhat*jhatstar*dhat';       @ recolor the estimate @

retp(jhatpw);
endp;




/*
** Truncated Procedure
*/
proc(2) = tru(m,vstar,ahat);
local mm,gam,jhatpw,rv,oldtrap,va;

   rv=rows(vstar);   gam=0;   mm=m;
   jhatpw=jpw(vstar,ahat,&truker,m,mm,rv);

    /* IF Jhatpw is not positive definite, set m=0 and recalculate Jhatpw. */
      oldtrap=trapchk(0xffff);
      trap 1;
         va=eigh(jhatpw);
      if (scalerr(va)/=0) or (minc(va)<=0);
         m=0; mm=m; gam=1;
         jhatpw=jpw(vstar,ahat,&truker,m,mm,rv);
      endif;
      trap oldtrap,0xffff;

retp(jhatpw,gam);
endp;




/*
** Quadratic Spectral Procedure
*/
proc (2) = qsbar(vstar,ahat,ww);
local s,gam,m,mm,jhatstar,jhatpw,rv;

   rv=rows(vstar);

   s=qsssmall(vstar,rv,ww);
   gam = qscgam * ( ( (s[qsq+1,.]/s[1,.])^2 ) ^ (1/(2*qsq+1)) ) ;
   m= gam * (rv^(1/(2*qsq+1))) ;
   mm=rv-1;

   jhatpw=jpw(vstar,ahat,&qsker,m,mm,rv);

retp(jhatpw,gam);
endp;




/*
** Bartlett Procedure
*/
proc (2) = btbar(p,vstar,ahat,ww);
local s,gam,m,mm,Jhatstar,Jhatpw,rv;

   rv=rows(vstar);

   s=ssmall(vstar,p,rv,ww);
   gam = btcgam * ( ( (s[btq+1,.]/s[1,.])^2 ) ^ (1/(2*btq+1)) ) ;
   m=trunc( gam * (rv^(1/(2*btq+1))) ) ;
   mm=(m>(rv-1))*(rv-1-m)+m;

   jhatpw=jpw(vstar,ahat,&btker,m,mm,rv);

retp(jhatpw,gam);
endp;




/* -----------------------------------------------------------------------
procedure: pwhite1
Purpose: computes VAR (order 1) coefficients and residuals.  It
         uses singular value decomposition to restrict the resulting
         VAR coefficients.  This procedure matches the prewhitening
         procedure in Andrews and Monahan(1992).

*
* INPUT:
*   singpar:= 2x1 vector to restrict the VAR parameters using singular value
*             decomposition.  Set singpar=0 if this restriction is not needed.
*      hcap:= TxK matrix.
*
* OUTPUT:
*    ahat:= KxK matrix, coefficients from VAR(1);
*   vstar:= (T-1)xK matrix, residuals from VAR(1);
--------------------------------------------------------------------------- */
proc (2) = pwhite1(singpar,hcap);
local xxx,yyy,ahat,vstar,j,rowhcap,colhcap,bhat,s,chat,
      deltahat,singpar1,singpar2;

  /**
   ** Prewhiten hcap with VAR(1), see Andrews and Monahan(1992)
   **/
         rowhcap=rows(hcap);
         colhcap=cols(hcap);
         xxx=hcap[1:rowhcap-1,.];

    /** Step 1: Compute VAR(1) coefficients **/
     ahat={};
     j=1; do until j>colhcap;
         yyy=hcap[2:rowhcap,j];
         ahat=ahat~(invpd(xxx'xxx)*xxx'yyy);
     j=j+1; endo;

    /** Step 2: Use singular value decomposition to restrict parameters **/
     if rows(singpar)==2;
         singpar1=minc(singpar); singpar2=maxc(singpar);
         {bhat,s,chat}=svd1(ahat');
         deltahat=substute(s,(s .>  singpar2), singpar2);
         deltahat=substute(deltahat,(deltahat .< singpar1), singpar1);
         ahat=bhat*deltahat*chat';
     endif;

    /** Step 3: Compute VAR(1) residuals using the _new_ Ahat **/
     vstar={};
     j=1; do until j>colhcap;
         yyy=hcap[2:rowhcap,j];
         vstar=vstar~(yyy-xxx*ahat[j,.]');
     j=j+1; endo;

retp(ahat,vstar);
endp;