%%
%Function to Estimate the forecast error variances according to Pesaran and
%Shin (1998)
% Input:
% -estimates: data structure containing estimates Ahat (squared VAR matrix) and Sighat( residual covariance matrix)
% -real: data structure containing real values A (squared VAR matrix) and Sig( residual covariance matrix)
% -H: Forecast horizon. 
% -target: the type of output matrix. Can take values 'FEVD', 'FEV',
% 'MSFE', 'IRF', 'IRV', 'structuralFEVD'

function D=FEVD(estimates,real,H,target)

Ahat=estimates.Ahat;
Sighat=estimates.Sighat;

Areg=estimates.Areg;
SigReg=estimates.SigReg;

[N,~,R,L]=size(Areg);
[~,~,~,~,F]=size(SigReg);


if isstruct(real)==0
    A=eye(N);
    Sig=eye(N);
    Struct=eye(N);
    
else
    A=real.A;
    Sig=real.Sig;
    Struct = real.Struct;
end
R2 = size(A,3);

%%
% Get MA matrices
for r=1:R
    for h=1:H
        MAhat(:,:,h,r)=Ahat(:,:,r)^(h-1);
    end
end

for r=1:R2
    for h=1:H
        MA(:,:,h,r)=A(:,:,r)^(h-1);
    end
end


%% Forecast Errors
%OLS
for r=1:R
    P = Sighat(:,:,r)*sparse(diag(diag( Sighat(:,:,r) ).^(-0.5) ));
    for h=1:H
        summands(:,:,h) = (MAhat(:,:,h,r)*P).^2;
        summands2(:,:,h) = (MAhat(:,:,h,r)*Sighat(:,:,r)*MAhat(:,:,h,r)');
    end
    forecastVarianceEstimated(:,:,r) = sum(summands,3);
    MSEestimated(:,:,r) = sum(summands2,3);
    VarianceDecomposition(:,:,r) = sparse(diag(diag(MSEestimated(:,:,r))))\forecastVarianceEstimated(:,:,r);
end

%Regularized
clear summands summands2
for r=1:R
    for l=1:L
        for f=1:F
            P = SigReg(:,:,r,l,f)*sparse(diag(diag(SigReg(:,:,r,l,f)).^(-0.5) ));
            for h=1:H
                MAreg = Areg(:,:,r,l)^(h-1);
                summands(:,:,h) = (MAreg*P).^2;
                summands2(:,:,h) = (MAreg*SigReg(:,:,r,l,f)*MAreg');
                for i=1:h
                    summandAux(:,:,i) = Areg(:,:,r,l)^(i-1);
                end
                summands3(:,:,h) = sum(summandAux,3)*SigReg(:,:,r,l,f)/sparse(diag(diag(SigReg(:,:,r,l,f))))*SigReg(:,:,r,l,f)*(sum(summandAux,3))';
                correctedSig = (SigReg(:,:,r,l,f));
                summands4(:,:,h) = (MAreg*correctedSig);
            end
            forecastVarianceEstimatedReg(:,:,r,l,f) = sum(summands,3);
            MSEestimatedReg(:,:,r,l,f) = sum(summands2,3);
            VarianceDecompositionReg(:,:,r,l,f) = sparse(diag(diag(MSEestimatedReg(:,:,r,l,f))))\forecastVarianceEstimatedReg(:,:,r,l,f);
            impulseVariance(:,:,r,l,f) = sum(summands4,3);
            impulseVariance2(:,:,r,l,f) = sum(summands3,3);
        end
    end
end

%% Forecast Errors for real matrices
clear summands summands2
for r=1:R2
    P = Sig(:,:,r)*sparse(diag(diag(Sig(:,:,r)).^(-0.5)));
    for h=1:H
        summands(:,:,h) = (MA(:,:,h,r)*P).^2;
        summands2(:,:,h) = (MA(:,:,h,r)*Sig(:,:,r)*MA(:,:,h,r)');
        summands3(:,:,h) = (MA(:,:,h,r)*Struct(:,:,r)).^2;
    end
    FEVreal(:,:,r) = sum(summands,3);
    MSE(:,:,r) = sum(summands2,3);
    FEVDreal(:,:,r) = sparse(diag(diag(MSE(:,:,r))))\FEVreal(:,:,r);
    FEVDrealStruct(:,:,r) = sparse(diag(diag(MSE(:,:,r))))\sum(summands3,3);
end

%% Switch targets

switch target
    case 'FEV'
        D.Reg=forecastVarianceEstimatedReg;
        D.estimated=forecastVarianceEstimated;
        D.real=FEVDreal;
    case 'MSFE'
        D.Reg=MSEestimatedReg;
        D.estimated=MSEestimated;
        D.real=MSE;
    case 'FEVD'
        D.Reg = VarianceDecompositionReg;
        D.estimated = VarianceDecomposition;
        D.real = FEVDreal;
        D.Struct = FEVDrealStruct;
    case 'IRV'
        D.Reg = impulseVariance;
        D.estimated=[];
        D.real=[];
    case 'IRV2'
        D.Reg = impulseVariance2;
        D.estimated=[];
        D.real=[];
    case 'structuralFEVD'
        D.Reg = VarianceDecompositionReg;
        D.estimated = VarianceDecomposition;
        D.real = FEVDrealStruct;
        D.Struct = FEVDreal;

end