diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m
index 8b9e8aae13eab750746be7a8d0162799fbd134df..e50f1adae8845c2e2dccdd31e15757ea5a6e9227 100644
--- a/matlab/compute_moments_varendo.m
+++ b/matlab/compute_moments_varendo.m
@@ -1,4 +1,4 @@
-function oo_ = compute_moments_varendo(type, options_, M_, oo_, var_list_)
+function oo_ = compute_moments_varendo(type, options_, M_, oo_, estim_params_, var_list_)
 
 % Computes the second order moments (autocorrelation function, covariance
 % matrix and variance decomposition) distributions for all the endogenous variables selected in
@@ -55,7 +55,7 @@ end
 
 if strcmpi(type,'posterior')
     posterior = 1;
-    if nargin==4
+    if nargin==5
         var_list_ = options_.varobs;
     end
     if isfield(oo_,'PosteriorTheoreticalMoments')
@@ -63,7 +63,7 @@ if strcmpi(type,'posterior')
     end
 elseif strcmpi(type,'prior')
     posterior = 0;
-    if nargin==4
+    if nargin==5
         var_list_ = options_.prior_analysis_endo_var_list;
         if isempty(var_list_)
             options_.prior_analysis_var_list = options_.varobs;
@@ -97,13 +97,13 @@ end
 if posterior
     for i=1:NumberOfEndogenousVariables
         for j=i:NumberOfEndogenousVariables
-            oo_ = posterior_analysis('variance', var_list_{i}, var_list_{j}, NumberOfLags, options_, M_, oo_);
+            oo_ = posterior_analysis('variance', var_list_{i}, var_list_{j}, NumberOfLags, options_, M_, oo_, estim_params_);
         end
     end
 else
     for i=1:NumberOfEndogenousVariables
         for j=i:NumberOfEndogenousVariables
-            oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, [], options_, M_, oo_);
+            oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, [], options_, M_, oo_, estim_params_);
         end
     end
 end
@@ -113,7 +113,7 @@ if posterior
     for h=NumberOfLags:-1:1
         for i=1:NumberOfEndogenousVariables
             for j=1:NumberOfEndogenousVariables
-                oo_ = posterior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_);
+                oo_ = posterior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_, estim_params_);
             end
         end
     end
@@ -121,7 +121,7 @@ else
     for h=NumberOfLags:-1:1
         for i=1:NumberOfEndogenousVariables
             for j=1:NumberOfEndogenousVariables
-                oo_ = prior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_);
+                oo_ = prior_analysis('correlation', var_list_{i}, var_list_{j}, h, options_, M_, oo_, estim_params_);
             end
         end
     end
@@ -135,7 +135,7 @@ if options_.order==1
             if posterior
                 for i=1:NumberOfEndogenousVariables
                     for j=1:NumberOfExogenousVariables
-                        oo_ = posterior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_);
+                        oo_ = posterior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_, estim_params_);
                         temp(i,j) = oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.Mean.(var_list_{i}).(M_.exo_names{j});
                     end
                 end
@@ -144,7 +144,7 @@ if options_.order==1
             else
                 for i=1:NumberOfEndogenousVariables
                     for j=1:NumberOfExogenousVariables
-                        oo_ = prior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_);
+                        oo_ = prior_analysis('decomposition', var_list_{i}, M_.exo_names{j}, [], options_, M_, oo_, estim_params_);
                         temp(i,j)=oo_.PriorTheoreticalMoments.dsge.VarianceDecomposition.Mean.(var_list_{i}).(M_.exo_names{j});
                     end
                 end
@@ -181,7 +181,7 @@ if options_.order==1
                             temp(i,j,:) = oo_.PosteriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(observable_name_requested_vars{i}).(M_.exo_names{j});
                         end
                         endo_index_varlist = strmatch(observable_name_requested_vars{i}, var_list_, 'exact');
-                        oo_ = posterior_analysis('decomposition', var_list_{endo_index_varlist}, 'ME', [], options_, M_, oo_);
+                        oo_ = posterior_analysis('decomposition', var_list_{endo_index_varlist}, 'ME', [], options_, M_, oo_, estim_params_);
                         temp(i,j+1,:) = oo_.PosteriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(observable_name_requested_vars{i}).('ME');
                     end
                     title='Posterior mean variance decomposition (in percent) with measurement error';
@@ -219,7 +219,7 @@ if options_.order==1
             if posterior
                 for i=1:NumberOfEndogenousVariables
                     for j=1:NumberOfExogenousVariables
-                        oo_ = posterior_analysis('conditional decomposition', var_list_{i}, M_.exo_names{j}, Steps, options_, M_, oo_);
+                        oo_ = posterior_analysis('conditional decomposition', var_list_{i}, M_.exo_names{j}, Steps, options_, M_, oo_, estim_params_);
                         temp(i,j,:) = oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean.(var_list_{i}).(M_.exo_names{j});
                     end
                 end
@@ -261,7 +261,7 @@ if options_.order==1
                                 temp(i,j,:) = oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(observable_name_requested_vars{i}).(M_.exo_names{j});
                             end
                             endo_index_varlist = strmatch(observable_name_requested_vars{i}, var_list_, 'exact');
-                            oo_ = posterior_analysis('conditional decomposition', var_list_{endo_index_varlist}, 'ME', Steps, options_, M_, oo_);
+                            oo_ = posterior_analysis('conditional decomposition', var_list_{endo_index_varlist}, 'ME', Steps, options_, M_, oo_, estim_params_);
                             temp(i,j+1,:) = oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(observable_name_requested_vars{i}).('ME');
                         end
                         title = 'Posterior mean conditional variance decomposition (in percent) with measurement error';
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 69a5e7fd538ddeefc73ba979743365289f910621..1e216266b40382372c5a57a0af884683b46538b0 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -497,7 +497,7 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
                        end
                    end
                 end
-                oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
+                oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
             end
             if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
                 if error_flag
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
index 2c8afc8468afbc0a5b33b2359d791af265e722d1..3b523c0a7c392bda8dbfacea7d3cb63c13bd9ff9 100644
--- a/matlab/posterior_analysis.m
+++ b/matlab/posterior_analysis.m
@@ -1,4 +1,4 @@
-function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_)
+function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_,estim_params_)
 
 % Copyright © 2008-2021 Dynare Team
 %
@@ -29,7 +29,7 @@ switch info
     if drsize*SampleSize>MaxMegaBytes
         drsize=0;
     end
-    SampleAddress = selec_posterior_draws(SampleSize,drsize);
+    selec_posterior_draws(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state,estim_params_,SampleSize,drsize); %save draws to disk
     oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
   case {4,5}
     oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
diff --git a/matlab/selec_posterior_draws.m b/matlab/selec_posterior_draws.m
index 8aa1bf1452afb78012e23dbcf572c7a804524e45..72b094ab7cf2d4732b1e653b6464e67dd67092b8 100644
--- a/matlab/selec_posterior_draws.m
+++ b/matlab/selec_posterior_draws.m
@@ -1,4 +1,4 @@
-function SampleAddress = selec_posterior_draws(SampleSize,drsize)
+function SampleAddress = selec_posterior_draws(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,estim_params_,SampleSize,drsize)
 % Selects a sample of draws from the posterior distribution and if nargin>1
 % saves the draws in _pdraws mat files (metropolis folder). If drsize>0
 % the dr structure, associated to the parameters, is also saved in _pdraws.
@@ -6,8 +6,14 @@ function SampleAddress = selec_posterior_draws(SampleSize,drsize)
 % _mh file cannot be opened twice.
 %
 % INPUTS
-%   o SampleSize     [integer]  Size of the sample to build.
-%   o drsize         [double]   structure dr is drsize megaoctets.
+%   o M_                    [structure]     Matlab's structure describing the model
+%   o options_              [structure]     Matlab's structure describing the current options
+%   o dr                    [structure]     Reduced form model.
+%   o endo_steady_state     [vector]        steady state value for endogenous variables
+%   o exo_steady_state      [vector]        steady state value for exogenous variables
+%   o exo_det_steady_state  [vector]        steady state value for exogenous deterministic variables                                    
+%   o SampleSize            [integer]       Size of the sample to build.
+%   o drsize                [double]        structure dr is drsize megaoctets.
 %
 % OUTPUTS
 %   o SampleAddress  [integer]  A (SampleSize*4) array, each line specifies the
@@ -37,8 +43,6 @@ function SampleAddress = selec_posterior_draws(SampleSize,drsize)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global M_ options_ estim_params_ oo_
-
 % Number of parameters:
 npar = estim_params_.nvx;
 npar = npar + estim_params_.nvn;
@@ -48,9 +52,9 @@ npar = npar + estim_params_.np;
 
 % Select one task:
 switch nargin
-  case 1
+  case 8
     info = 0;
-  case 2
+  case 9
     MAX_mega_bytes = 10;% Should be an option...
     if drsize>0
         info=2;
@@ -114,8 +118,8 @@ if info
             end
             pdraws(i,1) = {x2(SampleAddress(i,4),:)};
             if info==2
-                set_parameters(pdraws{i,1});
-                [dr,~,M_.params] =compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+                M_ = set_parameters_locally(M_,pdraws{i,1});
+                [dr,~,M_.params] =compute_decision_rules(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
                 pdraws(i,2) = { dr };
             end
             old_mhfile = mhfile;
@@ -141,7 +145,7 @@ if info
             end
             pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
             if info==2
-                set_parameters(pdraws{linee,1});
+                M_ = set_parameters_locally(M_,pdraws{i,1});
                 [dr,~,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
                 pdraws(linee,2) = { dr };
             end