diff --git a/matlab/nonlinear-filters/src/auxiliary_initialization.m b/matlab/nonlinear-filters/src/auxiliary_initialization.m
index 8355e75ab8c82185bd06d7b065fb0952368c64ae..3e279117b760d3188b636ea52dd1332bfe832f9a 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 2caa3e91914da4c290cb6e894b6ed0bf2b6a86bb..b84eca837ad672b81a30e67266a100dced5ca931 100644
--- a/matlab/nonlinear-filters/src/auxiliary_particle_filter.m
+++ b/matlab/nonlinear-filters/src/auxiliary_particle_filter.m
@@ -24,6 +24,8 @@ function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptio
if isempty(start)
start = 1;
end
+% Get perturbation order
+order = DynareOptions.order;
% Set flag for prunning
pruning = ParticleOptions.pruning;
@@ -36,6 +38,7 @@ state_variables_steady_state = ReducedForm.state_variables_steady_state;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
sample_size = size(Y,2);
+number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
number_of_particles = ParticleOptions.number_of_particles;
@@ -51,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
@@ -76,33 +88,42 @@ weights = ones(1,number_of_particles)/number_of_particles ;
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
%StateVectors = bsxfun(@plus,zeros(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
-nodes = zeros(1,number_of_structural_innovations) ;
-nodes_weights = ones(number_of_structural_innovations,1) ;
-
for t=1:sample_size
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
if pruning
- yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
- tmp = 0 ;
- tmp_ = 0 ;
- for i=1:size(nodes)
- [tmp1, tmp1_] = local_state_space_iteration_2(yhat,nodes(i,:)'*ones(1,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
- tmp = tmp + nodes_weights(i)*tmp1 ;
- tmp_ = tmp_ + nodes_weights(i)*tmp1_ ;
+ 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 > 3');
end
else
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, udr);
- end
+ tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations,number_of_particles), dr, Model, DynareOptions, udr);
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 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 > 3: use_k_order_solver should be set to true');
end
end
end
@@ -118,13 +139,25 @@ for t=1:sample_size
weights_stage_1 = weights(indx)./tau_tilde(indx) ;
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
if pruning
- [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,:);
+ 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
+ StateVectors_ = tmp_(mf0_,:);
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
StateVectors = tmp(mf0,:);
@@ -135,9 +168,8 @@ for t=1:sample_size
if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
if pruning
temp = resample([StateVectors' StateVectors_'],weights',ParticleOptions);
- number_of_state_variables=size(StateVectors,1);
StateVectors = temp(:,1:number_of_state_variables)';
- StateVectors_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)';
+ StateVectors_ = temp(:,number_of_state_variables+1:end)';
else
StateVectors = resample(StateVectors',weights',ParticleOptions)';
end
diff --git a/matlab/nonlinear-filters/src/conditional_filter_proposal.m b/matlab/nonlinear-filters/src/conditional_filter_proposal.m
index 47e53b5ac6eb73fa16fe00ff0b986e1092075ca0..45d7886ae281aa9288fe21ceefce746b558c97b1 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 6108ab21ddaf9343ac78d9038a58e0c447e5c32c..d7a649fa417bb9a09d2e0832de34b9c105db998c 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 <https://www.gnu.org/licenses/>.
+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 88d3628b89f67d42a7948c3107f6f10535a8cf42..69e3eb6137722237f558db0480eb57ad3ebf38d0 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 <https://www.gnu.org/licenses/>.
+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 4e11dd064e079419a38fb9bd5003bd7a44644dab..477cc1782cc451c45caa353717387b5233fd376f 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 <https://www.gnu.org/licenses/>.
+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 9decb32256a3beedb79549cf78d0604a487b2e11..d46c0cd7eab3697248e3c7964c856c123ebd76ac 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 bfb4383c57c84122dba28afb63619890ab5e111e..8d0745514ada4583a559a9c7d5de720583302c8b 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,:));
@@ -159,7 +197,7 @@ for t=1:sample_size
indx = resample(0, tau_tilde', DynareOptions.particle);
StateVectors = StateVectors(:,indx);
xparam = fore_xparam(:,indx);
- if DynareOptions.order>=3 && pruning
+ if pruning
StateVectors_ = StateVectors_(:,indx);
end
w_stage1 = weights(indx)./tau_tilde(indx);
@@ -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 ab841009c0c218af96b644eec8152af5edfc4510..4c9c8313e60393a02edc997653f436fea13d1f6c 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)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
if pruning
- temp = resample([tmp(mf0,:)' tmp_(mf0,:)'],weights',ParticleOptions);
+ temp = resample([tmp(mf0,:)' tmp_(mf0_,:)'],weights',ParticleOptions);
StateVectors = temp(:,1:number_of_state_variables)';
StateVectors_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)';
else
@@ -139,7 +174,7 @@ for t=1:sample_size
elseif ParticleOptions.resampling.status.none
StateVectors = tmp(mf0,:);
if pruning
- StateVectors_ = tmp_(mf0,:);
+ StateVectors_ = tmp_(mf0_,:);
end
end
end