diff --git a/matlab/+gsa/map_calibration.m b/matlab/+gsa/map_calibration.m
index aa68e09008c4c03115b9c05622778aa0a5722181..3c9757bd0d2f844d2df271ead766a123bf2e88d1 100644
--- a/matlab/+gsa/map_calibration.m
+++ b/matlab/+gsa/map_calibration.m
@@ -212,7 +212,7 @@ if ~isempty(indx_irf)
         indx_irf_matrix(:,plot_indx(ij)) = indx_irf_matrix(:,plot_indx(ij)) + indx_irf(:,ij);
         for ik=1:size(mat_irf{ij},2)
             [Mean,Median,Var,HPD,Distrib] = ...
-                posterior_moments(mat_irf{ij}(:,ik),0,options_.mh_conf_sig);
+                posterior_moments(mat_irf{ij}(:,ik),options_.mh_conf_sig);
             irf_mean{plot_indx(ij)} = [irf_mean{plot_indx(ij)}; Mean];
             irf_median{plot_indx(ij)} = [irf_median{plot_indx(ij)}; Median];
             irf_var{plot_indx(ij)} = [irf_var{plot_indx(ij)}; Var];
@@ -417,7 +417,7 @@ if ~isempty(indx_moment)
         indx_moment_matrix(:,plot_indx(ij)) = indx_moment_matrix(:,plot_indx(ij)) + indx_moment(:,ij);
         for ik=1:size(mat_moment{ij},2)
             [Mean,Median,Var,HPD,Distrib] = ...
-                posterior_moments(mat_moment{ij}(:,ik),0,options_.mh_conf_sig);
+                posterior_moments(mat_moment{ij}(:,ik),options_.mh_conf_sig);
             moment_mean{plot_indx(ij)} = [moment_mean{plot_indx(ij)}; Mean];
             moment_median{plot_indx(ij)} = [moment_median{plot_indx(ij)}; Median];
             moment_var{plot_indx(ij)} = [moment_var{plot_indx(ij)}; Var];
diff --git a/matlab/+gsa/reduced_form_mapping.m b/matlab/+gsa/reduced_form_mapping.m
index de575f4cb2a2e1fa082ba0dbe72e78f5191dd28e..5baca5682d0797b547b637c9aab905d6d4f029d4 100644
--- a/matlab/+gsa/reduced_form_mapping.m
+++ b/matlab/+gsa/reduced_form_mapping.m
@@ -720,7 +720,7 @@ function indmcf = redform_mcf(y0, x0, options_mcf, options_, fname, parnames, es
 
 hh_fig=dyn_figure(options_.nodisplay,'name',options_mcf.amcf_title);
 
-[~, ~, ~, ~, post_deciles] = posterior_moments(y0,1,0.9);
+[~, ~, ~, ~, post_deciles] = posterior_moments(y0,0.9);
 post_deciles = [-inf; post_deciles; inf];
 
 for jt=1:10
diff --git a/matlab/+identification/plot.m b/matlab/+identification/plot.m
index d20531f752f31c7192a18e0a84c2a8200f5f8f3a..14701e9157f22064cea07b9e9ae953af99481ab6 100644
--- a/matlab/+identification/plot.m
+++ b/matlab/+identification/plot.m
@@ -493,7 +493,7 @@ else
             post_median=NaN(1,nparam);
             hpd_interval=NaN(nparam,2);
             for i=1:nparam
-                [~, post_median(:,i), ~, hpd_interval(i,:)] = posterior_moments(VVV(:,i),0,0.9);
+                [~, post_median(:,i), ~, hpd_interval(i,:)] = posterior_moments(VVV(:,i),0.9);
             end
             bar(post_median)
             hold on, plot(hpd_interval,'--*r'),
diff --git a/matlab/+occbin/forecast.m b/matlab/+occbin/forecast.m
index 07c8a94f139f8a6b3c6d7a6202a1b240475cd512..3cc855223134928d0fe1a29f177ee47ea8f873f4 100644
--- a/matlab/+occbin/forecast.m
+++ b/matlab/+occbin/forecast.m
@@ -88,7 +88,7 @@ if opts.replic
 
     for i=1:M_.endo_nbr
         for j=1:forecast
-            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),1,options_.forecasts.conf_sig);
+            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),options_.forecasts.conf_sig);
         end
         oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})=post_mean;
         oo_.occbin.linear_forecast.Median.(M_.endo_names{i})=post_median;
@@ -99,7 +99,7 @@ if opts.replic
         oo_.occbin.linear_forecast.Min.(M_.endo_names{i})=zlinmin(:,i);
         oo_.occbin.linear_forecast.Max.(M_.endo_names{i})=zlinmax(:,i);
         for j=1:forecast
-            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),1,options_.forecasts.conf_sig);
+            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),options_.forecasts.conf_sig);
         end
         oo_.occbin.forecast.Mean.(M_.endo_names{i})=post_mean;
         oo_.occbin.forecast.Median.(M_.endo_names{i})=post_median;
diff --git a/matlab/estimation/GetPosteriorParametersStatistics.m b/matlab/estimation/GetPosteriorParametersStatistics.m
index bc20b4415e277ae4229fd809eac9ee5b22a8042f..cd791de522c7007e206bf6dd883783c9f2a369c9 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -105,7 +105,7 @@ if np
     for i=1:np
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_)
             draws = getalldraws(ip);
-            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
             name = bayestopt_.name{ip};
             oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
         else
@@ -114,7 +114,7 @@ if np
                 [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
             catch
                 draws = getalldraws(ip);
-                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
                 name = bayestopt_.name{ip};
                 oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
             end
@@ -148,7 +148,7 @@ if nvx
     for i=1:nvx
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
             draws = getalldraws(ip);
-            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
             k = estim_params_.var_exo(i,1);
             name = M_.exo_names{k};
             oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@@ -160,7 +160,7 @@ if nvx
                 [post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
             catch
                 draws = getalldraws(ip);
-                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
                 k = estim_params_.var_exo(i,1);
                 name = M_.exo_names{k};
                 oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@@ -191,7 +191,7 @@ if nvn
     for i=1:nvn
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
             draws = getalldraws(ip);
-            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
             name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
             oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
         else
@@ -200,7 +200,7 @@ if nvn
                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
             catch
                 draws = getalldraws(ip);
-                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,1,options_.mh_conf_sig);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,options_.mh_conf_sig);
                 name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
                 oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
             end
@@ -230,7 +230,7 @@ if ncx
     for i=1:ncx
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
             draws = getalldraws(ip);
-            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,1,options_.mh_conf_sig);
+            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,options_.mh_conf_sig);
             k1 = estim_params_.corrx(i,1);
             k2 = estim_params_.corrx(i,2);
             name = sprintf('%s,%s', M_.exo_names{k1}, M_.exo_names{k2});
@@ -247,7 +247,7 @@ if ncx
                 [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
             catch
                 draws = getalldraws(ip);
-                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
                 k1 = estim_params_.corrx(i,1);
                 k2 = estim_params_.corrx(i,2);
                 name = sprintf('%s,%s', M_.exo_names{k1}, M_.exo_names{k2});
@@ -281,7 +281,7 @@ if ncn
     for i=1:ncn
         if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
             draws = getalldraws(ip);
-            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
             k1 = estim_params_.corrn(i,1);
             k2 = estim_params_.corrn(i,2);
             name = sprintf('%s,%s', M_.endo_names{k1}, M_.endo_names{k2});
@@ -296,7 +296,7 @@ if ncn
                 [post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
             catch
                 draws = getalldraws(ip);
-                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
                 k1 = estim_params_.corrn(i,1);
                 k2 = estim_params_.corrn(i,2);
                 name = sprintf('%s,%s', M_.endo_names{k1}, M_.endo_names{k2});
diff --git a/matlab/estimation/PosteriorIRF.m b/matlab/estimation/PosteriorIRF.m
index b539d13e9402ccebb66cce5b912a7164881a50f5..d53bc344d17795840f4a8fe34d2bf9d7dc5c3e0f 100644
--- a/matlab/estimation/PosteriorIRF.m
+++ b/matlab/estimation/PosteriorIRF.m
@@ -283,7 +283,7 @@ for file = 1:NumberOfIRFfiles_dsge
             for k = 1:size(temp.STOCK_IRF_DSGE,1)
                 kk = k+kdx;
                 [MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),...
-                 DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(temp.STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
+                 DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(temp.STOCK_IRF_DSGE(k,j,i,:)),options_.mh_conf_sig);
             end
         end
     end
@@ -323,7 +323,7 @@ if MAX_nirfs_dsgevar
                     kk = k+kdx;
                     [MeanIRFdsgevar(kk,j,i),MedianIRFdsgevar(kk,j,i),VarIRFdsgevar(kk,j,i),...
                      HPDIRFdsgevar(kk,:,j,i),DistribIRFdsgevar(kk,:,j,i)] = ...
-                        posterior_moments(squeeze(temp.STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig);
+                        posterior_moments(squeeze(temp.STOCK_IRF_BVARDSGE(k,j,i,:)),options_.mh_conf_sig);
                 end
             end
         end
diff --git a/matlab/estimation/pm3.m b/matlab/estimation/pm3.m
index 5465b81ab287b82759ca7d3df8ace8fe5c099504..820c83ca4935b1e001c015f6469dea50146c94c8 100644
--- a/matlab/estimation/pm3.m
+++ b/matlab/estimation/pm3.m
@@ -193,10 +193,10 @@ if strcmp(var_type,'_trend_coeff') %two dimensional arrays
     for i = 1:nvar
         if options_.estimation.moments_posterior_density.indicator
             [Mean(1,i),Median(1,i),Var(1,i),HPD(:,1,i),Distrib(:,1,i),Density(:,:,1,i)] = ...
-                posterior_moments(squeeze(stock1(SelecVariables(i),:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
+                posterior_moments(squeeze(stock1(SelecVariables(i),:)),options_.mh_conf_sig,options_.estimation.moments_posterior_density);
         else
             [Mean(1,i),Median(1,i),Var(1,i),HPD(:,1,i),Distrib(:,1,i)] = ...
-                posterior_moments(squeeze(stock1(SelecVariables(i),:)),0,options_.mh_conf_sig);
+                posterior_moments(squeeze(stock1(SelecVariables(i),:)),options_.mh_conf_sig);
         end
     end
 else %three dimensional arrays
@@ -204,21 +204,21 @@ else %three dimensional arrays
         for j = 1:n2
             if options_.estimation.moments_posterior_density.indicator
                 [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i),Density(:,:,j,i)] = ...
-                    posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
+                    posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),options_.mh_conf_sig,options_.estimation.moments_posterior_density);
             else
                 [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ...
-                    posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig);
+                    posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),options_.mh_conf_sig);
             end
             if filter_step_ahead_indicator
                 if options_.estimation.moments_posterior_density.indicator
                     for K_step = 1:length(options_.filter_step_ahead)
                         [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j),Density_filter_step_ahead(:,:,K_step,i,j) ] = ...
-                            posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
+                            posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),options_.mh_conf_sig,options_.estimation.moments_posterior_density);
                     end
                 else
                     for K_step = 1:length(options_.filter_step_ahead)
                         [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ...
-                            posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig);
+                            posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),options_.mh_conf_sig);
                     end
                 end
             end
diff --git a/matlab/estimation/posterior_moments.m b/matlab/estimation/posterior_moments.m
index e4d647bfdafd2f474a085583d8aa8862d8161566..a12cc5ba0522996c9666d397f5626f2130a80830 100644
--- a/matlab/estimation/posterior_moments.m
+++ b/matlab/estimation/posterior_moments.m
@@ -1,10 +1,11 @@
-function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info,mh_conf_sig,kernel_options)
+function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,mh_conf_sig,kernel_options)
+% [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,mh_conf_sig,kernel_options)
 % Computes posterior mean, median, variance, HPD interval, deciles, and density from posterior draws.
 %
 % INPUTS
-%    xx            [double]    Vector of posterior draws (or prior draws ;-)
-%    info          [integer]   If equal to one the posterior density is estimated.
-%    mh_config_sig [double]    Scalar between 0 and 1 specifying the size of the HPD interval.
+%    xx                 [double]    Vector of posterior draws (or prior draws ;-)
+%    mh_config_sig      [double]    Scalar between 0 and 1 specifying the size of the HPD interval.
+%    kernel_options     [structure] options for kernel density estimate
 %
 %
 % OUTPUTS
@@ -21,7 +22,7 @@ function [post_mean, post_median, post_var, hpd_interval, post_deciles, density]
 %                                                   kernel_density_estimate.m.
 %
 
-% Copyright © 2005-2017 Dynare Team
+% Copyright © 2005-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -38,20 +39,9 @@ function [post_mean, post_median, post_var, hpd_interval, post_deciles, density]
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if nargin<4
-    number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.
-    bandwidth = 0;                    % Rule of thumb optimal bandwidth parameter.
-    kernel_function = 'gaussian';     % Gaussian kernel for Fast Fourrier Transform approximaton.
-else
-    number_of_grid_points = kernel_options.gridpoints;
-    bandwidth = kernel_options.bandwidth;
-    kernel_function = kernel_options.kernel_function;
-end
-
 xx = xx(:);
 xx = sort(xx);
 
-
 post_mean = mean(xx);
 post_median = median(xx);
 post_var = var(xx);
@@ -84,8 +74,18 @@ if length(xx)>9
 else
     post_deciles=NaN(9,1);
 end
+
 density = [];
-if info
+if nargout == 6
+    if nargin<3
+        number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.
+        bandwidth = 0;                    % Rule of thumb optimal bandwidth parameter.
+        kernel_function = 'gaussian';     % Gaussian kernel for Fast Fourrier Transform approximaton.
+    else
+        number_of_grid_points = kernel_options.gridpoints;
+        bandwidth = kernel_options.bandwidth;
+        kernel_function = kernel_options.kernel_function;
+    end
     if post_var > 1e-12
         optimal_bandwidth = mh_optimal_bandwidth(xx,number_of_draws,bandwidth,kernel_function);
         [density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,...
diff --git a/matlab/estimation/variance_decomposition_ME_mc_analysis.m b/matlab/estimation/variance_decomposition_ME_mc_analysis.m
index 9d2e7e85e383e2d658e39e6e8d4ae89d159794bb..a9a4cc30788249dc9ab810e98cc261703bbe157d 100644
--- a/matlab/estimation/variance_decomposition_ME_mc_analysis.m
+++ b/matlab/estimation/variance_decomposition_ME_mc_analysis.m
@@ -95,10 +95,10 @@ end
 
 if options_.estimation.moments_posterior_density.indicator
     [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 else
     [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
-        posterior_moments(tmp,0,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 end
 
 oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.Mean.(var).(exo) = p_mean;
diff --git a/matlab/estimation/variance_decomposition_mc_analysis.m b/matlab/estimation/variance_decomposition_mc_analysis.m
index dc4f6f35e51daee54fb90422602ac41f18d314b3..51bc45532a5e790f30b3edfa2112c009a588cfa0 100644
--- a/matlab/estimation/variance_decomposition_mc_analysis.m
+++ b/matlab/estimation/variance_decomposition_mc_analysis.m
@@ -93,10 +93,10 @@ end
 
 if options_.estimation.moments_posterior_density.indicator
     [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 else
     [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
-        posterior_moments(tmp,0,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 end
 
 oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.Mean.(var).(exo) = p_mean;
diff --git a/matlab/moments/conditional_variance_decomposition_ME_mc_analysis.m b/matlab/moments/conditional_variance_decomposition_ME_mc_analysis.m
index c4d4e349694592eaded997db30f36fa436a72ca1..061179eb87661297b12c5d37fab723b02075da65 100644
--- a/matlab/moments/conditional_variance_decomposition_ME_mc_analysis.m
+++ b/matlab/moments/conditional_variance_decomposition_ME_mc_analysis.m
@@ -113,11 +113,11 @@ p_hpdsup = NaN(1,length(Steps));
 for i=1:length(Steps)
     if options_.estimation.moments_posterior_density.indicator
         [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ...
-            posterior_moments(tmp(:,i),1,mh_conf_sig);
+            posterior_moments(tmp(:,i),mh_conf_sig);
         p_density(:,:,i) = pp_density;
     else
         [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ...
-            posterior_moments(tmp(:,i),0,mh_conf_sig);
+            posterior_moments(tmp(:,i),mh_conf_sig);
     end
     p_mean(i) = pp_mean;
     p_median(i) = pp_median;
diff --git a/matlab/moments/conditional_variance_decomposition_mc_analysis.m b/matlab/moments/conditional_variance_decomposition_mc_analysis.m
index 74d86fffd5558934484a4cf6860e8b8188cfd26f..81aeb855fbd6941d5e1905f9015069f6bb543391 100644
--- a/matlab/moments/conditional_variance_decomposition_mc_analysis.m
+++ b/matlab/moments/conditional_variance_decomposition_mc_analysis.m
@@ -104,11 +104,11 @@ p_hpdsup = NaN(1,length(Steps));
 for i=1:length(Steps)
     if options_.estimation.moments_posterior_density.indicator
         [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ...
-            posterior_moments(tmp(:,i),1,mh_conf_sig);
+            posterior_moments(tmp(:,i),mh_conf_sig);
         p_density(:,:,i) = pp_density;
     else
         [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ...
-            posterior_moments(tmp(:,i),0,mh_conf_sig);
+            posterior_moments(tmp(:,i),mh_conf_sig);
     end
     p_mean(i) = pp_mean;
     p_median(i) = pp_median;
diff --git a/matlab/moments/correlation_mc_analysis.m b/matlab/moments/correlation_mc_analysis.m
index ecf58c4a458b6f179f8ca7351c57f28f99793529..b6df2e223332c35439f08a8a039ad8a140632fdd 100644
--- a/matlab/moments/correlation_mc_analysis.m
+++ b/matlab/moments/correlation_mc_analysis.m
@@ -93,10 +93,10 @@ end
 name = [ var1 '.' var2 ];
 if options_.estimation.moments_posterior_density.indicator
     [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 else
     [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
-        posterior_moments(tmp,0,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 end
 if isfield(oo_,[ TYPE 'TheoreticalMoments'])
     temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
diff --git a/matlab/moments/covariance_mc_analysis.m b/matlab/moments/covariance_mc_analysis.m
index 5215a74a7da473dd1c6096cff5fe3361bd678bcd..9d7cab85a37b24492df1a5bf63a01e1f492f40d5 100644
--- a/matlab/moments/covariance_mc_analysis.m
+++ b/matlab/moments/covariance_mc_analysis.m
@@ -110,11 +110,11 @@ end
 
 if options_.estimation.moments_posterior_density.indicator
     [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
     oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.density.(var1).(var2) = density;
 else
     [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
-        posterior_moments(tmp,0,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 end
 oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Mean.(var1).(var2) = p_mean;
 oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Median.(var1).(var2) = p_median;
@@ -126,11 +126,11 @@ oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.deciles.(var1).(var2) = p_dec
 if options_.contemporaneous_correlation
     if options_.estimation.moments_posterior_density.indicator
         [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
-            posterior_moments(tmp_corr_mat,1,mh_conf_sig);
+            posterior_moments(tmp_corr_mat,mh_conf_sig);
         oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.density.(var1).(var2) = density;
     else
         [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
-            posterior_moments(tmp_corr_mat,0,mh_conf_sig);
+            posterior_moments(tmp_corr_mat,mh_conf_sig);
     end
     oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.Mean.(var1).(var2) = p_mean;
     oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.Median.(var1).(var2) = p_median;