diff --git a/matlab/GetPosteriorMeanVariance.m b/matlab/GetPosteriorMeanVariance.m
index 562f39486..5a2c44af2 100644
--- a/matlab/GetPosteriorMeanVariance.m
+++ b/matlab/GetPosteriorMeanVariance.m
@@ -1,6 +1,20 @@
-function [mean,variance] = GetPosteriorMeanVariance(M,drop)
+function [mean,variance] = GetPosteriorMeanVariance(M_,drop)
+%function [mean,variance] = GetPosteriorMeanVariance(M,drop)
+% Computes the posterior mean and variance
+% (+updates of oo_ & TeX output).
+%
+% INPUTS
+%   M_               [structure]    Dynare model structure
+%   drop             [double]       share of draws to drop
+%
+% OUTPUTS
+%   mean             [double]       mean parameter vector
+%   variance         [double]       variance 
+%
+% SPECIAL REQUIREMENTS
+%   None.
 
-% Copyright © 2012-2017 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -17,8 +31,8 @@ function [mean,variance] = GetPosteriorMeanVariance(M,drop)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-MetropolisFolder = CheckPath('metropolis',M.dname);
-FileName = M.fname;
+MetropolisFolder = CheckPath('metropolis',M_.dname);
+FileName = M_.fname;
 BaseName = [MetropolisFolder filesep FileName];
 record=load_last_mh_history_file(MetropolisFolder, FileName);
 NbrDraws = sum(record.MhDraws(:,1));
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index dee97ac76..b7713b66d 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -1,17 +1,18 @@
 function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
+% function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
 % This function prints and saves posterior estimates after the mcmc
 % (+updates of oo_ & TeX output).
 %
 % INPUTS
-%   estim_params_    [structure]
-%   M_               [structure]
-%   options_         [structure]
-%   bayestopt_       [structure]
-%   oo_              [structure]
-%   pnames           [cell]        cell of strings, names of the prior shapes available
+%   estim_params_    [structure]    Dynare estimation parameter structure
+%   M_               [structure]    Dynare model structure
+%   options_         [structure]    Dynare options structure
+%   bayestopt_       [structure]    Dynare structure describing priors
+%   oo_              [structure]    Dynare results structure
+%   pnames           [cell]         cell of strings, names of the prior shapes available
 %
 % OUTPUTS
-%   oo_              [structure]
+%   oo_              [structure]    Dynare results structure
 %
 % SPECIAL REQUIREMENTS
 %   None.
@@ -50,7 +51,6 @@ FileName = M_.fname;
 
 record=load_last_mh_history_file(MetropolisFolder,FileName);
 
-FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine = record.KeepedDraws.FirstLine;
 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
@@ -68,12 +68,10 @@ skipline(2)
 disp('ESTIMATION RESULTS')
 skipline()
 
-try
-    disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean))
-catch
-    [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
-    disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean))
+if ~isfield(oo_,'MarginalDensity') || ~isfield(oo_.MarginalDensity,'ModifiedHarmonicMean')
+    [~,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
 end
+fprintf('\nLog data density is %f.\n', oo_.MarginalDensity.ModifiedHarmonicMean);
 
 num_draws=NumberOfDraws*mh_nblck;
 hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
diff --git a/matlab/convergence_diagnostics/McMCDiagnostics.m b/matlab/convergence_diagnostics/McMCDiagnostics.m
index 1d0da5d61..cbe02ce8d 100644
--- a/matlab/convergence_diagnostics/McMCDiagnostics.m
+++ b/matlab/convergence_diagnostics/McMCDiagnostics.m
@@ -1,11 +1,12 @@
 function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
-% function McMCDiagnostics
+% function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
 % Computes convergence tests
 %
 % INPUTS
-%   options_         [structure]
-%   estim_params_    [structure]
-%   M_               [structure]
+%   options_         [structure]    Dynare options structure
+%   estim_params_    [structure]    Dynare estimation parameter structure
+%   M_               [structure]    Dynare model structure
+%   oo_              [structure]    Dynare results structure
 %
 % OUTPUTS
 %   oo_              [structure]
@@ -58,7 +59,6 @@ if issue_an_error_message
 end
 
 % compute inefficiency factor
-FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine = record.KeepedDraws.FirstLine;
 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
@@ -68,6 +68,7 @@ NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws
 param_name = {};
 param_name_tex = {};
 
+Ifac=NaN(nblck,npar);
 for jj = 1:npar
     if options_.TeX
         [par_name_temp, par_name_tex_temp] = get_the_name(jj, options_.TeX, M_,estim_params_, options_);
@@ -107,8 +108,6 @@ update_last_mh_history_file(MetropolisFolder, ModelName, record);
 
 
 PastDraws = sum(record.MhDraws,1);
-LastFileNumber = PastDraws(2);
-LastLineNumber = record.MhDraws(end,3);
 NumberOfDraws  = PastDraws(1);
 
 if NumberOfDraws<=2000
@@ -217,8 +216,6 @@ ALPHA = 0.2;                                % increase too much with the number
 time = 1:NumberOfDraws;
 xx = Origin:StepSize:NumberOfDraws;
 NumberOfLines = length(xx);
-tmp = zeros(NumberOfDraws*nblck,3);
-UDIAG = zeros(NumberOfLines,6,npar);
 
 if NumberOfDraws < Origin
     disp('Estimation::mcmc::diagnostics: The number of simulations is too small to compute the MCMC convergence diagnostics.')
@@ -256,7 +253,6 @@ if isnumeric(options_.parallel)
     clear fout
     % Parallel execution!
 else
-    ModelName = ModelName;
     if ~isempty(M_.bvar)
         ModelName = [ModelName '_bvar'];
     end
@@ -391,7 +387,7 @@ if reste
         fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(pages+1));
         fprintf(fidTeX,'\\end{figure}\n');
         fprintf(fidTeX,'\n');
-        fprintf(fidTeX,'% End Of TeX file.');
+        fprintf(fidTeX,'%% End Of TeX file.');
         fclose(fidTeX);
     end
 end % if reste > 0
diff --git a/matlab/get_name_of_the_last_mh_file.m b/matlab/get_name_of_the_last_mh_file.m
index 6c9ccca24..1ddf7a060 100644
--- a/matlab/get_name_of_the_last_mh_file.m
+++ b/matlab/get_name_of_the_last_mh_file.m
@@ -1,4 +1,5 @@
 function [mhname,info] = get_name_of_the_last_mh_file(M_)
+% function [mhname,info] = get_name_of_the_last_mh_file(M_)
 % This function returns the name of the last mh file and test if the metropolis was completed.
 %
 % INPUTS
@@ -28,7 +29,6 @@ function [mhname,info] = get_name_of_the_last_mh_file(M_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-mhname = [];
 info = 1;
 
 MetropolisFolder = CheckPath('metropolis',M_.dname);
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index 8eb56bed4..48014ba61 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -1,16 +1,16 @@
 function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
-% function marginal = marginal_density()
+% function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
 % Computes the marginal density
 %
 % INPUTS
-%   options_         [structure]
-%   estim_params_    [structure]
-%   M_               [structure]
-%   oo_              [structure]
+%   options_         [structure]    Dynare options structure
+%   estim_params_    [structure]    Dynare estimation parameter structure
+%   M_               [structure]    Dynare model structure
+%   oo_              [structure]    Dynare results structure
 %
 % OUTPUTS
-%   marginal:        [double]     marginal density (modified harmonic mean)
-%   oo_              [structure]
+%   marginal:        [double]       marginal density (modified harmonic mean)
+%   oo_              [structure]    Dynare results structure
 %
 % SPECIAL REQUIREMENTS
 %    none
@@ -42,10 +42,9 @@ record=load_last_mh_history_file(MetropolisFolder, ModelName);
 [nblck, npar] = size(record.LastParameters);
 
 FirstMhFile = record.KeepedDraws.FirstMhFile;
-FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
+FirstLine = record.KeepedDraws.FirstLine;
 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
 TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws);
 
 fprintf('Estimation::marginal density: I''m computing the posterior mean and covariance... ');
diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m
index b5dbe38e4..690dcac5d 100644
--- a/matlab/metropolis_draw.m
+++ b/matlab/metropolis_draw.m
@@ -1,5 +1,5 @@
 function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
-% function [xparams, logpost]=metropolis_draw(init)
+% function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
 % Builds draws from metropolis
 %
 % INPUTS:
@@ -62,8 +62,6 @@ if init
     record=load_last_mh_history_file(MetropolisFolder, FileName);
     FirstMhFile = record.KeepedDraws.FirstMhFile;
     FirstLine = record.KeepedDraws.FirstLine;
-    TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
-    LastMhFile = TotalNumberOfMhFiles;
     TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
     NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
     MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ?
diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m
index d16ff7540..fbc395895 100644
--- a/matlab/posterior_sampler.m
+++ b/matlab/posterior_sampler.m
@@ -3,21 +3,20 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun
 % Random Walk Metropolis-Hastings algorithm.
 %
 % INPUTS
-%   o TargetFun  [char]     string specifying the name of the objective
-%                           function (posterior kernel).
-%   o ProposalFun  [char]   string specifying the name of the proposal
-%                           density
-%   o xparam1    [double]   (p*1) vector of parameters to be estimated (initial values).
-%   o sampler_options       structure
-%            .invhess       [double]   (p*p) matrix, posterior covariance matrix (at the mode).
-%   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters.
-%   o dataset_              data structure
-%   o dataset_info          dataset info structure
-%   o options_              options structure
-%   o M_                    model structure
-%   o estim_params_         estimated parameters structure
-%   o bayestopt_            estimation options structure
-%   o oo_                   outputs structure
+%   o TargetFun         [char]      string specifying the name of the objective
+%                                   function (posterior kernel).
+%   o ProposalFun       [char]      string specifying the name of the proposal
+%                                   density
+%   o xparam1           [double]    (p*1) vector of parameters to be estimated (initial values).
+%   o sampler_options               structure
+%   o mh_bounds         [double]    (p*2) matrix defining lower and upper bounds for the parameters.
+%   o dataset_          [structure] data structure
+%   o dataset_info      [structure] dataset info structure
+%   o options_          [structure] options structure
+%   o M_                [structure] model structure
+%   o estim_params_     [structure] estimated parameters structure
+%   o bayestopt_        [structure] prior specification structure
+%   o oo_               [structure] output structure
 %
 % SPECIAL REQUIREMENTS
 %   None.
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index 9a8190269..2ac7f6963 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -58,16 +58,9 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa
 %Initialize outputs
 ix2 = [];
 ilogpo2 = [];
-ModelName = [];
-MetropolisFolder = [];
 FirstBlock = [];
 FirstLine = [];
-npar = [];
-NumberOfBlocks = [];
-nruns = [];
 NewFile = [];
-MAX_nruns = [];
-d = [];
 
 ModelName = M_.fname;
 if ~isempty(M_.bvar)
@@ -171,7 +164,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
         else
             bayestopt0 = load([PreviousFolder0 filesep 'prior' filesep 'definition.mat']);
         end
-        [common_parameters,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name);
+        [~,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name);
         new_estimated_parameters = ~ismember(bayestopt_.name,bayestopt0.bayestopt_.name);
         ix2 = zeros(NumberOfBlocks,npar);
         ilogpo2 = zeros(NumberOfBlocks,1);
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index 5f2b61db5..51ab6f02d 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -1,6 +1,5 @@
 function prior_posterior_statistics(type,dataset,dataset_info)
-
-% function prior_posterior_statistics(type,dataset)
+% function prior_posterior_statistics(type,dataset,dataset_info)
 % Computes Monte Carlo filter smoother and forecasts
 %
 % INPUTS
@@ -8,6 +7,7 @@ function prior_posterior_statistics(type,dataset,dataset_info)
 %                  prior
 %                  gsa
 %    dataset:      data structure
+%    dataset_info: dataset structure
 %
 % OUTPUTS
 %    none
@@ -208,13 +208,11 @@ localVars.ifil=ifil;
 localVars.DirectoryName = DirectoryName;
 
 if strcmpi(type,'posterior')
-    BaseName = [DirectoryName filesep 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));
-    LastMhFile = TotalNumberOfMhFiles;
     TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
     NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
     mh_nblck = options_.mh_nblck;
@@ -318,9 +316,9 @@ if strcmpi(type,'gsa')
     return
 end
 
-if ~isnumeric(options_.parallel),
+if ~isnumeric(options_.parallel)
     leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen;
-    if options_.parallel_info.leaveSlaveOpen == 0,
+    if options_.parallel_info.leaveSlaveOpen == 0
         % Commenting for testing!!!
         options_.parallel_info.leaveSlaveOpen = 1; % Force locally to leave open remote matlab sessions (repeated pm3 calls)
     end