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

Ref. dynare#1673

matlab/nonlinear-filters/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,20 +123,30 @@ 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).
+            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 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.
             z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ;
@@ -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,16 +177,25 @@ 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).
+                    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 ReducedForm.use_k_order_solver
+                        tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions);
+                    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);
@@ -184,6 +203,7 @@ for t=1:sample_size
                         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,:));
                     wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1)));

matlab/nonlinear-filters/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,10 +182,10 @@ 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);
@@ -189,7 +194,7 @@ if setinitialcondition
         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!')