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

Question

Vic 2024-3-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095951-optimizing-code-for-hac-weight-matrix-generation-how-to-be-faster

编辑： Alan Stevens 2024-3-20

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,

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2024-3-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095951-optimizing-code-for-hac-weight-matrix-generation-how-to-be-faster#answer_1427221

在 MATLAB Online 中打开

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.

2 个评论
显示无隐藏无

Vic 2024-3-20

在 MATLAB Online 中打开

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 2024-3-20

编辑：Alan Stevens 2024-3-20

在 MATLAB Online 中打开

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