diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index b7713b66d7eda050588186906be6860c5b64b879..1bec62ce62b695d69b35a10ac0f11a6ebb5c8dd3 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -46,7 +46,7 @@ ncn     = estim_params_.ncn;
 np      = estim_params_.np ;
 
 MetropolisFolder = CheckPath('metropolis',M_.dname);
-OutputFolder = CheckPath('Output',M_.dname);
+latexFolder = CheckPath('latex',M_.dname);
 FileName = M_.fname;
 
 record=load_last_mh_history_file(MetropolisFolder,FileName);
@@ -86,7 +86,7 @@ end
 if np
     type = 'parameters';
     if TeX
-        fid = TeXBegin(OutputFolder, M_.fname, 1, type);
+        fid = TeXBegin(latexFolder, M_.fname, 1, type);
     end
     skipline()
     disp(type)
@@ -129,7 +129,7 @@ end
 if nvx
     type = 'shocks_std';
     if TeX
-        fid = TeXBegin(OutputFolder, FileName,2, 'standard deviation of structural shocks');
+        fid = TeXBegin(latexFolder, FileName,2, 'standard deviation of structural shocks');
     end
     ip = 1;
     skipline()
@@ -173,7 +173,7 @@ end
 if nvn
     type = 'measurement_errors_std';
     if TeX
-        fid = TeXBegin(OutputFolder, FileName, 3, 'standard deviation of measurement errors');
+        fid = TeXBegin(latexFolder, FileName, 3, 'standard deviation of measurement errors');
     end
     skipline()
     disp('standard deviation of measurement errors')
@@ -212,7 +212,7 @@ end
 if ncx
     type = 'shocks_corr';
     if TeX
-        fid = TeXBegin(OutputFolder,FileName,4,'correlation of structural shocks');
+        fid = TeXBegin(latexFolder,FileName,4,'correlation of structural shocks');
     end
     skipline()
     disp('correlation of shocks')
@@ -262,7 +262,7 @@ end
 if ncn
     type = 'measurement_errors_corr';
     if TeX
-        fid = TeXBegin(OutputFolder, FileName, 5, 'correlation of measurement errors');
+        fid = TeXBegin(latexFolder, FileName, 5, 'correlation of measurement errors');
     end
     skipline()
     disp('correlation of measurement errors')
diff --git a/matlab/PlotPosteriorDistributions.m b/matlab/PlotPosteriorDistributions.m
index 5cb0a2d7bba77601a19f4d66917c6292a006d1ba..376810db224926dd0dea36b0127ca90b223aa56a 100644
--- a/matlab/PlotPosteriorDistributions.m
+++ b/matlab/PlotPosteriorDistributions.m
@@ -16,7 +16,7 @@ function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2005-2018 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -32,8 +32,8 @@ function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-OutputDirectoryName = CheckPath('Output',M_.dname);
+latexDirectoryName = CheckPath('latex',M_.dname);
+graphDirectoryName = CheckPath('graphs',M_.dname);
 
 TeX     = options_.TeX;
 nblck   = options_.mh_nblck;
@@ -50,7 +50,7 @@ nn = sqrt(MaxNumberOfPlotPerFigure);
 figurename = 'Priors and posteriors';
 
 if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-    fidTeX = fopen([OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors.tex'],'w');
+    fidTeX = fopen([latexDirectoryName filesep M_.fname '_PriorsAndPosteriors.tex'],'w');
     fprintf(fidTeX,'%% TeX eps-loader file generated by PlotPosteriorDistributions.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
     fprintf(fidTeX,' \n');
@@ -63,7 +63,7 @@ for i=1:npar
     subplotnum = subplotnum+1;
     if subplotnum == 1
         figunumber = figunumber+1;
-        hfig=dyn_figure(options_.nodisplay, 'Name', figurename);
+        hh_fig=dyn_figure(options_.nodisplay, 'Name', figurename);
     end
     [nam,texnam] = get_the_name(i, TeX, M_, estim_params_, options_);
     [x2, f2, abscissa, dens, binf2, bsup2] = draw_prior_density(i, bayestopt_);
@@ -128,8 +128,8 @@ for i=1:npar
     borneinf = min(binf1, binf2);
     bornesup = max(bsup1, bsup2);
     subplot(nn, nn, subplotnum)
-    hh = plot(x2, f2, '-k', 'linewidth', 2);
-    set(hh, 'color', [0.7 0.7 0.7]);
+    hh_fig = plot(x2, f2, '-k', 'linewidth', 2);
+    set(hh_fig, 'color', [0.7 0.7 0.7]);
     hold on;
     plot(x1, f1, '-k', 'linewidth', 2);
     if ~options_.mh_posterior_mode_estimation
@@ -145,12 +145,12 @@ for i=1:npar
     hold off
     drawnow
     if subplotnum == MaxNumberOfPlotPerFigure || i == npar
-        dyn_saveas(hfig,[OutputDirectoryName '/' M_.fname '_PriorsAndPosteriors' int2str(figunumber)], options_.nodisplay, options_.graph_format);
+        dyn_saveas(hh_fig,[graphDirectoryName filesep M_.fname '_PriorsAndPosteriors' int2str(figunumber)], options_.nodisplay, options_.graph_format);
         if TeX && any(strcmp('eps', cellstr(options_.graph_format)))
             fprintf(fidTeX, '\\begin{figure}[H]\n');
             fprintf(fidTeX, '\\centering\n');
             fprintf(fidTeX, '\\includegraphics[width=%2.2f\\textwidth]{%s/%s_PriorsAndPosteriors%s}\n', ...
-                    options_.figures.textwidth*min(subplotnum/nn,1), OutputDirectoryName, M_.fname, int2str(figunumber));
+                    options_.figures.textwidth*min(subplotnum/nn,1), graphDirectoryName, M_.fname, int2str(figunumber));
             fprintf(fidTeX,'\\caption{Priors and posteriors.}');
             fprintf(fidTeX,'\\label{Fig:PriorsAndPosteriors:%s}\n', int2str(figunumber));
             fprintf(fidTeX,'\\end{figure}\n');
diff --git a/matlab/check_posterior_sampler_options.m b/matlab/check_posterior_sampler_options.m
index d00d1913dca838866df440e9775452c6f0bf7153..7848cf6ce563fac6264dfe359195b9f0719b44c1 100644
--- a/matlab/check_posterior_sampler_options.m
+++ b/matlab/check_posterior_sampler_options.m
@@ -1,5 +1,5 @@
-function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_)
-% function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_)
+function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_,outputFolderName)
+% function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_,outputFolderName)
 % initialization of posterior samplers
 %
 % INPUTS
@@ -14,6 +14,7 @@ function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sam
 %   posterior_sampler_options:       checked posterior sampler options
 %   options_:       structure storing the options
 %   bayestopt_:     structure storing information about priors
+%   outputFolderName: string of folder to store mat files
 %
 % SPECIAL REQUIREMENTS
 %   none
@@ -35,6 +36,9 @@ function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sam
 % 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 < 7
+    outputFolderName = 'Output';
+end
 
 init=0;
 if isempty(posterior_sampler_options)
@@ -392,7 +396,7 @@ if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
 end
 
 if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix
-    [~, invhess] = compute_mh_covariance_matrix(bayestopt_,fname,dname);
+    [~, invhess] = compute_mh_covariance_matrix(bayestopt_,fname,dname,outputFolderName);
     posterior_sampler_options.invhess = invhess;
 end
 
diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m
index 3e1abab1568f98e339f6ca934778f5da817dbf7d..e7f0c6a5955d4bf6b259a7d3d41867c8974f32e0 100644
--- a/matlab/compute_mh_covariance_matrix.m
+++ b/matlab/compute_mh_covariance_matrix.m
@@ -1,5 +1,5 @@
-function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname)
-% function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname)
+function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname,outputFolderName)
+% function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname,outputFolderName)
 % Estimation of the posterior covariance matrix, posterior mean, posterior mode and evaluation of the posterior kernel at the
 % estimated mode, using the draws from a metropolis-hastings. The estimated posterior mode and covariance matrix are saved in
 % a file <fname>_mh_mode.mat.
@@ -8,6 +8,7 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at
 %   o  bayestopt_                    [struct]  characterizing priors
 %   o  fname                         [string]  name of model
 %   o  dname                         [string]  name of directory with metropolis folder
+%   o  outputFolderName              [string]  name of directory to store results
 %
 % OUTPUTS
 %   o  posterior_mean                [double]  (n*1) vector, posterior expectation of the parameters.
@@ -34,6 +35,9 @@ 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/>.
+if nargin < 4
+    outputFolderName = 'Output';
+end
 MetropolisFolder = CheckPath('metropolis',dname);
 BaseName = [MetropolisFolder filesep fname];
 
@@ -70,4 +74,4 @@ hh = inv(posterior_covariance);
 fval = posterior_kernel_at_the_mode;
 parameter_names = bayestopt_.name;
 
-save([dname filesep 'Output' filesep fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names');
\ No newline at end of file
+save([dname filesep outputFolderName filesep fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names');
\ No newline at end of file
diff --git a/matlab/convergence_diagnostics/mcmc_diagnostics.m b/matlab/convergence_diagnostics/mcmc_diagnostics.m
index c61de8fc6c64b5f2d36530030babee055b4da585..93e19b0bcd9001ec53bfbf8d305cc252f1608518 100644
--- a/matlab/convergence_diagnostics/mcmc_diagnostics.m
+++ b/matlab/convergence_diagnostics/mcmc_diagnostics.m
@@ -33,8 +33,8 @@ function oo_ = mcmc_diagnostics(options_, estim_params_, M_, oo_)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-OutputFolder = CheckPath('Output',M_.dname);
+graphFolder = CheckPath('graphs',M_.dname);
+latexFolder = CheckPath('latex',M_.dname);
 MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
 
@@ -228,7 +228,7 @@ if NumberOfDraws < Origin
 end
 
 if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-    fidTeX = fopen([OutputFolder '/' ModelName '_UnivariateDiagnostics.tex'],'w');
+    fidTeX = fopen([latexFolder '/' ModelName '_UnivariateDiagnostics.tex'],'w');
     fprintf(fidTeX,'%% TeX eps-loader file generated by mcmc_diagnostics.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
     fprintf(fidTeX,' \n');
@@ -277,7 +277,7 @@ clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp;
 pages = floor(npar/npardisp); % changed from 3 to npardisp
 k = 0;
 for i = 1:pages
-    h=dyn_figure(options_.nodisplay,'Name','MCMC univariate convergence diagnostic (Brooks and Gelman,1998)');
+    hh_fig = dyn_figure(options_.nodisplay,'Name','MCMC univariate convergence diagnostic (Brooks and Gelman,1998)');
     boxplot = 1;
     for j = 1:npardisp % Loop over parameters  %npardisp instead of 3
         k = k+1;
@@ -311,11 +311,11 @@ for i = 1:pages
             boxplot = boxplot + 1;
         end
     end
-    dyn_saveas(h,[OutputFolder '/' ModelName '_udiag' int2str(i)],options_.nodisplay,options_.graph_format);
+    dyn_saveas(hh_fig,[graphFolder '/' ModelName '_udiag' int2str(i)],options_.nodisplay,options_.graph_format);
     if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         fprintf(fidTeX,'\\begin{figure}[H]\n');
         fprintf(fidTeX,'\\centering \n');
-        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n',options_.figures.textwidth*min((boxplot-1)/3,1),[OutputFolder '/' ModelName],int2str(i));
+        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n',options_.figures.textwidth*min((boxplot-1)/3,1),[graphFolder '/' ModelName],int2str(i));
         fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n');
         fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n');
         fprintf(fidTeX,'the eighty percent interval, the second and third moments.}');
@@ -333,7 +333,7 @@ if reste
         nr = npardisp;
         nc = 3;
     end
-    h = dyn_figure(options_.nodisplay,'Name','MCMC univariate convergence diagnostic (Brooks and Gelman, 1998)');
+    hh_fig = dyn_figure(options_.nodisplay,'Name','MCMC univariate convergence diagnostic (Brooks and Gelman, 1998)');
     boxplot = 1;
     for j = 1:reste
         k = k+1;
@@ -375,11 +375,11 @@ if reste
             boxplot = boxplot + 1;
         end
     end
-    dyn_saveas(h,[ OutputFolder '/' ModelName '_udiag' int2str(pages+1)],options_.nodisplay,options_.graph_format);
+    dyn_saveas(hh_fig,[ graphFolder '/' ModelName '_udiag' int2str(pages+1)],options_.nodisplay,options_.graph_format);
     if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         fprintf(fidTeX,'\\begin{figure}[H]\n');
         fprintf(fidTeX,'\\centering \n');
-        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n',options_.figures.textwidth*min((boxplot-1)/nc,1),[OutputFolder '/' ModelName],int2str(pages+1));
+        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n',options_.figures.textwidth*min((boxplot-1)/nc,1),[graphFolder '/' ModelName],int2str(pages+1));
         if reste == 2
             fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n');
             fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n');
@@ -401,7 +401,7 @@ clear UDIAG;
 % Multivariate diagnostic.
 %
 if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-    fidTeX = fopen([OutputFolder '/' ModelName '_MultivariateDiagnostics.tex'],'w');
+    fidTeX = fopen([latexFolder '/' ModelName '_MultivariateDiagnostics.tex'],'w');
     fprintf(fidTeX,'%% TeX eps-loader file generated by mcmc_diagnostics.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
     fprintf(fidTeX,' \n');
@@ -445,7 +445,7 @@ for iter  = Origin:StepSize:NumberOfDraws
 end
 MDIAG(:,[2 4 6],:) = MDIAG(:,[2 4 6],:)/nblck;
 
-h = dyn_figure(options_.nodisplay,'Name','Multivariate convergence diagnostic');
+hh_fig = dyn_figure(options_.nodisplay,'Name','Multivariate convergence diagnostic');
 boxplot = 1;
 for crit = 1:3
     if crit == 1
@@ -470,12 +470,12 @@ for crit = 1:3
     title(namnam,'Interpreter','none');
     boxplot = boxplot + 1;
 end
-dyn_saveas(h,[ OutputFolder '/' ModelName '_mdiag'],options_.nodisplay,options_.graph_format);
+dyn_saveas(hh_fig,[ graphFolder '/' ModelName '_mdiag'],options_.nodisplay,options_.graph_format);
 
 if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
     fprintf(fidTeX,'\\begin{figure}[H]\n');
     fprintf(fidTeX,'\\centering \n');
-    fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_mdiag}\n',[OutputFolder '/' ModelName]);
+    fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_mdiag}\n',[graphFolder '/' ModelName]);
     fprintf(fidTeX,'\\caption{Multivariate convergence diagnostics for the Metropolis-Hastings.\n');
     fprintf(fidTeX,'The first, second and third rows are respectively the criteria based on\n');
     fprintf(fidTeX,'the eighty percent interval, the second and third moments. The different \n');
diff --git a/matlab/display_estimation_results_table.m b/matlab/display_estimation_results_table.m
index 4d0b7031a45460e27b5384d4456cae99d6e496ba..6078d50317506c955c0ff39556dd245febbf3495 100644
--- a/matlab/display_estimation_results_table.m
+++ b/matlab/display_estimation_results_table.m
@@ -20,7 +20,7 @@ function oo_=display_estimation_results_table(xparam1,stdh,M_,options_,estim_par
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright © 2014-2018 Dynare Team
+% Copyright © 2014-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -177,11 +177,11 @@ if any(xparam1(1:nvx+nvn)<0)
     warning(sprintf('Some estimated standard deviations are negative.\n         Dynare internally works with variances so that the sign does not matter.\n         Nevertheless, it is recommended to impose either prior restrictions (Bayesian Estimation)\n         or a lower bound (ML) to assure positive values.'))
 end
 
-OutputDirectoryName = CheckPath('Output',M_.dname);
+latexDirectoryName = CheckPath('latex',M_.dname);
 
 if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior mode) Latex output
     if np
-        filename = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_1.tex'];
+        filename = [latexDirectoryName '/' M_.fname '_Posterior_Mode_1.tex'];
         fidTeX = fopen(filename,'w');
         TeXBegin_Bayesian(fidTeX,1,'parameters')
         ip = nvx+nvn+ncx+ncn+1;
@@ -198,7 +198,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
         TeXEnd(fidTeX)
     end
     if nvx
-        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_2.tex'];
+        TeXfile = [latexDirectoryName '/' M_.fname '_Posterior_Mode_2.tex'];
         fidTeX = fopen(TeXfile,'w');
         TeXBegin_Bayesian(fidTeX,2,'standard deviation of structural shocks')
         ip = 1;
@@ -216,7 +216,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
         TeXEnd(fidTeX)
     end
     if nvn
-        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_3.tex'];
+        TeXfile = [latexDirectoryName '/' M_.fname '_Posterior_Mode_3.tex'];
         fidTeX  = fopen(TeXfile,'w');
         TeXBegin_Bayesian(fidTeX,3,'standard deviation of measurement errors')
         ip = nvx+1;
@@ -234,7 +234,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
         TeXEnd(fidTeX)
     end
     if ncx
-        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_4.tex'];
+        TeXfile = [latexDirectoryName '/' M_.fname '_Posterior_Mode_4.tex'];
         fidTeX = fopen(TeXfile,'w');
         TeXBegin_Bayesian(fidTeX,4,'correlation of structural shocks')
         ip = nvx+nvn+1;
@@ -253,7 +253,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
         TeXEnd(fidTeX)
     end
     if ncn
-        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_5.tex'];
+        TeXfile = [latexDirectoryName '/' M_.fname '_Posterior_Mode_5.tex'];
         fidTeX = fopen(TeXfile,'w');
         TeXBegin_Bayesian(fidTeX,5,'correlation of measurement errors')
         ip = nvx+nvn+ncx+1;
@@ -273,7 +273,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
     end
 elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
     if np
-        filename = [OutputDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_1.tex'];
+        filename = [latexDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_1.tex'];
         fidTeX = fopen(filename, 'w');
         TeXBegin_ML(fidTeX, 1, 'parameters', table_title, LaTeXtitle)
         ip = nvx+nvn+ncx+ncn+1;
@@ -288,7 +288,7 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
         TeXEnd(fidTeX)
     end
     if nvx
-        filename = [OutputDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_2.tex'];
+        filename = [latexDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_2.tex'];
         fidTeX = fopen(filename, 'w');
         TeXBegin_ML(fidTeX, 2, 'standard deviation of structural shocks', table_title, LaTeXtitle)
         ip = 1;
@@ -304,7 +304,7 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
         TeXEnd(fidTeX)
     end
     if nvn
-        filename = [OutputDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_3.tex'];
+        filename = [latexDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_3.tex'];
         fidTeX = fopen(filename, 'w');
         TeXBegin_ML(fidTeX, 3, 'standard deviation of measurement errors', table_title, LaTeXtitle)
         ip = nvx+1;
@@ -320,7 +320,7 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
         TeXEnd(fidTeX)
     end
     if ncx
-        filename = [OutputDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_4.tex'];
+        filename = [latexDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_4.tex'];
         fidTeX = fopen(filename, 'w');
         TeXBegin_ML(fidTeX, 4, 'correlation of structural shocks', table_title,LaTeXtitle)
         ip = nvx+nvn+1;
@@ -337,7 +337,7 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
         TeXEnd(fidTeX)
     end
     if ncn
-        filename = [OutputDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_5.tex'];
+        filename = [latexDirectoryName '/' M_.fname '_' LaTeXtitle '_Mode_5.tex'];
         fidTeX = fopen(filename, 'w');
         TeXBegin_ML(fidTeX, 5, 'correlation of measurement errors', table_title, LaTeXtitle)
         ip = nvx+nvn+ncx+1;
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index 2f8c45bd76d9565c84d26a57ec804fc0dd9d253e..dc6826e59587626c26e91bb60f24b079dcac4d20 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -1,5 +1,5 @@
-function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
-% function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
+function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_, outputFolderName)
+% function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_, outputFolderName)
 % Computes the marginal density
 %
 % INPUTS
@@ -7,6 +7,7 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay
 %   estim_params_    [structure]    Dynare estimation parameter structure
 %   M_               [structure]    Dynare model structure
 %   oo_              [structure]    Dynare results structure
+%   outputFolderName [string]       name of folder with results
 %
 % OUTPUTS
 %   marginal:        [double]       marginal density (modified harmonic mean)
@@ -31,7 +32,9 @@ function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bay
 %
 % 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 < 6
+    outputFolderName = 'Output';
+end
 
 MetropolisFolder = CheckPath('metropolis',M_.dname);
 ModelName = M_.fname;
@@ -47,8 +50,8 @@ TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
 TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws);
 
-fprintf('Estimation::marginal density: I''m computing the posterior mean and covariance... ');
-[posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname);
+fprintf('marginal density: I''m computing the posterior mean and covariance... ');
+[posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname,outputFolderName);
 
 MU = transpose(posterior_mean);
 SIGMA = posterior_covariance;
@@ -64,7 +67,7 @@ end
 % (usefull if the user wants to perform some computations using
 % the posterior mean instead of the posterior mode ==> ).
 parameter_names = bayestopt_.name;
-save([M_.dname filesep 'Output' filesep M_.fname '_mean.mat'],'xparam1','hh','parameter_names','SIGMA');
+save([M_.dname filesep outputFolderName filesep M_.fname '_mean.mat'],'xparam1','hh','parameter_names','SIGMA');
 
 fprintf('Estimation::marginal density: I''m computing the posterior log marginal density (modified harmonic mean)... ');
 try 
diff --git a/matlab/plot_priors.m b/matlab/plot_priors.m
index 35759f3093fb00de7e69264f4e4dc395a2f6bbec..9cfa5385ee37dfa25f4774665023c599e3c5933c 100644
--- a/matlab/plot_priors.m
+++ b/matlab/plot_priors.m
@@ -15,7 +15,7 @@ function plot_priors(bayestopt_,M_,estim_params_,options_,optional_title)
 % SPECIAL REQUIREMENTS
 %    None
 
-% Copyright © 2004-2020 Dynare Team
+% Copyright © 2004-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -31,6 +31,8 @@ function plot_priors(bayestopt_,M_,estim_params_,options_,optional_title)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+latexDirectoryName = CheckPath('latex',M_.dname);
+graphDirectoryName = CheckPath('graphs',M_.dname);
 
 TeX = options_.TeX;
 if nargin<5
@@ -41,17 +43,14 @@ end
 npar = length(bayestopt_.p1);
 [nbplt,nr,nc,lr,lc,nstar] = pltorg(npar);
 
-if ~exist([M_.dname '/graphs'],'dir')
-    mkdir(M_.dname,'graphs');
-end
 if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-    fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_Priors.tex'],'w');
+    fidTeX = fopen([latexDirectoryName filesep M_.fname '_Priors.tex'],'w');
     fprintf(fidTeX,'%% TeX eps-loader file generated by plot_priors.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
     fprintf(fidTeX,' \n');
 end
 for plt = 1:nbplt
-    hplt = dyn_figure(options_.nodisplay,'Name',figurename);
+    hh_fig = dyn_figure(options_.nodisplay,'Name',figurename);
     if TeX
         TeXNAMES = [];
         NAMES    = [];
@@ -63,8 +62,8 @@ for plt = 1:nbplt
         [x,f,abscissa,dens,binf,bsup] = draw_prior_density(i,bayestopt_);
         [nam,texnam] = get_the_name(i,TeX,M_,estim_params_,options_);
         subplot(nr,nc,index)
-        hh = plot(x,f,'-k','linewidth',2);
-        set(hh,'color',[0.7 0.7 0.7]);
+        hh_plt = plot(x,f,'-k','linewidth',2);
+        set(hh_plt,'color',[0.7 0.7 0.7]);
         box on
         if TeX
             title(texnam,'Interpreter','latex')
@@ -73,11 +72,11 @@ for plt = 1:nbplt
         end
         drawnow
     end
-    dyn_saveas(hplt,[M_.dname, '/graphs/' M_.fname '_Priors' int2str(plt)],options_.nodisplay,options_.graph_format);
+    dyn_saveas(hh_fig,[graphDirectoryName filesep M_.fname '_Priors' int2str(plt)],options_.nodisplay,options_.graph_format);
     if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         fprintf(fidTeX,'\\begin{figure}[H]\n');
         fprintf(fidTeX,'\\centering\n');
-        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_Priors%s}\n',options_.figures.textwidth*min(index/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
+        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_Priors%s}\n',options_.figures.textwidth*min(index/nc,1),[graphDirectoryName filesep M_.fname],int2str(plt));
         fprintf(fidTeX,'\\caption{Priors.}');
         fprintf(fidTeX,'\\label{Fig:Priors:%s}\n',int2str(plt));
         fprintf(fidTeX,'\\end{figure}\n');