diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index de04e839850365b5f1b7b9df97d4f62672a35564..93f5e2335f432b92331d6772e5961719f49adaf0 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -733,7 +733,7 @@ if do_bayesian_estimation_mcmc
                 oo_.mom.posterior.metropolis = oo_load_mh.oo_.mom.posterior.metropolis;
             end
         end
-        [error_flag,~,options_mom_]= metropolis_draw(1,options_mom_,estim_params_,M_);
+        [options_mom_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_mom_); %check whether number of sub_draws is feasible
         if ~(~isempty(options_mom_.sub_draws) && options_mom_.sub_draws==0)
             % THIS IS PROBABLY NOT USEFUL HERE AND CAN BE REMOVED (PREPROCESSOR: REMOVE bayesian_irf, moments_varendo)
             %if options_mom_.bayesian_irf
diff --git a/matlab/convergence_diagnostics/mcmc_diagnostics.m b/matlab/convergence_diagnostics/mcmc_diagnostics.m
index 6872074fbecf75cd1a2df41bde77c1e278760b2e..ad216be37f4bfd6eaee218142d97b1bf8bd2420d 100644
--- a/matlab/convergence_diagnostics/mcmc_diagnostics.m
+++ b/matlab/convergence_diagnostics/mcmc_diagnostics.m
@@ -59,16 +59,14 @@ if issue_an_error_message
 end
 
 % compute inefficiency factor
-FirstLine = record.KeepedDraws.FirstLine;
-TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-FirstMhFile = record.KeepedDraws.FirstMhFile;
 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
 
 param_name = {};
 param_name_tex = {};
 
 Ifac=NaN(nblck,npar);
+AllDraws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, 'all');
 for jj = 1:npar
     if options_.TeX
         [par_name_temp, par_name_tex_temp] = get_the_name(jj, options_.TeX, M_,estim_params_, options_.varobs);
@@ -79,8 +77,7 @@ for jj = 1:npar
         par_name_temp = get_the_name(jj, options_.TeX, M_, estim_params_, options_.varobs);
         param_name = vertcat(param_name, par_name_temp);
     end
-    Draws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
-    Draws = reshape(Draws, [NumberOfDraws nblck]);
+    Draws = reshape(AllDraws(:,jj), [NumberOfDraws nblck]);
     Nc = min(1000, NumberOfDraws/2);
     for ll = 1:nblck
         Ifac(ll,jj) = mcmc_ifac(Draws(:,ll), Nc);
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 95d932248fe8cb231a21026171de81146c0fac47..8de89f5f7e6c0b4c00bbae1bcf58fedbe061add9 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -659,7 +659,6 @@ options_.particleswarm = particleswarm;
 
 % prior analysis
 options_.prior_mc = 20000;
-options_.prior_analysis_endo_var_list = {};
 
 % did model undergo block decomposition + minimum feedback set computation ?
 options_.block = false;
diff --git a/matlab/estimation/CutSample.m b/matlab/estimation/CutSample.m
index 462aa070a4c3906dd6537c2063f0c3902796cb28..4cbbbc887712e44db34d183ce6e1240ccdeeb270 100644
--- a/matlab/estimation/CutSample.m
+++ b/matlab/estimation/CutSample.m
@@ -46,7 +46,7 @@ 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' ]);
 
-if ~length(mh_files)
+if isempty(mh_files)
     error('%s: I can''t find MH file to load here!',dispString)
 end
 
diff --git a/matlab/estimation/GetAllPosteriorDraws.m b/matlab/estimation/GetAllPosteriorDraws.m
index 9a4b67b3e44d31fb34931ce82c1ea4a5d84c9596..5c1847cf0283ecf87da59089dd78f843a7306ab4 100644
--- a/matlab/estimation/GetAllPosteriorDraws.m
+++ b/matlab/estimation/GetAllPosteriorDraws.m
@@ -1,26 +1,22 @@
-function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
-
+function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstLine , FirstMhFile, blck)
+% function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstLine , FirstMhFile, blck)
 % Gets all posterior draws.
 %
 % INPUTS
-% - options_               [struct]   Dynare's options.
-% - dname                  [char]     name of directory with results.
-% - fname                  [char]     name of mod file.
-% - column                 [integer]  scalar, parameter index.
-% - FirstMhFile            [integer]  scalar, first MH file.
-% - FirstLine              [integer]  scalar, first line in first MH file.
-% - TotalNumberOfMhFile    [integer]  scalar, total number of MH file.
-% - NumberOfDraws          [integer]  scalar, number of posterior draws.
-% - nblcks                 [integer]  scalar, total number of blocks.
-% - blck:                  [integer]  scalar, desired block to read.
+% - options_               [struct]             Dynare's options.
+% - dname                  [char]               name of directory with results.
+% - fname                  [char]               name of mod file.
+% - column                 [integer or string]  scalar, parameter index
+%                                               'all' if all parameters are 
+%                                               requested
+% - FirstMhFile            [integer]            optional scalar, first MH file to be used
+% - FirstLine              [integer]            optional scalar, first line in first MH file to be used
+% - blck:                  [integer]  scalar, desired block to read (optional)
 %
 % OUTPUTS
 % - draws:                 [double]   NumberOfDraws×1 vector, draws from posterior distribution.
-%
-% REMARKS
-% Only the first and third input arguments are required for SMC samplers.
 
-% Copyright © 2005-2023 Dynare Team
+% Copyright © 2005-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,71 +33,89 @@ function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstMhFil
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if ishssmc(options_)
-    % Load draws from the posterior distribution
-    pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
-    posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
-    if column==0
-        draws = posterior.tlogpostkernel;
-    else
-        draws = transpose(posterior.particles(column,:));
-    end
-elseif isdime(options_)
-    posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
-    tune = posterior.tune;
-    chains = posterior.chains(end-tune:end,:,:);
-    if column>0
-        chains = reshape(chains, [], size(chains, 3));
-        draws = chains(:,column);
+if isnumeric(column) && column<0
+    error('GetAllPosteriorDraws:: column cannot be negative');
+end
+
+if issmc(options_)
+    if ishssmc(options_)
+        % Load draws from the posterior distribution
+        pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
+        posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
+        if column>0 || strcmp(column,'all')
+            if strcmp(column,'all')
+                draws = transpose(posterior.particles);
+            else
+                draws = transpose(posterior.particles(column,:));
+            end
+        else
+            draws = posterior.tlogpostkernel;
+        end
+    elseif isdime(options_)
+        posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
+        tune = posterior.tune;
+        chains = posterior.chains(end-tune:end,:,:);
+        if column>0 || strcmp(column,'all')
+            chains = reshape(chains, [], size(chains, 3));
+            if strcmp(column,'all')
+                draws = chains;
+            else
+                draws = chains(:,column);
+            end
+        else
+            draws = posterior.lprobs;
+        end
     else
-        draws = posterior.lprobs;
+        error('GetAllPosteriorDraws:: case should not happen. Please contact the developers')
     end
 else
+    DirectoryName = CheckPath('metropolis',dname);
+    record=load_last_mh_history_file(DirectoryName,fname);
+    if nargin<5 || isempty(FirstLine)
+        FirstLine = record.KeepedDraws.FirstLine;
+    end
+    if nargin<6 || isempty(FirstMhFile)
+        FirstMhFile = record.KeepedDraws.FirstMhFile;
+    end
+    TotalNumberOfMhFile = sum(record.MhDraws(:,2));
+    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+    NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+    [nblcks, npar] = size(record.LastParameters);
     iline = FirstLine;
     linee = 1;
-    DirectoryName = CheckPath('metropolis',dname);
-    if nblcks>1 && nargin<10
-        draws = zeros(NumberOfDraws*nblcks,1);
-        iline0=iline;
-        if column>0
-            for blck = 1:nblcks
-                iline=iline0;
-                for file = FirstMhFile:TotalNumberOfMhFile
-                    load([DirectoryName '/'  fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
-                    NumberOfLines = size(x2(iline:end,:),1);
-                    draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
-                    linee = linee+NumberOfLines;
-                    iline = 1;
-                end
-            end
-        else
-            for blck = 1:nblcks
-                iline=iline0;
-                for file = FirstMhFile:TotalNumberOfMhFile
-                    load([DirectoryName '/'  fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2')
-                    NumberOfLines = size(logpo2(iline:end),1);
-                    draws(linee:linee+NumberOfLines-1) = logpo2(iline:end);
-                    linee = linee+NumberOfLines;
-                    iline = 1;
-                end
-            end
-        end
+    if nargin==7
+        blocks_to_load=blck;
+        nblcks=length(blck);
     else
-        if nblcks==1
-            blck=1;
-        end
-        if column>0
+        blocks_to_load=1:nblcks;
+    end
+    if strcmp(column,'all')
+        draws = zeros(NumberOfDraws*nblcks,npar);
+    else
+        draws = zeros(NumberOfDraws*nblcks,1);
+    end
+    iline0=iline;
+    if column>0
+        for blck = blocks_to_load
+            iline=iline0;
             for file = FirstMhFile:TotalNumberOfMhFile
                 load([DirectoryName '/'  fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
                 NumberOfLines = size(x2(iline:end,:),1);
-                draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
+                if strcmp(column,'all')
+                    draws(linee:linee+NumberOfLines-1,:) = x2(iline:end,:);
+                else
+                    draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
+                end
                 linee = linee+NumberOfLines;
                 iline = 1;
             end
-        else
+        end
+    else
+        for blck = blocks_to_load
+            iline=iline0;
             for file = FirstMhFile:TotalNumberOfMhFile
                 load([DirectoryName '/'  fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2')
-                NumberOfLines = size(logpo2(iline:end,:),1);
+                NumberOfLines = size(logpo2(iline:end),1);
                 draws(linee:linee+NumberOfLines-1) = logpo2(iline:end);
                 linee = linee+NumberOfLines;
                 iline = 1;
diff --git a/matlab/estimation/GetOneDraw.m b/matlab/estimation/GetOneDraw.m
deleted file mode 100644
index 8c92a84ba70dac74984995a35d658ca067dd4c54..0000000000000000000000000000000000000000
--- a/matlab/estimation/GetOneDraw.m
+++ /dev/null
@@ -1,47 +0,0 @@
-function [xparams, logpost] = GetOneDraw(distribution, M_, estim_params_, oo_, options_, bayestopt_)
-
-% draws one parameter vector and its posterior from MCMC or the prior
-%
-% INPUTS
-%    type:      [string]       'posterior': draw from MCMC draws
-%                              'prior': draw from prior
-%    M_         [structure]     Definition of the model
-%    estim_params_ [structure]  characterizing parameters to be estimated
-%    oo_         [structure]    Storage of results
-%    options_    [structure]    Options
-%    bayestopt_  [structure]    describing the priors
-%
-% OUTPUTS
-%    xparams:   vector of estimated parameters (drawn from posterior or prior distribution)
-%    logpost:   log of the posterior density of this parameter vector
-%
-% SPECIAL REQUIREMENTS
-%    none
-
-% Copyright © 2005-2023 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-if isprior(distribution)
-    xparams = distribution.draw();
-    if nargout>1
-        logpost = evaluate_posterior_kernel(xparams, M_, estim_params_, oo_, options_, bayestopt_);
-    end
-elseif ischar(distribution) && strcmpi(distribution, 'posterior')
-    [xparams, logpost] = metropolis_draw(0);
-else
-    error('GetOneDraw:: Wrong inputs.')
-end
diff --git a/matlab/estimation/GetPosteriorMeanVariance.m b/matlab/estimation/GetPosteriorMeanVariance.m
index 6cb5a7f4ab5eaaec2ee0e70ba183e93b7b401115..9a1c1205129f6cc24b7a0745b15e92ec94db039b 100644
--- a/matlab/estimation/GetPosteriorMeanVariance.m
+++ b/matlab/estimation/GetPosteriorMeanVariance.m
@@ -1,7 +1,6 @@
 function [mean, variance] = GetPosteriorMeanVariance(options_, M_)
 % [mean,variance] = GetPosteriorMeanVariance(options_, M_)
 % Computes the posterior mean and variance
-% (+updates of oo_ & TeX output).
 %
 % INPUTS
 % - options_         [struct]    Dynare's options.
@@ -28,14 +27,25 @@ function [mean, variance] = GetPosteriorMeanVariance(options_, M_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if ishssmc(options_)
-    % Load draws from the posterior distribution
-    pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
-    posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
-    % Compute the posterior mean
-    mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
-    % Compute the posterior covariance
-    variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
+if issmc(options_)
+    if ishssmc(options_)
+        % Load draws from the posterior distribution
+        pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+        posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+        % Compute the posterior mean
+        mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
+        % Compute the posterior covariance
+        variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
+    elseif isdime(options_)
+        posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+        tune = posterior.tune;
+        chains = transpose(reshape(posterior.chains(end-tune:end,:,:), [], size(posterior.chains, 3)));
+        mean = sum(chains, 2)/size(chains,1);
+        % Compute the posterior covariance
+        variance = (chains-mean)*(chains-mean)'/size(chains,1);
+    else
+        error('GetPosteriorMeanVariance:: case should not happen. Please contact the developers')
+    end
 else
     MetropolisFolder = CheckPath('metropolis',M_.dname);
     FileName = M_.fname;
diff --git a/matlab/estimation/GetPosteriorParametersStatistics.m b/matlab/estimation/GetPosteriorParametersStatistics.m
index 299ccdaca07cc4fafb3cc257e63b3c69336d3bc3..208a37779290ca18a3daa690cfb706c9f7e419b8 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -44,17 +44,7 @@ np      = estim_params_.np ;
 latexFolder = CheckPath('latex',M_.dname);
 FileName = M_.fname;
 
-if ~issmc(options_)
-    MetropolisFolder = CheckPath('metropolis',M_.dname);
-    record=load_last_mh_history_file(MetropolisFolder,FileName);
-    FirstLine = record.KeepedDraws.FirstLine;
-    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;
-end
+[~, ~, num_draws]=set_number_of_subdraws(M_,options_); %get number of draws
 
 header_width = row_header_width(M_, estim_params_, bayestopt_);
 hpd_interval=[num2str(options_.mh_conf_sig*100), '% HPD interval'];
@@ -67,31 +57,18 @@ skipline()
 
 if ishssmc(options_)
     dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
-    % Set function handle for GetAllPosteriorDraws
-    getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, [], i);
+    hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
 elseif isdime(options_)
-    getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, [], i);
+    hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
 else
     if ~isfield(oo_,'MarginalDensity') || (issmc(options_) && ~isfield(oo_.MarginalDensity,'ModifiedHarmonicMean'))
         [~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
     end
     fprintf('Log data density (Modified Harmonic Mean) is %f.', oo_.MarginalDensity.ModifiedHarmonicMean);
-    % Set function handle for GetAllPosteriordraws
-    getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, M_.fname, i, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
-end
-
-if ishssmc(options_)
-    num_draws = options_.posterior_sampler_options.hssmc.particles;
-    hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
-elseif isdime(options_)
-    nchain = options_.posterior_sampler_options.current_options.nchain;
-    tune = options_.posterior_sampler_options.current_options.tune;
-    num_draws = nchain*tune;
-    hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
-else
-    num_draws=NumberOfDraws*mh_nblck;
     hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
 end
+% Set function handle for GetAllPosteriorDraws
+getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, M_.fname, i);
 
 if hpd_draws<2
     fprintf('posterior_moments: There are not enough draws computes to compute HPD Intervals. Skipping their computation.\n')
diff --git a/matlab/estimation/PosteriorIRF.m b/matlab/estimation/PosteriorIRF.m
index d53bc344d17795840f4a8fe34d2bf9d7dc5c3e0f..e5a9861ef7a727c60dacbf60db77bfb3352fa226 100644
--- a/matlab/estimation/PosteriorIRF.m
+++ b/matlab/estimation/PosteriorIRF.m
@@ -25,7 +25,7 @@ function oo_=PosteriorIRF(type,options_,estim_params_,oo_,M_,bayestopt_,dataset_
 % functions associated with it(the _core1 and _core2).
 % See also the comments posterior_sampler.m funtion.
 
-% Copyright © 2006-2023 Dynare Team
+% Copyright © 2006-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -83,12 +83,13 @@ end
 
 DirectoryName = CheckPath('Output',M_.dname);
 if strcmpi(type,'posterior')
-    MhDirectoryName = CheckPath('metropolis',M_.dname);
+    folder_name=get_posterior_folder_name(options_);
+    MhDirectoryName = CheckPath(folder_name,M_.dname);
 elseif strcmpi(type,'gsa')
     if options_.opt_gsa.pprior
-        MhDirectoryName = CheckPath(['GSA' filesep 'prior'],M_.dname);
+        MhDirectoryName = CheckPath(['gsa' filesep 'prior'],M_.dname);
     else
-        MhDirectoryName = CheckPath(['GSA' filesep 'mc'],M_.dname);
+        MhDirectoryName = CheckPath(['gsa' filesep 'mc'],M_.dname);
     end
 else
     MhDirectoryName = CheckPath('prior',M_.dname);
@@ -156,10 +157,7 @@ localVars.npar = npar;
 
 localVars.type=type;
 if strcmpi(type,'posterior')
-    while b<B
-        b = b + 1;
-        x(b,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
-    end
+    x=get_posterior_subsample(M_,options_,B);
 end
 
 if ~strcmpi(type,'prior')
diff --git a/matlab/estimation/PosteriorIRF_core1.m b/matlab/estimation/PosteriorIRF_core1.m
index 096c12c78e19ed0ff4ee39397f9668e930fcb80b..fa29c9f4b23730e978ce6e2ab47334b863016c96 100644
--- a/matlab/estimation/PosteriorIRF_core1.m
+++ b/matlab/estimation/PosteriorIRF_core1.m
@@ -88,17 +88,7 @@ if whoiam
     Parallel=myinputs.Parallel;
 end
 
-if strcmpi(type,'posterior')
-    MhDirectoryName = CheckPath('metropolis',M_.dname);
-elseif strcmpi(type,'gsa')
-    if options_.opt_gsa.pprior
-        MhDirectoryName = CheckPath(['gsa' filesep 'prior'],M_.dname);
-    else
-        MhDirectoryName = CheckPath(['gsa' filesep 'mc'],M_.dname);
-    end
-else
-    MhDirectoryName = CheckPath('prior',M_.dname);
-end
+MhDirectoryName = myinputs.MhDirectoryName;
 
 RemoteFlag = 0;
 
diff --git a/matlab/estimation/ReshapeMatFiles.m b/matlab/estimation/ReshapeMatFiles.m
index 58107602eb15542bb828e7cbb244dc057b818580..a4051bfb416d4b92cc3b0b9832250bf2f13db18c 100644
--- a/matlab/estimation/ReshapeMatFiles.m
+++ b/matlab/estimation/ReshapeMatFiles.m
@@ -30,7 +30,7 @@ function ReshapeMatFiles(fname, dname, exo_nbr, endo_nbr, options_, type, type2)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2003-2023 Dynare Team
+% Copyright © 2003-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -48,10 +48,12 @@ function ReshapeMatFiles(fname, dname, exo_nbr, endo_nbr, options_, type, type2)
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if nargin==6
-    MhDirectoryName = [ CheckPath('metropolis',dname) filesep ];
+    folder_name=get_posterior_folder_name(options_);
+    MhDirectoryName = [CheckPath(folder_name,dname) filesep ];
 else
     if strcmpi(type2,'posterior')
-        MhDirectoryName = [CheckPath('metropolis',dname) filesep ];
+        folder_name=get_posterior_folder_name(options_);
+        MhDirectoryName = [CheckPath(folder_name,dname) filesep ];
     elseif strcmpi(type2,'gsa')
         if options_.opt_gsa.morris==1
             MhDirectoryName = [CheckPath('gsa/screen',dname) filesep ];
diff --git a/matlab/estimation/check_posterior_analysis_data.m b/matlab/estimation/check_posterior_analysis_data.m
index c8a8c71f2bb47d5c155a16305d5b435be6090a9d..6de65d2de26ed4380f63f3785b3ddc4dc5f5a467 100644
--- a/matlab/estimation/check_posterior_analysis_data.m
+++ b/matlab/estimation/check_posterior_analysis_data.m
@@ -1,11 +1,12 @@
-function [info,description] = check_posterior_analysis_data(type,M_)
-% function [info,description] = check_posterior_analysis_data(type,M_)
+function [info,description] = check_posterior_analysis_data(type,M_,options_)
+% function [info,description] = check_posterior_analysis_data(type,M_,options_)
 % Checks the status of posterior analysis and in particular if files need to be
 % created or updated; called by posterior_analysis.m
 %
 % Inputs:
 %   type        [string]        name of the posterior moment considered
 %   M_          [structure]     Dynare model structure
+%   options_    [structure]     Dynare options structure
 %
 % Outputs:
 %   info        [scalar]        return code
@@ -17,7 +18,7 @@ function [info,description] = check_posterior_analysis_data(type,M_)
 %                                   info = 6; % Ok (nothing to do ;-)
 %   description [string]        Message corresponding to info
 
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,7 +40,17 @@ if nargout>1
     description = '';
 end
 
-[MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
+if ~issmc(options_)
+    [MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
+else
+    if ishssmc(options_)
+        [MetropolisFolder, info] = CheckPath('hssmc',M_.dname);
+    elseif isdime(options_)
+        [MetropolisFolder, info] = CheckPath('dime',M_.dname);
+    else
+        error('check_posterior_analysis_data:: case should not happen. Please contact the developers')
+    end
+end
 
 % Get informations about mcmc files.
 if info
@@ -47,9 +58,6 @@ if info
     return
 end
 
-mhname = get_name_of_the_last_mh_file(M_);
-mhdate = get_date_of_a_file([MetropolisFolder filesep mhname]);
-
 % Get informations about _posterior_draws files.
 drawsinfo = dir([ MetropolisFolder filesep M_.fname '_posterior_draws*.mat']);
 if isempty(drawsinfo)
@@ -60,6 +68,18 @@ if isempty(drawsinfo)
     return
 else
     number_of_last_posterior_draws_file = length(drawsinfo);
+    if ~issmc(options_)
+        mhname = get_name_of_the_last_mh_file(M_);
+        mhdate = get_date_of_a_file([MetropolisFolder filesep mhname]);
+    else
+        if ishssmc(options_)
+            % Load draws from the posterior distribution
+            pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+            mhdate = get_date_of_a_file(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+        elseif isdime(options_)
+            mhdate = get_date_of_a_file(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+        end
+    end
     pddate = get_date_of_a_file([ MetropolisFolder filesep M_.fname '_posterior_draws' int2str(number_of_last_posterior_draws_file) '.mat']);
     if pddate<mhdate
         info = 2; % _posterior_draws files have to be updated.
diff --git a/matlab/estimation/check_prior_analysis_data.m b/matlab/estimation/check_prior_analysis_data.m
index 569f324d22e6cf840889c56b7dc328c8b0b86130..ee400509870971f196387a15b016014c31d056f5 100644
--- a/matlab/estimation/check_prior_analysis_data.m
+++ b/matlab/estimation/check_prior_analysis_data.m
@@ -47,7 +47,6 @@ if ~exist([ M_.dname '/prior/draws'],'dir')
 end
 
 prior_draws_info = dir([ M_.dname '/prior/draws/prior_draws*.mat']);
-date_of_the_last_prior_draw_file = prior_draws_info(end).datenum;
 
 %% Get informations about _posterior_draws files.
 if isempty(prior_draws_info)
@@ -57,6 +56,7 @@ if isempty(prior_draws_info)
     end
     return
 else
+    date_of_the_last_prior_draw_file = prior_draws_info(end).datenum;
     date_of_the_prior_definition = get_date_of_a_file([ M_.dname '/prior/definition.mat']);
     if date_of_the_prior_definition>date_of_the_last_prior_draw_file
         info = 2;
diff --git a/matlab/estimation/conditional_variance_decomposition_ME_mc_analysis.m b/matlab/estimation/conditional_variance_decomposition_ME_mc_analysis.m
index 6c787b0c80551fd78f1eb2982d5c37d9a7918876..f09dff723fdc1ced3923aad5ed54b3fc9b0a2b80 100644
--- a/matlab/estimation/conditional_variance_decomposition_ME_mc_analysis.m
+++ b/matlab/estimation/conditional_variance_decomposition_ME_mc_analysis.m
@@ -23,7 +23,7 @@ function oo_ = ...
 % OUTPUTS
 %   oo_          [structure]        Dynare structure where the results are saved.
 
-% Copyright © 2017-2023 Dynare Team
+% Copyright © 2017-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -41,8 +41,9 @@ function oo_ = ...
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep ];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
diff --git a/matlab/estimation/conditional_variance_decomposition_mc_analysis.m b/matlab/estimation/conditional_variance_decomposition_mc_analysis.m
index ab4f7e77d88e6082b18af8cd1d545bfbaac52b91..b84a7b64cd9793f8a537d0d1204b36f22b598bf2 100644
--- a/matlab/estimation/conditional_variance_decomposition_mc_analysis.m
+++ b/matlab/estimation/conditional_variance_decomposition_mc_analysis.m
@@ -23,7 +23,7 @@ function oo_ = ...
 % OUTPUTS
 %   oo_          [structure]        Dynare structure where the results are saved.
 
-% Copyright © 2009-2018 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -41,8 +41,9 @@ function oo_ = ...
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
diff --git a/matlab/estimation/correlation_mc_analysis.m b/matlab/estimation/correlation_mc_analysis.m
index 94b10bb74a2626edacdfd48c76715092045eb1c6..d4b8235ec271ccf2f2d4d677ce41486946f2b394 100644
--- a/matlab/estimation/correlation_mc_analysis.m
+++ b/matlab/estimation/correlation_mc_analysis.m
@@ -2,7 +2,7 @@ function oo_ = correlation_mc_analysis(SampleSize,type,dname,fname,vartan,nvar,v
 % This function analyses the (posterior or prior) distribution of the
 % endogenous variables correlation function.
 
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -20,8 +20,9 @@ function oo_ = correlation_mc_analysis(SampleSize,type,dname,fname,vartan,nvar,v
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
diff --git a/matlab/estimation/covariance_mc_analysis.m b/matlab/estimation/covariance_mc_analysis.m
index 92b7b451c9b9c7811e8b6e717b21e3bace3b7ae3..d0c4f77f7077392c655e3b28915d17cbb1bc6452 100644
--- a/matlab/estimation/covariance_mc_analysis.m
+++ b/matlab/estimation/covariance_mc_analysis.m
@@ -19,7 +19,7 @@ function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,varta
 % OUTPUTS
 %   oo_                     [structure]        Dynare structure where the results are saved.
 
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,8 +37,9 @@ function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,varta
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
diff --git a/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m
index 6e052b6ec6c506af9d502c62a11188088b793920..ef55dad23f8fddd4f01bfc64b92c1f2cf066ef1c 100644
--- a/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m
+++ b/matlab/estimation/dsge_simulated_theoretical_conditional_variance_decomposition.m
@@ -19,7 +19,7 @@ function [nvar,vartan,NumberOfConditionalDecompFiles] = ...
 %   vartan                           [char]     array of characters (with nvar rows).
 %   NumberOfConditionalDecompFiles   [integer]  scalar, number of prior or posterior data files (for covariance).
 
-% Copyright © 2009-2023 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,9 +37,11 @@ function [nvar,vartan,NumberOfConditionalDecompFiles] = ...
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 
+folder_name=get_posterior_folder_name(options_);
+
 % Get informations about the _posterior_draws files.
 if strcmpi(type,'posterior')
-    NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]));
+    NumberOfDrawsFiles = length(dir([M_.dname filesep folder_name filesep M_.fname '_' type '_draws*' ]));
     posterior = 1;
 elseif strcmpi(type,'prior')
     NumberOfDrawsFiles = length(dir([M_.dname '/prior/draws/' type '_draws*' ]));
@@ -51,27 +53,16 @@ end
 
 %delete old stale files before creating new ones
 if posterior
-    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecomposition*'])
-    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecompositionME*'])
+    delete_stale_file([M_.dname filesep folder_name filesep M_.fname '_PosteriorConditionalVarianceDecomposition*'])
+    delete_stale_file([M_.dname filesep folder_name filesep M_.fname '_PosteriorConditionalVarianceDecompositionME*'])
 else
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecomposition*'])
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecompositionME*'])
 end
 
 % Set varlist (vartan)
-if ~posterior
-    if isfield(options_,'varlist')
-        temp = options_.varlist;
-    end
-    options_.varlist = options_.prior_analysis_endo_var_list;
-end
 endo_names=options_.varlist;
 [ivar,vartan ] = get_variables_list(options_, M_);
-if ~posterior
-    if exist('temp','var')
-        options_.varlist = temp;
-    end
-end
 nvar = length(ivar);
 
 % Set the size of the auto-correlation function to zero.
@@ -118,9 +109,12 @@ linea = 0;
 linea_ME = 0;
 for file = 1:NumberOfDrawsFiles
     if posterior
-        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
+        temp=load([M_.dname filesep folder_name filesep M_.fname '_' type '_draws' num2str(file) ]);
     else
         temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        if columns(temp.pdraws)==3 && isstruct(temp.pdraws{1,3})
+            temp.pdraws(:,2)=[];
+        end
     end
     isdrsaved = columns(temp.pdraws)-1;
     NumberOfDraws = rows(temp.pdraws);
@@ -146,7 +140,7 @@ for file = 1:NumberOfDrawsFiles
             ConditionalDecompFileNumber = ConditionalDecompFileNumber + 1;
             linea = 0;
             if posterior
-                save([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ...
+                save([M_.dname filesep folder_name filesep M_.fname '_PosteriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ...
                      'Conditional_decomposition_array','endo_names');
             else
                 save([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ...
@@ -167,7 +161,7 @@ for file = 1:NumberOfDrawsFiles
                 ConditionalDecompFileNumber_ME = ConditionalDecompFileNumber_ME + 1;
                 linea_ME = 0;
                 if posterior
-                    save([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecompME' int2str(ConditionalDecompFileNumber_ME) '.mat' ], ...
+                    save([M_.dname filesep folder_name filesep M_.fname '_PosteriorConditionalVarianceDecompME' int2str(ConditionalDecompFileNumber_ME) '.mat' ], ...
                          'Conditional_decomposition_array_ME','endo_names');
                 else
                     save([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecompME' int2str(ConditionalDecompFileNumber_ME) '.mat' ], ...
diff --git a/matlab/estimation/dsge_simulated_theoretical_correlation.m b/matlab/estimation/dsge_simulated_theoretical_correlation.m
index d874cb54ba08fc185ebb176db849ec2b86259604..7fa195a2cdc4c1e850df072df660cae0bbca667e 100644
--- a/matlab/estimation/dsge_simulated_theoretical_correlation.m
+++ b/matlab/estimation/dsge_simulated_theoretical_correlation.m
@@ -18,7 +18,7 @@ function [nvar,vartan,CorrFileNumber] = dsge_simulated_theoretical_correlation(S
 %   vartan         [char]           array of characters (with nvar rows).
 %   CorrFileNumber [integer]        scalar, number of prior or posterior data files (for correlation).
 
-% Copyright © 2007-2021 Dynare Team
+% Copyright © 2007-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -38,9 +38,11 @@ function [nvar,vartan,CorrFileNumber] = dsge_simulated_theoretical_correlation(S
 [ivar,vartan] = get_variables_list(options_, M_);
 nvar = length(ivar);
 
+folder_name=get_posterior_folder_name(options_);
+
 % Get informations about the _posterior_draws files.
 if strcmpi(type,'posterior')
-    CorrFileNumber = length(dir([M_.dname '/metropolis/' M_.fname '_PosteriorCorrelations*']));
+    CorrFileNumber = length(dir([M_.dname filesep folder_name filesep M_.fname '_PosteriorCorrelations*']));
 elseif strcmpi(type,'prior')
     CorrFileNumber = length(dir([M_.dname '/prior/moments/' M_.fname '_PriorCorrelations*']));
 else
diff --git a/matlab/estimation/dsge_simulated_theoretical_covariance.m b/matlab/estimation/dsge_simulated_theoretical_covariance.m
index 31d5a1486d1bde18d5a08919e70948af9db42303..267f0d0e4d1ffebbc095736c703ea3aeb35e2265 100644
--- a/matlab/estimation/dsge_simulated_theoretical_covariance.m
+++ b/matlab/estimation/dsge_simulated_theoretical_covariance.m
@@ -16,7 +16,7 @@ function [nvar,vartan,CovarFileNumber] = dsge_simulated_theoretical_covariance(S
 %   vartan            [char]     array of characters (with nvar rows).
 %   CovarFileNumber   [integer]  scalar, number of prior or posterior data files (for covariance).
 
-% Copyright © 2007-2021 Dynare Team
+% Copyright © 2007-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -35,9 +35,11 @@ function [nvar,vartan,CovarFileNumber] = dsge_simulated_theoretical_covariance(S
 
 nodecomposition = 1;
 
+folder_name=get_posterior_folder_name(options_);
+    
 % Get informations about the _posterior_draws files.
 if strcmpi(type,'posterior')
-    NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]));
+    NumberOfDrawsFiles = length(dir([M_.dname filesep folder_name filesep M_.fname '_' type '_draws*' ]));
     posterior = 1;
 elseif strcmpi(type,'prior')
     NumberOfDrawsFiles = length(dir([M_.dname '/prior/draws/' type '_draws*' ]));
@@ -49,27 +51,16 @@ end
 
 %delete old stale files before creating new ones
 if posterior
-    delete_stale_file([M_.dname '/metropolis/' M_.fname '_Posterior2ndOrderMoments*'])
-    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorCorrelations*']);
+    delete_stale_file([M_.dname filesep folder_name filesep M_.fname '_Posterior2ndOrderMoments*'])
+    delete_stale_file([M_.dname filesep folder_name filesep M_.fname '_PosteriorCorrelations*']);
 else
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_Prior2ndOrderMoments*'])
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorCorrelations*']);
 end
 
 % Set varlist (vartan)
-if ~posterior
-    if isfield(options_,'varlist')
-        temp = options_.varlist;
-    end
-    options_.varlist = options_.prior_analysis_endo_var_list;
-end
 endo_names=options_.varlist;
 [ivar,vartan] = get_variables_list(options_,M_);
-if ~posterior
-    if exist('temp','var')
-        options_.varlist = temp;
-    end
-end
 nvar = length(ivar);
 
 if options_.pruning
@@ -112,9 +103,12 @@ linea_cov = 0;
 linea_corr = 0;
 for file = 1:NumberOfDrawsFiles
     if posterior
-        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
+        temp=load([M_.dname filesep folder_name filesep M_.fname '_' type '_draws' num2str(file) ]);
     else
         temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        if columns(temp.pdraws)==3 && isstruct(temp.pdraws{1,3})
+            temp.pdraws(:,2)=[];
+        end
     end
     NumberOfDraws = rows(temp.pdraws);
     isdrsaved = columns(temp.pdraws)-1;
@@ -148,7 +142,7 @@ for file = 1:NumberOfDrawsFiles
 
         if linea_cov == NumberOfCovarLines
             if posterior
-                save([ M_.dname '/metropolis/' M_.fname '_Posterior2ndOrderMoments' int2str(CovarFileNumber) '.mat' ],'Covariance_matrix','endo_names');
+                save([ M_.dname filesep folder_name filesep M_.fname '_Posterior2ndOrderMoments' int2str(CovarFileNumber) '.mat' ],'Covariance_matrix','endo_names');
             else
                 save([ M_.dname '/prior/moments/' M_.fname '_Prior2ndOrderMoments' int2str(CovarFileNumber) '.mat' ],'Covariance_matrix','endo_names');
             end
@@ -166,7 +160,7 @@ for file = 1:NumberOfDrawsFiles
         end
         if linea_corr == NumberOfCorrLines
             if posterior
-                save([ M_.dname '/metropolis/' M_.fname '_PosteriorCorrelations' int2str(CorrFileNumber) '.mat' ],'Correlation_array','endo_names');
+                save([ M_.dname filesep folder_name filesep M_.fname '_PosteriorCorrelations' int2str(CorrFileNumber) '.mat' ],'Correlation_array','endo_names');
             else
                 save([ M_.dname '/prior/moments/' M_.fname '_PriorCorrelations' int2str(CorrFileNumber) '.mat' ],'Correlation_array','endo_names');
             end
diff --git a/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m b/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m
index 880b27438fd24c358c264befc21dd5e356e41381..be017702ff941d73e7af247ff5434fe1477bdfdb 100644
--- a/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m
+++ b/matlab/estimation/dsge_simulated_theoretical_variance_decomposition.m
@@ -18,7 +18,7 @@ function [nvar,vartan,NumberOfDecompFiles] = ...
 %   vartan            [char]     array of characters (with nvar rows).
 %   CovarFileNumber   [integer]  scalar, number of prior or posterior data files (for covariance).
 
-% Copyright © 2007-2021 Dynare Team
+% Copyright © 2007-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,9 +37,11 @@ function [nvar,vartan,NumberOfDecompFiles] = ...
 
 nodecomposition = 0;
 
+folder_name=get_posterior_folder_name(options_);
+
 % Get informations about the _posterior_draws files.
 if strcmpi(type,'posterior')
-    NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]));
+    NumberOfDrawsFiles = length(dir([M_.dname filesep folder_name filesep M_.fname '_' type '_draws*' ]));
     posterior = 1;
 elseif strcmpi(type,'prior')
     NumberOfDrawsFiles = length(dir([M_.dname '/prior/draws/' type '_draws*' ]));
@@ -51,27 +53,16 @@ end
 
 %delete old stale files before creating new ones
 if posterior
-    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecomposition*']);
-    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecompME*']);
+    delete_stale_file([M_.dname filesep folder_name filesep M_.fname '_PosteriorVarianceDecomposition*']);
+    delete_stale_file([M_.dname filesep folder_name filesep M_.fname '_PosteriorVarianceDecompME*']);
 else
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecomposition*']);
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME*']);
 end
 
 % Set varlist (vartan)
-if ~posterior
-    if isfield(options_,'varlist')
-        temp = options_.varlist;
-    end
-    options_.varlist = options_.prior_analysis_endo_var_list;
-end
 [ivar,vartan,options_] = get_variables_list(options_,M_);
 endo_names=options_.varlist;
-if ~posterior
-    if exist('temp','var')
-        options_.varlist = temp;
-    end
-end
 nvar = length(ivar);
 
 % Set the size of the auto-correlation function to zero.
@@ -122,9 +113,12 @@ linea_ME = 0;
 only_non_stationary_vars=0;
 for file = 1:NumberOfDrawsFiles
     if posterior
-        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
+        temp=load([M_.dname filesep folder_name filesep M_.fname '_' type '_draws' num2str(file) ]);
     else
         temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        if columns(temp.pdraws)==3 && isstruct(temp.pdraws{1,3})
+            temp.pdraws(:,2)=[];
+        end
     end
     isdrsaved = columns(temp.pdraws)-1;
     NumberOfDraws = rows(temp.pdraws);
@@ -170,7 +164,7 @@ for file = 1:NumberOfDrawsFiles
         end
         if linea == NumberOfDecompLines
             if posterior
-                save([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecomposition' int2str(DecompFileNumber) '.mat' ],'Decomposition_array','endo_names');
+                save([M_.dname filesep folder_name filesep M_.fname '_PosteriorVarianceDecomposition' int2str(DecompFileNumber) '.mat' ],'Decomposition_array','endo_names');
             else
                 save([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecomposition' int2str(DecompFileNumber) '.mat' ],'Decomposition_array','endo_names');
             end
@@ -189,7 +183,7 @@ for file = 1:NumberOfDrawsFiles
         if ME_present
             if linea_ME == NumberOfDecompLines_ME
                 if posterior
-                    save([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecompME' int2str(DecompFileNumber_ME) '.mat' ],'Decomposition_array_ME','endo_names');
+                    save([M_.dname filesep folder_name filesep M_.fname '_PosteriorVarianceDecompME' int2str(DecompFileNumber_ME) '.mat' ],'Decomposition_array_ME','endo_names');
                 else
                     save([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME' int2str(DecompFileNumber_ME) '.mat' ],'Decomposition_array_ME','endo_names');
                 end
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index 8d1efefbc4b2c4cc48206352a3923d96b2522f4e..a2e9ffd8c100a8ce8a4180369dde10be3eeba310 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -483,9 +483,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
             end
             % Store posterior mean in a vector and posterior variance in
             % a matrix
-            if ~isdime(options_)
-                [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] = GetPosteriorMeanVariance(options_, M_);
-            end
+            [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] = GetPosteriorMeanVariance(options_, M_);
         elseif options_.load_mh_file && options_.load_results_after_load_mh
             % load fields from previous MCMC run stored in results-file
             field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
@@ -506,28 +504,22 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
                 oo_.posterior.metropolis=oo_load_mh.oo_.posterior.metropolis;
             end
         end
-        if ~issmc(options_)
-            [error_flag, ~, options_]= metropolis_draw(1, options_, estim_params_, M_);
-        else
-            error_flag=false;
-        end
+        [options_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_); %check whether number is feasible
         if ~(~isempty(options_.sub_draws) && options_.sub_draws==0)
             if options_.bayesian_irf
-                if ~issmc(options_)
-                    if error_flag
-                        error('%s: I cannot compute the posterior IRFs!',dispString)
-                    end
-                    oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
-                else
-                    fprintf('%s: SMC does not yet support the bayesian_irf option. Skipping computation.\n',dispString);
+                if error_flag
+                    error('%s: I cannot compute the posterior IRFs!',dispString)
                 end
+                oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
             end
             if options_.moments_varendo
-                if ~issmc(options_)
-                    if error_flag
-                        error('%s: I cannot compute the posterior moments for the endogenous variables!',dispString)
-                    end
-                    if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results
+                if error_flag
+                    error('%s: I cannot compute the posterior moments for the endogenous variables!',dispString)
+                end
+                if options_.load_mh_file
+                    if issmc(options_)
+                        error('%s: SMC does not yet support the load_mh_file option.\n',dispString);
+                    elseif options_.mh_replic==0 %user wants to recompute results for standard MCMC
                         [MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
                         if ~info
                             generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'};
@@ -546,23 +538,17 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
                             end
                         end
                     end
-                    oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
-                else
-                    fprintf('%s: SMC does not yet support the moments_varendo option. Skipping computation.\n',dispString);
                 end
+                oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
             end
             if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
-                if ~ishssmc(options_) && ~isdime(options_)
-                    if error_flag
-                        error('%s: I cannot compute the posterior statistics!',dispString)
-                    end
-                    if options_.order==1 && ~options_.particle.status
-                        oo_=prior_posterior_statistics('posterior',dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString); %get smoothed and filtered objects and forecasts
-                    else
-                        error('%s: Particle Smoothers are not yet implemented.',dispString)
-                    end
+                if error_flag
+                    error('%s: I cannot compute the posterior statistics!',dispString)
+                end
+                if options_.order==1 && ~options_.particle.status
+                    oo_=prior_posterior_statistics('posterior',dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString); %get smoothed and filtered objects and forecasts
                 else
-                    fprintf('%s: SMC does not yet support the smoother and forecast options. Skipping computation.\n',dispString);
+                    error('%s: Particle Smoothers are not yet implemented.',dispString)
                 end
             end
         else
diff --git a/matlab/estimation/execute_prior_posterior_function.m b/matlab/estimation/execute_prior_posterior_function.m
index f23d4664babdb2bb206837bc6671e1a6abeb5253..898cf63005952007353f40404ec49f83ee749613 100644
--- a/matlab/estimation/execute_prior_posterior_function.m
+++ b/matlab/estimation/execute_prior_posterior_function.m
@@ -17,7 +17,7 @@ function oo_=execute_prior_posterior_function(posterior_function_name,M_,options
 % OUTPUTS
 %   oo_          [structure]     Matlab/Octave structure gathering the results (initialized by dynare, see @ref{oo_}).
 
-% Copyright © 2013-2023 Dynare Team
+% Copyright © 2013-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -53,9 +53,8 @@ if strcmpi(type,'posterior')
     % Get informations about the _posterior_draws files.
     % discard first mh_drop percent of the draws:
     CutSample(M_, options_, 'prior_posterior_function');
-    % initialize metropolis draws
     options_.sub_draws = n_draws; % set draws for sampling; changed value is not returned to base workspace
-    [error_flag, ~, options_] = metropolis_draw(1, options_, estim_params_, M_);
+    [options_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_); %check whether number is feasible
     if error_flag
         error('EXECUTE_POSTERIOR_FUNCTION: The draws could not be initialized')
     end
@@ -69,7 +68,7 @@ elseif strcmpi(type,'prior')
             error('The prior distributions are not properly set up.')
         end
     end
-    if exist([M_.fname '_prior_restrictions.m'])
+    if exist([M_.fname '_prior_restrictions.m'],"file")
         warning('prior_function currently does not support endogenous prior restrictions. They will be ignored. Consider using a posterior_function with nobs=1.')
     end
     Prior = dprior(bayestopt_, options_.prior_trunc);
@@ -77,13 +76,11 @@ else
     error('EXECUTE_POSTERIOR_FUNCTION: Unknown type!')
 end
 
+
 if strcmpi(type, 'prior')
     parameter_mat = Prior.draws(n_draws);
 else
-    parameter_mat = NaN(length(bayestopt_.p6), n_draws);
-    for i = 1:n_draws
-        parameter_mat(:,i) = GetOneDraw(type, M_, estim_params_, oo_, options_, bayestopt_);
-    end
+    parameter_mat=get_posterior_subsample(M_,options_,n_draws)';
 end
 
 % Get output size
diff --git a/matlab/estimation/get_posterior_folder_name.m b/matlab/estimation/get_posterior_folder_name.m
new file mode 100644
index 0000000000000000000000000000000000000000..25b86f1c8b6375b691dfff125cbee03b0752aafd
--- /dev/null
+++ b/matlab/estimation/get_posterior_folder_name.m
@@ -0,0 +1,42 @@
+function folder_name=get_posterior_folder_name(options_)
+% function folder_name=get_posterior_folder_name(options_)
+% runs the estimation of the model
+%
+% INPUTS
+%   options_            [structure] describing the options
+%
+% OUTPUTS
+%   folder_name         [char]      output folder for different posterior
+%                                   samplers
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+if ~issmc(options_)
+    folder_name='metropolis';
+else
+    if ishssmc(options_)
+        folder_name='hssmc';
+    elseif isdime(options_)
+        folder_name='dime';
+    else
+        error('get_posterior_folder_name:: case should not happen. Please contact the developers')
+    end
+end
\ No newline at end of file
diff --git a/matlab/estimation/get_posterior_subsample.m b/matlab/estimation/get_posterior_subsample.m
new file mode 100644
index 0000000000000000000000000000000000000000..5df136f5a2f3f9bc71891037b91f00cadd02b58f
--- /dev/null
+++ b/matlab/estimation/get_posterior_subsample.m
@@ -0,0 +1,105 @@
+function [pdraws, log_posterior]=get_posterior_subsample(M_,options_,SampleSize)
+% function [pdraws, log_posterior]=get_posterior_subsample(M_,options_,SampleSize)
+% Checks the status of posterior analysis and in particular if files need to be
+% created or updated; called by posterior_analysis.m
+%
+% Inputs:
+%   M_              [structure]     Dynare model structure
+%   options_        [structure]     Dynare options structure
+%   SampleSize      [integer]       number of desired subdraws
+%
+% Outputs:
+%   pdraws          [double]        SampleSize by npar matrix of parameter draws
+%   log_posterior   [double]        SampleSize by 1 vector of log posterior
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+if ~issmc(options_)
+    % Get informations about the mcmc:
+    record=load_last_mh_history_file([M_.dname filesep 'metropolis'], M_.fname);
+    npar=size(record.InitialParameters,2);
+    FirstMhFile = record.KeepedDraws.FirstMhFile;
+    FirstLine = record.KeepedDraws.FirstLine;
+    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+    NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+    MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
+
+    % Randomly select draws in the posterior distribution:
+    SampleAddress = zeros(SampleSize,4);
+    for iter = 1:SampleSize
+        ChainNumber = ceil(rand*options_.mh_nblck);
+        DrawNumber  = ceil(rand*NumberOfDraws);
+        SampleAddress(iter,1) = DrawNumber;
+        SampleAddress(iter,2) = ChainNumber;
+        if DrawNumber <= MAX_nruns-FirstLine+1 % in first file where FirstLine may be bigger than 1
+            MhFileNumber = FirstMhFile;
+            MhLineNumber = FirstLine+DrawNumber-1;
+        else
+            DrawNumber  = DrawNumber-(MAX_nruns-FirstLine+1);
+            MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
+            MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
+        end
+        SampleAddress(iter,3) = MhFileNumber;
+        SampleAddress(iter,4) = MhLineNumber;
+    end
+    SampleAddress = sortrows(SampleAddress,[3 2]);
+
+    % Selected draws in the posterior distribution, and if drsize>0
+    % reduced form solutions, are saved on disk.
+    pdraws = NaN(SampleSize,npar);
+    if nargout>1
+        log_posterior = NaN(SampleSize,1);
+    end
+    old_mhfile = 0;
+    old_mhblck = 0;
+    for iter = 1:SampleSize
+        mhfile = SampleAddress(iter,3);
+        mhblck = SampleAddress(iter,2);
+        if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
+            temp=load([M_.dname filesep 'metropolis' filesep M_.fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2','logpo2');
+        end
+        pdraws(iter,:) = temp.x2(SampleAddress(iter,4),:);
+        if nargout>1
+            log_posterior(iter,1) = temp.logpo2(SampleAddress(iter,4));
+        end
+        old_mhfile = mhfile;
+        old_mhblck = mhblck;
+    end
+else
+    if ishssmc(options_)
+        % Load draws from the posterior distribution
+        pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+        posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+        NumberOfDraws=size(posterior.particles,2);
+        DrawNumbers  = ceil(rand(SampleSize,1)*NumberOfDraws);
+        pdraws = transpose(posterior.particles(:,DrawNumbers));
+        if nargout>1
+            log_posterior = posterior.tlogpostkernel(DrawNumbers);
+        end
+    elseif isdime(options_)
+        posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+        tune = posterior.tune;
+        chains = reshape(posterior.chains(end-tune:end,:,:), [], size(posterior.chains, 3));
+        NumberOfDraws=size(chains,1);
+        DrawNumbers  = ceil(rand(SampleSize,1)*NumberOfDraws);
+        pdraws = chains(DrawNumbers,:);
+        if nargout>1
+            log_posterior = posterior.lprobs(DrawNumbers);
+        end
+    end
+end
\ No newline at end of file
diff --git a/matlab/estimation/metropolis_draw.m b/matlab/estimation/metropolis_draw.m
deleted file mode 100644
index 505d1dbeb6b5deaa7964a366b882f5fabf2cd2a7..0000000000000000000000000000000000000000
--- a/matlab/estimation/metropolis_draw.m
+++ /dev/null
@@ -1,99 +0,0 @@
-function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
-% function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
-% Builds draws from metropolis
-%
-% INPUTS:
-%   init:              scalar equal to
-%                      1: first call to store the required information
-%                           on files, lines, and chains to be read
-%                           in persistent variables to make them available
-%                           for future calls
-%                      0: load a parameter draw
-% Additional Inputs required for initialization
-%   options_           [structure]     Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
-%   M_                 [structure]     Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
-%   estim_params_      [structure]     Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
-%
-% OUTPUTS:
-%   xparams:           if init==0: vector of estimated parameters
-%                      if init==1: error flaog
-%   logpost:           log of posterior density
-%   options_:          [structure]     Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
-%
-% SPECIAL REQUIREMENTS
-%
-%   Requires CutSample to be run before in order to set up mh_history-file
-
-% Copyright © 2003-2023 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns
-
-xparams = 0;
-logpost = 0;
-
-if init
-    %get number of parameters
-    nvx  = estim_params_.nvx;
-    nvn  = estim_params_.nvn;
-    ncx  = estim_params_.ncx;
-    ncn  = estim_params_.ncn;
-    np   = estim_params_.np ;
-    npar = nvx+nvn+ncx+ncn+np;
-    %get path of metropolis files
-    MetropolisFolder = CheckPath('metropolis',M_.dname);
-    FileName = M_.fname;
-    BaseName = [MetropolisFolder filesep FileName];
-    %load mh_history-file with info on what to load
-    record=load_last_mh_history_file(MetropolisFolder, FileName);
-    FirstMhFile = record.KeepedDraws.FirstMhFile;
-    FirstLine = record.KeepedDraws.FirstLine;
-    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 ?
-    mh_nblck = options_.mh_nblck;
-    % set sub_draws option if empty
-    if isempty(options_.sub_draws)
-        options_.sub_draws = min(options_.posterior_max_subsample_draws, ceil(.25*NumberOfDraws));
-    else
-        if options_.sub_draws>NumberOfDraws*mh_nblck
-            skipline()
-            disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*mh_nblck) ')!'])
-            disp('You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).')
-            skipline()
-            xparams = 1; % xparams is interpreted as an error flag
-        end
-    end
-    return
-else %not initialization, return one draw
-     %get random draw from random chain
-    ChainNumber = ceil(rand*mh_nblck);
-    DrawNumber  = ceil(rand*NumberOfDraws);
-
-    if DrawNumber <= MAX_nruns-FirstLine+1 %draw in first file, needs to account for first line
-        MhFilNumber = FirstMhFile;
-        MhLine = FirstLine+DrawNumber-1;
-    else %draw in other file
-        DrawNumber  = DrawNumber-(MAX_nruns-FirstLine+1);
-        MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
-        MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
-    end
-    %load parameters and posterior
-    load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2');
-    xparams = x2(MhLine,:);
-    logpost= logpo2(MhLine);
-end
\ No newline at end of file
diff --git a/matlab/estimation/mh_autocorrelation_function.m b/matlab/estimation/mh_autocorrelation_function.m
index 953b542efae722a1fd5db4f8fa27e7b853111277..04ef6fb2c0c848a25dfbdd10b2b2a62f9c6871cb 100644
--- a/matlab/estimation/mh_autocorrelation_function.m
+++ b/matlab/estimation/mh_autocorrelation_function.m
@@ -45,21 +45,8 @@ if isempty(column)
     return
 end
 
-% Get informations about the posterior draws:
-MetropolisFolder = CheckPath('metropolis',M_.dname);
-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(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck, blck);
+PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, [],[], blck);
 
 % Compute the autocorrelation function:
 [~,autocor] = sample_autocovariance(PosteriorDraws,options_.mh_autocorrelation_function_size);
diff --git a/matlab/estimation/posterior_analysis.m b/matlab/estimation/posterior_analysis.m
index 72f4713da3bdd09b24a987195124ba3f0566a261..b334bdce9ae6dcb7a281b525180c37987d335f74 100644
--- a/matlab/estimation/posterior_analysis.m
+++ b/matlab/estimation/posterior_analysis.m
@@ -17,28 +17,28 @@ function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_,estim_para
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-info = check_posterior_analysis_data(type,M_);
+info = check_posterior_analysis_data(type,M_,options_);
 SampleSize = options_.sub_draws;
 switch info
-  case 0
-    disp('check_posterior_analysis_data:: Can''t find any mcmc file!')
-    error('Check the options of the estimation command...')
-  case {1,2}
-    MaxMegaBytes = options_.MaximumNumberOfMegaBytes;
-    drsize = size_of_the_reduced_form_model(oo_.dr);
-    if drsize*SampleSize>MaxMegaBytes
-        drsize=0;
-    end
-    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_);
-  case 6
-    [ivar,vartan] = get_variables_list(options_,M_);
-    nvar = length(ivar);
-    oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
-  otherwise
-    error('posterior_analysis:: Check_posterior_analysis_data gave a meaningless output!')
+    case 0
+        disp('check_posterior_analysis_data:: Can''t find any mcmc file!')
+        error('Check the options of the estimation command...')
+    case {1,2} %need to get draws and store posterior
+        MaxMegaBytes = options_.MaximumNumberOfMegaBytes;
+        drsize = size_of_the_reduced_form_model(oo_.dr);
+        if drsize*SampleSize>MaxMegaBytes
+            drsize=0;
+        end
+        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} %process draws and save files
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
+    case 6 %just read out
+        [ivar,vartan] = get_variables_list(options_,M_);
+        nvar = length(ivar);
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
+    otherwise
+        error('posterior_analysis:: Check_posterior_analysis_data gave a meaningless output!')
 end
 
 
diff --git a/matlab/estimation/prior_analysis.m b/matlab/estimation/prior_analysis.m
index a247ba6f1caf908c92456eb89d258deb594b2651..e9065a4516d3e66519bd66240cdd49fba889abe4 100644
--- a/matlab/estimation/prior_analysis.m
+++ b/matlab/estimation/prior_analysis.m
@@ -19,26 +19,26 @@ function oo_ = prior_analysis(type,arg1,arg2,arg3,options_,M_,oo_,estim_params_)
 info = check_prior_analysis_data(type,M_);
 SampleSize = options_.prior_mc;
 switch info
-  case {0,1,2}
-    MaxMegaBytes = options_.MaximumNumberOfMegaBytes;
-    drsize = size_of_the_reduced_form_model(oo_.dr);
-    if drsize*SampleSize>MaxMegaBytes
-        drsave=0;
-    else
-        drsave=1;
-    end
-    load([M_.dname '/prior/definition.mat']);
-    prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_);
-    clear('bayestopt_');
-    oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
-  case {4,5}
-    oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
-  case 6
-    [ivar,vartan] = get_variables_list(options_,M_);
-    nvar = length(ivar);
-    oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
-  otherwise
-    error('prior_analysis:: Check_prior_analysis_data gave a meaningless output!')
+    case {0,1,2}
+        MaxMegaBytes = options_.MaximumNumberOfMegaBytes;
+        drsize = size_of_the_reduced_form_model(oo_.dr);
+        if drsize*SampleSize>MaxMegaBytes
+            drsave=0;
+        else
+            drsave=1;
+        end
+        load([M_.dname '/prior/definition.mat']);
+        prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_);
+        clear('bayestopt_');
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
+    case {4,5}
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
+    case 6
+        [ivar,vartan] = get_variables_list(options_,M_);
+        nvar = length(ivar);
+        oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
+    otherwise
+        error('prior_analysis:: Check_prior_analysis_data gave a meaningless output!')
 end
 
 
@@ -53,7 +53,7 @@ switch type
   case 'variance'
     if nargin==narg1
         [nvar,vartan] = ...
-            dsge_simulated_theoretical_covariance(SampleSize,M_,options_,oo_,'prior');
+            dsge_simulated_theoretical_covariance(SampleSize,arg3,M_,options_,oo_,'prior');
     end
     oo_ = covariance_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                  vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_,options_);
@@ -64,6 +64,17 @@ switch type
     end
     oo_ = variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                              M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+    if ~all(diag(M_.H)==0)
+        if strmatch(arg1,options_.varobs,'exact')
+            if isoctave && octave_ver_less_than('8.4') %Octave bug #60347
+                observable_name_requested_vars=intersect_stable(vartan,options_.varobs);
+            else
+                observable_name_requested_vars=intersect(vartan,options_.varobs,'stable');
+            end
+            oo_ = variance_decomposition_ME_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
+                [M_.exo_names;'ME'],arg2,observable_name_requested_vars,arg1,options_.mh_conf_sig,oo_,options_);
+        end
+    end
   case 'correlation'
     if nargin==narg1
         [nvar,vartan] = ...
@@ -79,9 +90,9 @@ switch type
     oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                                       arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
     if ~all(diag(M_.H)==0)
-        if strmatch(vartan(arg1,:),options_.varobs,'exact')
+        if strmatch(arg1,options_.varobs,'exact')
             oo_ = conditional_variance_decomposition_ME_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
-                                                              arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+                                                              arg3,[M_.exo_names;'ME'],arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
         end
     end
   otherwise
diff --git a/matlab/estimation/prior_posterior_statistics.m b/matlab/estimation/prior_posterior_statistics.m
index 6971f9465660c4f7e6b624b54786b382638b83bd..0768a4d30f653d033f4513c82da60a77ccbf12f4 100644
--- a/matlab/estimation/prior_posterior_statistics.m
+++ b/matlab/estimation/prior_posterior_statistics.m
@@ -22,7 +22,7 @@ function oo_=prior_posterior_statistics(type,dataset_,dataset_info,M_,oo_,option
 % See the comments in the posterior_sampler.m funtion.
 
 
-% Copyright © 2005-2023 Dynare Team
+% Copyright © 2005-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -62,7 +62,8 @@ end
 maxlag = M_.maximum_endo_lag;
 
 if strcmpi(type,'posterior')
-    DirectoryName = CheckPath('metropolis',M_.dname);
+    folder_name=get_posterior_folder_name(options_);
+    DirectoryName = CheckPath(folder_name,M_.dname);
     B = options_.sub_draws;
 elseif strcmpi(type,'gsa')
     RootDirectoryName = CheckPath('gsa',M_.dname);
@@ -202,25 +203,14 @@ localVars.bayestopt_=bayestopt_;
 
 
 if strcmpi(type,'posterior')
-    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));
-    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-    NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
-    mh_nblck = options_.mh_nblck;
-    if B==NumberOfDraws*mh_nblck
+    [~, ~, NumberOfDraws]=set_number_of_subdraws(M_,options_);
+
+    if B==NumberOfDraws
         % we load all retained MH runs !
-        logpost=GetAllPosteriorDraws(options_, M_.dname, M_.fname, 0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
-        for column=1:npar
-            x(:,column) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
-        end
+        logpost=GetAllPosteriorDraws(options_, M_.dname, M_.fname, 0);
+        x = GetAllPosteriorDraws(options_, M_.dname, M_.fname, 'all');
     else
-        logpost=NaN(B,1);
-        for b=1:B
-            [x(b,:), logpost(b)] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
-        end
+        [x, logpost]=get_posterior_subsample(M_,options_,B);
     end
     localVars.logpost=logpost;
 end
diff --git a/matlab/estimation/prior_posterior_statistics_core.m b/matlab/estimation/prior_posterior_statistics_core.m
index 65660dabb59324575ac2bbcee33e3d90c0aff31b..5cf8f364e5ff56e0b1de55f0842b981effb04afa 100644
--- a/matlab/estimation/prior_posterior_statistics_core.m
+++ b/matlab/estimation/prior_posterior_statistics_core.m
@@ -119,7 +119,8 @@ end
 
 % DirectoryName = myinputs.DirectoryName;
 if strcmpi(type,'posterior')
-    DirectoryName = CheckPath('metropolis',M_.dname);
+    folder_name=get_posterior_folder_name(options_);
+    DirectoryName = CheckPath(folder_name,M_.dname);
 elseif strcmpi(type,'gsa')
     if options_.opt_gsa.pprior
         DirectoryName = CheckPath(['gsa',filesep,'prior'],M_.dname);
diff --git a/matlab/estimation/prior_sampler.m b/matlab/estimation/prior_sampler.m
index 172190ce716fbc363cfb2c4a0b7b27304d7f2130..3521df61f3befa36cd0d34382244e1a562a10198 100644
--- a/matlab/estimation/prior_sampler.m
+++ b/matlab/estimation/prior_sampler.m
@@ -6,6 +6,8 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 %   M_          [structure]  Model description.
 %   bayestopt_  [structure]  Prior distribution description.
 %   options_    [structure]  Global options of Dynare.
+%   oo_         [structure]     Dynare structure where the results are saved.
+%   estim_params_       [structure] characterizing parameters to be estimated
 %
 % OUTPUTS:
 %   results     [structure]  Various statistics.
@@ -13,7 +15,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2009-2023 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,14 +35,11 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
 % Initialization.
 Prior = dprior(bayestopt_, options_.prior_trunc);
 PriorDirectoryName = CheckPath('prior/draws',M_.dname);
-work = ~drsave;
 iteration = 0;
 loop_indx = 0;
-file_indx = [];
 count_bk_indeterminacy = 0;
 count_bk_unstability = 0;
 count_bk_singularity = 0;
-count_static_var_def = 0;
 count_no_steadystate = 0;
 count_steadystate_file_exit = 0;
 count_dll_problem = 0;
diff --git a/matlab/estimation/selec_posterior_draws.m b/matlab/estimation/selec_posterior_draws.m
index 624c6e54f6e648553ec6d882860faba041cc983e..161f8e8e39f284442928c94723c02036e7ebd6c6 100644
--- a/matlab/estimation/selec_posterior_draws.m
+++ b/matlab/estimation/selec_posterior_draws.m
@@ -2,8 +2,7 @@ function SampleAddress = selec_posterior_draws(M_,options_,dr,endo_steady_state,
 % 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.
-% This routine is more efficient than metropolis_draw.m because here an
-% _mh file cannot be opened twice.
+% This routine assures an _mh file cannot be opened twice.
 %
 % INPUTS
 %   o M_                    [structure]     Matlab's structure describing the model
@@ -55,7 +54,7 @@ switch nargin
   case 8
     info = 0;
   case 9
-    MAX_mega_bytes = 10;% Should be an option...
+    MAX_mega_bytes = 100;% Should be an option...
     if drsize>0
         info=2;
     else
@@ -66,100 +65,55 @@ switch nargin
     error('selec_posterior_draws:: Unexpected number of input arguments!')
 end
 
-MetropolisFolder = CheckPath('metropolis',M_.dname);
-ModelName = M_.fname;
-BaseName = [MetropolisFolder filesep ModelName];
-
-% Get informations about the mcmc:
-record=load_last_mh_history_file(MetropolisFolder, ModelName);
-FirstMhFile = record.KeepedDraws.FirstMhFile;
-FirstLine = record.KeepedDraws.FirstLine;
-TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
-MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
-mh_nblck = options_.mh_nblck;
-
-% Randomly select draws in the posterior distribution:
-SampleAddress = zeros(SampleSize,4);
-for i = 1:SampleSize
-    ChainNumber = ceil(rand*mh_nblck);
-    DrawNumber  = ceil(rand*NumberOfDraws);
-    SampleAddress(i,1) = DrawNumber;
-    SampleAddress(i,2) = ChainNumber;
-    if DrawNumber <= MAX_nruns-FirstLine+1
-        MhFileNumber = FirstMhFile;
-        MhLineNumber = FirstLine+DrawNumber-1;
+if ~issmc(options_)
+    MetropolisFolder = CheckPath('metropolis',M_.dname);
+else
+    if ishssmc(options_)
+        [MetropolisFolder] = CheckPath('hssmc',M_.dname);
+    elseif isdime(options_)
+        [MetropolisFolder] = CheckPath('dime',M_.dname);
     else
-        DrawNumber  = DrawNumber-(MAX_nruns-FirstLine+1);
-        MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
-        MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
+        error('check_posterior_analysis_data:: case should not happen. Please contact the developers')
     end
-    SampleAddress(i,3) = MhFileNumber;
-    SampleAddress(i,4) = MhLineNumber;
 end
-SampleAddress = sortrows(SampleAddress,[3 2]);
+
+ModelName = M_.fname;
+BaseName = [MetropolisFolder filesep ModelName];
+
+parameter_draws=get_posterior_subsample(M_,options_,SampleSize);
 
 % Selected draws in the posterior distribution, and if drsize>0
 % reduced form solutions, are saved on disk.
 if info
     %delete old stale files before creating new ones
     delete_stale_file([BaseName '_posterior_draws*.mat'])
-    if  SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file.
-        pdraws = cell(SampleSize,info);
-        old_mhfile = 0;
-        old_mhblck = 0;
-        for i = 1:SampleSize
-            mhfile = SampleAddress(i,3);
-            mhblck = SampleAddress(i,2);
-            if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
-                load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
-            end
-            pdraws(i,1) = {x2(SampleAddress(i,4),:)};
-            if info==2
-                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;
-            old_mhblck = mhblck;
-        end
-        clear('x2')
-        save([BaseName '_posterior_draws1.mat'],'pdraws','estim_params_')
-    else% The posterior draws are saved in xx files.
-        NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
-        NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
-        NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
-        linee = 0;
-        fnum  = 1;
+    NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
+    NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
+    NumberOfLinesLastFile = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
+    linee = 0;
+    fnum  = 1;
+    if fnum < NumberOfFiles
         pdraws = cell(NumberOfDrawsPerFile,info);
-        old_mhfile = 0;
-        old_mhblck = 0;
-        for i=1:SampleSize
-            linee = linee+1;
-            mhfile = SampleAddress(i,3);
-            mhblck = SampleAddress(i,2);
-            if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
-                load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
-            end
-            pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
-            if info==2
-                M_ = set_parameters_locally(M_,pdraws{linee,1});
-                [dr,~,M_.params] = compute_decision_rules(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
-                pdraws(linee,2) = { dr };
-            end
-            old_mhfile = mhfile;
-            old_mhblck = mhblck;
-            if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile
-                linee = 0;
-                save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
-                fnum = fnum+1;
-                if fnum < NumberOfFiles
-                    pdraws = cell(NumberOfDrawsPerFile,info);
-                else
-                    pdraws = cell(NumberOfLines,info);
-                end
+    else
+        pdraws = cell(NumberOfLinesLastFile,info);
+    end
+    for i=1:SampleSize
+        linee = linee+1;
+        pdraws(i,1) = {parameter_draws(i,:)};
+        if info==2
+            M_ = set_parameters_locally(M_,pdraws{linee,1});
+            [dr,~,M_.params] = compute_decision_rules(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
+            pdraws(linee,2) = { dr };
+        end
+        if (fnum < NumberOfFiles && linee == NumberOfDrawsPerFile) || (fnum <= NumberOfFiles && i==SampleSize)
+            linee = 0;
+            save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
+            fnum = fnum+1;
+            if fnum < NumberOfFiles
+                pdraws = cell(NumberOfDrawsPerFile,info);
+            else
+                pdraws = cell(NumberOfLinesLastFile,info);
             end
         end
-        save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
     end
 end
\ No newline at end of file
diff --git a/matlab/estimation/set_number_of_subdraws.m b/matlab/estimation/set_number_of_subdraws.m
new file mode 100644
index 0000000000000000000000000000000000000000..c6dcf6911d635903353e46e964907fdcc218f070
--- /dev/null
+++ b/matlab/estimation/set_number_of_subdraws.m
@@ -0,0 +1,83 @@
+function [sub_draws, error_flag, NumberOfDrawsPerChain]=set_number_of_subdraws(M_,options_)
+% function [sub_draws, error_flag, NumberOfDrawsPerChain]=set_number_of_subdraws(M_,options_)
+% Set option field sub_draws based on size of available sample
+% Inputs:
+%   M_          [structure]     Dynare model structure
+%   options_    [structure]     Dynare options structure
+%
+% Outputs:
+%   sub_draws   [scalar]        number of sub-draws
+%   error_flag  [boolean]       error indicate of insufficient draws are
+%                               available
+%   NumberOfDrawsPerChain [integer]     number of available posterior draws
+%                                       per chain
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+
+error_flag=false;
+if ~issmc(options_)
+    record=load_last_mh_history_file([M_.dname filesep 'metropolis'], M_.fname);
+    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+    NumberOfDrawsPerChain = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+    NumberOfDraws=NumberOfDrawsPerChain*record.Nblck;
+    if isempty(options_.sub_draws)
+        sub_draws = min(options_.posterior_max_subsample_draws, ceil(NumberOfDraws));
+    else
+        if options_.sub_draws>NumberOfDraws*record.Nblck
+            skipline()
+            disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*options_.mh_nblck) ')!'])
+            disp('You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).')
+            skipline()
+            error_flag = true;
+            sub_draws=options_.sub_draws; %pass through
+        else
+            sub_draws=options_.sub_draws; %pass through
+        end
+    end
+else
+    if ishssmc(options_)
+        % Load draws from the posterior distribution
+        pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
+        posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+        NumberOfDraws = size(posterior.particles,2);
+        NumberOfDrawsPerChain=NumberOfDraws;
+    elseif isdime(options_)
+        posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
+        tune = posterior.tune;
+        chains = reshape(posterior.chains(end-tune:end,:,:), [], size(posterior.chains, 3));
+        NumberOfDraws=size(chains,1);
+        NumberOfDrawsPerChain=NumberOfDraws;
+    else
+        error('set_number_of_subdraws:: case should not happen. Please contact the developers')
+    end
+    if isempty(options_.sub_draws)
+        sub_draws = min(options_.posterior_max_subsample_draws, NumberOfDraws);
+    else
+        if options_.sub_draws>NumberOfDraws
+            skipline()
+            disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws) ')!'])
+            disp('You can either change the value of sub_draws or run another MCMC.')
+            skipline()
+            error_flag = true;
+            sub_draws=options_.sub_draws; %pass through
+        else
+            sub_draws=options_.sub_draws; %pass through
+        end
+    end
+end
\ No newline at end of file
diff --git a/matlab/estimation/trace_plot.m b/matlab/estimation/trace_plot.m
index a3bb401abb322d9175fead4981e1f1d0065e47b7..2d0dd7c2d6e2c6edbfaa6b2771d1e5157c56d4d1 100644
--- a/matlab/estimation/trace_plot.m
+++ b/matlab/estimation/trace_plot.m
@@ -53,17 +53,7 @@ if isempty(column)
 end
 
 if ~issmc(options_)
-    % Get informations about the posterior draws:
-    MetropolisFolder = CheckPath('metropolis',M_.dname);
-    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;
-
+    options_.mh_drop=0; % locally, do not drop as we want all draws
     n_nblocks_to_plot=length(blck);
 else
     if ishssmc(options_)
@@ -73,14 +63,11 @@ else
     end
 end
 
+[~, ~, TotalNumberOfMhDraws]=set_number_of_subdraws(M_,options_); 
+
 if n_nblocks_to_plot==1
 % Get all the posterior draws:
-    if ishssmc(options_)
-        PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,[],column);
-        TotalNumberOfMhDraws=length(PosteriorDraws);
-    else
-        PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,M_.fname,column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck);
-    end
+    PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,M_.fname,column, 1, 1, blck);
 else
     PosteriorDraws=NaN(TotalNumberOfMhDraws,n_nblocks_to_plot);
     save_string='';
@@ -88,7 +75,7 @@ else
         title_string_tex='';
     end
     for block_iter=1:n_nblocks_to_plot
-        PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck(block_iter));
+        PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, 1, 1, blck(block_iter));
         save_string=[save_string,'_',num2str(blck(block_iter))];
         if options_.TeX
             title_string_tex=[title_string_tex, ', ' num2str(blck(block_iter))];
diff --git a/matlab/estimation/variance_decomposition_ME_mc_analysis.m b/matlab/estimation/variance_decomposition_ME_mc_analysis.m
index 91939d23b351eb9fea52dee2358d0042ad91d28f..22c84544d5aa6b9600333fa432a0e42e39c80974 100644
--- a/matlab/estimation/variance_decomposition_ME_mc_analysis.m
+++ b/matlab/estimation/variance_decomposition_ME_mc_analysis.m
@@ -25,7 +25,7 @@ function oo_ = variance_decomposition_ME_mc_analysis(NumberOfSimulations,type,dn
 
 
 
-% Copyright © 2017 Dynare Team
+% Copyright © 2017-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -43,8 +43,9 @@ function oo_ = variance_decomposition_ME_mc_analysis(NumberOfSimulations,type,dn
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
diff --git a/matlab/estimation/variance_decomposition_mc_analysis.m b/matlab/estimation/variance_decomposition_mc_analysis.m
index d55de6989b1107d026ff453358b2686ccaa46724..a34a507cf50e7359a8024d628d68ccc4bb64eac0 100644
--- a/matlab/estimation/variance_decomposition_mc_analysis.m
+++ b/matlab/estimation/variance_decomposition_mc_analysis.m
@@ -25,7 +25,7 @@ function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname
 
 
 
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -43,8 +43,9 @@ function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
diff --git a/matlab/list_of_functions_to_be_cleared.m b/matlab/list_of_functions_to_be_cleared.m
index 90c7187b3c8c34dae908937f592ff520c492f735..dc1a5eadcf86bd551cfc23a038927d1b6f0056f3 100644
--- a/matlab/list_of_functions_to_be_cleared.m
+++ b/matlab/list_of_functions_to_be_cleared.m
@@ -1 +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', '+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', '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'};
diff --git a/matlab/moments/compute_moments_varendo.m b/matlab/moments/compute_moments_varendo.m
index 4b26fb17a6dd4c7890e2ca9d8983f780f168175b..874e51c9e47f9bcf9ce4c9d56cd2d7133f7e7e2c 100644
--- a/matlab/moments/compute_moments_varendo.m
+++ b/matlab/moments/compute_moments_varendo.m
@@ -64,10 +64,7 @@ if strcmpi(type,'posterior')
 elseif strcmpi(type,'prior')
     posterior = 0;
     if nargin==5
-        var_list_ = options_.prior_analysis_endo_var_list;
-        if isempty(var_list_)
-            options_.prior_analysis_var_list = options_.varobs;
-        end
+        var_list_ = options_.varobs;
     end
     if isfield(oo_,'PriorTheoreticalMoments')
         oo_=rmfield(oo_,'PriorTheoreticalMoments');
@@ -103,7 +100,7 @@ if posterior
 else
     for i=1:NumberOfEndogenousVariables
         for j=i:NumberOfEndogenousVariables
-            oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, [], options_, M_, oo_, estim_params_);
+            oo_ = prior_analysis('variance', var_list_{i}, var_list_{j}, NumberOfLags, options_, M_, oo_, estim_params_);
         end
     end
 end
diff --git a/tests/moments/fs2000_post_moments.mod b/tests/moments/fs2000_post_moments.mod
index f97968ca137b1663392da45507cecfa9d80c82e0..b3bc124c18f876c50dd254a559c5c59af63a4b23 100644
--- a/tests/moments/fs2000_post_moments.mod
+++ b/tests/moments/fs2000_post_moments.mod
@@ -285,6 +285,12 @@ for var_iter_1=1:nvars
     end
 end
 
+
+options_.prior_mc=100;
+options_.qz_criterium = 1+1e-6;
+oo_ = compute_moments_varendo('prior',options_,M_,oo_,estim_params_,M_.endo_names(1:M_.orig_endo_nbr));
+
+
 /*
  * The following lines were used to generate the data file. If you want to
  * generate another random data file, comment the "estimation" line and uncomment