function data = resample(w,data,range,ndraw,varargin)
% resample  Resample from a VAR object.
%
% Syntax
% =======
%
%     d = resample(w,d,range,ndraw,...)
%     d = resample(w,[],range,ndraw,...)
%
% Input arguments
% ================
%
% * `w` [ VAR ] - VAR object to resample from.
%
% * `input` [ struct | tseries ] - Input database or tseries used in
% bootstrap; not needed when `'method=' 'montecarlo'`.
%
% * `range` [ numeric ] - Range for which data will be returned.
%
% Output arguments
% =================
%
% * `d` [ struct | tseries ] - Resampled output database or tseries.
%
% Options
% ========
%
% * `'method='` [ 'bootstrap' | *'montecarlo'* | function_handle ] -
% Bootstrap from estimated residuals, resample from normal distribution, or
% use user-supplied sampler.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% * `'randomise='` [ `true` | *`false`* ] - Randomise or fix pre-sample initial
% condition.
%
% * `'wild='` [ `true` | *`false`* ] - Use wild bootstrap instead of standard
% Efron bootstrap when `'method=' 'bootstrap'`.
%
% Description
% ============
%
% Example
% ========
%

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

% Handle obsolete syntax.
throwWarning = false;
if nargin < 4
    % resample(w,data,ndraw)
    ndraw = range;
    range = Inf;
    throwWarning = true;
elseif ischar(ndraw)
    % resample(w,data,ndraw,...)
    varargin = [{ndraw},varargin];
    ndraw = range;
    range = Inf;
    throwWarning = true;
end

if throwWarning
    warning('iris:VAR', ...
        ['Calling VAR.resample with three input arguments is obsolete, ', ...
        'and will not be supported in future versions of IRIS.\n']);
end

% Parse required input arguments.
P = inputParser();
P.addRequired('w',@isvar);
P.addRequired('input',@(x) isnumeric(x) || isstruct(x) || istseries(x));
P.addRequired('range',@(x) isnumeric(x));
P.addRequired('ndraw',@(x) isnumericscalar(x));
P.parse(w,data,range,ndraw);

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

if ischar(options.method)
    options.method = lower(options.method);
end

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

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

% Works only with single parameterisation.
if nalt > 1
    VAR.error(19,'<<matlab:idoc(''VAR.resample'') resample>>');
end

if isempty(range)
    VAR.error(34);
elseif isequal(range,Inf)
    range = w.range(w.fitted);
end

xrange = range(1)-p : range(end);
% `nper` includes pre-sample.
nper = numel(xrange);
[outputformat,xrange,y,e] = VAR.datarequest(w,data,xrange,options);
ndata = size(y,3);

% Works only with a single data set.
if ndata > 1
    VAR.error(20,'<<matlab:idoc(''VAR.resample'') resample>>');
end

% Pre-allocate resampled data.
Y = nan([ny,nper,ndraw]);

% Fixed initial condition.
if isempty(data)
    % Asymptotic initial condition.
    [ans,temp] = mean(w,1);
    Y(:,1:p,:) = temp(:,:,ones([1,ndraw]));
else
    % Initial condition from pre-sample data.
    Y(:,1:p,:) = y(:,1:p,ones([1,ndraw]));
end

% TODO: randomise initial condition
%{
if options.randomise
else
end
%}

[A,B,K] = mysystem(w);

% Back out reduced-form residuals if needed. The B matrix is then
% discarded, and only the covariance matrix of reduced-form residuals is
% used.
if ~isempty(B)
    e = B*e;
end

if ~isequal(options.method,'bootstrap')
    % Safely factorise (chol/svd) the covariance matrix of reduced-form
    % residuals so that we can draw from multivariate normal.
    F = covfun.factorise(w.Omega);
    if isa(options.method,'function_handle')
        allSampleE = options.method(ny*(nper-p),ndraw);
    end
end

% Create a command-window progress bar.
if options.progress
    progress = progressbar('IRIS VAR.resample progress');
end

nanInit = false([1,ndraw]);
nanResid = false([1,ndraw]);
for idraw = 1 : ndraw
    Ei = zeros([ny,nper]);
    Ei(:,p+1:end) = dodrawresiduals();
    Yi = Y(:,:,idraw);
    if any(any(isnan(Yi(:,1:p))))
        nanInit(idraw) = true;
    end
    if any(isnan(Ei(:)))
        nanResid(idraw) = true;
    end
    for t = p+1 : nper
        Yilags = Yi(:,t-(1:p));
        Yi(:,t) = A*Yilags(:) + Ei(:,t) + K;
    end
    Y(:,:,idraw) = Yi;
    % Update the progress bar.
    if options.progress
        update(progress,idraw/ndraw);
    end
end

% Report NaNs in initial conditions.
if any(nanInit)
    temp = sprintf('#%g ',find(nanInit));
    VAR.warning(14,temp(1:end-1));
end

% Report NaNs in resampled residuals.
if any(nanResid)
    temp = sprintf('#%g ',find(nanInit));
    VAR.warning(15,temp(1:end-1));
end

% Return only endogenous variables, not shocks.
data = VAR.outputdata(w,outputformat,xrange,Y,[],w.ynames);

% Nested functions.

%**************************************************************************
    function E = dodrawresiduals()
        if isequal(options.method,'bootstrap')
            if options.wild
                % Wild bootstrap.
                % Setting draw = ones([1,nper-p]) would reproduce sample.
                draw = randn([1,nper-p]);
                E = e(:,p+1:end).*draw(ones([1,ny]),:);
            else
                % Standard Efron bootstrap.
                % Setting draw = 1 : nper-p would reproduce sample;
                % draw is uniform integer [1,nper-p].
                draw = randi([1,nper-p],[1,nper-p]);
                E = e(:,p+draw);
            end
        else
            if isa(options.method,'function_handle')
                thisSampleE = allSampleE(:,idraw);
                thisSampleE = reshape(thisSampleE,[ny,nper-p]);
            else
                thisSampleE = randn([ny,nper-p]);
            end
            E = F*thisSampleE;
        end
    end
% dodrawresiduals().

end

