% This file computes variance decompositions

% Need to know IRF, p, h, and q

function [DECOMP,TOT] = vardecomp_eric(IRF,h,p,q);

SIRF = IRF.^2; % This squares every element of the IRF

for i = 1:q^2
    SIRF(i,:) = cumsum(SIRF(i,:));
end

TOT = zeros(h,q);

    
    
for i = 1:q
    for j = 1:h
        TOT(j,i) = sum(SIRF(i:q:q^2,j));
    end
end

% The above calculates the total variance for each variable at each horizon


% Now we've got the cumulative sum

DECOMP = zeros(q^2,h);
       
       
for i = 1:h
    for j = 1:q
        for k = 1:q
            
            DECOMP((j-1)*q + k,i) = SIRF((j-1)*q + k,i)/TOT(i,k);
        end
    end
end

% Now do a check to make sure everything adds to 1:

CHECK = zeros(h,q);

for i = 1:q
    for j = 1:h
        CHECK(j,i) = sum(DECOMP(i:q:q^2,j));
    end
end

% and it checks out!

        
 