Merge branch 'pf_options' into 'master'

Particle filter: enable passing of option values and prepare for other samplers

Closes preprocessor#59

See merge request Dynare/dynare!1884

matlab/check_particle_filter_options.m

+function [particle_options] = check_particle_filter_options(particle_options)
+% function [particle_filter_options, options_] = check_particle_filter_options(particle_filter_options_string, options_)
+% initialization of particle filter options
+%
+% INPUTS
+%   particle_filter_options:                structure storing the options
+
+% OUTPUTS
+%   particle_filter_options:                checked particle filter options
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2021 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 <http://www.gnu.org/licenses/>.
+
+if strcmpi(particle_options.filter_algorithm, 'sis')
+    particle_options.algorithm = 'sequential_importance_particle_filter';
+elseif strcmpi(particle_options.filter_algorithm, 'apf')
+    particle_options.algorithm = 'auxiliary_particle_filter';
+elseif strcmpi(particle_options.filter_algorithm, 'gf')
+    particle_options.algorithm = 'gaussian_filter';
+elseif strcmpi(particle_options.filter_algorithm,  'gmf')
+    particle_options.algorithm = 'gaussian_mixture_filter';
+elseif strcmpi(particle_options.filter_algorithm, 'cpf')
+    particle_options.algorithm = 'conditional_particle_filter';
+elseif strcmpi(particle_options.filter_algorithm, 'nlkf')
+    particle_options.algorithm = 'nonlinear_kalman_filter';
+else
+    error(['Estimation: Unknown filter ' particle_options.filter_algorithm])
+end
+
+if ~isempty(particle_options.particle_filter_options)
+    % set default options and user defined options
+    options_list = read_key_value_string(particle_options.particle_filter_options);
+    for i=1:rows(options_list)
+        switch options_list{i,1}
+            case 'posterior_sampler'
+                if ~(strcmpi(options_list{i,2}, 'Herbst_Schorfheide') || ...
+                        strcmpi(options_list{i,2}, 'DSMH'))
+                    error(['check_particle_filter_options:: the proposal_distribution option to estimation takes either ' ...
+                        'Herbst_Schorfheide or Herbst_Schorfheide as options']);
+                else
+                    particle_options.posterior_sampler=options_list{i,2};
+                end
+            case 'initial_state_prior_std'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the initial_state_prior_std option takes a positive argument');
+                else
+                    particle_options.initial_state_prior_std=options_list{i,2};
+                end
+            case 'pruning'
+                if ~islogical(options_list{i,2})
+                    error('check_particle_filter_options:: the pruning options takes only true or false');
+                else
+                    particle_options.pruning=options_list{i,2};
+                end
+            case 'unscented_alpha'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the unscented_alpha option takes a positive argument');
+                else
+                    particle_options.unscented.alpha=options_list{i,2};
+                end
+            case 'unscented_beta'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the unscented_beta option takes a positive argument');
+                else
+                    particle_options.unscented.beta=options_list{i,2};
+                end
+            case 'unscented_kappa'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the unscented_kappa option takes a positive argument');
+                else
+                    particle_options.unscented.kappa=options_list{i,2};
+                end
+            case 'mixture_state_variables'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the mixture_state_variables option takes a positive integer');
+                else
+                    particle_options.mixture_state_variables=options_list{i,2};
+                end
+            case 'mixture_structural_shocks'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the mixture_structural_shocks option takes a positive integer');
+                else
+                    particle_options.mixture_structural_shocks=options_list{i,2};
+                end
+            case 'mixture_measurement_shocks'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the mixture_measurement_shocks option takes a positive integer');
+                else
+                    particle_options.mixture_measurement_shocks=options_list{i,2};
+                end
+            case 'liu_west_delta'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the liu_west_delta option takes a positive argument');
+                else
+                    particle_options.liu_west_delta=options_list{i,2};
+                end
+            case 'liu_west_chol_sigma_bar'
+                if options_list{i,2} <= 0
+                    error('check_particle_filter_options:: the liu_west_chol_sigma_bar option takes a positive argument');
+                else
+                    particle_options.liu_west_chol_sigma_bar=options_list{i,2};
+                end
+                
+            otherwise
+                warning(['check_particle_filter_options: Unknown option (' options_list{i,1} ')!'])
+        end
+    end
+end

matlab/default_option_values.m

@@ -246,8 +246,6 @@ particle.status = false;
 % How do we initialize the states?
 particle.initialization = 1;
 particle.initial_state_prior_std = .1;
-% Set the default order of approximation of the model (perturbation).
-particle.perturbation = 2;
 % Set the default number of particles.
 particle.number_of_particles = 5000;
 % Set the default approximation order (Smolyak)
@@ -292,7 +290,8 @@ particle.liu_west_chol_sigma_bar = .01 ;
 % Options for setting the weights in conditional particle filters.
 particle.cpf_weights_method.amisanotristani = true;
 particle.cpf_weights_method.murrayjonesparslow = false;
-% Copy ep structure in options_ global structure
+particle.particle_filter_options ='';
+% Copy particle structure in options_ global structure
 options_.particle = particle;
 options_.rwgmh.init_scale = 1e-4 ;
 options_.rwgmh.scale_chain = 1 ;

matlab/dynare_estimation_1.m

@@ -85,21 +85,7 @@ end
 if ~options_.dsge_var
     if options_.particle.status
         objective_function = str2func('non_linear_dsge_likelihood');
-        if strcmpi(options_.particle.filter_algorithm, 'sis')
-            options_.particle.algorithm = 'sequential_importance_particle_filter';
-        elseif strcmpi(options_.particle.filter_algorithm, 'apf')
-            options_.particle.algorithm = 'auxiliary_particle_filter';
-        elseif strcmpi(options_.particle.filter_algorithm, 'gf')
-            options_.particle.algorithm = 'gaussian_filter';
-        elseif strcmpi(options_.particle.filter_algorithm,  'gmf')
-            options_.particle.algorithm = 'gaussian_mixture_filter';
-        elseif strcmpi(options_.particle.filter_algorithm, 'cpf')
-            options_.particle.algorithm = 'conditional_particle_filter';
-        elseif strcmpi(options_.particle.filter_algorithm, 'nlkf')
-            options_.particle.algorithm = 'nonlinear_kalman_filter';
-        else
-            error(['Estimation: Unknown filter ' options_.particle.filter_algorithm])
-        end
+        [options_.particle] = check_particle_filter_options(options_.particle);
     else
         if options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter
             objective_function = str2func('occbin.IVF_posterior');
@@ -353,6 +339,14 @@ if ~options_.cova_compute
     stdh = NaN(length(xparam1),1);
 end
 
+if options_.particle.status && isfield(options_.particle,'posterior_sampler')
+    if strcmpi(options_.particle.posterior_sampler,'Herbst_Schorfheide')
+        Herbst_Schorfheide_sampler(objective_function,xparam1,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
+    elseif strcmpi(options_.particle.posterior_sampler,'DSMH')
+        DSMH_sampler(objective_function,xparam1,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
+    end
+end
+        
 if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
     % display results table and store parameter estimates and standard errors in results
     oo_ = display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, prior_dist_names, 'Posterior', 'posterior');

tests/particle/dsge_base2.mod

 // DGP
 @#ifndef RISKY_CALIBRATION
-@#define RISKY_CALIBRATION = 1
+    @#define RISKY_CALIBRATION = 1
 @#endif
 @#ifndef EXTREME_CALIBRATION
-@#define EXTREME_CALIBRATION = 0
+    @#define EXTREME_CALIBRATION = 0
 @#endif
 @#ifndef BENCHMARK_CALIBRATION
-@#define BENCHMARK_CALIBRATION = 0
+    @#define BENCHMARK_CALIBRATION = 0
 @#endif
 
 // ALGORITHM
 @#ifndef LINEAR_KALMAN
-@#define LINEAR_KALMAN = 0
+    @#define LINEAR_KALMAN = 0
 @#endif
 @#ifndef NON_LINEAR_KALMAN
-@#define NON_LINEAR_KALMAN = 1
+    @#define NON_LINEAR_KALMAN = 1
 @#endif
 @#ifndef ALGO_SIR
-@#define ALGO_SIR = 0
+    @#define ALGO_SIR = 0
 @#endif
 @#ifndef ALGO_SISmoothR
-@#define ALGO_SISmoothR = 0
+    @#define ALGO_SISmoothR = 0
 @#endif
 @#ifndef ALGO_APF
-@#define ALGO_APF = 0
+    @#define ALGO_APF = 0
 @#endif
 @#ifndef ALGO_CPF
-@#define ALGO_CPF = 0
+    @#define ALGO_CPF = 0
 @#endif
 @#ifndef ALGO_GPF
-@#define ALGO_GPF = 0
+    @#define ALGO_GPF = 0
 @#endif
 @#ifndef ALGO_GCF
-@#define ALGO_GCF = 0
+    @#define ALGO_GCF = 0
 @#endif
 @#ifndef ALGO_GUF
-@#define ALGO_GUF = 0
+    @#define ALGO_GUF = 0
 @#endif
 @#ifndef ALGO_GMPF
-@#define ALGO_GMPF = 0
+    @#define ALGO_GMPF = 0
 @#endif
 @#ifndef ALGO_GMCF
-@#define ALGO_GMCF = 0
+    @#define ALGO_GMCF = 0
 @#endif
 @#ifndef ALGO_ONLINE_1
-@#define ALGO_ONLINE_1 = 0
+    @#define ALGO_ONLINE_1 = 0
 @#endif
 @#ifndef ALGO_ONLINE_2
-@#define ALGO_ONLINE_2 = 0
+    @#define ALGO_ONLINE_2 = 0
 @#endif
 @#ifndef MCMC
-@#define MCMC = 0
+    @#define MCMC = 0
 @#endif
 
 var k A c l i y;
@@ -84,8 +84,6 @@ steady_state_model;
 
 end;
 
-
-
 shocks;
 var e_a; stderr 0.035;
 end;
@@ -207,13 +205,11 @@ options_.threads.local_state_space_iteration_2 = 4;
 @#endif
 
 @#if ALGO_ONLINE_2
-  options_.particle.liu_west_delta = 0.9 ;
-  estimation(order=2,number_of_particles=1000,mode_compute=11,mh_replic=0);
+  estimation(order=2,number_of_particles=1000,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
 @#endif
 
 @#if ALGO_ONLINE_1
-  options_.particle.liu_west_delta = 0.9 ;
-  estimation(order=1,number_of_particles=1000,mode_compute=11,mh_replic=0);
+  estimation(order=1,number_of_particles=1000,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
 @#endif
 
 @#if MCMC