diff --git a/matlab/SMC_samplers_initialization.m b/matlab/SMC_samplers_initialization.m
index a4a90acbf..f144ab66b 100644
--- a/matlab/SMC_samplers_initialization.m
+++ b/matlab/SMC_samplers_initialization.m
@@ -85,7 +85,7 @@ fprintf(fidlog,'%% Session 1.\n');
 fprintf(fidlog,' \n');
 prior_draw(bayestopt_,options_.prior_trunc);
 % Find initial values for the NumberOfParticles chains...
-set_dynare_seed('default');
+options_=set_dynare_seed_local_options(options_,'default');
 fprintf(fidlog,['  Initial values of the parameters:\n']);
 disp('Estimation::dsmh: Searching for initial values...');
 ix2 = zeros(npar,NumberOfParticles);
diff --git a/matlab/backward/simul_backward_model_init.m b/matlab/backward/simul_backward_model_init.m
index 504593fb6..99a771127 100644
--- a/matlab/backward/simul_backward_model_init.m
+++ b/matlab/backward/simul_backward_model_init.m
@@ -155,7 +155,7 @@ if nargin<6 || isempty(innovations)
     covariance_matrix_upper_cholesky = chol(covariance_matrix);
     % Set seed to its default state.
     if DynareOptions.bnlms.set_dynare_seed_to_default
-        set_dynare_seed('default');
+        DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
     end
     % Simulate structural innovations.
     switch DynareOptions.bnlms.innovation_distribution
diff --git a/matlab/ep/extended_path_initialization.m b/matlab/ep/extended_path_initialization.m
index 6f41c4d48..4417f7b5a 100644
--- a/matlab/ep/extended_path_initialization.m
+++ b/matlab/ep/extended_path_initialization.m
@@ -97,7 +97,7 @@ end
 
 % Set seed.
 if ep.set_dynare_seed_to_default
-    set_dynare_seed('default');
+    DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 end
 
 % hybrid correction
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 99fce8305..f76c01d42 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -142,7 +142,7 @@ priordens([],[],[],[],[],[],1);
 dyn_first_order_solver();
 
 % Set dynare random generator and seed.
-set_dynare_seed('default');
+options_=set_dynare_seed_local_options(options_,'default');
 
 % Load user configuration file.
 if isfield(options_, 'global_init_file')
diff --git a/matlab/nonlinear-filters/auxiliary_initialization.m b/matlab/nonlinear-filters/auxiliary_initialization.m
index e05741052..48224451e 100644
--- a/matlab/nonlinear-filters/auxiliary_initialization.m
+++ b/matlab/nonlinear-filters/auxiliary_initialization.m
@@ -83,7 +83,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
 %end
 
 % Set seed for randn().
-set_dynare_seed('default');
+DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 
 % Initialization of the likelihood.
 const_lik = log(2*pi)*number_of_observed_variables;
diff --git a/matlab/nonlinear-filters/auxiliary_particle_filter.m b/matlab/nonlinear-filters/auxiliary_particle_filter.m
index 21caa73bd..93daa5214 100644
--- a/matlab/nonlinear-filters/auxiliary_particle_filter.m
+++ b/matlab/nonlinear-filters/auxiliary_particle_filter.m
@@ -77,7 +77,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
 Q_lower_triangular_cholesky = chol(Q)';
 
 % Set seed for randn().
-set_dynare_seed('default');
+DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 
 % Initialization of the likelihood.
 const_lik = log(2*pi)*number_of_observed_variables+log(det(H));
diff --git a/matlab/nonlinear-filters/conditional_particle_filter.m b/matlab/nonlinear-filters/conditional_particle_filter.m
index db204dc85..7a752b827 100644
--- a/matlab/nonlinear-filters/conditional_particle_filter.m
+++ b/matlab/nonlinear-filters/conditional_particle_filter.m
@@ -78,7 +78,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot, 2);
 Q_lower_triangular_cholesky = chol(Q)';
 
 % Set seed for randn().
-set_dynare_seed('default');
+DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 
 % Initialization of the likelihood.
 lik = NaN(T, 1);
diff --git a/matlab/nonlinear-filters/gaussian_filter.m b/matlab/nonlinear-filters/gaussian_filter.m
index 9b1b29e4c..afd9c87c9 100644
--- a/matlab/nonlinear-filters/gaussian_filter.m
+++ b/matlab/nonlinear-filters/gaussian_filter.m
@@ -74,7 +74,7 @@ else
 end
 
 if ParticleOptions.distribution_approximation.montecarlo
-    set_dynare_seed('default');
+    DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 end
 
 % Get covariance matrices
diff --git a/matlab/nonlinear-filters/gaussian_mixture_filter.m b/matlab/nonlinear-filters/gaussian_mixture_filter.m
index e06461b85..e28d578db 100644
--- a/matlab/nonlinear-filters/gaussian_mixture_filter.m
+++ b/matlab/nonlinear-filters/gaussian_mixture_filter.m
@@ -79,7 +79,7 @@ else
 end
 
 if ParticleOptions.distribution_approximation.montecarlo
-    set_dynare_seed('default');
+    DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 end
 
 % Get covariance matrices
diff --git a/matlab/nonlinear-filters/nonlinear_kalman_filter.m b/matlab/nonlinear-filters/nonlinear_kalman_filter.m
index ccae0731d..f741cbd77 100644
--- a/matlab/nonlinear-filters/nonlinear_kalman_filter.m
+++ b/matlab/nonlinear-filters/nonlinear_kalman_filter.m
@@ -105,7 +105,7 @@ else
 end
 
 if ParticleOptions.distribution_approximation.montecarlo
-    set_dynare_seed('default');
+    DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 end
 
 % Get covariance matrices
diff --git a/matlab/nonlinear-filters/online_auxiliary_filter.m b/matlab/nonlinear-filters/online_auxiliary_filter.m
index f81469e6c..c9ab692d7 100644
--- a/matlab/nonlinear-filters/online_auxiliary_filter.m
+++ b/matlab/nonlinear-filters/online_auxiliary_filter.m
@@ -39,7 +39,7 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 % Set seed for randn().
-set_dynare_seed('default');
+DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 pruning = DynareOptions.particle.pruning;
 second_resample = DynareOptions.particle.resampling.status.systematic;
 variance_update = true;
diff --git a/matlab/nonlinear-filters/sequential_importance_particle_filter.m b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
index c69ded02f..323e2cd96 100644
--- a/matlab/nonlinear-filters/sequential_importance_particle_filter.m
+++ b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
@@ -101,7 +101,7 @@ state_variance_rank = size(StateVectorVarianceSquareRoot,2);
 Q_lower_triangular_cholesky = chol(Q)';
 
 % Set seed for randn().
-set_dynare_seed('default');
+DynareOptions=set_dynare_seed_local_options(DynareOptions,'default');
 
 % Initialization of the weights across particles.
 weights = ones(1,number_of_particles)/number_of_particles ;
diff --git a/matlab/parallel/InitializeComputationalEnvironment.m b/matlab/parallel/InitializeComputationalEnvironment.m
index 15372da7c..9d52cf015 100644
--- a/matlab/parallel/InitializeComputationalEnvironment.m
+++ b/matlab/parallel/InitializeComputationalEnvironment.m
@@ -64,7 +64,7 @@ isHybridMatlabOctave = isHybridMatlabOctave && ~isoctave;
 options_.parallel_info.isHybridMatlabOctave = isHybridMatlabOctave;
 if isHybridMatlabOctave
     % Reset dynare random generator and seed.
-    set_dynare_seed('default');
+    options_=set_dynare_seed_local_options(options_,'default');
 end
 
 
diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m
index 9b751c558..405c89c88 100644
--- a/matlab/posterior_sampler_core.m
+++ b/matlab/posterior_sampler_core.m
@@ -126,15 +126,15 @@ for curr_block = fblck:nblck
         %
         % Set the random number generator type (the seed is useless but needed by the function)
         if ~isoctave
-            set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
+            options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
         else
-            set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
+            options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
         end
         % Set the state of the RNG
         set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);
     catch
         % If the state set by master is incompatible with the slave, we only reseed
-        set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
+        options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
     end
     mh_recover_flag=0;
     if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index ab9ffe063..db70a69da 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -176,7 +176,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
     end
     % Find initial values for the NumberOfBlocks chains...
     if NumberOfBlocks > 1 || options_.mh_initialize_from_previous_mcmc.status% Case 1: multiple chains
-        set_dynare_seed('default');
+        options_=set_dynare_seed_local_options(options_,'default');
         fprintf(fidlog,['  Initial values of the parameters:\n']);
         fprintf('%s: Searching for initial values...\n', dispString);
         if ~options_.mh_initialize_from_previous_mcmc.status
@@ -295,7 +295,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
     else
         for j=1:NumberOfBlocks
             % we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
-            set_dynare_seed(options_.DynareRandomStreams.seed+j);
+            options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+j);
             % record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
             [record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state();
         end
diff --git a/matlab/set_dynare_seed.m b/matlab/set_dynare_seed.m
index 9925d12cc..9791bded4 100644
--- a/matlab/set_dynare_seed.m
+++ b/matlab/set_dynare_seed.m
@@ -1,8 +1,10 @@
-function set_dynare_seed(a,b)
-% Set seeds depending on matlab (octave) version. This routine is called in dynare_config and can be called by the
+function set_dynare_seed(varargin)
+% Set seeds depending on matlab (octave) version. This routine is called is
+% a wrapper for set_dynare_seed_local_options
+% dynare_config and can be called by the
 % user in the mod file.
 %
-% Copyright © 2010-2020 Dynare Team
+% Copyright © 2010-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -18,92 +20,11 @@ function set_dynare_seed(a,b)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
 global options_
 
-if ~nargin
+if nargin<1
     error('set_dynare_seed:: I need at least one input argument!')
 end
 
-matlab_random_streams = ~(isoctave || options_.parallel_info.isHybridMatlabOctave);
-
-if matlab_random_streams% Use new matlab interface.
-    if nargin==1
-        if ischar(a) && strcmpi(a,'default')
-            options_.DynareRandomStreams.algo = 'mt19937ar';
-            options_.DynareRandomStreams.seed = 0;
-            s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
-            reset(RandStream.setGlobalStream(s));
-            return
-        end
-        if ischar(a) && strcmpi(a,'reset')
-            s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
-            reset(RandStream.setGlobalStream(s));
-            return
-        end
-        if ~ischar(a) || (ischar(a) && strcmpi(a, 'clock'))
-            options_.DynareRandomStreams.algo = 'mt19937ar';
-            if ischar(a)
-                options_.DynareRandomStreams.seed = rem(floor(now*24*60*60), 2^32);
-            else
-                options_.DynareRandomStreams.seed = a;
-            end
-            s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
-            reset(RandStream.setGlobalStream(s));
-            return
-        end
-        error('set_dynare_seed:: something is wrong in the calling sequence!')
-    elseif nargin==2
-        if ~ischar(a) || ~( strcmpi(a,'mcg16807') || ...
-                            strcmpi(a,'mlfg6331_64') || ...
-                            strcmpi(a,'mrg32k3a') || ...
-                            strcmpi(a,'mt19937ar') || ...
-                            strcmpi(a,'shr3cong') || ...
-                            strcmpi(a,'swb2712') )
-            disp('set_dynare_seed:: First argument must be string designing the uniform random number algorithm!')
-            RandStream.list
-            skipline()
-            disp('set_dynare_seed:: Change the first input accordingly...')
-            skipline()
-            error(' ')
-        end
-        if ~isint(b)
-            error('set_dynare_seed:: The second input argument must be an integer!')
-        end
-        options_.DynareRandomStreams.algo = a;
-        options_.DynareRandomStreams.seed = b;
-        s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
-        reset(RandStream.setGlobalStream(s));
-    end
-else% Use old matlab interface.
-    if nargin==1
-        if ischar(a) && strcmpi(a,'default')
-            if isoctave
-                options_.DynareRandomStreams.algo = 'state';
-            else
-                options_.DynareRandomStreams.algo = 'twister';
-            end
-            options_.DynareRandomStreams.seed = 0;
-            rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
-            randn('state',options_.DynareRandomStreams.seed);
-            return
-        end
-        if ischar(a) && strcmpi(a,'reset')
-            rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
-            randn('state',options_.DynareRandomStreams.seed);
-            return
-        end
-        if (~ischar(a) && isint(a)) || (ischar(a) && strcmpi(a,'clock'))
-            if ischar(a)
-                options_.DynareRandomStreams.seed = floor(now*24*60*60);
-            else
-                options_.DynareRandomStreams.seed = a;
-            end
-            rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
-            randn('state',options_.DynareRandomStreams.seed);
-            return
-        end
-        error('set_dynare_seed:: Something is wrong in the calling sequence!')
-    else
-        error('set_dynare_seed:: Cannot use more than one input argument with your version of Matlab/Octave!')
-    end
-end
+options_=set_dynare_seed_local_options(options_,varargin{:});
\ No newline at end of file
diff --git a/matlab/set_dynare_seed_local_options.m b/matlab/set_dynare_seed_local_options.m
new file mode 100644
index 000000000..266a2a724
--- /dev/null
+++ b/matlab/set_dynare_seed_local_options.m
@@ -0,0 +1,121 @@
+function options_=set_dynare_seed_local_options(options_,a,b)
+% options_=set_dynare_seed_local_options(options_,a,b)
+% Set seeds depending on Matlab (octave) version
+% Inputs:
+%   o options_              options structure
+%   o a                     first input argument, 
+%                           for single argument input, either
+%                               [number]       seed
+%                               [string]        'default'
+%                                               'reset'
+%       	                                    'clock'
+%                           for two inputs: 
+%                               [string]        indicating feasible RNG algorithm (default: 'mt19937ar')
+%   o b                     second input argument
+%                               [number]        seed for algorithm
+%                                               specified with first input
+
+% Copyright © 2010-2023 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 <https://www.gnu.org/licenses/>.
+
+if nargin<2
+    error('set_dynare_seed:: I need at least two input arguments!')
+end
+
+matlab_random_streams = ~(isoctave || options_.parallel_info.isHybridMatlabOctave);
+
+if matlab_random_streams% Use new matlab interface.
+    if nargin==2
+        if ischar(a) && strcmpi(a,'default')
+            options_.DynareRandomStreams.algo = 'mt19937ar';
+            options_.DynareRandomStreams.seed = 0;
+            s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
+            reset(RandStream.setGlobalStream(s));
+            return
+        end
+        if ischar(a) && strcmpi(a,'reset')
+            s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
+            reset(RandStream.setGlobalStream(s));
+            return
+        end
+        if ~ischar(a) || (ischar(a) && strcmpi(a, 'clock'))
+            options_.DynareRandomStreams.algo = 'mt19937ar';
+            if ischar(a)
+                options_.DynareRandomStreams.seed = rem(floor(now*24*60*60), 2^32);
+            else
+                options_.DynareRandomStreams.seed = a;
+            end
+            s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
+            reset(RandStream.setGlobalStream(s));
+            return
+        end
+        error('set_dynare_seed:: something is wrong in the calling sequence!')
+    elseif nargin==3
+        if ~ischar(a) || ~( strcmpi(a,'mcg16807') || ...
+                            strcmpi(a,'mlfg6331_64') || ...
+                            strcmpi(a,'mrg32k3a') || ...
+                            strcmpi(a,'mt19937ar') || ...
+                            strcmpi(a,'shr3cong') || ...
+                            strcmpi(a,'swb2712') )
+            disp('set_dynare_seed:: First argument must be string designing the uniform random number algorithm!')
+            RandStream.list
+            skipline()
+            disp('set_dynare_seed:: Change the first input accordingly...')
+            skipline()
+            error(' ')
+        end
+        if ~isint(b)
+            error('set_dynare_seed:: The second input argument must be an integer!')
+        end
+        options_.DynareRandomStreams.algo = a;
+        options_.DynareRandomStreams.seed = b;
+        s = RandStream(options_.DynareRandomStreams.algo,'Seed',options_.DynareRandomStreams.seed);
+        reset(RandStream.setGlobalStream(s));
+    end
+else% Use old matlab interface.
+    if nargin==2
+        if ischar(a) && strcmpi(a,'default')
+            if isoctave
+                options_.DynareRandomStreams.algo = 'state';
+            else
+                options_.DynareRandomStreams.algo = 'twister';
+            end
+            options_.DynareRandomStreams.seed = 0;
+            rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
+            randn('state',options_.DynareRandomStreams.seed);
+            return
+        end
+        if ischar(a) && strcmpi(a,'reset')
+            rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
+            randn('state',options_.DynareRandomStreams.seed);
+            return
+        end
+        if (~ischar(a) && isint(a)) || (ischar(a) && strcmpi(a,'clock'))
+            if ischar(a)
+                options_.DynareRandomStreams.seed = floor(now*24*60*60);
+            else
+                options_.DynareRandomStreams.seed = a;
+            end
+            rand(options_.DynareRandomStreams.algo,options_.DynareRandomStreams.seed);
+            randn('state',options_.DynareRandomStreams.seed);
+            return
+        end
+        error('set_dynare_seed:: Something is wrong in the calling sequence!')
+    else
+        error('set_dynare_seed:: Cannot use more than one input argument with your version of Matlab/Octave!')
+    end
+end
diff --git a/matlab/simul_static_model.m b/matlab/simul_static_model.m
index 51d417a7c..bf3978a81 100644
--- a/matlab/simul_static_model.m
+++ b/matlab/simul_static_model.m
@@ -53,7 +53,7 @@ if nargin<2 || isempty(innovations)
     covariance_matrix_upper_cholesky = chol(covariance_matrix);
     % Set seed to its default state.
     if options_.bnlms.set_dynare_seed_to_default
-        set_dynare_seed('default');
+        options_=set_dynare_seed_local_options(options_,'default');
     end
     % Simulate structural innovations.
     switch options_.bnlms.innovation_distribution