diff --git a/matlab/estimation/online/online_auxiliary_filter.m b/matlab/estimation/online/online_auxiliary_filter.m
index 9e006d6e807489f8fa9ba441537f067f9df8c6bc..f4945f491f57cf98ff706546f32519b8e976d73f 100644
--- a/matlab/estimation/online/online_auxiliary_filter.m
+++ b/matlab/estimation/online/online_auxiliary_filter.m
@@ -41,7 +41,7 @@ SimulationFolder = CheckPath('online', M_.dname);
 bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
 
 % initialization of state particles
-[~, M_, options_, oo_, ReducedForm] = solve_model_for_online_filter(true, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+[~, ~, ReducedForm] = solve_model_for_online_filter(true, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
 
 order = options_.order;
 mf0 = ReducedForm.mf0;
@@ -80,7 +80,7 @@ for i=1:number_of_particles
     info = 12042009;
     while info
         candidate = Prior.draw();
-        [info, M_, options_, oo_] = solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+        [info] = solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
         if ~info
             xparam(:,i) = candidate(:);
         end
@@ -117,8 +117,8 @@ for t=1:sample_size
     tau_tilde = zeros(1,number_of_particles);
     for i=1:number_of_particles
         % model resolution
-        [info, M_, options_, oo_, ReducedForm] = ...
-            solve_model_for_online_filter(false, fore_xparam(:,i), dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_);
+        [info, M_, ReducedForm] = ...
+            solve_model_for_online_filter(false, fore_xparam(:,i), dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
         if ~info(1)
             steadystate = ReducedForm.steadystate;
             state_variables_steady_state = ReducedForm.state_variables_steady_state;
@@ -211,8 +211,8 @@ for t=1:sample_size
             candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
             if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
                 % model resolution for new parameters particles
-                [info, M_, options_, oo_, ReducedForm] = ...
-                    solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_) ;
+                [info, M_, ReducedForm] = ...
+                    solve_model_for_online_filter(false, candidate, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
                 if ~info(1)
                     xparam(:,i) = candidate ;
                     steadystate = ReducedForm.steadystate;
@@ -359,7 +359,6 @@ for t=1:sample_size
         end 
     end
     save(sprintf('%s%sparameters_particles-%u.mat', SimulationFolder, filesep(), t), 'xparam');
-
     dyntable(options_,'', {'Parameter'; 'Lower Bound (95%)';'Mean';'Upper Bound (95%)'}, bayestopt_.name, [lb95_xparam(:,t), mean_xparam(:,t), ub95_xparam(:,t)], cellofchararraymaxlength(bayestopt_.name)+2, 10, 6)
     disp('')
 
diff --git a/matlab/nonlinear-filters/solve_model_for_online_filter.m b/matlab/nonlinear-filters/solve_model_for_online_filter.m
index 988bcb59cfc7e34ccba1381d21f10595193acf66..f93eb112508fa101c10894a911559ec4a9d97039 100644
--- a/matlab/nonlinear-filters/solve_model_for_online_filter.m
+++ b/matlab/nonlinear-filters/solve_model_for_online_filter.m
@@ -1,7 +1,7 @@
-function [info, M_, options_, oo_, ReducedForm] = ...
-    solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_)
-% [info, M_, options_, oo_, ReducedForm] = ...
-%     solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, oo_)
+function [info, M_, ReducedForm] = ...
+    solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
+% [info, M_, ReducedForm] = ...
+%     solve_model_for_online_filter(setinitialcondition, xparam1, dataset_, options_, M_, estim_params_, bayestopt_, bounds, dr , endo_steady_state, exo_steady_state, exo_det_steady_state)
 
 % Solves the dsge model for an particular parameters set.
 %
@@ -14,13 +14,14 @@ function [info, M_, options_, oo_, ReducedForm] = ...
 % - estim_params_            [struct]     Estimated parameters.
 % - bayestopt_               [struct]     Prior definition.
 % - bounds                   [struct]     Prior bounds.
-% - oo_                      [struct]     Dynare results.
+% - dr                       [structure]  Reduced form model.
+% - endo_steady_state        [vector]     steady state value for endogenous variables
+% - exo_steady_state         [vector]     steady state value for exogenous variables
+% - exo_det_steady_state     [vector]     steady state value for exogenous deterministic variables
 %
 % OUTPUTS
 % - info                     [integer]    scalar, nonzero if any problem occur when computing the reduced form.
 % - M_                       [struct]     M_ description.
-% - options_                 [struct]     Dynare options.
-% - oo_                      [struct]     Dynare results.
 % - ReducedForm              [struct]     Reduced form model.
 
 % Copyright © 2013-2024 Dynare Team
@@ -40,21 +41,23 @@ function [info, M_, options_, oo_, ReducedForm] = ...
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-info = zeros(4,1);
+info        = zeros(4,1);
 
 %----------------------------------------------------
 % 1. Get the structural parameters & define penalties
 %----------------------------------------------------
 
+% Ensure that xparam1 is a column vector.
+% (Don't do the transformation if xparam1 is empty, otherwise it would become a
+%  0×1 matrix, which create issues with older MATLABs when comparing with [] in
+%  check_bounds_and_definiteness_estimation)
 if ~isempty(xparam1)
     xparam1 = xparam1(:);
 end
 
 M_ = set_all_parameters(xparam1,estim_params_,M_);
 
-
 [~,info,~,Q,H]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, bounds);
-
 if info(1)    
     return
 end
@@ -68,20 +71,17 @@ M_.H = H;
 %------------------------------------------------------------------------------
 
 warning('off', 'MATLAB:nearlySingularMatrix')
-[oo_.dr, info, M_.params] = ...
-    compute_decision_rules(M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+[dr, info, M_.params] = ...
+    compute_decision_rules(M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
 warning('on', 'MATLAB:nearlySingularMatrix')
 
 if info(1)~=0
-    if nargout==5
+    if nargout==3
         ReducedForm = 0;
     end
     return
 end
 
-% Get decision rules and transition equations.
-dr = oo_.dr;
-
 % Set persistent variables (first call).
 mf0 = bayestopt_.mf0;
 mf1 = bayestopt_.mf1;
@@ -91,7 +91,7 @@ number_of_state_variables = length(mf0);
 
 
 % Return reduced form model.
-if nargout>4
+if nargout>2
     ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
     ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
     ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
@@ -133,16 +133,11 @@ if setinitialcondition
         StateVectorVariance = StateVectorVariance(mf0,mf0);
       case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
         StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
-        old_DynareOptionsperiods = options_.periods;
         options_.periods = 5000;
-        old_DynareOptionspruning =  options_.pruning;
         options_.pruning = options_.particle.pruning;
         y_ = simult(dr.ys, dr, M_, options_);
         y_ = y_(dr.order_var(state_variables_idx),2001:options_.periods);
         StateVectorVariance = cov(y_');
-        options_.periods = old_DynareOptionsperiods;
-        options_.pruning = old_DynareOptionspruning;
-        clear('old_DynareOptionsperiods','y_');
       case 3% Initial state vector covariance is a diagonal matrix.
         StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
         StateVectorVariance = options_.particle.initial_state_prior_std*eye(number_of_state_variables);