diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 6a0ef6e857e85bd4ec26693b5680a70fe57626ff..4558edc862df31613f6d59ea49a21ad38ccf2c02 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -288,6 +288,8 @@ switch minimizer_algorithm
newratflag = options_.newrat.hess; %default
end
nit=options_.newrat.maxiter;
+ epsilon = options_.gradient_epsilon;
+ crit_rel = NaN;
Verbose = options_.newrat.verbosity;
Save_files = options_.newrat.Save_files;
if ~isempty(options_.optim_opt)
@@ -303,8 +305,12 @@ switch minimizer_algorithm
else
newratflag=flag;
end
+ case 'NumgradEpsilon'
+ epsilon = options_list{i,2};
case 'TolFun'
crit = options_list{i,2};
+ case 'TolFunRel'
+ crit_rel = options_list{i,2};
case 'verbosity'
Verbose = options_list{i,2};
case 'SaveFiles'
@@ -318,9 +324,10 @@ switch minimizer_algorithm
Save_files = 0;
Verbose = 0;
end
+ hess_info.crit_rel=crit_rel;
hess_info.gstep=options_.gstep;
hess_info.htol = 1.e-4;
- hess_info.h1=options_.gradient_epsilon*ones(n_params,1);
+ hess_info.h1=epsilon*ones(n_params,1);
% here we force 7th input argument (flagg) to be 0, since outer product
% gradient Hessian is handled in dynare_estimation_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,prior_information.p2,options_.gradient_epsilon,parameter_names,varargin{:}); %hessian_mat is the plain outer product gradient Hessian
diff --git a/matlab/optimization/mr_gstep.m b/matlab/optimization/mr_gstep.m
index a7f39a604152b0f2a2d0fc16f138780793031751..425a495bcb56ace2d9ddcf583312a5b604886acb 100644
--- a/matlab/optimization/mr_gstep.m
+++ b/matlab/optimization/mr_gstep.m
@@ -72,15 +72,77 @@ while i<n
gg(i)=(f1(i)'-f_1(i)')./(2.*h1(i));
hh(i) = 1/max(1.e-9,abs( (f1(i)+f_1(i)-2*f0)./(h1(i)*h1(i)) ));
if gg(i)*(hh(i)*gg(i))/2 > htol(i)
- [f0, x, ~, ~] = csminit1(func0,x,penalty,f0,gg,0,diag(hh),Verbose,varargin{:});
+ [ff, xx,fcount,retcode] = csminit1(func0,x,penalty,f0,gg,0,diag(hh),Verbose,varargin{:});
+ if retcode
+% do_nothing=true;
+ xa=transpose(linspace(x(i)/2, x(i)*3/2, 7));
+ for k=1:7
+ xh1(i)=xa(k);
+ fa(k,1) = penalty_objective_function(xh1,func0,penalty,varargin{:});
+ end
+ b=[ones(7,1) xa xa.*xa./2]\fa;
+ gg(i)=x(i)*b(3)+b(2);
+ hh(i)=1/b(3);
+ [ff2, xx2, fcount2, retcode2] = csminit1(func0,x,penalty,f0,gg,0,diag(hh),Verbose,varargin{:});
+ if ff2<ff
+ ff=ff2;
+ xx=xx2;
+ end
+ if min(fa)<ff
+ [ff, im]=min(fa);
+ xx(i)=xa(im);
+ end
+ end
ig(i)=1;
+ if not(isequal(xx , check_bounds(xx,bounds)))
+ xx = check_bounds(xx,bounds);
+ if xx(i)<x(i)
+ % lower bound
+ xx(i) = min(xx(i)+h1(i), 0.5*(xx(i)+x(i)));
+ else
+ % upper bound
+ xx(i) = max(xx(i)-h1(i), 0.5*(xx(i)+x(i)));
+ end
+ [ff,exit_flag]=penalty_objective_function(xx,func0,penalty,varargin{:});
+ if exit_flag~=1
+ disp('last step exited with bad status!')
+ elseif ff<f0
+ f0=ff;
+ x=xx;
+ end
+ else
+ % check improvement wrt predicted one
+ if abs(f0-ff) < abs(gg(i)*(hh(i)*gg(i))/2/100) || abs(x(i)-xx(i))<1.e-10
+ [ff1, xx1, fcount, retcode] = csminit1(func0,x,penalty,f0,-gg,0,diag(hh),Verbose,varargin{:});
+ if not(isequal(xx1 , check_bounds(xx1,bounds)))
+ xx1 = check_bounds(xx1,bounds);
+ if xx1(i)<x(i)
+ % lower bound
+ xx1(i) = min(xx1(i)+h1(i), 0.5*(xx1(i)+x(i)));
+ else
+ % upper bound
+ xx1(i) = max(xx1(i)-h1(i), 0.5*(xx1(i)+x(i)));
+ end
+ [ff1,exit_flag]=penalty_objective_function(xx1,func0,penalty,varargin{:});
+ if exit_flag~=1
+ disp('last step exited with bad status!')
+ end
+ end
+ if ff1<ff
+ ff=ff1;
+ xx=xx1;
+ end
+ end
+ f0=ff;
+ x=xx;
+ end
if Verbose
- fprintf(['Done for param %s = %8.4f\n'],parameter_names{i},x(i))
+ fprintf(['Done for param %s = %8.4f; f = %8.4f\n'],parameter_names{i},x(i),f0)
end
end
xh1=x;
end
- x = check_bounds(x,bounds);
+% penalty=f0;
if Save_files
save('gstep.mat','x','h1','f0')
end
@@ -96,10 +158,10 @@ function x = check_bounds(x,bounds)
inx = find(x>=bounds(:,2));
if ~isempty(inx)
- x(inx) = bounds(inx,2)-eps;
+ x(inx) = bounds(inx,2)-1.e-10;
end
inx = find(x<=bounds(:,1));
if ~isempty(inx)
- x(inx) = bounds(inx,1)+eps;
+ x(inx) = bounds(inx,1)+1.e-10;
end
diff --git a/matlab/optimization/mr_hessian.m b/matlab/optimization/mr_hessian.m
index 4cd15a90cff0d053eabe09a9f1da727cae57d7fc..03dbc548ec0f79cbb3a285c91579c6c334d0688d 100644
--- a/matlab/optimization/mr_hessian.m
+++ b/matlab/optimization/mr_hessian.m
@@ -78,11 +78,16 @@ else
end
-hess_info.h1 = min(hess_info.h1,0.5.*hmax);
+hess_info.h1 = min(hess_info.h1,0.9.*hmax);
if htol0<hess_info.htol
hess_info.htol=htol0;
end
+if not(isnan(hess_info.crit_rel))
+ htol1=hess_info.htol;
+ hess_info.htol=abs(hess_info.crit_rel*f0);
+end
+
xh1=x;
f1=zeros(size(f0,1),n);
f_1=f1;
@@ -106,22 +111,39 @@ while i<n
ic=0;
icount = 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
+ istoobig=false;
+ istoosmall=false;
+ give_it_up=false;
+ while (abs(dx(it))<0.5*hess_info.htol || abs(dx(it))>(3*hess_info.htol)) && icount<10 && ic==0 && hmax(i)>2*1.e-10 && not(give_it_up)
icount=icount+1;
+ istoobig(it)=false;
+ istoosmall(it)=false;
if abs(dx(it))<0.5*hess_info.htol
+ istoosmall(it)=true;
+ if hess_info.h1(i)==0.9*hmax(i)% || hess_info.h1(i)==0.3*abs(x(i))
+ give_it_up=true;
+ end
if abs(dx(it)) ~= 0
- 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));
+% htmp=min(max(1.e-10,0.3*abs(x(i))), 0.9*hess_info.htol/abs(dx(it))*hess_info.h1(i));
+ htmp=0.9*hess_info.htol/abs(dx(it))*hess_info.h1(i);
else
- hess_info.h1(i)=2.1*hess_info.h1(i);
+ htmp=2.1*hess_info.h1(i);
end
- 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);
+ htmp = min(htmp,0.9*hmax(i));
+ if any(h0(istoobig(1:it))) %&& htmp>=min(h0(istoobig(1:it)))
+ htmp = 0.5*(min(h0(istoobig))+max(h0(istoosmall)));
+ end
+ hess_info.h1(i) = max(htmp,1.e-10);
xh1(i)=x(i)+hess_info.h1(i);
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
end
if abs(dx(it))>(3*hess_info.htol)
- hess_info.h1(i)= hess_info.htol/abs(dx(it))*hess_info.h1(i);
- hess_info.h1(i) = max(hess_info.h1(i),1e-10);
+ istoobig(it)=true;
+ htmp= hess_info.htol/abs(dx(it))*hess_info.h1(i);
+ if any(h0(istoosmall(1:it))) %&& htmp<=max(h0(istoosmall(1:it)))
+ htmp = 0.5*(min(h0(istoobig))+max(h0(istoosmall)));
+ end
+ hess_info.h1(i) = max(htmp,1e-10);
xh1(i)=x(i)+hess_info.h1(i);
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
iter=0;
@@ -136,11 +158,27 @@ while i<n
it=it+1;
dx(it)=(fx-f0);
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))
+ if (hess_info.h1(i)<1.e-12*min(1,h2(i)) && hess_info.h1(i)<0.9*hmax(i))
ic=1;
hcheck=1;
end
end
+ if icount == 10 || hess_info.h1(i)==1.e-10
+ istoobig(it)=false;
+ istoosmall(it)=false;
+ if abs(dx(it))<0.5*hess_info.htol
+ istoosmall(it)=true;
+ end
+ if abs(dx(it))>(3*hess_info.htol)
+ istoobig(it)=true;
+ end
+ if any(istoobig) && (istoosmall(it) || istoobig(it))
+ % always better to be wrong from above, to avoid numerical noise
+ [ddx, ij]=min(dx(istoobig));
+ fx=f0+ddx;
+ hess_info.h1(i)=h0(ij);
+ end
+ end
f1(:,i)=fx;
if outer_product_gradient
if any(isnan(ffx)) || isempty(ffx)
@@ -265,4 +303,6 @@ else
hh1 = [];
end
-htol1=hhtol;
+if isnan(hess_info.crit_rel)
+ htol1=hhtol;
+end
\ No newline at end of file
diff --git a/matlab/optimization/newrat.m b/matlab/optimization/newrat.m
index a3c07db051bb281f97f0f922b8ffb453a3e4086f..9832930da8942285677fffda02d9e90bb99cba49 100644
--- a/matlab/optimization/newrat.m
+++ b/matlab/optimization/newrat.m
@@ -94,6 +94,7 @@ fval=fval0;
outer_product_gradient=1;
if isempty(hh)
+ penalty=fval0;
[dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(x,func0,penalty,flagit,htol,hess_info,bounds,prior_std,Save_files,varargin{:});
if isempty(dum)
outer_product_gradient=0;
@@ -176,9 +177,25 @@ while norm(gg)>gtol && check==0 && jit<nit
iggx(ig_pos,ig_pos) = inv( hhx(ig_pos,ig_pos) );
[fvala,x0,fc,retcode] = csminit1(func0,x0,penalty,fval,ggx,0,iggx,Verbose,varargin{:});
end
- x0 = check_bounds(x0,bounds);
+ if not(isequal(x0 , check_bounds(x0,bounds)))
+ x0 = check_bounds(x0,bounds);
+ [fvala,exit_flag]=penalty_objective_function(x0,func0,penalty,varargin{:});
+ if exit_flag==1
+ penalty=fvala;
+ else
+ disp('last step exited with bad status!')
+ end
+ end
[fvala, x0, ig] = mr_gstep(h1,x0,bounds,func0,penalty,htol0,Verbose,Save_files,gradient_epsilon, parameter_names,varargin{:});
- x0 = check_bounds(x0,bounds);
+ if not(isequal(x0 , check_bounds(x0,bounds)))
+ x0 = check_bounds(x0,bounds);
+ [fvala,exit_flag]=penalty_objective_function(x0,func0,penalty,varargin{:});
+ if exit_flag==1
+ penalty=fvala;
+ else
+ disp('last step exited with bad status!')
+ end
+ end
if Save_files
nig=[nig ig];
end