diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index e24ee6d6ff4582032a17873c73d988c73f0846e7..cda0ecb02694b30bb7b6385b62746f61a12d671f 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -274,7 +274,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(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,[bounds.lb bounds.ub],bayestopt_.p2,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx);
+                    hh = reshape(mr_hessian(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,[bounds.lb bounds.ub],bayestopt_.p2,0,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx);
                 end
                 options_.kalman_algo = kalman_algo0;
             end
diff --git a/matlab/optimization/csminit1.m b/matlab/optimization/csminit1.m
index 22391027ebcbce9f5dcb4d3350c94474eb8aeb3f..ec088edc724ec08ee1461ff22b0e08c6896ad034 100644
--- a/matlab/optimization/csminit1.m
+++ b/matlab/optimization/csminit1.m
@@ -129,7 +129,6 @@ else
     done=0;
     factor=3;
     shrink=1;
-    lambdaMin=0;
     lambdaMax=inf;
     lambdaPeak=0;
     fPeak=f0;
diff --git a/matlab/optimization/mr_hessian.m b/matlab/optimization/mr_hessian.m
index d1ab01e7ba68f8fee6c8c8008f8abefddd93b197..1012ee2fdc88f672528a5e2ea8f2eecfe38a4809 100644
--- a/matlab/optimization/mr_hessian.m
+++ b/matlab/optimization/mr_hessian.m
@@ -1,15 +1,15 @@
-function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,func,penalty,hflag,htol0,hess_info,bounds,prior_std,varargin)
-% function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,func,penalty,hflag,htol0,hess_info,bounds,prior_std,varargin)
+function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,func,penalty,hflag,htol0,hess_info,bounds,prior_std,Save_files,varargin)
+% function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,func,penalty,hflag,htol0,hess_info,bounds,prior_std,Save_files,varargin)
 %  numerical gradient and Hessian, with 'automatic' check of numerical
 %  error
 %
 % adapted from Michel Juillard original routine hessian.m
 %
 % Inputs:
+%  - x                  parameter values
 %  - 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
@@ -26,6 +26,7 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,f
 %                           computation of Hessian
 %  - bounds                 prior bounds of parameters 
 %  - prior_std              prior standard devation of parameters (can be NaN)
+%  - Save_files             indicator whether files should be saved
 %  - varargin               other inputs
 %                           e.g. in dsge_likelihood
 %                           varargin{1} --> DynareDataset
@@ -99,11 +100,7 @@ while i<n
     h10=hess_info.h1(i);
     hcheck=0;
     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
+    [fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
     it=1;
     dx=(fx-f0);
     ic=0;
@@ -120,21 +117,13 @@ while i<n
             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
+            [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);
             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
+            [fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
             iter=0;
             while (fx-f0)==0 && iter<50
                 hess_info.h1(i)= hess_info.h1(i)*2;
@@ -188,7 +177,7 @@ gg=(f1'-f_1')./(2.*hess_info.h1);
 
 if outer_product_gradient
     if hflag==2
-        gg=(f1'-f_1')./(2.*hess_info.h1);
+        % full numerical Hessian
         hessian_mat = zeros(size(f0,1),n*n);
         for i=1:n
             if i > 1
@@ -209,19 +198,17 @@ if outer_product_gradient
                 xh1(j)=x(j);
                 xh_1(i)=x(i);
                 xh_1(j)=x(j);
-                j=j+1;
             end
-            i=i+1;
         end
     elseif hflag==1
+        % full numerical 2nd order derivs only in diagonal
         hessian_mat = zeros(size(f0,1),n*n);
         for i=1:n
             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
-                hessian_mat(:,(i-1)*n+i)=max(eps, gg(i)^2);
-            end
+            hessian_mat(:,(i-1)*n+i)=dum;
+            if any(dum<=eps)
+                 hessian_mat(dum<=eps,(i-1)*n+i)=max(eps, gg(i)^2);
+            end                       
         end
     end
 
@@ -230,26 +217,27 @@ if outer_product_gradient
     hh_mat0=ggh'*ggh;  % outer product hessian
     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
+    ihh=A'*(hh_mat\A);  % inverted outer product hessian (based on rescaling)
     if hflag>0 && min(eig(reshape(hessian_mat,n,n)))>0
         hh0 = A*reshape(hessian_mat,n,n)*A';  %rescaled second order derivatives
-        hh = reshape(hessian_mat,n,n);  %rescaled second order derivatives
+        hh = reshape(hessian_mat,n,n);  %second order derivatives
         sd0=sqrt(diag(hh0));   %rescaled 'standard errors' using second order derivatives
         sd=sqrt(diag(hh_mat));  %rescaled 'standard errors' using outer product
         hh_mat=hh_mat./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
-        igg=inv(hh_mat);   % rescaled outer product hessian with 'true' std's
-        ihh=A'*(hh_mat\A);  % inverted outer product hessian
-        hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with 'true' std's
+        ihh=A'*(hh_mat\A);  % update inverted outer product hessian with 'true' std's 
         sd=sqrt(diag(ihh));   %standard errors
         sdh=sqrt(1./diag(hh));   %diagonal standard errors
         for j=1:length(sd)
+            % some heuristic normalizations of the standard errors that
+            % avoid numerical issues in outer product
             sd0(j,1)=min(prior_std(j), sd(j));  %prior std
             sd0(j,1)=10^(0.5*(log10(sd0(j,1))+log10(sdh(j,1))));
         end
+        inv_A=inv(A);
         ihh=ihh./(sd*sd').*(sd0*sd0');  %inverse outer product with modified std's
-        igg=inv(A)'*ihh*inv(A);  % inverted rescaled outer product hessian with modified std's
-        hh_mat=inv(igg);   % outer product rescaled hessian with modified std's
-        hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with modified std's
+        igg=inv_A'*ihh*inv_A;  % inverted rescaled outer product hessian with modified std's
+        % hh_mat=inv(igg);   % outer product rescaled hessian with modified std's
+        hh_mat0=inv_A'/igg*inv_A;  % outer product hessian with modified std's
                                         %     sd0=sqrt(1./diag(hh0));   %rescaled 'standard errors' using second order derivatives
                                         %     sd=sqrt(diag(igg));  %rescaled 'standard errors' using outer product
                                         %     igg=igg./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
@@ -267,7 +255,9 @@ if outer_product_gradient
         hessian_mat=hh_mat0(:);
     end
     hh1=hess_info.h1;
-    save hess.mat hessian_mat
+    if Save_files
+        save('hess.mat','hessian_mat')
+    end
 else
     hessian_mat=[];
     ihh=[];
diff --git a/matlab/optimization/newrat.m b/matlab/optimization/newrat.m
index e94d565540aef541026fa7203ed5df127c5f68fb..0d7ece4fee471b34a1619010dffb171713e86dfd 100644
--- a/matlab/optimization/newrat.m
+++ b/matlab/optimization/newrat.m
@@ -70,7 +70,7 @@ nx=length(x);
 xparam1=x;
 %ftol0=1.e-6;
 htol_base = max(1.e-7, ftol0);
-flagit=0;  % mode of computation of hessian in each iteration
+flagit=0;  % mode of computation of hessian in each iteration; hard-coded outer-product of gradients as it performed best in tests
 ftol=ftol0;
 gtol=1.e-3;
 htol=htol_base;
@@ -84,13 +84,17 @@ end
 % func0 = str2func([func2str(func0),'_hh']);
 % func0 = func0;
 [fval0,exit_flag,gg,hh]=penalty_objective_function(x,func0,penalty,varargin{:});
+if ~exit_flag
+    disp_verbose('Bad initial parameter.',Verbose)
+    return
+end
 fval=fval0;
 
 % initialize mr_gstep and mr_hessian
 
 outer_product_gradient=1;
 if isempty(hh)
-    [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(x,func0,penalty,flagit,htol,hess_info,bounds,prior_std,varargin{:});
+    [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;
         igg = 1e-4*eye(nx);
@@ -117,15 +121,16 @@ else
     h1=[];
 end
 H = igg;
-disp_verbose(['Gradient norm ',num2str(norm(gg))],Verbose)
-ee=eig(hh);
-disp_verbose(['Minimum Hessian eigenvalue ',num2str(min(ee))],Verbose)
-disp_verbose(['Maximum Hessian eigenvalue ',num2str(max(ee))],Verbose)
+if Verbose
+    disp_eigenvalues_gradient(gg,hh);   
+end
 g=gg;
 check=0;
-if max(eig(hh))<0
-    disp_verbose('Negative definite Hessian! Local maximum!',Verbose)
-    pause
+if Verbose
+    if max(eig(hh))<0
+        disp('Negative definite Hessian! Local maximum!')
+        pause
+    end
 end
 if Save_files
     save('m1.mat','x','hh','g','hhg','igg','fval0')
@@ -135,7 +140,9 @@ igrad=1;
 igibbs=1;
 inx=eye(nx);
 jit=0;
-nig=[];
+if Save_files
+    nig=[];
+end
 ig=ones(nx,1);
 ggx=zeros(nx,1);
 while norm(gg)>gtol && check==0 && jit<nit
@@ -156,22 +163,25 @@ while norm(gg)>gtol && check==0 && jit<nit
         fval=fval1;
         x0=x01;
     end
-    if length(find(ig))<nx
-        ggx=ggx*0;
-        ggx(find(ig))=gg(find(ig));
+    ig_pos=find(ig);
+    if length(ig_pos)<nx
+        ggx=ggx*0;        
+        ggx(ig_pos)=gg(ig_pos);
         if analytic_derivation || ~outer_product_gradient
             hhx=hh;
         else
             hhx = reshape(dum,nx,nx);
         end
         iggx=eye(length(gg));
-        iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
+        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);
     [fvala, x0, ig] = mr_gstep(h1,x0,bounds,func0,penalty,htol0,Verbose,Save_files,gradient_epsilon, parameter_names,varargin{:});
     x0 = check_bounds(x0,bounds);
-    nig=[nig ig];
+    if Save_files
+        nig=[nig ig];
+    end
     disp_verbose('Sequence of univariate steps!!',Verbose)
     fval=fvala;
     if (fval0(icount)-fval)<ftol && flagit==0
@@ -208,7 +218,7 @@ while norm(gg)>gtol && check==0 && jit<nit
             if flagit==2
                 hh=hh0;
             elseif flagg>0
-                [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagg,ftol0,hess_info,bounds,prior_std,varargin{:});
+                [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagg,ftol0,hess_info,bounds,prior_std,Save_files,varargin{:});
                 if flagg==2
                     hh = reshape(dum,nx,nx);
                     ee=eig(hh);
@@ -220,15 +230,14 @@ while norm(gg)>gtol && check==0 && jit<nit
                 end
             end
         end
-        disp_verbose(['Actual dxnorm ',num2str(norm(x(:,end)-x(:,end-1)))],Verbose)
-        disp_verbose(['FVAL          ',num2str(fval)],Verbose)
-        disp_verbose(['Improvement   ',num2str(fval0(icount)-fval)],Verbose)
-        disp_verbose(['Ftol          ',num2str(ftol)],Verbose)
-        disp_verbose(['Htol          ',num2str(max(htol0))],Verbose)
-        disp_verbose(['Gradient norm  ',num2str(norm(gg))],Verbose)
-        ee=eig(hh);
-        disp_verbose(['Minimum Hessian eigenvalue ',num2str(min(ee))],Verbose)
-        disp_verbose(['Maximum Hessian eigenvalue ',num2str(max(ee))],Verbose)
+        if Verbose
+            disp(['Actual dxnorm ',num2str(norm(x(:,end)-x(:,end-1)))])
+            disp(['FVAL          ',num2str(fval)])
+            disp(['Improvement   ',num2str(fval0(icount)-fval)])
+            disp(['Ftol          ',num2str(ftol)])
+            disp(['Htol          ',num2str(max(htol0))])
+            disp_eigenvalues_gradient(gg,hh);
+        end
         g(:,icount+1)=gg;
     else
         df = fval0(icount)-fval;
@@ -248,7 +257,7 @@ while norm(gg)>gtol && check==0 && jit<nit
                     save('m1.mat','x','fval0','nig')
                 end
             end
-            [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagit,htol,hess_info,bounds,prior_std,varargin{:});
+            [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagit,htol,hess_info,bounds,prior_std,Save_files,varargin{:});
             if isempty(dum)
                 outer_product_gradient=0;
             end
@@ -280,13 +289,11 @@ while norm(gg)>gtol && check==0 && jit<nit
             hhg=hh;
             H = inv(hh);
         end
-        disp_verbose(['Gradient norm  ',num2str(norm(gg))],Verbose)
-        ee=eig(hh);
-        disp_verbose(['Minimum Hessian eigenvalue ',num2str(min(ee))],Verbose)
-        disp_verbose(['Maximum Hessian eigenvalue ',num2str(max(ee))],Verbose)
-        if max(eig(hh))<0
-            disp_verbose('Negative definite Hessian! Local maximum!',Verbose)
-            pause(1)
+        if Verbose
+            if max(eig(hh))<0
+                disp('Negative definite Hessian! Local maximum!')
+                pause(1)
+            end
         end
         t=toc(tic1);
         disp_verbose(['Elapsed time for iteration ',num2str(t),' s.'],Verbose)
@@ -334,3 +341,10 @@ inx = find(x<=bounds(:,1));
 if ~isempty(inx)
     x(inx) = bounds(inx,1)+eps;
 end
+
+function ee=disp_eigenvalues_gradient(gg,hh)
+
+disp(['Gradient norm  ',num2str(norm(gg))])
+ee=eig(hh);
+disp(['Minimum Hessian eigenvalue ',num2str(min(ee))])
+disp(['Maximum Hessian eigenvalue ',num2str(max(ee))])