% function W = nw(ZV, T)    
% 
% L=ceil(1.3*T^(1/2)); % Newey-West (1987) with truncation parameter recommended by Lazarus, Lewis, Stock, Watson (JBES 2018)
%     for j=0:L
%         if j>0
%             acv(:,:,j+1)=ZV(j+1:end,:)'*ZV(1:end-j,:)/T + ZV(1:end-j,:)'*ZV(j+1:end,:)/T;
%             w_l=1-j/(L);
%             W=W+w_l*acv(:,:,j+1);
%         else
%             acv(:,:,j+1)=ZV'*ZV/T;
%             W=acv(:,:,j+1);
%         end
%     end
% 
% end

function W = nw(ZV, T)
% Compute the Newey-West variance-covariance matrix for ZV (T x K)

    L = ceil(1.3 * sqrt(T));  % Truncation lag, Lazarus et al. (2018)
    [T, K] = size(ZV);        % dimensions of ZV
    W = zeros(K, K);          % initialize variance-covariance matrix

    for j = 0:L
        weight = 1 - j / (L + 1);  % Bartlett weight (classic NW)

        if j == 0
            Gamma0 = (ZV' * ZV) / T;
            W = W + Gamma0;
        else
            Gamma_j = (ZV((j+1):end, :)' * ZV(1:end-j, :)) / T;
            Gamma_j_sym = Gamma_j + Gamma_j';  % symmetrize
            W = W + weight * Gamma_j_sym;
        end
    end
end
