From ccce0cae1f4838047c049abf328c30ac7186ae8d Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 23 Sep 2021 13:04:19 +0200 Subject: [PATCH] Replace k-order solver by Fortran routine --- src/auxiliary_particle_filter.m | 5 +++-- src/conditional_filter_proposal.m | 2 +- src/gaussian_filter_bank.m | 2 +- src/gaussian_mixture_filter_bank.m | 2 +- src/measurement_equations.m | 2 +- src/nonlinear_kalman_filter.m | 3 ++- src/online_auxiliary_filter.m | 5 +++-- src/sequential_importance_particle_filter.m | 3 ++- 8 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/auxiliary_particle_filter.m b/src/auxiliary_particle_filter.m index 6a4e89f..228aeb3 100644 --- a/src/auxiliary_particle_filter.m +++ b/src/auxiliary_particle_filter.m @@ -42,6 +42,7 @@ number_of_particles = ParticleOptions.number_of_particles; if ReducedForm.use_k_order_solver dr = ReducedForm.dr; + [udr] = folded_to_unfolded_dr(dr, DynareModel, DynareOptions); else % Set local state space model (first order approximation). ghx = ReducedForm.ghx; @@ -96,7 +97,7 @@ for t=1:sample_size 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); + tmp = tmp + nodes_weights(i)*local_state_space_iteration_k(yhat, nodes(i,:)'*ones(1,number_of_particles), dr, Model, DynareOptions, udr); end else tmp = 0; @@ -121,7 +122,7 @@ for t=1:sample_size StateVectors_ = tmp_(mf0,:); else if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); else tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); end diff --git a/src/conditional_filter_proposal.m b/src/conditional_filter_proposal.m index 3b6b655..b816d13 100644 --- a/src/conditional_filter_proposal.m +++ b/src/conditional_filter_proposal.m @@ -79,7 +79,7 @@ epsilon = Q_lower_triangular_cholesky*nodes'; yhat = repmat(StateVectors-state_variables_steady_state, 1, size(epsilon, 2)); if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions,folded_to_unfolded_dr(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 diff --git a/src/gaussian_filter_bank.m b/src/gaussian_filter_bank.m index d736120..0731eca 100644 --- a/src/gaussian_filter_bank.m +++ b/src/gaussian_filter_bank.m @@ -81,7 +81,7 @@ StateVectors = sigma_points(1:number_of_state_variables,:); epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:); yhat = bsxfun(@minus, StateVectors, state_variables_steady_state); if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, folded_to_unfolded_dr(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 diff --git a/src/gaussian_mixture_filter_bank.m b/src/gaussian_mixture_filter_bank.m index 74b9dea..963b3eb 100644 --- a/src/gaussian_mixture_filter_bank.m +++ b/src/gaussian_mixture_filter_bank.m @@ -77,7 +77,7 @@ epsilon = bsxfun(@plus, StructuralShocksSqrtP*nodes3(:,number_of_state_variables StateVectors = bsxfun(@plus, StateSqrtP*nodes3(:,1:number_of_state_variables)', StateMu); yhat = bsxfun(@minus, StateVectors, state_variables_steady_state); if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions,folded_to_unfolded_dr(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 diff --git a/src/measurement_equations.m b/src/measurement_equations.m index 87edc75..45cee63 100644 --- a/src/measurement_equations.m +++ b/src/measurement_equations.m @@ -32,7 +32,7 @@ state_variables_steady_state = ReducedForm.state_variables_steady_state; number_of_structural_innovations = length(ReducedForm.Q); yhat = bsxfun(@minus, StateVectors, state_variables_steady_state); if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions, folded_to_unfolded_dr(dr, Model, DynareOptions)); measure = tmp(mf1,:); else measure = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); diff --git a/src/nonlinear_kalman_filter.m b/src/nonlinear_kalman_filter.m index 322eb6e..594c751 100644 --- a/src/nonlinear_kalman_filter.m +++ b/src/nonlinear_kalman_filter.m @@ -56,6 +56,7 @@ end if ReducedForm.use_k_order_solver dr = ReducedForm.dr; + [udr] = folded_to_unfolded_dr(dr, Model, DynareOptions); else % Set local state space model (first-order approximation). ghx = ReducedForm.ghx; @@ -116,7 +117,7 @@ for t=1:sample_size epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:); yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions,udr); else tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); end diff --git a/src/online_auxiliary_filter.m b/src/online_auxiliary_filter.m index 91c51f5..a7b79fa 100644 --- a/src/online_auxiliary_filter.m +++ b/src/online_auxiliary_filter.m @@ -125,6 +125,7 @@ for t=1:sample_size % Set local state space model (second-order approximation). if ReducedForm.use_k_order_solver dr = ReducedForm.dr; + [udr] = folded_to_unfolded_dr(dr, Model, DynareOptions); else constant = ReducedForm.constant; % Set local state space model (first-order approximation). @@ -138,7 +139,7 @@ for t=1:sample_size % particle likelihood contribution yhat = bsxfun(@minus, StateVectors(:,i), state_variables_steady_state); if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions, udr); else if pruning yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state); @@ -194,7 +195,7 @@ for t=1:sample_size % compute particles likelihood contribution yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state); if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); else if pruning yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state); diff --git a/src/sequential_importance_particle_filter.m b/src/sequential_importance_particle_filter.m index 975ac12..365bbc2 100644 --- a/src/sequential_importance_particle_filter.m +++ b/src/sequential_importance_particle_filter.m @@ -51,6 +51,7 @@ end if ReducedForm.use_k_order_solver dr = ReducedForm.dr; + [udr] = folded_to_unfolded_dr(dr, Model, DynareOptions); else % Set local state space model (first order approximation). ghx = ReducedForm.ghx; @@ -106,7 +107,7 @@ 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); else if ReducedForm.use_k_order_solver - tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions); + tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); else tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); end -- GitLab