function d = simulate(this,d,range,varargin)
    
opt = passvalopt('rhsmodel.simulate',varargin{:});
o = rhsmodel.optimset(opt);

if opt.ndraw > 1 && ~isa(opt.residuals,'function_handle')
    utils.error('bkwmodel', ...
        'Options ''ndraw'' and ''residuals'' are set inconsistently.');
end

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

range = range(1) : range(end);
nper = length(range);
maxlag = size(this.occur,3) - 1;
xrange = range(1)-maxlag : range(end);
nxper = length(xrange);
n = length(this.name);
ne = sum(this.nametype == 3);
nalt = size(this.Assign,3);
first = maxlag + 1;
ny = sum(this.nametype == 1);

% Get and check input data. 3rd dimension is adjusted to match `nalt`
% within the `db2dp` function so that `ndata` >= `nalt`.
[A,X] = db2dp(this,d,xrange);
ndata = size(X,3);

%{
% Index of endogenous data points in the en-bloc part of the simulation.
endoga = false(n,nxper);
endoga(name1index,first:nxper) = true;

% Get endogenised/exogenised anchors.
isplan = isa(opt.plan,'plan');
if isplan
    % Check number of exogenised and endogenised data points.
    nexog = nnzexog(opt.plan);
    nendog = nnzendog(opt.plan);
    if nexog ~= nendog
        utils.warning('model', ...
            ['Number of exogenised data points (%g) does not match ', ...
            'number of endogenised data points (%g).'], ...
            nexog(opt.plan),nendog(opt.plan));
    end
    if nendog == 0
        isplan = false;
    else
        [planexoga,planendoga] = myanchors(this,opt.plan,range);
        planexoga = [false(n,maxlag),planexoga];
        planendoga = [false(n,maxlag),planendoga];
        lastanchor = find(any([planexoga;planendoga],1),1,'last');
        if isequal(opt.enbloc,false)
            opt.enbloc = lastanchor;
        end
        endoga(planexoga) = false;
        endoga(planendoga) = true;
    end
end
endoga(:,first:end) = false;

if (opt.ndraw > 1 && nalt > 1) ...
        || (opt.ndraw > 1 && ndata > 1)
    utils.error('bkwmodel', ...
        ['Option ''ndraw'' can be used only with a single parameterisation ', ...
        'and a single data set.']);
elseif opt.ndraw > 1
    % Expand the number of datasets to match `ndraw`.
    X(:,:,end+1:opt.ndraw) = X(:,:,end*ones([1,opt.ndraw-end]));
    ndata = opt.ndraw;
end

% Draw random residuals if requested by the user.
if isa(opt.residuals,'function_handle')
    drawresiduals_();
end
X(this.nametype == 3,1:maxlag,:) = NaN;
%}

% Check for NaN initial conditions.
do_chkinit();

% Check for NaN, Inf, imag exogenous variables.
% chkexog_();
% Check for exogenised NaNs or Infs.
%if isplan
%     chkplan_();
%end

nloop = ndata;

if opt.progress
    progress = progressbar('IRIS rhsmodel.simulate progress');
end

for iloop = 1 : nloop
    
    x = X(:,:,iloop);
    a = A(:,:,iloop);
    
    for t = first : nxper
        switch opt.solvefor
            case 'endogenous'
                a = prediction(this,a,x,t,o);
            case 'residuals'
                x = residuals(this,a,x,t,o);
        end
    end
    
    X(:,:,iloop) = x;
    A(:,:,iloop) = a;
    
    if opt.progress
        update(progress,iloop/nloop);
    end
    
end

d = dp2db(this,A,X,xrange);

% Nested functions follow.

    % @ *******************************************************************
    function drawresiduals_()
        index = this.nametype == 3;
        C = covfun.stdcorr2cov(permute(this.stdcorr,[2,3,1]),ne);
        if nalt == 1
            F = chol(C).';
            X(index,maxlag+(1:nper),:) = reshape( ...
                F*opt.residuals([ne,nper*ndata]), ...
                [ne,nper,ndata]);
        else
            X(index,maxlag+(1:nper),:) = reshape( ...
                opt.residuals([ne,nper*ndata]), ...
                [ne,nper,ndata]);
            for idata = 1 : ndata
                if idata <= nalt
                    F = chol(C(:,:,idata)).';
                end
                X(index,maxlag+(1:nper),idata) = ...
                    F*X(index,maxlag+(1:nper),idata);
            end
        end
    end
    % @ drawresiduals_().
    
    % @ *******************************************************************
    % Nested function.
    function do_chkinit()
        % Check for NaN initial conditions in endogenous variables.
        nanindex = any(isnan(A(:,first-1,:)),3);
        if any(nanindex)
            index = unique(real(this.systemid(nanindex)));
            utils.error('bkwmodel', ...
                'Initial condition not available for this variable: ''%s''.', ...
                this.name{index});
        end
    end
    % @ chkinit_().
    
    % @ *******************************************************************
    function chkexog_()
        nanindex = any(isnan(X),3);
        infimagindex = any(isinf(X) | imag(X) ~= 0,3);
        nanexog = false(size(this.name));
        infimagexog = false(size(this.name));
        t = maxlag+1 : nper;
        for i = find(this.nametype == 2)
            lags = -(find(any(this.occur(:,i,:),1)) - 1);
            lags = lags(:).';
            pos = [];
            for j = maxlag+ (1 : nper)
                pos = [pos,j+lags];
            end
            pos = unique(pos);
            nanexog(i) = any(nanindex(i,pos));
            infimagexog(i) = any(infimagindex(i,pos));
        end
        if any(nanexog)
            utils.error('bkwmodel', ...
                ['This exogenous variables has NaNs ', ...
                'within the simulation range: ''%s''.'], ...
                this.name{nanexog});
        end
        if any(infimagexog)
            utils.error('bkwmodel', ...
                ['This exogenous variables has results in Inf or imag ', ...
                'within the simulation range: ''%s''.'], ...
                this.name{infimagexog});
        end
    end
    % @ chkexog_().
    
    % @ *******************************************************************
    function chkplan_()
        xnan = any(~isfinite(X),3);
        if ~any(xnan(planexoga))
            return
        end
        missing = false(1,n);
        for j = 1 : ny
            xnanj = xnan(j,first:nper);
            exogaj = exoga(j,first:nper);
            missing(j) = any(xnanj(exogaj));
        end
        utils.error('bkwmodel', ...
            'This variable is exogenised to NaNs or Infs: ''%s''.', ...
            this.name{missing});
    end
    % @ chkplan_().
    
end