diff --git a/.gitignore b/.gitignore
index 9556e47b5b1c43eba4d48a3eab42b49a671deaf4..144f6d24cf8c050d634a591ff8ad22868f87fb35 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,7 @@
 *.pdf
 
-tex/auto
\ No newline at end of file
+tex/auto
+
+model/+fs2000/
+model/fs2000/
+model/fs2000.log
diff --git a/matlab/load_all_distributions.m b/matlab/load_all_distributions.m
new file mode 100644
index 0000000000000000000000000000000000000000..f468abec9fcd076224051af4dbfccd32a323561d
--- /dev/null
+++ b/matlab/load_all_distributions.m
@@ -0,0 +1,12 @@
+function particles = load_all_distributions(modelname, nparticles, idparam)
+
+% (C) Stéphane Adjemian 2024
+
+n = number_of_smc_files(modelname);
+
+particles = NaN(nparticles, n);
+
+for i=1:n
+    p = load(sprintf('%s%shssmc%sparticles-%u-%u.mat', modelname, filesep, filesep, i, n));
+    particles(:,i) = transpose(p.particles(idparam,:));
+end
diff --git a/matlab/load_posterior_draws.m b/matlab/load_posterior_draws.m
new file mode 100644
index 0000000000000000000000000000000000000000..a37d12d7dc10f9e0855904fca37d76a63a478c14
--- /dev/null
+++ b/matlab/load_posterior_draws.m
@@ -0,0 +1,7 @@
+function p = load_posterior_draws(modelname)
+
+% (C) Stéphane Adjemian 2024
+
+n = number_of_smc_files(modelname);
+
+p = load(sprintf('%s%shssmc%sparticles-%u-%u.mat', modelname, filesep, filesep, n, n));
diff --git a/matlab/number_of_smc_files.m b/matlab/number_of_smc_files.m
new file mode 100644
index 0000000000000000000000000000000000000000..a56c4f0b664c22268aba61d10fb40a9c44d1bedc
--- /dev/null
+++ b/matlab/number_of_smc_files.m
@@ -0,0 +1,5 @@
+function n = number_of_smc_files(modelname)
+
+% (C) Stéphane Adjemian 2024
+
+n = length(dir(sprintf('%s%shssmc%sparticles-*-*.mat', modelname, filesep, filesep)));
diff --git a/matlab/plot_all_distributions.m b/matlab/plot_all_distributions.m
new file mode 100644
index 0000000000000000000000000000000000000000..11f672a1d5ad937ba9f3dc0bf7f2266d7c8b70c2
--- /dev/null
+++ b/matlab/plot_all_distributions.m
@@ -0,0 +1,24 @@
+function densities = plot_all_distributions(modelname, nparticles, idparam)
+
+% (C) Stéphane Adjemian, 2024
+
+    particles = load_all_distributions(modelname, nparticles, idparam);
+
+    [N,n] = size(particles);
+
+    densities = cell(n,2);
+
+    for i=1:n
+        bandwidth = mh_optimal_bandwidth(particles(:,i), N, 0, 'gaussian');
+        [x, f] = kernel_density_estimate(particles(:,i), 512, N, bandwidth, 'gaussian');
+        densities{i,1} = x;
+        densities{i,2} = f;
+    end
+
+    plot3(densities{1,1}, ones(512,1), densities{1,2}, '-k')
+    hold on
+    for i=2:25, plot3(densities{i,1}, i*ones(512,1), densities{i,2}, '-k'), end
+    axis tight
+    grid on
+    hold off
+end
diff --git a/model/fs2000.mod b/model/fs2000.mod
index b81f9a510b4d9d227ef0e01518d8a8374b80cc6c..17b676d4a0c0895e31081ee9ae98a60eb7bc6bca 100644
--- a/model/fs2000.mod
+++ b/model/fs2000.mod
@@ -89,3 +89,6 @@ estimation(order=1, datafile='fsdat_simul.m', nobs=192, loglinear,
                                       'particles', 20000,
                                       'scale',.5,
                                       'target', .25));
+
+addpath('../matlab');
+plot_all_distributions('fs2000', 20000, 7); // autoregressive parameter (rho)