diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index d656ac6b69a12603077883eb86ce21404882e4d3..8ecb60e92377d485894c026353cd276014faafab 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -513,7 +513,7 @@ newrat.tolerance.f=1e-5;
 newrat.tolerance.f_analytic=1e-7;
 newrat.maxiter=1000;
 newrat.verbosity=1;
-newrat.Save_files=1;
+newrat.Save_files=0;
 
 options_.newrat=newrat;
 
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index d92c141c9ecd26a367a5c642f6f401230f3f2b80..ca97cd30f569daec4b0c021a452904f265f7dad0 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -265,7 +265,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,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,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 13d1eae6e0d5c80b2699ef5da3eb819e7e192d22..d8ffbb4b1a334882161cb9bb78f770567d2a7a1f 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -222,6 +222,9 @@ switch minimizer_algorithm
         csminwel1(objective_function, start_par_value, H0, analytic_grad, crit, nit, numgrad, epsilon, Verbose, Save_files, varargin{:});
     hessian_mat=inv(inverse_hessian_mat);
   case 5
+    if isempty(prior_information) %mr_hessian requires it, but can be NaN
+        prior_information.p2=NaN(n_params,1);
+    end
     if options_.analytic_derivation==-1 %set outside as code for use of analytic derivation
         analytic_grad=1;
         crit = options_.newrat.tolerance.f_analytic;
@@ -265,9 +268,11 @@ switch minimizer_algorithm
     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
+    [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,varargin{:});    %hessian_mat is the plain outer product gradient Hessian
   case 6
+    if isempty(prior_information) %Inf will be reset
+        prior_information.p2=Inf(n_params,1);
+    end
     [opt_par_values, hessian_mat, Scale, fval] = gmhmaxlik(objective_function, start_par_value, ...
                                                       Initial_Hessian, options_.mh_jscale, bounds, prior_information.p2, options_.gmhmaxlik, options_.optim_opt, varargin{:});
   case 7
diff --git a/matlab/optimization/mr_hessian.m b/matlab/optimization/mr_hessian.m
index ec8ac1e900781028153a65e5c8a5594208df96e7..d1ab01e7ba68f8fee6c8c8008f8abefddd93b197 100644
--- a/matlab/optimization/mr_hessian.m
+++ b/matlab/optimization/mr_hessian.m
@@ -1,5 +1,5 @@
-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)
+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)
 %  numerical gradient and Hessian, with 'automatic' check of numerical
 %  error
 %
@@ -24,7 +24,10 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,f
 %                       derivatives
 %  - hess_info              structure storing the step sizes for
 %                           computation of Hessian
-%  - varargin               other inputs:
+%  - bounds                 prior bounds of parameters 
+%  - prior_std              prior standard devation of parameters (can be NaN)
+%  - varargin               other inputs
+%                           e.g. in dsge_likelihood
 %                           varargin{1} --> DynareDataset
 %                           varargin{2} --> DatasetInfo
 %                           varargin{3} --> DynareOptions
@@ -64,9 +67,9 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1, hess_info] = mr_hessian(x,f
 n=size(x,1);
 
 [f0,exit_flag, ff0]=penalty_objective_function(x,func,penalty,varargin{:});
-h2=varargin{7}.ub-varargin{7}.lb;
-hmax=varargin{7}.ub-x;
-hmax=min(hmax,x-varargin{7}.lb);
+h2=bounds(:,2)-bounds(:,1);
+hmax=bounds(:,2)-x;
+hmax=min(hmax,x-bounds(:,1));
 if isempty(ff0)
     outer_product_gradient=0;
 else
@@ -240,7 +243,7 @@ if outer_product_gradient
         sd=sqrt(diag(ihh));   %standard errors
         sdh=sqrt(1./diag(hh));   %diagonal standard errors
         for j=1:length(sd)
-            sd0(j,1)=min(varargin{6}.p2(j), sd(j));  %prior std
+            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
         ihh=ihh./(sd*sd').*(sd0*sd0');  %inverse outer product with modified std's
diff --git a/matlab/optimization/newrat.m b/matlab/optimization/newrat.m
index 9f14d9f1c7b0322e0a0a7a2d0334e1b64f597c0c..e1ec4c1ed064ce1103b80a18344613fe3d9bac7a 100644
--- a/matlab/optimization/newrat.m
+++ b/matlab/optimization/newrat.m
@@ -1,4 +1,4 @@
-function [xparam1, hh, gg, fval, igg, hess_info] = newrat(func0, x, bounds, analytic_derivation, ftol0, nit, flagg, Verbose, Save_files, hess_info, varargin)
+function [xparam1, hh, gg, fval, igg, hess_info] = newrat(func0, x, bounds, analytic_derivation, ftol0, nit, flagg, Verbose, Save_files, hess_info, prior_std, 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
@@ -8,6 +8,7 @@ function [xparam1, hh, gg, fval, igg, hess_info] = newrat(func0, x, bounds, anal
 %  - 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
+%  - bounds                 prior bounds of parameters 
 %  - analytic_derivation    1 if analytic derivatives, 0 otherwise
 %  - ftol0                  termination criterion for function change
 %  - nit                    maximum number of iterations
@@ -21,7 +22,10 @@ function [xparam1, hh, gg, fval, igg, hess_info] = newrat(func0, x, bounds, anal
 %  - Save_files             1 if intermediate output is to be saved
 %  - hess_info              structure storing the step sizes for
 %                           computation of Hessian
-%  - varargin               other inputs:
+%  - prior_std              prior standard devation of parameters (can be NaN); 
+%                           passed to mr_hessian
+%  - varargin               other inputs
+%                           e.g. in dsge_likelihood and others:
 %                           varargin{1} --> DynareDataset
 %                           varargin{2} --> DatasetInfo
 %                           varargin{3} --> DynareOptions
@@ -69,7 +73,6 @@ ftol=ftol0;
 gtol=1.e-3;
 htol=htol_base;
 htol0=htol_base;
-gibbstol=length(varargin{6}.pshape)/50; %25;
 
 % force fcn, grad to function handle
 if ischar(func0)
@@ -85,7 +88,7 @@ fval=fval0;
 
 outer_product_gradient=1;
 if isempty(hh)
-    [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(x,func0,penalty,flagit,htol,hess_info,varargin{:});
+    [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(x,func0,penalty,flagit,htol,hess_info,bounds,prior_std,varargin{:});
     if isempty(dum)
         outer_product_gradient=0;
         igg = 1e-4*eye(nx);
@@ -203,7 +206,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,varargin{:});
+                [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagg,ftol0,hess_info,bounds,prior_std,varargin{:});
                 if flagg==2
                     hh = reshape(dum,nx,nx);
                     ee=eig(hh);
@@ -243,7 +246,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,varargin{:});
+            [dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(xparam1,func0,penalty,flagit,htol,hess_info,bounds,prior_std,varargin{:});
             if isempty(dum)
                 outer_product_gradient=0;
             end