function [X,Y,x,y] = fevd(this,time)
% fevd  Forecast error variance decomposition for SVAR variables.
%
% Syntax
% =======
%
%     [X,Y,x,y] = fevd(V,NPER)
%     [X,Y,x,y] = fevd(V,RANGE)
%
% Input arguments
% ================
%
% * `V` [ VAR ] - Structural VAR model.
%
% * `NPER` [ numeric ] - Number of periods.
%
% * `RANGE` [ numeric ] - Date range.
%
% Output arguments
% =================
%
% * `X` [ numeric ] - Forecast error variance decomposition into absolute
% contributions of residuals; absolute contributions sum up to the total
% variance.
%
% * `Y` [ numeric ] - Forecast error variance decomposition into relative
% contributions of residuals; relative contributions sum up to 1.
%
% * `x` [ tseries ] - `X` converted to a tseries object.
%
% * `y` [ tseries ] - `Y` converted to a tseries object.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.
    
    % Tell whether time is nper or range
    if length(time) == 1 && round(time) == time && time > 0
        range = 1 : time;
    else
        range = time(1) : time(end);
    end
    nper = length(range);
    
%**************************************************************************
    
    [ny,p,nalt] = sizeof(this);
    
    Phi = timedom.var2vma(this.A,this.B,nper);
    X = cumsum(Phi.^2,3);
    Y = nan(size(X));
    varvec = this.std .^ 2;
    for ialt = 1 : nalt
        for t = 1 : nper
            if varvec(ialt) ~= 1
                X(:,:,t,ialt) = X(:,:,t,ialt) .* varvec(ialt);
            end
            Xsum = sum(X(:,:,t,ialt),2);
            Xsum = Xsum(:,ones(1,ny));
            Y(:,:,t,ialt) = X(:,:,t,ialt) ./ Xsum;
        end
    end
    
    if nargout > 2
        x = tseries(range,permute(X,[3,1,2,4]));
    end
    
    if nargout > 3
        y = tseries(range,permute(Y,[3,1,2,4]));
    end
    
end
