function [PSTAR,POS,COV,HESS,THIS,V,void1,void2,DELTA,PDELTA] = ...
    estimate(THIS,data,range,D,varargin)
% estimate  Estimate model parameters by optimising selected objective function.
%
% Syntax
% =======
%
%     [P,POS,COV,HESS,M,V,~,~,DELTA,PDELTA] = estimate(M,D,RANGE,E,...)
%
% Input arguments
% ================
%
% * `M` [ struct ] - Model object.
%
% * `D` [ struct | cell ] - Input database or datapack from which the
% measurement variables will be taken.
%
% * `RANGE` [ struct ] - Date range.
%
% * `E` [ struct ] - Database with the list of paremeters that will be
% estimated, and the estimation specifications for each of them.
%
% Output arguments
% =================
%
% * `P` [ struct ] - Database with point estimates of requested
% parameters.
%
% * `POS` [ poster ] - Posterior, [`poster`](poster/Contents), object; this
% object also gives you access to the value of the objective function at
% optimum or at any point in the parameter space, see the
% [`poster/eval`](poster/eval) function.
%
% * `COV` [ numeric ] - Approximate covariance matrix for the estimates of
% parameters with slack bounds based on the asymptotic Fisher information
% matrix (not on the Hessian returned from the optimisation routine).
%
% * `HESS` [ cell ] - `HESS{1}` is the total hessian of the objective
% function; `HESS{2}` is the contributions of the priors to the hessian.
%
% * `M` [ model ] - Model object solved with the estimated parameters
% (including out-of-likelihood parameters and common variance factor).
%
% The remaining five output arguments, `V`, `F`, `PE`, `DELTA`, `PDELTA`,
% are the same as the [`model/loglik`](model/loglik) output arguments of the same names.
%
% A tilda, ~, denotes a void (unused) output argument that is included for
% backward compatibility.
%
% Options
% ========
%
% * `'exclude='` [ cellstr | *empty* ] - List of measurement variables that
% are to be excluded from the likelihood function (treated as
% deterministic).
%
% * `'filter='` [ cell | *empty* ] - Cell array of options that will be
% passed on to the Kalman filter; see help on [`model/filter`](model/filter) for the options
% available.
%
% * `'maxIter='` [ numeric | *`500`* ] - Maximum number of iterations allowed.
%
% * `'maxFunEvals='` [ numeric | *`2000`* ] - Maximum number of objective
% function calls allowed.
%
% * `'noSolution='` [ *`'error'`* | `'penalty'` ] - Specifies what happens if
% solution or steady state fails to solve in an iteration: `'error='` stops
% the execution with an error message, `'penalty='` returns an extremely low
% value of the likelihood.
%
% * `'objective='` [ *`'-loglik'`* | `'prederr'` ] - Objective function
% optimised; it can be minus the log likelihood function or weighted sum of
% prediction errors, either adjusted for prior distributions.
%
% * `'objectiveSample='` [ numeric | *`Inf`* ] - The objective function will
% be computed on this subrange only; `Inf` means the entire range will be used.
%
% * `'optimSet='` [ cell | *empty* ] - Cell array used to create the
% Optimization Toolbox options structure; work only with
% `'solver='` `'default'`.
%
% * `'refresh='` [ *`true`* | `false` ] - Refresh dynamic links in each
% iteration.
%
% * `'solve='` [ *`true`* | `false` ] - Re-compute solution in each iteration.
%
% * `'solver='` [ *'default'* | cell | function_handle ] - Minimisation
% procedure; `'default'` means the Optimization Toolbox `fminunc` or
% `fmincon` functions. You can enter a function handle to your own
% optimisation procedure, or a cell array with a function handle and
% additional input arguments -- see below.
%
% * `'sstate='` [ `true` | *`false`* | cell | function_handle ] - Re-compute
% steady state in each iteration. You can specify a cell array with options
% for the `sstate` function, or a function handle whose behaviour is
% described below.
%
% * `'tolFun='` [ numeric | *`1e-6`* ] - Termination tolerance on the
% objective function.
%
% * `'tolX='` [ numeric | *`1e-6`* ] - Termination tolerance on the estimated
% parameters.
%
% Description
% ============
%
% In the input parameter database, `E`, you can provide the following
% four specifications for each parameter:
%
%     E.parameter_name = {start,lower,upper,logprior}
%
% where `start` is the value from which the numerical
% optimisation will start, `lower` is the lower bound, `upper` is the upper
% bound, and `logprior` is a function handle expected to return the log of
% the prior density. You can use the [`logdist`](logdist/Contents) package
% to create function handles for some of the basic prior distributions.
%
% You can use `NaN` for `start` if you wish to use the value currently
% assigned in the model object. You can use `-Inf` and `Inf` for the
% bounds, or leave the bounds empty or not specify them at all. You can
% leave the prior distribution empty or not specify it at all.
%
% --User-supplied optimisation (minimisation) routine
%
% You can supply a function handle to your own minimisation routine through
% the option `'solver='`. This routine will be used instead of the Optim
% Tbx's `fminunc` or `fmincon` functions. The user-supplied function is
% expected to take at least five input arguments and return three output
% arguments:
%
%     [POPT,OBJOPT,HESS] = yourminfunc(F,P0,PLOW,PHIGH,OPT)
%
% with the following input arguments:
%
% * `F` is a function handle to the function minimised;
% * `P0` is a 1-by-N vector of initial parameter values;
% * `PLOW` is a 1-by-N vector of lower bounds (with `-Inf` indicating
% no lower
% bound);
% * `PHIGH` is a 1-by-N vector of upper bounds (with `Inf` indicating no
% upper bounds);
% * `OPT` is an Optim Tbx style struct with the optimisation settings
% (tolerance, number of iterations, etc); of course you may simply ignore
% this information and leave the input argument unused.
%
% and the following output arguments:
%
% * `POPT` is a 1-by-N vector of estimated parameters;
% * `OBJOPT` is the value of the objective function at optimum;
% * `HESS` is a N-by-N approximate hessian matrix at optimum.
%
% If you need to use extra input arguments in your minimisation function,
% enter a cell array instead of a plain function handle:
%
%     {@yourminfunc,ARG1,ARG2,...}
%
% In that case, the solver will be called the following way:
%
%     [POPT,OBJOPT,HESS] = yourminfunc(F,P0,PLOW,PHIGH,OPT,ARG1,ARG2,...)
%
% --User-supplied steady-state solver
%
% You can supply a function handle to your own steady state solver (i.e. a
% function that finds the steady state for given parameters) through the
% `'sstate='` option.
%
% The function is expected to take one input argument, the model object
% with newly assigned parameters, and return at least two output arguments,
% the model object with a new steady state (or balanced-growth path) and a
% success flag. The flag is `true` if the steady state has been successfully
% computed, and `false` if not:
%
%     [M,SUCCESS] = yoursstatesolver(M)
%
% It is your responsibility to add the growth characteristics if some of
% the model variables drift over time. In other words, you need to take
% care of the imaginary parts of the steady state values in the model
% object returned by the solver.
%
% Example
% ========
%

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

% Validate required input arguments.
IP = inputParser();
IP.addRequired('model',@ismodel);
IP.addRequired('data',@(x) isstruct(x) || iscell(x) || isempty(x));
IP.addRequired('range',@isnumeric);
IP.addRequired('parameters',@isstruct);
IP.parse(THIS,data,range,D);

% Process the `estimate` function's own options; what's left is the
% kalman filter options. This is for bkw compatibility only. The kalman
% filter options are entered through a separate `'filter='` option as a
% cell array.
[opt,varargin] = passvalopt('estimateobj.myestimate',varargin{:});
[op1,varargin] = passvalopt('model.estimate',varargin{:});
opt = dbmerge(opt,op1);

% Initialise persistent calls for non-lin steady-state solver.
opt.sstate = mysstateopt(THIS,opt.sstate);
clear('model/mysstatenonlin');

% If the user does not specify the `'filter='` option, create one using
% any unused option.
if isempty(opt.filter)
    opt.filter = varargin;
end

if isempty(THIS.Refresh)
    opt.refresh = false;
end

% Process Kalman filter options; `mypreloglik` also expands solution
% forward if needed for tunes on the mean of shocks.
[loglikopt,THIS,mllfunc] = ...
    mypreloglik(THIS,range,opt.domain,[],opt.filter{:});
loglikopt.ahead = 1;

% Get the first column of input data.
data = datarequest('y',THIS,data,range,1,loglikopt);

% Void output argument for bkw compatibility.
void1 = [];
void2 = [];

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

% Multiple parameterisations not allowed.
nalt = size(THIS.Assign,3);
if nalt > 1
    utils.error('model', ...
        ['Cannot run ESTIMATE on objects ', ...
        'with multiple parameterisations or multiple data sets.']);
end

% Retrieve names of parameters to be estimated, initial values, lower &
% upper bounds, penalties, and prior distributions.
[s,loglikopt] = ...
    myestimstruct(THIS,D,opt.penalty,loglikopt,'warning',opt.initval);
np = length(s.plist);

% Set up objective functions and handles to anonymous functions. This
% can be done only after data have been converted.
s.minusloglik = mllfunc;

% Populate the objective function struct supplying all that is needed to
% evaluate the objective function; also merge `s` with `opt`.
dopopulatestruct();

% BACKEND MYESTIMATE
% ===================

[THIS,pstar,objstar,HESS,boundhit,s] = myestimate(THIS,s);
    
% INFORMATION MATRIX AND COVARIANCE
% ==================================

COV = nan(np);
if np > 0
    COV = myproposalcov(THIS,HESS,boundhit,range,s,loglikopt);
end

% SET UP POSTERIOR OBJECT
% ========================
%
% Set up posterior object before we assign out-of-liks and scale std errors
% in the model object.

POS = poster();
dopopulateposterobj();

% RERUN LOGLIK FOR OUT-OF-LIK PARAMS
% ===================================

% Re-run the Kalman filter or FD likelihood to get the estimates of V and
% out-of-lik parameters.
V = 1;
if nargout >= 5 || loglikopt.relative
    Y = struct('Pdelta',[]);
    [~,V,delta,Y] = mllfunc(THIS,s.data,Y,loglikopt);
    PDELTA = Y.Pdelta;
    DELTA = struct();
    % Assign out-of-lik parameters.
    if ~isempty(loglikopt.outoflik)
        for i = 1 : length(loglikopt.outoflik)
            namepos = loglikopt.outoflik(i);
            name = THIS.name{i};
            THIS.Assign(1,namepos,:) = delta(1,i,:);
            DELTA.(name) = permute(delta(1,i,:),[1,3,2]);
        end
    end
    % 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.
    %
    % Furthermore, there is no need to re-compute the steady state or solution,
    % or refresh the dynamic links because neither std devs nor out-of-lik
    % params are allowed on the RHSs of dynamic links.
    if V ~= 1
        ne = sum(THIS.nametype == 3);
        THIS.stdcorr(1:ne) = THIS.stdcorr(1:ne)*sqrt(V);
    end
end

% Database with point estimates.
PSTAR = cell2struct(num2cell(pstar),s.plist,2);

% Clean-up.
clear('model/mysstatenonlin');

% Nested functions.

%**************************************************************************
    function dopopulatestruct()
        s.data = data;
        s.loglikopt = loglikopt;
        s.assign = THIS.Assign;
        s.stdcorr = THIS.stdcorr;
        optlist = fieldnames(opt);
        for ii = 1 : length(optlist)
            s.(optlist{ii}) = opt.(optlist{ii});
        end
    end
% dopopulatestruct().

%**************************************************************************
    function dopopulateposterobj()
        % Make sure that draws that fail to solve do not cause an error and
        % hence do not interupt the posterior simulator.
        s.nosolution = 'penalty';
        
        POS.paramList = s.plist;
        POS.minusLogPostFunc = @objfunc;
        POS.minusLogPostFuncArgs = {THIS,s};
        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 = s.pl;
        POS.upperBounds = s.pu;
    end
% dopopulateposterobj().

end