diff --git a/.gitmodules b/.gitmodules
index 2cc13d27ea15221701b4f6fa32ee74a9dc23eec8..c0399c3cf74e3ad4e746a7dd34e42b10fcc3bcd4 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -10,9 +10,6 @@
[submodule "matlab/utilities/tests"]
path = matlab/utilities/tests
url = ../../Dynare/m-unit-tests.git
-[submodule "matlab/particles"]
- path = matlab/particles
- url = ../../Dynare/particles.git
[submodule "matlab/modules/dseries"]
path = matlab/modules/dseries
url = ../../Dynare/dseries.git
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index af43656e73a86e972fd4e9b82274da74e9e6d1a1..11d4a0831ac98093a46146dda83d3c286c69fdbf 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -53,7 +53,7 @@ p = {'/distributions/' ; ...
'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ...
'/../contrib/jsonlab/' ; ...
'/parallel/' ; ...
- '/particles/src' ; ...
+ '/nonlinear-filters/' ; ...
'/gsa/' ; ...
'/ep/' ; ...
'/backward/' ; ...
diff --git a/matlab/nonlinear-filters/auxiliary_initialization.m b/matlab/nonlinear-filters/auxiliary_initialization.m
new file mode 100644
index 0000000000000000000000000000000000000000..e05741052a346b0e512e782158dc1b90dcc90345
--- /dev/null
+++ b/matlab/nonlinear-filters/auxiliary_initialization.m
@@ -0,0 +1,117 @@
+function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
+
+% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
+
+% Copyright © 2011-2022 Dynare Team
+%
+% This file is part of Dynare (particles module).
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare particles module is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+persistent init_flag mf0 mf1 number_of_particles
+persistent number_of_observed_variables number_of_structural_innovations
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+% Set flag for prunning
+%pruning = ParticleOptions.pruning;
+
+% Get steady state and mean.
+%steadystate = ReducedForm.steadystate;
+constant = ReducedForm.constant;
+ss = ReducedForm.ys;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+% Set persistent variables.
+if isempty(init_flag)
+ mf0 = ReducedForm.mf0;
+ mf1 = ReducedForm.mf1;
+ number_of_observed_variables = length(mf1);
+ number_of_structural_innovations = length(ReducedForm.Q);
+ number_of_particles = ParticleOptions.number_of_particles;
+ init_flag = 1;
+end
+
+order = DynareOptions.order;
+
+% 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;
+ghs2 = ReducedForm.ghs2;
+if (order == 3)
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+end
+
+% Get covariance matrices
+Q = ReducedForm.Q;
+H = ReducedForm.H;
+if isempty(H)
+ H = 0;
+end
+
+% Get initial condition for the state vector.
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+state_variance_rank = size(StateVectorVarianceSquareRoot,2);
+%Q_lower_triangular_cholesky = chol(Q)';
+%if pruning
+% StateVectorMean_ = StateVectorMean;
+% StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
+%end
+
+% Set seed for randn().
+set_dynare_seed('default');
+
+% Initialization of the likelihood.
+const_lik = log(2*pi)*number_of_observed_variables;
+
+% Initialization of the weights across particles.
+weights = ones(1,number_of_particles)/number_of_particles ;
+StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
+%if pruning
+% StateVectors_ = StateVectors;
+%end
+yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
+%if pruning
+% yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
+% [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
+%else
+if (order == 2)
+ tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
+elseif (order == 3)
+ tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, ss, options_.threads.local_state_space_iteration_3, false);
+else
+ error('Orders > 3 not allowed');
+end
+%end
+PredictedObservedMean = weights*(tmp(mf1,:)');
+PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
+dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean');
+PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + H;
+wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
+tau_tilde = weights.*wtilde ;
+tau_tilde = tau_tilde/sum(tau_tilde);
+initial_distribution = resample(StateVectors',tau_tilde',ParticleOptions)' ;
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/auxiliary_particle_filter.m b/matlab/nonlinear-filters/auxiliary_particle_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..21caa73bd14362b156b565a2f5c7d071567299c1
--- /dev/null
+++ b/matlab/nonlinear-filters/auxiliary_particle_filter.m
@@ -0,0 +1,183 @@
+function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model)
+
+% Evaluates the likelihood of a nonlinear model with the auxiliary particle filter
+% allowing eventually resampling.
+%
+% Copyright © 2011-2022 Dynare Team
+%
+% This file is part of Dynare (particles module).
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare particles module is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+% Get perturbation order
+order = DynareOptions.order;
+
+% Set flag for prunning
+pruning = ParticleOptions.pruning;
+
+% Get steady state and mean.
+steadystate = ReducedForm.steadystate;
+constant = ReducedForm.constant;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+mf0 = ReducedForm.mf0;
+mf1 = ReducedForm.mf1;
+sample_size = size(Y,2);
+number_of_state_variables = length(mf0);
+number_of_observed_variables = length(mf1);
+number_of_structural_innovations = length(ReducedForm.Q);
+number_of_particles = ParticleOptions.number_of_particles;
+
+if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ udr = ReducedForm.udr;
+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;
+ ghs2 = ReducedForm.ghs2;
+ if (order == 3)
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+end
+
+% Get covariance matrices
+Q = ReducedForm.Q;
+H = ReducedForm.H;
+
+% Get initial condition for the state vector.
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
+state_variance_rank = size(StateVectorVarianceSquareRoot,2);
+Q_lower_triangular_cholesky = chol(Q)';
+
+% Set seed for randn().
+set_dynare_seed('default');
+
+% Initialization of the likelihood.
+const_lik = log(2*pi)*number_of_observed_variables+log(det(H));
+lik = NaN(sample_size,1);
+LIK = NaN;
+
+% Initialization of the weights across particles.
+weights = ones(1,number_of_particles)/number_of_particles ;
+StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
+%StateVectors = bsxfun(@plus,zeros(state_variance_rank,number_of_particles),StateVectorMean);
+if pruning
+ if order == 2
+ StateVectors_ = StateVectors;
+ state_variables_steady_state_ = state_variables_steady_state;
+ mf0_ = mf0;
+ elseif order == 3
+ StateVectors_ = repmat(StateVectors,3,1);
+ state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
+ mf0_ = repmat(mf0,1,3);
+ mask2 = number_of_state_variables+1:2*number_of_state_variables;
+ mask3 = 2*number_of_state_variables+1:3*number_of_state_variables;
+ mf0_(mask2) = mf0_(mask2)+size(ghx,1);
+ mf0_(mask3) = mf0_(mask3)+2*size(ghx,1);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+end
+
+for t=1:sample_size
+ yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
+ if pruning
+ yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_);
+ if order == 2
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ [tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+ else
+ if ReducedForm.use_k_order_solver
+ tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations,number_of_particles), dr, Model, DynareOptions, udr);
+ else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations,number_of_particles), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+ end
+ end
+ PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
+ z = sum(PredictionError.*(H\PredictionError),1) ;
+ tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ;
+ tau_tilde = tau_tilde/sum(tau_tilde) ;
+ indx = resample(0,tau_tilde',ParticleOptions);
+ if pruning
+ yhat_ = yhat_(:,indx) ;
+ end
+ yhat = yhat(:,indx) ;
+ weights_stage_1 = weights(indx)./tau_tilde(indx) ;
+ epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
+ if pruning
+ if order == 2
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ [tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+ StateVectors_ = tmp_(mf0_,:);
+ else
+ if ReducedForm.use_k_order_solver
+ tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
+ else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+ end
+ end
+ StateVectors = tmp(mf0,:);
+ PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
+ weights_stage_2 = weights_stage_1.*(exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) + 1e-99) ;
+ lik(t) = log(mean(weights_stage_2)) ;
+ weights = weights_stage_2/sum(weights_stage_2);
+ if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
+ if pruning
+ temp = resample([StateVectors' StateVectors_'],weights',ParticleOptions);
+ StateVectors = temp(:,1:number_of_state_variables)';
+ StateVectors_ = temp(:,number_of_state_variables+1:end)';
+ else
+ StateVectors = resample(StateVectors',weights',ParticleOptions)';
+ end
+ weights = ones(1,number_of_particles)/number_of_particles;
+ end
+end
+
+LIK = -sum(lik(start:end));
diff --git a/matlab/nonlinear-filters/conditional_filter_proposal.m b/matlab/nonlinear-filters/conditional_filter_proposal.m
new file mode 100644
index 0000000000000000000000000000000000000000..4d60a443586a21f21692065f9ecf44c9ab0d4cb9
--- /dev/null
+++ b/matlab/nonlinear-filters/conditional_filter_proposal.m
@@ -0,0 +1,164 @@
+function [ProposalStateVector, Weights, flag] = conditional_filter_proposal(ReducedForm, y, StateVectors, SampleWeights, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ...
+ H, ParticleOptions, ThreadsOptions, DynareOptions, Model)
+
+% Computes the proposal for each past particle using Gaussian approximations
+% for the state errors and the Kalman filter
+%
+% INPUTS
+% - ReducedForm [structure] Matlab's structure describing the reduced form model.
+% - y [double] p×1 vector, current observation (p is the number of observed variables).
+% - StateVectors
+% - SampleWeights
+% - Q_lower_triangular_cholesky
+% - H_lower_triangular_cholesky
+% - H
+% - ParticleOptions
+% - ThreadsOptions
+% - DynareOptions
+% - Model
+%
+% OUTPUTS
+% - ProposalStateVector
+% - Weights
+% - flag
+
+% Copyright © 2012-2022 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+flag = false;
+
+order = DynareOptions.order;
+
+if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ udr = ReducedForm.udr;
+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;
+ ghs2 = ReducedForm.ghs2;
+ if order == 3
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+end
+
+constant = ReducedForm.constant;
+steadystate = ReducedForm.steadystate;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+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);
+
+if ParticleOptions.proposal_approximation.montecarlo
+ nodes = randn(ParticleOptions.number_of_particles/10, number_of_structural_innovations);
+ weights = 1.0/ParticleOptions.number_of_particles;
+ weights_c = weights;
+elseif ParticleOptions.proposal_approximation.cubature
+ [nodes, weights] = spherical_radial_sigma_points(number_of_structural_innovations);
+ weights_c = weights;
+elseif ParticleOptions.proposal_approximation.unscented
+ [nodes, weights, weights_c] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions);
+else
+ error('Estimation: This approximation for the proposal is not implemented or unknown!')
+end
+
+epsilon = Q_lower_triangular_cholesky*nodes';
+yhat = repmat(StateVectors-state_variables_steady_state, 1, size(epsilon, 2));
+
+if ReducedForm.use_k_order_solver
+ tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
+else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+end
+
+PredictedStateMean = tmp(mf0,:)*weights;
+PredictedObservedMean = tmp(mf1,:)*weights;
+
+if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
+ PredictedStateMean = sum(PredictedStateMean, 2);
+ PredictedObservedMean = sum(PredictedObservedMean, 2);
+ dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean)'.*sqrt(weights);
+ dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)'.*sqrt(weights);
+ PredictedStateVariance = dState*dState';
+ big_mat = [dObserved dState; H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)];
+ [~, mat] = qr2(big_mat,0);
+ mat = mat';
+ 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));
+ Error = y-PredictedObservedMean;
+ StateVectorMean = PredictedStateMean+(CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*Error;
+ if ParticleOptions.cpf_weights_method.amisanotristani
+ Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1), PredictedObservedVarianceSquareRoot, Error);
+ end
+else
+ dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean);
+ dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean);
+ PredictedStateVariance = dState*diag(weights_c)*dState';
+ PredictedObservedVariance = dObserved*diag(weights_c)*dObserved'+H;
+ PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved';
+ KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance;
+ Error = y-PredictedObservedMean;
+ StateVectorMean = PredictedStateMean+KalmanFilterGain*Error;
+ StateVectorVariance = PredictedStateVariance-KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
+ StateVectorVariance = 0.5*(StateVectorVariance+StateVectorVariance');
+ [StateVectorVarianceSquareRoot, p] = chol(StateVectorVariance, 'lower') ;
+ if p
+ flag = true;
+ ProposalStateVector = zeros(number_of_state_variables, 1);
+ Weights = 0.0;
+ return
+ end
+ if ParticleOptions.cpf_weights_method.amisanotristani
+ Weights = SampleWeights.*probability2(zeros(number_of_observed_variables, 1), chol(PredictedObservedVariance)', Error);
+ end
+end
+
+ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot, 2), 1)+StateVectorMean;
+if ParticleOptions.cpf_weights_method.murrayjonesparslow
+ PredictedStateVariance = 0.5*(PredictedStateVariance+PredictedStateVariance');
+ [PredictedStateVarianceSquareRoot, p] = chol(PredictedStateVariance, 'lower');
+ if p
+ flag = true;
+ ProposalStateVector = zeros(number_of_state_variables,1);
+ Weights = 0.0;
+ return
+ end
+ Prior = probability2(PredictedStateMean, PredictedStateVarianceSquareRoot, ProposalStateVector);
+ Posterior = probability2(StateVectorMean, StateVectorVarianceSquareRoot, ProposalStateVector);
+ Likelihood = probability2(y, H_lower_triangular_cholesky, measurement_equations(ProposalStateVector, ReducedForm, ThreadsOptions, DynareOptions, Model));
+ Weights = SampleWeights.*Likelihood.*(Prior./Posterior);
+end
diff --git a/matlab/nonlinear-filters/conditional_particle_filter.m b/matlab/nonlinear-filters/conditional_particle_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..db204dc855b7e7d757a53da9bcdf3f5261bb4216
--- /dev/null
+++ b/matlab/nonlinear-filters/conditional_particle_filter.m
@@ -0,0 +1,110 @@
+function [LIK,lik] = conditional_particle_filter(ReducedForm, Y, s, ParticleOptions, ThreadsOptions, DynareOptions, Model)
+
+% Evaluates the likelihood of a non-linear model with a particle filter
+%
+% INPUTS
+% - ReducedForm [structure] Matlab's structure describing the reduced form model.
+% - Y [double] p×T matrix of (detrended) data, where p is the number of observed variables.
+% - s [integer] scalar, likelihood evaluation starts at s (has to be smaller than T, the sample length provided in Y).
+% - ParticlesOptions [struct]
+% - ThreadsOptions [struct]
+% - DynareOptions [struct]
+% - Model [struct]
+%
+% OUTPUTS
+% - LIK [double] scalar, likelihood
+% - lik [double] (T-s+1)×1 vector, density of observations in each period.
+%
+% REMARKS
+% - The proposal is built using the Kalman updating step for each particle.
+% - we need draws in the errors distributions
+% Whether we use Monte-Carlo draws from a multivariate gaussian distribution
+% as in Amisano & Tristani (JEDC 2010).
+% Whether we use multidimensional Gaussian sparse grids approximations:
+% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak
+% operator (ref: Winschel & Kratzig, 2010).
+% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009a,2009b).
+% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der
+% Merwe & Wan 2003).
+%
+% Pros:
+% - Allows using current observable information in the proposal
+% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach
+% Cons:
+% - The use of the Kalman updating step may biais the proposal distribution since
+% it has been derived in a linear context and is implemented in a nonlinear
+% context. That is why particle resampling is performed.
+
+% Copyright © 2009-2020 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% Set default for third input argument.
+if isempty(s)
+ s = 1;
+end
+
+T = size(Y,2);
+p = length(ReducedForm.mf1);
+n = ParticleOptions.number_of_particles;
+
+% Get covariance matrices
+Q = ReducedForm.Q;
+H = ReducedForm.H;
+if isempty(H)
+ H = 0;
+ H_lower_triangular_cholesky = 0;
+else
+ H_lower_triangular_cholesky = chol(H)';
+end
+
+% Get initial condition for the state vector.
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
+state_variance_rank = size(StateVectorVarianceSquareRoot, 2);
+Q_lower_triangular_cholesky = chol(Q)';
+
+% Set seed for randn().
+set_dynare_seed('default');
+
+% Initialization of the likelihood.
+lik = NaN(T, 1);
+ks = 0;
+StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, n), StateVectorMean);
+SampleWeights = ones(1, n)/n;
+
+for t=1:T
+ flags = false(n, 1);
+ for i=1:n
+ [StateParticles(:,i), SampleWeights(i), flags(i)] = ...
+ conditional_filter_proposal(ReducedForm, Y(:,t), StateParticles(:,i), SampleWeights(i), Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ParticleOptions, ThreadsOptions, DynareOptions, Model);
+ end
+ if any(flags)
+ LIK = -Inf;
+ lik(t) = -Inf;
+ return
+ end
+ SumSampleWeights = sum(SampleWeights);
+ lik(t) = log(SumSampleWeights);
+ SampleWeights = SampleWeights./SumSampleWeights;
+ if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*T) || ParticleOptions.resampling.status.systematic
+ ks = ks + 1;
+ StateParticles = resample(StateParticles', SampleWeights', ParticleOptions)';
+ SampleWeights = ones(1, n)/n;
+ end
+end
+
+LIK = -sum(lik(s:end));
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/fit_gaussian_mixture.m b/matlab/nonlinear-filters/fit_gaussian_mixture.m
new file mode 100644
index 0000000000000000000000000000000000000000..9836fa59e6ffe33cf49fdf77388e5891d6f14c63
--- /dev/null
+++ b/matlab/nonlinear-filters/fit_gaussian_mixture.m
@@ -0,0 +1,51 @@
+function [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(X,X_weights,StateMu,StateSqrtP,StateWeights,crit,niters,check)
+
+% Copyright © 2013-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+[dim,Ndata] = size(X);
+M = size(StateMu,2) ;
+if check % Ensure that covariances don't collapse
+ MIN_COVAR_SQRT = sqrt(eps);
+ init_covars = StateSqrtP;
+end
+eold = -Inf;
+for n=1:niters
+ % Calculate posteriors based on old parameters
+ [prior,likelihood,marginal,posterior] = probability3(StateMu,StateSqrtP,StateWeights,X,X_weights);
+ e = sum(log(marginal));
+ if (n > 1 && abs((e - eold)/eold) < crit)
+ return;
+ else
+ eold = e;
+ end
+ new_pr = (sum(posterior,2))';
+ StateWeights = new_pr/Ndata;
+ StateMu = bsxfun(@rdivide,(posterior*X')',new_pr);
+ for j=1:M
+ diffs = bsxfun(@minus,X,StateMu(:,j));
+ tpost = (1/sqrt(new_pr(j)))*sqrt(posterior(j,:));
+ diffs = bsxfun(@times,diffs,tpost);
+ [foo,tcov] = qr2(diffs',0);
+ StateSqrtP(:,:,j) = tcov';
+ if check
+ if min(abs(diag(StateSqrtP(:,:,j)))) < MIN_COVAR_SQRT
+ StateSqrtP(:,:,j) = init_covars(:,:,j);
+ end
+ end
+ end
+end
diff --git a/matlab/nonlinear-filters/gaussian_densities.m b/matlab/nonlinear-filters/gaussian_densities.m
new file mode 100644
index 0000000000000000000000000000000000000000..727d083c49c6f5285cb31eab41cc9c9cd57b1f3c
--- /dev/null
+++ b/matlab/nonlinear-filters/gaussian_densities.m
@@ -0,0 +1,50 @@
+function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sqr_Pss_t_t_1,particles,H,normconst,weigths1,weigths2,ReducedForm,ThreadsOptions,DynareOptions, Model)
+%
+% Elements to calculate the importance sampling ratio
+%
+% INPUTS
+% reduced_form_model [structure] Matlab's structure describing the reduced form model.
+% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
+% reduced_form_model.state.dr [structure] output of resol.m.
+% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
+% start [integer] scalar, likelihood evaluation starts at 'start'.
+% smolyak_accuracy [integer] scalar.
+%
+% OUTPUTS
+% LIK [double] scalar, likelihood
+% lik [double] vector, density of observations in each period.
+%
+% REFERENCES
+%
+% NOTES
+% The vector "lik" is used to evaluate the jacobian of the likelihood.
+
+% Copyright © 2009-2019 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% proposal density
+proposal = probability2(mut_t, sqr_Pss_t_t, particles);
+
+% prior density
+prior = probability2(st_t_1, sqr_Pss_t_t_1, particles);
+
+% likelihood
+yt_t_1_i = measurement_equations(particles, ReducedForm, ThreadsOptions, DynareOptions, Model);
+likelihood = probability2(obs, sqrt(H), yt_t_1_i);
+
+IncrementalWeights = likelihood.*prior./proposal;
diff --git a/matlab/nonlinear-filters/gaussian_filter.m b/matlab/nonlinear-filters/gaussian_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..9b1b29e4c8cefc4e4bcd8217d7a2a558dc16fc26
--- /dev/null
+++ b/matlab/nonlinear-filters/gaussian_filter.m
@@ -0,0 +1,135 @@
+function [LIK,lik] = gaussian_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model)
+
+% Evaluates the likelihood of a non-linear model approximating the
+% predictive (prior) and filtered (posterior) densities for state variables
+% by gaussian distributions.
+% Gaussian approximation is done by:
+% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009).
+% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1995)
+% - Monte-Carlo draws from a multivariate gaussian distribution.
+% First and second moments of prior and posterior state densities are computed
+% from the resulting nodes/particles and allows to generate new distributions at the
+% following observation.
+% Pros: The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles
+% filters since it treats a lesser number of particles. Furthermore, in all cases, there is no need
+% of resampling.
+% Cons: estimations may be biaised if the model is truly non-gaussian
+% since predictive and filtered densities are unimodal.
+%
+% INPUTS
+% Reduced_Form [structure] Matlab's structure describing the reduced form model.
+% Y [double] matrix of original observed variables.
+% start [double] structural parameters.
+% ParticleOptions [structure] Matlab's structure describing options concerning particle filtering.
+% ThreadsOptions [structure] Matlab's structure.
+%
+% OUTPUTS
+% LIK [double] scalar, likelihood
+% lik [double] vector, density of observations in each period.
+%
+% REFERENCES
+%
+% NOTES
+% The vector "lik" is used to evaluate the jacobian of the likelihood.
+
+% Copyright © 2009-2019 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+mf0 = ReducedForm.mf0;
+mf1 = ReducedForm.mf1;
+sample_size = size(Y,2);
+number_of_state_variables = length(mf0);
+number_of_observed_variables = length(mf1);
+number_of_particles = ParticleOptions.number_of_particles;
+
+% compute gaussian quadrature nodes and weights on states and shocks
+if ParticleOptions.distribution_approximation.cubature
+ [nodes2, weights2] = spherical_radial_sigma_points(number_of_state_variables);
+ weights_c2 = weights2;
+elseif ParticleOptions.distribution_approximation.unscented
+ [nodes2, weights2, weights_c2] = unscented_sigma_points(number_of_state_variables,ParticleOptions);
+else
+ if ~ParticleOptions.distribution_approximation.montecarlo
+ error('This approximation for the proposal is unknown!')
+ end
+end
+
+if ParticleOptions.distribution_approximation.montecarlo
+ set_dynare_seed('default');
+end
+
+% Get covariance matrices
+Q = ReducedForm.Q;
+H = ReducedForm.H;
+if isempty(H)
+ H = 0;
+ H_lower_triangular_cholesky = 0;
+else
+ H_lower_triangular_cholesky = reduced_rank_cholesky(H)';
+end
+
+% Get initial condition for the state vector.
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+state_variance_rank = size(StateVectorVarianceSquareRoot,2);
+Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)';
+
+% Initialization of the likelihood.
+const_lik = (2*pi)^(number_of_observed_variables/2) ;
+lik = NaN(sample_size,1);
+LIK = NaN;
+
+for t=1:sample_size
+ [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
+ gaussian_filter_bank(ReducedForm, Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, ...
+ H, ParticleOptions, ThreadsOptions, DynareOptions, Model);
+ if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented
+ StateParticles = bsxfun(@plus, StateVectorMean, StateVectorVarianceSquareRoot*nodes2');
+ IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
+ PredictedStateVarianceSquareRoot, StateParticles, H, const_lik, ...
+ weights2, weights_c2, ReducedForm, ThreadsOptions, ...
+ DynareOptions, Model);
+ SampleWeights = weights2.*IncrementalWeights;
+ else
+ StateParticles = bsxfun(@plus, StateVectorVarianceSquareRoot*randn(state_variance_rank, number_of_particles), StateVectorMean) ;
+ IncrementalWeights = gaussian_densities(Y(:,t), StateVectorMean, StateVectorVarianceSquareRoot, PredictedStateMean, ...
+ PredictedStateVarianceSquareRoot,StateParticles,H,const_lik, ...
+ 1/number_of_particles,1/number_of_particles,ReducedForm,ThreadsOptions, ...
+ DynareOptions, Model);
+ SampleWeights = IncrementalWeights/number_of_particles;
+ end
+ SampleWeights = SampleWeights + 1e-6*ones(size(SampleWeights, 1), 1);
+ SumSampleWeights = sum(SampleWeights);
+ lik(t) = log(SumSampleWeights);
+ SampleWeights = SampleWeights./SumSampleWeights;
+ if not(ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented)
+ if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
+ StateParticles = resample(StateParticles', SampleWeights, ParticleOptions)';
+ SampleWeights = ones(number_of_particles, 1)/number_of_particles;
+ end
+ end
+ StateVectorMean = StateParticles*SampleWeights;
+ temp = bsxfun(@minus, StateParticles, StateVectorMean);
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky(bsxfun(@times,SampleWeights',temp)*temp')';
+end
+
+LIK = -sum(lik(start:end));
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/gaussian_filter_bank.m b/matlab/nonlinear-filters/gaussian_filter_bank.m
new file mode 100644
index 0000000000000000000000000000000000000000..17ae110cc63bf063241ca63790439e2581bcef8c
--- /dev/null
+++ b/matlab/nonlinear-filters/gaussian_filter_bank.m
@@ -0,0 +1,139 @@
+function [PredictedStateMean, PredictedStateVarianceSquareRoot, StateVectorMean, StateVectorVarianceSquareRoot] = ...
+ gaussian_filter_bank(ReducedForm, obs, StateVectorMean, StateVectorVarianceSquareRoot, Q_lower_triangular_cholesky, H_lower_triangular_cholesky, H, ...
+ ParticleOptions, ThreadsOptions, DynareOptions, Model)
+%
+% Computes the proposal with a gaussian approximation for importance
+% sampling
+% This proposal is a gaussian distribution calculated à la Kalman
+%
+% INPUTS
+% reduced_form_model [structure] Matlab's structure describing the reduced form model.
+% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
+% reduced_form_model.state.dr [structure] output of resol.m.
+% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
+%
+% OUTPUTS
+% LIK [double] scalar, likelihood
+% lik [double] vector, density of observations in each period.
+%
+% REFERENCES
+%
+% NOTES
+% The vector "lik" is used to evaluate the jacobian of the likelihood.
+
+% Copyright © 2009-2022 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+order = DynareOptions.order;
+
+if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ udr = ReducedForm.udr;
+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;
+ ghs2 = ReducedForm.ghs2;
+ if order == 3
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+end
+
+constant = ReducedForm.constant;
+steadystate = ReducedForm.steadystate;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+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);
+
+if ParticleOptions.proposal_approximation.montecarlo
+ nodes = randn(ParticleOptions.number_of_particles, number_of_state_variables+number_of_structural_innovations) ;
+ weights = 1/ParticleOptions.number_of_particles ;
+ weights_c = weights ;
+elseif ParticleOptions.proposal_approximation.cubature
+ [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations) ;
+ weights_c = weights ;
+elseif ParticleOptions.proposal_approximation.unscented
+ [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations, ParticleOptions);
+else
+ error('This approximation for the proposal is not implemented or unknown!')
+end
+
+xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1)] ;
+sqr_Px = [ StateVectorVarianceSquareRoot, zeros(number_of_state_variables, number_of_structural_innovations);
+ zeros(number_of_structural_innovations, number_of_state_variables) Q_lower_triangular_cholesky];
+sigma_points = bsxfun(@plus, xbar, sqr_Px*(nodes'));
+StateVectors = sigma_points(1:number_of_state_variables,:);
+epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
+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, udr);
+else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+end
+
+PredictedStateMean = tmp(mf0,:)*weights;
+PredictedObservedMean = tmp(mf1,:)*weights;
+
+if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
+ PredictedStateMean = sum(PredictedStateMean, 2);
+ PredictedObservedMean = sum(PredictedObservedMean, 2);
+ dState = bsxfun(@minus,tmp(mf0,:), PredictedStateMean)'.*sqrt(weights);
+ dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean)'.*sqrt(weights);
+ PredictedStateVarianceSquareRoot = chol(dState'*dState)';
+ big_mat = [dObserved, dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables,number_of_state_variables)];
+ [~, mat] = qr2(big_mat, 0);
+ mat = mat';
+ 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));
+ PredictionError = obs - PredictedObservedMean;
+ StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError;
+else
+ dState = bsxfun(@minus, tmp(mf0,:), PredictedStateMean);
+ dObserved = bsxfun(@minus, tmp(mf1,:), PredictedObservedMean);
+ PredictedStateVariance = dState*diag(weights_c)*dState';
+ PredictedObservedVariance = dObserved*diag(weights_c)*dObserved' + H;
+ PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved';
+ PredictedStateVarianceSquareRoot = chol(PredictedStateVariance)';
+ PredictionError = obs - PredictedObservedMean;
+ KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance;
+ StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError;
+ StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
+ StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
+end
diff --git a/matlab/nonlinear-filters/gaussian_mixture_densities.m b/matlab/nonlinear-filters/gaussian_mixture_densities.m
new file mode 100644
index 0000000000000000000000000000000000000000..0cf9dd79e5aa0372af0ed9a8ddf475c8e08db933
--- /dev/null
+++ b/matlab/nonlinear-filters/gaussian_mixture_densities.m
@@ -0,0 +1,55 @@
+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
+% reduced_form_model [structure] Matlab's structure describing the reduced form model.
+% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
+% reduced_form_model.state.dr [structure] output of resol.m.
+% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
+% start [integer] scalar, likelihood evaluation starts at 'start'.
+% smolyak_accuracy [integer] scalar.
+%
+% OUTPUTS
+% LIK [double] scalar, likelihood
+% lik [double] vector, density of observations in each period.
+%
+% REFERENCES
+%
+% NOTES
+% The vector "lik" is used to evaluate the jacobian of the likelihood.
+
+% Copyright © 2009-2019 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% Compute the density of particles under the prior distribution
+[~, ~, prior] = probability(StateMuPrior, StateSqrtPPrior, StateWeightsPrior, StateParticles);
+prior = prior';
+
+% Compute the density of particles under the proposal distribution
+[~, ~, 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, DynareOptions, Model);
+
+% likelihood
+likelihood = probability2(obs, sqrt(H), yt_t_1_i);
+IncrementalWeights = likelihood.*prior./proposal;
diff --git a/matlab/nonlinear-filters/gaussian_mixture_filter.m b/matlab/nonlinear-filters/gaussian_mixture_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..e06461b85acbbf61b0ac98a4cd2815aaab4ea5a0
--- /dev/null
+++ b/matlab/nonlinear-filters/gaussian_mixture_filter.m
@@ -0,0 +1,226 @@
+function [LIK, lik] = gaussian_mixture_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model)
+
+% Evaluates the likelihood of a non-linear model approximating the state
+% variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing
+% a wide variety of generalized distributions (when multimodal for instance).
+% Each gaussian distribution is obtained whether
+% - with a radial-spherical cubature
+% - with scaled unscented sigma-points
+% A Sparse grid Kalman Filter is implemented on each component of the mixture,
+% which confers it a weight about current information.
+% Information on the current observables is then embodied in the proposal
+% distribution in which we draw particles, which allows
+% - reaching a greater precision relatively to a standard particle filter,
+% - reducing the number of particles needed,
+% - still being faster.
+%
+%
+% INPUTS
+% reduced_form_model [structure] Matlab's structure describing the reduced form model.
+% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
+% reduced_form_model.state.dr [structure] output of resol.m.
+% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
+% start [integer] scalar, likelihood evaluation starts at 'start'.
+%
+% OUTPUTS
+% LIK [double] scalar, likelihood
+% lik [double] vector, density of observations in each period.
+%
+% REFERENCES
+%
+% Van der Meerwe & Wan, Gaussian Mixture Sigma-Point Particle Filters for Sequential
+% Probabilistic Inference in Dynamic State-Space Models.
+% Heiss & Winschel, 2010, Journal of Applied Economics.
+% Winschel & Kratzig, 2010, Econometrica.
+%
+% NOTES
+% The vector "lik" is used to evaluate the jacobian of the likelihood.
+% Copyright © 2009-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+mf0 = ReducedForm.mf0;
+mf1 = ReducedForm.mf1;
+sample_size = size(Y,2);
+number_of_state_variables = length(mf0);
+number_of_observed_variables = length(mf1);
+number_of_structural_innovations = length(ReducedForm.Q);
+G = ParticleOptions.mixture_state_variables; % number of GM components in state
+number_of_particles = ParticleOptions.number_of_particles;
+
+% compute gaussian quadrature nodes and weights on states and shocks
+if ParticleOptions.distribution_approximation.cubature
+ [nodes, weights] = spherical_radial_sigma_points(number_of_state_variables);
+elseif ParticleOptions.distribution_approximation.unscented
+ [nodes, weights] = unscented_sigma_points(number_of_state_variables, ParticleOptions);
+else
+ if ~ParticleOptions.distribution_approximation.montecarlo
+ error('This approximation for the proposal is unknown!')
+ end
+end
+
+if ParticleOptions.distribution_approximation.montecarlo
+ set_dynare_seed('default');
+end
+
+% Get covariance matrices
+Q = ReducedForm.Q;
+H = ReducedForm.H;
+if isempty(H)
+ H = 0;
+ H_lower_triangular_cholesky = 0;
+else
+ H_lower_triangular_cholesky = reduced_rank_cholesky(H)';
+end
+Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)';
+
+% Initialize mixtures
+StateWeights = ones(1, G)/G;
+StateMu = ReducedForm.StateVectorMean;
+StateSqrtP = zeros(number_of_state_variables, number_of_state_variables, G);
+temp = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+StateMu = bsxfun(@plus, StateMu, bsxfun(@times,diag(temp), (-(G-1)/2:1:(G-1)/2))/10);
+for g=1:G
+ StateSqrtP(:,:,g) = temp/sqrt(G) ;
+end
+
+if ~ParticleOptions.mixture_structural_shocks
+ StructuralShocksMu = zeros(1, number_of_structural_innovations);
+ StructuralShocksWeights = 1;
+ I = 1;
+ StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu';
+ StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I);
+ StructuralShocksSqrtP(:,:,1) = Q_lower_triangular_cholesky;
+elseif ParticleOptions.mixture_structural_shocks==1
+ if ParticleOptions.proposal_approximation.cubature
+ [StructuralShocksMu, StructuralShocksWeights] = spherical_radial_sigma_points(number_of_structural_innovations);
+ StructuralShocksWeights = ones(size(StructuralShocksMu, 1), 1)*StructuralShocksWeights;
+ elseif ParticleOptions.proposal_approximation.unscented
+ [StructuralShocksMu, StructuralShocksWeights] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions);
+ else
+ if ~ParticleOptions.distribution_approximation.montecarlo
+ error('This approximation for the proposal is unknown!')
+ end
+ end
+ I = size(StructuralShocksWeights, 1);
+ StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu';
+ StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I);
+ for i=1:I
+ StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky;
+ end
+else
+ if ParticleOptions.proposal_approximation.cubature
+ [StructuralShocksMu, StructuralShocksWeights] = spherical_radial_sigma_points(number_of_structural_innovations);
+ StructuralShocksWeights = ones(size(StructuralShocksMu, 1), 1)*StructuralShocksWeights ;
+ elseif ParticleOptions.proposal_approximation.unscented
+ [StructuralShocksMu, StructuralShocksWeights] = unscented_sigma_points(number_of_structural_innovations, ParticleOptions);
+ else
+ if ~ParticleOptions.distribution_approximation.montecarlo
+ error('This approximation for the proposal is unknown!')
+ end
+ end
+ I = size(StructuralShocksWeights, 1);
+ StructuralShocksMu = Q_lower_triangular_cholesky*StructuralShocksMu';
+ StructuralShocksSqrtP = zeros(number_of_structural_innovations, number_of_structural_innovations, I);
+ for i=1:I
+ StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky/sqrt(StructuralShocksWeights(i));
+ end
+end
+
+ObservationShocksWeights = 1;
+J = 1 ;
+
+Gprime = G*I;
+Gsecond = G*I*J;
+SampleWeights = ones(Gsecond, 1)/Gsecond;
+
+StateWeightsPrior = zeros(1,Gprime);
+StateMuPrior = zeros(number_of_state_variables,Gprime);
+StateSqrtPPrior = zeros(number_of_state_variables, number_of_state_variables, Gprime);
+
+StateWeightsPost = zeros(1, Gsecond);
+StateMuPost = zeros(number_of_state_variables, Gsecond);
+StateSqrtPPost = zeros(number_of_state_variables, number_of_state_variables, Gsecond);
+
+const_lik = (2*pi)^(.5*number_of_observed_variables);
+
+lik = NaN(sample_size, 1);
+LIK = NaN;
+for t=1:sample_size
+ % Build the proposal joint quadratures of Gaussian on states, structural
+ % shocks and observation shocks based on each combination of mixtures
+ for i=1:I
+ for j=1:J
+ for g=1:G
+ gprime = g + (i-1)*G;
+ gsecond = gprime + (j-1)*Gprime;
+ [StateMuPrior(:,gprime), StateSqrtPPrior(:,:,gprime), StateWeightsPrior(1,gprime), ...
+ StateMuPost(:,gsecond), StateSqrtPPost(:,:,gsecond), StateWeightsPost(1,gsecond)] = ...
+ gaussian_mixture_filter_bank(ReducedForm,Y(:,t), StateMu(:,g), StateSqrtP(:,:,g), StateWeights(g),...
+ StructuralShocksMu(:,i), StructuralShocksSqrtP(:,:,i), StructuralShocksWeights(i),...
+ ObservationShocksWeights(j), H, H_lower_triangular_cholesky, const_lik, ...
+ ParticleOptions, ThreadsOptions, DynareOptions, Model);
+ end
+ end
+ end
+
+ % Normalize weights
+ StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior, 2);
+ StateWeightsPost = StateWeightsPost/sum(StateWeightsPost, 2);
+
+ if ParticleOptions.distribution_approximation.cubature || ParticleOptions.distribution_approximation.unscented
+ for i=1:Gsecond
+ StateParticles = bsxfun(@plus, StateMuPost(:,i), StateSqrtPPost(:,:,i)*nodes');
+ IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
+ StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
+ ReducedForm, ThreadsOptions, DynareOptions, Model);
+ SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights);
+ end
+ SumSampleWeights = sum(SampleWeights);
+ lik(t) = log(SumSampleWeights);
+ SampleWeights = SampleWeights./SumSampleWeights;
+ [~, SortedRandomIndx] = sort(rand(1,Gsecond));
+ SortedRandomIndx = SortedRandomIndx(1:G);
+ indx = resample(0,SampleWeights,ParticleOptions);
+ indx = indx(SortedRandomIndx);
+ StateMu = StateMuPost(:,indx);
+ StateSqrtP = StateSqrtPPost(:,:,indx);
+ StateWeights = ones(1,G)/G;
+ else
+ % Sample particle in the proposal distribution, ie the posterior state GM
+ StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles);
+ IncrementalWeights = gaussian_mixture_densities(Y(:,t), StateMuPrior, StateSqrtPPrior, StateWeightsPrior, ...
+ StateMuPost, StateSqrtPPost, StateWeightsPost, StateParticles, H, ...
+ ReducedForm, ThreadsOptions, DynareOptions, Model);
+ SampleWeights = IncrementalWeights/number_of_particles;
+ SumSampleWeights = sum(SampleWeights,1);
+ SampleWeights = SampleWeights./SumSampleWeights;
+ lik(t) = log(SumSampleWeights);
+ if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
+ StateParticles = resample(StateParticles',SampleWeights',ParticleOptions)';
+ SampleWeights = ones(number_of_particles,1)/number_of_particles;
+ end
+ [StateMu, StateSqrtP, StateWeights] = fit_gaussian_mixture(StateParticles, SampleWeights', StateMu, StateSqrtP, StateWeights, 0.001, 10, 1);
+ end
+end
+
+LIK = -sum(lik(start:end));
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/gaussian_mixture_filter_bank.m b/matlab/nonlinear-filters/gaussian_mixture_filter_bank.m
new file mode 100644
index 0000000000000000000000000000000000000000..acbe77f2b07dfdbd1082d1644efe939555a1b277
--- /dev/null
+++ b/matlab/nonlinear-filters/gaussian_mixture_filter_bank.m
@@ -0,0 +1,146 @@
+function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPPost,StateWeightsPost] =...
+ 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
+%
+% INPUTS
+% reduced_form_model [structure] Matlab's structure describing the reduced form model.
+% reduced_form_model.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% reduced_form_model.state.Q [double] (qq x qq) variance matrix of state errors.
+% reduced_form_model.state.dr [structure] output of resol.m.
+% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
+%
+% OUTPUTS
+% LIK [double] scalar, likelihood
+% lik [double] vector, density of observations in each period.
+%
+% REFERENCES
+%
+% NOTES
+% The vector "lik" is used to evaluate the jacobian of the likelihood.
+
+% Copyright © 2009-2022 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+order = DynareOptions.order;
+
+if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ udr = ReducedForm.udr;
+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;
+ ghs2 = ReducedForm.ghs2;
+ if order == 3
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+end
+
+constant = ReducedForm.constant;
+steadystate = ReducedForm.steadystate;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+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;
+
+if ParticleOptions.proposal_approximation.cubature
+ [nodes3, weights3] = spherical_radial_sigma_points(numb);
+ weights_c3 = weights3;
+elseif ParticleOptions.proposal_approximation.unscented
+ [nodes3, weights3, weights_c3] = unscented_sigma_points(numb, ParticleOptions);
+else
+ 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);
+if ReducedForm.use_k_order_solver
+ tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
+else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+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);
+ PredictedStateVariance = dState'*dState;
+ big_mat = [dObserved, dState ; H_lower_triangular_cholesky, zeros(number_of_observed_variables, number_of_state_variables)];
+ [~, mat] = qr2(big_mat, 0);
+ mat = mat';
+ 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);
+ PredictedStateVariance = dState*diag(weights_c3)*dState';
+ PredictedObservedVariance = dObserved*diag(weights_c3)*dObserved' + H;
+ PredictedStateAndObservedCovariance = dState*diag(weights_c3)*dObserved';
+ sqrdet = sqrt(det(PredictedObservedVariance));
+ iPredictedObservedVariance = inv(PredictedObservedVariance);
+ PredictionError = obs - PredictedObservedMean;
+ KalmanFilterGain = PredictedStateAndObservedCovariance*iPredictedObservedVariance;
+ StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError;
+ StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
+ StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
+end
+
+data_lik_GM_g = exp(-0.5*PredictionError'*iPredictedObservedVariance*PredictionError)/abs(normfactO*sqrdet) + 1e-99;
+StateMuPrior = PredictedStateMean;
+StateSqrtPPrior = reduced_rank_cholesky(PredictedStateVariance)';
+StateWeightsPrior = StateWeights*StructuralShocksWeights;
+StateMuPost = StateVectorMean;
+StateSqrtPPost = StateVectorVarianceSquareRoot;
+StateWeightsPost = StateWeightsPrior*ObservationShocksWeights*data_lik_GM_g;
diff --git a/matlab/nonlinear-filters/importance_sampling.m b/matlab/nonlinear-filters/importance_sampling.m
new file mode 100644
index 0000000000000000000000000000000000000000..d477f9f50d6f477e9b7103696535f1a328c58e8c
--- /dev/null
+++ b/matlab/nonlinear-filters/importance_sampling.m
@@ -0,0 +1,28 @@
+function State_Particles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost,numP)
+
+% Copyright © 2013-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+[Xdim,Gsecond] = size(StateMuPost) ;
+u = rand(numP,1);
+[Nc,comp] = histc(u, cumsum([0; StateWeightsPost]));
+State_Particles = zeros(Xdim,numP);
+for k=1:Gsecond
+ idx = bsxfun(@eq,comp,k*ones(size(comp))) ;
+ State_Particles(:,idx) = StateSqrtPPost(:,:,k)*randn(Xdim,Nc(k));
+end
+State_Particles= State_Particles + StateMuPost(:,comp);
diff --git a/matlab/nonlinear-filters/measurement_equations.m b/matlab/nonlinear-filters/measurement_equations.m
new file mode 100644
index 0000000000000000000000000000000000000000..8ed06349d97652ab032d6330bb0d7f2f1c4bc9a8
--- /dev/null
+++ b/matlab/nonlinear-filters/measurement_equations.m
@@ -0,0 +1,57 @@
+function measure = measurement_equations(StateVectors,ReducedForm,ThreadsOptions, DynareOptions, Model)
+
+% Copyright © 2013-2022 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+order = DynareOptions.order;
+mf1 = ReducedForm.mf1;
+if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ udr = ReducedForm.udr;
+else
+ ghx = ReducedForm.ghx(mf1,:);
+ ghu = ReducedForm.ghu(mf1,:);
+ ghxx = ReducedForm.ghxx(mf1,:);
+ ghuu = ReducedForm.ghuu(mf1,:);
+ ghxu = ReducedForm.ghxu(mf1,:);
+ ghs2 = ReducedForm.ghs2(mf1,:);
+ if order == 3
+ ghxxx = ReducedForm.ghxxx(mf1,:);
+ ghuuu = ReducedForm.ghuuu(mf1,:);
+ ghxxu = ReducedForm.ghxxu(mf1,:);
+ ghxuu = ReducedForm.ghxuu(mf1,:);
+ ghxss = ReducedForm.ghxss(mf1,:);
+ ghuss = ReducedForm.ghuss(mf1,:);
+ end
+end
+steadystate = ReducedForm.steadystate(mf1,:);
+constant = ReducedForm.constant(mf1,:);
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+number_of_structural_innovations = length(ReducedForm.Q);
+yhat = bsxfun(@minus, StateVectors, state_variables_steady_state);
+if ReducedForm.use_k_order_solver
+ tmp = local_state_space_iteration_k(yhat, zeros(number_of_structural_innovations, size(yhat,2)), dr, Model, DynareOptions, udr);
+ measure = tmp(mf1,:);
+else
+ if order == 2
+ measure = local_state_space_iteration_2(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ measure = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, size(yhat,2)), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+end
diff --git a/matlab/nonlinear-filters/multivariate_smooth_resampling.m b/matlab/nonlinear-filters/multivariate_smooth_resampling.m
new file mode 100644
index 0000000000000000000000000000000000000000..a83d348c2e08a5cfe5bfe90a93ad2c8256711edb
--- /dev/null
+++ b/matlab/nonlinear-filters/multivariate_smooth_resampling.m
@@ -0,0 +1,72 @@
+function new_particles = multivariate_smooth_resampling(particles,weights)
+% Smooth Resampling of the particles.
+
+%@info:
+%! @deftypefn {Function File} {@var{new_particles} =} multivariate_smooth_resampling (@var{weights}, @var{particles}, @var{number_of_new_particles}, @var{number_of_partitions})
+%! @anchor{particle/multivariate_smooth_resampling}
+%! @sp 1
+%! Smooth Resampling of the particles (multivariate version).
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item weights
+%! n*1 vector of doubles, particles' weights.
+%! @item particles
+%! n*1 vector of doubles, particles.
+%! @item number_of_new_particles
+%! Integer scalar.
+%! @item number_of_partitions
+%! Integer scalar.
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item indx
+%! number_of_new_particles*1 vector of doubles, new particles.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 1
+%! @ref{particle/sequantial_importance_particle_filter}
+%! @sp 2
+%! @strong{This function calls:}
+%! @sp 1
+%! @ref{particle/univariate_smooth_resampling}
+%! @sp 2
+%! @end deftypefn
+%@eod:
+
+% Copyright © 2012-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr
+% stephane DOT adjemian AT univ DASH lemans DOT fr
+
+number_of_particles = length(weights);
+number_of_states = size(particles,2);
+[P,D] = eig(particles'*(bsxfun(@times,1/number_of_particles,particles))) ;
+D = diag(D) ;
+vectors = bsxfun(@times,P,sqrt(D)') ;
+orthogonalized_particles = bsxfun(@rdivide,particles*vectors,D') ;
+new_particles = zeros(number_of_particles,number_of_states) ;
+for j=1:number_of_states
+ tout = sortrows( [orthogonalized_particles(:,j) weights],1) ;
+ new_particles(:,j) = univariate_smooth_resampling(tout(:,2),tout(:,1),number_of_particles) ;
+end
+new_particles = new_particles*(vectors') ;
diff --git a/matlab/nonlinear-filters/mykmeans.m b/matlab/nonlinear-filters/mykmeans.m
new file mode 100644
index 0000000000000000000000000000000000000000..e7c424b00646f8ccf3dbcee70395b218c632cd54
--- /dev/null
+++ b/matlab/nonlinear-filters/mykmeans.m
@@ -0,0 +1,53 @@
+function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod)
+
+% Copyright © 2013-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+[n,m] = size(x) ;
+indold = zeros(1,m) ;
+if cod==0
+ d = transpose(sum(bsxfun(@power,bsxfun(@minus,x,mean(x)),2)));
+ d = sortrows( [transpose(1:m) d],2) ;
+ d = d((1+(0:1:g-1))*m/g,1) ;
+ c = x(:,d);
+else
+ c = init ;
+end
+for iter=1:300
+ dist = zeros(g,m) ;
+ for i=1:g
+ dist(i,:) = sum(bsxfun(@power,bsxfun(@minus,x,c(:,i)),2));
+ end
+ [rien,ind] = min(dist) ;
+ if isequal(ind,indold)
+ break ;
+ end
+ indold = ind ;
+ for i=1:g
+ lin = bsxfun(@eq,ind,i.*ones(1,m)) ;
+ h = x(:,lin) ;
+ c(:,i) = mean(h,2) ;
+ end
+end
+SqrtVariance = zeros(n,n,g) ;
+Weights = zeros(1,g) ;
+for i=1:g
+ temp = x(:,bsxfun(@eq,ind,i*ones(1,m))) ;
+ u = bsxfun(@minus,temp,mean(temp,2)); %temp-mean(temp,1)' ;
+ SqrtVariance(:,:,i) = chol( (u*u')/size(temp,2) )' ;
+ Weights(i) = size(temp,2)/m ;
+end
diff --git a/matlab/nonlinear-filters/neff.m b/matlab/nonlinear-filters/neff.m
new file mode 100644
index 0000000000000000000000000000000000000000..bd428fcd1bdea591f25433dcf387990eda8f7e9f
--- /dev/null
+++ b/matlab/nonlinear-filters/neff.m
@@ -0,0 +1,21 @@
+function n = neff(w)
+% Evaluates the criterion for resampling
+
+% Copyright © 2013-2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+n = 1/dot(w,w);
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/nonlinear_kalman_filter.m b/matlab/nonlinear-filters/nonlinear_kalman_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..ccae0731dbb3f7d654a89f305b319eef50a01842
--- /dev/null
+++ b/matlab/nonlinear-filters/nonlinear_kalman_filter.m
@@ -0,0 +1,182 @@
+function [LIK,lik] = nonlinear_kalman_filter(ReducedForm, Y, start, ParticleOptions, ThreadsOptions, DynareOptions, Model)
+
+% Evaluates the likelihood of a non-linear model approximating the
+% predictive (prior) and filtered (posterior) densities for state variables
+% by a Kalman filter.
+% Gaussian distribution approximation is done by:
+% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009).
+% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1995)
+% - Monte-Carlo draws from a multivariate gaussian distribution.
+% First and second moments of prior and posterior state densities are computed
+% from the resulting nodes/particles and allows to generate new distributions at the
+% following observation.
+% Pros: The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles
+% filters since it treats a lesser number of particles.
+% Cons: 1. Application a linear projection formulae in a nonlinear context.
+% 2. Parameter estimations may be biaised if the model is truly non-gaussian since predictive and
+% filtered densities are unimodal.
+%
+% INPUTS
+% Reduced_Form [structure] Matlab's structure describing the reduced form model.
+% Y [double] matrix of original observed variables.
+% start [double] structural parameters.
+% ParticleOptions [structure] Matlab's structure describing options concerning particle filtering.
+% ThreadsOptions [structure] Matlab's structure.
+%
+% OUTPUTS
+% LIK [double] scalar, likelihood
+% lik [double] vector, density of observations in each period.
+%
+% REFERENCES
+%
+% NOTES
+% The vector "lik" is used to evaluate the jacobian of the likelihood.
+
+% Copyright © 2009-2022 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+order = DynareOptions.order;
+
+if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ udr = ReducedForm.udr;
+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;
+ ghs2 = ReducedForm.ghs2;
+ if (order == 3)
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+end
+
+constant = ReducedForm.constant;
+steadystate = ReducedForm.steadystate;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+mf0 = ReducedForm.mf0;
+mf1 = ReducedForm.mf1;
+sample_size = size(Y,2);
+number_of_state_variables = length(mf0);
+number_of_observed_variables = length(mf1);
+number_of_structural_innovations = length(ReducedForm.Q);
+
+% compute gaussian quadrature nodes and weights on states and shocks
+if ParticleOptions.proposal_approximation.montecarlo
+ nodes = randn(ParticleOptions.number_of_particles,number_of_state_variables+number_of_structural_innovations);
+ weights = 1/ParticleOptions.number_of_particles;
+ weights_c = weights;
+elseif ParticleOptions.proposal_approximation.cubature
+ [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations);
+ weights_c = weights;
+elseif ParticleOptions.proposal_approximation.unscented
+ [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,ParticleOptions);
+else
+ error('Estimation: This approximation for the proposal is not implemented or unknown!')
+end
+
+if ParticleOptions.distribution_approximation.montecarlo
+ set_dynare_seed('default');
+end
+
+% Get covariance matrices
+H = ReducedForm.H;
+H_lower_triangular_cholesky = chol(H)' ;
+Q_lower_triangular_cholesky = chol(ReducedForm.Q)';
+
+% Get initial condition for the state vector.
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
+
+% Initialization of the likelihood.
+lik = NaN(sample_size,1);
+LIK = NaN;
+
+for t=1:sample_size
+ xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1) ] ;
+ sqr_Px = [StateVectorVarianceSquareRoot zeros(number_of_state_variables,number_of_structural_innovations);
+ zeros(number_of_structural_innovations,number_of_state_variables) Q_lower_triangular_cholesky];
+ sigma_points = bsxfun(@plus,xbar,sqr_Px*(nodes'));
+ StateVectors = sigma_points(1:number_of_state_variables,:);
+ epsilon = sigma_points(number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations,:);
+ 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, udr);
+ else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, false);
+ end
+ end
+ PredictedStateMean = tmp(mf0,:)*weights ;
+ PredictedObservedMean = tmp(mf1,:)*weights;
+ if ParticleOptions.proposal_approximation.cubature || ParticleOptions.proposal_approximation.montecarlo
+ PredictedStateMean = sum(PredictedStateMean,2);
+ PredictedObservedMean = sum(PredictedObservedMean,2);
+ dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean)'.*sqrt(weights);
+ dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean)'.*sqrt(weights);
+ big_mat = [dObserved dState ; [H_lower_triangular_cholesky zeros(number_of_observed_variables,number_of_state_variables)] ];
+ [~, mat] = qr2(big_mat,0);
+ mat = mat';
+ 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));
+ PredictionError = Y(:,t) - PredictedObservedMean;
+ StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*PredictionError;
+ else
+ dState = bsxfun(@minus,tmp(mf0,:),PredictedStateMean);
+ dObserved = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
+ PredictedStateVariance = dState*diag(weights_c)*dState';
+ PredictedObservedVariance = dObserved*diag(weights_c)*dObserved' + H;
+ PredictedStateAndObservedCovariance = dState*diag(weights_c)*dObserved';
+ PredictionError = Y(:,t) - PredictedObservedMean;
+ KalmanFilterGain = PredictedStateAndObservedCovariance/PredictedObservedVariance;
+ StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError;
+ StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
+ [StateVectorVarianceSquareRoot, p]= chol(StateVectorVariance,'lower');
+ if p
+ LIK=-Inf;
+ lik(t)=-Inf;
+ return
+ end
+ [~, p]= chol(PredictedObservedVariance,'lower');
+ if p
+ LIK=-Inf;
+ lik(t)=-Inf;
+ return
+ end
+ end
+ lik(t) = log( sum(probability2(Y(:,t),H_lower_triangular_cholesky,tmp(mf1,:)).*weights,1) ) ;
+end
+
+LIK = -sum(lik(start:end));
diff --git a/matlab/nonlinear-filters/online_auxiliary_filter.m b/matlab/nonlinear-filters/online_auxiliary_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..bbfb21e3e5f9c6b1059f9c93207d2584d4ae01bd
--- /dev/null
+++ b/matlab/nonlinear-filters/online_auxiliary_filter.m
@@ -0,0 +1,455 @@
+function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, DynareResults)
+
+% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
+%
+% INPUTS
+% - xparam1 [double] n×1 vector, Initial condition for the estimated parameters.
+% - DynareDataset [dseries] Sample used for estimation.
+% - dataset_info [struct] Description of the sample.
+% - DynareOptions [struct] Option values (options_).
+% - Model [struct] Description of the model (M_).
+% - EstimatedParameters [struct] Description of the estimated parameters (estim_params_).
+% - BayesInfo [struct] Prior definition (bayestopt_).
+% - DynareResults [struct] Results (oo_).
+%
+% OUTPUTS
+% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters).
+% - pmode [double] n×1 vector, mode of the particles at the end of the sample (for the parameters).
+% - pmedian [double] n×1 vector, median of the particles at the end of the sample (for the parameters).
+% - pstdev [double] n×1 vector, st. dev. of the particles at the end of the sample (for the parameters).
+% - p025 [double] n×1 vector, 2.5 percent of the particles are below p025(i) for i=1,…,n.
+% - p975 [double] n×1 vector, 97.5 percent of the particles are below p975(i) for i=1,…,n.
+% - covariance [double] n×n matrix, covariance of the particles at the end of the sample.
+
+% Copyright © 2013-2022 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% Set seed for randn().
+set_dynare_seed('default');
+pruning = DynareOptions.particle.pruning;
+second_resample = DynareOptions.particle.resampling.status.systematic;
+variance_update = true;
+
+bounds = prior_bounds(BayesInfo, DynareOptions.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
+
+% initialization of state particles
+[~, Model, DynareOptions, DynareResults, ReducedForm] = solve_model_for_online_filter(true, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
+
+order = DynareOptions.order;
+mf0 = ReducedForm.mf0;
+mf1 = ReducedForm.mf1;
+number_of_particles = DynareOptions.particle.number_of_particles;
+number_of_parameters = size(xparam1,1);
+Y = DynareDataset.data;
+sample_size = size(Y,1);
+number_of_observed_variables = length(mf1);
+number_of_structural_innovations = length(ReducedForm.Q);
+liu_west_delta = DynareOptions.particle.liu_west_delta;
+
+% Get initial conditions for the state particles
+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 order == 2
+ StateVectors_ = StateVectors;
+ elseif order == 3
+ StateVectors_ = repmat(StateVectors,3,1);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+end
+
+% parameters for the Liu & West filter
+small_a = (3*liu_west_delta-1)/(2*liu_west_delta);
+b_square = 1-small_a*small_a;
+
+% Initialization of parameter particles
+xparam = zeros(number_of_parameters,number_of_particles);
+Prior = dprior(BayesInfo, DynareOptions.prior_trunc);
+for i=1:number_of_particles
+ info = 12042009;
+ while info
+ candidate = Prior.draw();
+ [info, Model, DynareOptions, DynareResults] = solve_model_for_online_filter(false, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
+ if ~info
+ xparam(:,i) = candidate(:);
+ end
+ end
+end
+
+% Initialization of the weights of particles.
+weights = ones(1,number_of_particles)/number_of_particles;
+
+% Initialization of the likelihood.
+const_lik = log(2*pi)*number_of_observed_variables;
+mean_xparam = zeros(number_of_parameters,sample_size);
+mode_xparam = zeros(number_of_parameters,sample_size);
+median_xparam = zeros(number_of_parameters,sample_size);
+std_xparam = zeros(number_of_parameters,sample_size);
+lb95_xparam = zeros(number_of_parameters,sample_size);
+ub95_xparam = zeros(number_of_parameters,sample_size);
+
+%% The Online filter
+for t=1:sample_size
+ if t>1
+ fprintf('\nSubsample with %s first observations.\n\n', int2str(t))
+ else
+ fprintf('\nSubsample with only the first observation.\n\n')
+ end
+ % Moments of parameters particles distribution
+ m_bar = xparam*(weights');
+ temp = bsxfun(@minus,xparam,m_bar);
+ sigma_bar = (bsxfun(@times,weights,temp))*(temp');
+ if variance_update
+ chol_sigma_bar = chol(b_square*sigma_bar)';
+ end
+ % Prediction (without shocks)
+ fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam);
+ tau_tilde = zeros(1,number_of_particles);
+ for i=1:number_of_particles
+ % model resolution
+ [info, Model, DynareOptions, DynareResults, ReducedForm] = ...
+ solve_model_for_online_filter(false, fore_xparam(:,i), DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults);
+ if ~info(1)
+ 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;
+ udr = ReducedForm.udr;
+ else
+ steadystate = ReducedForm.steadystate;
+ 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;
+ ghs2 = ReducedForm.ghs2;
+ if (order == 3)
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+ if pruning
+ if order == 2
+ state_variables_steady_state_ = state_variables_steady_state;
+ elseif order == 3
+ state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+ end
+
+ 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, udr);
+ else
+ if pruning
+ yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state_);
+ if order == 2
+ [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);
+ elseif order == 3
+ [tmp, tmp_] = local_state_space_iteration_3(yhat_, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+ else
+ if order == 2
+ 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);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, zeros(number_of_structural_innovations, 1), ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+ 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) ;
+ tau_tilde(i) = weights(i).*(tpdf(z, 3*ones(size(z)))+1e-99) ;
+ end
+ end
+ % particles selection
+ tau_tilde = tau_tilde/sum(tau_tilde);
+ indx = resample(0, tau_tilde', DynareOptions.particle);
+ StateVectors = StateVectors(:,indx);
+ xparam = fore_xparam(:,indx);
+ if pruning
+ StateVectors_ = StateVectors_(:,indx);
+ end
+ w_stage1 = weights(indx)./tau_tilde(indx);
+ % draw in the new distributions
+ wtilde = zeros(1, number_of_particles);
+ for i=1:number_of_particles
+ info = 12042009;
+ counter=0;
+ while info(1) && counter <DynareOptions.particle.liu_west_max_resampling_tries
+ counter=counter+1;
+ candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
+ if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
+ % model resolution for new parameters particles
+ [info, Model, DynareOptions, DynareResults, ReducedForm] = ...
+ solve_model_for_online_filter(false, candidate, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults) ;
+ if ~info(1)
+ xparam(:,i) = candidate ;
+ 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;
+ udr = ReducedForm.udr;
+ 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;
+ ghs2 = ReducedForm.ghs2;
+ if (order == 3)
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+ if pruning
+ if order == 2
+ state_variables_steady_state_ = state_variables_steady_state;
+ mf0_ = mf0;
+ elseif order == 3
+ state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
+ mf0_ = repmat(mf0,1,3);
+ mask2 = number_of_state_variables+1:2*number_of_state_variables;
+ mask3 = 2*number_of_state_variables+1:number_of_state_variables;
+ mf0_(mask2) = mf0_(mask2)+size(ghx,1);
+ mf0_(mask3) = mf0_(mask3)+2*size(ghx,1);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+ end
+ 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, udr);
+ else
+ if pruning
+ yhat_ = bsxfun(@minus,StateVectors_(:,i), state_variables_steady_state_);
+ if order == 2
+ [tmp, tmp_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, yhat_, steadystate, DynareOptions.threads.local_state_space_iteration_2);
+ elseif order == 3
+ [tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+ StateVectors_(:,i) = tmp_(mf0_,:);
+ else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, DynareOptions.threads.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, DynareOptions.threads.local_state_space_iteration_3, pruning);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+ 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)));
+ end
+ end
+ if counter==DynareOptions.particle.liu_west_max_resampling_tries
+ fprintf('\nLiu & West particle filter: I haven''t been able to solve the model in %u tries.\n',DynareOptions.particle.liu_west_max_resampling_tries)
+ fprintf('Liu & West particle filter: The last error message was: %s\n',get_error_message(info))
+ fprintf('Liu & West particle filter: You can try to increase liu_west_max_resampling_tries, but most\n')
+ fprintf('Liu & West particle filter: likely there is an issue with the model.\n')
+ error('Liu & West particle filter: unable to solve the model.')
+ end
+ end
+ end
+ % normalization
+ weights = wtilde/sum(wtilde);
+ if variance_update && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size)
+ variance_update = false;
+ end
+ % final resampling (not advised)
+ if second_resample
+ [~, idmode] = max(weights);
+ mode_xparam(:,t) = xparam(:,idmode);
+ indx = resample(0, weights,DynareOptions.particle);
+ StateVectors = StateVectors(:,indx) ;
+ if pruning
+ StateVectors_ = StateVectors_(:,indx);
+ end
+ xparam = xparam(:,indx);
+ weights = ones(1, number_of_particles)/number_of_particles;
+ mean_xparam(:,t) = mean(xparam, 2);
+ mat_var_cov = bsxfun(@minus, xparam, mean_xparam(:,t));
+ mat_var_cov = (mat_var_cov*mat_var_cov')/(number_of_particles-1);
+ std_xparam(:,t) = sqrt(diag(mat_var_cov));
+ for i=1:number_of_parameters
+ temp = sortrows(xparam(i,:)');
+ lb95_xparam(i,t) = temp(0.025*number_of_particles);
+ median_xparam(i,t) = temp(0.5*number_of_particles);
+ ub95_xparam(i,t) = temp(0.975*number_of_particles);
+ end
+ end
+ if second_resample
+ [~, idmode] = max(weights);
+ mode_xparam(:,t) = xparam(:,idmode);
+ mean_xparam(:,t) = xparam*(weights');
+ mat_var_cov = bsxfun(@minus, xparam,mean_xparam(:,t));
+ mat_var_cov = mat_var_cov*(bsxfun(@times, mat_var_cov, weights)');
+ std_xparam(:,t) = sqrt(diag(mat_var_cov));
+ for i=1:number_of_parameters
+ temp = sortrows([xparam(i,:)' weights'], 1);
+ cumulated_weights = cumsum(temp(:,2));
+ pass1 = false;
+ pass2 = false;
+ pass3 = false;
+ for j=1:number_of_particles
+ if ~pass1 && cumulated_weights(j)>=0.025
+ lb95_xparam(i,t) = temp(j,1);
+ pass1 = true;
+ end
+ if ~pass2 && cumulated_weights(j)>=0.5
+ median_xparam(i,t) = temp(j,1);
+ pass2 = true;
+ end
+ if ~pass3 && cumulated_weights(j)>=0.975
+ ub95_xparam(i,t) = temp(j,1);
+ pass3 = true;
+ end
+ end
+ end
+ end
+ str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)');
+ for l=1:size(xparam,1)
+ str = sprintf('%s\n %5.4f \t\t %7.5f \t\t %5.4f', str, lb95_xparam(l,t), mean_xparam(l,t), ub95_xparam(l,t));
+ end
+ disp(str)
+ disp('')
+end
+
+pmean = xparam(:,sample_size);
+pmode = mode_xparam(:,sample_size);
+pstdev = std_xparam(:,sample_size) ;
+p025 = lb95_xparam(:,sample_size) ;
+p975 = ub95_xparam(:,sample_size) ;
+pmedian = median_xparam(:,sample_size) ;
+covariance = mat_var_cov;
+
+%% Plot parameters trajectory
+TeX = DynareOptions.TeX;
+
+nr = ceil(sqrt(number_of_parameters)) ;
+nc = floor(sqrt(number_of_parameters));
+nbplt = 1 ;
+
+if TeX
+ fidTeX = fopen([Model.fname '_param_traj.tex'],'w');
+ fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
+ fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+ fprintf(fidTeX,' \n');
+end
+
+for plt = 1:nbplt
+ hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories');
+ for k=1:length(pmean)
+ subplot(nr,nc,k)
+ [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
+ % Draw the surface for an interval containing 95% of the particles.
+ area(1:sample_size, ub95_xparam(k,:), 'FaceColor', [.9 .9 .9], 'BaseValue', min(lb95_xparam(k,:)));
+ hold on
+ area(1:sample_size, lb95_xparam(k,:), 'FaceColor', [1 1 1], 'BaseValue', min(lb95_xparam(k,:)));
+ % Draw the mean of particles.
+ plot(1:sample_size, mean_xparam(k,:), '-k', 'linewidth', 2)
+ if TeX
+ title(texname,'interpreter','latex')
+ else
+ title(name,'interpreter','none')
+ end
+ hold off
+ axis tight
+ drawnow
+ end
+ dyn_saveas(hh, [Model.fname '_param_traj' int2str(plt)], DynareOptions.nodisplay, DynareOptions.graph_format);
+ if TeX
+ % TeX eps loader file
+ fprintf(fidTeX,'\\begin{figure}[H]\n');
+ fprintf(fidTeX,'\\centering \n');
+ fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',Model.fname,int2str(plt));
+ fprintf(fidTeX,'\\caption{Parameters trajectories.}');
+ fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt));
+ fprintf(fidTeX,'\\end{figure}\n');
+ fprintf(fidTeX,' \n');
+ end
+end
+
+% Plot Parameter Densities
+number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
+bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
+kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
+for plt = 1:nbplt
+ hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities');
+ for k=1:length(pmean)
+ subplot(nr,nc,k)
+ [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
+ optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function);
+ [density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ...
+ number_of_particles, optimal_bandwidth, kernel_function);
+ plot(density(:,1), density(:,2));
+ hold on
+ if TeX
+ title(texname,'interpreter','latex')
+ else
+ title(name,'interpreter','none')
+ end
+ hold off
+ axis tight
+ drawnow
+ end
+ dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
+ if TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format)))
+ % TeX eps loader file
+ fprintf(fidTeX, '\\begin{figure}[H]\n');
+ fprintf(fidTeX,'\\centering \n');
+ fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/nc,1),M_.fname,int2str(plt));
+ fprintf(fidTeX,'\\caption{Parameter densities based on the Liu/West particle filter.}');
+ fprintf(fidTeX,'\\label{Fig:ParameterDensities:%s}\n',int2str(plt));
+ fprintf(fidTeX,'\\end{figure}\n');
+ fprintf(fidTeX,' \n');
+ end
+end
diff --git a/matlab/nonlinear-filters/probability.m b/matlab/nonlinear-filters/probability.m
new file mode 100644
index 0000000000000000000000000000000000000000..a934b48d8077b95cb0bbedb2d46a400d495ec652
--- /dev/null
+++ b/matlab/nonlinear-filters/probability.m
@@ -0,0 +1,39 @@
+function [prior,likelihood,C,posterior] = probability(mu,sqrtP,prior,X)
+
+% Copyright © 2013-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+[dim,nov] = size(X);
+M = size(mu,2) ;
+if nargout>1
+ likelihood = zeros(M,nov);
+ normfact = (2*pi)^(dim/2);
+ for k=1:M
+ XX = bsxfun(@minus,X,mu(:,k));
+ S = sqrtP(:,:,k);
+ foo = S \ XX;
+ likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S))));
+ end
+end
+likelihood = likelihood + 1e-99;
+if nargout>2
+ C = prior*likelihood + 1e-99;
+end
+if nargout>3
+ posterior = bsxfun(@rdivide,bsxfun(@times,prior',likelihood),C) + 1e-99 ;
+ posterior = bsxfun(@rdivide,posterior,sum(posterior,1));
+end
diff --git a/matlab/nonlinear-filters/probability2.m b/matlab/nonlinear-filters/probability2.m
new file mode 100644
index 0000000000000000000000000000000000000000..e37cfe3d101c9383712b6341a74567ee814516d4
--- /dev/null
+++ b/matlab/nonlinear-filters/probability2.m
@@ -0,0 +1,37 @@
+function [density] = probability2(mu,S,X)
+%
+% Multivariate gaussian density
+%
+% INPUTS
+% n [integer] scalar, number of variables.
+%
+% OUTPUTS
+% nodes [double] nodes of the cubature
+% weigths [double] associated weigths
+%
+% REFERENCES
+%
+%
+% NOTES
+%
+% Copyright © 2009-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+dim = size(X,1) ;
+normfact = bsxfun(@power,(2*pi),(dim/2)) ;
+foo = S\(bsxfun(@minus,X,mu)) ;
+density = exp(-0.5*sum(foo.*foo)')./abs((normfact*prod(diag(S)))) + 1e-99 ;
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/probability3.m b/matlab/nonlinear-filters/probability3.m
new file mode 100644
index 0000000000000000000000000000000000000000..18dc442b8bc32d9ddfe53e1b53ae523011aefb12
--- /dev/null
+++ b/matlab/nonlinear-filters/probability3.m
@@ -0,0 +1,39 @@
+function [prior,likelihood,C,posterior] = probability3(mu,sqrtP,prior,X,X_weights)
+
+% Copyright © 2013-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+[dim,nov] = size(X);
+M = size(mu,2) ;
+if nargout>1
+ likelihood = zeros(M,nov);
+ normfact = (2*pi)^(dim/2);
+ for k=1:M
+ XX = bsxfun(@minus,X,mu(:,k));
+ S = sqrtP(:,:,k);
+ foo = S \ XX;
+ likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S))));
+ end
+end
+wlikelihood = bsxfun(@times,X_weights,likelihood) + 1e-99;
+if nargout>2
+ C = prior*wlikelihood + 1e-99;
+end
+if nargout>3
+ posterior = bsxfun(@rdivide,bsxfun(@times,prior',wlikelihood),C) + 1e-99 ;
+ posterior = bsxfun(@rdivide,posterior,sum(posterior,1));
+end
diff --git a/matlab/nonlinear-filters/resample.m b/matlab/nonlinear-filters/resample.m
new file mode 100644
index 0000000000000000000000000000000000000000..fc50173c28eabac2df1b8aa4549233fbb7b8f357
--- /dev/null
+++ b/matlab/nonlinear-filters/resample.m
@@ -0,0 +1,79 @@
+function resampled_output = resample(particles,weights,ParticleOptions)
+% Resamples particles.
+% if particles = 0, returns the resampling index (except for smooth resampling)
+% Otherwise, returns the resampled particles set.
+
+%@info:
+%! @deftypefn {Function File} {@var{indx} =} resample (@var{weights}, @var{method})
+%! @anchor{particle/resample}
+%! @sp 1
+%! Resamples particles.
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item weights
+%! n*1 vector of doubles, particles' weights.
+%! @item method
+%! string equal to 'residual' or 'traditional'.
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item indx
+%! n*1 vector of intergers, indices.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 1
+%! @ref{particle/sequantial_importance_particle_filter}
+%! @sp 2
+%! @strong{This function calls:}
+%! @sp 1
+%! @ref{residual_resampling}, @ref{traditional_resampling}
+%! @sp 2
+%! @end deftypefn
+%@eod:
+
+% Copyright © 2011-2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+defaultmethod = 1; % For residual based method set this variable equal to 0.
+
+if defaultmethod
+ if ParticleOptions.resampling.method.kitagawa
+ resampled_output = traditional_resampling(particles,weights,rand);
+ elseif ParticleOptions.resampling.method.stratified
+ resampled_output = traditional_resampling(particles,weights,rand(size(weights)));
+ elseif ParticleOptions.resampling.method.smooth
+ if particles==0
+ error('Particle = 0 is incompatible with this resampling method!')
+ end
+ resampled_output = multivariate_smooth_resampling(particles,weights);
+ else
+ error('Unknown sampling method!')
+ end
+else
+ if ParticleOptions.resampling.method.kitagawa
+ resampled_output = residual_resampling(particles,weights,rand);
+ elseif ParticleOptions.resampling.method.stratified
+ resampled_output = residual_resampling(particles,weights,rand(size(weights)));
+ else
+ error('Unknown sampling method!')
+ end
+end
diff --git a/matlab/nonlinear-filters/residual_resampling.m b/matlab/nonlinear-filters/residual_resampling.m
new file mode 100644
index 0000000000000000000000000000000000000000..ee0ccf9ee2104370fbfbaf60404e7955b786fe98
--- /dev/null
+++ b/matlab/nonlinear-filters/residual_resampling.m
@@ -0,0 +1,143 @@
+function return_resample = residual_resampling(particles,weights,noise)
+% Resamples particles.
+
+%@info:
+%! @deftypefn {Function File} {@var{indx} =} residual_resampling (@var{weights})
+%! @anchor{particle/residual_resampling}
+%! @sp 1
+%! Resamples particles.
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item weights
+%! n*1 vector of doubles, particles' weights.
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item indx
+%! n*1 vector of intergers, indices.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 1
+%! @ref{particle/resample}
+%! @sp 2
+%! @strong{This function calls:}
+%! @sp 2
+%! @end deftypefn
+%@eod:
+
+% Copyright © 2011-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr
+% stephane DOT adjemian AT univ DASH lemans DOT fr
+
+% What is the number of particles?
+number_of_particles = length(weights);
+
+switch length(noise)
+ case 1
+ kitagawa_resampling = 1;
+ case number_of_particles
+ kitagawa_resampling = 0;
+ otherwise
+ error(['particle::resampling: Unknown method! The size of the second argument (' inputname(3) ') is wrong.'])
+end
+
+% Set vectors of indices.
+jndx = 1:number_of_particles;
+indx = zeros(1,number_of_particles);
+
+% Multiply the weights by the number of particles.
+WEIGHTS = number_of_particles*weights;
+
+% Compute the integer part of the normalized weights.
+iWEIGHTS = fix(WEIGHTS);
+
+% Compute the number of resample
+number_of_trials = number_of_particles-sum(iWEIGHTS);
+
+if number_of_trials
+ WEIGHTS = (WEIGHTS-iWEIGHTS)/number_of_trials;
+ EmpiricalCDF = cumsum(WEIGHTS);
+ if kitagawa_resampling
+ u = (transpose(1:number_of_trials)-1+noise(:))/number_of_trials;
+ else
+ u = fliplr(cumprod(noise(1:number_of_trials).^(1./(number_of_trials:-1:1))));
+ end
+ j=1;
+ for i=1:number_of_trials
+ while (u(i)>EmpiricalCDF(j))
+ j=j+1;
+ end
+ iWEIGHTS(j)=iWEIGHTS(j)+1;
+ if kitagawa_resampling==0
+ j=1;
+ end
+ end
+end
+
+k=1;
+for i=1:number_of_particles
+ if (iWEIGHTS(i)>0)
+ for j=k:k+iWEIGHTS(i)-1
+ indx(j) = jndx(i);
+ end
+ end
+ k = k + iWEIGHTS(i);
+end
+
+if particles==0
+ return_resample = indx ;
+else
+ return_resample = particles(indx,:) ;
+end
+%@test:1
+%$ % Define the weights
+%$ weights = randn(2000,1).^2;
+%$ weights = weights/sum(weights);
+%$ % Initialize t.
+%$ t = ones(1,1);
+%$
+%$ try
+%$ indx1 = residual_resampling(weights);
+%$ catch
+%$ t(1) = 0;
+%$ end
+%$
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ % Define the weights
+%$ weights = exp(randn(2000,1));
+%$ weights = weights/sum(weights);
+%$ % Initialize t.
+%$ t = ones(1,1);
+%$
+%$ try
+%$ indx1 = residual_resampling(weights);
+%$ catch
+%$ t(1) = 0;
+%$ end
+%$
+%$ T = all(t);
+%@eof:2
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/sequential_importance_particle_filter.m b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..c69ded02f27f1d1d5eca6d7659b6962e84646089
--- /dev/null
+++ b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
@@ -0,0 +1,185 @@
+function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions, DynareOptions, Model)
+
+% Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
+
+% Copyright © 2011-2022 Dynare Team
+%
+% This file is part of Dynare (particles module).
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare particles module is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+persistent init_flag
+persistent mf0 mf1
+persistent number_of_particles number_of_state_variables
+persistent sample_size number_of_observed_variables number_of_structural_innovations
+
+% Set default value for start
+if isempty(start)
+ start = 1;
+end
+
+% Set flag for prunning
+pruning = ParticleOptions.pruning;
+
+% Get steady state and mean.
+steadystate = ReducedForm.steadystate;
+constant = ReducedForm.constant;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+order = DynareOptions.order;
+
+% Set persistent variables (if needed).
+if isempty(init_flag)
+ mf0 = ReducedForm.mf0;
+ mf1 = ReducedForm.mf1;
+ sample_size = size(Y,2);
+ number_of_state_variables = length(mf0);
+ number_of_observed_variables = length(mf1);
+ number_of_structural_innovations = length(ReducedForm.Q);
+ number_of_particles = ParticleOptions.number_of_particles;
+ init_flag = 1;
+end
+
+if ReducedForm.use_k_order_solver
+ dr = ReducedForm.dr;
+ udr = ReducedForm.udr;
+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;
+ ghs2 = ReducedForm.ghs2;
+ if order == 3
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
+end
+
+% Get covariance matrices.
+Q = ReducedForm.Q; % Covariance matrix of the structural innovations.
+H = ReducedForm.H; % Covariance matrix of the measurement errors.
+if isempty(H)
+ H = 0;
+end
+
+% Initialization of the likelihood.
+const_lik = log(2*pi)*number_of_observed_variables +log(det(H)) ;
+lik = NaN(sample_size,1);
+
+% Get initial condition for the state vector.
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+if pruning
+ StateVectorMean_ = StateVectorMean;
+ StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
+end
+
+% Get the rank of StateVectorVarianceSquareRoot
+state_variance_rank = size(StateVectorVarianceSquareRoot,2);
+
+% Factorize the covariance matrix of the structural innovations
+Q_lower_triangular_cholesky = chol(Q)';
+
+% Set seed for randn().
+set_dynare_seed('default');
+
+% Initialization of the weights across particles.
+weights = ones(1,number_of_particles)/number_of_particles ;
+StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
+if pruning
+ if order == 2
+ StateVectors_ = StateVectors;
+ state_variables_steady_state_ = state_variables_steady_state;
+ mf0_ = mf0;
+ elseif order == 3
+ StateVectors_ = repmat(StateVectors,3,1);
+ state_variables_steady_state_ = repmat(state_variables_steady_state,3,1);
+ mf0_ = repmat(mf0,1,3);
+ mask2 = number_of_state_variables+1:2*number_of_state_variables;
+ mask3 = 2*number_of_state_variables+1:3*number_of_state_variables;
+ mf0_(mask2) = mf0_(mask2)+size(ghx,1);
+ mf0_(mask3) = mf0_(mask3)+2*size(ghx,1);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+end
+
+% Loop over observations
+for t=1:sample_size
+ yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
+ epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
+ if pruning
+ yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state_);
+ if order == 2
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ [tmp, tmp_] = local_state_space_iteration_3(yhat_, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
+ else
+ error('Pruning is not available for orders > 3');
+ end
+ else
+ if ReducedForm.use_k_order_solver
+ tmp = local_state_space_iteration_k(yhat, epsilon, dr, Model, DynareOptions, udr);
+ else
+ if order == 2
+ tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
+ elseif order == 3
+ tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
+ else
+ error('Order > 3: use_k_order_solver should be set to true');
+ end
+ end
+ end
+ %PredictedObservedMean = tmp(mf1,:)*transpose(weights);
+ PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
+ %dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
+ %PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H;
+ %PredictedObservedVariance = H;
+ if rcond(H) > 1e-16
+ lnw = -.5*(const_lik+sum(PredictionError.*(H\PredictionError),1));
+ else
+ LIK = NaN;
+ return
+ end
+ dfac = max(lnw);
+ wtilde = weights.*exp(lnw-dfac);
+ lik(t) = log(sum(wtilde))+dfac;
+ weights = wtilde/sum(wtilde);
+ if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
+ if pruning
+ temp = resample([tmp(mf0,:)' tmp_(mf0_,:)'],weights',ParticleOptions);
+ StateVectors = temp(:,1:number_of_state_variables)';
+ StateVectors_ = temp(:,number_of_state_variables+1:end)';
+ else
+ StateVectors = resample(tmp(mf0,:)',weights',ParticleOptions)';
+ end
+ weights = ones(1,number_of_particles)/number_of_particles;
+ elseif ParticleOptions.resampling.status.none
+ StateVectors = tmp(mf0,:);
+ if pruning
+ StateVectors_ = tmp_(mf0_,:);
+ end
+ end
+end
+
+LIK = -sum(lik(start:end));
diff --git a/matlab/nonlinear-filters/solve_model_for_online_filter.m b/matlab/nonlinear-filters/solve_model_for_online_filter.m
new file mode 100644
index 0000000000000000000000000000000000000000..92cbf4c74042ee8d729a069e6fb7241796d58953
--- /dev/null
+++ b/matlab/nonlinear-filters/solve_model_for_online_filter.m
@@ -0,0 +1,215 @@
+function [info, Model, DynareOptions, DynareResults, ReducedForm] = ...
+ solve_model_for_online_filter(setinitialcondition, xparam1, DynareDataset, DynareOptions, Model, EstimatedParameters, BayesInfo, bounds, DynareResults)
+
+% Solves the dsge model for an particular parameters set.
+%
+% INPUTS
+% - setinitialcondition [logical] return initial condition if true.
+% - xparam1 [double] n×1 vector, parameter values.
+% - DynareDataset [struct] Dataset for estimation (dataset_).
+% - DynareOptions [struct] Dynare options (options_).
+% - Model [struct] Model description (M_).
+% - EstimatedParameters [struct] Estimated parameters (estim_params_).
+% - BayesInfo [struct] Prior definition (bayestopt_).
+% - DynareResults [struct] Dynare results (oo_).
+%
+% OUTPUTS
+% - info [integer] scalar, nonzero if any problem occur when computing the reduced form.
+% - Model [struct] Model description (M_).
+% - DynareOptions [struct] Dynare options (options_).
+% - DynareResults [struct] Dynare results (oo_).
+% - ReducedForm [struct] Reduced form model.
+
+% Copyright © 2013-2022 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+persistent init_flag restrict_variables_idx state_variables_idx mf0 mf1 number_of_state_variables
+
+info = 0;
+
+%----------------------------------------------------
+% 1. Get the structural parameters & define penalties
+%----------------------------------------------------
+
+% Test if some parameters are smaller than the lower bound of the prior domain.
+if any(xparam1<bounds.lb)
+ info = 41;
+ return
+end
+
+% Test if some parameters are greater than the upper bound of the prior domain.
+if any(xparam1>bounds.ub)
+ info = 42;
+ return
+end
+
+% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
+Q = Model.Sigma_e;
+H = Model.H;
+for i=1:EstimatedParameters.nvx
+ k =EstimatedParameters.var_exo(i,1);
+ Q(k,k) = xparam1(i)*xparam1(i);
+end
+offset = EstimatedParameters.nvx;
+if EstimatedParameters.nvn
+ for i=1:EstimatedParameters.nvn
+ H(i,i) = xparam1(i+offset)*xparam1(i+offset);
+ end
+ offset = offset+EstimatedParameters.nvn;
+else
+ H = zeros(size(DynareDataset.data, 2));
+end
+
+% Get the off-diagonal elements of the covariance matrix for the structural innovations. Test if Q is positive definite.
+if EstimatedParameters.ncx
+ for i=1:EstimatedParameters.ncx
+ k1 =EstimatedParameters.corrx(i,1);
+ k2 =EstimatedParameters.corrx(i,2);
+ Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
+ Q(k2,k1) = Q(k1,k2);
+ end
+ % Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
+ [~, testQ] = chol(Q);
+ if testQ
+ % The variance-covariance matrix of the structural innovations is not definite positive.
+ info = 43;
+ return
+ end
+ offset = offset+EstimatedParameters.ncx;
+end
+
+% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
+if EstimatedParameters.ncn
+ corrn_observable_correspondence = EstimatedParameters.corrn_observable_correspondence;
+ for i=1:EstimatedParameters.ncn
+ k1 = corrn_observable_correspondence(i,1);
+ k2 = corrn_observable_correspondence(i,2);
+ H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
+ H(k2,k1) = H(k1,k2);
+ end
+ % Try to compute the cholesky decomposition of H (possible iff H is positive definite)
+ [~, testH] = chol(H);
+ if testH
+ % The variance-covariance matrix of the measurement errors is not definite positive.
+ info = 44;
+ return
+ end
+ offset = offset+EstimatedParameters.ncn;
+end
+
+% Update estimated structural parameters in Mode.params.
+if EstimatedParameters.np > 0
+ Model.params(EstimatedParameters.param_vals(:,1)) = xparam1(offset+1:end);
+end
+
+% Update Model.Sigma_e and Model.H.
+Model.Sigma_e = Q;
+Model.H = H;
+
+%------------------------------------------------------------------------------
+% 2. call model setup & reduction program
+%------------------------------------------------------------------------------
+
+warning('off', 'MATLAB:nearlySingularMatrix')
+[~, ~, ~, info, Model, DynareResults] = ...
+ dynare_resolve(Model, DynareOptions, DynareResults, 'restrict');
+warning('on', 'MATLAB:nearlySingularMatrix')
+
+if info(1)~=0
+ if nargout==5
+ ReducedForm = 0;
+ end
+ return
+end
+
+% Get decision rules and transition equations.
+dr = DynareResults.dr;
+
+% Set persistent variables (first call).
+if isempty(init_flag)
+ mf0 = BayesInfo.mf0;
+ mf1 = BayesInfo.mf1;
+ restrict_variables_idx = dr.restrict_var_list;
+ state_variables_idx = restrict_variables_idx(mf0);
+ number_of_state_variables = length(mf0);
+ init_flag = true;
+end
+
+
+% Return reduced form model.
+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==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);
+ ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
+ elseif DynareOptions.order>=3
+ ReducedForm.use_k_order_solver = true;
+ ReducedForm.dr = dr;
+ ReducedForm.udr = folded_to_unfolded_dr(dr, Model, DynareOptions);
+ else
+ n_states=size(dr.ghx,2);
+ n_shocks=size(dr.ghu,2);
+ ReducedForm.use_k_order_solver = false;
+ ReducedForm.ghxx = zeros(size(restrict_variables_idx,1),n_states^2);
+ ReducedForm.ghuu = zeros(size(restrict_variables_idx,1),n_shocks^2);
+ ReducedForm.ghxu = zeros(size(restrict_variables_idx,1),n_states*n_shocks);
+ ReducedForm.constant = ReducedForm.steadystate;
+ end
+ ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
+ ReducedForm.Q = Q;
+ ReducedForm.H = H;
+ ReducedForm.mf0 = mf0;
+ ReducedForm.mf1 = mf1;
+end
+
+% Set initial condition
+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.state_variables_steady_state;%.constant(mf0);
+ [A,B] = kalman_transition_matrix(dr,dr.restrict_var_list,dr.restrict_columns,Model.exo_nbr);
+ StateVectorVariance2 = 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);
+ StateVectorVariance = lyapunov_symm(A, B*ReducedForm.Q*B', DynareOptions.lyapunov_fixed_point_tol, ...
+ DynareOptions.qz_criterium, DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
+ StateVectorVariance = StateVectorVariance(mf0,mf0);
+ 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.state_variables_steady_state;%.constant(mf0);
+ old_DynareOptionsperiods = DynareOptions.periods;
+ DynareOptions.periods = 5000;
+ old_DynareOptionspruning = DynareOptions.pruning;
+ DynareOptions.pruning = DynareOptions.particle.pruning;
+ y_ = simult(dr.ys, dr, Model, DynareOptions, DynareResults);
+ y_ = y_(dr.order_var(state_variables_idx),2001:DynareOptions.periods);
+ StateVectorVariance = cov(y_');
+ DynareOptions.periods = old_DynareOptionsperiods;
+ DynareOptions.pruning = old_DynareOptionspruning;
+ clear('old_DynareOptionsperiods','y_');
+ case 3% Initial state vector covariance is a diagonal matrix.
+ 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!')
+ end
+ ReducedForm.StateVectorMean = StateVectorMean;
+ ReducedForm.StateVectorVariance = StateVectorVariance;
+end
diff --git a/matlab/nonlinear-filters/spherical_radial_sigma_points.m b/matlab/nonlinear-filters/spherical_radial_sigma_points.m
new file mode 100644
index 0000000000000000000000000000000000000000..756cac28f0fa862a55cdfdfcc903e18dbaa032a9
--- /dev/null
+++ b/matlab/nonlinear-filters/spherical_radial_sigma_points.m
@@ -0,0 +1,36 @@
+function [nodes,weights] = spherical_radial_sigma_points(n)
+%
+% Computes nodes and weigths from a third-degree spherical-radial cubature
+% rule.
+% INPUTS
+% n [integer] scalar, number of variables.
+%
+% OUTPUTS
+% nodes [double] nodes of the cubature
+% weigths [double] associated weigths
+%
+% REFERENCES
+%
+% Arasaratnam & Haykin 2008,2009.
+%
+% NOTES
+%
+% Copyright © 2009-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+nodes = (sqrt(n)*([eye(n) -eye(n)]))' ;
+weights = (1/(2*n)) ;
diff --git a/matlab/nonlinear-filters/traditional_resampling.m b/matlab/nonlinear-filters/traditional_resampling.m
new file mode 100644
index 0000000000000000000000000000000000000000..ab9089c8c6b92b1defcd05b177b4b64393d73ed9
--- /dev/null
+++ b/matlab/nonlinear-filters/traditional_resampling.m
@@ -0,0 +1,237 @@
+function return_resample = traditional_resampling(particles,weights,noise)
+% Resamples particles.
+
+%@info:
+%! @deftypefn {Function File} {@var{indx} =} traditional_resampling (@var{weights},@var{noise})
+%! @anchor{particle/traditional_resampling}
+%! @sp 1
+%! Resamples particles (Resampling à la Kitagawa or stratified resampling).
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item weights
+%! n*1 vector of doubles, particles' weights.
+%! @item noise
+%! n*1 vector of doubles sampled from a [0,1] uniform distribution (stratified resampling) or scalar double
+%! sampled from a [0,1] uniform distribution (Kitagawa resampling).
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item indx
+%! n*1 vector of intergers, indices.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 1
+%! @ref{particle/resample}
+%! @sp 2
+%! @strong{This function calls:}
+%! @sp 2
+%! @end deftypefn
+%@eod:
+
+% Copyright © 2011-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr
+% stephane DOT adjemian AT univ DASH lemans DOT fr
+
+% What is the number of particles?
+number_of_particles = length(weights);
+
+% Initialize the returned argument.
+indx = ones(number_of_particles,1);
+
+% Select method.
+switch length(noise)
+ case 1
+ kitagawa_resampling = 1;
+ case number_of_particles
+ kitagawa_resampling = 0;
+ otherwise
+ error(['particle::resampling: Unknown method! The size of the second argument (' inputname(3) ') is wrong.'])
+end
+
+% Get the empirical CDF.
+c = cumsum(weights);
+
+% Draw a starting point.
+if kitagawa_resampling
+ randvec = (transpose(1:number_of_particles)-1+noise(:))/number_of_particles ;
+else
+ randvec = fliplr(cumprod(noise.^(1./(number_of_particles:-1:1))));
+end
+
+% Start at the bottom of the CDF
+if kitagawa_resampling
+ j = 1;
+ for i=1:number_of_particles
+ while (randvec(i)>c(j))
+ j = j+1;
+ end
+ indx(i) = j;
+ end
+else
+ for i=1:number_of_particles
+ indx(i) = sum(randvec(i)>c);
+ end
+ % Matlab's indices start at 1...
+ indx = indx+1;
+end
+
+if particles==0
+ return_resample = indx ;
+else
+ return_resample = particles(indx,:) ;
+end
+
+%@test:1
+%$ % Define the weights
+%$ weights = randn(2000,1).^2;
+%$ weights = weights/sum(weights);
+%$ % Initialize t.
+%$ t = ones(2,1);
+%$
+%$ % First, try the stratified resampling.
+%$ try
+%$ indx1 = traditional_resampling(weights,rand(2000,1));
+%$ catch
+%$ t(1) = 0;
+%$ end
+%$
+%$ % Second, try the Kitagawa resampling.
+%$ try
+%$ indx2 = traditional_resampling(weights,rand);
+%$ catch
+%$ t(2) = 0;
+%$ end
+%$
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ % Define the weights
+%$ weights = exp(randn(2000,1));
+%$ weights = weights/sum(weights);
+%$ % Initialize t.
+%$ t = ones(2,1);
+%$
+%$ % First, try the stratified resampling.
+%$ try
+%$ indx1 = traditional_resampling(weights,rand(2000,1));
+%$ catch
+%$ t(1) = 0;
+%$ end
+%$
+%$ % Second, try the Kitagawa resampling.
+%$ try
+%$ indx2 = traditional_resampling(weights,rand);
+%$ catch
+%$ t(2) = 0;
+%$ end
+%$
+%$ T = all(t);
+%@eof:2
+
+%@test:3
+%$ % Set the number of particles.
+%$ number_of_particles = 20000;
+%$
+%$ show_plot = 0;
+%$ show_time = 1;
+%$
+%$ % Define the weights
+%$ weights = randn(number_of_particles,1).^2;
+%$ weights = weights/sum(weights);
+%$
+%$ % Compute the empirical CDF
+%$ c = cumsum(weights);
+%$
+%$ % Stratified resampling.
+%$ noise = rand(number_of_particles,1);
+%$
+%$ if show_time
+%$ disp('Stratified resampling timing:')
+%$ tic
+%$ end
+%$
+%$ indx1 = traditional_resampling(weights,noise);
+%$
+%$ if show_time
+%$ toc
+%$ tic
+%$ end
+%$
+%$ indx1_ = zeros(number_of_particles,1);
+%$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles;
+%$ for i=1:number_of_particles
+%$ j = 1;
+%$ while (randvec(i)>c(j))
+%$ j = j + 1;
+%$ end
+%$ indx1_(i) = j;
+%$ end
+%$
+%$ if show_time
+%$ toc
+%$ end
+%$
+%$ % Kitagawa's resampling.
+%$ noise = rand;
+%$
+%$ if show_time
+%$ disp('Kitagawa''s resampling timing:')
+%$ tic
+%$ end
+%$
+%$ indx2 = traditional_resampling(weights,noise);
+%$
+%$ if show_time
+%$ toc
+%$ tic
+%$ end
+%$
+%$ indx2_ = zeros(number_of_particles,1);
+%$ randvec = (transpose(1:number_of_particles)-1+noise)/number_of_particles;
+%$ j = 1;
+%$ for i=1:number_of_particles
+%$ while (randvec(i)>c(j))
+%$ j = j + 1;
+%$ end
+%$ indx2_(i) = j;
+%$ end
+%$
+%$ if show_time
+%$ toc
+%$ end
+%$
+%$ % REMARK
+%$ % Note that the alternative code used in this test is sensibly faster than the code proposed
+%$ % in the routine for the resampling scheme à la Kitagawa...
+%$
+%$ if show_plot
+%$ plot(randvec,c,'-r'), hold on, plot([randvec(1),randvec(end)],[c(1),c(end)],'-k'), hold off, axis tight, box on
+%$ end
+%$
+%$ % Check results.
+%$ t(1) = dassert(indx1,indx1_);
+%$ t(2) = dassert(indx2,indx2_);
+%$ T = all(t);
+%@eof:3
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/univariate_smooth_resampling.m b/matlab/nonlinear-filters/univariate_smooth_resampling.m
new file mode 100644
index 0000000000000000000000000000000000000000..6661d8f0d6209d4838cf63e00ec266a927145cc0
--- /dev/null
+++ b/matlab/nonlinear-filters/univariate_smooth_resampling.m
@@ -0,0 +1,82 @@
+function new_particles = univariate_smooth_resampling(weights,particles,number_of_new_particles)
+% Smooth Resampling of the particles.
+
+%@info:
+%! @deftypefn {Function File} {@var{new_particles} =} univariate_smooth_resampling (@var{weights}, @var{number_of_new_particles})
+%! @anchor{particle/univariate_smooth_resampling}
+%! @sp 1
+%! Smooth Resampling of the particles (univariate version).
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item weights
+%! n*1 vector of doubles, particles' weights.
+%! @item particles
+%! n*1 vector of doubles, particles.
+%! @item number_of_new_particles
+%! Integer scalar.
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item indx
+%! number_of_new_particles*1 vector of doubles, new particles.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 1
+%! @ref{particle/sequantial_importance_particle_filter}
+%! @sp 2
+%! @strong{This function calls:}
+%! @sp 2
+%! @end deftypefn
+%@eod:
+
+% Copyright © 2012-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr
+% stephane DOT adjemian AT univ DASH lemans DOT fr
+
+M = length(particles) ;
+lambda_tilde = [ (.5*weights(1)) ;
+ (.5*(weights(1:M-1)+weights(2:M))) ;
+ (.5*weights(M)) ] ;
+lambda_bar = cumsum(lambda_tilde) ;
+u = rand(1,1) ;
+new_particles = zeros(number_of_new_particles,1) ;
+rj = 0 ;
+i = 1 ;
+j = 1 ;
+while i<=number_of_new_particles
+ u_j = ( i-1 + u)/number_of_new_particles ;
+ while u_j>lambda_bar(j)
+ rj = j ;
+ j = j+1 ;
+ end
+ if rj==0
+ new_particles(i) = particles(1) ;
+ elseif rj==M
+ new_particles(i) = particles(M) ;
+ else
+ u_star = (u_j - lambda_bar(rj))./lambda_tilde(rj+1) ;
+ new_particles(i) = (particles(rj+1) - particles(rj))*u_star + particles(rj) ;
+ end
+ i = i+1 ;
+end
\ No newline at end of file
diff --git a/matlab/nonlinear-filters/unscented_sigma_points.m b/matlab/nonlinear-filters/unscented_sigma_points.m
new file mode 100644
index 0000000000000000000000000000000000000000..9ff7a27e8c3add343078e367ba7ddc5be447075f
--- /dev/null
+++ b/matlab/nonlinear-filters/unscented_sigma_points.m
@@ -0,0 +1,40 @@
+function [nodes,W_m,W_c] = unscented_sigma_points(n,ParticleOptions)
+%
+% Computes nodes and weigths for a scaled unscented transform cubature
+% INPUTS
+% n [integer] scalar, number of variables.
+%
+% OUTPUTS
+% nodes [double] nodes of the cubature
+% weigths [double] associated weigths
+%
+% REFERENCES
+%
+%
+%
+% NOTES
+%
+% Copyright © 2009-2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+lambda = (ParticleOptions.unscented.alpha^2)*(n+ParticleOptions.unscented.kappa) - n ;
+nodes = [ zeros(n,1) ( sqrt(n+lambda).*([ eye(n) -eye(n)]) ) ]' ;
+W_m = lambda/(n+lambda) ;
+W_c = W_m + (1-ParticleOptions.unscented.alpha^2+ParticleOptions.unscented.beta) ;
+temp = ones(2*n,1)/(2*(n+lambda)) ;
+W_m = [W_m ; temp] ;
+W_c = [W_c ; temp] ;
diff --git a/matlab/particles b/matlab/particles
deleted file mode 160000
index 4a679b26eef0463169705cde460079feddc14819..0000000000000000000000000000000000000000
--- a/matlab/particles
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 4a679b26eef0463169705cde460079feddc14819