function [X,ynames,d] = fmse(this,time,varargin)
% fmse  Forecast mean square error matrices.
%
% Syntax
% =======
%
%     [M,x] = fmse(this,nper)
%     [M,x] = fmse(this,range)
%
% Input arguments
% ================
%
% * `v` [ VAR ] - VAR 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.
%
% * `x` [ dbase | tseries ] - Database or tseries with the std deviations
% of individual variables, i.e. the square roots of the diagonal elements
% of `M`.
%
% Options
% ========
%
% * `'output='` [ 'dbase' | *'tseries'* ] - Format of output data.
%
% Description
% ============
%
% Example
% ========
%
%}

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

options = passvalopt('VAR.fmse',varargin{:});

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

% Orthonormalise residuals so that we do not have to multiply the VMA
% representation by Omega.
B = covfun.factorise(this.Omega);

% Get VMA representation.
X = timedom.var2vma(this.A,B,nper);

% Compute FMSE matrices.
for ialt = 1 : nalt
   for t = 1 : nper
      X(:,:,t,ialt) = X(:,:,t,ialt)*transpose(X(:,:,t,ialt));
   end
end
X = cumsum(X,3);

% Return std devs for individual series.
template = tseries();
if nargout > 1
   x = nan([nper,ny,nalt]);
   for i = 1 : ny
      x(:,i,:) = sqrt(permute(X(i,i,:,:),[3,1,4,2]));
   end
   ynames = get(this,'yNames');
   if strcmpi(options.output,'dbase')
      d = struct();
      for i = 1 : ny
         tmp = x(:,i,:);
         d.(ynames{i}) = replace(template,tmp(:,:),range(1));
      end
   else
      d = replace(template,x,range(1),ynames(1,:,ones([1,nalt])));
   end
end

end
