From 18eede9d057dd086e41b3d30745f84e7109c3a3f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Karam=C3=A9?=
 <frederic.karame@univ-lemans.fr>
Date: Sat, 25 May 2024 16:21:39 +0200
Subject: [PATCH] Bug fix and simplification of the Student's pdf.

dynare was crashing without the statistical toolbox from Mathworks.
---
 matlab/nonlinear-filters/auxiliary_particle_filter.m | 4 +++-
 matlab/nonlinear-filters/online_auxiliary_filter.m   | 4 +++-
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/matlab/nonlinear-filters/auxiliary_particle_filter.m b/matlab/nonlinear-filters/auxiliary_particle_filter.m
index 713864308..f0c67c34a 100644
--- a/matlab/nonlinear-filters/auxiliary_particle_filter.m
+++ b/matlab/nonlinear-filters/auxiliary_particle_filter.m
@@ -132,7 +132,9 @@ for t=1:sample_size
     end
     PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
     z = sum(PredictionError.*(H\PredictionError),1) ;
-    tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ;
+%    tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ;
+    ddl = 3 ;
+    tau_tilde = weights.*(exp(gammaln((ddl + 1) / 2) - gammaln(ddl/2))./(sqrt(ddl*pi).*(1 + (z.^2)./ddl).^((ddl + 1)/2))+1e-99) ;
     tau_tilde = tau_tilde/sum(tau_tilde) ;
     indx = resample(0,tau_tilde',ParticleOptions);
     if pruning
diff --git a/matlab/nonlinear-filters/online_auxiliary_filter.m b/matlab/nonlinear-filters/online_auxiliary_filter.m
index 1c98b2bd3..f8e3c54d9 100644
--- a/matlab/nonlinear-filters/online_auxiliary_filter.m
+++ b/matlab/nonlinear-filters/online_auxiliary_filter.m
@@ -191,7 +191,9 @@ for t=1:sample_size
             PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
             % Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails.
             z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ;
-            tau_tilde(i) = weights(i).*(tpdf(z, 3*ones(size(z)))+1e-99) ;
+            ddl = 3 ;
+            %tau_tilde(i) = weights(i).*(tpdf(z,3*ones(size(z)))+1e-99) ;
+            tau_tilde(i) = weights(i)*(exp(gammaln((ddl + 1) / 2) - gammaln(ddl/2))/(sqrt(ddl*pi)*(1 + (z^2)/ddl)^((ddl + 1)/2))+1e-99) ;
         else
             tau_tilde(i) = 0 ;
         end
-- 
GitLab