stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Amends particle filters to use the local_state_space_iteration_3 mex

rm-particles^2
Normann Rion 2022-09-19 15:43:10 +02:00
parent 25401abeb8
commit 80d48b7745
9 changed files with 254 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/auxiliary_initialization.m
View File

@ -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,:));

33
src/auxiliary_particle_filter.m
View File

@ -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

19
src/conditional_filter_proposal.m
View File

@ -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;

19
src/gaussian_filter_bank.m
View File

@ -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;

19
src/gaussian_mixture_filter_bank.m
View File

@ -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;

17
src/measurement_equations.m
View File

@ -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

17
src/nonlinear_kalman_filter.m
View File

@ -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;

90
src/online_auxiliary_filter.m
View File

@ -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,:);

47
src/sequential_importance_particle_filter.m
View File

@ -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