function [d,exitflag,DISCREP] = simulate(this,d,range,varargin)

options = passvalopt('bkwmodel.simulate',varargin{:});

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

if iscell(options.optimset)
    options.optimset = optimset(options.optimset{:});
elseif isempty(options.optimset)
    options.optimset = optimset( ...
        'display',options.display, ...
        'maxiter',options.maxiter, ...
        'maxfunevals',options.maxfunevals, ...
        'tolx',options.tolx, ...
        'tolfun',options.tolfun);
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);
name1index = this.nametype == 1;
ne = sum(this.nametype == 3);
nalt = size(this.Assign,3);
first = maxlag + 1;

% Get and check input data. 3rd dimension is adjusted to match `nalt`
% within the `db2dp` function so that `ndata` >= `nalt`.
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(options.plan,'plan');
if isplan
    % Check number of exogenised and endogenised data points.
    nexog = nnzexog(options.plan);
    nendog = nnzendog(options.plan);
    if nexog ~= nendog
        utils.warning('model', ...
            ['Number of exogenised data points (%g) does not match ', ...
            'number of endogenised data points (%g).'], ...
            nexog(options.plan),nendog(options.plan));
    end
    if nendog == 0
        isplan = false;
    else
        [planexoga,planendoga] = myanchors(this,options.plan,range);
        planexoga = [false(n,maxlag),planexoga];
        planendoga = [false(n,maxlag),planendoga];
        lastanchor = find(any([planexoga;planendoga],1),1,'last');
        if isequal(options.enbloc,false)
            options.enbloc = lastanchor;
        end
        endoga(planexoga) = false;
        endoga(planendoga) = true;
    end
end

if islogical(options.enbloc)
    if options.enbloc
        enbloc = nxper;
    else
        enbloc = first-1;
    end
else
    enbloc = options.enbloc;
end
endoga(:,enbloc+1:end) = false;

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

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

if any(this.log)
    X(this.log,:,:) = log(X(this.log,:,:));
end

% Check for NaN initial conditions.
dochkinit();
% Check for NaN, Inf, imag exogenous variables.
dochkexog();
% Check for exogenised NaNs or Infs.
if isplan
    dochkplan();
end

nloop = ndata;

neqtn = length(this.eqtnEval);
exitflag = false([1,nloop]);
discrep = nan([neqtn,nxper,nloop]);

if options.deviation
    X0 = X;
    X0(this.nametype == 3,maxlag+1:end,:) = 0;
end

if options.progress
    progress = progressbar('IRIS bkwmodel.simulate progress');
elseif options.esttime
    eta = esttime('IRIS bkwmodel.simulate is running');
end

for idev = 0 : double(options.deviation)
    
    for iloop = 1 : nloop
        
        if enbloc >= first
            if idev == 0
                x = X(:,:,iloop);
            else
                x = X0(:,:,iloop);
            end
            % Set up initial condition.
            x0 = x;
            for t = first : enbloc
                y0 = x0(this.nametype <= 2,t);
                index = ~isfinite(y0);
                if any(index) && t > 1
                    y1 = x0(this.nametype == 1,t-1);
                    y0(index) = y1(index);
                    index = ~isfinite(y0);
                end
                y0(index) = 0;
                x0(this.nametype <= 2,t) = y0;
                e0 = x0(this.nametype == 3,t);
                index = ~isfinite(e0);
                e0(index) = 0;
                x0(this.nametype == 3,t) = e0;
            end
            x0 = x0(endoga);
            t = first : enbloc;
            if isequal(options.solver,'lsqnonlin')
                [xstar,ans,ans,thisexitflag] = ...
                    lsqnonlin(@doobjfunc,x0,[],[],options.optimset,endoga);
            else
                [xstar,fval,thisexitflag] = ...
                    fsolve(@doobjfunc,x0,options.optimset,endoga);
            end
            exitflag(iloop) = thisexitflag > 0;
            if nargout > 2
                discrep(:,maxlag+1:end,iloop) = doobjfunc(xstar,endoga);
            end
            x(endoga) = xstar;
            if idev == 0
                X(:,:,iloop) = x;
            else
                X0(:,:,iloop) = x;
            end
        end
        
        if enbloc < nxper
            if idev == 0
                x = X(:,:,iloop);
            else
                x = X0(:,:,iloop);
            end
            % Last pre-sample period.
            thisexitflag = nan([1,nxper]);
            for t = enbloc+1 : nxper
                % Set up initial conditions.
                x0 = x(name1index,t-1);
                x0(~isfinite(x0)) = 0;
                if isequal(options.solver,'lsqnonlin')
                    [thisxstar,ans,ans,thisexitflag(t)] = ...
                        lsqnonlin(@doobjfunc,x0,[],[],options.optimset);
                else
                    [thisxstar,fval,thisexitflag(t)] = ...
                        fsolve(@doobjfunc,x0,options.optimset);
                end
                x(name1index,t) = thisxstar;
                if nargout > 2
                    discrep(:,t,iloop) = doobjfunc(thisxstar);
                end
            end
            exitflag(iloop) = all(thisexitflag(enbloc+1:nxper) > 0);
            xstar = x(name1index,enbloc+1:nxper);
            if idev == 0
                X(name1index,enbloc+1:nxper,iloop) = xstar;
            else
                X0(name1index,enbloc+1:nxper,iloop) = xstar;
            end
        end
        
        if options.progress
            update(progress,iloop/nloop);
        elseif options.esttime
            update(eta,iloop/nloop);
        end
        
    end
    
end

if options.deviation
    X = X - X0;
end

if any(this.log)
    X(this.log,:,:) = exp(X(this.log,:,:));
end

if any(~exitflag)
    index = num2cell(find(~exitflag));
    utils.warning('bkwmodel', ...
        'This simulation run did not converge: #%g.',index{:});
end

d = dp2db(this,X,xrange);

if nargout > 2
    DISCREP = ...
        tseries(range,permute(discrep(:,maxlag+1:end,:),[2,1,3]), ...
        this.eqtnlabel(1,:,ones([1,nloop])));
end

% Nested functions.

%***********************************************************************
    function obj = doobjfunc(p,endoga)
        xx = x;
        if nargin == 1
            xx(this.nametype == 1,t) = p;
        else
            xx(endoga) = p;
        end
        if any(this.log)
            xx(this.log,:) = exp(xx(this.log,:));
        end
        obj = this.eqtnEvalAll(xx,t);
    end
% doobjfunc().

%***********************************************************************
    function dodrawresiduals()
        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*options.residuals([ne,nper*ndata]),[ne,nper,ndata]);
        else
            X(index,maxlag+(1:nper),:) = ...
                reshape(options.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
% dodrawresiduals().

%***********************************************************************
    function dochkinit()
        % Check for NaN initial conditions in endogenous variables.
        nanIndex = any(isnan(X(:,1:maxlag,:)),3);
        nanInit = false(size(this.name));
        for i = find(this.nametype == 1)
            lagIndex = any(this.occur(:,i,2:end),1);
            lagIndex = lagIndex(:).';
            lagIndex = lagIndex(end:-1:1);
            nanInit(i) = any(nanIndex(i,:) & lagIndex);
        end
        if any(nanInit)
            utils.error('bkwmodel', ...
                'Initial condition not available for this variable: ''%s''.', ...
                this.name{nanInit});
        end
    end
% dochkinit().

%***********************************************************************
    function dochkexog()
        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
% dochkexog().

%***********************************************************************
    function dochkplan()
        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
% dochkplan().

end