function [d,c,f,u,e] = forecast(this,data,range,j,varargin)
% forecast  Forecast FAVAR factors and observables.
%
% Syntax
% =======
%
%     [D,CC,F,U,E] = forecast(A,D,RANGE,J,...)
%
% Input arguments
% ================
%
% * `A` [ FAVAR ] - FAVAR object.
%
% * `D` [ struct | tseries ] - Input data with initial condition for the
% FAVAR factors.
%
% * `RANGE` [ numeric ] - Forecast range.
%
% * `J` [ struct | tseries] - Conditioning data with hard tunes on the
% FAVAR observables.
%
% Output arguments
% =================
%
% * `D` [ struct ] - Output database or tseries object with the FAVAR
% observables.
%
% * `CC` [ struct | tseries ] - Projection of common components in the
% observables.
%
% * `F` [ tseries ] - Projection of common factors.
%
% * `U` [ tseries ] - Conditional idiosyncratic residuals.
%
% * `E` [ tseries ] - Conditional structural residuals.
%
% Options
% ========
%
% See help on [`FAVAR/filter`](FAVAR/filter) for options available.
%
% Description
% ============
%
% Example
% ========
%

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

% Parse required input arguments.
p = inputParser();
p.addRequired('f',@isfavar);
p.addRequired('init',@(x) istseries(x) || isstruct(x));
p.addRequired('range',@isnumeric);
p.addRequired('j',@(x) isempty(x) || istseries(x) || isstruct(x));
p.parse(this,data,range,j);

% Parse options.
options = passvalopt('FAVAR.forecast',varargin{:});

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

ny = size(this.C,1);
nx = size(this.C,2);
p = size(this.A,2)/nx;
range = range(1) : range(end);

if isstruct(data) ...
      && ~isfield(data,'init') ...
      && isfield(data,'mean') ...
      && istseries(data.mean)
   data = data.mean;
end

if istseries(data)
   % Only mean tseries supplied. No uncertainty in initial condition.
   [~,~,x0] = VAR.datarequest(this,data,range(1)-p:range(1)-1);
   x0 = x0(:,end:-1:1,:,:);
   x0 = x0(:);
   P0 = 0;
else
   % Complete description of initial conditions.
   index = abs(range(1)-1 - data.init{3}) <= 0.01;
   if isempty(index) || ~any(index)
      % Initial condition not available.
      utils.error('FAVAR', ...
         'Initial condition for factors not available from input data.');
   end
   x0 = data.init{1}(:,index,:,:);
   P0 = data.init{2}(:,:,index,:,:);
end
nper = length(range);
ndata = size(x0,3);

if ~isempty(j)
   if isstruct(j) && isfield(j,'mean')
      j = j.mean;
   end
   [outputformat,range,y] = VAR.datarequest(this,j,range,options);
   [this,y] = standardise(this,y);
else
   y = nan([ny,nper,ndata]);
   outputformat = options.output;
   if strcmpi(outputformat,'auto')
      if isempty(this.ynames)
         outputformat = 'tseries';
      else
         outputformat = 'dbase';
      end
   end   
end

Sigma = this.Sigma;
% Reduce or zero off-diagonal elements in the cov matrix of idiosyncratic
% residuals if requested.
options.cross = double(options.cross);
if options.cross < 1
   index = logical(eye(size(Sigma)));
   Sigma(~index) = options.cross*Sigma(~index);
end

if isequal(options.invfunc,'auto')
   if this.cross == 1 && options.cross == 1
      invfunc = @pinv;
   else
      invfunc = @inv;
   end
else
   invfunc = options.invfunc;
end

% Run Kalman smoother to re-estimate the common factors taking the
% coefficient matrices as given. If `allobserved` is true then the VAR
% smoother makes an assumption that the factors are uniquely determined
% whenever all observables are available; this is only true if the
% idiosyncratic covariance matrix is not scaled down.
allobserved = this.cross == 1 && options.cross == 1;
[x,Px,e,u,y,Py,yindex] = timedom.varsmoother( ...
   this.A,this.B,[],this.C,[],1,Sigma,y,[],x0,P0,invfunc,allobserved,...
   options.tolerance,options.persist);

if options.meanonly
   Px = Px(:,:,[],[]);
   Py = [];
end

[y,Py] = FAVAR.destandardise(this.Mean,this.Std,y,Py);

ynames = get(this,'ynames');
d = VAR.outputdata(this,outputformat,range,y,Py,ynames);

if nargout > 1
   % Common components.
   [c,Pc] = FAVAR.cc(this.C,x(1:nx,:,:),Px(1:nx,1:nx,:,:));
   [c,Pc] = FAVAR.destandardise(this.Mean,this.Std,c,Pc);
   c = VAR.outputdata(this,outputformat,range,c,Pc,ynames);   
end

if nargout > 2
   f = VAR.outputdata(this,'tseries',range,x(1:nx,:,:),Px(1:nx,1:nx,:,:));
   if ~options.meanonly
      % Means and MSEs that can be used as initial condition.
      lastobs = find(any(yindex ~= 0,1),1,'last');
      if isempty(lastobs)
         index = 1 : nper;
      else
         index = lastobs : nper;
      end
      f.init = {x(:,index,:),Px(:,:,index,:),range(index)};
   end
end

if nargout > 3
   u = FAVAR.destandardise(0,this.Std,u);
   u = VAR.outputdata(this,outputformat,range,u,NaN,ynames);
end

if nargout > 5
   template = tseries();
   e = replace(template,permute(e,[2,1,3]),range(1));
   if ~options.meanonly
   e = struct( ...
      'mean',e, ...
      'std',replace(template,zeros([0,size(e,1),size(e,3)]),range(1)) ...
      );
   end
end

end
