diff --git a/nonlinear-filters/src/auxiliary_initialization.m b/nonlinear-filters/src/auxiliary_initialization.m
new file mode 100644
index 000000000..784550000
--- /dev/null
+++ b/nonlinear-filters/src/auxiliary_initialization.m
@@ -0,0 +1,116 @@
+function initial_distribution = auxiliary_initialization(ReducedForm,Y,start,DynareOptions)
+% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
+%
+% INPUTS
+% ReducedForm [structure] Matlab's structure describing the reduced form model.
+% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors.
+% ReducedForm.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'.
+% mf [integer] pp*1 vector of indices.
+% number_of_particles [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 (C) 2013 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 .
+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 = DynareOptions.particle.pruning;
+
+% Get steady state and mean.
+%steadystate = ReducedForm.steadystate;
+constant = ReducedForm.constant;
+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 = DynareOptions.particle.number_of_particles;
+ init_flag = 1;
+end
+
+% 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;
+
+% 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,DynareOptions.threads.local_state_space_iteration_2);
+%else
+ tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+%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',DynareOptions)' ;
\ No newline at end of file
diff --git a/nonlinear-filters/src/auxiliary_particle_filter.m b/nonlinear-filters/src/auxiliary_particle_filter.m
new file mode 100644
index 000000000..e869484fe
--- /dev/null
+++ b/nonlinear-filters/src/auxiliary_particle_filter.m
@@ -0,0 +1,159 @@
+function [LIK,lik] = auxiliary_particle_filter(ReducedForm,Y,start,DynareOptions)
+% Evaluates the likelihood of a nonlinear model with a particle filter allowing eventually resampling.
+%
+% INPUTS
+% ReducedForm [structure] Matlab's structure describing the reduced form model.
+% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors.
+% ReducedForm.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'.
+% mf [integer] pp*1 vector of indices.
+% number_of_particles [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 (C) 2011-2013 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 .
+persistent init_flag mf0 mf1 number_of_particles
+persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+% Set flag for prunning
+pruning = DynareOptions.particle.pruning;
+
+% Get steady state and mean.
+steadystate = ReducedForm.steadystate;
+constant = ReducedForm.constant;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+% Set persistent variables.
+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 = DynareOptions.particle.number_of_particles;
+ init_flag = 1;
+end
+
+% 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;
+
+% 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;
+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);
+if pruning
+ StateVectors_ = StateVectors;
+end
+for t=1:sample_size
+ 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,DynareOptions.threads.local_state_space_iteration_2);
+ else
+ tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+ 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 ;
+ sum_tau_tilde = sum(tau_tilde) ;
+ %var_wtilde = wtilde-sum_tau_tilde ;
+ %var_wtilde = var_wtilde'*var_wtilde/(number_of_particles-1) ;
+ lik(t) = log(sum_tau_tilde) ; %+ .5*var_wtilde/(number_of_particles*(sum_tau_tilde*sum_tau_tilde)) ;
+ tau_tilde = tau_tilde/sum_tau_tilde;
+ if pruning
+ temp = resample([yhat' yhat_'],tau_tilde',DynareOptions);
+ yhat = temp(:,1:number_of_state_variables)' ;
+ yhat_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)' ;
+ else
+ yhat = resample(yhat',tau_tilde',DynareOptions)' ;
+ end
+ if pruning
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
+ else
+ tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+ 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))) ;
+ epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
+ if pruning
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
+ StateVectors_ = tmp_(mf0,:);
+ else
+ tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+ end
+ StateVectors = tmp(mf0,:);
+ PredictedObservedMean = mean(tmp(mf1,:),2);
+ PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
+ dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
+ PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles + H;
+ lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)));
+ wtilde = lnw./wtilde;
+ weights = wtilde/sum(wtilde);
+end
+
+LIK = -sum(lik(start:end));
\ No newline at end of file
diff --git a/nonlinear-filters/src/conditional_filter_proposal.m b/nonlinear-filters/src/conditional_filter_proposal.m
new file mode 100644
index 000000000..aa82a449a
--- /dev/null
+++ b/nonlinear-filters/src/conditional_filter_proposal.m
@@ -0,0 +1,123 @@
+function [ProposalStateVector,Weights] = conditional_filter_proposal(ReducedForm,obs,StateVectors,SampleWeights,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions,normconst2)
+%
+% Computes the proposal for each past particle using Gaussian approximations
+% for the state errors and the Kalman filter
+%
+% 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 (C) 2012-2013 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 .
+%
+% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr
+% stephane DOT adjemian AT univ DASH lemans DOT fr
+
+persistent init_flag2 mf0 mf1
+persistent number_of_state_variables number_of_observed_variables
+persistent number_of_structural_innovations
+
+% Set local state space model (first-order approximation).
+ghx = ReducedForm.ghx;
+ghu = ReducedForm.ghu;
+% Set local state space model (second-order approximation).
+ghxx = ReducedForm.ghxx;
+ghuu = ReducedForm.ghuu;
+ghxu = ReducedForm.ghxu;
+
+if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ...
+ any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ...
+ any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4))
+ ghx
+ ghu
+ ghxx
+ ghuu
+ ghxu
+end
+
+constant = ReducedForm.constant;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+% Set persistent variables.
+if isempty(init_flag2)
+ mf0 = ReducedForm.mf0;
+ mf1 = ReducedForm.mf1;
+ number_of_state_variables = length(mf0);
+ number_of_observed_variables = length(mf1);
+ number_of_structural_innovations = length(ReducedForm.Q);
+ init_flag2 = 1;
+end
+
+if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.proposal_approximation.montecarlo
+ [nodes,weights] = spherical_radial_sigma_points(number_of_structural_innovations);
+ weights_c = weights ;
+elseif DynareOptions.particle.proposal_approximation.unscented
+ [nodes,weights,weights_c] = unscented_sigma_points(number_of_structural_innovations,DynareOptions);
+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)) ;
+
+tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+
+PredictedStateMean = tmp(mf0,:)*weights ;
+PredictedObservedMean = tmp(mf1,:)*weights;
+
+if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.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)] ];
+ [mat1,mat] = qr2(big_mat,0);
+ mat = mat';
+ clear('mat1');
+ PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables);
+ CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
+ StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
+ StateVectorMean = PredictedStateMean + (CovarianceObservedStateSquareRoot/PredictedObservedVarianceSquareRoot)*(obs - PredictedObservedMean);
+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 ;
+ StateVectorMean = PredictedStateMean + KalmanFilterGain*(obs - PredictedObservedMean);
+ StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
+ StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';
+end
+
+ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ;
+ypred = measurement_equations(ProposalStateVector,ReducedForm,DynareOptions) ;
+foo = H_lower_triangular_cholesky \ (obs - ypred) ;
+likelihood = exp(-0.5*(foo)'*foo)/normconst2 + 1e-99 ;
+Weights = SampleWeights.*likelihood;
diff --git a/nonlinear-filters/src/conditional_particle_filter.m b/nonlinear-filters/src/conditional_particle_filter.m
new file mode 100644
index 000000000..3f365f85c
--- /dev/null
+++ b/nonlinear-filters/src/conditional_particle_filter.m
@@ -0,0 +1,125 @@
+function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,DynareOptions)
+%
+% Evaluates the likelihood of a non-linear model with a particle filter
+% - 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.
+%
+% 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 (C) 2009-2010 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 .
+
+% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr
+% stephane DOT adjemian AT univ DASH lemans DOT fr
+
+persistent init_flag mf0 mf1
+persistent number_of_particles
+persistent sample_size number_of_state_variables number_of_observed_variables
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+% Set persistent variables.
+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);
+ init_flag = 1;
+ number_of_particles = DynareOptions.particle.number_of_particles ;
+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)';
+
+% Set seed for randn().
+set_dynare_seed('default');
+
+% Initialization of the likelihood.
+normconst2 = log(2*pi)*number_of_observed_variables*prod(diag(H_lower_triangular_cholesky)) ;
+lik = NaN(sample_size,1);
+LIK = NaN;
+
+ks = 0 ;
+
+StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
+SampleWeights = ones(1,number_of_particles)/number_of_particles ;
+for t=1:sample_size
+ for i=1:number_of_particles
+ [StateParticles(:,i),SampleWeights(i)] = ...
+ conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions,normconst2) ;
+ end
+ SumSampleWeights = sum(SampleWeights) ;
+ lik(t) = log(SumSampleWeights) ;
+ SampleWeights = SampleWeights./SumSampleWeights ;
+ if (DynareOptions.particle.resampling.status.generic && neff(SampleWeights).
+
+[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] = probability(StateMu,StateSqrtP,StateWeights,X);
+ 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/nonlinear-filters/src/gaussian_densities.m b/nonlinear-filters/src/gaussian_densities.m
new file mode 100644
index 000000000..3a770118d
--- /dev/null
+++ b/nonlinear-filters/src/gaussian_densities.m
@@ -0,0 +1,52 @@
+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,DynareOptions)
+%
+% 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 (C) 2009-2010 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 .
+
+% 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,DynareOptions) ;
+eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
+yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
+tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
+Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
+sqr_det = sqrt(det(Pyy)) ;
+foo = (eta_t_i/Pyy).*eta_t_i ;
+likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;
+IncrementalWeights = likelihood.*prior./proposal ;
diff --git a/nonlinear-filters/src/gaussian_filter.m b/nonlinear-filters/src/gaussian_filter.m
new file mode 100644
index 000000000..da3840c74
--- /dev/null
+++ b/nonlinear-filters/src/gaussian_filter.m
@@ -0,0 +1,162 @@
+function [LIK,lik] = gaussian_filter(ReducedForm,Y,start,DynareOptions)
+% 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 Kronrod-Paterson gaussian quadrature with a limited number of nodes.
+% Multidimensional quadrature is obtained by the Smolyak operator (ref: Winschel & Kratzig, 2010).
+% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2008,2009).
+% - a scaled unscented transform cubature (ref: )
+% - 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.
+% => The use of nodes is much faster than Monte-Carlo Gaussian particle and standard particles
+% filters since it treats a lesser number of particles and there is no need
+% of resampling.
+% However, estimations may reveal biaised if the model is truly non-gaussian
+% since predictive and filtered densities are unimodal.
+%
+% 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 (C) 2009-2013 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 .
+
+persistent init_flag mf0 mf1
+persistent nodes2 weights2 weights_c2 number_of_particles
+persistent sample_size number_of_state_variables number_of_observed_variables
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+% Set persistent variables.
+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_particles = DynareOptions.particle.number_of_particles;
+ init_flag = 1;
+end
+
+% compute gaussian quadrature nodes and weights on states and shocks
+if isempty(nodes2)
+ if DynareOptions.particle.distribution_approximation.cubature
+ [nodes2,weights2] = spherical_radial_sigma_points(number_of_state_variables);
+ weights_c2 = weights2;
+ elseif DynareOptions.particle.distribution_approximation.unscented
+ [nodes2,weights2,weights_c2] = unscented_sigma_points(number_of_state_variables,DynareOptions);
+ else
+ if ~DynareOptions.particle.distribution_approximation.montecarlo
+ error('Estimation: This approximation for the proposal is not implemented or unknown!')
+ end
+ end
+end
+
+if DynareOptions.particle.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;
+
+SampleWeights = 1/number_of_particles ;
+ks = 0 ;
+%Estimate = zeros(number_of_state_variables,sample_size) ;
+%V_Estimate = zeros(number_of_state_variables,number_of_state_variables,sample_size) ;
+for t=1:sample_size
+ % build the proposal
+ [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = ...
+ gaussian_filter_bank(ReducedForm,Y(:,t),StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions) ;
+ %Estimate(:,t) = PredictedStateMean ;
+ %V_Estimate(:,:,t) = PredictedStateVarianceSquareRoot ;
+ if DynareOptions.particle.distribution_approximation.cubature || DynareOptions.particle.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,DynareOptions) ;
+ SampleWeights = weights2.*IncrementalWeights ;
+ SumSampleWeights = sum(SampleWeights) ;
+ lik(t) = log(SumSampleWeights) ;
+ SampleWeights = SampleWeights./SumSampleWeights ;
+ else % Monte-Carlo draws
+ 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,DynareOptions) ;
+ SampleWeights = SampleWeights.*IncrementalWeights ;
+ SumSampleWeights = sum(SampleWeights) ;
+ %VarSampleWeights = IncrementalWeights-SumSampleWeights ;
+ %VarSampleWeights = VarSampleWeights*VarSampleWeights'/(number_of_particles-1) ;
+ lik(t) = log(SumSampleWeights) ; %+ .5*VarSampleWeights/(number_of_particles*(SumSampleWeights*SumSampleWeights)) ;
+ SampleWeights = SampleWeights./SumSampleWeights ;
+ Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
+ if (Neff<.5*sample_size && DynareOptions.particle.resampling.status.generic) || DynareOptions.particle.resampling.status.systematic
+ ks = ks + 1 ;
+ StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ;
+ StateVectorMean = mean(StateParticles,2) ;
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/(number_of_particles-1) - StateVectorMean*(StateVectorMean') )';
+ SampleWeights = 1/number_of_particles ;
+ elseif DynareOptions.particle.resampling.status.none
+ StateVectorMean = (sampleWeights*StateParticles)' ;
+ temp = sqrt(SampleWeights').*StateParticles ;
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )';
+ end
+ end
+end
+
+LIK = -sum(lik(start:end));
diff --git a/nonlinear-filters/src/gaussian_filter_bank.m b/nonlinear-filters/src/gaussian_filter_bank.m
new file mode 100644
index 000000000..6094dbcf1
--- /dev/null
+++ b/nonlinear-filters/src/gaussian_filter_bank.m
@@ -0,0 +1,122 @@
+function [PredictedStateMean,PredictedStateVarianceSquareRoot,StateVectorMean,StateVectorVarianceSquareRoot] = gaussian_filter_bank(ReducedForm,obs,StateVectorMean,StateVectorVarianceSquareRoot,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,DynareOptions)
+%
+% 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 (C) 2009-2012 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 .
+
+persistent init_flag2 mf0 mf1
+persistent number_of_state_variables number_of_observed_variables
+persistent number_of_structural_innovations
+
+% Set local state space model (first-order approximation).
+ghx = ReducedForm.ghx;
+ghu = ReducedForm.ghu;
+% Set local state space model (second-order approximation).
+ghxx = ReducedForm.ghxx;
+ghuu = ReducedForm.ghuu;
+ghxu = ReducedForm.ghxu;
+
+if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ...
+ any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ...
+ any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4))
+ ghx
+ ghu
+ ghxx
+ ghuu
+ ghxu
+end
+
+constant = ReducedForm.constant;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+% Set persistent variables.
+if isempty(init_flag2)
+ mf0 = ReducedForm.mf0;
+ mf1 = ReducedForm.mf1;
+ number_of_state_variables = length(mf0);
+ number_of_observed_variables = length(mf1);
+ number_of_structural_innovations = length(ReducedForm.Q);
+ init_flag2 = 1;
+end
+
+if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.proposal_approximation.montecarlo
+ [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables+number_of_structural_innovations) ;
+ weights_c = weights ;
+elseif DynareOptions.particle.proposal_approximation.unscented
+ [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables+number_of_structural_innovations,DynareOptions);
+else
+ error('Estimation: 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);
+tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+PredictedStateMean = tmp(mf0,:)*weights ;
+PredictedObservedMean = tmp(mf1,:)*weights;
+
+if DynareOptions.particle.proposal_approximation.cubature || DynareOptions.particle.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)] ];
+ [mat1,mat] = qr2(big_mat,0);
+ mat = mat';
+ clear('mat1');
+ PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables);
+ CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
+ StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
+ 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
\ No newline at end of file
diff --git a/nonlinear-filters/src/gaussian_mixture_densities.m b/nonlinear-filters/src/gaussian_mixture_densities.m
new file mode 100644
index 000000000..9a5b28e05
--- /dev/null
+++ b/nonlinear-filters/src/gaussian_mixture_densities.m
@@ -0,0 +1,58 @@
+function IncrementalWeights = gaussian_mixture_densities(obs,StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
+ StateMuPost,StateSqrtPPost,StateWeightsPost,...
+ StateParticles,H,normconst,weigths1,weigths2,ReducedForm,DynareOptions)
+
+%
+% 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 (C) 2009-2012 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 .
+
+% Compute the density of particles under the prior distribution
+[ras,ras,prior] = probability(StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateParticles) ;
+prior = prior' ;
+% Compute the density of particles under the proposal distribution
+[ras,ras,proposal] = probability(StateMuPost,StateSqrtPPost,StateWeightsPost,StateParticles) ;
+proposal = proposal' ;
+% Compute the density of the current observation conditionally to each particle
+yt_t_1_i = measurement_equations(StateParticles,ReducedForm,DynareOptions) ;
+eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
+yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
+tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
+Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
+sqr_det = sqrt(det(Pyy)) ;
+foo = (eta_t_i/Pyy).*eta_t_i ;
+likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;
+IncrementalWeights = likelihood.*prior./proposal ;
+
diff --git a/nonlinear-filters/src/gaussian_mixture_filter.m b/nonlinear-filters/src/gaussian_mixture_filter.m
new file mode 100644
index 000000000..5c60a81fb
--- /dev/null
+++ b/nonlinear-filters/src/gaussian_mixture_filter.m
@@ -0,0 +1,225 @@
+function [LIK,lik] = gaussian_mixture_filter(ReducedForm,Y,start,DynareOptions)
+% 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 Smolyak quadrature à la Kronrod & Paterson (Heiss & Winschel 2010, Winschel & Kratzig 2010).
+% - 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 (C) 2009-2013 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 .
+
+
+persistent init_flag mf0 mf1 Gprime Gsecond
+persistent nodes weights weights_c I J G number_of_particles
+persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
+
+% Set default
+if isempty(start)
+ start = 1;
+end
+
+% Set persistent variables.
+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);
+ G = DynareOptions.particle.mixture_state_variables; % number of GM components in state
+ I = DynareOptions.particle.mixture_structural_shocks ; % number of GM components in structural noise
+ J = DynareOptions.particle.mixture_measurement_shocks ; % number of GM components in observation noise
+ Gprime = G*I ;
+ Gsecond = G*I*J ;
+ number_of_particles = DynareOptions.particle.number_of_particles;
+ init_flag = 1;
+end
+
+SampleWeights = ones(Gsecond,1)/Gsecond ;
+
+% compute gaussian quadrature nodes and weights on states and shocks
+if isempty(nodes)
+ if DynareOptions.particle.distribution_approximation.cubature
+ [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables);
+ weights_c = weights;
+ elseif DynareOptions.particle.distribution_approximation.unscented
+ [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables,DynareOptions);
+ else
+ if ~DynareOptions.particle.distribution_approximation.montecarlo
+ error('Estimation: This approximation for the proposal is not implemented or unknown!')
+ end
+ end
+end
+
+if DynareOptions.particle.distribution_approximation.montecarlo
+ set_dynare_seed('default');
+ SampleWeights = 1/number_of_particles ;
+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 all matrices
+StateWeights = ones(1,G)/G ;
+StateMu = ReducedForm.StateVectorMean*ones(1,G) ;
+StateSqrtP = zeros(number_of_state_variables,number_of_state_variables,G) ;
+for g=1:G
+ StateSqrtP(:,:,g) = reduced_rank_cholesky(ReducedForm.StateVectorVariance)' ;
+end
+
+StructuralShocksWeights = ones(1,I)/I ;
+StructuralShocksMu = zeros(number_of_structural_innovations,I) ;
+StructuralShocksSqrtP = zeros(number_of_structural_innovations,number_of_structural_innovations,I) ;
+for i=1:I
+ StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky ;
+end
+
+ObservationShocksWeights = ones(1,J)/J ;
+ObservationShocksMu = zeros(number_of_observed_variables,J) ;
+ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ;
+for j=1:J
+ ObservationShocksSqrtP(:,:,j) = H_lower_triangular_cholesky ;
+end
+
+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) ;
+
+%estimate = zeros(sample_size,number_of_state_variables,3) ;
+const_lik = (2*pi)^(.5*number_of_observed_variables) ;
+
+ks = 0 ;
+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 ;
+ a = g + (j-1)*G ;
+ b = a + (i-1)*Gprime ;
+ [StateMuPrior(:,a),StateSqrtPPrior(:,:,a),StateWeightsPrior(1,a),...
+ StateMuPost(:,b),StateSqrtPPost(:,:,b),StateWeightsPost(1,b)] =...
+ gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(1,g),...
+ StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(1,i),...
+ ObservationShocksMu(:,j),ObservationShocksSqrtP(:,:,j),ObservationShocksWeights(1,j),...
+ H,H_lower_triangular_cholesky,const_lik,DynareOptions) ;
+ end
+ end
+ end
+
+ % Normalize weights
+ StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior,2) ;
+ StateWeightsPost = StateWeightsPost/sum(StateWeightsPost,2) ;
+
+ if DynareOptions.particle.distribution_approximation.cubature || DynareOptions.particle.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,const_lik,weights,weights_c,ReducedForm,DynareOptions) ;
+ SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ;
+ end
+ SumSampleWeights = sum(SampleWeights) ;
+ lik(t) = log(SumSampleWeights) ;
+ SampleWeights = SampleWeights./SumSampleWeights ;
+ [ras,SortedRandomIndx] = sort(rand(1,Gsecond));
+ SortedRandomIndx = SortedRandomIndx(1:G);
+ indx = index_resample(0,SampleWeights,DynareOptions) ;
+ 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) ;
+ % Compute prior, proposal and likelihood of particles
+ IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
+ StateMuPost,StateSqrtPPost,StateWeightsPost,...
+ StateParticles,H,const_lik,1/number_of_particles,...
+ 1/number_of_particles,ReducedForm,DynareOptions) ;
+ % calculate importance weights of particles
+ SampleWeights = SampleWeights.*IncrementalWeights ;
+ SumSampleWeights = sum(SampleWeights,1) ;
+ SampleWeights = SampleWeights./SumSampleWeights ;
+ lik(t) = log(SumSampleWeights) ;
+ % First possible state point estimates
+ %estimate(t,:,1) = SampleWeights*StateParticles' ;
+ % Resampling if needed of required
+ Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
+ if (DynareOptions.particle.resampling.status.generic && Neff<.5*sample_size) || DynareOptions.particle.resampling.status.systematic
+ ks = ks + 1 ;
+ StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ;
+ StateVectorMean = mean(StateParticles,2) ;
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/number_of_particles - StateVectorMean*(StateVectorMean') )';
+ SampleWeights = 1/number_of_particles ;
+ elseif DynareOptions.particle.resampling.status.none
+ StateVectorMean = StateParticles*sampleWeights ;
+ temp = sqrt(SampleWeights').*StateParticles ;
+ StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';
+ end
+ % Use the information from particles to update the gaussian mixture on state variables
+ [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(StateParticles,StateMu,StateSqrtP,StateWeights,0.001,10,1) ;
+ %estimate(t,:,3) = StateWeights*StateMu' ;
+ end
+end
+
+LIK = -sum(lik(start:end)) ;
\ No newline at end of file
diff --git a/nonlinear-filters/src/gaussian_mixture_filter_bank.m b/nonlinear-filters/src/gaussian_mixture_filter_bank.m
new file mode 100644
index 000000000..3a3caa776
--- /dev/null
+++ b/nonlinear-filters/src/gaussian_mixture_filter_bank.m
@@ -0,0 +1,138 @@
+function [StateMuPrior,StateSqrtPPrior,StateWeightsPrior,StateMuPost,StateSqrtPPost,StateWeightsPost] =...
+ gaussian_mixture_filter_bank(ReducedForm,obs,StateMu,StateSqrtP,StateWeights,...
+ StructuralShocksMu,StructuralShocksSqrtP,StructuralShocksWeights,...
+ ObservationShocksMu,ObservationShocksSqrtP,ObservationShocksWeights,...
+ H,H_lower_triangular_cholesky,normfactO,DynareOptions)
+%
+% 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 (C) 2009-2012 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 .
+
+
+
+persistent init_flag2 mf0 mf1 %nodes3 weights3 weights_c3
+persistent number_of_state_variables number_of_observed_variables
+persistent number_of_structural_innovations
+
+% Set local state space model (first-order approximation).
+ghx = ReducedForm.ghx;
+ghu = ReducedForm.ghu;
+% Set local state space model (second-order approximation).
+ghxx = ReducedForm.ghxx;
+ghuu = ReducedForm.ghuu;
+ghxu = ReducedForm.ghxu;
+
+if any(any(isnan(ghx))) || any(any(isnan(ghu))) || any(any(isnan(ghxx))) || any(any(isnan(ghuu))) || any(any(isnan(ghxu))) || ...
+ any(any(isinf(ghx))) || any(any(isinf(ghu))) || any(any(isinf(ghxx))) || any(any(isinf(ghuu))) || any(any(isinf(ghxu))) ...
+ any(any(abs(ghx)>1e4)) || any(any(abs(ghu)>1e4)) || any(any(abs(ghxx)>1e4)) || any(any(abs(ghuu)>1e4)) || any(any(abs(ghxu)>1e4))
+ ghx
+ ghu
+ ghxx
+ ghuu
+ ghxu
+end
+
+constant = ReducedForm.constant;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+% Set persistent variables.
+if isempty(init_flag2)
+ mf0 = ReducedForm.mf0;
+ mf1 = ReducedForm.mf1;
+ number_of_state_variables = length(mf0);
+ number_of_observed_variables = length(mf1);
+ number_of_structural_innovations = length(ReducedForm.Q);
+ init_flag2 = 1;
+end
+
+numb = number_of_state_variables+number_of_structural_innovations ;
+
+if DynareOptions.particle.proposal_approximation.cubature
+ [nodes3,weights3] = spherical_radial_sigma_points(numb);
+ weights_c3 = weights3;
+elseif DynareOptions.particle.proposal_approximation.unscented
+ [nodes3,weights3,weights_c3] = unscented_sigma_points(numb,DynareOptions);
+else
+ error('Estimation: This approximation for the proposal is not implemented or unknown!')
+end
+
+epsilon = bsxfun(@plus,StructuralShocksSqrtP*nodes3(:,number_of_state_variables+1:number_of_state_variables+number_of_structural_innovations)',StructuralShocksMu) ;
+StateVectors = bsxfun(@plus,StateSqrtP*nodes3(:,1:number_of_state_variables)',StateMu);
+yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
+tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+PredictedStateMean = tmp(mf0,:)*weights3;
+PredictedObservedMean = tmp(mf1,:)*weights3;
+
+if DynareOptions.particle.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)] ];
+ [mat1,mat] = qr2(big_mat,0);
+ mat = mat';
+ clear('mat1');
+ PredictedObservedVarianceSquareRoot = mat(1:number_of_observed_variables,1:number_of_observed_variables);
+ CovarianceObservedStateSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),1:number_of_observed_variables);
+ StateVectorVarianceSquareRoot = mat(number_of_observed_variables+(1:number_of_state_variables),number_of_observed_variables+(1:number_of_state_variables));
+ 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/nonlinear-filters/src/importance_sampling.m b/nonlinear-filters/src/importance_sampling.m
new file mode 100644
index 000000000..403b6beb3
--- /dev/null
+++ b/nonlinear-filters/src/importance_sampling.m
@@ -0,0 +1,29 @@
+function State_Particles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost,numP)
+
+% Copyright (C) 2013 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 .
+
+[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/nonlinear-filters/src/index_resample.m b/nonlinear-filters/src/index_resample.m
new file mode 100644
index 000000000..e91645605
--- /dev/null
+++ b/nonlinear-filters/src/index_resample.m
@@ -0,0 +1,72 @@
+function resampling_index = index_resample(particles,weights,DynareOptions)
+% Resamples particles.
+
+%@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 (C) 2013 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 .
+
+defaultmethod = 1; % For residual based method set this variable equal to 0.
+
+if defaultmethod
+ if DynareOptions.particle.resampling.method.kitagawa
+ resampling_index = traditional_resampling(particles,weights,rand);
+ elseif DynareOptions.particle.resampling.method.stratified
+ resampling_index = traditional_resampling(particles,weights,rand(size(weights)));
+ else
+ error('Unknow sampling method!')
+ end
+else
+ if DynareOptions.particle.resampling.method.kitagawa
+ resampled_particles = residual_resampling(particles,weights,rand);
+ elseif DynareOptions.particle.resampling.method.stratified
+ resampled_particles = residual_resampling(particles,weights,rand(size(weights)));
+ else
+ error('Unknown sampling method!')
+ end
+end
\ No newline at end of file
diff --git a/nonlinear-filters/src/local_state_space_iteration/local_state_space_iteration_2.m b/nonlinear-filters/src/local_state_space_iteration/local_state_space_iteration_2.m
new file mode 100644
index 000000000..0c6a52d07
--- /dev/null
+++ b/nonlinear-filters/src/local_state_space_iteration/local_state_space_iteration_2.m
@@ -0,0 +1,275 @@
+function [y,y_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,a,b,c)%yhat_,ss)
+
+%@info:
+%! @deftypefn {Function File} {@var{y}, @var{y_} =} local_state_equation_2 (@var{yhat},@var{epsilon}, @var{ghx}, @var{ghu}, @var{constant}, @var{ghxx}, @var{ghuu}, @var{ghxu}, @var{yhat_}, @var{ss})
+%! @anchor{particle/local_state_space_iteration_2}
+%! @sp 1
+%! Given the states (y) and structural innovations (epsilon), this routine computes the level of selected endogenous variables when the
+%! model is approximated by an order two taylor expansion around the deterministic steady state. Depending on the number of input/output
+%! argument the pruning algorithm advocated by C. Sims is used or not (this version should not be used if the selected endogenous variables
+%! are not the states of the model).
+%!
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item yhat
+%! n*1 vector of doubles, initial condition, where n is the number of state variables.
+%! @item epsilon
+%! q*1 vector of doubles, structural innovations.
+%! @item ghx
+%! m*n matrix of doubles, restricted dr.ghx where we only consider the lines corresponding to a subset of endogenous variables.
+%! @item ghu
+%! m*q matrix of doubles, restricted dr.ghu where we only consider the lines corresponding to a subset of endogenous variables.
+%! @item constant
+%! m*1 vector of doubles, deterministic steady state plus second order correction for a subset of endogenous variables.
+%! @item ghxx
+%! m*n² matrix of doubles, restricted dr.ghxx where we only consider the lines corresponding to a subset of endogenous variables.
+%! @item ghuu
+%! m*q² matrix of doubles, restricted dr.ghuu where we only consider the lines corresponding to a subset of endogenous variables.
+%! @item ghxu
+%! m*(nq) matrix of doubles, subset of dr.ghxu where we only consider the lines corresponding to a subset of endogenous variables.
+%! @item yhat_
+%! n*1 vector of doubles, spurious states for the pruning version.
+%! @item ss
+%! n*1 vector of doubles, steady state for the states.
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item y
+%! m*1 vector of doubles, selected endogenous variables.
+%! @item y_
+%! m*1 vector of doubles, update of the latent variables needed for the pruning version (first order update).
+%! @end table
+%! @sp 2
+%! @strong{Remarks}
+%! @sp 1
+%! [1] If the function has 10 input arguments then it must have 2 output arguments (pruning version).
+%! @sp 1
+%! [2] If the function has 08 input arguments then it must have 1 output argument.
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 2
+%! @strong{This function calls:}
+%!
+%!
+%! @end deftypefn
+%@eod:
+
+% Copyright (C) 2011-2012 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 .
+
+% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
+% frederic DOT karame AT univ DASH evry DOT fr
+
+if nargin==9
+ pruning = 0; numthreads = a;
+ if nargout>1
+ error('local_state_space_iteration_2:: Numbers of input and output argument are inconsistent!')
+ end
+elseif nargin==11
+ pruning = 1; numthreads = c; yhat_ = a; ss = b;
+ if nargout~=2
+ error('local_state_space_iteration_2:: Numbers of input and output argument are inconsistent!')
+ end
+else
+ error('local_state_space_iteration_2:: Wrong number of input arguments!')
+end
+
+number_of_threads = numthreads;
+
+switch pruning
+ case 0
+ for i =1:size(yhat,2)
+ y(:,i) = constant + ghx*yhat(:,i) + ghu*epsilon(:,i) ...
+ + A_times_B_kronecker_C(.5*ghxx,yhat(:,i),number_of_threads) ...
+ + A_times_B_kronecker_C(.5*ghuu,epsilon(:,i),number_of_threads) ...
+ + A_times_B_kronecker_C(ghxu,yhat(:,i),epsilon(:,i),number_of_threads);
+ end
+ case 1
+ for i =1:size(yhat,2)
+ y(:,i) = constant + ghx*yhat(:,i) + ghu*epsilon(:,i) ...
+ + A_times_B_kronecker_C(.5*ghxx,yhat_(:,i),number_of_threads) ...
+ + A_times_B_kronecker_C(.5*ghuu,epsilon(:,i),number_of_threads) ...
+ + A_times_B_kronecker_C(ghxu,yhat_(:,i),epsilon(:,i),number_of_threads);
+ end
+ y_ = ghx*yhat_+ghu*epsilon;
+ y_ = bsxfun(@plus,y_,ss);
+end
+
+%@test:1
+%$ n = 2;
+%$ q = 3;
+%$
+%$ yhat = zeros(n,1);
+%$ epsilon = zeros(q,1);
+%$ ghx = rand(n,n);
+%$ ghu = rand(n,q);
+%$ constant = ones(n,1);
+%$ ghxx = rand(n,n*n);
+%$ ghuu = rand(n,q*q);
+%$ ghxu = rand(n,n*q);
+%$ yhat_ = zeros(n,1);
+%$ ss = ones(n,1);
+%$
+%$ % Call the tested routine.
+%$ for i=1:10
+%$ y1 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1);
+%$ [y2,y2_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1);
+%$ end
+%$
+%$ % Check the results.
+%$ t(1) = dassert(y1,ones(n,1));
+%$ t(2) = dassert(y2,ones(n,1));
+%$ t(3) = dassert(y2_,ones(n,1));
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ old_path = pwd;
+%$ cd([fileparts(which('dynare')) '/../tests/']);
+%$ dynare('dsge_base2');
+%$ load dsge_base2;
+%$ cd(old_path);
+%$ dr = oo_.dr;
+%$ clear('oo_','options_','M_');
+%$ delete([fileparts(which('dynare')) '/../tests/dsge_base2.mat']);
+%$ istates = dr.nstatic+(1:dr.npred);
+%$ n = dr.npred;
+%$ q = size(dr.ghu,2);
+%$ yhat = zeros(n,1);
+%$ epsilon = zeros(q,1);
+%$ ghx = dr.ghx(istates,:);
+%$ ghu = dr.ghu(istates,:);
+%$ constant = dr.ys(istates,:)+dr.ghs2(istates,:);
+%$ ghxx = dr.ghxx(istates,:);
+%$ ghuu = dr.ghuu(istates,:);
+%$ ghxu = dr.ghxu(istates,:);
+%$ yhat_ = zeros(n,1);
+%$ ss = dr.ys(istates,:);
+%$
+%$ t = ones(2,1);
+%$
+%$ % Call the tested routine.
+%$ try
+%$ y1 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1);
+%$ catch
+%$ t(1) = 0;
+%$ end
+%$ try
+%$ [y2,y2_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1);
+%$ catch
+%$ t(2) = 0;
+%$ end
+%$
+%$ % Check the results.
+%$ T = all(t);
+%@eof:2
+
+%@test:3
+%$ Bohrbug = 1; % A bug that manifests reliably under a possibly unknown but well-defined set of conditions.
+%$ if ~Bohrbug
+%$ n = 2;
+%$ q = 3;
+%$
+%$ yhat = .01*randn(n,1);
+%$ epsilon = .001*randn(q,1);
+%$ ghx = rand(n,n);
+%$ ghu = rand(n,q);
+%$ constant = ones(n,1);
+%$ ghxx = rand(n,n*n);
+%$ ghuu = rand(n,q*q);
+%$ ghxu = rand(n,n*q);
+%$ yhat_ = zeros(n,1);
+%$ ss = ones(n,1);
+%$
+%$ % Call the tested routine (mex version).
+%$ y1a = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1);
+%$ [y2a,y2a_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1);
+%$
+%$ % Call the tested routine (matlab version)
+%$ path_to_mex = fileparts(which(['qmc_sequence.' mexext]));
+%$ where_am_i_coming_from = pwd;
+%$ cd(path_to_mex);
+%$ tar('local_state_space_iteration_2.tar',['local_state_space_iteration_2.' mexext]);
+%$ cd(where_am_i_coming_from);
+%$ dynare_config([],0);
+%$ y1b = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1);
+%$ [y2b,y2b_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,ss,1);
+%$ cd(path_to_mex);
+%$ untar('local_state_space_iteration_2.tar');
+%$ delete('local_state_space_iteration_2.tar');
+%$ cd(where_am_i_coming_from);
+%$ dynare_config([],0);
+%$ % Check the results.
+%$ t(1) = dassert(y1a,y1b);
+%$ t(2) = dassert(y2a,y2b);
+%$ t(3) = dassert(y2a_,y2b_);
+%$ T = all(t);
+%$ else
+%$ t(1) = 1;
+%$ T = all(t);
+%$ end
+%@eof:3
+
+
+%@test:4
+%$ % TIMING TEST (parallelization with openmp)
+%$ old_path = pwd;
+%$ cd([fileparts(which('dynare')) '/../tests/']);
+%$ dynare('dsge_base2');
+%$ load dsge_base2;
+%$ cd(old_path);
+%$ dr = oo_.dr;
+%$ clear('oo_','options_','M_');
+%$ delete([fileparts(which('dynare')) '/../tests/dsge_base2.mat']);
+%$ istates = dr.nstatic+(1:dr.npred);
+%$ n = dr.npred;
+%$ q = size(dr.ghu,2);
+%$ yhat = zeros(n,10000000);
+%$ epsilon = zeros(q,10000000);
+%$ ghx = dr.ghx(istates,:);
+%$ ghu = dr.ghu(istates,:);
+%$ constant = dr.ys(istates,:)+dr.ghs2(istates,:);
+%$ ghxx = dr.ghxx(istates,:);
+%$ ghuu = dr.ghuu(istates,:);
+%$ ghxu = dr.ghxu(istates,:);
+%$ yhat_ = zeros(n,10000000);
+%$ ss = dr.ys(istates,:);
+%$
+%$ t = NaN(4,1);
+%$ tic, for i=1:10, y1 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,1); end
+%$ t1 = toc;
+%$ tic, for i=1:10, y2 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,2); end
+%$ t2 = toc;
+%$ t(1) = dassert(y1,y2,1e-15); clear('y1');
+%$ tic, for i=1:10, y3 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,3); end
+%$ t3 = toc;
+%$ t(2) = dassert(y2,y3,1e-15); clear('y2');
+%$ tic, for i=1:10, y4 = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,4); end
+%$ t4 = toc;
+%$ t(3) = dassert(y4,y3,1e-15); clear('y3','y4');
+%$ t(4) = (t1>t2) && (t2>t3) && (t3>t4);
+%$ if ~t(4)
+%$ disp('Timmings:')
+%$ [t1, t2, t3, t4]
+%$ end
+%$ % Check the results.
+%$ T = all(t);
+%@eof:4
diff --git a/nonlinear-filters/src/measurement_equations.m b/nonlinear-filters/src/measurement_equations.m
new file mode 100644
index 000000000..c00d8ebc4
--- /dev/null
+++ b/nonlinear-filters/src/measurement_equations.m
@@ -0,0 +1,33 @@
+function measure = measurement_equations(StateVectors,ReducedForm,DynareOptions)
+
+% Copyright (C) 2013 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 .
+
+mf1 = ReducedForm.mf1;
+ghx = ReducedForm.ghx(mf1,:);
+ghu = ReducedForm.ghu(mf1,:);
+ghxx = ReducedForm.ghxx(mf1,:);
+ghuu = ReducedForm.ghuu(mf1,:);
+ghxu = ReducedForm.ghxu(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) ;
+measure = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,size(yhat,2)),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+
+
+
diff --git a/nonlinear-filters/src/multivariate_smooth_resampling.m b/nonlinear-filters/src/multivariate_smooth_resampling.m
new file mode 100644
index 000000000..09ca287b9
--- /dev/null
+++ b/nonlinear-filters/src/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 (C) 2012-2013 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 .
+
+% 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/nonlinear-filters/src/mykmeans.m b/nonlinear-filters/src/mykmeans.m
new file mode 100644
index 000000000..56af5390b
--- /dev/null
+++ b/nonlinear-filters/src/mykmeans.m
@@ -0,0 +1,53 @@
+function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod)
+
+% Copyright (C) 2013 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 .
+
+[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/nonlinear-filters/src/neff.m b/nonlinear-filters/src/neff.m
new file mode 100644
index 000000000..5251b9f17
--- /dev/null
+++ b/nonlinear-filters/src/neff.m
@@ -0,0 +1,21 @@
+function n = neff(w)
+% Evaluates the criterion for resampling
+
+% Copyright (C) 2013 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 .
+
+n = dot(w,w);
diff --git a/nonlinear-filters/src/online_auxiliary_filter.m b/nonlinear-filters/src/online_auxiliary_filter.m
new file mode 100644
index 000000000..4a8243a2d
--- /dev/null
+++ b/nonlinear-filters/src/online_auxiliary_filter.m
@@ -0,0 +1,373 @@
+function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
+
+% Carvalho & Lopes particle filter = auxiliary particle filter including Liu & West filter on parameters.
+%
+% INPUTS
+% ReducedForm [structure] Matlab's structure describing the reduced form model.
+% ReducedForm.measurement.H [double] (pp x pp) variance matrix of measurement errors.
+% ReducedForm.state.Q [double] (qq x qq) variance matrix of state errors.
+% ReducedForm.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'.
+% mf [integer] pp*1 vector of indices.
+% number_of_particles [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 (C) 2013 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 .
+persistent Y init_flag mf0 mf1 bounds number_of_particles number_of_parameters liu_west_delta liu_west_chol_sigma_bar
+persistent start_param sample_size number_of_observed_variables number_of_structural_innovations
+
+% Set seed for randn().
+set_dynare_seed('default') ;
+pruning = DynareOptions.particle.pruning;
+second_resample = 1 ;
+variance_update = 1 ;
+
+% initialization of state particles
+[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
+ solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
+
+% Set persistent variables.
+if isempty(init_flag)
+ 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,2);
+ number_of_observed_variables = length(mf1);
+ number_of_structural_innovations = length(ReducedForm.Q);
+ liu_west_delta = DynareOptions.particle.liu_west_delta ;
+ %liu_west_chol_sigma_bar = DynareOptions.particle.liu_west_chol_sigma_bar*eye(number_of_parameters) ;
+ start_param = xparam1 ;
+ %liu_west_chol_sigma_bar = sqrt(bsxfun(@times,eye(number_of_parameters),BayesInfo.p2)) ;
+ %start_param = BayesInfo.p1 ;
+ bounds = [BayesInfo.lb BayesInfo.ub] ;
+ init_flag = 1;
+end
+
+% Get initial conditions for the state particles
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
+state_variance_rank = size(StateVectorVarianceSquareRoot,2);
+StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
+if pruning
+ StateVectors_ = StateVectors;
+end
+
+% parameters for the Liu & West filter
+h_square = (3*liu_west_delta-1)/(2*liu_west_delta) ;
+h_square = 1-h_square*h_square ;
+small_a = sqrt(1-h_square) ;
+
+% Initialization of parameter particles
+xparam = zeros(number_of_parameters,number_of_particles) ;
+stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/100 ;
+stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/50 ;
+i = 1 ;
+while i<=number_of_particles
+ %candidate = start_param + .001*liu_west_chol_sigma_bar*randn(number_of_parameters,1) ;
+ candidate = start_param + bsxfun(@times,stderr,randn(number_of_parameters,1)) ;
+ if all(candidate(:) >= bounds(:,1)) && all(candidate(:) <= bounds(:,2))
+ xparam(:,i) = candidate(:) ;
+ i = i+1 ;
+ end
+end
+
+%xparam = bsxfun(@plus,bounds(:,1),bsxfun(@times,(bounds(:,2)-bounds(:,1)),rand(number_of_parameters,number_of_particles))) ;
+
+% 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) ;
+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
+ disp(t)
+ % 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==1
+ chol_sigma_bar = chol(h_square*sigma_bar)' ;
+ end
+ % Prediction (without shocks)
+ ObservedVariables = zeros(number_of_observed_variables,number_of_particles) ;
+ for i=1:number_of_particles
+ % model resolution
+ [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
+ solve_model_for_online_filter(t,xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
+ steadystate = ReducedForm.steadystate;
+ state_variables_steady_state = ReducedForm.state_variables_steady_state;
+ % Set local state space model (second-order approximation).
+ constant = ReducedForm.constant;
+ ghx = ReducedForm.ghx;
+ ghu = ReducedForm.ghu;
+ ghxx = ReducedForm.ghxx;
+ ghuu = ReducedForm.ghuu;
+ ghxu = ReducedForm.ghxu;
+ % particle likelihood contribution
+ yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
+ if pruning
+ yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
+ else
+ tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+ end
+ ObservedVariables(:,i) = tmp(mf1,:) ;
+ end
+ PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
+ PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables);
+ dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
+ PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
+ wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
+ % unormalized weights and observation likelihood contribution
+ tau_tilde = weights.*wtilde ;
+ sum_tau_tilde = sum(tau_tilde) ;
+ % particles selection
+ tau_tilde = tau_tilde/sum_tau_tilde ;
+ indx = index_resample(0,tau_tilde',DynareOptions);
+ StateVectors = StateVectors(:,indx) ;
+ if pruning
+ StateVectors_ = StateVectors_(:,indx) ;
+ end
+ xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ;
+ xparam = xparam(:,indx) ;
+ wtilde = wtilde(indx) ;
+ % draw in the new distributions
+ i = 1 ;
+ while i<=number_of_particles
+ candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters,1) ;
+ if all(candidate >= bounds(:,1)) && all(candidate <= bounds(:,2))
+ xparam(:,i) = candidate ;
+ % model resolution for new parameters particles
+ [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
+ solve_model_for_online_filter(t,xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
+ steadystate = ReducedForm.steadystate;
+ state_variables_steady_state = ReducedForm.state_variables_steady_state;
+ % Set local state space model (second order approximation).
+ constant = ReducedForm.constant;
+ ghx = ReducedForm.ghx;
+ ghu = ReducedForm.ghu;
+ ghxx = ReducedForm.ghxx;
+ ghuu = ReducedForm.ghuu;
+ ghxu = ReducedForm.ghxu;
+ % Get covariance matrices and structural shocks
+ epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations,1) ;
+ % compute particles likelihood contribution
+ yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
+ if pruning
+ yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
+ StateVectors_(:,i) = tmp_(mf0,:) ;
+ else
+ tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+ end
+ StateVectors(:,i) = tmp(mf0,:) ;
+ ObservedVariables(:,i) = tmp(mf1,:) ;
+ i = i+1 ;
+ end
+ end
+ PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
+ PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables);
+ dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
+ PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
+ lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)));
+ % importance ratio
+ wtilde = lnw./wtilde ;
+ % normalization
+ weights = wtilde/sum(wtilde);
+ if (variance_update==1) && (neff(weights)0.025 && pass1==1
+ lb95_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ;
+ pass1 = 2 ;
+ end
+ if cumulated_weights(j)>0.5 && pass2==1
+ median_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ;
+ pass2 = 2 ;
+ end
+ if cumulated_weights(j)>0.975 && pass3==1
+ ub95_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ;
+ pass3 = 2 ;
+ end
+ end
+ end
+ end
+ disp([lb95_xparam(:,t) mean_xparam(:,t) ub95_xparam(:,t)])
+end
+distrib_param = xparam ;
+xparam = mean_xparam(:,sample_size) ;
+std_param = std_xparam(:,sample_size) ;
+lb_95 = lb95_xparam(:,sample_size) ;
+ub_95 = ub95_xparam(:,sample_size) ;
+median_param = median_xparam(:,sample_size) ;
+
+%% Plot parameters trajectory
+TeX = DynareOptions.TeX;
+
+[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters);
+
+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
+
+z = 1:1:sample_size ;
+
+for plt = 1:nbplt,
+ if TeX
+ NAMES = [];
+ TeXNAMES = [];
+ end
+ hh = dyn_figure(DynareOptions,'Name','Parameters Trajectories');
+ for k=1:min(nstar,length(xparam)-(plt-1)*nstar)
+ subplot(nr,nc,k)
+ kk = (plt-1)*nstar+k;
+ [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
+ if TeX
+ if isempty(NAMES)
+ NAMES = name;
+ TeXNAMES = texname;
+ else
+ NAMES = char(NAMES,name);
+ TeXNAMES = char(TeXNAMES,texname);
+ end
+ end
+ y = [mean_xparam(kk,:)' median_xparam(kk,:)' lb95_xparam(kk,:)' ub95_xparam(kk,:)' xparam(kk)*ones(sample_size,1)] ;
+ plot(z,y);
+ hold on
+ title(name,'interpreter','none')
+ hold off
+ axis tight
+ drawnow
+ end
+ dyn_saveas(hh,[ Model.fname '_param_traj' int2str(plt) ],DynareOptions);
+ if TeX
+ % TeX eps loader file
+ fprintf(fidTeX,'\\begin{figure}[H]\n');
+ for jj = 1:min(nstar,length(x)-(plt-1)*nstar)
+ fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
+ end
+ 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,
+ if TeX
+ NAMES = [];
+ TeXNAMES = [];
+ end
+ hh = dyn_figure(DynareOptions,'Name','Parameters Densities');
+ for k=1:min(nstar,length(xparam)-(plt-1)*nstar)
+ subplot(nr,nc,k)
+ kk = (plt-1)*nstar+k;
+ [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
+ if TeX
+ if isempty(NAMES)
+ NAMES = name;
+ TeXNAMES = texname;
+ else
+ NAMES = char(NAMES,name);
+ TeXNAMES = char(TeXNAMES,texname);
+ end
+ end
+ optimal_bandwidth = mh_optimal_bandwidth(distrib_param(kk,:)',number_of_particles,bandwidth,kernel_function);
+ [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(kk,:)',number_of_grid_points,...
+ number_of_particles,optimal_bandwidth,kernel_function);
+ plot(density(:,1),density(:,2));
+ hold on
+ title(name,'interpreter','none')
+ hold off
+ axis tight
+ drawnow
+ end
+ dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions);
+ if TeX
+ % TeX eps loader file
+ fprintf(fidTeX,'\\begin{figure}[H]\n');
+ for jj = 1:min(nstar,length(x)-(plt-1)*nstar)
+ fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
+ end
+ fprintf(fidTeX,'\\centering \n');
+ fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParametersDensities%s}\n',Model.fname,int2str(plt));
+ fprintf(fidTeX,'\\caption{ParametersDensities.}');
+ fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
+ fprintf(fidTeX,'\\end{figure}\n');
+ fprintf(fidTeX,' \n');
+ end
+end
+
\ No newline at end of file
diff --git a/nonlinear-filters/src/probability.m b/nonlinear-filters/src/probability.m
new file mode 100644
index 000000000..5d76e49ee
--- /dev/null
+++ b/nonlinear-filters/src/probability.m
@@ -0,0 +1,39 @@
+function [prior,likelihood,C,posterior] = probability(mu,sqrtP,prior,X)
+
+% Copyright (C) 2013 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 .
+
+[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/nonlinear-filters/src/probability2.m b/nonlinear-filters/src/probability2.m
new file mode 100644
index 000000000..57ab65f8f
--- /dev/null
+++ b/nonlinear-filters/src/probability2.m
@@ -0,0 +1,38 @@
+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 (C) 2009-2012 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 .
+
+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 ;
+
diff --git a/nonlinear-filters/src/resample.m b/nonlinear-filters/src/resample.m
new file mode 100644
index 000000000..a72894428
--- /dev/null
+++ b/nonlinear-filters/src/resample.m
@@ -0,0 +1,74 @@
+function resampled_particles = resample(particles,weights,DynareOptions)
+% Resamples particles.
+
+%@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 (C) 2011-2013 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 .
+
+defaultmethod = 1; % For residual based method set this variable equal to 0.
+
+if defaultmethod
+ if DynareOptions.particle.resampling.method.kitagawa
+ resampled_particles = traditional_resampling(particles,weights,rand);
+ elseif DynareOptions.particle.resampling.method.stratified
+ resampled_particles = traditional_resampling(particles,weights,rand(size(weights)));
+ elseif DynareOptions.particle.resampling.method.smooth
+ resampled_particles = multivariate_smooth_resampling(particles,weights);
+ else
+ error('Unknow sampling method!')
+ end
+else
+ if DynareOptions.particle.resampling.method.kitagawa
+ resampled_particles = residual_resampling(particles,weights,rand);
+ elseif DynareOptions.particle.resampling.method.stratified
+ resampled_particles = residual_resampling(particles,weights,rand(size(weights)));
+ else
+ error('Unknown sampling method!')
+ end
+end
diff --git a/nonlinear-filters/src/residual_resampling.m b/nonlinear-filters/src/residual_resampling.m
new file mode 100644
index 000000000..c2f09e3b8
--- /dev/null
+++ b/nonlinear-filters/src/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 (C) 2011-2013 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 .
+
+% 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/nonlinear-filters/src/sequential_importance_particle_filter.m b/nonlinear-filters/src/sequential_importance_particle_filter.m
new file mode 100644
index 000000000..eadaeba41
--- /dev/null
+++ b/nonlinear-filters/src/sequential_importance_particle_filter.m
@@ -0,0 +1,178 @@
+function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,DynareOptions)
+% Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
+% Standard Sequential Monte Carlo approach with
+% - the usual proposal (the state transition distribution)
+% - options on resampling: none, adaptive or systematic
+%@info:
+%! @deftypefn {Function File} {@var{y}, @var{y_} =} sequential_importance_particle_filter (@var{ReducedForm},@var{Y}, @var{start}, @var{DynareOptions})
+%! @anchor{particle/sequential_importance_particle_filter}
+%! @sp 1
+%! Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
+%!
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item ReducedForm
+%! Structure describing the state space model (built in @ref{non_linear_dsge_likelihood}).
+%! @item Y
+%! p*smpl matrix of doubles (p is the number of observed variables), the (detrended) data.
+%! @item start
+%! Integer scalar, likelihood evaluation starts at observation 'start'.
+%! @item DynareOptions
+%! Structure specifying Dynare's options.
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item LIK
+%! double scalar, value of (minus) the logged likelihood.
+%! @item lik
+%! smpl*1 vector of doubles, density of the observations at each period.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @ref{non_linear_dsge_likelihood}
+%! @sp 2
+%! @strong{This function calls:}
+%!
+%! @end deftypefn
+%@eod:
+
+% Copyright (C) 2011-2013 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 .
+
+% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr
+% stephane DOT adjemian AT univ DASH lemans DOT fr
+
+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 = DynareOptions.particle.pruning;
+
+% Get steady state and mean.
+steadystate = ReducedForm.steadystate;
+constant = ReducedForm.constant;
+state_variables_steady_state = ReducedForm.state_variables_steady_state;
+
+% 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 = DynareOptions.particle.number_of_particles;
+ init_flag = 1;
+end
+
+% 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;
+
+% 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;
+lik = NaN(sample_size,1);
+
+% Get initial condition for the state vector.
+StateVectorMean = ReducedForm.StateVectorMean;
+StateVectorVarianceSquareRoot = 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
+ StateVectors_ = StateVectors;
+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);
+ [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
+ else
+ tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
+ 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;
+ if rcond(PredictedObservedVariance) > 1e-16
+ lnw = -.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\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 (DynareOptions.particle.resampling.status.generic && neff(weights).
+
+% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
+% frederic DOT karame AT univ DASH lemans DOT fr
+
+%global objective_function_penalty_base
+% Declaration of the penalty as a persistent variable.
+persistent init_flag
+persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1
+persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations
+
+% Initialization of the returned arguments.
+fval = [];
+ys = [];
+trend_coeff = [];
+exit_flag = 1;
+
+% Set the number of observed variables
+nvobs = DynareDataset.info.nvobs;
+
+%------------------------------------------------------------------------------
+% 1. Get the structural parameters & define penalties
+%------------------------------------------------------------------------------
+
+% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
+%if (DynareOptions.mode_compute~=1) && any(xparam1BayesInfo.ub)
+% k = find(xparam1(:)>BayesInfo.ub);
+% fval = objective_function_penalty_base+sum((xparam1(k)-BayesInfo.ub(k)).^2);
+% exit_flag = 0;
+% 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(nvobs);
+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)
+% [CholQ,testQ] = chol(Q);
+% if testQ
+ % The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
+% a = diag(eig(Q));
+% k = find(a < 0);
+% if k > 0
+% fval = objective_function_penalty_base+sum(-a(k));
+% exit_flag = 0;
+% info = 43;
+% return
+% end
+% 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)
+% [CholH,testH] = chol(H);
+% if testH
+ % The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
+% a = diag(eig(H));
+% k = find(a < 0);
+% if k > 0
+% fval = objective_function_penalty_base+sum(-a(k));
+% exit_flag = 0;
+% info = 44;
+% return
+% end
+% 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
+%------------------------------------------------------------------------------
+
+% Linearize the model around the deterministic sdteadystate and extract the matrices of the state equation (T and R).
+[T,R,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+
+%if info(1) == 1 || info(1) == 2 || info(1) == 5
+% fval = objective_function_penalty_base+1;
+% exit_flag = 0;
+% return
+%elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21
+% fval = objective_function_penalty_base+info(2);
+% exit_flag = 0;
+% return
+%end
+
+% Define a vector of indices for the observed variables. Is this really usefull?...
+BayesInfo.mf = BayesInfo.mf1;
+
+% Define the deterministic linear trend of the measurement equation.
+if DynareOptions.noconstant
+ constant = zeros(nvobs,1);
+else
+ if DynareOptions.loglinear
+ constant = log(SteadyState(BayesInfo.mfys));
+ else
+ constant = SteadyState(BayesInfo.mfys);
+ end
+end
+
+% Define the deterministic linear trend of the measurement equation.
+if BayesInfo.with_trend
+ trend_coeff = zeros(DynareDataset.info.nvobs,1);
+ t = DynareOptions.trend_coeffs;
+ for i=1:length(t)
+ if ~isempty(t{i})
+ trend_coeff(i) = evalin('base',t{i});
+ end
+ end
+ trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs];
+else
+ trend = repmat(constant,1,DynareDataset.info.ntobs);
+end
+
+% Get needed informations for kalman filter routines.
+start = DynareOptions.presample+1;
+np = size(T,1);
+mf = BayesInfo.mf;
+Y = transpose(DynareDataset.rawdata);
+
+%------------------------------------------------------------------------------
+% 3. Initial condition of the Kalman filter
+%------------------------------------------------------------------------------
+
+% 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 = BayesInfo.restrict_var_list;
+ observed_variables_idx = restrict_variables_idx(mf1);
+ state_variables_idx = restrict_variables_idx(mf0);
+ sample_size = size(Y,2);
+ number_of_state_variables = length(mf0);
+ number_of_observed_variables = length(mf1);
+ number_of_structural_innovations = length(Q);
+ init_flag = 1;
+end
+
+ReducedForm.ghx = dr.ghx(restrict_variables_idx,:);
+ReducedForm.ghu = dr.ghu(restrict_variables_idx,:);
+ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
+ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
+ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
+ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
+ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
+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;
+
+% Set initial condition for t=1
+if observation_number==1
+ switch DynareOptions.particle.initialization
+ case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model.
+ StateVectorMean = ReducedForm.constant(mf0);
+ StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',1e-12,1e-12,[],[],DynareOptions.debug);
+ case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model).
+ StateVectorMean = ReducedForm.constant(mf0);
+ old_DynareOptionsperiods = DynareOptions.periods;
+ DynareOptions.periods = 5000;
+ y_ = simult(oo_.steady_state, dr,Model,DynareOptions,DynareResults);
+ y_ = y_(state_variables_idx,2001:5000);
+ StateVectorVariance = cov(y_');
+ DynareOptions.periods = old_DynareOptionsperiods;
+ clear('old_DynareOptionsperiods','y_');
+ case 3% Initial state vector covariance is a diagonal matrix.
+ StateVectorMean = ReducedForm.constant(mf0);
+ 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/nonlinear-filters/src/spherical_radial_sigma_points.m b/nonlinear-filters/src/spherical_radial_sigma_points.m
new file mode 100644
index 000000000..754caacfa
--- /dev/null
+++ b/nonlinear-filters/src/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 (C) 2009-2012 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 .
+
+nodes = (sqrt(n)*([eye(n) -eye(n)]))' ;
+weights = (1/(2*n)) ;
diff --git a/nonlinear-filters/src/traditional_resampling.m b/nonlinear-filters/src/traditional_resampling.m
new file mode 100644
index 000000000..a9d86eb36
--- /dev/null
+++ b/nonlinear-filters/src/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 (C) 2011-2013 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 .
+
+% 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/nonlinear-filters/src/univariate_smooth_resampling.m b/nonlinear-filters/src/univariate_smooth_resampling.m
new file mode 100644
index 000000000..8238902a1
--- /dev/null
+++ b/nonlinear-filters/src/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 (C) 2012 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 .
+
+% 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/nonlinear-filters/src/unscented_sigma_points.m b/nonlinear-filters/src/unscented_sigma_points.m
new file mode 100644
index 000000000..26fcc4b5c
--- /dev/null
+++ b/nonlinear-filters/src/unscented_sigma_points.m
@@ -0,0 +1,42 @@
+function [nodes,W_m,W_c] = unscented_sigma_points(n,DynareOptions)
+%
+% 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 (C) 2009-2012 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 .
+
+lambda = (DynareOptions.particle.unscented.alpha^2)*(n+DynareOptions.particle.unscented.kappa) - n ;
+nodes = [ zeros(n,1) ( sqrt(n+lambda).*([ eye(n) -eye(n)]) ) ]' ;
+W_m = lambda/(n+lambda) ;
+W_c = W_m + (1-DynareOptions.particle.unscented.alpha^2+DynareOptions.particle.unscented.beta) ;
+temp = ones(2*n,1)/(2*(n+lambda)) ;
+W_m = [W_m ; temp] ;
+W_c = [W_c ; temp] ;
+
+