Merge branch 'analytic_derivation' of git.dynare.org:JohannesPfeifer/dynare
commit
f9437f89c2
|
@ -6451,9 +6451,11 @@ block decomposition of the model (see :opt:`block`).
|
||||||
|
|
||||||
.. option:: analytic_derivation
|
.. option:: analytic_derivation
|
||||||
|
|
||||||
Triggers estimation with analytic gradient. The final hessian
|
Triggers estimation with analytic gradient at ``order=1``.
|
||||||
is also computed analytically. Only works for stationary models
|
The final hessian at the mode is also computed analytically.
|
||||||
without missing observations, i.e. for ``kalman_algo<3``.
|
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
|
.. option:: ar = INTEGER
|
||||||
|
|
||||||
|
|
|
@ -185,6 +185,11 @@ if DynareOptions.debug
|
||||||
end
|
end
|
||||||
DynareOptions.analytic_derivation=ana_deriv;
|
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 DynareOptions.mode_compute==5
|
||||||
% if ~strcmp(func2str(objective_function),'dsge_likelihood')
|
% if ~strcmp(func2str(objective_function),'dsge_likelihood')
|
||||||
% error('Options mode_compute=5 is not compatible with non linear filters or Dsge-VAR models!')
|
% error('Options mode_compute=5 is not compatible with non linear filters or Dsge-VAR models!')
|
||||||
|
|
|
@ -8,12 +8,13 @@ function [ldens,Dldens,D2ldens] = lpdfgam(x,a,b)
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% 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
|
% SPECIAL REQUIREMENTS
|
||||||
% none
|
% none
|
||||||
|
|
||||||
% Copyright (C) 2003-2017 Dynare Team
|
% Copyright (C) 2003-2021 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -42,6 +43,7 @@ end
|
||||||
|
|
||||||
|
|
||||||
if nargout >1
|
if nargout >1
|
||||||
|
Dldens = ldens ;
|
||||||
if length(a)==1
|
if length(a)==1
|
||||||
Dldens(idx) = (a-1)./(x(idx)) - ones(length(idx),1)/b ;
|
Dldens(idx) = (a-1)./(x(idx)) - ones(length(idx),1)/b ;
|
||||||
else
|
else
|
||||||
|
@ -50,6 +52,7 @@ if nargout >1
|
||||||
end
|
end
|
||||||
|
|
||||||
if nargout == 3
|
if nargout == 3
|
||||||
|
D2ldens = ldens ;
|
||||||
if length(a)==1
|
if length(a)==1
|
||||||
D2ldens(idx) = -(a-1)./(x(idx)).^2;
|
D2ldens(idx) = -(a-1)./(x(idx)).^2;
|
||||||
else
|
else
|
||||||
|
|
|
@ -2,7 +2,7 @@ function [ldens,Dldens,D2ldens] = lpdfgbeta(x,a,b,aa,bb)
|
||||||
% Evaluates the logged BETA PDF at x.
|
% Evaluates the logged BETA PDF at x.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% 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,
|
% a [double] m*n matrix of First BETA distribution parameters,
|
||||||
% b [double] m*n matrix of Second BETA distribution parameters,
|
% b [double] m*n matrix of Second BETA distribution parameters,
|
||||||
% aa [double] m*n matrix of lower bounds for (generalized) distribution,
|
% aa [double] m*n matrix of lower bounds for (generalized) distribution,
|
||||||
|
@ -10,11 +10,13 @@ function [ldens,Dldens,D2ldens] = lpdfgbeta(x,a,b,aa,bb)
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% 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
|
% SPECIAL REQUIREMENTS
|
||||||
% none
|
% none
|
||||||
|
|
||||||
% Copyright (C) 2003-2017 Dynare Team
|
% Copyright (C) 2003-2021 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -42,6 +44,7 @@ end
|
||||||
|
|
||||||
|
|
||||||
if nargout >1
|
if nargout >1
|
||||||
|
Dldens = ldens ;
|
||||||
if length(a)==1
|
if length(a)==1
|
||||||
Dldens(idx) = (a-1)./(x(idx)-aa) - (b-1)./(bb-x(idx)) ;
|
Dldens(idx) = (a-1)./(x(idx)-aa) - (b-1)./(bb-x(idx)) ;
|
||||||
else
|
else
|
||||||
|
@ -51,6 +54,7 @@ end
|
||||||
|
|
||||||
|
|
||||||
if nargout == 3
|
if nargout == 3
|
||||||
|
D2ldens = ldens ;
|
||||||
if length(a)==1
|
if length(a)==1
|
||||||
D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 - (b-1)./(bb-x(idx)).^2 ;
|
D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 - (b-1)./(bb-x(idx)).^2 ;
|
||||||
else
|
else
|
||||||
|
|
|
@ -13,11 +13,13 @@ function [ldens,Dldens,D2ldens] = lpdfig1(x,s,nu)
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% ldens [double] m*n matrix of logged INVERSE-GAMMA-1 densities evaluated at x.
|
% 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
|
% SPECIAL REQUIREMENTS
|
||||||
% none
|
% none
|
||||||
|
|
||||||
% Copyright (C) 2004-2017 Dynare Team
|
% Copyright (C) 2004-2021 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -44,6 +46,7 @@ else
|
||||||
end
|
end
|
||||||
|
|
||||||
if nargout >1
|
if nargout >1
|
||||||
|
Dldens = ldens ;
|
||||||
if length(s)==1
|
if length(s)==1
|
||||||
Dldens(idx) = - (nu+1)./(x(idx)) + s./(x(idx).^3) ;
|
Dldens(idx) = - (nu+1)./(x(idx)) + s./(x(idx).^3) ;
|
||||||
else
|
else
|
||||||
|
@ -52,6 +55,7 @@ if nargout >1
|
||||||
end
|
end
|
||||||
|
|
||||||
if nargout == 3
|
if nargout == 3
|
||||||
|
D2ldens = ldens ;
|
||||||
if length(s)==1
|
if length(s)==1
|
||||||
D2ldens(idx) = (nu+1)./(x(idx).^2) - 3*s(idx)./(x(idx).^4) ;
|
D2ldens(idx) = (nu+1)./(x(idx).^2) - 3*s(idx)./(x(idx).^4) ;
|
||||||
else
|
else
|
||||||
|
|
|
@ -13,11 +13,13 @@ function [ldens,Dldens,D2ldens] = lpdfig2(x,s,nu)
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% ldens [double] m*n matrix of logged INVERSE-GAMMA-2 densities evaluated at x.
|
% 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
|
% SPECIAL REQUIREMENTS
|
||||||
% none
|
% none
|
||||||
|
|
||||||
% Copyright (C) 2004-2017 Dynare Team
|
% Copyright (C) 2004-2021 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -44,6 +46,7 @@ else
|
||||||
end
|
end
|
||||||
|
|
||||||
if nargout >1
|
if nargout >1
|
||||||
|
Dldens = ldens;
|
||||||
if length(s)==1
|
if length(s)==1
|
||||||
Dldens(idx) = - .5*(nu+2)./(x(idx)) + .5*s./x(idx).^2;
|
Dldens(idx) = - .5*(nu+2)./(x(idx)) + .5*s./x(idx).^2;
|
||||||
else
|
else
|
||||||
|
@ -52,6 +55,7 @@ if nargout >1
|
||||||
end
|
end
|
||||||
|
|
||||||
if nargout == 3
|
if nargout == 3
|
||||||
|
D2ldens = ldens;
|
||||||
if length(s)==1
|
if length(s)==1
|
||||||
D2ldens(idx) = .5*(nu+2)./(x(idx)).^2 - s./x(idx).^3;
|
D2ldens(idx) = .5*(nu+2)./(x(idx)).^2 - s./x(idx).^3;
|
||||||
else
|
else
|
||||||
|
|
|
@ -8,12 +8,13 @@ function [ldens,Dldens,D2ldens] = lpdfnorm(x,a,b)
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% 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
|
% SPECIAL REQUIREMENTS
|
||||||
% none
|
% none
|
||||||
|
|
||||||
% Copyright (C) 2003-2017 Dynare Team
|
% Copyright (C) 2003-2021 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
|
|
@ -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)
|
% - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
|
||||||
% - [ ] improve check for duplicate moments by using the cellfun and unique functions
|
% - [ ] improve check for duplicate moments by using the cellfun and unique functions
|
||||||
% - [ ] dirname option to save output to different directory not yet implemented
|
% - [ ] 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
|
% The TeX option crashes MATLAB R2014a run with "-nodisplay" option
|
||||||
% (as is done from the testsuite).
|
% (as is done from the testsuite).
|
||||||
|
@ -183,9 +182,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
|
||||||
if options_mom_.order > 3
|
if options_mom_.order > 3
|
||||||
error('method_of_moments: perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM.\n');
|
error('method_of_moments: perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM.\n');
|
||||||
end
|
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
|
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 )
|
% 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;
|
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
|
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
|
% 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
|
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
|
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;
|
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
|
else
|
||||||
options_mom_.mom.compute_derivs = false;
|
options_mom_.mom.compute_derivs = false;
|
||||||
objective_function1 = objective_function;
|
|
||||||
end
|
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_);
|
Bounds, oo_, estim_params_, M_, options_mom_);
|
||||||
if options_mom_.vector_output
|
if options_mom_.vector_output
|
||||||
fval = fval'*fval;
|
fval = fval'*fval;
|
||||||
|
|
|
@ -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_)
|
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, junk1, junk2, 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
|
% 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 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 info: vector storing error code and penalty
|
||||||
% o exit_flag: 0 if error, 1 if no error
|
% o exit_flag: 0 if error, 1 if no error
|
||||||
% o junk1: empty matrix required for optimizer interface
|
% o df: analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only)
|
||||||
% o junk2: empty matrix required for optimizer interface
|
% o junk1: empty matrix required for optimizer interface (Hessian would go here)
|
||||||
% o oo_: structure containing the results with the following updated fields:
|
% o oo_: structure containing the results with the following updated fields:
|
||||||
% - mom.model_moments [numMom x 1] vector with model moments
|
% - mom.model_moments [numMom x 1] vector with model moments
|
||||||
% - mom.Q value of the quadratic form of the moment difference
|
% - mom.Q value of the quadratic form of the moment difference
|
||||||
% o M_: Matlab's structure describing the model
|
% 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 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
|
% This function is called by
|
||||||
% o method_of_moments.m
|
% 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...
|
% 0. Initialization of the returned variables and others...
|
||||||
%------------------------------------------------------------------------------
|
%------------------------------------------------------------------------------
|
||||||
if options_mom_.vector_output == 1
|
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
|
||||||
|
if options_mom_.vector_output == 1
|
||||||
if options_mom_.mom.penalized_estimator
|
if options_mom_.mom.penalized_estimator
|
||||||
df = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1));
|
df = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1));
|
||||||
else
|
else
|
||||||
df = nan(size(oo_.mom.data_moments,1),length(xparam1));
|
df = nan(size(oo_.mom.data_moments,1),length(xparam1));
|
||||||
end
|
end
|
||||||
else
|
else
|
||||||
df = nan(1,length(xparam1));
|
df = nan(1,length(xparam1));
|
||||||
|
end
|
||||||
|
else
|
||||||
|
df=[]; %required to be empty by e.g. newrat
|
||||||
end
|
end
|
||||||
junk1 = [];
|
junk1 = [];
|
||||||
junk2 = [];
|
junk2 = [];
|
||||||
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
|
||||||
% -------------------------------------------------------------------------
|
|
||||||
% 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
|
|
||||||
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
[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
|
|
@ -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
|
% http://sims.princeton.edu/yftp/optimize/mfiles/csminwel.m
|
||||||
%
|
%
|
||||||
% Copyright (C) 1993-2007 Christopher Sims
|
% Copyright (C) 1993-2007 Christopher Sims
|
||||||
% Copyright (C) 2006-2017 Dynare Team
|
% Copyright (C) 2006-2021 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -79,9 +79,9 @@ if ischar(fcn)
|
||||||
end
|
end
|
||||||
if ischar(grad)
|
if ischar(grad)
|
||||||
grad = str2func(grad);
|
grad = str2func(grad);
|
||||||
grad_fun_provided=1;
|
grad_fun_provided = true;
|
||||||
else
|
else
|
||||||
grad_fun_provided=0;
|
grad_fun_provided = false;
|
||||||
end
|
end
|
||||||
%tailstr = ')';
|
%tailstr = ')';
|
||||||
%stailstr = [];
|
%stailstr = [];
|
||||||
|
@ -228,16 +228,16 @@ while ~done
|
||||||
f2=f;f3=f;f1=f;retcode2=retcode1;retcode3=retcode1;
|
f2=f;f3=f;f1=f;retcode2=retcode1;retcode3=retcode1;
|
||||||
end
|
end
|
||||||
%how to pick gh and xh
|
%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;
|
ih=3;
|
||||||
fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
|
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;
|
ih=2;
|
||||||
fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
|
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;
|
ih=1;
|
||||||
fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
|
fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
|
||||||
else
|
else % stuck or bad gradient
|
||||||
[fh,ih] = min([f1,f2,f3]);
|
[fh,ih] = min([f1,f2,f3]);
|
||||||
%disp_verbose(sprintf('ih = %d',ih))
|
%disp_verbose(sprintf('ih = %d',ih))
|
||||||
%eval(['xh=x' num2str(ih) ';'])
|
%eval(['xh=x' num2str(ih) ';'])
|
||||||
|
@ -253,12 +253,10 @@ while ~done
|
||||||
%eval(['retcodeh=retcode' num2str(ih) ';'])
|
%eval(['retcodeh=retcode' num2str(ih) ';'])
|
||||||
retcodei=[retcode1,retcode2,retcode3];
|
retcodei=[retcode1,retcode2,retcode3];
|
||||||
retcodeh=retcodei(ih);
|
retcodeh=retcodei(ih);
|
||||||
if exist('gh')
|
|
||||||
nogh=isempty(gh);
|
nogh=isempty(gh);
|
||||||
else
|
badgh=1;
|
||||||
nogh=1;
|
if nogh %recompute gradient
|
||||||
end
|
|
||||||
if nogh
|
|
||||||
if NumGrad
|
if NumGrad
|
||||||
[gh, badgh]=get_num_grad(method,fcn,penalty,fh,xh,epsilon,varargin{:});
|
[gh, badgh]=get_num_grad(method,fcn,penalty,fh,xh,epsilon,varargin{:});
|
||||||
elseif grad_fun_provided
|
elseif grad_fun_provided
|
||||||
|
@ -268,7 +266,6 @@ while ~done
|
||||||
badgh = ~cost_flag;
|
badgh = ~cost_flag;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
badgh=1;
|
|
||||||
end
|
end
|
||||||
%end of picking
|
%end of picking
|
||||||
stuck = (abs(fh-f) < crit);
|
stuck = (abs(fh-f) < crit);
|
||||||
|
|
|
@ -70,25 +70,34 @@ switch minimizer_algorithm
|
||||||
error('Optimization algorithm 1 requires the Optimization Toolbox')
|
error('Optimization algorithm 1 requires the Optimization Toolbox')
|
||||||
end
|
end
|
||||||
% Set default optimization options for fmincon.
|
% Set default optimization options for fmincon.
|
||||||
optim_options = optimset('display','iter', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
|
optim_options = optimoptions('fmincon','display','iter', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
|
||||||
if ~isoctave
|
|
||||||
optim_options = optimset(optim_options, 'LargeScale','off');
|
|
||||||
end
|
|
||||||
if isoctave
|
if isoctave
|
||||||
% Under Octave, use the active-set (i.e. octave_sqp of nonlin_min)
|
% Under Octave, use the active-set (i.e. octave_sqp of nonlin_min)
|
||||||
% algorithm. On a simple example (fs2000.mod), the default algorithm
|
% algorithm. On a simple example (fs2000.mod), the default algorithm
|
||||||
% is not able to even move away from the initial point.
|
% 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
|
end
|
||||||
if ~isempty(options_.optim_opt)
|
if ~isempty(options_.optim_opt)
|
||||||
eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
|
eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
|
||||||
end
|
end
|
||||||
if options_.silent_optimizer
|
if options_.silent_optimizer
|
||||||
optim_options = optimset(optim_options,'display','off');
|
optim_options = optimoptions(optim_options,'display','off');
|
||||||
end
|
end
|
||||||
if options_.analytic_derivation
|
if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1) %use wrapper
|
||||||
optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7);
|
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
|
end
|
||||||
|
else
|
||||||
if ~isoctave
|
if ~isoctave
|
||||||
[opt_par_values,fval,exitflag,output,lamdba,grad,hessian_mat] = ...
|
[opt_par_values,fval,exitflag,output,lamdba,grad,hessian_mat] = ...
|
||||||
fmincon(objective_function,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options,varargin{:});
|
fmincon(objective_function,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options,varargin{:});
|
||||||
|
@ -99,6 +108,8 @@ switch minimizer_algorithm
|
||||||
[opt_par_values,fval,exitflag,output] = ...
|
[opt_par_values,fval,exitflag,output] = ...
|
||||||
fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
|
fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
|
||||||
end
|
end
|
||||||
|
end
|
||||||
|
|
||||||
case 2
|
case 2
|
||||||
%simulating annealing
|
%simulating annealing
|
||||||
sa_options = options_.saopt;
|
sa_options = options_.saopt;
|
||||||
|
@ -157,16 +168,24 @@ switch minimizer_algorithm
|
||||||
error('Optimization algorithm 3 requires the Optimization Toolbox')
|
error('Optimization algorithm 3 requires the Optimization Toolbox')
|
||||||
end
|
end
|
||||||
% Set default optimization options for fminunc.
|
% 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)
|
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_.analytic_derivation
|
|
||||||
optim_options = optimset(optim_options,'GradObj','on');
|
|
||||||
end
|
end
|
||||||
if options_.silent_optimizer
|
if options_.silent_optimizer
|
||||||
optim_options = optimset(optim_options,'display','off');
|
optim_options = optimoptions(optim_options,'display','off');
|
||||||
end
|
end
|
||||||
|
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
|
||||||
if ~isoctave
|
if ~isoctave
|
||||||
[opt_par_values,fval,exitflag] = fminunc(objective_function,start_par_value,optim_options,varargin{:});
|
[opt_par_values,fval,exitflag] = fminunc(objective_function,start_par_value,optim_options,varargin{:});
|
||||||
else
|
else
|
||||||
|
@ -174,6 +193,7 @@ switch minimizer_algorithm
|
||||||
func = @(x) objective_function(x,varargin{:});
|
func = @(x) objective_function(x,varargin{:});
|
||||||
[opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
|
[opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
|
||||||
end
|
end
|
||||||
|
end
|
||||||
case 4
|
case 4
|
||||||
% Set default options.
|
% Set default options.
|
||||||
H0 = 1e-4*eye(n_params);
|
H0 = 1e-4*eye(n_params);
|
||||||
|
@ -480,6 +500,34 @@ switch minimizer_algorithm
|
||||||
% Minimize the penalized objective (note that the penalty is not updated).
|
% 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, fval, exitflag, output] = particleswarm(objfun, length(start_par_value), LB, UB, particleswarmOptions);
|
||||||
opt_par_values = opt_par_values(:);
|
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
|
case 101
|
||||||
solveoptoptions = options_.solveopt;
|
solveoptoptions = options_.solveopt;
|
||||||
if ~isempty(options_.optim_opt)
|
if ~isempty(options_.optim_opt)
|
||||||
|
@ -506,7 +554,12 @@ switch minimizer_algorithm
|
||||||
if options_.silent_optimizer
|
if options_.silent_optimizer
|
||||||
solveoptoptions.verbosity = 0;
|
solveoptoptions.verbosity = 0;
|
||||||
end
|
end
|
||||||
|
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{:});
|
[opt_par_values,fval]=solvopt(start_par_value,objective_function,[],[],[],solveoptoptions,varargin{:});
|
||||||
|
end
|
||||||
case 102
|
case 102
|
||||||
if isoctave
|
if isoctave
|
||||||
error('Optimization algorithm 2 is not available under Octave')
|
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')
|
error('Optimization algorithm 2 requires the Global Optimization Toolbox')
|
||||||
end
|
end
|
||||||
% Set default optimization options for simulannealbnd.
|
% 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)
|
if ~isempty(options_.optim_opt)
|
||||||
eval(['optim_options = saoptimset(optim_options,' options_.optim_opt ');']);
|
eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
|
||||||
end
|
end
|
||||||
if options_.silent_optimizer
|
if options_.silent_optimizer
|
||||||
optim_options = optimset(optim_options,'display','off');
|
optim_options = optimoptions(optim_options,'display','off');
|
||||||
end
|
end
|
||||||
func = @(x)objective_function(x,varargin{:});
|
func = @(x)objective_function(x,varargin{:});
|
||||||
[opt_par_values,fval,exitflag,output] = simulannealbnd(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
|
[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
|
otherwise
|
||||||
if ischar(minimizer_algorithm)
|
if ischar(minimizer_algorithm)
|
||||||
if exist(minimizer_algorithm)
|
if exist(minimizer_algorithm)
|
||||||
|
|
|
@ -200,7 +200,7 @@ end
|
||||||
epsnorm=1.e-15;
|
epsnorm=1.e-15;
|
||||||
epsnorm2=1.e-30; % epsilon & epsilon^2
|
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
|
cnteps=optim.TolXConstraint; % Max. admissible residual
|
||||||
else
|
else
|
||||||
h1=sign(optim.minimizer_indicator); % Minimize resp. maximize a function
|
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 ---{
|
%Display control ---{
|
||||||
if optim.verbosity<=0, dispdata=0
|
if optim.verbosity<=0, dispdata=0;
|
||||||
if optim.verbosity==-1
|
if optim.verbosity==-1
|
||||||
dispwarn=0;
|
dispwarn=0;
|
||||||
else
|
else
|
||||||
|
@ -320,7 +320,7 @@ elseif abs(f)==Inf
|
||||||
end
|
end
|
||||||
xrec=x; frec=f; % record point and function value
|
xrec=x; frec=f; % record point and function value
|
||||||
% Constrained problem
|
% Constrained problem
|
||||||
if constr, fp=f; kless=0
|
if constr, fp=f; kless=0;
|
||||||
if trx
|
if trx
|
||||||
fc=feval(func,x');
|
fc=feval(func,x');
|
||||||
else
|
else
|
||||||
|
@ -541,7 +541,7 @@ while 1
|
||||||
xx=x;
|
xx=x;
|
||||||
end
|
end
|
||||||
if kj==0
|
if kj==0
|
||||||
kd=4
|
kd=4;
|
||||||
end
|
end
|
||||||
kj=kj+1; w=-.9; h=h*2; % use large coef. of space dilation
|
kj=kj+1; w=-.9; h=h*2; % use large coef. of space dilation
|
||||||
if kj>2*kd
|
if kj>2*kd
|
||||||
|
@ -574,7 +574,7 @@ while 1
|
||||||
% RESETTING ----{
|
% RESETTING ----{
|
||||||
if kcheck>1
|
if kcheck>1
|
||||||
idx=find(abs(g)>ZeroGrad); numelem=size(idx,2);
|
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 all(abs(g1(idx))<=abs(g(idx))*grbnd) || nrmz==0
|
||||||
if dispwarn
|
if dispwarn
|
||||||
disp(wrnmes)
|
disp(wrnmes)
|
||||||
|
@ -801,7 +801,7 @@ while 1
|
||||||
n_grad_evals=n_grad_evals+1;
|
n_grad_evals=n_grad_evals+1;
|
||||||
end
|
end
|
||||||
if size(g,2)==1
|
if size(g,2)==1
|
||||||
g=g'
|
g=g';
|
||||||
end
|
end
|
||||||
ng=norm(g);
|
ng=norm(g);
|
||||||
if isnan(ng)
|
if isnan(ng)
|
||||||
|
@ -926,7 +926,7 @@ while 1
|
||||||
% DISPLAY THE CURRENT VALUES ----{
|
% DISPLAY THE CURRENT VALUES ----{
|
||||||
if k==ld
|
if k==ld
|
||||||
disp('Iter.# ..... Function ... Step Value ... Gradient Norm ');
|
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;
|
ld=k+dispdata;
|
||||||
end
|
end
|
||||||
%----}
|
%----}
|
||||||
|
@ -941,7 +941,7 @@ while 1
|
||||||
termflag=0;
|
termflag=0;
|
||||||
end
|
end
|
||||||
if knan
|
if knan
|
||||||
termflag=0
|
termflag=0;
|
||||||
end
|
end
|
||||||
if kc>=mxtc
|
if kc>=mxtc
|
||||||
termflag=0;
|
termflag=0;
|
||||||
|
|
|
@ -83,8 +83,8 @@ if nargin > 6 && initialization
|
||||||
end
|
end
|
||||||
|
|
||||||
logged_prior_density = 0.0;
|
logged_prior_density = 0.0;
|
||||||
dlprior = 0.0;
|
dlprior = zeros(1,length(x));
|
||||||
d2lprior = 0.0;
|
d2lprior = dlprior;
|
||||||
|
|
||||||
if tt1
|
if tt1
|
||||||
logged_prior_density = logged_prior_density + sum(lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1))) ;
|
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
|
return
|
||||||
end
|
end
|
||||||
if nargout==2
|
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
|
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
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
|
@ -61,29 +61,74 @@ var e_a; stderr 0.014;
|
||||||
var e_m; stderr 0.005;
|
var e_m; stderr 0.005;
|
||||||
end;
|
end;
|
||||||
|
|
||||||
steady;
|
stoch_simul(order=1,periods=200, irf=0,nomoments,noprint);
|
||||||
|
save('my_data.mat','gp_obs','gy_obs');
|
||||||
check;
|
|
||||||
|
|
||||||
estimated_params;
|
estimated_params;
|
||||||
alp, beta_pdf, 0.356, 0.02;
|
alp, 0.356;
|
||||||
bet, beta_pdf, 0.993, 0.002;
|
rho, 0.129;
|
||||||
gam, normal_pdf, 0.0085, 0.003;
|
psi, 0.65;
|
||||||
mst, normal_pdf, 1.0002, 0.007;
|
stderr e_a, 0.035449;
|
||||||
rho, beta_pdf, 0.129, 0.1;
|
stderr e_m, 0.008862;
|
||||||
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;
|
|
||||||
end;
|
end;
|
||||||
|
|
||||||
varobs gp_obs gy_obs;
|
varobs gp_obs gy_obs;
|
||||||
|
|
||||||
options_.solve_tolf = 1e-12;
|
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=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);
|
||||||
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);
|
if (isoctave && user_has_octave_forge_package('optim', '1.6')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
|
||||||
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=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
|
||||||
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);
|
%,optim = ('DerivativeCheck', 'on','FiniteDifferenceType','central')
|
||||||
//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);
|
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
|
||||||
|
|
|
@ -89,7 +89,7 @@ end;
|
||||||
estimated_params_init(use_calibration);
|
estimated_params_init(use_calibration);
|
||||||
end;
|
end;
|
||||||
|
|
||||||
@#for optimizer in [1, 3, 13]
|
@#for optimizer in [1, 3, 4, 101, 13]
|
||||||
@#if estimParams == 2 && optimizer == 13
|
@#if estimParams == 2 && optimizer == 13
|
||||||
%skip due to buggy behavior in Octave
|
%skip due to buggy behavior in Octave
|
||||||
if ~isoctave
|
if ~isoctave
|
||||||
|
@ -104,7 +104,7 @@ end;
|
||||||
, nograph
|
, nograph
|
||||||
, mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance
|
, mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance
|
||||||
%, additional_optimizer_steps = [1 3 13]
|
%, additional_optimizer_steps = [1 3 13]
|
||||||
% , optim = ('TolFun', 1e-6
|
% , optim = ('DerivativeCheck', 'on','FiniteDifferenceType','central'
|
||||||
% ,'TolX', 1e-6
|
% ,'TolX', 1e-6
|
||||||
% ,'MaxIter', 3000
|
% ,'MaxIter', 3000
|
||||||
% ,'MaxFunEvals', 1D6
|
% ,'MaxFunEvals', 1D6
|
||||||
|
|
Loading…
Reference in New Issue