From 79b6184d876d2a43107c67dd4094351b8c34462e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Argos=29?=
<stepan@adjemian.eu>
Date: Mon, 24 Mar 2025 12:30:07 +0100
Subject: [PATCH] Make SIS (particle filter) compatible with linear models.
---
matlab/estimation/dynare_estimation_1.m | 62 +++++++++----------
.../estimation/non_linear_dsge_likelihood.m | 29 +++++----
.../sequential_importance_particle_filter.m | 51 ++++++++-------
3 files changed, 73 insertions(+), 69 deletions(-)
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index e82e4389dd..1c42b50992 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -12,7 +12,7 @@ function dynare_estimation_1(var_list_,dname)
% SPECIAL REQUIREMENTS
% none
-% Copyright © 2003-2024 Dynare Team
+% Copyright © 2003-2025 Dynare Team
%
% This file is part of Dynare.
%
@@ -68,35 +68,31 @@ else
end
% Set particle filter flag.
-if options_.order > 1
- if options_.particle.status
- skipline()
- disp('Estimation using a non linear filter!')
- skipline()
- if ~options_.nointeractive && ismember(options_.mode_compute,[1,3,4]) && ~strcmpi(options_.particle.filter_algorithm,'gf')% Known gradient-based optimizers
- disp('You are using a gradient-based mode-finder. Particle filtering introduces discontinuities in the')
- disp('objective function w.r.t the parameters. Thus, should use a non-gradient based optimizer.')
- fprintf('\nPlease choose a mode-finder:\n')
- fprintf('\t 0 - Continue using gradient-based method (it is most likely that you will no get any sensible result).\n')
- fprintf('\t 6 - Monte Carlo based algorithm\n')
- fprintf('\t 7 - Nelder-Mead simplex based optimization routine (Matlab optimization toolbox required)\n')
- fprintf('\t 8 - Nelder-Mead simplex based optimization routine (Dynare''s implementation)\n')
- fprintf('\t 9 - CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm\n')
- choice = [];
- while isempty(choice)
- choice = input('Please enter your choice: ');
- if isnumeric(choice) && isint(choice) && ismember(choice,[0 6 7 8 9])
- if choice
- options_.mode_compute = choice;
- end
- else
- fprintf('\nThis is an invalid choice (you have to choose between 0, 6, 7, 8 and 9).\n')
- choice = [];
+if options_.particle.status
+ skipline()
+ disp('Estimation using a non linear filter!')
+ skipline()
+ if ~options_.nointeractive && ismember(options_.mode_compute,[1,3,4]) && ~strcmpi(options_.particle.filter_algorithm,'gf')% Known gradient-based optimizers
+ disp('You are using a gradient-based mode-finder. Particle filtering introduces discontinuities in the')
+ disp('objective function w.r.t the parameters. Thus, should use a non-gradient based optimizer.')
+ fprintf('\nPlease choose a mode-finder:\n')
+ fprintf('\t 0 - Continue using gradient-based method (it is most likely that you will no get any sensible result).\n')
+ fprintf('\t 6 - Monte Carlo based algorithm\n')
+ fprintf('\t 7 - Nelder-Mead simplex based optimization routine (Matlab optimization toolbox required)\n')
+ fprintf('\t 8 - Nelder-Mead simplex based optimization routine (Dynare''s implementation)\n')
+ fprintf('\t 9 - CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm\n')
+ choice = [];
+ while isempty(choice)
+ choice = input('Please enter your choice: ');
+ if isnumeric(choice) && isint(choice) && ismember(choice,[0 6 7 8 9])
+ if choice
+ options_.mode_compute = choice;
end
+ else
+ fprintf('\nThis is an invalid choice (you have to choose between 0, 6, 7, 8 and 9).\n')
+ choice = [];
end
end
- else
- error('For estimating the model with a second order approximation using a non linear filter, one should have options_.particle.status=true;')
end
end
@@ -579,14 +575,14 @@ end
%Run and store classical smoother if needed
if options_.smoother && ... %Bayesian smoother requested before
- (any(bayestopt_.pshape > 0) && options_.mh_replic || ... % Bayesian with MCMC run
- any(bayestopt_.pshape > 0) && options_.load_mh_file) % Bayesian with loaded MCMC
- % nothing to do
+ (any(bayestopt_.pshape > 0) && options_.mh_replic || ... % Bayesian with MCMC run
+ any(bayestopt_.pshape > 0) && options_.load_mh_file) % Bayesian with loaded MCMC
+ % nothing to do
elseif options_.partial_information ||...
- options_.order>1 %no particle smoother
- % smoothing not yet supported
+ options_.order>1 %no particle smoother
+ % smoothing not yet supported
else
- %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
+ % ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
end
diff --git a/matlab/estimation/non_linear_dsge_likelihood.m b/matlab/estimation/non_linear_dsge_likelihood.m
index a7b13c1798..0ee495f431 100644
--- a/matlab/estimation/non_linear_dsge_likelihood.m
+++ b/matlab/estimation/non_linear_dsge_likelihood.m
@@ -120,7 +120,10 @@ state_variables_idx = restrict_variables_idx(mf0);
number_of_state_variables = length(mf0);
ReducedForm.steadystate = dr.ys(dr.order_var(restrict_variables_idx));
-ReducedForm.constant = ReducedForm.steadystate + .5*dr.ghs2(restrict_variables_idx);
+ReducedForm.constant = ReducedForm.steadystate;
+if options_.order==2
+ ReducedForm.constant = ReducedForm.constant + .5*dr.ghs2(restrict_variables_idx);
+end
ReducedForm.state_variables_steady_state = dr.ys(dr.order_var(state_variables_idx));
ReducedForm.Q = Q;
ReducedForm.H = H;
@@ -138,17 +141,19 @@ else
ReducedForm.use_k_order_solver = false;
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.ghs2 = dr.ghs2(restrict_variables_idx,:);
- if options_.order==3
- ReducedForm.ghxxx = dr.ghxxx(restrict_variables_idx,:);
- ReducedForm.ghuuu = dr.ghuuu(restrict_variables_idx,:);
- ReducedForm.ghxxu = dr.ghxxu(restrict_variables_idx,:);
- ReducedForm.ghxuu = dr.ghxuu(restrict_variables_idx,:);
- ReducedForm.ghxss = dr.ghxss(restrict_variables_idx,:);
- ReducedForm.ghuss = dr.ghuss(restrict_variables_idx,:);
+ if options_.order>1
+ ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:);
+ ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:);
+ ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:);
+ ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:);
+ if options_.order==3
+ ReducedForm.ghxxx = dr.ghxxx(restrict_variables_idx,:);
+ ReducedForm.ghuuu = dr.ghuuu(restrict_variables_idx,:);
+ ReducedForm.ghxxu = dr.ghxxu(restrict_variables_idx,:);
+ ReducedForm.ghxuu = dr.ghxuu(restrict_variables_idx,:);
+ ReducedForm.ghxss = dr.ghxss(restrict_variables_idx,:);
+ ReducedForm.ghuss = dr.ghuss(restrict_variables_idx,:);
+ end
end
end
diff --git a/matlab/nonlinear-filters/sequential_importance_particle_filter.m b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
index 0ac02e8f7d..2bfb281676 100644
--- a/matlab/nonlinear-filters/sequential_importance_particle_filter.m
+++ b/matlab/nonlinear-filters/sequential_importance_particle_filter.m
@@ -2,7 +2,7 @@ function [LIK,lik] = sequential_importance_particle_filter(ReducedForm,Y,start,P
% Evaluates the likelihood of a nonlinear model with a particle filter (optionally with resampling).
-% Copyright © 2011-2022 Dynare Team
+% Copyright © 2011-2025 Dynare Team
%
% This file is part of Dynare (particles module).
%
@@ -47,20 +47,21 @@ else
% Set local state space model (first order approximation).
ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu;
-
% Set local state space model (second order approximation).
- ghxx = ReducedForm.ghxx;
- ghuu = ReducedForm.ghuu;
- ghxu = ReducedForm.ghxu;
- ghs2 = ReducedForm.ghs2;
- if order == 3
- % Set local state space model (third order approximation).
- ghxxx = ReducedForm.ghxxx;
- ghuuu = ReducedForm.ghuuu;
- ghxxu = ReducedForm.ghxxu;
- ghxuu = ReducedForm.ghxuu;
- ghxss = ReducedForm.ghxss;
- ghuss = ReducedForm.ghuss;
+ if order>1
+ ghxx = ReducedForm.ghxx;
+ ghuu = ReducedForm.ghuu;
+ ghxu = ReducedForm.ghxu;
+ ghs2 = ReducedForm.ghs2;
+ if order == 3
+ % Set local state space model (third order approximation).
+ ghxxx = ReducedForm.ghxxx;
+ ghuuu = ReducedForm.ghuuu;
+ ghxxu = ReducedForm.ghxxu;
+ ghxuu = ReducedForm.ghxuu;
+ ghxss = ReducedForm.ghxss;
+ ghuss = ReducedForm.ghuss;
+ end
end
end
@@ -126,20 +127,22 @@ for t=1:sample_size
if ReducedForm.use_k_order_solver
tmp = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr);
else
- if order == 2
+ switch order
+ case 1
+ tmp = NaN(rows(ghx), columns(yhat));
+ parfor i=1:columns(yhat)
+ tmp(:,i) = constant + ghx*yhat(:,i) + ghu*epsilon(:,i)
+ end
+ case 2
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
- elseif order == 3
+ case 3
tmp = local_state_space_iteration_3(yhat, epsilon, ghx, ghu, ghxx, ghuu, ghxu, ghs2, ghxxx, ghuuu, ghxxu, ghxuu, ghxss, ghuss, steadystate, ThreadsOptions.local_state_space_iteration_3, pruning);
- else
+ otherwise
error('Order > 3: use_k_order_solver should be set to true');
end
end
end
- %PredictedObservedMean = tmp(mf1,:)*transpose(weights);
PredictionError = bsxfun(@minus,Y(:,t),tmp(ReducedForm.mf1,:));
- %dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
- %PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H;
- %PredictedObservedVariance = H;
if rcond(H) > 1e-16
lnw = -.5*(const_lik+sum(PredictionError.*(H\PredictionError),1));
else
@@ -150,13 +153,13 @@ for t=1:sample_size
wtilde = weights.*exp(lnw-dfac);
lik(t) = log(sum(wtilde))+dfac;
weights = wtilde/sum(wtilde);
- if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
+ if (ParticleOptions.resampling.status.generic && neff(weights)<ParticleOptions.resampling.threshold*number_of_particles) || ParticleOptions.resampling.status.systematic
if pruning
- temp = resample([tmp(ReducedForm.mf0,:)' tmp_(mf0_,:)'],weights',ParticleOptions);
+ temp = resample([tmp(ReducedForm.mf0,:)' tmp_(mf0_,:)'], weights', ParticleOptions);
StateVectors = temp(:,1:number_of_state_variables)';
StateVectors_ = temp(:,number_of_state_variables+1:end)';
else
- StateVectors = resample(tmp(ReducedForm.mf0,:)',weights',ParticleOptions)';
+ StateVectors = resample(tmp(ReducedForm.mf0,:)', weights', ParticleOptions)';
end
weights = ones(1,number_of_particles)/number_of_particles;
elseif ParticleOptions.resampling.status.none
--
GitLab