function OUTP = resample(THIS,INP,RANGE,NDRAW,varargin)
% resample  Resample from the model implied distribution.
%
% Syntax
% =======
%
%     OUTP = resample(M,INP,RANGE,NDRAW,...)
%     OUTP = resample(M,INP,RANGE,NDRAW,J,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Solved model object.
%
% * `INP` [ struct | empty ] - Input data (if needed) for the
% distributions of initial condition and/or empirical shocks; if not
% bootstrap is invovled
%
% * `RANGE` [ numeric ] - Resampling date range.
%
% * `NDRAW` [ numeric ] - Number of draws.
%
% * `J` [ struct | empty ] - Database user-supplied (time-varying) tunes
% on std devs, corr coeffs, and/or means of shocks.
%
% Output arguments
% =================
%
% * `OUTP` [ struct | cell ] - Output database with resampled data.
%
% Options
% ========
%
% * `'deviation='` [ `true` | *`false`* ] - Treat input and output data as
% deviations from balanced-growth path.
%
% * `'dtrends='` [ *`'auto'`* | `true` | `false` ] - Add deterministic trends to
% measurement variables.
%
% * `'method='` [ `'bootstrap'` | *`'montecarlo'`* ] - Method of
% randomising shocks and initial condition.
%
% * `'output='` [ *`'auto'`* | `'dbase'` | `'dpack'` ] - Format of output
% data.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the
% command window.
%
% * `'randomise='` [ *`true`* | `false` | numeric ] - Randomise initial
% condition.
%
% * `'stateVector='` [ *`'alpha'`* | `'x'` ] - When resampling initial
% condition, use the transformed state vector, `'alpha'`, or the vector of
% original variables, `'x'`; this option is meant to guarantee
% replicability of results.
%
% * `'svdOnly='` [ `true` | *`false`* ] - Do not attempt Cholesky and only
% use SVD to factorize the covariance matrix when resampling initial
% condition; only applies when `'randomise=' true`.
%
% * `'wild='` [ `true` | *`false`* ] - Use wild bootstrap instead of Efron
% bootstrap; only applies when `'method=' 'bootstrap'`.
%
% Description
% ============
%
% When you use wild bootstrap for resampling the initial condition, note
% that the results are based on an assumption that the mean of the initial
% condition is the asymptotic mean implied by the model (i.e. the steady
% state).
%
% Example
% ========
%

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

% Called `resample(m,source,range,ndraw,tune,...)'
tune = [];
if ~isempty(varargin) && isstruct(varargin{1})
    % Tunes std devs of shocks.
    % filter(this,data,range,tunes,options...)
    % TODO: tunes on means of shocks.
    tune = varargin{1};
    varargin(1) = [];
end

% Parse required input arguments.
P = inputParser();
P.addRequired('m',@ismodel);
P.addRequired('source',@(x) isnumeric(x) || isstruct(x) || istseries(x));
P.addRequired('range',@(x) isnumeric(x));
P.addRequired('ndraw',@(x) isnumericscalar(x));
P.addRequired('tune',@(x) isempty(x) || isstruct(x));
P.parse(THIS,INP,RANGE,NDRAW,tune);

% Parse options.
opt = passvalopt('model.resample',varargin{:});
inputFormat = model.dataformat(INP,true);
outputFormat = 'dbase';

% If `'dtrends='` option is `'auto'` switch on/off the dtrends according to
% `'deviation='`.
if ischar(opt.dtrends)
    opt.dtrends = ~opt.deviation;
end

% `ninit` is the number of pre-sample periods used to resample the initial
% condition if user does not wish to factorise the covariance matrix.
ninit = 0;
if isnumeric(opt.randomise)
    if isequal(opt.method,'bootstrap')
        utils.error('model', ...
            'Cannot simulate initial conditions in bootstrap resampling.');
    else
        ninit = round(opt.randomise);
        opt.randomise = false;
    end
end

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

% TODO: Allow to combine montecarloed initial condtion and bootstrapped
% shocks, and vice versa.

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

RANGE = RANGE(1) : RANGE(end);
nper = length(RANGE);
realsmall = getrealsmall();
nalt = size(THIS.Assign,3);

% Cannot RESAMPLE from multiple parameterisations.
if nalt > 1
    utils.error('model', ...
        ['Cannot RESAMPLE from model objects ', ...
        'with multiple parameterisations.']);
end

% Check if solution is available.
if isnan(THIS,'solution')
    utils.warning('model', ...
        '#Solution_not_available',' #1');
    OUTP = struct();
    return
end

ny = sum(THIS.nametype == 1);
nx = size(THIS.solution{1},1);
nb = size(THIS.solution{1},2);
nf = nx - nb;
ne = sum(THIS.nametype == 3);

% Combine user-supplied stdcorr with model stdcorr.
usrstdcorr = mytune2stdcorr(THIS,RANGE,tune,opt);
usrstdcorrix = ~isnan(usrstdcorr);

% Get tunes on the mean of shocks.
if ~isempty(tune)
    usrshock = datarequest('e',THIS,tune,RANGE);
    anyusrshock = any(usrshock(:) ~= 0);
end

[T,R,K,Z,H,D,U,Omg] = mysspace(THIS,1,false);
nunit = sum(abs(abs(THIS.eigval)-1) <= realsmall);
nstable = nb - nunit;
Ta = T(nf+1:end,:);
Ra = R(nf+1:end,:);
Ta2 = Ta(nunit+1:end,nunit+1:end);
Ra2 = Ra(nunit+1:end,:);

if opt.dtrends
    W = mydtrendsrequest(THIS,'range',RANGE,1);
end

% Describe the distribution of initial conditions.
if isequal(opt.method,'bootstrap')
    % (1) Bootstrap.
    if ~opt.wild
        % (1a) Efron boostrap.
        sourceAlpha = datarequest('alpha',THIS,INP,RANGE);
    else
        % (1b) Wild bootstrap.
        sourceAlpha0 = datarequest('init',THIS,INP,RANGE);
        Ea = douncmean();
    end
else
    % (2) Monte Carlo or user-supplied sampler.
    if strcmpi(inputFormat,'struct')
        % (2a) User-supplied distribution.
        [Ea,~,~,Pa] = datarequest('init',THIS,INP,RANGE);
        if strcmpi(opt.statevector,'alpha')
            % (2ai) Resample `alpha` vector.
            Fa = covfun.factorise(Pa,opt.svdonly);
        else
            % (2aii) Resample original `x` vector.
            Ex = U*Ea;
            Px = U*Pa*U.';
            Ui = inv(U);
            Fx = covfun.factorise(Px,opt.svdonly);
        end
    else
        % (2b) Asymptotic distribution.
        Ea = douncmean();
        Fa = zeros(nb);
        Pa = zeros(nb);
        Pa(nunit+1:end,nunit+1:end) = covfun.acovf(Ta2,Ra2, ...
            [],[],[],[],[],Omg,THIS.eigval(nunit+1:end),0);
        if strcmpi(opt.statevector,'alpha')
            % (2bi) Resample `alpha` vector.
            Fa(nunit+1:end,nunit+1:end) = covfun.factorise( ...
                Pa(nunit+1:end,nunit+1:end),opt.svdonly);
        else
            % (2bii) Resample original `x` vector.
            Ex = U*Ea;
            Px = U*Pa*U.';
            Ui = inv(U);
            Fx = covfun.factorise(Px,opt.svdonly);
        end
    end
end

% Describe the distribution of shocks.
if isequal(opt.method,'bootstrap')
    % (1) Bootstrap.
    sourceE = datarequest('e',THIS,INP,RANGE);
else
    % (2) Monte Carlo.
    % TODO: Use `mycombinestdcorr` instead.
    stdcorr = permute(THIS.stdcorr,[2,3,1]);
    stdcorr = stdcorr(:,ones(1,nper));
    % Combine the model object stdevs with the user-supplied stdevs.
    if any(usrstdcorrix(:))
        stdcorr(usrstdcorrix) = usrstdcorr(usrstdcorrix);
    end
    % Add model-object std devs for pre-sample if random initial conditions
    % are obtained by simulation.
    if ninit > 0
        stdcorr = [stdcorr(:,ones(1,ninit)),stdcorr];
    end
    
    % Periods in which corr coeffs are all zero. In these periods, we simply
    % mutliply the standard normal shocks by std devs one by one. In all
    % other periods, we need to compute and factorize the entire cov matrix.
    zerocorr = all(stdcorr(ne+1:end,:) == 0, 1);
    if any(~zerocorr)
        Pe = nan(ne,ne,ninit+nper);
        Fe = nan(ne,ne,ninit+nper);
        Pe(:,:,~zerocorr) = ...
            covfun.stdcorr2cov(stdcorr(:,~zerocorr),ne);
        Fe(:,:,~zerocorr) = covfun.factorise(Pe(:,:,~zerocorr));
    end
    
    % If user supplies sampler, sample all shocks and inital conditions at
    % once. This allows for advanced user-supplied simulation methods, e.g.
    % latin hypercube.
    if isa(opt.method,'function_handle')
        allSampleE = opt.method(ne*(ninit+nper),NDRAW);
        if opt.randomise
            allSampleInit = opt.method(nb,NDRAW);
        end
    end
end

% Pre-allocate output data.
OUTP = { ...
    nan(ny,nper,NDRAW), ...
    nan(nx,nper,NDRAW), ...
    nan(ne,nper,NDRAW), ...
    RANGE, ...
    };

% Pre-llocate initial conditions.
init = nan(nx,1,NDRAW);

if ~opt.deviation
    W = mydtrendsrequest(THIS,'range',RANGE);
end

% Distinguish between transition and measurement residuals.
rindex = any(abs(R(:,1:ne)) > 0,1);
hindex = any(abs(H(:,1:ne)) > 0,1);

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

for idraw = 1 : NDRAW
    e = dodrawshocks();
    if ~isempty(tune) && anyusrshock
        e = e + usrshock;
    end
    a0 = dodrawinitcond();
    % Transition variables.
    a = nan(nx,ninit+nper);
    a(:,1) = T*a0 + R(:,rindex)*e(rindex,1);
    if ~opt.deviation
        a(:,1) = a(:,1) + K;
    end
    for t = 2 : ninit+nper
        a(:,t) = T*a(nf+1:end,t-1) + R(:,rindex)*e(rindex,t);
        if ~opt.deviation
            a(:,t) = a(:,t) + K;
        end
    end
    % Measurement variables.
    y = Z*a(nf+1:end,ninit+1:end) + H(:,hindex)*e(hindex,ninit+1:end);
    if ~opt.deviation
        y = y + D(:,ones(1,nper));
    end
    if opt.dtrends
        y = y + W;
    end
    % Store initial condition.
    if ninit == 0
        init(nf+1:end,1,idraw) = a0;
    else
        init(nf+1:end,1,idraw) = a(nf+1:end,ninit);
    end
    % Store this draw.
    OUTP{1}(:,:,idraw) = y;
    OUTP{2}(:,:,idraw) = a(:,ninit+1:end);
    OUTP{3}(:,:,idraw) = e(:,ninit+1:end);
    % Update the progress bar.
    if opt.progress
        update(progress,idraw/NDRAW);
    end
end

% Add pre-sample init cond.
OUTP{1} = [nan(ny,1,NDRAW),OUTP{1}];
OUTP{2} = [init,OUTP{2}];
OUTP{3} = [nan(ne,1,NDRAW),OUTP{3}];
OUTP{4} = [RANGE(1)-1,RANGE];

% Convert alpha to xb;
OUTP = myalpha2xb(THIS,OUTP);

% Convert datapack to database if requested.
if strcmpi(outputFormat,'dbase')
    OUTP = dp2db(THIS,OUTP);
end

% Nested functions.

%**************************************************************************
% Nested function.
    function Ea = douncmean()
        Ea = zeros(nb,1);
        if ~opt.deviation
            Ka2 = K(nf+nunit+1:end);
            Ea(nunit+1:end) = (eye(nstable) - Ta2) \ Ka2;
        end
    end
% douncmean().

%**************************************************************************
    function e = dodrawshocks()
        % Resample residuals.
        if isequal(opt.method,'bootstrap')
            % In boostrap, `ninit` is always zero.
            if opt.wild
                % Wild bootstrap.
                draw = randn(1,nper);
                % To reproduce input sample: draw = ones(1,nper);
                e = sourceE.*draw(ones(1,ne),:);
            else
                % Standard Efron bootstrap.
                % draw is uniform on [1,nper].
                draw = randi([1,nper],[1,nper]);
                % To reproduce input sample: draw = 0 : nper-1;
                e = sourceE(:,draw);
            end
        else
            if isa(opt.method,'function_handle')
                % Fetch and reshape the presampled shocks.
                thisSampleE = allSampleE(:,idraw);
                thisSampleE = reshape(thisSampleE,[ne,ninit+nper]);
            else
                % Draw shocks from standardised normal.
                thisSampleE = randn(ne,ninit+nper);
            end
            % Scale standardised normal by the std devs.
            e = zeros(ne,ninit+nper);
            e(:,zerocorr) = ...
                stdcorr(1:ne,zerocorr) .* thisSampleE(:,zerocorr);
            if any(~zerocorr)
                % Some corrs are non-zero.
                for i = find(~zerocorr)
                    e(:,i) = Fe(:,:,i)*thisSampleE(:,i);
                end
            end
        end
    end
% dodrawshocks().

%**************************************************************************
    function a0 = dodrawinitcond()
        % Randomise initial condition for stable alpha.
        if isequal(opt.method,'bootstrap')
            % Bootstrap from empirical distribution.
            if opt.randomise
                if opt.wild
                    % Wild-bootstrap initial condition for alpha from given
                    % sample initial condition. This assumes that the mean is
                    % the unconditional distribution.
                    a0 = [ ...
                        sourceAlpha0(1:nunit,1); ...
                        Ea2 + randn()*(sourceAlpha0(nunit+1:end,1) - Ea2); ...
                        ];
                else
                    % Efron-bootstrap init cond for alpha from sample.
                    draw = randi([1,nper],[1,nper]);
                    a0 = sourceAlpha(:,draw);
                end
            else
                % Fix init cond to given pre-sample init cond.
                a0 = sourceAlpha0;
            end
        else
            % Gaussian Monte Carlo from theoretical distribution.
            if opt.randomise
                if isa(opt.method,'function_handle')
                    % Fetch the pre-sampled initial conditions.
                    thisSampleInit = allSampleInit(:,idraw);
                else
                    % Draw from standardised normal.
                    thisSampleInit = randn(nb,1);
                end
                if strcmpi(opt.statevector,'alpha')
                    a0 = Ea + Fa*thisSampleInit;
                else
                    x0 = Ex + Fx*thisSampleInit;
                    a0 = Ui*x0;
                end
            else
                % Fix initial conditions to mean.
                a0 = Ea;
            end
        end
    end
% dodrawinitcond().

end