diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index c795d33bc..f46c63d26 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -6451,9 +6451,11 @@ block decomposition of the model (see :opt:`block`).
.. option:: analytic_derivation
- Triggers estimation with analytic gradient. The final hessian
- is also computed analytically. Only works for stationary models
- without missing observations, i.e. for ``kalman_algo<3``.
+ Triggers estimation with analytic gradient at ``order=1``.
+ The final hessian at the mode is also computed analytically.
+ Only works for stationary models without missing observations,
+ i.e. for ``kalman_algo<3``. Optimizers that rely on analytic
+ gradients are ``mode_compute=1,3,4,5,101``.
.. option:: ar = INTEGER
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index fd1c7de98..43b000987 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -185,6 +185,11 @@ if DynareOptions.debug
end
DynareOptions.analytic_derivation=ana_deriv;
+if DynareOptions.mode_compute==13
+ error('Options mode_compute=13 is only compatible with quadratic objective functions')
+end
+
+
% if DynareOptions.mode_compute==5
% if ~strcmp(func2str(objective_function),'dsge_likelihood')
% error('Options mode_compute=5 is not compatible with non linear filters or Dsge-VAR models!')
diff --git a/matlab/lpdfgam.m b/matlab/lpdfgam.m
index 7dbcf2579..4a9230bb0 100644
--- a/matlab/lpdfgam.m
+++ b/matlab/lpdfgam.m
@@ -7,13 +7,14 @@ function [ldens,Dldens,D2ldens] = lpdfgam(x,a,b)
% b [double] m*n matrix or scalar, Second GAMMA distribution parameters (scale),
%
% OUTPUTS
-% ldens [double] m*n matrix of logged GAMMA densities evaluated at x.
-%
+% ldens [double] m*n matrix of logged GAMMA densities evaluated at x.
+% Dldens [double] m*n matrix of first derivatives of logged GAMMA densities.
+% D2ldens [double] m*n matrix of second derivatives of logged matrix of logged GAMMA densities.
%
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2003-2017 Dynare Team
+% Copyright (C) 2003-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -42,6 +43,7 @@ end
if nargout >1
+ Dldens = ldens ;
if length(a)==1
Dldens(idx) = (a-1)./(x(idx)) - ones(length(idx),1)/b ;
else
@@ -50,6 +52,7 @@ if nargout >1
end
if nargout == 3
+ D2ldens = ldens ;
if length(a)==1
D2ldens(idx) = -(a-1)./(x(idx)).^2;
else
diff --git a/matlab/lpdfgbeta.m b/matlab/lpdfgbeta.m
index 6f0dd0fd8..7d8e189a2 100644
--- a/matlab/lpdfgbeta.m
+++ b/matlab/lpdfgbeta.m
@@ -2,19 +2,21 @@ function [ldens,Dldens,D2ldens] = lpdfgbeta(x,a,b,aa,bb)
% Evaluates the logged BETA PDF at x.
%
% INPUTS
-% x [double] m*n matrix of loactions,
+% x [double] m*n matrix of locations,
% a [double] m*n matrix of First BETA distribution parameters,
% b [double] m*n matrix of Second BETA distribution parameters,
% aa [double] m*n matrix of lower bounds for (generalized) distribution,
% bb [double] m*n matrix of upper bounds for (generalized) distribution
%
% OUTPUTS
-% ldens [double] m*n matrix of logged (generalized) BETA densities.
+% ldens [double] m*n matrix of logged (generalized) BETA densities.
+% Dldens [double] m*n matrix of first derivatives of logged (generalized) BETA densities.
+% D2ldens [double] m*n matrix of second derivatives of logged matrix of logged (generalized) BETA densities.
%
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2003-2017 Dynare Team
+% Copyright (C) 2003-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -42,6 +44,7 @@ end
if nargout >1
+ Dldens = ldens ;
if length(a)==1
Dldens(idx) = (a-1)./(x(idx)-aa) - (b-1)./(bb-x(idx)) ;
else
@@ -51,6 +54,7 @@ end
if nargout == 3
+ D2ldens = ldens ;
if length(a)==1
D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 - (b-1)./(bb-x(idx)).^2 ;
else
diff --git a/matlab/lpdfig1.m b/matlab/lpdfig1.m
index 2b14c8fab..8de3dbb7a 100644
--- a/matlab/lpdfig1.m
+++ b/matlab/lpdfig1.m
@@ -13,11 +13,13 @@ function [ldens,Dldens,D2ldens] = lpdfig1(x,s,nu)
%
% OUTPUTS
% ldens [double] m*n matrix of logged INVERSE-GAMMA-1 densities evaluated at x.
+% Dldens [double] m*n matrix of first derivatives of logged INVERSE-GAMMA-1 densities.
+% D2ldens [double] m*n matrix of second derivatives of logged matrix of logged INVERSE-GAMMA-1 densities.
%
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2004-2017 Dynare Team
+% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -44,6 +46,7 @@ else
end
if nargout >1
+ Dldens = ldens ;
if length(s)==1
Dldens(idx) = - (nu+1)./(x(idx)) + s./(x(idx).^3) ;
else
@@ -52,6 +55,7 @@ if nargout >1
end
if nargout == 3
+ D2ldens = ldens ;
if length(s)==1
D2ldens(idx) = (nu+1)./(x(idx).^2) - 3*s(idx)./(x(idx).^4) ;
else
diff --git a/matlab/lpdfig2.m b/matlab/lpdfig2.m
index e1b3ca2a1..162bc3eaf 100644
--- a/matlab/lpdfig2.m
+++ b/matlab/lpdfig2.m
@@ -13,11 +13,13 @@ function [ldens,Dldens,D2ldens] = lpdfig2(x,s,nu)
%
% OUTPUTS
% ldens [double] m*n matrix of logged INVERSE-GAMMA-2 densities evaluated at x.
+% Dldens [double] m*n matrix of first derivatives of logged INVERSE-GAMMA-2 densities.
+% D2ldens [double] m*n matrix of second derivatives of logged matrix of logged INVERSE-GAMMA-2 densities.
%
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2004-2017 Dynare Team
+% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -44,6 +46,7 @@ else
end
if nargout >1
+ Dldens = ldens;
if length(s)==1
Dldens(idx) = - .5*(nu+2)./(x(idx)) + .5*s./x(idx).^2;
else
@@ -52,6 +55,7 @@ if nargout >1
end
if nargout == 3
+ D2ldens = ldens;
if length(s)==1
D2ldens(idx) = .5*(nu+2)./(x(idx)).^2 - s./x(idx).^3;
else
diff --git a/matlab/lpdfnorm.m b/matlab/lpdfnorm.m
index 2c1249996..280771d7b 100644
--- a/matlab/lpdfnorm.m
+++ b/matlab/lpdfnorm.m
@@ -7,13 +7,14 @@ function [ldens,Dldens,D2ldens] = lpdfnorm(x,a,b)
% b [double] m*n matrix or scalar, Second GAUSSIAN distribution parameters (standard deviation).
%
% OUTPUTS
-% ldens [double] m*n matrix of logged GAUSSIAN densities evaluated at x.
-%
+% ldens [double] m*n matrix of logged GAUSSIAN densities evaluated at x.
+% Dldens [double] m*n matrix of first derivatives of logged GAUSSIAN densities.
+% D2ldens [double] m*n matrix of second derivatives of logged matrix of GAUSSIAN densities.
%
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2003-2017 Dynare Team
+% Copyright (C) 2003-2021 Dynare Team
%
% This file is part of Dynare.
%
diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m
index c83bae21b..87eb6f1a4 100644
--- a/matlab/method_of_moments/method_of_moments.m
+++ b/matlab/method_of_moments/method_of_moments.m
@@ -93,7 +93,6 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
% - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
% - [ ] improve check for duplicate moments by using the cellfun and unique functions
% - [ ] dirname option to save output to different directory not yet implemented
-% - [ ] add analytic_jacobian option for mode_compute 4 and 101
% The TeX option crashes MATLAB R2014a run with "-nodisplay" option
% (as is done from the testsuite).
@@ -183,9 +182,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
if options_mom_.order > 3
error('method_of_moments: perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM.\n');
end
- options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false); % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
- options_mom_.mom = set_default_option(options_mom_.mom,'analytic_jacobian',false); % use analytic Jacobian in optimization, only available for GMM and gradient-based optimizers
end
+options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false); % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
+options_mom_.mom = set_default_option(options_mom_.mom,'analytic_jacobian',false); % use analytic Jacobian in optimization, only available for GMM and gradient-based optimizers
% initialize flag to compute derivs in objective function (needed for GMM with either analytic_standard_errors or analytic_jacobian )
options_mom_.mom.compute_derivs = false;
@@ -352,7 +351,7 @@ options_mom_.vector_output= false; % specifies whether the objective f
optimizer_vec=[options_mom_.mode_compute;num2cell(options_mom_.additional_optimizer_steps)]; % at each stage one can possibly use different optimizers sequentially
-analytic_jacobian_optimizers = [1, 3, 13]; %these are currently supported, see to-do list
+analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; %these are currently supported, see to-do list
% -------------------------------------------------------------------------
% Step 1d: Other options that need to be initialized
@@ -914,24 +913,11 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
end
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(optimizer_vec{optim_iter},analytic_jacobian_optimizers) %do this only for gradient-based optimizers
options_mom_.mom.compute_derivs = true;
- objective_function1 = str2func('method_of_moments_objective_function_gradient_helper');
- switch optimizer_vec{optim_iter}
- case 1
- options_mom_.optim_opt = [optim_opt0, ',''GradObj'',''on''']; % make sure GradObj option is on for fmincon
- case 3
- options_mom_.optim_opt = [optim_opt0, ',''GradObj'',''on''']; % make sure GradObj option is on for fmincon
- case 13
- options_mom_.optim_opt = [optim_opt0, ',''Jacobian'',''on''']; % make sure Jacobian option is on for lsqnonlin
- end
- if strcmp(options_mom_.optim_opt(1),',')
- options_mom_.optim_opt(1) = []; %remove the comma if optim_opt was empty
- end
else
options_mom_.mom.compute_derivs = false;
- objective_function1 = objective_function;
end
- [xparam1, fval, exitflag] = dynare_minimize_objective(objective_function1, xparam0, optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],...
+ [xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],...
Bounds, oo_, estim_params_, M_, options_mom_);
if options_mom_.vector_output
fval = fval'*fval;
diff --git a/matlab/method_of_moments/method_of_moments_objective_function.m b/matlab/method_of_moments/method_of_moments_objective_function.m
index 3a8548e76..9aecc87a2 100644
--- a/matlab/method_of_moments/method_of_moments_objective_function.m
+++ b/matlab/method_of_moments/method_of_moments_objective_function.m
@@ -1,5 +1,5 @@
-function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
-% [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+% [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
% -------------------------------------------------------------------------
% This function evaluates the objective function for GMM/SMM estimation
% =========================================================================
@@ -15,14 +15,13 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = meth
% o fval: value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
% o info: vector storing error code and penalty
% o exit_flag: 0 if error, 1 if no error
-% o junk1: empty matrix required for optimizer interface
-% o junk2: empty matrix required for optimizer interface
+% o df: analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only)
+% o junk1: empty matrix required for optimizer interface (Hessian would go here)
% o oo_: structure containing the results with the following updated fields:
% - mom.model_moments [numMom x 1] vector with model moments
% - mom.Q value of the quadratic form of the moment difference
% o M_: Matlab's structure describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
-% o df: analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only)
% -------------------------------------------------------------------------
% This function is called by
% o method_of_moments.m
@@ -59,14 +58,18 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = meth
%------------------------------------------------------------------------------
% 0. Initialization of the returned variables and others...
%------------------------------------------------------------------------------
-if options_mom_.vector_output == 1
- if options_mom_.mom.penalized_estimator
- df = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1));
+if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
+ if options_mom_.vector_output == 1
+ if options_mom_.mom.penalized_estimator
+ df = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1));
+ else
+ df = nan(size(oo_.mom.data_moments,1),length(xparam1));
+ end
else
- df = nan(size(oo_.mom.data_moments,1),length(xparam1));
+ df = nan(1,length(xparam1));
end
else
- df = nan(1,length(xparam1));
+ df=[]; %required to be empty by e.g. newrat
end
junk1 = [];
junk2 = [];
diff --git a/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m b/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m
deleted file mode 100644
index f1c88af49..000000000
--- a/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m
+++ /dev/null
@@ -1,55 +0,0 @@
-function [out1, out2] = method_of_moments_objective_function_gradient_helper(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
-% [out1, out2] = method_of_moments_objective_function_gradient_helper(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
-% -------------------------------------------------------------------------
-% This helper function evaluates the objective function for GMM/SMM estimation and
-% outputs the function value fval at first and the gradient df at second place,
-% needed for gradient-based optimizers if analytic_jacobian option is set
-% =========================================================================
-% INPUTS
-% o xparam1: current value of estimated parameters as returned by set_prior()
-% o Bounds: structure containing parameter bounds
-% o oo_: structure for results
-% o estim_params_: structure describing the estimated_parameters
-% o M_ structure describing the model
-% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
-% -------------------------------------------------------------------------
-% OUTPUTS: dependent on the optimizer calling this function, see below
-% -------------------------------------------------------------------------
-% This function calls
-% o method_of_moments_objective_function.m
-% =========================================================================
-% Copyright (C) 2020-2021 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see .
-% -------------------------------------------------------------------------
-% This function is called by
-% o method_of_moments.m
-% o dynare_minimize_objective.m
-% -------------------------------------------------------------------------
-% This function calls
-% o method_of_moments_objective_function
-% -------------------------------------------------------------------------
-% Author(s):
-% o Willi Mutschler (willi@mutschler.eu)
-% =========================================================================
-[fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
-
-% mode_compute=1|3|13 require the gradient as second output argument
-out1=fval;
-out2=df;
-
-end%main function end
-
diff --git a/matlab/optimization/analytic_gradient_wrapper.m b/matlab/optimization/analytic_gradient_wrapper.m
new file mode 100644
index 000000000..8ef31deaf
--- /dev/null
+++ b/matlab/optimization/analytic_gradient_wrapper.m
@@ -0,0 +1,37 @@
+function [fval, grad, hess, exit_flag]=analytic_gradient_wrapper(x, fcn, varargin)
+%function [fval, grad, hess, exitflag]=analytic_gradient_wrapper(x, fcn, varargin)
+% Encapsulates an objective function to be minimized for use with Matlab
+% optimizers
+%
+% INPUTS
+% - x [double] n*1 vector of instrument values.
+% - fcn [fhandle] objective function.
+% - varagin [cell] additional parameters for fcn.
+%
+% OUTPUTS
+% - fval [double] scalar, value of the objective function at x.
+% - grad gradient of the objective function
+% - hess Hessian of the objective function
+% - exit_flag [integer] scalar, flag returned by
+
+% Copyright (C) 2021 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+[fval, info, exit_flag, grad, hess] = fcn(x, varargin{:});
+if size(grad,2)==1
+ grad=grad'; %should be row vector for Matlab; exception lsqnonlin where Jacobian is required
+end
\ No newline at end of file
diff --git a/matlab/optimization/csminwel1.m b/matlab/optimization/csminwel1.m
index a5b94e4ec..9cc6cddc2 100644
--- a/matlab/optimization/csminwel1.m
+++ b/matlab/optimization/csminwel1.m
@@ -42,7 +42,7 @@ function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,m
% http://sims.princeton.edu/yftp/optimize/mfiles/csminwel.m
%
% Copyright (C) 1993-2007 Christopher Sims
-% Copyright (C) 2006-2017 Dynare Team
+% Copyright (C) 2006-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -79,9 +79,9 @@ if ischar(fcn)
end
if ischar(grad)
grad = str2func(grad);
- grad_fun_provided=1;
+ grad_fun_provided = true;
else
- grad_fun_provided=0;
+ grad_fun_provided = false;
end
%tailstr = ')';
%stailstr = [];
@@ -228,16 +228,16 @@ while ~done
f2=f;f3=f;f1=f;retcode2=retcode1;retcode3=retcode1;
end
%how to pick gh and xh
- if f3 < f - crit && badg3==0 && f3 < f2 && f3 < f1
+ if f3 < f - crit && badg3==0 && f3 < f2 && f3 < f1 %f3 has improved function, gradient is good and it is smaller than the other two
ih=3;
fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
- elseif f2 < f - crit && badg2==0 && f2 < f1
+ elseif f2 < f - crit && badg2==0 && f2 < f1 %f2 has improved function, gradient is good and it is smaller than f2
ih=2;
fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
- elseif f1 < f - crit && badg1==0
+ elseif f1 < f - crit && badg1==0 %f1 has improved function, gradient is good
ih=1;
fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
- else
+ else % stuck or bad gradient
[fh,ih] = min([f1,f2,f3]);
%disp_verbose(sprintf('ih = %d',ih))
%eval(['xh=x' num2str(ih) ';'])
@@ -253,12 +253,10 @@ while ~done
%eval(['retcodeh=retcode' num2str(ih) ';'])
retcodei=[retcode1,retcode2,retcode3];
retcodeh=retcodei(ih);
- if exist('gh')
- nogh=isempty(gh);
- else
- nogh=1;
- end
- if nogh
+
+ nogh=isempty(gh);
+ badgh=1;
+ if nogh %recompute gradient
if NumGrad
[gh, badgh]=get_num_grad(method,fcn,penalty,fh,xh,epsilon,varargin{:});
elseif grad_fun_provided
@@ -268,7 +266,6 @@ while ~done
badgh = ~cost_flag;
end
end
- badgh=1;
end
%end of picking
stuck = (abs(fh-f) < crit);
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 512f3c8dc..740534d56 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -70,35 +70,46 @@ switch minimizer_algorithm
error('Optimization algorithm 1 requires the Optimization Toolbox')
end
% Set default optimization options for fmincon.
- optim_options = optimset('display','iter', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
- if ~isoctave
- optim_options = optimset(optim_options, 'LargeScale','off');
- end
+ optim_options = optimoptions('fmincon','display','iter', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
if isoctave
% Under Octave, use the active-set (i.e. octave_sqp of nonlin_min)
% algorithm. On a simple example (fs2000.mod), the default algorithm
% is not able to even move away from the initial point.
- optim_options = optimset(optim_options, 'Algorithm','active-set');
+ optim_options = optimoptions(optim_options, 'Algorithm','active-set');
+ end
+ if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1)
+ optim_options = optimoptions(optim_options,'GradObj','on','TolX',1e-7); %alter default TolX
end
if ~isempty(options_.optim_opt)
- eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
+ eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
end
if options_.silent_optimizer
- optim_options = optimset(optim_options,'display','off');
+ optim_options = optimoptions(optim_options,'display','off');
end
- if options_.analytic_derivation
- optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7);
- end
- if ~isoctave
- [opt_par_values,fval,exitflag,output,lamdba,grad,hessian_mat] = ...
- fmincon(objective_function,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options,varargin{:});
+ if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1) %use wrapper
+ func = @(x) analytic_gradient_wrapper(x,objective_function,varargin{:});
+ if ~isoctave
+ [opt_par_values,fval,exitflag,output,lamdba,grad,hessian_mat] = ...
+ fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
+ else
+ % Under Octave, use a wrapper, since fmincon() does not have an 11th
+ % arg. Also, only the first 4 output arguments are available.
+ [opt_par_values,fval,exitflag,output] = ...
+ fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
+ end
else
- % Under Octave, use a wrapper, since fmincon() does not have an 11th
- % arg. Also, only the first 4 output arguments are available.
- func = @(x) objective_function(x,varargin{:});
- [opt_par_values,fval,exitflag,output] = ...
- fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
+ if ~isoctave
+ [opt_par_values,fval,exitflag,output,lamdba,grad,hessian_mat] = ...
+ fmincon(objective_function,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options,varargin{:});
+ else
+ % Under Octave, use a wrapper, since fmincon() does not have an 11th
+ % arg. Also, only the first 4 output arguments are available.
+ func = @(x) objective_function(x,varargin{:});
+ [opt_par_values,fval,exitflag,output] = ...
+ fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
+ end
end
+
case 2
%simulating annealing
sa_options = options_.saopt;
@@ -157,22 +168,31 @@ switch minimizer_algorithm
error('Optimization algorithm 3 requires the Optimization Toolbox')
end
% Set default optimization options for fminunc.
- optim_options = optimset('display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
+ optim_options = optimoptions('fminunc','display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
if ~isempty(options_.optim_opt)
- eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
- end
- if options_.analytic_derivation
- optim_options = optimset(optim_options,'GradObj','on');
+ eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
end
if options_.silent_optimizer
- optim_options = optimset(optim_options,'display','off');
+ optim_options = optimoptions(optim_options,'display','off');
end
- if ~isoctave
- [opt_par_values,fval,exitflag] = fminunc(objective_function,start_par_value,optim_options,varargin{:});
+ if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1)
+ optim_options = optimoptions(optim_options,'GradObj','on');
+ if ~isoctave
+ func = @(x) analytic_gradient_wrapper(x,objective_function,varargin{:});
+ [opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
+ else
+ % Under Octave, use a wrapper, since fminunc() does not have a 4th arg
+ func = @(x) analytic_gradient_wrapper(x,objective_function,varargin{:});
+ [opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
+ end
else
- % Under Octave, use a wrapper, since fminunc() does not have a 4th arg
- func = @(x) objective_function(x,varargin{:});
- [opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
+ if ~isoctave
+ [opt_par_values,fval,exitflag] = fminunc(objective_function,start_par_value,optim_options,varargin{:});
+ else
+ % Under Octave, use a wrapper, since fminunc() does not have a 4th arg
+ func = @(x) objective_function(x,varargin{:});
+ [opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
+ end
end
case 4
% Set default options.
@@ -480,6 +500,34 @@ switch minimizer_algorithm
% Minimize the penalized objective (note that the penalty is not updated).
[opt_par_values, fval, exitflag, output] = particleswarm(objfun, length(start_par_value), LB, UB, particleswarmOptions);
opt_par_values = opt_par_values(:);
+ case 13
+ % Matlab's lsqnonlin (Optimization toolbox needed).
+ if isoctave && ~user_has_octave_forge_package('optim')
+ error('Option mode_compute=13 requires the optim package')
+ elseif ~isoctave && ~user_has_matlab_license('optimization_toolbox')
+ error('Option mode_compute=13 requires the Optimization Toolbox')
+ end
+ optim_options = optimoptions('lsqnonlin','display','iter','MaxFunEvals',5000,'MaxIter',5000,'TolFun',1e-6,'TolX',1e-6);
+ if ~isempty(options_.optim_opt)
+ eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
+ end
+ if options_.silent_optimizer
+ optim_options.Display='off';
+ end
+ if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1)
+ optim_options.SpecifyObjectiveGradient=true;
+ func = @(x) analytic_gradient_wrapper(x,objective_function,varargin{:});
+ [opt_par_values,Resnorm,fval,exitflag,OUTPUT,LAMBDA,JACOB] = ...
+ lsqnonlin(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
+ else
+ if ~isoctave
+ [opt_par_values,Resnorm,fval,exitflag,OUTPUT,LAMBDA,JACOB] = lsqnonlin(objective_function,start_par_value,bounds(:,1),bounds(:,2),optim_options,varargin{:});
+ else
+ % Under Octave, use a wrapper, since lsqnonlin() does not have a 6th arg
+ func = @(x)objective_function(x,varargin{:});
+ [opt_par_values,Resnorm,fval,exitflag,OUTPUT,LAMBDA,JACOB] = lsqnonlin(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
+ end
+ end
case 101
solveoptoptions = options_.solveopt;
if ~isempty(options_.optim_opt)
@@ -506,7 +554,12 @@ switch minimizer_algorithm
if options_.silent_optimizer
solveoptoptions.verbosity = 0;
end
- [opt_par_values,fval]=solvopt(start_par_value,objective_function,[],[],[],solveoptoptions,varargin{:});
+ if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1)
+ func = @(x) analytic_gradient_wrapper(x,objective_function,varargin{:});
+ [opt_par_values,fval]=solvopt(start_par_value,func,1,[],[],solveoptoptions);
+ else
+ [opt_par_values,fval]=solvopt(start_par_value,objective_function,[],[],[],solveoptoptions,varargin{:});
+ end
case 102
if isoctave
error('Optimization algorithm 2 is not available under Octave')
@@ -514,36 +567,15 @@ switch minimizer_algorithm
error('Optimization algorithm 2 requires the Global Optimization Toolbox')
end
% Set default optimization options for simulannealbnd.
- optim_options = saoptimset('display','iter','TolFun',1e-8);
+ optim_options = optimoptions('simulannealbnd','display','iter','TolFun',1e-8);
if ~isempty(options_.optim_opt)
- eval(['optim_options = saoptimset(optim_options,' options_.optim_opt ');']);
+ eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
end
if options_.silent_optimizer
- optim_options = optimset(optim_options,'display','off');
+ optim_options = optimoptions(optim_options,'display','off');
end
func = @(x)objective_function(x,varargin{:});
[opt_par_values,fval,exitflag,output] = simulannealbnd(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
- case 13
- % Matlab's lsqnonlin (Optimization toolbox needed).
- if isoctave && ~user_has_octave_forge_package('optim')
- error('Option mode_compute=13 requires the optim package')
- elseif ~isoctave && ~user_has_matlab_license('optimization_toolbox')
- error('Option mode_compute=13 requires the Optimization Toolbox')
- end
- optim_options = optimset('display','iter','MaxFunEvals',5000,'MaxIter',5000,'TolFun',1e-6,'TolX',1e-6);
- if ~isempty(options_.optim_opt)
- eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
- end
- if options_.silent_optimizer
- optim_options = optimset(optim_options,'display','off');
- end
- if ~isoctave
- [opt_par_values,Resnorm,fval,exitflag,OUTPUT,LAMBDA,JACOB] = lsqnonlin(objective_function,start_par_value,bounds(:,1),bounds(:,2),optim_options,varargin{:});
- else
- % Under Octave, use a wrapper, since lsqnonlin() does not have a 6th arg
- func = @(x)objective_function(x,varargin{:});
- [opt_par_values,Resnorm,fval,exitflag,OUTPUT,LAMBDA,JACOB] = lsqnonlin(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
- end
otherwise
if ischar(minimizer_algorithm)
if exist(minimizer_algorithm)
diff --git a/matlab/optimization/solvopt.m b/matlab/optimization/solvopt.m
index 56050b074..fb23371f6 100644
--- a/matlab/optimization/solvopt.m
+++ b/matlab/optimization/solvopt.m
@@ -200,7 +200,7 @@ end
epsnorm=1.e-15;
epsnorm2=1.e-30; % epsilon & epsilon^2
-if constr, h1=-1 % NLP: restricted to minimization
+if constr, h1=-1; % NLP: restricted to minimization
cnteps=optim.TolXConstraint; % Max. admissible residual
else
h1=sign(optim.minimizer_indicator); % Minimize resp. maximize a function
@@ -219,7 +219,7 @@ knorms=0; gnorms=zeros(1,10); % Gradient norms stored
%---}
%Display control ---{
-if optim.verbosity<=0, dispdata=0
+if optim.verbosity<=0, dispdata=0;
if optim.verbosity==-1
dispwarn=0;
else
@@ -320,7 +320,7 @@ elseif abs(f)==Inf
end
xrec=x; frec=f; % record point and function value
% Constrained problem
-if constr, fp=f; kless=0
+if constr, fp=f; kless=0;
if trx
fc=feval(func,x');
else
@@ -541,7 +541,7 @@ while 1
xx=x;
end
if kj==0
- kd=4
+ kd=4;
end
kj=kj+1; w=-.9; h=h*2; % use large coef. of space dilation
if kj>2*kd
@@ -574,7 +574,7 @@ while 1
% RESETTING ----{
if kcheck>1
idx=find(abs(g)>ZeroGrad); numelem=size(idx,2);
- if numelem>0, grbnd=epsnorm*numelem^2
+ if numelem>0, grbnd=epsnorm*numelem^2;
if all(abs(g1(idx))<=abs(g(idx))*grbnd) || nrmz==0
if dispwarn
disp(wrnmes)
@@ -801,7 +801,7 @@ while 1
n_grad_evals=n_grad_evals+1;
end
if size(g,2)==1
- g=g'
+ g=g';
end
ng=norm(g);
if isnan(ng)
@@ -926,7 +926,7 @@ while 1
% DISPLAY THE CURRENT VALUES ----{
if k==ld
disp('Iter.# ..... Function ... Step Value ... Gradient Norm ');
- disp(sprintf('%5i %13.5e %13.5e %13.5e',k,f,dx,ng));
+ fprintf('%5i %13.5e %13.5e %13.5e\n',k,f,dx,ng);
ld=k+dispdata;
end
%----}
@@ -941,7 +941,7 @@ while 1
termflag=0;
end
if knan
- termflag=0
+ termflag=0;
end
if kc>=mxtc
termflag=0;
diff --git a/matlab/priordens.m b/matlab/priordens.m
index b988e7ea1..1d402397a 100644
--- a/matlab/priordens.m
+++ b/matlab/priordens.m
@@ -83,8 +83,8 @@ if nargin > 6 && initialization
end
logged_prior_density = 0.0;
-dlprior = 0.0;
-d2lprior = 0.0;
+dlprior = zeros(1,length(x));
+d2lprior = dlprior;
if tt1
logged_prior_density = logged_prior_density + sum(lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1))) ;
@@ -181,9 +181,9 @@ if tt8
return
end
if nargout==2
- [tmp, dlprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8))
+ [tmp, dlprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8));
elseif nargout==3
- [tmp, dlprior(id8), ds2lprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8))
+ [tmp, dlprior(id8), d2lprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8));
end
end
diff --git a/tests/analytic_derivatives/fs2000_analytic_derivation.mod b/tests/analytic_derivatives/fs2000_analytic_derivation.mod
index 1d514405c..677b0f2e1 100644
--- a/tests/analytic_derivatives/fs2000_analytic_derivation.mod
+++ b/tests/analytic_derivatives/fs2000_analytic_derivation.mod
@@ -61,29 +61,74 @@ var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
-steady;
-
-check;
+stoch_simul(order=1,periods=200, irf=0,nomoments,noprint);
+save('my_data.mat','gp_obs','gy_obs');
estimated_params;
-alp, beta_pdf, 0.356, 0.02;
-bet, beta_pdf, 0.993, 0.002;
-gam, normal_pdf, 0.0085, 0.003;
-mst, normal_pdf, 1.0002, 0.007;
-rho, beta_pdf, 0.129, 0.1;
-psi, beta_pdf, 0.65, 0.05;
-del, beta_pdf, 0.01, 0.005;
-stderr e_a, inv_gamma_pdf, 0.035449, inf;
-stderr e_m, inv_gamma_pdf, 0.008862, inf;
+alp, 0.356;
+rho, 0.129;
+psi, 0.65;
+stderr e_a, 0.035449;
+stderr e_m, 0.008862;
end;
varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
-estimation(order=1,mode_compute=9,analytic_derivation,kalman_algo=1,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8);
-estimation(order=1,mode_compute=5,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8);
-estimation(order=1,mode_compute=4,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8);
-estimation(order=1,mode_compute=4,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8);
-//estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=3,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8);
-//estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=4,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8);
\ No newline at end of file
+estimation(order=1,mode_compute=9,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0,prior_trunc=0);
+if (isoctave && user_has_octave_forge_package('optim', '1.6')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
+ estimation(order=1,mode_compute=1,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0
+ %,optim = ('DerivativeCheck', 'on','FiniteDifferenceType','central')
+ );
+ estimation(order=1,mode_compute=3,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+ estimation(order=1,mode_compute=101,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+end
+estimation(order=1,mode_compute=5,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+estimation(order=1,mode_compute=4,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+estimation(order=1,mode_compute=4,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+options_.debug=1;
+estimation(order=1,mode_compute=0,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,plot_priors=0);
+fval_ML_1=oo_.likelihood_at_initial_parameters;
+estimation(order=1,mode_compute=0,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=my_data,nobs=192,mh_replic=0,plot_priors=0);
+fval_ML_2=oo_.likelihood_at_initial_parameters;
+options_.analytic_derivation=0;
+estimation(order=1,mode_compute=0,mode_file=fs2000_analytic_derivation_mode,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,plot_priors=0);
+fval_ML_3=oo_.likelihood_at_initial_parameters;
+
+if abs(fval_ML_1-fval_ML_2)>1e-5 || abs(fval_ML_1-fval_ML_3)>1e-5
+ error('Likelihood does not match')
+end
+options_.debug=0;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+rho, beta_pdf, 0.129, 0.100;
+psi, beta_pdf, 0.65, 0.05;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+estimation(order=1,mode_compute=9,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0,prior_trunc=0);
+if (isoctave && user_has_octave_forge_package('optim', '1.6')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
+ estimation(order=1,mode_compute=1,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0
+ %,optim = ('DerivativeCheck', 'on','FiniteDifferenceType','central')
+ );
+ estimation(order=1,mode_compute=3,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+ estimation(order=1,mode_compute=101,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+end
+estimation(order=1,mode_compute=5,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+estimation(order=1,mode_compute=4,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+estimation(order=1,mode_compute=4,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=my_data,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8,plot_priors=0);
+options_.debug=1;
+estimation(order=1,mode_compute=0,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,plot_priors=0);
+fval_Bayes_1=oo_.likelihood_at_initial_parameters;
+estimation(order=1,mode_compute=0,mode_file=fs2000_analytic_derivation_mode,analytic_derivation,kalman_algo=2,datafile=my_data,nobs=192,mh_replic=0,plot_priors=0);
+fval_Bayes_2=oo_.likelihood_at_initial_parameters;
+options_.analytic_derivation=0;
+estimation(order=1,mode_compute=0,mode_file=fs2000_analytic_derivation_mode,kalman_algo=1,datafile=my_data,nobs=192,mh_replic=0,plot_priors=0);
+fval_Bayes_3=oo_.likelihood_at_initial_parameters;
+
+if abs(fval_Bayes_1-fval_Bayes_2)>1e-5 || abs(fval_Bayes_1-fval_Bayes_3)>1e-5
+ error('Likelihood does not match')
+end
diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod
index 39a58b674..f58c9de04 100644
--- a/tests/estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod
+++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod
@@ -89,7 +89,7 @@ end;
estimated_params_init(use_calibration);
end;
- @#for optimizer in [1, 3, 13]
+ @#for optimizer in [1, 3, 4, 101, 13]
@#if estimParams == 2 && optimizer == 13
%skip due to buggy behavior in Octave
if ~isoctave
@@ -104,7 +104,7 @@ end;
, nograph
, mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance
%, additional_optimizer_steps = [1 3 13]
-% , optim = ('TolFun', 1e-6
+% , optim = ('DerivativeCheck', 'on','FiniteDifferenceType','central'
% ,'TolX', 1e-6
% ,'MaxIter', 3000
% ,'MaxFunEvals', 1D6