function [PSTAR,POS,COV,HESS,THIS,VSTAR] = ...
        estimate(THIS,d,RANGE,E,varargin)

% Validate required inputs.
IP = inputParser();
IP.addRequired('MODEL',@(x) isa(x,'rhsmodel'));
IP.addRequired('DATA',@isstruct);
IP.addRequired('RANGE',@(x) isnumeric(x) && ~isempty(x));
IP.addRequired('ESTPAR',@(x) isstruct(E) ...
    && ~isempty(E) && ~isempty(fieldnames(E)))
IP.parse(THIS,d,RANGE,E);

% Parse estimate-specific options.
[opt,varargin] = passvalopt('rhsmodel.estimate',varargin{:});

% Create optimset.
o = optimset( ...
        'algorithm','interior-point', ...
        'display',opt.display, ...
        'maxIter',opt.maxiter, ...
        'maxFunEvals',opt.maxfunevals, ...
        'tolX',opt.tolx, ...
        'tolFun',opt.tolfun);

if iscell(opt.optimset)
    o = optimset(o,opt.optimset{:});
end

% Parse kalman-specific options.
if ~isempty(opt.filter)
    kalmanopt = passvalopt('rhsmodel.kalman',opt.filter{:});
else
    kalmanopt = passvalopt('rhsmodel.kalman',varargin{:});
end

% Create optimset for Kalman filter.
okalman = rhsmodel.optimset(kalmanopt);

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

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

maxlag = size(THIS.occur,3) - 1;
ne = sum(THIS.nametype == 3);

% Get input data. The input `RANGE` does not include pre-sample.
if ~isinf(RANGE(1))
    xrange = [RANGE(1)-maxlag,RANGE(end)];
else
    xrange = RANGE;
end
[A,X,~,~,sy] = db2dp(THIS,d,xrange);

% Get information about parameters to be estimated.
p = myparamstruct(THIS,E,0,'warning',opt.initval);
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);

if np == 0
    utils.error('rhsmodel', ...
        'No parameters to estimate.');
end

passignposnan = isnan(passignpos);
pstdcorrposnan = isnan(pstdcorrpos);
passignpos = passignpos(~passignposnan);
pstdcorrpos = pstdcorrpos(~pstdcorrposnan);

HESS = {zeros(np),zeros(np)};
VSTAR = 1;

% Populate the obj function struct.
s = struct();
dopopulatestruct();

%chkobj = objfunc(p0,THIS,s,kalmanopt);
keyboard

[pstar,objstar,EXITFLAG,OUTPUT,lambda,GRAD,HESS{1}] = ...
    fmincon(@objfunc,p0,[],[],[],[],pl,pu,[], ...
    o,THIS,s,kalmanopt);

[chkobj,~,~,VSTAR] = objfunc(pstar,THIS,s,kalmanopt);

if abs(objstar-chkobj) > 1e-10
    utils.warning('internal', ...
        'Large discrepany between true OBJFUNC and the value returned from Optim Tbx.');
    keyboard
end
    
% Parameter bounds hit.
boundhit = lambda.lower ~= 0 | lambda.upper ~= 0;
boundhit = boundhit(:).';

% Contribution of priors to the overall hessian.
HESS{2} = estimateobj.diffprior(prior,pstar);

% Covariance matrix of parameter estimates.
COV = HESS{1};
COV(~boundhit,~boundhit) = inv(COV(~boundhit,~boundhit));
COV(boundhit,boundhit) = NaN;

% Update parameters in the output model.
pstar = pstar(:).';
THIS.Assign(1,passignpos) = pstar(1,~passignposnan);
THIS.stdcorr(1,pstdcorrpos) = pstar(1,~pstdcorrposnan);

% Create the `PSTAR` struct before we adjust stddevs.
PSTAR = struct();
for i = 1 : np
    PSTAR.(plist{i}) = pstar(i);
end

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

% Adjust stdevs after we saved the model object in `logPostFuncArgs`. This
% is to keep the model object consistent with `initParam` -- the posterior
% simulator is initialised with the unadjusted stdevs so that the
% lower/upper bounds and the cov matrix remain valid.
if VSTAR ~= 1
    STDSTAR = sqrt(VSTAR);
    THIS.stdcorr(1:ne) = THIS.stdcorr(1:ne)*STDSTAR;
end

% Nested functions.

%**************************************************************************
    function dopopulatestruct()
        s.isprior = isprior;
        s.priorindex = priorindex;
        s.prior = prior;
        s.passignpos = passignpos;
        s.pstdcorrpos = pstdcorrpos;
        s.passignposnan = passignposnan;
        s.pstdcorrposnan = pstdcorrposnan;
        s.o = okalman;
        s.A = A;
        s.X = X;
        s.sy = sy;
    end
% dopopulatestruct().
    
%**************************************************************************  
    function dopopulateposterobj()
        POS.paramList = plist;
        POS.minusLogPostFunc = @objfunc;
        POS.minusLogPostFuncArgs = {THIS,s,kalmanopt};
        POS.initLogPost = -objstar;
        POS.initParam = pstar;
        try
            C = COV;
            % Set the rows and columns corresponding to constraints hit to
            % zeros.
            if any(boundhit)
                for ii = find(boundhit)
                    C(ii,:) = 0;
                    C(:,ii) = 0;
                    C(ii,ii) = 1/HESS{1}(ii,ii);
                end
            end
            POS.initProposalCov = C;
        catch Error
            utils.warning('model', ...
                ['Posterior simulator object cannot be initialised.', ...
                '\nThe following error occurs:\n\n%s'], ...
                Error.message);
        end
        POS.lowerBounds = pl;
        POS.upperBounds = pu;        
    end
% do_populateposterobj().
    
end