function [data,idata] = forecast(this,data,range,varargin)
% forecast  Unconditional or conditional forecasts.
%
% Syntax
% =======
%
%     [OUTP,INST] = forecast(V,INP,RANGE,J,...)
%
% Input arguments
% ================
%
% * `V` [ VAR ] - VAR object.
%
% * `INP` [ struct | tseries ] - Input data that will be used to set up
% initial condition for the forecast.
%
% * `RANGE` [ numeric ] - Forecast range.
%
% * `J` [ struct | tseries ] - Conditioning database with structural
% tunes on the mean of residuals, reduced-form tunes on the endogenous
% variables, and tunes on the conditioning instruments.
%
% Output arguments
% =================
%
% * `OUTP` [ struct | tseries ] - Output database or output tseries
% object.
%
% * `INST` [ struct | empty ] - Output database with paths for conditioning
% instruments.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright 2007-2012 Jaromir Benes.

jdata = [];
if ~isempty(varargin) && ~ischar(varargin{1})
    jdata = varargin{1};
    varargin(1) = [];
end

% Parse input arguments.
p = inputParser();
p.addRequired('w',@isvar);
p.addRequired('input',@(x) isnumeric(x) || isstruct(x) || istseries(x));
p.addRequired('range',@isnumeric);
p.addRequired('cond',@(x) isempty(x) || isstruct(x) || istseries(x));
p.parse(this,data,range,jdata);

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

% Bkw compatibility only.
if ~isempty(options.returnmse)
    warning('iris:obsolete', ...
        ['\n*** The ''returnMSE'' option is obsolete and will not be supported in future versions. ', ...
        'Use ''meanOnly'' instead.']);
    options.meanonly = ~options.returnmse;
end

[ny,p,nalt,ni] = sizeof(this);
if istseries(data) && ~any(size(data,2) == [2*ny,ny])
    VAR.error(16,'input time series');
end

if istseries(jdata) && ~any(size(jdata,2) == [2*ny,ny])
    VAR.error(16,'conditions');
end

if isequal(range,Inf)
    VAR.error(27,'<<matlab:idoc(''VAR.forecast'') forecast>>');
end

retinstr = nargout > 1;

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

if isempty(range)
    warning('iris:VAR', ...
        '\n*** Forecast range is empty.');
    if options.meanonly
        data = [];
    else
        data = struct();
        data.mean = [];
        data.std = [];
    end
end

if (range(1) > range(end))
    % Go backward in time.
    backcast = true;
    this = backward(this);
    xrange = range(end) : range(1)+p;
    range = range(end) : range(1);
else
    backcast = false;
    xrange = range(1)-p : range(end);
end

% Include pre-sample.
[outputformat,xrange,x0,e] = VAR.datarequest(this,data,xrange,options);
e(isnan(e)) = 0;

% Do not include pre-sample.
[ans,ans,jx,ans,ji] = VAR.datarequest(this,jdata,range); %#ok<NOANS,ASGLU>

if backcast
    x0 = x0(:,end:-1:1,:,:);
    e = e(:,end:-1:1,:,:);
    jx = jx(:,end:-1:1,:,:);
    ji = ji(:,end:-1:1,:,:);
end

nper = length(xrange);
x0 = x0(:,1:p,:);
e = e(:,p+1:end,:);

ndata = size(x0,3);
ncond = size(jx,3);
ninst = size(ji,3);
nloop = max([nalt,ndata,ncond,ninst]);

% Stack initial conditions.
tempx0 = x0(:,p:-1:1,:);
x0 = reshape(tempx0(:),[ny*p,nloop]);

X = nan(ny,nper,nloop);
E = nan(ny,nper,nloop);
P = zeros(ny,ny,nper,nloop);
if retinstr
    I = nan(ni,nper,nloop);
end

Zi = this.Zi;
if isempty(Zi)
    Zi = zeros(0,1+ny*p);
end

for iloop = 1 : nloop
    if iloop <= nalt
        [Ai,Bi,Ki,Omgi] = mysystem(this,iloop);
        if ~isempty(options.omega)
            Omgi(:,:) = options.omega;
        end
        % Use the allobserved option in `varsmoother` only if the cov matrix
        % is full rank. Otherwise, we'd get singularity warnings.
        allobserved = rank(Omgi) == ny;
        % Remove the constant vector if this is a deviation simulation.
        if options.deviation
            Ki = [];
        end
        % Reduce or zero off-diagonal elements in the cov matrix of residuals
        % if requested. This only matters in VARs, not SVARs.
        if double(options.cross) < 1
            index = logical(eye(size(Omgi)));
            Omgi(~index) = double(options.cross)*Omgi(~index);
        end
    end
    if iloop <= ndata
        x0i = x0(:,iloop);
        ei = e(:,:,iloop);
    end
    if iloop <= ncond
        jxi = jx(:,:,iloop);
    end
    if iloop <= ninst && ~isempty(ji)
        index = any(~isnan(ji(:,:,iloop)),2);
        jii = ji(index,:,iloop);
        Cii = Zi(index,1);
        Zii = Zi(index,2:end);
    else
        jii = [];
        Cii = [];
        Zii = [];
    end
    if ~isempty(jii)
        jxi = [jxi;jii]; %#ok<AGROW>
        Z = [eye(ny,ny*p);Zii];
        D = [zeros(ny,1);Cii];
        allobserved = false;
    else
        Z = eye(ny);
        D = zeros(ny,1);
    end
    % Add pre-sample initial condition.
    X(:,p:-1:1,iloop) = reshape(x0i,[ny,p]);
    % Run Kalman filter and smoother.
    [Xi,Pi,E(:,p+1:end,iloop)] = timedom.varsmoother( ...
        Ai,Bi,Ki,Z,D,Omgi,0,jxi,ei,x0i,0,@inv,allobserved);
    X(:,p+1:end,iloop) = Xi(1:ny,:);
    P(:,:,p+1:end,iloop) = Pi(1:ny,1:ny,:);
    % Evaluate conditioning instruments.
    if retinstr
        I(:,p+1:end,iloop) = bsxfun(@plus,Zi(:,1),Zi(:,2:end)*Xi);
    end
    % Set MSE matrix entries to zero for observed data points.
    index = ~isnan(jxi);
    for t = find(any(index,1))
        P(index(:,t),:,p+t,iloop) = 0;
        P(:,index(:,t),p+t,iloop) = 0;
    end
end

if backcast
    X = X(:,end:-1:1,:);
    E = E(:,end:-1:1,:);
    I = I(:,end:-1:1,:);
    P = P(:,:,end:-1:1,:);
end

names = get(this,'names');
if ~options.returnresiduals
    E = zeros(0,size(E,2),size(E,3));
    if ~isempty(names)
        names = names(1:ny);
    end
end

% Output data for endougenous variables and residuals.
if options.meanonly
    data = VAR.outputdata(this,outputformat,xrange,[X;E],[],names);
else
    data = VAR.outputdata(this,outputformat,xrange,[X;E],P,names);
end

% Output data for conditioning instruments.
if retinstr
    inames = this.inames;
    idata = VAR.outputdata(this,outputformat,xrange,I,[],inames);
else
    idata = [];
end

end