From 7721ea041c5cbe66f38cacf5d9a2b0d1fed5020b Mon Sep 17 00:00:00 2001 From: Normann Rion Date: Mon, 19 Sep 2022 15:43:10 +0200 Subject: [PATCH] Amends particle filters to use the local_state_space_iteration_3 mex --- .../src/auxiliary_initialization.m | 20 ++++- .../src/auxiliary_particle_filter.m | 33 +++++-- .../src/conditional_filter_proposal.m | 19 +++- .../src/gaussian_filter_bank.m | 19 +++- .../src/gaussian_mixture_filter_bank.m | 19 +++- .../src/measurement_equations.m | 17 +++- .../src/nonlinear_kalman_filter.m | 17 +++- .../src/online_auxiliary_filter.m | 90 +++++++++++++++++-- .../sequential_importance_particle_filter.m | 47 ++++++++-- 9 files changed, 254 insertions(+), 27 deletions(-) diff --git a/matlab/nonlinear-filters/src/auxiliary_initialization.m b/matlab/nonlinear-filters/src/auxiliary_initialization.m index 8355e75ab..3e279117b 100644 --- a/matlab/nonlinear-filters/src/auxiliary_initialization.m +++ b/matlab/nonlinear-filters/src/auxiliary_initialization.m @@ -2,7 +2,7 @@ function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,Par % Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling. -% Copyright © 2011-2017 Dynare Team +% Copyright © 2011-2022 Dynare Team % % This file is part of Dynare (particles module). % @@ -45,6 +45,8 @@ if isempty(init_flag) init_flag = 1; end +order = DynareOptions.order; + % Set local state space model (first order approximation). ghx = ReducedForm.ghx; ghu = ReducedForm.ghu; @@ -52,6 +54,14 @@ ghu = ReducedForm.ghu; ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; +if (order == 3) + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; +end % Get covariance matrices Q = ReducedForm.Q; @@ -87,7 +97,13 @@ yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); % yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state); % [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); %else -tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); +if (order == 2) + tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); +elseif (order == 3) + tmp = local_state_space_iteration_3(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,ThreadsOptions.local_state_space_iteration_3); +else + error('Orders > 3 not allowed'); +end %end PredictedObservedMean = weights*(tmp(mf1,:)'); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); diff --git a/matlab/nonlinear-filters/src/auxiliary_particle_filter.m b/matlab/nonlinear-filters/src/auxiliary_particle_filter.m index 2d3b00112..b84eca837 100644 --- a/matlab/nonlinear-filters/src/auxiliary_particle_filter.m +++ b/matlab/nonlinear-filters/src/auxiliary_particle_filter.m @@ -54,6 +54,15 @@ else ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end end % Get covariance matrices @@ -83,8 +92,14 @@ if pruning StateVectors_ = StateVectors; state_variables_steady_state_ = state_variables_steady_state; mf0_ = mf0; + elseif order == 3 + StateVectors_ = repmat(StateVectors,2,1); + state_variables_steady_state_ = repmat(state_variables_steady_state,2,1); + mf0_ = repmat(mf0,1,2); + mask = number_of_state_variables+1:2*number_of_state_variables; + mf0_(mask) = mf0_(mask)+size(ghx,1); else - error('Pruning is not available for orders > 2'); + error('Pruning is not available for orders > 3'); end end @@ -94,8 +109,10 @@ for t=1:sample_size yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_); if order == 2 [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_3); else - error('Pruning is not available for orders > 2'); + error('Pruning is not available for orders > 3'); end else if ReducedForm.use_k_order_solver @@ -103,8 +120,10 @@ for t=1:sample_size else if order == 2 tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,ThreadsOptions.local_state_space_iteration_3); else - error('Order > 2: use_k_order_solver should be set to true'); + error('Order > 3: use_k_order_solver should be set to true'); end end end @@ -122,8 +141,10 @@ for t=1:sample_size if pruning if order == 2 [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_3); else - error('Pruning is not available for orders > 2'); + error('Pruning is not available for orders > 3'); end StateVectors_ = tmp_(mf0_,:); else @@ -132,8 +153,10 @@ for t=1:sample_size else if order == 2 tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ThreadsOptions.local_state_space_iteration_3); else - error('Order > 2: use_k_order_solver should be set to true'); + error('Order > 3: use_k_order_solver should be set to true'); end end end diff --git a/matlab/nonlinear-filters/src/conditional_filter_proposal.m b/matlab/nonlinear-filters/src/conditional_filter_proposal.m index 47e53b5ac..45d7886ae 100644 --- a/matlab/nonlinear-filters/src/conditional_filter_proposal.m +++ b/matlab/nonlinear-filters/src/conditional_filter_proposal.m @@ -41,6 +41,8 @@ function [ProposalStateVector, Weights, flag] = conditional_filter_proposal(Redu flag = false; +order = DynareOptions.order; + if ReducedForm.use_k_order_solver dr = ReducedForm.dr; udr = ReducedForm.udr; @@ -52,6 +54,15 @@ else ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end end constant = ReducedForm.constant; @@ -82,7 +93,13 @@ 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, udr); else - tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ThreadsOptions.local_state_space_iteration_3); + else + error('Order > 3: use_k_order_solver should be set to true'); + end end PredictedStateMean = tmp(mf0,:)*weights; diff --git a/matlab/nonlinear-filters/src/gaussian_filter_bank.m b/matlab/nonlinear-filters/src/gaussian_filter_bank.m index 6108ab21d..d7a649fa4 100644 --- a/matlab/nonlinear-filters/src/gaussian_filter_bank.m +++ b/matlab/nonlinear-filters/src/gaussian_filter_bank.m @@ -39,6 +39,8 @@ function [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +order = DynareOptions.order; + if ReducedForm.use_k_order_solver dr = ReducedForm.dr; udr = ReducedForm.udr; @@ -50,6 +52,15 @@ else ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end end constant = ReducedForm.constant; @@ -84,7 +95,13 @@ 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, udr); else - tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ThreadsOptions.local_state_space_iteration_3); + else + error('Order > 3: use_k_order_solver should be set to true'); + end end PredictedStateMean = tmp(mf0,:)*weights; diff --git a/matlab/nonlinear-filters/src/gaussian_mixture_filter_bank.m b/matlab/nonlinear-filters/src/gaussian_mixture_filter_bank.m index 88d3628b8..69e3eb613 100644 --- a/matlab/nonlinear-filters/src/gaussian_mixture_filter_bank.m +++ b/matlab/nonlinear-filters/src/gaussian_mixture_filter_bank.m @@ -41,6 +41,8 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +order = DynareOptions.order; + if ReducedForm.use_k_order_solver dr = ReducedForm.dr; udr = ReducedForm.udr; @@ -52,6 +54,15 @@ else ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end end constant = ReducedForm.constant; @@ -80,7 +91,13 @@ 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, udr); else - tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ThreadsOptions.local_state_space_iteration_3); + else + error('Order > 3: use_k_order_solver should be set to true'); + end end PredictedStateMean = tmp(mf0,:)*weights3; PredictedObservedMean = tmp(mf1,:)*weights3; diff --git a/matlab/nonlinear-filters/src/measurement_equations.m b/matlab/nonlinear-filters/src/measurement_equations.m index 4e11dd064..477cc1782 100644 --- a/matlab/nonlinear-filters/src/measurement_equations.m +++ b/matlab/nonlinear-filters/src/measurement_equations.m @@ -17,6 +17,7 @@ function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +order = DynareOptions.order; mf1 = ReducedForm.mf1; if ReducedForm.use_k_order_solver dr = ReducedForm.dr; @@ -27,6 +28,14 @@ else ghxx = ReducedForm.ghxx(mf1,:); ghuu = ReducedForm.ghuu(mf1,:); ghxu = ReducedForm.ghxu(mf1,:); + if order == 3 + ghxxx = ReducedForm.ghxxx(mf1,:); + ghuuu = ReducedForm.ghuuu(mf1,:); + ghxxu = ReducedForm.ghxxu(mf1,:); + ghxuu = ReducedForm.ghxuu(mf1,:); + ghxss = ReducedForm.ghxss(mf1,:); + ghuss = ReducedForm.ghuss(mf1,:); + end end constant = ReducedForm.constant(mf1,:); state_variables_steady_state = ReducedForm.state_variables_steady_state; @@ -36,5 +45,11 @@ if ReducedForm.use_k_order_solver tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions, udr); 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); + if order == 2 + 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); + elseif order == 3 + measure = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ThreadsOptions.local_state_space_iteration_3); + else + error('Order > 3: use_k_order_solver should be set to true'); + end end diff --git a/matlab/nonlinear-filters/src/nonlinear_kalman_filter.m b/matlab/nonlinear-filters/src/nonlinear_kalman_filter.m index 9decb3225..d46c0cd7e 100644 --- a/matlab/nonlinear-filters/src/nonlinear_kalman_filter.m +++ b/matlab/nonlinear-filters/src/nonlinear_kalman_filter.m @@ -54,6 +54,8 @@ if isempty(start) start = 1; end +order = DynareOptions.order; + if ReducedForm.use_k_order_solver dr = ReducedForm.dr; udr = ReducedForm.udr; @@ -65,6 +67,15 @@ else ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end end constant = ReducedForm.constant; @@ -119,7 +130,11 @@ for t=1:sample_size if ReducedForm.use_k_order_solver 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); + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ThreadsOptions.local_state_space_iteration_3); + end end PredictedStateMean = tmp(mf0,:)*weights ; PredictedObservedMean = tmp(mf1,:)*weights; diff --git a/matlab/nonlinear-filters/src/online_auxiliary_filter.m b/matlab/nonlinear-filters/src/online_auxiliary_filter.m index eefd163bb..8d0745514 100644 --- a/matlab/nonlinear-filters/src/online_auxiliary_filter.m +++ b/matlab/nonlinear-filters/src/online_auxiliary_filter.m @@ -49,6 +49,7 @@ bounds = prior_bounds(BayesInfo, DynareOptions.prior_trunc); % Reset bounds as l % initialization of state particles [~, Model, DynareOptions, DynareResults, ReducedForm] = solve_model_for_online_filter(true, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults); +order = DynareOptions.order; mf0 = ReducedForm.mf0; mf1 = ReducedForm.mf1; number_of_particles = DynareOptions.particle.number_of_particles; @@ -64,8 +65,14 @@ StateVectorMean = ReducedForm.StateVectorMean; StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)'; state_variance_rank = size(StateVectorVarianceSquareRoot,2); StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); -if DynareOptions.order<3 && pruning - StateVectors_ = StateVectors; +if pruning + if order == 2 + StateVectors_ = StateVectors; + elseif order == 3 + StateVectors_ = repmat(StateVectors,2,1); + else + error('Pruning is not available for orders > 3'); + end end % parameters for the Liu & West filter @@ -135,6 +142,25 @@ for t=1:sample_size ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end + if pruning + if order == 2 + state_variables_steady_state_ = state_variables_steady_state; + elseif order == 3 + state_variables_steady_state_ = repmat(state_variables_steady_state,2,1); + else + error('Pruning is not available for orders > 3'); + end + end + end % particle likelihood contribution yhat = bsxfun(@minus, StateVectors(:,i), state_variables_steady_state); @@ -142,10 +168,22 @@ for t=1:sample_size 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); - [tmp, ~] = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); + yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state_); + if order == 2 + [tmp, ~] = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + [tmp, ~] = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_3); + else + error('Pruning is not available for orders > 3'); + end else - tmp = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); + if order == 2 + tmp = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, DynareOptions.threads.local_state_space_iteration_3); + else + error('Order > 3: use_k_order_solver should be set to true'); + end end end PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:)); @@ -192,6 +230,28 @@ for t=1:sample_size ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if (order == 3) + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end + if pruning + if order == 2 + state_variables_steady_state_ = state_variables_steady_state; + mf0_ = mf0; + elseif order == 3 + state_variables_steady_state_ = repmat(state_variables_steady_state,2,1); + mf0_ = repmat(mf0,1,2); + mask = number_of_state_variables+1:2*number_of_state_variables; + mf0_(mask) = mf0_(mask)+size(ghx,1); + else + error('Pruning is not available for orders > 3'); + end + end end % Get covariance matrices and structural shocks epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations, 1); @@ -201,11 +261,23 @@ for t=1:sample_size tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr); else if pruning - yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state); - [tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); - StateVectors_(:,i) = tmp_(mf0,:); + yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state_); + if order == 2 + [tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_3); + else + error('Pruning is not available for orders > 3'); + end + StateVectors_(:,i) = tmp_(mf0_,:); else - tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); + if order == 2 + tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, DynareOptions.threads.local_state_space_iteration_3); + else + error('Order > 3: use_k_order_solver should be set to true'); + end end end StateVectors(:,i) = tmp(mf0,:); diff --git a/matlab/nonlinear-filters/src/sequential_importance_particle_filter.m b/matlab/nonlinear-filters/src/sequential_importance_particle_filter.m index ab841009c..4c9c8313e 100644 --- a/matlab/nonlinear-filters/src/sequential_importance_particle_filter.m +++ b/matlab/nonlinear-filters/src/sequential_importance_particle_filter.m @@ -37,6 +37,8 @@ steadystate = ReducedForm.steadystate; constant = ReducedForm.constant; state_variables_steady_state = ReducedForm.state_variables_steady_state; +order = DynareOptions.order; + % Set persistent variables (if needed). if isempty(init_flag) mf0 = ReducedForm.mf0; @@ -61,6 +63,15 @@ else ghxx = ReducedForm.ghxx; ghuu = ReducedForm.ghuu; ghxu = ReducedForm.ghxu; + if order == 3 + % Set local state space model (third order approximation). + ghxxx = ReducedForm.ghxxx; + ghuuu = ReducedForm.ghuuu; + ghxxu = ReducedForm.ghxxu; + ghxuu = ReducedForm.ghxuu; + ghxss = ReducedForm.ghxss; + ghuss = ReducedForm.ghuss; + end end % Get covariance matrices. @@ -95,7 +106,19 @@ set_dynare_seed('default'); weights = ones(1,number_of_particles)/number_of_particles ; StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean); if pruning - StateVectors_ = StateVectors; + if order == 2 + StateVectors_ = StateVectors; + state_variables_steady_state_ = state_variables_steady_state; + mf0_ = mf0; + elseif order == 3 + StateVectors_ = repmat(StateVectors,2,1); + state_variables_steady_state_ = repmat(state_variables_steady_state,2,1); + mf0_ = repmat(mf0,1,2); + mask = number_of_state_variables+1:2*number_of_state_variables; + mf0_(mask) = mf0_(mask)+size(ghx,1); + else + error('Pruning is not available for orders > 3'); + end end % Loop over observations @@ -103,13 +126,25 @@ for t=1:sample_size yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles); if pruning - yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state); - [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); + yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_); + if order == 2 + [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + [tmp, tmp_] = local_state_space_iteration_3(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_3); + else + error('Pruning is not available for orders > 3'); + end else if ReducedForm.use_k_order_solver 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); + if order == 2 + tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); + elseif order == 3 + tmp = local_state_space_iteration_3(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,ThreadsOptions.local_state_space_iteration_3); + else + error('Order > 3: use_k_order_solver should be set to true'); + end end end %PredictedObservedMean = tmp(mf1,:)*transpose(weights); @@ -129,7 +164,7 @@ for t=1:sample_size weights = wtilde/sum(wtilde); if (ParticleOptions.resampling.status.generic && neff(weights)