function [X,LIST,D] = fmse(THIS,TIME,varargin)
% fmse  Forecast mean square error matrices.
%
% Syntax
% =======
%
%     [M,LIST,D] = fmse(M,NPER,...)
%     [M,LIST,D] = fmse(M,RANGE,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Model object for which the forecast MSE matrices will
% be computed.
%
% * `NPER` [ numeric ] - Number of periods.
%
% * `RANGE` [ numeric ] - Date range.
%
% Output arguments
% =================
%
% * `M` [ numeric ] - Forecast MSE matrices.
%
% * `LIST` [ cellstr ] - List of variables in rows and columns of `M`.
%
% * `D` [ dbase ] - Database with the std deviations of
% individual variables, i.e. the square roots of the diagonal elements of
% `M`.
%
% Options
% ========
%
% * `'output='` [ *'namedmat'* | numeric ] - Output matrix `M` will be
% either a namedmat object or a plain numeric array; if the option
% `'select='` is used, `'output='` is always `'namedmat'`.
%
% * `'select='` [ cellstr | *`Inf`* ] - Return the FMSE matrices for
% selected variables only; `Inf` means all variables. The option does not
% apply to the database output `D`.
%
% Description
% ============
%
% Example
% ========
%

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

opt = passvalopt('model.fmse',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 = length(THIS.solutionid{1});
nx = length(THIS.solutionid{2});
nalt = size(THIS.Assign,3);
X = zeros(ny+nx,ny+nx,nper,nalt);

% Compute FMSE for all available parameterisations.
[flag,index] = isnan(THIS,'solution');
for ialt = find(~index)
    [T,R,K,Z,H,D,U,Omg] = mysspace(THIS,ialt,false);
    X(:,:,:,ialt) = timedom.fmse(T,R,K,Z,H,D,U,Omg,nper);
end

% Some solution(s) not available.
if flag
    utils.warning('model', ...
        '#Solution_not_available', ...
        sprintf(' #%g',find(index)));
end

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

% Database of std deviations.
if nargout > 2
    % Select only contemporaneous variables.
    id = [THIS.solutionid{1:2}];
    D = struct();
    for i = find(imag(id) == 0)
        name = THIS.name{id(i)};
        D.(name) = tseries(range,sqrt(permute(X(i,i,:,:),[3,4,1,2])));
    end
    for j = find(THIS.nametype == 4)
        D.(THIS.name{j}) = permute(THIS.Assign(1,j,:),[1,3,2]);
    end
end

% Convert output matrix to namedmat object.
if isnamedmat
    X = namedmat(X,LIST,LIST);
end

% Select variables.
if isselect
    X = select(X,opt.select);
end

end