function DATA = simulate(THIS,DATA,RANGE,varargin)
% simulate  Simulate VAR model.
%
% Syntax
% =======
%
%     OUTP = simulate(V,INP,RANGE,...)
%
% Input arguments
% ================
%
% * `V` [ VAR ] - VAR object that will be simulated.
%
% * `INP` [ tseries | struct ] - Input data from which the initial
% condtions and residuals will be taken.
%
% * `RANGE` [ numeric ] - Simulation range.
%
% Output arguments
% =================
%
% * `OUTP` [ tseries ] - Simulated output data.
%
% Options
% ========
%
% * `'contributions='` [ `true` | *`false`* ] - Decompose the simulated
% paths into contributions of individual residuals.
%
% * `'deviation='` [ `true` | *`false`* ] - Treat input and output data as
% deviations from unconditional mean.
%
% * `'output='` [ *`'auto'`* | `'dbase'` | `'tseries'` ] - Format of output
% data.
%
% Description
% ============
%
% Example
% ========
%

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

% Parse input arguments.
p = inputParser();
p.addRequired('w',@isvar);
p.addRequired('input',@(x) isnumeric(x) || istseries(x) || isstruct(x));
p.addRequired('range',@isnumeric);
p.parse(THIS,DATA,RANGE);

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

[ny,p,nalt] = sizeof(THIS);

if istseries(DATA) && ~any(size(DATA,2) == [2*ny,ny])
    VAR.error(16,'input time series');
end

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

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

if isempty(RANGE)
    return
end

if RANGE(1) > RANGE(end)
    backcast = true;
    THIS = backward(THIS);
    RANGE = RANGE(end) : RANGE(1)+p;
else
    backcast = false;
    RANGE = RANGE(1)-p : RANGE(end);
end

[outputformat,RANGE,x,e] = VAR.datarequest(THIS,DATA,RANGE,options);
e(isnan(e)) = 0;

if backcast
    x = x(:,end:-1:1,:,:);
    e = e(:,end:-1:1,:,:);
end

e(:,1:p,:) = NaN;
nper = length(RANGE);
ndata = size(x,3);
nloop = max([nalt,ndata]);

if options.contributions
    if nloop > 1
        % Cannot run contributions for multiple data sets or params.
        VAR.error(17,'SIMULATE');
    else
        % Simulation of contributions.
        nloop = ny + 1;
    end
end

if ndata < nloop
    expand = ones(1,nloop-ndata);
    x = cat(3,x,x(:,:,end*expand));
    e = cat(3,e,e(:,:,end*expand));
end

for iloop = 1 : nloop
    if iloop <= nalt
        [Ai,Bi,Ki] = mysystem(THIS,iloop);
    end
    constant = ~options.deviation;
    if options.contributions
        if iloop <= ny
            index = true(1,ny);
            index(iloop) = false;
            e(index,:,iloop) = 0;
            x(:,1:p,iloop) = 0;
            constant = false;
        else
            e(:,:,iloop) = 0;
        end
    end
    if isempty(Bi)
        Bei = e(:,:,iloop);
    else
        Bei = Bi*e(:,:,iloop);
    end
    xi = x(:,:,iloop);
    for t = p + 1 : nper
        xilags = xi(:,t-(1:p));
        xi(:,t) = Ai*xilags(:) + Bei(:,t);
        if constant
            xi(:,t) = xi(:,t) + Ki;
        end
    end
    x(:,:,iloop) = xi;
end

if backcast
    x = x(:,end:-1:1,:,:);
    e = e(:,end:-1:1,:,:);
end

names = get(THIS,'names');
if ~options.returnresiduals
    e = [];
    if ~isempty(names)
        names = names(1:ny);
    end
end

% Output data.
DATA = VAR.outputdata(THIS,outputformat,RANGE,[x;e],[],names);

end
