From a6009908e44e0ceb3204ba21c245d26eab56c51c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Scylla=29?=
<stephane.adjemian@univ-lemans.fr>
Date: Thu, 28 Aug 2014 17:41:38 +0200
Subject: [PATCH] Added option resampling (estimation command, non linear
filters).
---
doc/dynare.texi | 20 +++++++++++++++++++
matlab/global_initialization.m | 5 ++++-
matlab/particle/conditional_particle_filter.m | 3 +--
matlab/particle/gaussian_filter.m | 5 ++---
matlab/particle/gaussian_mixture_filter.m | 4 ++--
.../sequential_importance_particle_filter.m | 5 ++---
preprocessor/DynareBison.yy | 6 +++++-
preprocessor/DynareFlex.ll | 6 ++++++
8 files changed, 42 insertions(+), 12 deletions(-)
diff --git a/doc/dynare.texi b/doc/dynare.texi
index 6a81446893..82b0dfa3bf 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -5301,10 +5301,30 @@ variables. This is equivalent to manually listing all the observed
variables after the @code{estimation} command.
@item number_of_particles = @var{INTEGER}
+@anchor{number_of_particles}
Number of particles used when evaluating the likelihood of a non linear state space model. Default: @code{1000}.
+@item resampling = @var{OPTION}
+@anchor{resampling}
+Determines if resampling of the particles is done. Possible values for @var{OPTION} are:
+
+@table @code
+
+@item none
+No resampling.
+
+@item systematic
+Resampling at each iteration, this is the default value.
+
+@item generic
+Resampling if and only if the effective sample size is below a certain level defined by @code{resampling_threshold*number_of_particles}.
+
@end table
+
+@end table
+
+
@customhead{Note}
If no @code{mh_jscale} parameter is used in estimated_params, the
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 1a1840c915..f7afd262a4 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -236,7 +236,10 @@ particle.unscented.alpha = .01;
particle.unscented.beta = 2;
particle.unscented.kappa = 0;
% Configuration of resampling in case of particles
-particle.resampling.status = 'systematic'; % 'none', 'generic', 'smoothed'
+particle.resampling.status.systematic = 1;
+particle.resampling.status.none = 0;
+particle.resampling.status.generic = 0;
+
particle.resampling.neff_threshold = .5;
% Choice of the resampling method
particle.resampling.method1 = 'traditional' ;
diff --git a/matlab/particle/conditional_particle_filter.m b/matlab/particle/conditional_particle_filter.m
index 941e01331e..36244809b0 100644
--- a/matlab/particle/conditional_particle_filter.m
+++ b/matlab/particle/conditional_particle_filter.m
@@ -113,8 +113,7 @@ for t=1:sample_size
SumSampleWeights = sum(SampleWeights) ;
lik(t) = log(SumSampleWeights) ;
SampleWeights = SampleWeights./SumSampleWeights ;
- if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(SampleWeights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ...
- strcmp(DynareOptions.particle.resampling.status,'systematic')
+ if (DynareOptions.particle.resampling.status.generic && neff(SampleWeights)<DynareOptions.particle.resampling.neff_threshold*sample_size) || DynareOptions.particle.resampling.status.systematic
ks = ks + 1 ;
StateParticles = resample(StateParticles',SampleWeights',DynareOptions)';
SampleWeights = ones(1,number_of_particles)/number_of_particles ;
diff --git a/matlab/particle/gaussian_filter.m b/matlab/particle/gaussian_filter.m
index ebe2b1d31a..04edbcc5d0 100644
--- a/matlab/particle/gaussian_filter.m
+++ b/matlab/particle/gaussian_filter.m
@@ -146,14 +146,13 @@ for t=1:sample_size
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 && strcmpi(DynareOptions.particle.resampling.status,'generic')) || ...
- strcmpi(DynareOptions.particle.resampling.status,'systematic')
+ 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 strcmpi(DynareOptions.particle.resampling.status,'none')
+ elseif DynareOptions.particle.resampling.status.none
StateVectorMean = (sampleWeights*StateParticles)' ;
temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp'*temp - StateVectorMean*(StateVectorMean') )';
diff --git a/matlab/particle/gaussian_mixture_filter.m b/matlab/particle/gaussian_mixture_filter.m
index 8f05ea390c..0506b3c507 100644
--- a/matlab/particle/gaussian_mixture_filter.m
+++ b/matlab/particle/gaussian_mixture_filter.m
@@ -202,13 +202,13 @@ for t=1:sample_size
%estimate(t,:,1) = SampleWeights*StateParticles' ;
% Resampling if needed of required
Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
- if (Neff<.5*sample_size && strcmpi(DynareOptions.particle.resampling.status,'generic')) || strcmpi(DynareOptions.particle.resampling.status,'systematic')
+ 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 strcmpi(DynareOptions.particle.resampling.status,'none')
+ elseif DynareOptions.particle.resampling.status.none
StateVectorMean = StateParticles*sampleWeights ;
temp = sqrt(SampleWeights').*StateParticles ;
StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';
diff --git a/matlab/particle/sequential_importance_particle_filter.m b/matlab/particle/sequential_importance_particle_filter.m
index 6474f6cdab..3c86f17462 100644
--- a/matlab/particle/sequential_importance_particle_filter.m
+++ b/matlab/particle/sequential_importance_particle_filter.m
@@ -158,8 +158,7 @@ for t=1:sample_size
wtilde = weights.*exp(lnw-dfac);
lik(t) = log(sum(wtilde))+dfac;
weights = wtilde/sum(wtilde);
- if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ...
- strcmp(DynareOptions.particle.resampling.status,'systematic')
+ if (DynareOptions.particle.resampling.status.generic && neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size) || DynareOptions.particle.resampling.status.systematic
if pruning
temp = resample([tmp(mf0,:)' tmp_(mf0,:)'],weights',DynareOptions);
StateVectors = temp(:,1:number_of_state_variables)';
@@ -168,7 +167,7 @@ for t=1:sample_size
StateVectors = resample(tmp(mf0,:)',weights',DynareOptions)';
end
weights = ones(1,number_of_particles)/number_of_particles;
- elseif strcmp(DynareOptions.particle.resampling.status,'none')
+ elseif DynareOptions.particle.resampling.status.none
StateVectors = tmp(mf0,:);
if pruning
StateVectors_ = tmp_(mf0,:);
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 44b42e37ae..2362e0f1ec 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -107,7 +107,7 @@ class ParsingDriver;
%token MFS MH_CONF_SIG MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS
%token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL MCMC_JUMPING_COVARIANCE MOMENT_CALIBRATION
-%token NUMBER_OF_PARTICLES
+%token NUMBER_OF_PARTICLES RESAMPLING SYSTEMATIC GENERIC
%token <string_val> NAME
%token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
%token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS
@@ -1678,6 +1678,7 @@ estimation_options : o_datafile
| o_consider_all_endogenous
| o_consider_only_observed
| o_number_of_particles
+ | o_resampling
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -2641,6 +2642,9 @@ o_bvar_prior_train : BVAR_PRIOR_TRAIN EQUAL INT_NUMBER { driver.option_num("bvar
o_bvar_replic : BVAR_REPLIC EQUAL INT_NUMBER { driver.option_num("bvar_replic", $3); };
o_number_of_particles : NUMBER_OF_PARTICLES EQUAL INT_NUMBER { driver.option_num("particle.number_of_particles", $3); };
+o_resampling : RESAMPLING EQUAL SYSTEMATIC
+ | RESAMPLING EQUAL NONE {driver.option_num("particle.resampling.status.systematic", "0"); driver.option_num("particle.resampling.status.none", "1"); }
+ | RESAMPLING EQUAL GENERIC {driver.option_num("particle.resampling.status.systematic", "0"); driver.option_num("particle.resampling.status.generic", "1"); };
o_gsa_identification : IDENTIFICATION EQUAL INT_NUMBER { driver.option_num("identification", $3); }; /*not in doc */
o_gsa_morris : MORRIS EQUAL INT_NUMBER { driver.option_num("morris", $3); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 37386530a0..f1fab51561 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -352,7 +352,13 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>use_univariate_filters_if_singularity_is_detected {return token::USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED;}
<DYNARE_STATEMENT>hybrid {return token::HYBRID;}
<DYNARE_STATEMENT>default {return token::DEFAULT;}
+
<DYNARE_STATEMENT>number_of_particles {return token::NUMBER_OF_PARTICLES;}
+<DYNARE_STATEMENT>resampling {return token::RESAMPLING;}
+<DYNARE_STATEMENT>systematic {return token::SYSTEMATIC;}
+<DYNARE_STATEMENT>none {return token::NONE;}
+<DYNARE_STATEMENT>generic {return token::GENERIC;}
+
<DYNARE_STATEMENT>alpha {
yylval->string_val = new string(yytext);
return token::ALPHA;
--
GitLab