diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index d656ac6b69a12603077883eb86ce21404882e4d3..835e64ffdac919de0b56c87b6e82f1d2c205448a 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -425,7 +425,7 @@ options_.smoother = false;
 options_.posterior_max_subsample_draws = 1200;
 options_.sub_draws = [];
 options_.ME_plot_tol=1e-6;
-% options_.use_mh_covariance_matrix = 0;
+options_.use_mh_covariance_matrix = false;
 options_.gradient_method = 2; %used by csminwel and newrat
 options_.gradient_epsilon = 1e-6; %used by csminwel and newrat
 options_.posterior_sampler_options.sampling_opt = []; %extended set of options for individual posterior samplers
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 7ba2347852bd4b6c0c5471f03786aefa4e37f990..963c9d4ecf459b143aa0fe3409e61ac6425491e8 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -140,7 +140,11 @@ SteadyState = [];
 trend_coeff = [];
 exit_flag   = 1;
 info        = zeros(4,1);
-DLIK        = [];
+if DynareOptions.analytic_derivation
+    DLIK        = NaN(1,length(xparam1));
+else
+    DLIK        = [];
+end
 Hess        = [];
 
 % Ensure that xparam1 is a column vector.
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index ef6b895686ba66c18479f132ffac97d3297b473c..54f9ea25a77a2761745fa263cc013d92997a7089 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -177,16 +177,22 @@ catch % if check fails, provide info on using calibration if present
 end
 
 if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
-    if options_.smoother
-        [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
-        [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
-        if options_.forecast > 0
-            oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info);
+    if options_.order==1 && ~options_.particle.status
+        if options_.smoother
+            [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
+            [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+            if options_.forecast > 0
+                oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info);
+            end
+        end
+        %reset qz_criterium
+        options_.qz_criterium=qz_criterium_old;
+        return
+    else %allow to continue, e.g. with MCMC_jumping_covariance
+        if options_.smoother
+            error('Estimation:: Particle Smoothers are not yet implemented.')
         end
     end
-    %reset qz_criterium
-    options_.qz_criterium=qz_criterium_old;
-    return
 end
 
 %% Estimation of the posterior mode or likelihood mode
@@ -533,13 +539,26 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
                 if error_flag
                     error('Estimation::mcmc: I cannot compute the posterior moments for the endogenous variables!')
                 end
+                if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results
+                   [MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
+                   if ~info
+                       generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'};
+                       for ii=1:length(generic_post_data_file_name)
+                           delete_stale_file([MetropolisFolder filesep M_.fname '_' generic_post_data_file_name{1,ii} '*']);
+                       end
+                   end
+                end
                 oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
             end
             if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
                 if error_flag
                     error('Estimation::mcmc: I cannot compute the posterior statistics!')
                 end
-                prior_posterior_statistics('posterior',dataset_,dataset_info);
+                if options_.order==1 && ~options_.particle.status
+                    prior_posterior_statistics('posterior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts
+                else
+                    error('Estimation::mcmc: Particle Smoothers are not yet implemented.')
+                end
             end
         else
             fprintf('Estimation:mcmc: sub_draws was set to 0. Skipping posterior computations.')
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index 39a45d155fb4c2c6cbe88071b4de86689fea0912..ffdc9ad7af287fecd34f7129961616ba9d60a390 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -51,6 +51,12 @@ if DynareOptions.order>1
     if BayesInfo.with_trend
         error('initial_estimation_checks:: particle filtering does not support trends')
     end
+    if DynareOptions.order>2 && DynareOptions.particle.pruning==1
+        error('initial_estimation_checks:: the particle filter with order>2 does not support pruning')
+    end
+    if DynareOptions.particle.pruning~=DynareOptions.pruning
+        fprintf('initial_estimation_checks:: the pruning settings differ between the particle filter and the one used for IRFs/simulations. Make sure this is intended.\n')
+    end
 end
 
 non_zero_ME=length(EstimatedParameters.H_entries_to_check_for_positive_definiteness);
@@ -163,6 +169,11 @@ if DynareOptions.debug
 end
 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 ~strcmp(func2str(objective_function),'dsge_likelihood')
 %         error('Options mode_compute=5 is not compatible with non linear filters or Dsge-VAR models!')
diff --git a/matlab/kalman/likelihood/kalman_filter.m b/matlab/kalman/likelihood/kalman_filter.m
index ec26524eb55139e3074a2d929361b9d1b79b90d5..14ace968ddd12545936f360be9b07a02a00a7c9e 100644
--- a/matlab/kalman/likelihood/kalman_filter.m
+++ b/matlab/kalman/likelihood/kalman_filter.m
@@ -53,7 +53,7 @@ function [LIK, LIKK, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_t
 %! @table @ @var
 %! @item LIK
 %! Double scalar, value of (minus) the likelihood.
-%! @item likk
+%! @item LIKK
 %! Column vector of doubles, values of the density of each observation.
 %! @item a
 %! Vector (@var{mm}*1) of doubles, mean of the state vector at the end of the (sub)sample.
@@ -254,9 +254,9 @@ end
 if t <= last
     if analytic_derivation
         if analytic_derivation==2
-            [tmp, tmp2] = kalman_filter_ss(Y, t, last, a, T, K, iF, dF, Z, pp, Zflag, analytic_derivation, Da, DT, DYss, D2a, D2T, D2Yss);
+            [tmp, tmp2] = kalman_filter_ss(Y, t, last, a, T, K, iF, log_dF, Z, pp, Zflag, analytic_derivation, Da, DT, DYss, D2a, D2T, D2Yss);
         else
-            [tmp, tmp2] = kalman_filter_ss(Y, t, last, a, T, K, iF, dF, Z, pp, Zflag, analytic_derivation, Da, DT, DYss, asy_hess);
+            [tmp, tmp2] = kalman_filter_ss(Y, t, last, a, T, K, iF, log_dF, Z, pp, Zflag, analytic_derivation, Da, DT, DYss, asy_hess);
         end
         likk(s+1:end) = tmp2{1};
         dlikk(s+1:end,:) = tmp2{2};
diff --git a/matlab/kalman/likelihood/kalman_filter_fast.m b/matlab/kalman/likelihood/kalman_filter_fast.m
index db3dec3b755d406d55bf89005c80b6946b3a49e6..5b463323a390c91b6d702a60fe086e7e0650ac8c 100644
--- a/matlab/kalman/likelihood/kalman_filter_fast.m
+++ b/matlab/kalman/likelihood/kalman_filter_fast.m
@@ -209,10 +209,10 @@ end
 if t <= last
     if analytic_derivation
         if analytic_derivation==2
-            [tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
+            [tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,log(dF),Z,pp,Zflag, ...
                                            analytic_derivation,Da,DT,DYss,D2a,D2T,D2Yss);
         else
-            [tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
+            [tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,log(dF),Z,pp,Zflag, ...
                                            analytic_derivation,Da,DT,DYss,asy_hess);
         end
         likk(s+1:end)=tmp2{1};
diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m
index 42d563a9bc58adc3a0c4717c7bad0adbc812cb89..8f2f96023e9b745cad8242b6e1e99ff929bc1ff7 100644
--- a/matlab/non_linear_dsge_likelihood.m
+++ b/matlab/non_linear_dsge_likelihood.m
@@ -183,7 +183,7 @@ ReducedForm.H = H;
 ReducedForm.mf0 = mf0;
 ReducedForm.mf1 = mf1;
 
-if DynareOptions.k_order_solver
+if DynareOptions.k_order_solver && ~(DynareOptions.particle.pruning && DynareOptions.order==2)
     ReducedForm.use_k_order_solver = true;
     ReducedForm.dr = dr;
 else
@@ -205,10 +205,13 @@ switch DynareOptions.particle.initialization
     StateVectorMean = ReducedForm.constant(mf0);
     old_DynareOptionsperiods = DynareOptions.periods;
     DynareOptions.periods = 5000;
+    old_DynareOptionspruning =  DynareOptions.pruning;
+    DynareOptions.pruning = DynareOptions.particle.pruning;
     y_ = simult(DynareResults.steady_state, dr,Model,DynareOptions,DynareResults);
     y_ = y_(state_variables_idx,2001:5000);
     StateVectorVariance = cov(y_');
     DynareOptions.periods = old_DynareOptionsperiods;
+    DynareOptions.pruning = old_DynareOptionspruning;
     clear('old_DynareOptionsperiods','y_');
   case 3% Initial state vector covariance is a diagonal matrix (to be used
         % if model has stochastic trends).
diff --git a/matlab/optimization/analytic_gradient_wrapper.m b/matlab/optimization/analytic_gradient_wrapper.m
new file mode 100644
index 0000000000000000000000000000000000000000..8ef31deaff85f1815695bab04a9f6396d969a7a8
--- /dev/null
+++ b/matlab/optimization/analytic_gradient_wrapper.m
@@ -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
\ No newline at end of file
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 1f8fef8afe0e389c044deddf5cb769137feab1df..f70e3687a26589200032112db95391afe9249055 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -70,9 +70,10 @@ switch minimizer_algorithm
         error('Optimization algorithm 1 requires the Optimization Toolbox')
     end
     % Set default optimization options for fmincon.
-    optim_options = optimset('display','iter', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
-    if ~isoctave
-        optim_options = optimset(optim_options, 'LargeScale','off');
+    if isoctave || matlab_ver_less_than('8.1') 
+        optim_options = optimset('display','iter', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
+    else
+        optim_options = optimoptions('fmincon','display','iter', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
     end
     if isoctave
         % Under Octave, use the active-set (i.e. octave_sqp of nonlin_min)
@@ -80,25 +81,51 @@ switch minimizer_algorithm
         % is not able to even move away from the initial point.
         optim_options = optimset(optim_options, 'Algorithm','active-set');
     end
+    if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1)
+        if isoctave || matlab_ver_less_than('8.1') 
+            optim_options = optimset('GradObj','on','TolX',1e-7);
+        else            
+            optim_options = optimoptions(optim_options,'GradObj','on','TolX',1e-7); %alter default TolX
+        end
+    end
     if ~isempty(options_.optim_opt)
-        eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
+        if isoctave || matlab_ver_less_than('8.1') 
+            eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
+        else           
+            eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
+        end
     end
     if options_.silent_optimizer
-        optim_options = optimset(optim_options,'display','off');
-    end
-    if options_.analytic_derivation
-        optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7);
+        if isoctave || matlab_ver_less_than('8.1') 
+            optim_options = optimset(optim_options,'display','off');
+        else           
+            optim_options = optimoptions(optim_options,'display','off');
+        end
     end
-    if ~isoctave
-        [opt_par_values,fval,exitflag,output,lamdba,grad,hessian_mat] = ...
-            fmincon(objective_function,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options,varargin{:});
+    if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1) %use wrapper
+        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
     else
-        % Under Octave, use a wrapper, since fmincon() does not have an 11th
-        % arg. Also, only the first 4 output arguments are available.
-        func = @(x) objective_function(x,varargin{:});
-        [opt_par_values,fval,exitflag,output] = ...
-            fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
-    end
+        if ~isoctave
+            [opt_par_values,fval,exitflag,output,lamdba,grad,hessian_mat] = ...
+                fmincon(objective_function,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options,varargin{:});
+        else
+            % Under Octave, use a wrapper, since fmincon() does not have an 11th
+            % arg. Also, only the first 4 output arguments are available.
+            func = @(x) objective_function(x,varargin{:});
+            [opt_par_values,fval,exitflag,output] = ...
+                fmincon(func,start_par_value,[],[],[],[],bounds(:,1),bounds(:,2),[],optim_options);
+        end    
+    end
+    
   case 2
     %simulating annealing
     sa_options = options_.saopt;
@@ -157,22 +184,47 @@ switch minimizer_algorithm
         error('Optimization algorithm 3 requires the Optimization Toolbox')
     end
     % Set default optimization options for fminunc.
-    optim_options = optimset('display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
-    if ~isempty(options_.optim_opt)
-        eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
+    if isoctave || matlab_ver_less_than('8.1') 
+        optim_options = optimset('display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
+    else
+        optim_options = optimoptions('fminunc','display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
     end
-    if options_.analytic_derivation
-        optim_options = optimset(optim_options,'GradObj','on');
+    if ~isempty(options_.optim_opt)
+        if isoctave || matlab_ver_less_than('8.1') 
+             eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
+        else
+            eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
+        end
     end
     if options_.silent_optimizer
-        optim_options = optimset(optim_options,'display','off');
+        if isoctave || matlab_ver_less_than('8.1') 
+            optim_options = optimset(optim_options,'display','off');
+        else
+            optim_options = optimoptions(optim_options,'display','off');
+        end
     end
-    if ~isoctave
-        [opt_par_values,fval,exitflag] = fminunc(objective_function,start_par_value,optim_options,varargin{:});
+    if options_.analytic_derivation || (isfield(options_,'mom') && options_.mom.analytic_jacobian==1)
+         if isoctave || matlab_ver_less_than('8.1') 
+            optim_options = optimset(optim_options,'GradObj','on');
+        else           
+            optim_options = optimoptions(optim_options,'GradObj','on');
+         end
+         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
-        % Under Octave, use a wrapper, since fminunc() does not have a 4th arg
-        func = @(x) objective_function(x,varargin{:});
-        [opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
+        if ~isoctave
+            [opt_par_values,fval,exitflag] = fminunc(objective_function,start_par_value,optim_options,varargin{:});
+        else
+            % Under Octave, use a wrapper, since fminunc() does not have a 4th arg
+            func = @(x) objective_function(x,varargin{:});
+            [opt_par_values,fval,exitflag] = fminunc(func,start_par_value,optim_options);
+        end
     end
   case 4
     % Set default options.
@@ -501,7 +553,12 @@ switch minimizer_algorithm
     if options_.silent_optimizer
         solveoptoptions.verbosity = 0;
     end
-    [opt_par_values,fval]=solvopt(start_par_value,objective_function,[],[],[],solveoptoptions,varargin{:});
+    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{:});
+    end
   case 102
     if isoctave
         error('Optimization algorithm 2 is not available under Octave')
@@ -509,12 +566,24 @@ switch minimizer_algorithm
         error('Optimization algorithm 2 requires the Global Optimization Toolbox')
     end
     % Set default optimization options for simulannealbnd.
-    optim_options = saoptimset('display','iter','TolFun',1e-8);
+    if isoctave || matlab_ver_less_than('8.1') 
+        optim_options = saoptimset('display','iter','TolFun',1e-8);
+    else
+        optim_options = optimoptions('simulannealbnd','display','iter','TolFun',1e-8);
+    end   
     if ~isempty(options_.optim_opt)
-        eval(['optim_options = saoptimset(optim_options,' options_.optim_opt ');']);
+        if isoctave || matlab_ver_less_than('8.1') 
+            eval(['optim_options = saoptimset(optim_options,' options_.optim_opt ');']);
+        else
+            eval(['optim_options = optimoptions(optim_options,' options_.optim_opt ');']);
+        end
     end
     if options_.silent_optimizer
-        optim_options = optimset(optim_options,'display','off');
+        if isoctave || matlab_ver_less_than('8.1') 
+            optim_options = optimset(optim_options,'display','off');
+        else
+            optim_options = optimoptions(optim_options,'display','off');
+        end
     end
     func = @(x)objective_function(x,varargin{:});
     [opt_par_values,fval,exitflag,output] = simulannealbnd(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
diff --git a/matlab/optimization/solvopt.m b/matlab/optimization/solvopt.m
index 04bf44db1ad8ea194e04e8c5c4484c597f36b746..56050b07400f615efca75f94856d956c40f1eb9e 100644
--- a/matlab/optimization/solvopt.m
+++ b/matlab/optimization/solvopt.m
@@ -14,9 +14,8 @@ function [x,f,exitflag,n_f_evals,n_grad_evals,n_constraint_evals,n_constraint_gr
 % fun       name of an M-file (M-function) which computes the value
 %           of the objective function <fun> at a point x,
 %           synopsis: f=fun(x)
-% grad      name of an M-file (M-function) which computes the gradient
+% grad      indicator whether objective function provides the gradient
 %           vector of the function <fun> at a point x,
-%           synopsis: g=grad(x)
 % func      name of an M-file (M-function) which computes the MAXIMAL
 %           RESIDUAL(!) for a set of constraints at a point x,
 %           synopsis: fc=func(x)
@@ -126,7 +125,7 @@ if nargin<2           % Function and/or starting point are not specified
 end
 if nargin<3
     app=1;             % No user-supplied gradients
-elseif isempty(grad)
+elseif isempty(grad) || grad==0
     app=1;
 else
     app=0;                     % Exact gradients are supplied
@@ -271,9 +270,19 @@ stopf=0;
 
 % COMPUTE THE FUNCTION  ( FIRST TIME ) ----{
 if trx
-    f=feval(fun,x',varargin{:});
+    if app
+        f=feval(fun,x',varargin{:});
+    else
+        [f,g]=feval(fun,x',varargin{:});
+        n_grad_evals=n_grad_evals+1;
+    end    
 else
-    f=feval(fun,x,varargin{:});
+    if app
+        f=feval(fun,x,varargin{:});
+    else
+        [f,g]=feval(fun,x,varargin{:});
+        n_grad_evals=n_grad_evals+1;
+    end
 end
 n_f_evals=n_f_evals+1;
 if isempty(f)
@@ -378,12 +387,7 @@ if app
     end
     n_f_evals=n_f_evals+n;
 else
-    if trx
-        g=feval(grad,x',varargin{:});
-    else
-        g=feval(grad,x,varargin{:});
-    end
-    n_grad_evals=n_grad_evals+1;
+    %done above
 end
 if size(g,2)==1, g=g'; end
 ng=norm(g);
@@ -674,6 +678,9 @@ while 1
                 end
                 if ksm || kc>=mxtc
                     exitflag=-3;
+                    % don't return with NaN or Inf despite error code
+                    x=x1;
+                    f=f1;                  
                     if trx
                         x=x';
                     end
@@ -786,10 +793,10 @@ while 1
             end
             n_f_evals=n_f_evals+n;
         else
-            if trx
-                g=feval(grad,x',varargin{:});
+            if trx                
+                [~,g]=feval(fun,x',varargin{:});
             else
-                g=feval(grad,x,varargin{:});
+                [~,g]=feval(fun,x,varargin{:});
             end
             n_grad_evals=n_grad_evals+1;
         end
@@ -1084,7 +1091,7 @@ while 1
                     elseif isnan(f)
                         if dispwarn
                             disp(errmes)
-                            disp(error32)
+                            disp(error31)
                         end
                         exitflag=-3;
                         if trx
@@ -1105,9 +1112,9 @@ while 1
                         n_f_evals=n_f_evals+n;
                     else
                         if trx
-                            g=feval(grad,x',varargin{:});
+                            [~,g]=feval(fun,x',varargin{:});
                         else
-                            g=feval(grad,x,varargin{:});
+                            [~,g]=feval(fun,x,varargin{:});
                         end
                         n_grad_evals=n_grad_evals+1;
                     end
@@ -1210,9 +1217,9 @@ while 1
                         n_f_evals=n_f_evals+n;
                     else
                         if trx
-                            gt=feval(grad,x1',varargin{:});
+                            [~,gt]=feval(fun,x1',varargin{:});
                         else
-                            gt=feval(grad,x1,varargin{:});
+                            [~,gt]=feval(fun,x1,varargin{:});
                         end
                         n_grad_evals=n_grad_evals+1;
                     end
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index 8f0fb2afce2bdddba7af5daa36057c0e0a188a41..912a51656a2480981d80cf870f6480ba501a3213 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -225,7 +225,6 @@ if ~options_.load_mh_file && ~options_.mh_recover
     record.LastFileNumber = AnticipatedNumberOfFiles ;
     record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
     record.MCMCConcludedSuccessfully = 0;
-    record.MCMC_sampler=options_.posterior_sampler_options.posterior_sampling_method;
     fprintf('Ok!\n');
     id = write_mh_history_file(MetropolisFolder, ModelName, record);
     disp(['Estimation::mcmc: Details about the MCMC are available in ' BaseName '_mh_history_' num2str(id) '.mat'])
@@ -285,7 +284,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
     end
     ilogpo2 = record.LastLogPost;
     ix2 = record.LastParameters;
-    [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d);
+    [d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d);
     FirstBlock = 1;
     NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
     fprintf('Estimation::mcmc: I am writing a new mh-history file... ');
@@ -367,7 +366,7 @@ elseif options_.mh_recover
         LastFileFullIndicator=1;
     end
     if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice')
-        [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d);
+        [d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d);
     end
     %% Now find out what exactly needs to be redone
     % 1. Check if really something needs to be done
@@ -488,25 +487,32 @@ elseif options_.mh_recover
     FirstBlock = find(FBlock==1,1);
 end
 
-function [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d)
-if isfield(record,'ProposalCovariance') && isfield(record,'ProposalCovariance')
-    if isfield(record,'MCMC_sampler')
-        if ~strcmp(record.MCMC_sampler,options_.posterior_sampler_options.posterior_sampling_method)
-            error(fprintf('Estimation::mcmc: The posterior_sampling_method differs from the one of the original chain. Please reset it to %s',record.MCMC_sampler))
+function [d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d)
+if ~options_.use_mh_covariance_matrix
+    if isfield(record,'ProposalCovariance') && isfield(record,'ProposalScaleVec')
+        fprintf('Estimation::mcmc: Recovering the previous proposal density.\n')
+        d=record.ProposalCovariance;
+        bayestopt_.jscale=record.ProposalScaleVec;
+    else
+        if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice')
+            % pass through input d unaltered
+            if options_.mode_compute~=0
+                fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by mode_compute\n.');
+            elseif ~isempty(options_.mode_file)
+                fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by the mode_file\n.');
+            else
+                error('Estimation::mcmc: No stored previous proposal density found, no mode-finding conducted, and no mode-file provided. I don''t know how to continue')
+            end
         end
     end
-    fprintf('Estimation::mcmc: Recovering the previous proposal density.\n')
-    d=record.ProposalCovariance;
-    bayestopt_.jscale=record.ProposalScaleVec;
 else
-    if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice')
-        % pass through input d unaltered
-        if options_.mode_compute~=0
-            fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by mode_compute\n.');
-        elseif ~isempty(options_.mode_file)
-            fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by the mode_file\n.');
-        else
-            error('Estimation::mcmc: No stored previous proposal density found, no mode-finding conducted, and no mode-file provided. I don''t know how to continue')
-        end
+    % pass through input d unaltered
+    fprintf('Estimation::mcmc: use_mh_covariance_matrix specified, continuing with proposal density implied by the previous MCMC run.\n.');
+end
+
+if isfield(record,'Sampler')
+    if ~strcmp(record.Sampler,options_.posterior_sampler_options.posterior_sampling_method)
+        warning('Estimation::mcmc: The posterior_sampling_method %s selected differs from the %s of the original chain. This may create problems with the convergence diagnostics.',options_.posterior_sampler_options.posterior_sampling_method,record.Sampler)
+        record.Sampler=options_.posterior_sampler_options.posterior_sampling_method; %update sampler used
     end
 end
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 7081d500b50a979ddbfcb930766da4a39dd21619..fde38b5ea76623b7a7ebee322b5f46ed36ec2a8a 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -26,6 +26,10 @@ if M_.exo_nbr==0
     error('stoch_simul:: does not support having no varexo in the model. As a workaround you could define a dummy exogenous variable.')
 end
 
+if ismember(options_.solve_algo,[10,11])
+    error('stoch_simul:: perturbation solutions are not compatible with mixed complementarity problems and their solvers')
+end
+
 test_for_deep_parameters_calibration(M_);
 
 dr = oo_.dr;