diff --git a/src/online_auxiliary_filter.m b/src/online_auxiliary_filter.m
index 18f762ec92aea71f05d5d8a6f9755cf42cfbf1ae..6acadc6080c6a7512934ded68bf6cfebac3b1d70 100644
--- a/src/online_auxiliary_filter.m
+++ b/src/online_auxiliary_filter.m
@@ -64,7 +64,7 @@ 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 pruning
+if DynareOptions.order<3 && pruning
StateVectors_ = StateVectors;
end
@@ -123,19 +123,29 @@ for t=1:sample_size
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set local state space model (second-order approximation).
- constant = ReducedForm.constant;
- ghx = ReducedForm.ghx;
- ghu = ReducedForm.ghu;
- ghxx = ReducedForm.ghxx;
- ghuu = ReducedForm.ghuu;
- ghxu = ReducedForm.ghxu;
+ if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ else
+ constant = ReducedForm.constant;
+ % 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
% particle likelihood contribution
- yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
- 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 ReducedForm.use_k_order_solver
+ tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions);
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 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);
+ 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);
+ end
end
PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
% Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails.
@@ -148,7 +158,7 @@ for t=1:sample_size
indx = resample(0, tau_tilde', DynareOptions.particle);
StateVectors = StateVectors(:,indx);
xparam = fore_xparam(:,indx);
- if pruning
+ if DynareOptions.order>=3 && pruning
StateVectors_ = StateVectors_(:,indx);
end
w_stage1 = weights(indx)./tau_tilde(indx);
@@ -167,22 +177,32 @@ for t=1:sample_size
steadystate = ReducedForm.steadystate;
state_variables_steady_state = ReducedForm.state_variables_steady_state;
% Set local state space model (second order approximation).
- constant = ReducedForm.constant;
- ghx = ReducedForm.ghx;
- ghu = ReducedForm.ghu;
- ghxx = ReducedForm.ghxx;
- ghuu = ReducedForm.ghuu;
- ghxu = ReducedForm.ghxu;
+ if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ else
+ constant = ReducedForm.constant;
+ % 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 and structural shocks
epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations, 1);
% compute particles likelihood contribution
yhat = bsxfun(@minus,StateVectors(:,i), state_variables_steady_state);
- 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,:);
+ 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, DynareOptions.threads.local_state_space_iteration_2);
+ 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,:);
+ else
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2);
+ end
end
StateVectors(:,i) = tmp(mf0,:);
PredictionError = bsxfun(@minus,Y(t,:)', tmp(mf1,:));
diff --git a/src/solve_model_for_online_filter.m b/src/solve_model_for_online_filter.m
index 4df219d8728659725e695ecbf72b5a6be6572095..332378d5a930aafe81a154e7ffa325d2545842ec 100644
--- a/src/solve_model_for_online_filter.m
+++ b/src/solve_model_for_online_filter.m
@@ -155,12 +155,17 @@ if nargout>4
ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
- if DynareOptions.order>1
+ if DynareOptions.order==2
+ ReducedForm.use_k_order_solver = false;
ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
+ elseif DynareOptions.order>=3
+ ReducedForm.use_k_order_solver = true;
+ ReducedForm.dr = dr;
else
+ ReducedForm.use_k_order_solver = false;
ReducedForm.ghxx = zeros(size(restrict_variables_idx,1),size(dr.kstate,2));
ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),size(dr.ghu,2));
ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),size(dr.ghx,2));
@@ -177,19 +182,19 @@ end
if setinitialcondition
switch DynareOptions.particle.initialization
case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
- StateVectorMean = ReducedForm.constant(mf0);
+ StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
- StateVectorMean = ReducedForm.constant(mf0);
+ StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
old_DynareOptionsperiods = DynareOptions.periods;
DynareOptions.periods = 5000;
- y_ = simult(oo_.steady_state, dr,Model,DynareOptions,DynareResults);
+ y_ = simult(oo_.steady_state, dr, Model, DynareOptions, DynareResults);
y_ = y_(state_variables_idx,2001:5000);
StateVectorVariance = cov(y_');
DynareOptions.periods = old_DynareOptionsperiods;
clear('old_DynareOptionsperiods','y_');
case 3% Initial state vector covariance is a diagonal matrix.
- StateVectorMean = ReducedForm.constant(mf0);
+ StateVectorMean = ReducedForm.state_variables_steady_state;%.constant(mf0);
StateVectorVariance = DynareOptions.particle.initial_state_prior_std*eye(number_of_state_variables);
otherwise
error('Unknown initialization option!')