From 8eaccada384548fe52ff1e0591955fae65dfe403 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Sat, 21 Dec 2019 10:37:48 +0100 Subject: [PATCH] Allow k order approximation in Auxiliary Particle Filter. Ref. dynare#1673 --- src/auxiliary_particle_filter.m | 79 +++++++++++++++------------------ 1 file changed, 35 insertions(+), 44 deletions(-) diff --git a/src/auxiliary_particle_filter.m b/src/auxiliary_particle_filter.m index 61a43bc13..905a569d2 100644 --- a/src/auxiliary_particle_filter.m +++ b/src/auxiliary_particle_filter.m @@ -1,9 +1,9 @@ -function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions) +function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model) % Evaluates the likelihood of a nonlinear model with the auxiliary particle filter % allowing eventually resampling. % -% Copyright (C) 2011-2017 Dynare Team +% Copyright (C) 2011-2019 Dynare Team % % This file is part of Dynare (particles module). % @@ -20,9 +20,6 @@ function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptio % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -persistent init_flag mf0 mf1 number_of_particles -persistent sample_size number_of_observed_variables number_of_structural_innovations - % Set default if isempty(start) start = 1; @@ -36,24 +33,24 @@ steadystate = ReducedForm.steadystate; constant = ReducedForm.constant; state_variables_steady_state = ReducedForm.state_variables_steady_state; -% Set persistent variables. -if isempty(init_flag) - mf0 = ReducedForm.mf0; - mf1 = ReducedForm.mf1; - sample_size = size(Y,2); - number_of_observed_variables = length(mf1); - number_of_structural_innovations = length(ReducedForm.Q); - number_of_particles = ParticleOptions.number_of_particles; - init_flag = 1; -end +mf0 = ReducedForm.mf0; +mf1 = ReducedForm.mf1; +sample_size = size(Y,2); +number_of_observed_variables = length(mf1); +number_of_structural_innovations = length(ReducedForm.Q); +number_of_particles = ParticleOptions.number_of_particles; -% Set local state space model (first order approximation). -ghx = ReducedForm.ghx; -ghu = ReducedForm.ghu; -% Set local state space model (second order approximation). -ghxx = ReducedForm.ghxx; -ghuu = ReducedForm.ghuu; -ghxu = ReducedForm.ghxu; +if ReducedForm.use_k_order_solver + dr = ReducedForm.dr; +else + % Set local state space model (first order approximation). + ghx = ReducedForm.ghx; + ghu = ReducedForm.ghu; + % Set local state space model (second order approximation). + ghxx = ReducedForm.ghxx; + ghuu = ReducedForm.ghuu; + ghxu = ReducedForm.ghxu; +end % Get covariance matrices Q = ReducedForm.Q; @@ -81,19 +78,6 @@ if pruning StateVectors_ = StateVectors; end -% Uncomment for building the mean average predictions based on a sparse -% grids of structural shocks. Otherwise, all shocks are set to 0 in the -% prediction. -% if ParticleOptions.proposal_approximation.cubature -% [nodes,nodes_weights] = spherical_radial_sigma_points(number_of_structural_innovations) ; -% nodes_weights = ones(size(nodes,1),1)*nodes_weights ; -% elseif ParticleOptions.proposal_approximation.unscented -% [nodes,nodes_weights,nodes_weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions); -% else -% error('Estimation: This approximation for the proposal is not implemented or unknown!') -% end -% nodes = (Q_lower_triangular_cholesky*(nodes'))' ; - nodes = zeros(1,number_of_structural_innovations) ; nodes_weights = ones(number_of_structural_innovations,1) ; @@ -109,15 +93,19 @@ for t=1:sample_size tmp_ = tmp_ + nodes_weights(i)*tmp1_ ; end else - tmp = 0 ; - for i=1:size(nodes) - tmp = tmp + nodes_weights(i)*local_state_space_iteration_2(yhat,nodes(i,:)'*ones(1,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); + if ReducedForm.use_k_order_solver + tmp = 0; + for i=1:size(nodes) + tmp = tmp + nodes_weights(i)*local_state_space_iteration_k(yhat, nodes(i,:)'*ones(1,number_of_particles), dr, Model, DynareOptions); + end + else + tmp = 0; + for i=1:size(nodes) + tmp = tmp + nodes_weights(i)*local_state_space_iteration_2(yhat,nodes(i,:)'*ones(1,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); + end end end PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); - %tau_tilde = weights.*(exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) + 1e-99) ; - % Replace Gaussian density with a Student density with 3 degrees of - % freedom for fat tails. z = sum(PredictionError.*(H\PredictionError),1) ; tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ; tau_tilde = tau_tilde/sum(tau_tilde) ; @@ -132,7 +120,11 @@ for t=1:sample_size [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); StateVectors_ = tmp_(mf0,:); else - tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); + if ReducedForm.use_k_order_solver + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + else + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + end end StateVectors = tmp(mf0,:); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); @@ -151,5 +143,4 @@ for t=1:sample_size end end -%plot(lik) ; -LIK = -sum(lik(start:end)); +LIK = -sum(lik(start:end)); \ No newline at end of file