diff --git a/src/online_auxiliary_filter.m b/src/online_auxiliary_filter.m index 18f762ec9..6acadc608 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; - % 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,:)); diff --git a/src/solve_model_for_online_filter.m b/src/solve_model_for_online_filter.m index 4df219d87..332378d5a 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!')