diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m
index 68d90928c792d8a40b942f3c855dde8a7ba5c615..4b3153e9318b77efd49ae7bc24e52219f352512c 100644
--- a/matlab/+occbin/set_default_options.m
+++ b/matlab/+occbin/set_default_options.m
@@ -67,7 +67,7 @@ end
 
 if ismember(flag,{'likelihood','all'})
     options_occbin_.likelihood.curb_retrench = false;
-    options_occbin_.likelihood.first_period_occbin_update = true;
+    options_occbin_.likelihood.first_period_occbin_update = 1;
     options_occbin_.likelihood.full_output = false;
     options_occbin_.likelihood.IF_likelihood = false;
     options_occbin_.likelihood.init_regime_history = [];
@@ -187,7 +187,7 @@ if ismember(flag,{'smoother','all'})
     options_occbin_.smoother.curb_retrench = false;
     options_occbin_.smoother.debug = false;
     options_occbin_.smoother.fast = false;
-    options_occbin_.smoother.first_period_occbin_update = true;
+    options_occbin_.smoother.first_period_occbin_update = 1;
     options_occbin_.smoother.full_output = false;
 %     options.occbin.smoother.init_mode = 1; % 0 = standard;  1 = unconditional frcsts zero shocks+smoothed states in each period
     options_occbin_.smoother.init_regime_history = [];
diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index f5f0c1446c79ff6edf9b1c829e8b98c90f2216ca..2f7834df9d589ea14b7c7f497f30411d432eaad2 100644
--- a/matlab/+occbin/solve_one_constraint.m
+++ b/matlab/+occbin/solve_one_constraint.m
@@ -105,7 +105,7 @@ end
 SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods);
 SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods);
 SS_out.C = nan(DM.n_vars,n_shocks_periods);
-if ~exist('regime_history_','var') || isempty(regime_history_guess)
+if isempty(regime_history_guess)
     regime_history = struct();
     guess_history = false;
 else
@@ -152,9 +152,7 @@ for shock_period = 1:n_shocks_periods
         if length(binding_indicator)<(nperiods_0 + 1)
             binding_indicator=[binding_indicator; false(nperiods_0 + 1-length(binding_indicator),1)];
         end
-        
-        binding_indicator_history{iter}=binding_indicator;
-        
+               
         if iter==1 && guess_history_it
             regime = regime_history_guess(shock_period).regime;
             regime_start = regime_history_guess(shock_period).regimestart;
@@ -164,6 +162,7 @@ for shock_period = 1:n_shocks_periods
             end
             nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required
         end
+        binding_indicator_history{iter}=binding_indicator;
         % analyze when each regime starts based on current guess
         [regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug);
         regime_history(shock_period).regime = regime;
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index 80744ac50943fa24b2e3ee285cb7ad3d3c587152..eab9c269b6614321f7e9cbd40a36567957834c52 100644
--- a/matlab/+occbin/solve_two_constraints.m
+++ b/matlab/+occbin/solve_two_constraints.m
@@ -115,7 +115,7 @@ end
 SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods);
 SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods);
 SS_out.C = NaN(DM.n_vars,n_shocks_periods);
-if ~exist('regime_history_','var') || isempty(regime_history_guess)
+if isempty(regime_history_guess)
     regime_history = struct();
     guess_history = false;
 else
@@ -164,7 +164,6 @@ for shock_period = 1:n_shocks_periods
             % unconstrained regime (nperiods_0)
             binding_indicator=[binding_indicator; false(nperiods_0 + 1-size(binding_indicator,1),2)];
         end
-        binding_indicator_history{iter}=binding_indicator;
         
         if iter==1 && guess_history_it
             regime_1 = regime_history_guess(shock_period).regime1;
@@ -181,6 +180,8 @@ for shock_period = 1:n_shocks_periods
             end
             nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required
         end
+        binding_indicator_history{iter}=binding_indicator;
+        
         % analyse violvec and isolate contiguous periods in the other regime.
         [regime_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug);
         regime_history(shock_period).regime1 = regime_1;
diff --git a/matlab/+occbin/solver.m b/matlab/+occbin/solver.m
index 11ac70eafa860dc806753c2fa63a3e9bb7b2d1fc..e38293bcf01d637ce9b6e0829573cb6281b5ba6a 100644
--- a/matlab/+occbin/solver.m
+++ b/matlab/+occbin/solver.m
@@ -54,6 +54,9 @@ if solve_dr
     if isempty(options_.qz_criterium)
         options_.qz_criterium = 1+1e-6;
     end
+    if options_.order>1
+        options_.order = 1;
+    end
 
     [dr,error_flag,M_,oo_] = resol(0,M_,options_,oo_);
     out.error_flag=error_flag;
diff --git a/matlab/evaluate_likelihood.m b/matlab/evaluate_likelihood.m
index f48385a843e286c0e2ed0c0cecf1710514ac50f1..692588145c71fe566499197d5f8729ee3bf21f30 100644
--- a/matlab/evaluate_likelihood.m
+++ b/matlab/evaluate_likelihood.m
@@ -47,24 +47,24 @@ end
 
 if ischar(parameters)
     switch parameters
-      case 'posterior mode'
-        parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
-      case 'posterior mean'
-        parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
-      case 'posterior median'
-        parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
-      case 'prior mode'
-        parameters = bayestopt_.p5(:);
-      case 'prior mean'
-        parameters = bayestopt_.p1;
-      otherwise
-        disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:')
-        disp('                   ''posterior mode'', ')
-        disp('                   ''posterior mean'', ')
-        disp('                   ''posterior median'', ')
-        disp('                   ''prior mode'' or')
-        disp('                   ''prior mean''.')
-        error
+        case 'posterior mode'
+            parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
+        case 'posterior mean'
+            parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
+        case 'posterior median'
+            parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
+        case 'prior mode'
+            parameters = bayestopt_.p5(:);
+        case 'prior mean'
+            parameters = bayestopt_.p1;
+        otherwise
+            disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:')
+            disp('                   ''posterior mode'', ')
+            disp('                   ''posterior mean'', ')
+            disp('                   ''posterior median'', ')
+            disp('                   ''prior mode'' or')
+            disp('                   ''prior mean''.')
+            error
     end
 end
 
@@ -73,10 +73,23 @@ if isempty(dataset)
 end
 options_=select_qz_criterium_value(options_);
 
+if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
+    % Plot prior densities.
+    % Set prior bounds
+    bounds = prior_bounds(bayestopt_, options_.prior_trunc);
+else  % estimated parameters but no declared priors
+    % No priors are declared so Dynare will estimate the model by
+    % maximum likelihood with inequality constraints for the parameters.
+    [~,~,~,lb,ub] = set_prior(estim_params_,M_,options_);
+    bounds.lb = lb;
+    bounds.ub = ub;
+end
+
+
 if options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter
-    llik = -occbin.IVF_posterior(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_);
-else    
-    llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_);
+    llik = -occbin.IVF_posterior(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
+else
+    llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
 end
 ldens = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_);
 llik = llik - ldens;
\ No newline at end of file