diff --git a/matlab/CutSample.m b/matlab/CutSample.m
index 409f65cbfa552d6e8c2bb00a4e36269601a4dc3f..e29f010d5b2156eb74c8074e08700c866b3c8f32 100644
--- a/matlab/CutSample.m
+++ b/matlab/CutSample.m
@@ -16,7 +16,7 @@ function CutSample(M_, options_, estim_params_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2005-2017 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,8 +33,6 @@ function CutSample(M_, options_, estim_params_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
-
 % Get the path to the metropolis files.
 MetropolisFolder = CheckPath('metropolis',M_.dname);
 
@@ -42,7 +40,8 @@ MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
 
 % Load the last mh-history file.
-load_last_mh_history_file(MetropolisFolder, M_.fname);
+record=load_last_mh_history_file(MetropolisFolder, M_.fname);
+npar=size(record.LastParameters,2);
 
 % Get the list of files where the mcmc draw are saved.
 mh_files = dir([ MetropolisFolder ,filesep, M_.fname '_mh*.mat' ]);
diff --git a/matlab/GetAllPosteriorDraws.m b/matlab/GetAllPosteriorDraws.m
index bc6a21c37e1b3af7ce0e867a8dd2e5b28c587d8e..f8cde8ef29d406f3a3f87f4dcbaaddf560554474 100644
--- a/matlab/GetAllPosteriorDraws.m
+++ b/matlab/GetAllPosteriorDraws.m
@@ -1,22 +1,23 @@
-function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, blck)
-
-% function Draws = GetAllPosteriorDraws(column,FirstMhFile,FirstLine,TotalNumberOfMhFile,NumberOfDraws)
+function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
+% function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
 % Gets all posterior draws
 %
 % INPUTS
-%    column:               column
+%    column:               column of desired parameter in draw matrix
 %    FirstMhFile:          first mh file
 %    FirstLine:            first line
 %    TotalNumberOfMhFile:  total number of mh file
 %    NumberOfDraws:        number of draws
-
+%    nblcks:               total number of blocks 
+%    blck:                 desired block to read
+%
 % OUTPUTS
 %    Draws:                draws from posterior distribution
 %
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2005-2017 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,19 +34,17 @@ function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumbe
 % 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_
-
-nblck = options_.mh_nblck;
+global M_ 
 
 iline = FirstLine;
 linee = 1;
 DirectoryName = CheckPath('metropolis',M_.dname);
 
-if nblck>1 && nargin<6
-    Draws = zeros(NumberOfDraws*nblck,1);
+if nblcks>1 && nargin<7
+    Draws = zeros(NumberOfDraws*nblcks,1);
     iline0=iline;
     if column>0
-        for blck = 1:nblck
+        for blck = 1:nblcks
             iline=iline0;
             for file = FirstMhFile:TotalNumberOfMhFile
                 load([DirectoryName '/'  M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
@@ -56,7 +55,7 @@ if nblck>1 && nargin<6
             end
         end
     else
-        for blck = 1:nblck
+        for blck = 1:nblcks
             iline=iline0;
             for file = FirstMhFile:TotalNumberOfMhFile
                 load([DirectoryName '/'  M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2')
@@ -68,7 +67,7 @@ if nblck>1 && nargin<6
         end
     end
 else
-    if nblck==1
+    if nblcks==1
         blck=1;
     end
     if column>0
diff --git a/matlab/GetPosteriorMeanVariance.m b/matlab/GetPosteriorMeanVariance.m
index fa095c6aa032223f0aee3ac40e9c27f80ab3e496..562f39486b1430bb3052a88d1eb4fd3736c0c176 100644
--- a/matlab/GetPosteriorMeanVariance.m
+++ b/matlab/GetPosteriorMeanVariance.m
@@ -20,7 +20,7 @@ function [mean,variance] = GetPosteriorMeanVariance(M,drop)
 MetropolisFolder = CheckPath('metropolis',M.dname);
 FileName = M.fname;
 BaseName = [MetropolisFolder filesep FileName];
-load_last_mh_history_file(MetropolisFolder, FileName);
+record=load_last_mh_history_file(MetropolisFolder, FileName);
 NbrDraws = sum(record.MhDraws(:,1));
 NbrFiles = sum(record.MhDraws(:,2));
 NbrBlocks = record.Nblck;
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index e33fae487493f8bc52e3d1b5afa5530f64078ecc..dee97ac76af8e69d526e97e786b279740e0a58e8 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -16,7 +16,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright © 2006-2018 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -38,19 +38,17 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
 %end
 
 TeX     = options_.TeX;
-nblck   = options_.mh_nblck;
 nvx     = estim_params_.nvx;
 nvn     = estim_params_.nvn;
 ncx     = estim_params_.ncx;
 ncn     = estim_params_.ncn;
 np      = estim_params_.np ;
-nx      = nvx+nvn+ncx+ncn+np;
 
 MetropolisFolder = CheckPath('metropolis',M_.dname);
 OutputFolder = CheckPath('Output',M_.dname);
 FileName = M_.fname;
 
-load_last_mh_history_file(MetropolisFolder,FileName);
+record=load_last_mh_history_file(MetropolisFolder,FileName);
 
 FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine = record.KeepedDraws.FirstLine;
@@ -58,6 +56,7 @@ TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
 FirstMhFile = record.KeepedDraws.FirstMhFile;
 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+mh_nblck = size(record.LastParameters,1);
 clear record;
 
 header_width = row_header_width(M_, estim_params_, bayestopt_);
@@ -76,7 +75,7 @@ catch
     disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean))
 end
 
-num_draws=NumberOfDraws*options_.mh_nblck;
+num_draws=NumberOfDraws*mh_nblck;
 hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
 
 if hpd_draws<2
@@ -97,7 +96,7 @@ if np
     ip = nvx+nvn+ncx+ncn+1;
     for i=1:np
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
-            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
             name = bayestopt_.name{ip};
             oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@@ -106,7 +105,7 @@ if np
                 name = bayestopt_.name{ip};
                 [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
             catch
-                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
                 name = bayestopt_.name{ip};
                 oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@@ -140,7 +139,7 @@ if nvx
     disp(tit2)
     for i=1:nvx
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
-            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
             k = estim_params_.var_exo(i,1);
             name = M_.exo_names{k};
@@ -152,7 +151,7 @@ if nvx
                 name = M_.exo_names{k};
                 [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
             catch
-                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
                     posterior_moments(Draws, 1, options_.mh_conf_sig);
                 k = estim_params_.var_exo(i,1);
@@ -184,7 +183,7 @@ if nvn
     ip = nvx+1;
     for i=1:nvn
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
-            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
             name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
             oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@@ -193,7 +192,7 @@ if nvn
                 name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
             catch
-                Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+                Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws, mh_nblck);
                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1,options_.mh_conf_sig);
                 name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
                 oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
@@ -223,7 +222,7 @@ if ncx
     ip = nvx+nvn+1;
     for i=1:ncx
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
-            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+            Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1,options_.mh_conf_sig);
             k1 = estim_params_.corrx(i,1);
             k2 = estim_params_.corrx(i,2);
@@ -240,7 +239,7 @@ if ncx
                 NAME = sprintf('%s_%s', M_.exo_names{k1}, M_.exo_names{k2});
                 [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
             catch
-                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
                 k1 = estim_params_.corrx(i,1);
                 k2 = estim_params_.corrx(i,2);
@@ -273,7 +272,7 @@ if ncn
     ip = nvx+nvn+ncx+1;
     for i=1:ncn
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
-            Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+            Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws, mh_nblck);
             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
             k1 = estim_params_.corrn(i,1);
             k2 = estim_params_.corrn(i,2);
@@ -288,7 +287,7 @@ if ncn
                 NAME = sprintf('%s_%s', M_.endo_names{k1}, M_.endo_names{k2});
                 [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
             catch
-                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+                Draws = GetAllPosteriorDraws(ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
                 k1 = estim_params_.corrn(i,1);
                 k2 = estim_params_.corrn(i,2);
diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m
index b4435190a596695e49d606ae252a56cda7ed5cb1..c8f96e270c08a4748352721c1ac0f240b0cf7e43 100644
--- a/matlab/compute_mh_covariance_matrix.m
+++ b/matlab/compute_mh_covariance_matrix.m
@@ -15,7 +15,7 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -32,27 +32,20 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at
 % 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_ bayestopt_
-
-
-n = estim_params_.np + ...
-    estim_params_.nvn+ ...
-    estim_params_.ncx+ ...
-    estim_params_.ncn+ ...
-    estim_params_.nvx;
-
-nblck = options_.mh_nblck;
+global M_ bayestopt_
 
 MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
 BaseName = [MetropolisFolder filesep ModelName];
 
-load_last_mh_history_file(MetropolisFolder, ModelName);
+record=load_last_mh_history_file(MetropolisFolder, ModelName);
 
 FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine   = record.KeepedDraws.FirstLine;
 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 
+[nblck, n] = size(record.LastParameters);
+
 posterior_kernel_at_the_mode = -Inf;
 posterior_mean = zeros(n,1);
 posterior_mode = NaN(n,1);
diff --git a/matlab/convergence_diagnostics/McMCDiagnostics.m b/matlab/convergence_diagnostics/McMCDiagnostics.m
index 16f58ed852a83f1048c91b40c9d48cea36b66fd5..1d0da5d61364b76a16508e9f0fe98711b82dbf15 100644
--- a/matlab/convergence_diagnostics/McMCDiagnostics.m
+++ b/matlab/convergence_diagnostics/McMCDiagnostics.m
@@ -16,7 +16,7 @@ function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
 % PARALLEL CONTEXT
 % See the comment in posterior_sampler.m funtion.
 
-% Copyright © 2005-2018 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -38,25 +38,12 @@ MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
 
 TeX = options_.TeX;
-nblck = options_.mh_nblck;
-npar = estim_params_.nvx;
-npar = npar + estim_params_.nvn;
-npar = npar + estim_params_.ncx;
-npar = npar + estim_params_.ncn;
-npar = npar + estim_params_.np ;
-MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
-
-load_last_mh_history_file(MetropolisFolder, ModelName);
+
+record=load_last_mh_history_file(MetropolisFolder, ModelName);
 NumberOfMcFilesPerBlock = record.LastFileNumber;
+[nblck, npar] = size(record.LastParameters);
 npardisp = options_.convergence.brooksgelman.plotrows;
 
-% Check that the declared number of blocks is consistent with informations saved in mh-history files.
-if ~isequal(nblck,record.Nblck)
-    disp(['Estimation::mcmc::diagnostics: The number of MCMC chains you declared (' num2str(nblck) ') is inconsistent with the information available in the mh-history files (' num2str(record.Nblck) ' chains)!'])
-    disp(['                               I reset the number of MCMC chains to ' num2str(record.Nblck) '.'])
-    nblck = record.Nblck;
-end
-
 % check if all the mh files are available.
 issue_an_error_message = 0;
 for b = 1:nblck
@@ -91,7 +78,7 @@ for jj = 1:npar
         par_name_temp = get_the_name(jj, options_.TeX, M_, estim_params_, options_);
         param_name = vertcat(param_name, par_name_temp);
     end
-    Draws = GetAllPosteriorDraws(jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+    Draws = GetAllPosteriorDraws(jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
     Draws = reshape(Draws, [NumberOfDraws nblck]);
     Nc = min(1000, NumberOfDraws/2);
     for ll = 1:nblck
diff --git a/matlab/generate_trace_plots.m b/matlab/generate_trace_plots.m
index 2fcbda4f3f2297cc8a363e0a269aed6e75349d11..a241cc07617b37a68e9d0f9e8047472c2c00cde0 100644
--- a/matlab/generate_trace_plots.m
+++ b/matlab/generate_trace_plots.m
@@ -11,7 +11,7 @@ function generate_trace_plots(chain_number)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2016-2018 Dynare Team
+% Copyright © 2016-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,7 +33,7 @@ global M_ options_ estim_params_
 
 % Get informations about the posterior draws:
 MetropolisFolder = CheckPath('metropolis', M_.dname);
-load_last_mh_history_file(MetropolisFolder, M_.fname);
+record=load_last_mh_history_file(MetropolisFolder, M_.fname);
 if chain_number>record.Nblck
     error('generate_trace_plots:: chain number is bigger than existing number of chains')
 end
diff --git a/matlab/get_name_of_the_last_mh_file.m b/matlab/get_name_of_the_last_mh_file.m
index 5d49a780ba3afc8646f4ece3b6c4172886ab1cb7..6c9ccca24015839b763a32237ae691145c3b4dd8 100644
--- a/matlab/get_name_of_the_last_mh_file.m
+++ b/matlab/get_name_of_the_last_mh_file.m
@@ -11,7 +11,7 @@ function [mhname,info] = get_name_of_the_last_mh_file(M_)
 %                          file. Otherwise info is equal to zero (a likely reason for this is
 %                          that the mcmc simulations were not completed).
 
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -35,7 +35,7 @@ MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
 BaseName = [MetropolisFolder filesep ModelName];
 
-load_last_mh_history_file(MetropolisFolder, ModelName);
+record=load_last_mh_history_file(MetropolisFolder, ModelName);
 
 mh_number = record.LastFileNumber ;
 bk_number = record.Nblck ;
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index 8bca13e6681c7204709d7dcc4d12e4983356b9f8..8eb56bed4a05c85400f3cfe41ed8fb2bb64aa36f 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -15,7 +15,7 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2005-2018 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,14 +33,13 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 
-npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
-nblck = options_.mh_nblck;
-
 MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
 BaseName = [MetropolisFolder filesep ModelName];
 
-load_last_mh_history_file(MetropolisFolder, ModelName);
+record=load_last_mh_history_file(MetropolisFolder, ModelName);
+
+[nblck, npar] = size(record.LastParameters);
 
 FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m
index 533d79ad36793d13e298435bb0a83ec310fbbb84..b5dbe38e401dde5716e1d64f7ce53142cdab1028 100644
--- a/matlab/metropolis_draw.m
+++ b/matlab/metropolis_draw.m
@@ -24,7 +24,7 @@ function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params
 %
 %   Requires CutSample to be run before in order to set up mh_history-file
 
-% Copyright © 2003-2017 Dynare Team
+% Copyright © 2003-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -59,7 +59,7 @@ if init
     FileName = M_.fname;
     BaseName = [MetropolisFolder filesep FileName];
     %load mh_history-file with info on what to load
-    load_last_mh_history_file(MetropolisFolder, FileName);
+    record=load_last_mh_history_file(MetropolisFolder, FileName);
     FirstMhFile = record.KeepedDraws.FirstMhFile;
     FirstLine = record.KeepedDraws.FirstLine;
     TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
diff --git a/matlab/mh_autocorrelation_function.m b/matlab/mh_autocorrelation_function.m
index 16524a190b4c8b96e15a5971c30898277bc770aa..8c07143bd1cd5e3289fedf18457d0bc737552d27 100644
--- a/matlab/mh_autocorrelation_function.m
+++ b/matlab/mh_autocorrelation_function.m
@@ -18,7 +18,7 @@ function mh_autocorrelation_function(options_,M_,estim_params_,type,blck,name1,n
 %
 % SPECIAL REQUIREMENTS
 
-% Copyright © 2003-2017 Dynare Team
+% Copyright © 2003-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -48,17 +48,19 @@ end
 
 % Get informations about the posterior draws:
 MetropolisFolder = CheckPath('metropolis',M_.dname);
-load_last_mh_history_file(MetropolisFolder, M_.fname);
+record=load_last_mh_history_file(MetropolisFolder, M_.fname);
 
 FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+nblck = size(record.LastParameters,1);
+
 clear record;
 
 % Get all the posterior draws:
-PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, blck);
+PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck, blck);
 
 % Compute the autocorrelation function:
 [autocov,autocor] = sample_autocovariance(PosteriorDraws,options_.mh_autocorrelation_function_size);
diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m
index 2942a61295743461991dfdcfafdde603a1d86877..d16ff754016daab9ec06453bd895d23396390a46 100644
--- a/matlab/posterior_sampler.m
+++ b/matlab/posterior_sampler.m
@@ -37,7 +37,7 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun
 % Then the comments write here can be used for all the other pairs of
 % parallel functions and also for management functions.
 
-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -62,7 +62,7 @@ vv = sampler_options.invhess;
 InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
 
 % Load last mh history file
-load_last_mh_history_file(MetropolisFolder, ModelName);
+record=load_last_mh_history_file(MetropolisFolder, ModelName);
 
 % Only for test parallel results!!!
 
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index 60a909d2c6d6993707dcf33bd667aa21c01a2f26..9a819026906a0bb5ad419d73f17c6d1ec6a83a9d 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -340,7 +340,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
 elseif options_.load_mh_file && ~options_.mh_recover
     % Here we consider previous mh files (previous mh did not crash).
     disp('Estimation::mcmc: I am loading past Metropolis-Hastings simulations...')
-    load_last_mh_history_file(MetropolisFolder, ModelName);
+    record=load_last_mh_history_file(MetropolisFolder, ModelName);
     if ~isnan(record.MCMCConcludedSuccessfully) && ~record.MCMCConcludedSuccessfully
         error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.')
     end
@@ -400,7 +400,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
 elseif options_.mh_recover
     % The previous metropolis-hastings crashed before the end! I try to recover the saved draws...
     disp('Estimation::mcmc: Recover mode!')
-    load_last_mh_history_file(MetropolisFolder, ModelName);
+    record=load_last_mh_history_file(MetropolisFolder, ModelName);
     NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains.
     options_.mh_nblck = NumberOfBlocks;
 
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index e93387afb6e46c78e52f70cecf7bfbbe8c237702..5f2b61db5116b2391522f43c602f76234876549f 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -19,7 +19,7 @@ function prior_posterior_statistics(type,dataset,dataset_info)
 % See the comments in the posterior_sampler.m funtion.
 
 
-% Copyright © 2005-2020 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -209,7 +209,8 @@ localVars.DirectoryName = DirectoryName;
 
 if strcmpi(type,'posterior')
     BaseName = [DirectoryName filesep M_.fname];
-    load_last_mh_history_file(DirectoryName, M_.fname);
+    record=load_last_mh_history_file(DirectoryName, M_.fname);
+    [nblck, npar] = size(record.LastParameters);
     FirstMhFile = record.KeepedDraws.FirstMhFile;
     FirstLine = record.KeepedDraws.FirstLine;
     TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
@@ -219,9 +220,9 @@ if strcmpi(type,'posterior')
     mh_nblck = options_.mh_nblck;
     if B==NumberOfDraws*mh_nblck
         % we load all retained MH runs !
-        logpost=GetAllPosteriorDraws(0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+        logpost=GetAllPosteriorDraws(0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
         for column=1:npar
-            x(:,column) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws);
+            x(:,column) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
         end
     else
         logpost=NaN(B,1);
diff --git a/matlab/selec_posterior_draws.m b/matlab/selec_posterior_draws.m
index 0a10dab82758751c0511fb06500e22b06edab2a7..0b5dadcc5ba1b38134b581167ccb211fb7b6dfe3 100644
--- a/matlab/selec_posterior_draws.m
+++ b/matlab/selec_posterior_draws.m
@@ -67,7 +67,7 @@ ModelName = M_.fname;
 BaseName = [MetropolisFolder filesep ModelName];
 
 % Get informations about the mcmc:
-load_last_mh_history_file(MetropolisFolder, ModelName);
+record=load_last_mh_history_file(MetropolisFolder, ModelName);
 FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine = record.KeepedDraws.FirstLine;
 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
diff --git a/matlab/trace_plot.m b/matlab/trace_plot.m
index af4a6b876bcfb75a845700e714f90b06a3154160..8deecb8174c4372d34ec499d41187671363a85fb 100644
--- a/matlab/trace_plot.m
+++ b/matlab/trace_plot.m
@@ -53,16 +53,17 @@ end
 
 % Get informations about the posterior draws:
 MetropolisFolder = CheckPath('metropolis',M_.dname);
-load_last_mh_history_file(MetropolisFolder, M_.fname);
+record=load_last_mh_history_file(MetropolisFolder, M_.fname);
 
 FirstMhFile = 1;
 FirstLine = 1;
 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+[mh_nblck] = size(record.LastParameters,2);
 clear record;
 
 % Get all the posterior draws:
-PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, blck);
+PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck);
 
 
 % Plot the posterior draws: