Allow k order approximation in Gaussian Mixture Filter (gmf).

Ref. dynare#1673

src/gaussian_mixture_densities.m

-function  IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
-                                                  StateMuPost,StateSqrtPPost,StateWeightsPost,...
-                                                  StateParticles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions)
-%
+function  IncrementalWeights = gaussian_mixture_densities(obs, StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
+                                                      StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
+                                                      ReducedForm, ThreadsOptions, DynareOptions, Model)
+
 % Elements to calculate the importance sampling ratio
 %
 % INPUTS
@@ -21,7 +21,8 @@ function  IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,State
 %
 % NOTES
 %   The vector "lik" is used to evaluate the jacobian of the likelihood.
-% Copyright (C) 2009-2017 Dynare Team
+
+% Copyright (C) 2009-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,19 +40,16 @@ function  IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,State
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 % Compute the density of particles under the prior distribution
-[ras,ras,prior] = probability(StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateParticles) ;
-prior = prior' ;
+[~, ~, prior] = probability(StateMuPrior, StateSqrtPPrior, StateWeightsPrior, StateParticles);
+prior = prior';
+
 % Compute the density of particles under the proposal distribution
-[ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ;
-proposal = proposal' ;
+[~, ~, proposal] = probability(StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles);
+proposal = proposal';
+
 % Compute the density of the current observation conditionally to each particle
-yt_t_1_i = measurement_equations(StateParticles,ReducedForm,ThreadsOptions) ;
-%eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
-%yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
-%tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
-%Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
-%sqr_det = sqrt(det(Pyy)) ;
-%foo = (eta_t_i/Pyy).*eta_t_i ;
-%likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;
-likelihood = probability2(obs,sqrt(H),yt_t_1_i) ;
-IncrementalWeights = likelihood.*prior./proposal ;
+yt_t_1_i = measurement_equations(StateParticles, ReducedForm, ThreadsOptions, DynareOptions, Model);
+
+% likelihood
+likelihood = probability2(obs, sqrt(H), yt_t_1_i);
+IncrementalWeights = likelihood.*prior./proposal;

src/gaussian_mixture_filter_bank.m

 function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPPost,StateWeightsPost] =...
-    gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,...
-                                 StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,...
-                                 ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,...
-                                 H,H_lower_triangular_cholesky,normfactO,ParticleOptions,ThreadsOptions)
-%
+    gaussian_mixture_filter_bank(ReducedForm, obs, StateMu, StateSqrtP, StateWeights, ...
+                                 StructuralShocksMu, StructuralShocksSqrtP, StructuralShocksWeights, ...
+                                 ObservationShocksWeights, H, H_lower_triangular_cholesky, normfactO, ...
+                                 ParticleOptions, ThreadsOptions, DynareOptions, Model)
+
 % Computes the proposal with a gaussian approximation for importance
 % sampling
 % This proposal is a gaussian distribution calculated à la Kalman
@@ -23,7 +23,8 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP
 %
 % NOTES
 %   The vector "lik" is used to evaluate the jacobian of the likelihood.
-% Copyright (C) 2009-2017 Dynare Team
+
+% Copyright (C) 2009-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -40,86 +41,73 @@ function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPP
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-
-
-persistent init_flag2 mf0 mf1 %nodes3 weights3 weights_c3
-persistent number_of_state_variables number_of_observed_variables
-persistent number_of_structural_innovations
-
-% 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;
-
-if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ...
-        any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ...
-        any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4))
-    ghx
-    ghu
-    ghxx
-    ghuu
-    ghxu
+if ReducedForm.use_k_order_solver
+    dr = ReducedForm.dr;
+else
+    % 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
 
 constant = ReducedForm.constant;
 state_variables_steady_state = ReducedForm.state_variables_steady_state;
 
-% Set persistent variables.
-if isempty(init_flag2)
-    mf0 = ReducedForm.mf0;
-    mf1 = ReducedForm.mf1;
-    number_of_state_variables = length(mf0);
-    number_of_observed_variables = length(mf1);
-    number_of_structural_innovations = length(ReducedForm.Q);
-    init_flag2 = 1;
-end
+mf0 = ReducedForm.mf0;
+mf1 = ReducedForm.mf1;
+number_of_state_variables = length(mf0);
+number_of_observed_variables = length(mf1);
+number_of_structural_innovations = length(ReducedForm.Q);
 
-numb = number_of_state_variables+number_of_structural_innovations ;
+numb = number_of_state_variables+number_of_structural_innovations;
 
 if ParticleOptions.proposal_approximation.cubature
-    [nodes3,weights3] = spherical_radial_sigma_points(numb);
+    [nodes3, weights3] = spherical_radial_sigma_points(numb);
     weights_c3 = weights3;
 elseif ParticleOptions.proposal_approximation.unscented
-    [nodes3,weights3,weights_c3] = unscented_sigma_points(numb,ParticleOptions);
+    [nodes3, weights3, weights_c3] = unscented_sigma_points(numb, ParticleOptions);
 else
-    error('Estimation: This approximation for the proposal is not implemented or unknown!')
+    error('This approximation for the proposal is unknown!')
 end
 
-epsilon =  bsxfun(@plus,StructuralShocksSqrtP*nodes3(:,number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations)',StructuralShocksMu) ;
-StateVectors = bsxfun(@plus,StateSqrtP*nodes3(:,1:number_of_state_variables)',StateMu);
-yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
-tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
+epsilon = bsxfun(@plus, StructuralShocksSqrtP*nodes3(:,number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations)', StructuralShocksMu);
+StateVectors = bsxfun(@plus, StateSqrtP*nodes3(:,1:number_of_state_variables)', StateMu);
+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);
+else
+    tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
+end
 PredictedStateMean = tmp(mf0,:)*weights3;
 PredictedObservedMean = tmp(mf1,:)*weights3;
 
 if ParticleOptions.proposal_approximation.cubature
-    PredictedStateMean = sum(PredictedStateMean,2);
-    PredictedObservedMean = sum(PredictedObservedMean,2);
-    dState = (bsxfun(@minus,tmp(mf0,:),PredictedStateMean)').*sqrt(weights3);
-    dObserved = (bsxfun(@minus,tmp(mf1,:),PredictedObservedMean)').*sqrt(weights3);
+    PredictedStateMean = sum(PredictedStateMean, 2);
+    PredictedObservedMean = sum(PredictedObservedMean, 2);
+    dState = (bsxfun(@minus, tmp(mf0,:), PredictedStateMean)').*sqrt(weights3);
+    dObserved = (bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)').*sqrt(weights3);
     PredictedStateVariance = dState'*dState;
-    big_mat = [dObserved  dState ; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ];
-    [mat1,mat] = qr2(big_mat,0);
+    big_mat = [dObserved,  dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables, number_of_state_variables)];
+    [~, mat] = qr2(big_mat, 0);
     mat = mat';
-    clear('mat1');
-    PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables);
-    CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
-    StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
+    PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables, 1:number_of_observed_variables);
+    CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables), 1:number_of_observed_variables);
+    StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables), number_of_observed_variables+(1:number_of_state_variables));
     iPredictedObservedVarianceSquareRoot = inv(PredictedObservedVarianceSquareRoot);
     iPredictedObservedVariance = iPredictedObservedVarianceSquareRoot'*iPredictedObservedVarianceSquareRoot;
     sqrdet = 1/sqrt(det(iPredictedObservedVariance));
     PredictionError = obs - PredictedObservedMean;
     StateVectorMean = PredictedStateMean + CovarianceObservedStateSquareRoot*iPredictedObservedVarianceSquareRoot*PredictionError;
 else
-    dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean);
-    dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
+    dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean);
+    dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean);
     PredictedStateVariance = dState*diag(weights_c3)*dState';
     PredictedObservedVariance = dObserved*diag(weights_c3)*dObserved' + H;
     PredictedStateAndObservedCovariance = dState*diag(weights_c3)*dObserved';
-    sqrdet = sqrt(det(PredictedObservedVariance)) ;
+    sqrdet = sqrt(det(PredictedObservedVariance));
     iPredictedObservedVariance = inv(PredictedObservedVariance);
     PredictionError = obs - PredictedObservedMean;
     KalmanFilterGain = PredictedStateAndObservedCovariance*iPredictedObservedVariance;
@@ -130,9 +118,9 @@ else
 end
 
 data_lik_GM_g = exp(-0.5*PredictionError'*iPredictedObservedVariance*PredictionError)/abs(normfactO*sqrdet) + 1e-99;
-StateMuPrior = PredictedStateMean ;
+StateMuPrior = PredictedStateMean;
 StateSqrtPPrior = reduced_rank_cholesky(PredictedStateVariance)';
 StateWeightsPrior = StateWeights*StructuralShocksWeights;
 StateMuPost = StateVectorMean;
 StateSqrtPPost = StateVectorVarianceSquareRoot;
-StateWeightsPost = StateWeightsPrior*ObservationShocksWeights*data_lik_GM_g ;
+StateWeightsPost = StateWeightsPrior*ObservationShocksWeights*data_lik_GM_g;
\ No newline at end of file