Added option resampling (estimation command, non linear filters).

time-shift
Stéphane Adjemian (Scylla) 2014-08-28 17:41:38 +02:00
parent 77a29e95c1
commit a6009908e4
8 changed files with 42 additions and 12 deletions

View File

@ -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

View File

@ -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' ;

View File

@ -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 ;

View File

@ -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') )';

View File

@ -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') )';

View File

@ -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,:);

View File

@ -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); };

View File

@ -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;