diff --git a/matlab/estimation/posterior_sampler_core.m b/matlab/estimation/posterior_sampler_core.m
index bf269a1410a075023d49224bc3b5151bed71d948..61492fd4fe09369cfa62f66a698076cf9bc1f2b9 100644
--- a/matlab/estimation/posterior_sampler_core.m
+++ b/matlab/estimation/posterior_sampler_core.m
@@ -86,12 +86,6 @@ endo_steady_state = myinputs.endo_steady_state;
 exo_steady_state=myinputs.exo_steady_state;
 exo_det_steady_state=myinputs.exo_det_steady_state;
 
-% Necessary only for remote computing!
-if whoiam
-    % initialize persistent variables in priordens()
-    priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
-end
-
 MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
 BaseName = [MetropolisFolder filesep ModelName];
diff --git a/matlab/estimation/priordens.m b/matlab/estimation/priordens.m
index f285606510ed9d0ffccc5ca8326360feac98a134..e620d232efeaed629f59f7f790424a69d6b8d35b 100644
--- a/matlab/estimation/priordens.m
+++ b/matlab/estimation/priordens.m
@@ -1,4 +1,4 @@
-function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4, initialization)
+function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4)
 % Computes a prior density for the structural parameters of DSGE models
 %
 % INPUTS
@@ -8,14 +8,13 @@ function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape,
 %    p7:            [double]      vector with n elements, second parameter of the prior distribution (bayestopt_.p7).
 %    p3:            [double]      vector with n elements, lower bounds of the untruncated standard or generalized distribution
 %    p4:            [double]      vector with n elements, upper bound of the untruncated standard or generalized distribution
-%    initialization [integer]     if 1: initialize persistent variables
 %
 % OUTPUTS
 %    logged_prior_density  [double]  scalar, log of the prior density evaluated at x.
 %    info                  [double]  error code for index of Inf-prior parameter
 %
 
-% Copyright © 2003-2023 Dynare Team
+% Copyright © 2003-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -32,54 +31,49 @@ function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape,
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-persistent id1 id2 id3 id4 id5 id6 id8
-persistent tt1 tt2 tt3 tt4 tt5 tt6 tt8
-
 info=0;
 
-if nargin > 6  && initialization
-    % Beta indices.
-    tt1 = true;
-    id1 = find(pshape==1);
-    if isempty(id1)
-        tt1 = false;
-    end
-    % Gamma indices.
-    tt2 = true;
-    id2 = find(pshape==2);
-    if isempty(id2)
-        tt2 = false;
-    end
-    % Gaussian indices.
-    tt3 = true;
-    id3 = find(pshape==3);
-    if isempty(id3)
-        tt3 = false;
-    end
-    % Inverse-Gamma-1 indices.
-    tt4 = true;
-    id4 = find(pshape==4);
-    if isempty(id4)
-        tt4 = false;
-    end
-    % Uniform indices.
-    tt5 = true;
-    id5 = find(pshape==5);
-    if isempty(id5)
-        tt5 = false;
-    end
-    % Inverse-Gamma-2 indices.
-    tt6 = true;
-    id6 = find(pshape==6);
-    if isempty(id6)
-        tt6 = false;
-    end
-    % Weibull indices.
-    tt8 = true;
-    id8 = find(pshape==8);
-    if isempty(id8)
-        tt8 = false;
-    end
+% Beta indices.
+tt1 = true;
+id1 = find(pshape==1);
+if isempty(id1)
+    tt1 = false;
+end
+% Gamma indices.
+tt2 = true;
+id2 = find(pshape==2);
+if isempty(id2)
+    tt2 = false;
+end
+% Gaussian indices.
+tt3 = true;
+id3 = find(pshape==3);
+if isempty(id3)
+    tt3 = false;
+end
+% Inverse-Gamma-1 indices.
+tt4 = true;
+id4 = find(pshape==4);
+if isempty(id4)
+    tt4 = false;
+end
+% Uniform indices.
+tt5 = true;
+id5 = find(pshape==5);
+if isempty(id5)
+    tt5 = false;
+end
+% Inverse-Gamma-2 indices.
+tt6 = true;
+id6 = find(pshape==6);
+if isempty(id6)
+    tt6 = false;
+end
+% Weibull indices.
+tt8 = true;
+id8 = find(pshape==8);
+if isempty(id8)
+    tt8 = false;
 end
 
 logged_prior_density = 0.0;
@@ -257,8 +251,8 @@ end
 
 % Call the tested routine
 try
-    % Initialization of priordens.
-    lpdstar = priordens(p5, p0, p6, p7, p3, p4, true);
+    % Prior (logged) density at the mode
+    lpdstar = priordens(p5, p0, p6, p7, p3, p4);
     % Do simulations in a loop and estimate recursively the mean and teh variance.
     LPD = NaN(10000,1);
     for i = 1:10000
diff --git a/matlab/estimation/set_prior.m b/matlab/estimation/set_prior.m
index e6739f4686f7433ef53e285cb02df6890626c5b1..3069a8a09598ace5870d106aee51f254b80af1d9 100644
--- a/matlab/estimation/set_prior.m
+++ b/matlab/estimation/set_prior.m
@@ -314,7 +314,3 @@ if exist([ M_.dname '/prior/definition.mat'],'file')
 else
     save([M_.dname '/prior/definition.mat'],'bayestopt_');
 end
-
-% initialize persistent variables in priordens()
-priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
-          bayestopt_.p3,bayestopt_.p4,1);
\ No newline at end of file
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 3853e28d30b422315083329c6bd7dc424db8334b..00da3b260b06a791d3c616956cd674208a712688 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -136,8 +136,6 @@ else
     options_ = default_option_values(M_);
 end
 
-% initialize persistent variables in priordens()
-priordens([],[],[],[],[],[],1);
 % initialize persistent variables in dyn_first_order_solver()
 dyn_first_order_solver();
 
diff --git a/matlab/list_of_functions_to_be_cleared.m b/matlab/list_of_functions_to_be_cleared.m
index 200478bb2b0060eff7f8945ff55f8cc79472d13f..90c7187b3c8c34dae908937f592ff520c492f735 100644
--- a/matlab/list_of_functions_to_be_cleared.m
+++ b/matlab/list_of_functions_to_be_cleared.m
@@ -1,2 +1 @@
-list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', '+gsa/prior_draw.m', '+identification/analysis.m', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', 'priordens',...
-    '+occbin/solver.m','+occbin/mkdatap_anticipated_dyn.m','+occbin/mkdatap_anticipated_2constraints_dyn.m','+occbin/match_function.m','+occbin/solve_one_constraint.m','+occbin/solve_two_constraint.m','+occbin/plot/shock_decomposition.m'};
+list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', '+gsa/prior_draw.m', '+identification/analysis.m', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', '+occbin/solver.m','+occbin/mkdatap_anticipated_dyn.m','+occbin/mkdatap_anticipated_2constraints_dyn.m','+occbin/match_function.m','+occbin/solve_one_constraint.m','+occbin/solve_two_constraint.m','+occbin/plot/shock_decomposition.m'};