Optimizing code for HAC Weight matrix generation - How to be faster?

Hi everyone,

I am re-creating Newey-West procedure for a heteroskedasticity and autocorrelation consistent standard errors (HAC) from scratch.

In order to compute the sandwich matrix V(b), I need to generate the weight matrix W

Since I have to run that test thousands of times, I need to optimize the code. Right now, I have this:

%% Weight matrix generation for HAC

T = length(Dependant_Variable);
k = length(Coeff);
df = T-k;

Optimal_Lags = floor(4*(T/100)^(2/9));

Residual_correlation = Residual*Residual';
Weight_correlation = zeros(size(Residual_correlation));

for i = 1:size(Weight_correlation,1)
    for j = 1:size(Weight_correlation,2)
        if abs(i-j)>Optimal_Lags
            Weight_correlation(i,j)=0;
        else
            Weight_correlation(i,j)=(T/df)*Residual_correlation(i,j)*(1-abs(j-i)/(Optimal_Lags+1));
        end
    end
end

Question:

How can I get the same result without using loops?

Thanks and best regards,

Answer 1

[Alan Stevens]

Try replacing the loops with something like this:

i = (1:size(Weight_correlation,1))';
j = 1:size(Weight_correlation,2);
c = ones(1,size(Weight_correlation,2));
r = ones(size(Weight_correlation,1),1);
lags = abs(i*c-r*j);
ind1 = lags>Optimal_Lags;
ind2 = lags<=Optimal_Lags;
Weight_correlation(ind1)=0;
Weight_correlation(ind2)=(T/df)*Residual_correlation(ind2).*(1-lags(ind2)/(Optimal_Lags+1));

Check carefully, as I can't really test it properly, not having your data.

[Vic]

You are right. This does not work.

close all; clearvars;clc;

Variables = rand(10,10);
Companies = rand(10,1);

Coeff = (Variables'*Variables)\Variables'*Companies;
Expected = Variables*Coeff;
Residual = Companies-Expected;

T = length(Companies);
k = length(Coeff);
df = T-k;

Optimal_Lags = floor(4*(T/100)^(2/9));

Residual_correlation = Residual*Residual';
Weight_correlation = zeros(size(Residual_correlation));

for i = 1:size(Weight_correlation,1)
    for j = 1:size(Weight_correlation,2)
        if abs(i-j)>Optimal_Lags
            Weight_correlation(i,j)=0;
        else
            Weight_correlation(i,j)=Residual_correlation(i,j)*(1-abs(j-i)/(Optimal_Lags+1));
        end
    end
end

i = (1:size(Weight_correlation,1))';
j = 1:size(Weight_correlation,2);
c = ones(1,size(Weight_correlation,2));
r = ones(size(Weight_correlation,1),1);
lags = abs(i*c-r*j);
ind1 = lags>Optimal_Lags;
ind2 = lags<=Optimal_Lags;
Weight_corr(ind1)=0;
Weight_corr(ind2)=(T/df)*Residual_correlation(ind2).*(1-lags(ind2)/(Optimal_Lags+1));

Check = Weight_correlation - Weight_corr;

Arrays have incompatible sizes for this operation.

Weight_correlation is (10x10) and Weight_corr is (1x100). I reshaped the outcome of this code but here is the result:

Any idea?

[Alan Stevens]

Try this:

Variables = rand(10,10);
Companies = rand(10,1);

Coeff = (Variables'*Variables)\Variables'*Companies;
Expected = Variables*Coeff;
Residual = Companies-Expected;

T = length(Companies);
k = length(Coeff);
df = T-k;

Optimal_Lags = floor(4*(T/100)^(2/9));

Residual_correlation = Residual*Residual';
Weight_correlation = zeros(size(Residual_correlation));

for i = 1:size(Weight_correlation,1)
    for j = 1:size(Weight_correlation,2)
        if abs(i-j)>Optimal_Lags
            Weight_correlation(i,j)=0;
        else
            Weight_correlation(i,j)=Residual_correlation(i,j)*(1-abs(j-i)/(Optimal_Lags+1));
        end
    end
end

i = (1:size(Weight_correlation,1))';
j = 1:size(Weight_correlation,2);
c = ones(1,size(Weight_correlation,2));
r = ones(size(Weight_correlation,1),1);
lags = abs(i*c-r*j);
ind2 = lags<=Optimal_Lags; 
Weight_corr = ind2.*Residual_correlation.*(1-lags/(Optimal_Lags+1));
Check = Weight_correlation - Weight_corr;
disp(Check)

     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0

(I note that you removed the T/df term which gave the inf values!).