/*************   Newey-West (Bartolini & Kramer, Jun95) *******************/
/*                                                                        */
/*                       NEWEY - WEST errors                              */
/*                                                                        */
/**************************************************************************/
/*
                 authors: Leonardo Bartolini
                          Charles Kramer

                          Research Department, IMF
                          700, 19th St. NW,
                          Washington, D.C. 20431, USA

                 correspondence to: LBARTOLINI@IMF.ORG
/**************************************************************************/

Notes:

 The routine is written as a procedure NW(y,X,b,lag);
 It arguments are:

       - a vector y[n,1] of dependent variables;
       - a matrix X[n,k] of independent variables;
       - a vector b[k,1] of the associated OLS coefficients,
       - an integer lag-truncation ("lag");

 The procedures returns a matrix (k x 2), with

       - in the first column the Newey-West errors;
       - in the second column the OLS errors.

 The program is written following the treatment and notation of
 Gallant, R. (1987), "Nonlinear Statistical Models," pp.137-139.

 This program is for public, non commercial use.
 It was succesfully tested against the NW errors provided by TSP.
 Nevertheless, the authors disclaim any responsibility for its use.

                                                                     */
/*********************************************************************/

PROC NW(y,X,b,lag);

    LOCAL sse,n,yhat,e,G,w,a,t,ga,V,F,nwerr,olserr,k,za,hhat;

    n=ROWS(X); k=ROWS(B);
    yhat=X*b; e=y-yhat;
    hhat=e'.*x';
    G=ZEROS(k,k); w=ZEROS(2*lag+1,1);

    a=0;

    DO UNTIL a==lag+1;
        ga=ZEROS(ROWS(b),ROWS(b));
        w[lag+1+a]=(lag+1-a)/(lag+1);
        za=hhat[.,(a+1):n]*hhat[.,1:n-a]';
        IF a==0; ga=ga+za; ELSE; ga=ga+za+za'; ENDIF;
        G=G+w[lag+1+a]*ga;
        a=a+1;
    ENDO;

        F=X'*X;
        V=INV(F)*G*INV(F);
        nwerr= (DIAG(V))^.5;
        olserr=(DIAG(INV(X'X)*e'e/(n-k)))^.5;
    RETP(nwerr~olserr);

ENDP;


