diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 6e71390b9dd78a0eb30ae616bd61ecfc05fae832..1442ab9ff91c864fef2f0d26984f507d3daa9766 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -242,7 +242,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1),
% with flag==0, we force to use the hessian from outer
% product gradient of optimizer 5
- hh = reshape(hessian(objective_function,xparam1, ...
+ hh = reshape(penalty_hessian(objective_function,xparam1,fval, ...
options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_),nx,nx);
elseif isnumeric(options_.mode_compute) && isequal(options_.mode_compute,5)
% other numerical hessian options available with optimizer 5
@@ -269,7 +269,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
if compute_hessian,
crit = options_.newrat.tolerance.f;
newratflag = newratflag>0;
- hh = reshape(mr_hessian(0,xparam1,objective_function,newratflag,crit,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx);
+ hh = reshape(mr_hessian(xparam1,objective_function,newratflag,crit,new_rat_hess_info,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx);
end
options_.kalman_algo = kalman_algo0;
end
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 87ad5445445cb52cd8f4abf3b708ff0f512d670a..bf6904f7a69376567a5dbe8764557d9932a0f073 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -1,5 +1,5 @@
-function [opt_par_values,fval,exitflag,hessian_mat,options_,Scale]=dynare_minimize_objective(objective_function,start_par_value,minimizer_algorithm,options_,bounds,parameter_names,prior_information,Initial_Hessian,varargin)
-
+function [opt_par_values,fval,exitflag,hessian_mat,options_,Scale,new_rat_hess_info]=dynare_minimize_objective(objective_function,start_par_value,minimizer_algorithm,options_,bounds,parameter_names,prior_information,Initial_Hessian,varargin)
+% function [opt_par_values,fval,exitflag,hessian_mat,options_,Scale,new_rat_hess_info]=dynare_minimize_objective(objective_function,start_par_value,minimizer_algorithm,options_,bounds,parameter_names,prior_information,Initial_Hessian,new_rat_hess_info,varargin)
% Calls a minimizer
%
% INPUTS
@@ -11,6 +11,7 @@ function [opt_par_values,fval,exitflag,hessian_mat,options_,Scale]=dynare_minimi
% parameter_names [n_params by 1] cell array strings containing the parameters names
% prior_information [matlab structure] Dynare prior information structure (bayestopt_) provided for algorithm 6
% Initial_Hessian [n_params by n_params] matrix initial hessian matrix provided for algorithm 6
+% new_rat_hess_info [matlab structure] step size info used by algorith 5
% varargin [cell array] Input arguments for objective function
%
% OUTPUTS
@@ -25,7 +26,7 @@ function [opt_par_values,fval,exitflag,hessian_mat,options_,Scale]=dynare_minimi
% none.
%
%
-% Copyright (C) 2014-2015 Dynare Team
+% Copyright (C) 2014-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -59,7 +60,7 @@ Scale=[];
exitflag=1;
fval=NaN;
opt_par_values=NaN(size(start_par_value));
-
+new_rat_hess_info=[];
switch minimizer_algorithm
case 1
@@ -244,7 +245,10 @@ switch minimizer_algorithm
Save_files = 0;
Verbose = 0;
end
- [opt_par_values,hessian_mat,gg,fval,invhess] = newrat(objective_function,start_par_value,bounds,analytic_grad,crit,nit,0,Verbose, Save_files,varargin{:});
+ hess_info.gstep=options_.gstep;
+ hess_info.htol = 1.e-4;
+ hess_info.h1=options_.gradient_epsilon*ones(n_params,1);
+ [opt_par_values,hessian_mat,gg,fval,invhess,new_rat_hess_info] = newrat(objective_function,start_par_value,bounds,analytic_grad,crit,nit,0,Verbose, Save_files,hess_info,varargin{:});
%hessian_mat is the plain outer product gradient Hessian
case 6
[opt_par_values, hessian_mat, Scale, fval] = gmhmaxlik(objective_function, start_par_value, ...
diff --git a/matlab/optimization/mr_hessian.m b/matlab/optimization/mr_hessian.m
index 096584e90b34c18fbfa7bbc75064879c6a221c71..eb1c0752e6e594ed9314c47b83ea67a9bb2dff9f 100644
--- a/matlab/optimization/mr_hessian.m
+++ b/matlab/optimization/mr_hessian.m
@@ -1,37 +1,48 @@
-function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,penalty,hflag,htol0,varargin)
-% [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,penalty,hflag,htol0,varargin)
-%
+function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,func,penalty,hflag,htol0,hess_info,varargin)
+% function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,func,penalty,hflag,htol0,hess_info,varargin)
% numerical gradient and Hessian, with 'automatic' check of numerical
% error
%
% adapted from Michel Juillard original routine hessian.m
%
-% func = function handle. The function must give two outputs:
-% - the log-likelihood AND the single contributions at times t=1,...,T
-% of the log-likelihood to compute outer product gradient
-% x = parameter values
-% hflag = 0, Hessian computed with outer product gradient, one point
-% increments for partial derivatives in gradients
-% hflag = 1, 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
-% with correlation structure as from outer product gradient;
-% two point evaluation of derivatives for partial derivatives
-% in gradients
-% hflag = 2, full numerical Hessian, computes second order partial derivatives
-% uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27
-% p. 884.
-% htol0 = 'precision' of increment of function values for numerical
-% derivatives
-%
-% varargin{1} --> DynareDataset
-% varargin{2} --> DatasetInfo
-% varargin{3} --> DynareOptions
-% varargin{4} --> Model
-% varargin{5} --> EstimatedParameters
-% varargin{6} --> BayesInfo
-% varargin{7} --> BayesInfo
-% varargin{8} --> DynareResults
-
-
+% Inputs:
+% - func function handle. The function must give two outputs:
+% the log-likelihood AND the single contributions at times t=1,...,T
+% of the log-likelihood to compute outer product gradient
+% - x parameter values
+% - penalty penalty due to error code
+% - hflag 0: Hessian computed with outer product gradient, one point
+% increments for partial derivatives in gradients
+% 1: 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
+% with correlation structure as from outer product gradient;
+% two point evaluation of derivatives for partial derivatives
+% in gradients
+% 2: full numerical Hessian, computes second order partial derivatives
+% uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27
+% p. 884.
+% - htol0 'precision' of increment of function values for numerical
+% derivatives
+% - hess_info structure storing the step sizes for
+% computation of Hessian
+% - varargin other inputs:
+% varargin{1} --> DynareDataset
+% varargin{2} --> DatasetInfo
+% varargin{3} --> DynareOptions
+% varargin{4} --> Model
+% varargin{5} --> EstimatedParameters
+% varargin{6} --> BayesInfo
+% varargin{7} --> Bounds
+% varargin{8} --> DynareResults
+%
+% Outputs
+% - hessian_mat hessian
+% - gg Jacobian
+% - htol1 updated 'precision' of increment of function values for numerical
+% derivatives
+% - ihh inverse outer product with modified std's
+% - hh_mat0 outer product hessian with modified std's
+% - hh1 updated hess_info.h1
+% - hess_info structure with updated step length
% Copyright (C) 2004-2016 Dynare Team
%
@@ -50,15 +61,7 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,pe
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-persistent h1 htol
-
n=size(x,1);
-if init
- gstep_=varargin{3}.gstep;
- htol = 1.e-4;
- h1=varargin{3}.gradient_epsilon*ones(n,1);
- return
-end
[f0,exit_flag, ff0]=penalty_objective_function(x,func,penalty,varargin{:});
h2=varargin{7}.ub-varargin{7}.lb;
@@ -71,10 +74,10 @@ else
end
-h1 = min(h1,0.5.*hmax);
+hess_info.h1 = min(hess_info.h1,0.5.*hmax);
-if htol0<htol
- htol=htol0;
+if htol0<hess_info.htol
+ hess_info.htol=htol0;
end
xh1=x;
f1=zeros(size(f0,1),n);
@@ -86,13 +89,13 @@ if outer_product_gradient
end
i=0;
-hhtol=htol*ones(n,1);
+hhtol=hess_info.htol*ones(n,1);
while i<n
i=i+1;
- htol=hhtol(i);
- h10=h1(i);
+ hess_info.htol=hhtol(i);
+ h10=hess_info.h1(i);
hcheck=0;
- xh1(i)=x(i)+h1(i);
+ xh1(i)=x(i)+hess_info.h1(i);
try
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
catch
@@ -102,43 +105,43 @@ while i<n
dx=(fx-f0);
ic=0;
icount = 0;
- h0=h1(i);
- while (abs(dx(it))<0.5*htol || abs(dx(it))>(3*htol)) && icount<10 && ic==0
+ h0=hess_info.h1(i);
+ while (abs(dx(it))<0.5*hess_info.htol || abs(dx(it))>(3*hess_info.htol)) && icount<10 && ic==0
icount=icount+1;
- if abs(dx(it))<0.5*htol
+ if abs(dx(it))<0.5*hess_info.htol
if abs(dx(it)) ~= 0,
- h1(i)=min(max(1.e-10,0.3*abs(x(i))), 0.9*htol/abs(dx(it))*h1(i));
+ hess_info.h1(i)=min(max(1.e-10,0.3*abs(x(i))), 0.9*hess_info.htol/abs(dx(it))*hess_info.h1(i));
else
- h1(i)=2.1*h1(i);
+ hess_info.h1(i)=2.1*hess_info.h1(i);
end
- h1(i) = min(h1(i),0.5*hmax(i));
- h1(i) = max(h1(i),1.e-10);
- xh1(i)=x(i)+h1(i);
+ hess_info.h1(i) = min(hess_info.h1(i),0.5*hmax(i));
+ hess_info.h1(i) = max(hess_info.h1(i),1.e-10);
+ xh1(i)=x(i)+hess_info.h1(i);
try
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
catch
fx=1.e8;
end
end
- if abs(dx(it))>(3*htol)
- h1(i)= htol/abs(dx(it))*h1(i);
- xh1(i)=x(i)+h1(i);
+ if abs(dx(it))>(3*hess_info.htol)
+ hess_info.h1(i)= hess_info.htol/abs(dx(it))*hess_info.h1(i);
+ xh1(i)=x(i)+hess_info.h1(i);
try
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
catch
fx=1.e8;
end
while (fx-f0)==0
- h1(i)= h1(i)*2;
- xh1(i)=x(i)+h1(i);
+ hess_info.h1(i)= hess_info.h1(i)*2;
+ xh1(i)=x(i)+hess_info.h1(i);
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
ic=1;
end
end
it=it+1;
dx(it)=(fx-f0);
- h0(it)=h1(i);
- if (h1(i)<1.e-12*min(1,h2(i)) && h1(i)<0.5*hmax(i))
+ h0(it)=hess_info.h1(i);
+ if (hess_info.h1(i)<1.e-12*min(1,h2(i)) && hess_info.h1(i)<0.5*hmax(i))
ic=1;
hcheck=1;
end
@@ -151,7 +154,7 @@ while i<n
ff1=ffx;
end
end
- xh1(i)=x(i)-h1(i);
+ xh1(i)=x(i)-hess_info.h1(i);
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
f_1(:,i)=fx;
if outer_product_gradient,
@@ -160,42 +163,42 @@ while i<n
else
ff_1=ffx;
end
- ggh(:,i)=(ff1-ff_1)./(2.*h1(i));
+ ggh(:,i)=(ff1-ff_1)./(2.*hess_info.h1(i));
end
xh1(i)=x(i);
- if hcheck && htol<1
- htol=min(1,max(min(abs(dx))*2,htol*10));
- h1(i)=h10;
- hhtol(i) = htol;
+ if hcheck && hess_info.htol<1
+ hess_info.htol=min(1,max(min(abs(dx))*2,hess_info.htol*10));
+ hess_info.h1(i)=h10;
+ hhtol(i) = hess_info.htol;
i=i-1;
end
end
-h_1=h1;
+h_1=hess_info.h1;
xh1=x;
xh_1=xh1;
-gg=(f1'-f_1')./(2.*h1);
+gg=(f1'-f_1')./(2.*hess_info.h1);
if outer_product_gradient,
if hflag==2
- gg=(f1'-f_1')./(2.*h1);
+ gg=(f1'-f_1')./(2.*hess_info.h1);
hessian_mat = zeros(size(f0,1),n*n);
for i=1:n
if i > 1
k=[i:n:n*(i-1)];
hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
end
- hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
+ hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(hess_info.h1(i)*h_1(i));
temp=f1+f_1-f0*ones(1,n);
for j=i+1:n
- xh1(i)=x(i)+h1(i);
+ xh1(i)=x(i)+hess_info.h1(i);
xh1(j)=x(j)+h_1(j);
- xh_1(i)=x(i)-h1(i);
+ xh_1(i)=x(i)-hess_info.h1(i);
xh_1(j)=x(j)-h_1(j);
temp1 = penalty_objective_function(xh1,func,penalty,varargin{:});
temp2 = penalty_objective_function(xh_1,func,penalty,varargin{:});
- hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
+ hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*hess_info.h1(i)*h_1(j));
xh1(i)=x(i);
xh1(j)=x(j);
xh_1(i)=x(i);
@@ -207,7 +210,7 @@ if outer_product_gradient,
elseif hflag==1
hessian_mat = zeros(size(f0,1),n*n);
for i=1:n
- dum = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
+ dum = (f1(:,i)+f_1(:,i)-2*f0)./(hess_info.h1(i)*h_1(i));
if dum>eps
hessian_mat(:,(i-1)*n+i)=dum;
else
@@ -216,10 +219,10 @@ if outer_product_gradient,
end
end
- gga=ggh.*kron(ones(size(ff1)),2.*h1'); % re-scaled gradient
+ gga=ggh.*kron(ones(size(ff1)),2.*hess_info.h1'); % re-scaled gradient
hh_mat=gga'*gga; % rescaled outer product hessian
hh_mat0=ggh'*ggh; % outer product hessian
- A=diag(2.*h1); % rescaling matrix
+ A=diag(2.*hess_info.h1); % rescaling matrix
% igg=inv(hh_mat); % inverted rescaled outer product hessian
ihh=A'*(hh_mat\A); % inverted outer product hessian
if hflag>0 && min(eig(reshape(hessian_mat,n,n)))>0
@@ -257,7 +260,7 @@ if outer_product_gradient,
ihh=hh_mat0;
hessian_mat=hh_mat0(:);
end
- hh1=h1;
+ hh1=hess_info.h1;
save hess.mat hessian_mat
else
hessian_mat=[];
diff --git a/matlab/optimization/newrat.m b/matlab/optimization/newrat.m
index d21ddd87e0985228849c8e928ff998477402126e..8085316d057f8d280720a51263c1f485fdd3d680 100644
--- a/matlab/optimization/newrat.m
+++ b/matlab/optimization/newrat.m
@@ -1,34 +1,43 @@
-function [xparam1, hh, gg, fval, igg] = newrat(func0, x, bounds, analytic_derivation, ftol0, nit, flagg, Verbose, Save_files, varargin)
-% [xparam1, hh, gg, fval, igg] = newrat(func0, x, bounds, analytic_derivation, ftol0, nit, flagg, Verbose, Save_files, varargin)
+function [xparam1, hh, gg, fval, igg, hess_info] = newrat(func0, x, bounds, analytic_derivation, ftol0, nit, flagg, Verbose, Save_files, hess_info, varargin)
+% [xparam1, hh, gg, fval, igg, hess_info] = newrat(func0, x, bounds, analytic_derivation, ftol0, nit, flagg, Verbose, Save_files, hess_info, varargin)
%
% Optimiser with outer product gradient and with sequences of univariate steps
% uses Chris Sims subroutine for line search
%
-% func0 = name of the function
-% there must be a version of the function called [func0,'_hh.m'], that also
-% gives as second OUTPUT the single contributions at times t=1,...,T
-% of the log-likelihood to compute outer product gradient
-%
-% x = starting guess
-% analytic_derivation = 1 if analytic derivs
-% ftol0 = ending criterion for function change
-% nit = maximum number of iterations
-%
-% In each iteration, Hessian is computed with outer product gradient.
-% for final Hessian (to start Metropolis):
-% flagg = 0, final Hessian computed with outer product gradient
-% flagg = 1, final 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
-% with correlation structure as from outer product gradient,
-% flagg = 2, full numerical Hessian
-%
-% varargin{1} --> DynareDataset
-% varargin{2} --> DatasetInfo
-% varargin{3} --> DynareOptions
-% varargin{4} --> Model
-% varargin{5} --> EstimatedParameters
-% varargin{6} --> BayesInfo
-% varargin{1} --> DynareResults
-
+% Inputs:
+% - func0 name of the function that also outputs the single contributions at times t=1,...,T
+% of the log-likelihood to compute outer product gradient
+% - x starting guess
+% - analytic_derivation 1 if analytic derivatives, 0 otherwise
+% - ftol0 termination criterion for function change
+% - nit maximum number of iterations
+% - flagg Indicator how to compute final Hessian (In each iteration, Hessian is computed with outer product gradient)
+% 0: final Hessian computed with outer product gradient
+% 1: final 'mixed' Hessian: diagonal elements computed with
+% numerical second order derivatives with correlation structure
+% as from outer product gradient
+% 2: full numerical Hessian
+% - Verbose 1 if explicit output is requested
+% - Save_files 1 if intermediate output is to be saved
+% - hess_info structure storing the step sizes for
+% computation of Hessian
+% - varargin other inputs:
+% varargin{1} --> DynareDataset
+% varargin{2} --> DatasetInfo
+% varargin{3} --> DynareOptions
+% varargin{4} --> Model
+% varargin{5} --> EstimatedParameters
+% varargin{6} --> BayesInfo
+% varargin{7} --> Bounds
+% varargin{8} --> DynareResults
+%
+% Outputs
+% - xparam1 parameter vector at optimum
+% - hh hessian
+% - gg gradient
+% - fval function value
+% - igg inverted outer product hessian
+% - hess_info structure with updated step length
% Copyright (C) 2004-2016 Dynare Team
%
@@ -76,8 +85,7 @@ fval=fval0;
outer_product_gradient=1;
if isempty(hh)
- mr_hessian(1,x,[],[],[],[],varargin{:});
- [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,penalty,flagit,htol,varargin{:});
+ [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(x,func0,penalty,flagit,htol,hess_info,varargin{:});
if isempty(dum),
outer_product_gradient=0;
igg = 1e-4*eye(nx);
@@ -195,7 +203,7 @@ while norm(gg)>gtol && check==0 && jit<nit
if flagit==2
hh=hh0;
elseif flagg>0
- [dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func0,penalty,flagg,ftol0,varargin{:});
+ [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagg,ftol0,hess_info,varargin{:});
if flagg==2
hh = reshape(dum,nx,nx);
ee=eig(hh);
@@ -235,7 +243,7 @@ while norm(gg)>gtol && check==0 && jit<nit
save('m1.mat','x','fval0','nig')
end
end
- [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func0,penalty,flagit,htol,varargin{:});
+ [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagit,htol,hess_info,varargin{:});
if isempty(dum),
outer_product_gradient=0;
end