/* proc for OLS */ 

proc(4)=OLS(Y,X);
local b_ols,e_ols,u_ols,s2,V0,V,aic,sq,q,hac_v;
X=ones(rows(X),1)~X;
b_ols=invpd(X'X)*X'Y;
e_ols=Y-X*b_ols;
u_ols=X.*e_ols;
s2=(e_ols'e_ols)/(rows(X)-cols(Y));
sq=(e_ols'e_ols)/rows(X);
aic=ln(sq)+2*cols(X)/rows(X);
V0=s2*invpd(X'X);
V=invpd(X'X)*(u_ols'u_ols)*invpd(X'X);
q=round(4*(rows(X)/100)^(2/9));
{hac_v}=NWHAC(Y,X,b_ols,q);


retp(b_ols,V0,V,hac_v);
endp;


proc(1)=NWHAC(y,X,b,q);  // Modified from Bartolini & Kramer (1995)

    local n,yhat,e,G,w,j,t,ga,VHAC,F,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*q+1,1);

    j=0;

    do while j < q+1;
        ga=ZEROS(ROWS(b),ROWS(b));
        w[q+1+j]=(q+1-j)/(q+1);
        za=hhat[.,(j+1):n]*hhat[.,1:n-j]';
        if j==0; ga=ga+za; else; ga=ga+za+za'; endif;
        G=G+w[q+1+j]*ga;
        j=j+1;
    ENDO;

        F=X'X;
        VHAC=INVPD(F)*G*INVPD(F);
        
    retp(VHAC);

ENDP;