function [PSTAR,DSTAR,POS,GRAD,HESS,BSTAR,VSTAR,CSTAR,range] = ...
    estimate(this,d,range,E,varargin)

% Validate required inputs.
IP = inputParser();
IP.addRequired('b',@(x) isa(x,'bkwmodel'));
IP.addRequired('d',@isstruct);
IP.addRequired('range',@(x) isnumeric(x) && ~isempty(x));
IP.addRequired('E',@(x) isstruct(E) ...
    && ~isempty(E) && ~isempty(fieldnames(E)))
IP.parse(this,d,range,E);

% Parese options and optimset for estimation.
options = passvalopt('bkwmodel.estimate',varargin{:});
if iscell(options.optimset)
    options.optimset = optimset(options.optimset{:});
elseif isempty(options.optimset)
    options.optimset = optimset( ...
        'algorithm','active-set', ...
        'display',options.display, ...
        'maxiter',options.maxiter, ...
        'maxfunevals',options.maxfunevals, ...
        'tolx',options.tolx, ...
        'tolfun',options.tolfun);
end

% Parse options and optimset for the Kalman filter.
options.filter = passvalopt('bkwmodel.filter',options.filter{:});
if iscell(options.filter.optimset)
    options.filter.optimset = optimset(options.filter.optimset{:});
elseif isempty(options.filter.optimset)
    options.filter.optimset = optimset( ...
        'display',options.filter.display, ...
        'maxiter',options.filter.maxiter, ...
        'maxfunevals',options.filter.maxfunevals, ...
        'tolx',options.filter.tolx, ...
        'tolfun',options.filter.tolfun);
end

if ischar(options.exclude)
    options.exclude = regexp(options.exclude,'\w+','match');
end
options.exclude(cellfun(@isempty,options.exclude)) = [];

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

if size(this.Assign,3) > 1
    utils.error('bkwmodel', ...
        'Cannot run `estimate` with multiple parameterisations.');
end

maxlag = size(this.occur,3) - 1;
n = sum(this.nametype == 1);
ne = sum(this.nametype == 3);
isidentity = ~isnan(this.identity);
identitynamepos = this.identity(isidentity);
isidentityname = false(1,n);
isidentityname(identitynamepos) = true;

% Get input data.
if ~isinf(range(1))
    xrange = [range(1)-maxlag,range(end)];
else
    xrange = range;
end
[X,xrange] = db2dp(this,d,xrange);

nper = length(xrange);
first = maxlag + 1;
t = first : nper;
range = xrange(t);

% Get variables to be treated as unobserved.
[u,E] = myunobsstruct(this,E);
Y = db2dp(this,u.tseries,xrange);
Y(n+1:end,:) = NaN;
% Index of unobserved data points.
Ya = ~isnan(Y);
isunobs = any(Ya(:));

% Index of residuals and identities to exclude from the objective function.
isexcluderes = false;
excluderes = false(1,ne);
excludeeqtn = false(1,n);
excludey = false(1,n);
excludeident = false(1,n);
if ~isempty(options.exclude)
    exclude_();
end

% Get information about parameters to be estimated.
p = myestimstruct(this,E,0,'warning');
plist = p.plist;
passignpos = p.assignpos;
pstdcorrpos = p.stdcorrpos;
p0 = p.p0;
pl = p.pl;
pu = p.pu;
prior = p.prior;
np = length(plist);
priorindex = p.priorindex;
isprior = any(priorindex);
isparam = np > 0;

filteronly = isunobs && ~isparam;

% Cannot treat identity variables as unobserved.
temp = any(Ya(1:n,:),2) & isidentityname.';
if any(temp)
    utils.error('bkwmodel', ...
        'Cannot treat this identity variable as unobserved: ''%s''.', ...
        this.name{temp});
end

S = permute(this.stdcorr,[2,3,1]);
nobs = sum(~isidentity)*(nper - maxlag);

GRAD = {zeros(1,np),zeros(1,np)};
HESS = {zeros(np),zeros(np)};
if filteronly
    C = covfun.stdcorr2cov(S,ne);
    if isexcluderes
        C = C(~excluderes,~excluderes);
    end
    C = (C+C.')/2;
    pstar = [];
    % Call `dataobjfunc_` which includes a call to `filter_`.
    [objstar,XSTAR,SSTAR,CSTAR,VSTAR] = dataobjfunc_(pstar,'runfilter');
    LAMBDA = struct('lower',zeros(np,1),'upper',zeros(np,1));
elseif isparam
    % Reset persistent reporting variables.
    objfunc_();
    [pstar,objstar,EXITFLAG,OUTPUT,LAMBDA,GRAD{1},HESS{1}] = ...
        fmincon(@objfunc_,p0,[],[],[],[],pl,pu,[],options.optimset);
    [ans,XSTAR,SSTAR,CSTAR,VSTAR] = dataobjfunc_(pstar,'runfilter');
    boundhit = LAMBDA.lower ~= 0 | LAMBDA.upper ~= 0;
    boundhit = boundhit(:).';
    HESS{2} = estimateobj.diffprior(prior,pstar);
end

if options.relative
    temp = SSTAR(1:ne,:);
    temp(~excluderes,:) = temp(~excluderes,:)*sqrt(VSTAR);
    SSTAR(1:ne,:) = temp;
end

BSTAR = this;
if options.relative || any(~isnan(pstdcorrpos))
    BSTAR.stdcorr(1,:,:) = permute(SSTAR,[3,1,2]);
end
if any(~isnan(passignpos))
    BSTAR.Assign(1,passignpos(~isnan(passignpos)),:) = ...
        pstar(~isnan(passignpos));
end

% Posterior simulator object.
POS = poster();
populateposterobj_();

PSTAR = struct();
DSTAR = struct();
%Hess{3} = struct();
%stdvec = sqrt(diag(inv(Hess{1})));
if isparam
    for i = 1 : np
        name = plist{i};
        PSTAR.(name) = pstar(i);
        %Hess{3}.(this.name{assignpos(i)}) = stdvec(i);
    end
end

% Create output database.
template = tseries();
% Convert endogenous and exogenous variables to database tseries.
template.start = xrange(1);
for i = find(this.nametype <= 2)
    name = this.name{i};
    DSTAR.(name) = template;
    DSTAR.(name).data = XSTAR(i,:).';
    DSTAR.(name).Comment = this.namelabel(i);
end
% Convert residuals to database tseries.
template.start = xrange(first);
for i = find(this.nametype == 3)
    name = this.name{i};
    DSTAR.(name) = template;
    DSTAR.(name).data = XSTAR(i,first:nper).';
    DSTAR.(name).Comment = this.namelabel(i);
end

res = XSTAR(this.nametype == 3,first:nper);
CSTAR = res*res.' / (nper - maxlag);

% Nested functions follow.
    
    %***********************************************************************
    % Nested function.
    function [x,res,RESNORM] = filter_(x,C)
        x(Ya) = Y(Ya);
        logx = x;
        logx(this.log,:) = log(x(this.log,:));
        % `C` has been adjusted for the excluded residuals.
        F = chol(C).';
        [logp1,RESNORM,RESIDUAL,EXITFLAG] = lsqnonlin( ...
            @filterobjfunc_,logx(Ya),[],[],options.filter.optimset,logx,F);
        % `res` is adjusted for the excluded residuals.
        [objres,x,res] = filterobjfunc_(logp1,logx,F);
    end
    % End of nested function filter_().
    
    %***********************************************************************
    % Nested function.
    function [objres,x,res] = filterobjfunc_(logp,logx,F)
        logx(Ya) = logp;
        x = logx;
        x(this.log,:) = exp(x(this.log,:));
        % Evaluate identities.
        if any(isidentity)
            % First, set all identity variables except those excluded by the
            % user to NaN to catch recursive identities that are not ordered
            % properly.
            x(this.identity(isidentity & ~excludeident),first:nper) = NaN;
            % Second, evaluate the identities one after another in
            % user-specified order.
            for ieq = find(isidentity & ~excludeident)
                % If there is a lag of the identity variable on the RHS, we
                % need to run the identity period by period. Otherwise, we do
                % it en bloc for all periods.
                islag = any(this.occur(ieq,this.identity(ieq),2:end),3);
                if ~islag
                    aux = [nan(size(x,1),maxlag),x];
                    value = ...
                        this.eqtnEvalIdent{ieq}(aux,maxlag+(1:nper));
                    nanindex = isnan(value);
                    value(1,nanindex) = X(this.identity(ieq),nanindex);
                    x(this.identity(ieq),:) = value;
                else
                    xx = [nan(size(x,1),maxlag),x];
                    for tt = 1 : nper
                        value = this.eqtnEvalIdent{ieq}(xx,maxlag+tt);
                        if ~isnan(value) || tt >= first
                            xx(this.identity(ieq),maxlag+tt) = value;
                        end
                    end
                    x = xx(:,maxlag+1:end);
                end
            end
        end
        % Evaluate residuals.
        res = this.eqtnEvalAllRes(x,t);
        if nargout > 1
            % Assign all residuals before we adjust `res` for the excluded
            % residuals.
            x(this.nametype == 3,first:nper) = res;
        end
        if isexcluderes
            res = res(~excluderes,:);
        end
        % Calculate independent standardised residuals. `F` has been adjusted
        % for the excluded residuals.
        res(isnan(res)) = options.nanresid;
        if any(~isfinite(res(:)))
            chkfiniteres_(res);
        end
        objres = F\res;
    end
    % End of nested function filterobjfunc_().
    
    %***********************************************************************
    % Nested function dataobjfunc_().
    function [L,x,s,C,v] = dataobjfunc_(p,varargin)
        x = X;
        s = S;
        p = p(:);
        index = ~isnan(passignpos);
        if any(index)
            x(passignpos(index),:) = p(index,ones([1,nper]));
        end
        index = ~isnan(pstdcorrpos);
        if any(index)
            s(pstdcorrpos(index)) = p(index);
        end
        C = covfun.stdcorr2cov(s,ne);
        if isexcluderes
            C = C(~excluderes,~excluderes);
        end
        C = (C+C.')/2;
        
        if isunobs && strcmp(varargin{1},'runfilter')
            [x,res,L] = filter_(x,C);
        else
            res = this.eqtnEvalAllRes(x,t);
            if nargout > 1
                % Assign all residuals before we adjust them for the excluded
                % residuals.
                x(this.nametype == 3,first:nper) = res;
            end
            res = res(~excluderes,:);
            res(isnan(res)) = options.nanresid;
            if any(~isfinite(res(:)))
                chkfiniteres_(res);
            end
        end
        
        % `C` and `res` have now the excluded residuals removed.
        W = sum(diag(C\(res*res.')));
        
        % Calculate the common variance scale.
        if options.relative
            v = W / nobs;
            W = nobs;
        else
            v = 1;
        end
        
        if filteronly
            % The objective function value is RESNORM from `filter_` if this
            % is only a filter run.
            return
        end
        
        logdetC = (nper-maxlag)*log(det(C)) + nobs*log(v);
        L = (logdetC + W)/2;
        
        if ~options.logdeta
            return
        end
        
        for j = maxlag+1 : nper
            A = system(this,x,j,true);
            A = A(1:n,1:n);
            if any(isidentity)
                % Substitute for identity variables:
                % Convert [A1,B1;A2,B2]*[X,I] = ... [E;0];
                % to (A1 - B1*(B2\A2))*X = E;
                A1 = A(~isidentity & ~excludeeqtn, ...
                    ~isidentityname & ~excludey);
                A2 = A(isidentity,~isidentityname & ~excludey);
                B1 = A(~isidentity & ~excludeeqtn,isidentityname);
                B2 = A(isidentity,isidentityname);
                A = A1 - B1*(B2\A2);
            else
                if isexcluderes
                    A = A(~excludeeqtn,~excludey);
                end
            end
            if any(~isfinite(A(:)))
                L = Inf;
                return
            end
            logdetA = log(det(A));
            if imag(logdetA) ~= 0
                L = Inf;
                return
            end
            L = L - logdetA;
        end
        
        if ~isfinite(L) || imag(L) ~= 0
            L = Inf;
        end
        
    end
    % End of ested function dataobjfunc_().
    
    %***********************************************************************
    % Nested function priorfunc_().
    function L = priorfunc_(p)
        L = 0;
        for j = find(priorindex)
            L = L - prior{j}(p(j));
            if isinf(L)
                break
            end
        end
    end
    % End of nested function priorfunc_().
    
    %***********************************************************************
    % Nested function objfunc_().
    function L = objfunc_(p)
        persistent COUNT SUB P;
        if nargin == 0
            COUNT = [];
            return
        end
        if ~isinf(options.report)
            if isempty(COUNT)
                COUNT = 0;
                SUB = ceil(sqrt(length(p)));
                figure();
            end
            COUNT = COUNT + 1;
            P = [P,p(:)];
            if mod(COUNT,options.report) == 0
                for j = 1 : length(p)
                    ax = subplot(SUB,SUB,j);
                    plot(1:COUNT,P(j,:));
                    set(ax,'xlim',[1,COUNT]);
                    title(plist{j},'interpreter','none');
                end
            end
            drawnow();
        end
        p = p(:);
        L = 0;
        if isprior
            L = L + priorfunc_(p);
        end
        if ~isinf(L)
            L = L + dataobjfunc_(p,'runfilter');
        end
        if ~isfinite(L)
            L = 1e10;
        end
    end
    % End of nested function objfunc_().
    
    %***********************************************************************
    % Nested function accessfunc_().
    function populateposterobj_()
        POS.paramList = plist;
        POS.minusLogLikFunc = @dataobjfunc_;
        POS.minusLogLikFuncArgs = {};
        POS.logPriorFunc = prior;
        POS.initParam = pstar;
        POS.initProposalCov = inv(HESS{1});
        POS.initLogPost = objstar;
        POS.lowerBounds = pl;
        POS.upperBounds = pu;
    end
    % End of nested function accessfunc_().
    
    %***********************************************************************
    % Nested function.
    function chkfiniteres_(res)
        index = any(~isfinite(res),2);
        if any(index)
            index = index(:).';
            name = this.name(this.nametype == 3);
            message = {};
            for i = find(index)
                message{end+1} = name{i}; %#ok<AGROW>
                dates = xrange(first) + find(~isfinite(res(i,:))) - 1;
                dates = dat2str(dates);
                message{end+1} = sprintf(' %s',dates{:}); %#ok<AGROW>
            end
            utils.error('bkwmodel', ...
                'Residual ''%s'' evaluates to NaN or Inf in these periods:%s.', ...
                message{:});
        end
    end
    % End of nested function chkfiniteres_().
    
    %***********************************************************************
    % Nested function.
    function exclude_()
        found = true(1,length(options.exclude));
        elist = this.name(this.nametype == 3);
        for i = 1 : length(options.exclude)
            % Try to match residual names.
            index = strcmp(options.exclude{i},elist);
            found = any(index);
            excluderes(index) = true;
            % Try to match identity variable names.
            if ~found
                for ieq = find(isidentity)
                    if strcmp(options.exclude{i},this.name{this.identity(ieq)})
                        found = true;
                        excludeident(ieq) = true;
                        break
                    end
                end
            end
        end
        if any(~found)
            utils.error('bkwmodel', ...
                ['This name in the ''exclude'' option is not ', ...
                'a residual or identity variable name: ''%s''.'], ...
                options.exclude{~found});
        end
        isexcluderes = any(excluderes);
        if isexcluderes
            % Index of equations to be excluded.
            excludeeqtn(~isidentity) = excluderes;
            % Try to find a unique index of the endogenous variables to be
            % excluded.
            occur = this.occur(:,this.nametype == 1,1);
            occur = occur(~excludeeqtn,:);
            excludey = all(~occur,1);
            if sum(excludeeqtn) ~= sum(excludey)
                utils.error('bkwmodel', ...
                    ['Equation exclusion does not produce ', ...
                    'a unique square system matrix.'])
            end
        end
    end
    % End of nested function exclude_().
    
end
