From b438c177d7a6ec7b5437ed44dc7d3459736374f6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Argos=29?=
 <stepan@adjemian.eu>
Date: Wed, 18 Sep 2024 07:55:15 +0200
Subject: [PATCH] Add plot (SMC distributions).

---
 .gitignore                      |  6 +++++-
 matlab/load_all_distributions.m | 12 ++++++++++++
 matlab/load_posterior_draws.m   |  7 +++++++
 matlab/number_of_smc_files.m    |  5 +++++
 matlab/plot_all_distributions.m | 24 ++++++++++++++++++++++++
 model/fs2000.mod                |  3 +++
 6 files changed, 56 insertions(+), 1 deletion(-)
 create mode 100644 matlab/load_all_distributions.m
 create mode 100644 matlab/load_posterior_draws.m
 create mode 100644 matlab/number_of_smc_files.m
 create mode 100644 matlab/plot_all_distributions.m

diff --git a/.gitignore b/.gitignore
index 9556e47..144f6d2 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 0000000..f468abe
--- /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 0000000..a37d12d
--- /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 0000000..a56c4f0
--- /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 0000000..11f672a
--- /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 b81f9a5..17b676d 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)
-- 
GitLab