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

9 次查看（过去 30 天）
Vic
Vic 2024-3-18
编辑： 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,

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

采纳的回答

Alan Stevens
Alan Stevens 2024-3-18
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
Vic 2024-3-20
Ran in:
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
Alan Stevens 2024-3-20
编辑：Alan Stevens 2024-3-20
Ran in:
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!).

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Nonspherical Models 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by