function [X,Y,LIST,x,y] = fevd(THIS,TIME,varargin)
% fevd  Forecast error variance decomposition for model variables.
%
% Syntax
% =======
%
%     [X,Y,LIST,A,B] = fevd(M,RANGE,...)
%     [X,Y,LIST,A,B] = fevd(M,NPER,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Model object for which the decomposition will be
% computed.
%
% * `RANGE` [ numeric ] - Decomposition date range with the first date
% beign the first forecast period.
%
% * `NPER` [ numeric ] - Number of periods for which the decomposition will
% be computed.
%
% Output arguments
% =================
%
% * `X` [ namedmat | numeric ] - Array with the absolute contributions of
% individual shocks to total variance of each variables.
%
% * `Y` [ namedmat | numeric ] - Array with the relative contributions of
% individual shocks to total variance of each variables.
%
% * `LIST` [ cellstr ] - List of variables in rows of the `X` an `Y`
% arrays, and shocks in columns of the `X` and `Y` arrays.
%
% * `A` [ struct ] - Database with the absolute contributions converted to
% time series.
%
% * `B` [ struct ] - Database with the relative contributions converted to
% time series.
%
% Options
% ========
%
% * `'output='` [ *'namedmat'* | numeric ] - Output matrices `X` and `Y`
% will be either namedmat objects or plain numeric arrays; if the option
% `'select='` is used, `'output='` is always `'namedmat'`.
%
% * `'select='` [ char | cellstr ] - Return the decomposition matrices, `X`
% and `Y`, for selected variables and/or shocks only; `Inf` means all
% variables. This option does not apply to the database outputs, `A` and
% `B`.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

% Parse options.
opt = passvalopt('model.fevd',varargin{:});

if ischar(opt.select)
    opt.select = regexp(opt.select,'\w+','match');
end

% 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);

isselect = iscellstr(opt.select);
isnamedmat = strcmpi(opt.output,'namedmat') || isselect;

%**************************************************************************

ny = sum(THIS.nametype == 1);
nx = length(THIS.solutionid{2});
ne = sum(THIS.nametype == 3);
nalt = size(THIS.Assign,3);
X = nan([ny+nx,ne,nper,nalt]);
Y = nan([ny+nx,ne,nper,nalt]);

% Calculate FEVD for all solved parameterisations that with no
% cross-correlated shocks.
[ans,nosolution] = isnan(THIS,'solution'); %#ok<NOANS,ASGLU>
nonzerocorr = permute(any(THIS.stdcorr(1,ne+1:end,:),2),[1,3,2]);
for ialt = find(~nosolution & ~nonzerocorr)
    [T,R,K,Z,H,D,Za,Omega] = mysspace(THIS,ialt,false);
    [Xi,Yi] = timedom.fevd(T,R,K,Z,H,D,Za,Omega,nper);
    X(:,:,:,ialt) = Xi;
    Y(:,:,:,ialt) = Yi;
end

% Report solution(s) not available.
if any(nosolution)
    utils.warning('model', ...
        '#Solution_not_available', ...
        sprintf(' #%g',find(nosolution)));
end

% Report parameterisations with non-zero cross-correlations.
if any(nonzerocorr)
    temp = sprintf(' #%g',find(nonzerocorr));
    utils.warning('model', ...
        ['Cannot compute FEVD for parameterisations with ', ...
        'non-zero cross-correlations:',temp,'.']);
end

LIST = {[THIS.solutionvector{1:2}],THIS.solutionvector{3}};

% Convert arrays to tseries databases.
if nargout > 3
    % Select only current dated variables.
    id = [THIS.solutionid{1:2}];
    name = THIS.name(real(id));
    ename = THIS.name(THIS.nametype == 3);
    x = struct();
    y = struct();
    for i = find(imag(id) == 0)
        comm = regexprep(ename,'.*',[LIST{1}{i},' <- $0']);
        x.(name{i}) = tseries(range,permute(X(i,:,:,:),[3,2,4,1]),comm);
        y.(name{i}) = tseries(range,permute(Y(i,:,:,:),[3,2,4,1]),comm);
    end
    % Add parameter database.
    for j = find(THIS.nametype == 4)
        x.(THIS.name{j}) = permute(THIS.Assign(1,j,:),[1,3,2]);
        y.(THIS.name{j}) = x.(THIS.name{j});
    end
end

% Convert output matrices to namedmat objects.
if isnamedmat
    X = namedmat(X,LIST{1},LIST{2});
    Y = namedmat(Y,LIST{1},LIST{2});
end

% Select variables; selection only applies to the matrix outputs, `X`
% and `Y`, and not to the database outputs, `x` and `y`.
if isselect
    [X,index] = select(X,opt.select);
    if nargout > 1
        Y = Y(index{1},index{2},:,:);
    end
end

end