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

Allow k order approximation in online filter (mode_compute=11). 

Ref. dynare#1673
rm-particles^2
Stéphane Adjemian (Charybdis) 2020-02-15 16:50:32 +01:00
parent 2ed534e52e
commit cef2015162
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
2 changed files with 55 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
src/online_auxiliary_filter.m
View File

@ -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;
% 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);
if ReducedForm.use_k_order_solver
dr = ReducedForm.dr;
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);
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 ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, 1), dr, Model, DynareOptions);
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);
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,:));

15
src/solve_model_for_online_filter.m
View File

@ -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!')