diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 9422d030165c1b6e5d7714cc1fdd7db7f2c9ef65..a1d501d980702871ec3c5b698f2e98cda5b330e1 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -289,9 +289,6 @@ particle.resampling.number_of_partitions = 200;
particle.mixture_state_variables = 5 ;
particle.mixture_structural_shocks = 1 ;
particle.mixture_measurement_shocks = 1 ;
-% Online approach
-particle.liu_west_delta = 0.99 ;
-particle.liu_west_max_resampling_tries = 5000;
% Options for setting the weights in conditional particle filters.
particle.cpf_weights_method.amisanotristani = true;
particle.cpf_weights_method.murrayjonesparslow = false;
@@ -517,6 +514,11 @@ options_.posterior_sampler_options.dsmh.particles = 20000 ;
options_.posterior_sampler_options.dsmh.alpha0 = 0.2 ;
options_.posterior_sampler_options.dsmh.alpha1 = 0.3 ;
options_.posterior_sampler_options.dsmh.tau = 10 ;
+% Liu and West online Sampler
+options_.posterior_sampler_options.online.particles= 5000 ;
+options_.posterior_sampler_options.online.liu_west_delta = 0.99 ;
+options_.posterior_sampler_options.online.liu_west_max_resampling_tries = 5000 ;
+options_.posterior_sampler_options.online.second_resampling = false;
options_.trace_plot_ma = 200;
options_.mh_autocorrelation_function_size = 30;
@@ -660,6 +662,12 @@ particleswarm.UseParallel = false;
particleswarm.UseVectorized = false;
options_.particleswarm = particleswarm;
+% One step
+one_step.gradient = [] ;
+one_step.hessian = [] ;
+one_step.outer_product = 0 ;
+options_.one_step= one_step ;
+
% prior analysis
options_.prior_mc = 20000;
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 5b2f208c85d563087a0011c48df2302b77b827a9..8d44e1caaa2756125d7f5d7fbc3b279cfe53d778 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -53,6 +53,7 @@ p = {'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ...
'/ep/' ; ...
'/estimation/'; ...
'/estimation/smc/'; ...
+ '/estimation/online/'; ...
'/estimation/resampler/'; ...
'/kalman/' ; ...
'/kalman/likelihood' ; ...
diff --git a/matlab/estimation/GetAllPosteriorDraws.m b/matlab/estimation/GetAllPosteriorDraws.m
index 5c1847cf0283ecf87da59089dd78f843a7306ab4..fdea14c93f56029116a2a2887379f0d34a63a3e2 100644
--- a/matlab/estimation/GetAllPosteriorDraws.m
+++ b/matlab/estimation/GetAllPosteriorDraws.m
@@ -51,6 +51,18 @@ if issmc(options_)
else
draws = posterior.tlogpostkernel;
end
+ elseif isonline(options_)
+ % Load draws from the posterior distribution
+ posterior = load(sprintf('%s/online/parameters_particles_final.mat', dname));
+ if column>0 || strcmp(column,'all')
+ if strcmp(column,'all')
+ draws = transpose(posterior.param);
+ else
+ draws = transpose(posterior.param(column,:));
+ end
+ else
+ draws=NaN(size(posterior.param,2),1);
+ end
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
diff --git a/matlab/estimation/GetPosteriorMeanVariance.m b/matlab/estimation/GetPosteriorMeanVariance.m
index 9a1c1205129f6cc24b7a0745b15e92ec94db039b..7f73899cc863c5b8d51f88a388f5fa755cc816aa 100644
--- a/matlab/estimation/GetPosteriorMeanVariance.m
+++ b/matlab/estimation/GetPosteriorMeanVariance.m
@@ -36,6 +36,12 @@ if issmc(options_)
mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
% Compute the posterior covariance
variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
+ elseif isonline(options_)
+ posterior = load(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
+ % Compute the posterior mean
+ mean = sum(posterior.param, 2)/size(posterior.param, 2);
+ % Compute the posterior covariance
+ variance = (posterior.param-mean)*(posterior.param-mean)'/size(posterior.param, 2);
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
diff --git a/matlab/estimation/GetPosteriorParametersStatistics.m b/matlab/estimation/GetPosteriorParametersStatistics.m
index 208a37779290ca18a3daa690cfb706c9f7e419b8..ebaf93666523353a3a2df5c5d54f5f594743395d 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -55,9 +55,11 @@ skipline(2)
disp('ESTIMATION RESULTS')
skipline()
-if ishssmc(options_)
+if ishssmc(options_)
dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
+elseif isonline(options_)
+ hpd_draws = num_draws;
elseif isdime(options_)
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
else
@@ -87,7 +89,7 @@ if np
disp(tit2)
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
- if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_)
+ if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_) || isonline(options_)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
name = bayestopt_.name{ip};
diff --git a/matlab/estimation/check_posterior_analysis_data.m b/matlab/estimation/check_posterior_analysis_data.m
index 6de65d2de26ed4380f63f3785b3ddc4dc5f5a467..6b9008d4422899c063c90cb77b20a86ab5d467cf 100644
--- a/matlab/estimation/check_posterior_analysis_data.m
+++ b/matlab/estimation/check_posterior_analysis_data.m
@@ -45,6 +45,8 @@ if ~issmc(options_)
else
if ishssmc(options_)
[MetropolisFolder, info] = CheckPath('hssmc',M_.dname);
+ elseif isonline(options_)
+ [MetropolisFolder, info] = CheckPath('online',M_.dname);
elseif isdime(options_)
[MetropolisFolder, info] = CheckPath('dime',M_.dname);
else
@@ -76,6 +78,10 @@ else
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
mhdate = get_date_of_a_file(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+ elseif isonline(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
+ mhdate = get_date_of_a_file(sprintf('%s/online/parameters_particles_final.mat', M_.dname, length(pfiles), length(pfiles)));
elseif isdime(options_)
mhdate = get_date_of_a_file(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
end
diff --git a/matlab/estimation/check_posterior_sampler_options.m b/matlab/estimation/check_posterior_sampler_options.m
index e93510075baf1f6956c47beb1a5c148f9eefc3a7..fe93f77a69ec5436c5ee84ed80b28f28c25e25e3 100644
--- a/matlab/estimation/check_posterior_sampler_options.m
+++ b/matlab/estimation/check_posterior_sampler_options.m
@@ -374,125 +374,150 @@ if init
posterior_sampler_options.WR=[];
end
- case 'hssmc'
+ case 'hssmc'
- % default options
- posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
+ % default options
+ posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
- % user defined options
- if ~isempty(options_.posterior_sampler_options.sampling_opt)
- options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
- for i=1:rows(options_list)
- switch options_list{i,1}
+ % user defined options
+ if ~isempty(options_.posterior_sampler_options.sampling_opt)
+ options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
+ for i=1:rows(options_list)
+ switch options_list{i,1}
case 'target'
- posterior_sampler_options.target = options_list{i,2};
+ posterior_sampler_options.target = options_list{i,2};
case 'steps'
- posterior_sampler_options.steps = options_list{i,2};
+ posterior_sampler_options.steps = options_list{i,2};
case 'scale'
- posterior_sampler_options.scale = options_list{i,2};
+ posterior_sampler_options.scale = options_list{i,2};
case 'particles'
- posterior_sampler_options.particles = options_list{i,2};
+ posterior_sampler_options.particles = options_list{i,2};
case 'lambda'
- posterior_sampler_options.lambda = options_list{i,2};
+ posterior_sampler_options.lambda = options_list{i,2};
otherwise
- warning(['hssmc: Unknown option (' options_list{i,1} ')!'])
- end
- end
- end
-
- options_.mode_compute = 0;
- options_.cova_compute = 0;
- options_.mh_replic = 0;
- options_.mh_posterior_mode_estimation = false;
-
- case 'dime_mcmc'
-
- % default options
- posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dime);
-
- % user defined options
- if ~isempty(options_.posterior_sampler_options.sampling_opt)
- options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
- for i=1:rows(options_list)
- switch options_list{i,1}
+ warning(['hssmc: Unknown option (' options_list{i,1} ')!'])
+ end
+ end
+ end
+
+ options_.mode_compute = 0;
+ options_.cova_compute = 0;
+ options_.mh_replic = 0;
+ options_.mh_posterior_mode_estimation = false;
+
+ case 'dime_mcmc'
+
+ % default options
+ posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dime);
+
+ if ~isempty(options_.posterior_sampler_options.sampling_opt)
+ options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
+ for i=1:rows(options_list)
+ switch options_list{i,1}
case 'nchain'
- posterior_sampler_options.nchain = options_list{i,2};
+ posterior_sampler_options.nchain = options_list{i,2};
case 'niter'
- posterior_sampler_options.niter = options_list{i,2};
+ posterior_sampler_options.niter = options_list{i,2};
case 'parallel'
- posterior_sampler_options.parallel = options_list{i,2};
+ posterior_sampler_options.parallel = options_list{i,2};
case 'tune'
- posterior_sampler_options.tune = options_list{i,2};
+ posterior_sampler_options.tune = options_list{i,2};
case 'aimh_prob'
- posterior_sampler_options.aimh_prob = options_list{i,2};
+ posterior_sampler_options.aimh_prob = options_list{i,2};
case 'sigma'
- posterior_sampler_options.sigma = options_list{i,2};
+ posterior_sampler_options.sigma = options_list{i,2};
case 'gamma'
- posterior_sampler_options.gamma = options_list{i,2};
+ posterior_sampler_options.gamma = options_list{i,2};
case 'df_proposal_dist'
- posterior_sampler_options.df_proposal_dist = options_list{i,2};
+ posterior_sampler_options.df_proposal_dist = options_list{i,2};
case 'rho'
- posterior_sampler_options.rho = options_list{i,2};
+ posterior_sampler_options.rho = options_list{i,2};
otherwise
- warning(['dime: Unknown option (' options_list{i,1} ')!'])
- end
- end
- end
-
- if ~isfield(posterior_sampler_options,'tune')
- posterior_sampler_options.tune = ceil(50000 / posterior_sampler_options.nchain);
- dprintf('dime: setting number of ensemble iterations to keep (`tune`) to %d.', posterior_sampler_options.tune)
- end
- options_.mode_compute = 0;
- options_.cova_compute = 0;
- options_.mh_replic = 0;
- options_.mh_posterior_mode_estimation = false;
+ warning(['dime: Unknown option (' options_list{i,1} ')!'])
+ end
+ end
+ end
+
+ if ~isfield(posterior_sampler_options,'tune')
+ posterior_sampler_options.tune = ceil(50000 / posterior_sampler_options.nchain);
+ dprintf('dime: setting number of ensemble iterations to keep (`tune`) to %d.', posterior_sampler_options.tune)
+ end
+
+ case 'online'
+
+ % default options
+ posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.online);
+
+ % user defined options
+ if ~isempty(options_.posterior_sampler_options.sampling_opt)
+ options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
+ for i=1:rows(options_list)
+ switch options_list{i,1}
+ case 'particles'
+ posterior_sampler_options.particles = options_list{i,2};
+ case 'liu_west_delta'
+ posterior_sampler_options.liu_west_delta = options_list{i,2};
+ case 'liu_west_max_resampling_tries'
+ posterior_sampler_options.liu_west_max_resampling_tries = options_list{i,2};
+ case 'second_resampling'
+ posterior_sampler_options.second_resampling = options_list{i,2};
+ otherwise
+ warning(['online: Unknown option (' options_list{i,1} ')!'])
+ end
+ end
+ end
+
+ options_.mode_compute = 0;
+ options_.cova_compute = 0;
+ options_.mh_replic = 0;
+ options_.mh_posterior_mode_estimation = false;
case 'dsmh'
- % default options
- posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dsmh);
-
- % user defined options
- if ~isempty(options_.posterior_sampler_options.sampling_opt)
- options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
- for i=1:rows(options_list)
- switch options_list{i,1}
- case 'proposal_distribution'
- if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
- strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
- error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
- 'rand_multivariate_student or rand_multivariate_normal as options']);
- else
- posterior_sampler_options.proposal_distribution=options_list{i,2};
- end
- case 'student_degrees_of_freedom'
- if options_list{i,2} <= 0
- error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
- else
- posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
- end
- case 'save_tmp_file'
- posterior_sampler_options.save_tmp_file = options_list{i,2};
- case 'number_of_particles'
- posterior_sampler_options.particles = options_list{i,2};
- otherwise
- warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
- end
- end
- end
-
- options_.mode_compute = 0;
- options_.cova_compute = 0;
- options_.mh_replic = 0;
- options_.mh_posterior_mode_estimation = true;
+ % default options
+ posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dsmh);
+
+ % user defined options
+ if ~isempty(options_.posterior_sampler_options.sampling_opt)
+ options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
+ for i=1:rows(options_list)
+ switch options_list{i,1}
+ case 'proposal_distribution'
+ if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
+ strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
+ error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
+ 'rand_multivariate_student or rand_multivariate_normal as options']);
+ else
+ posterior_sampler_options.proposal_distribution=options_list{i,2};
+ end
+ case 'student_degrees_of_freedom'
+ if options_list{i,2} <= 0
+ error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
+ else
+ posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
+ end
+ case 'save_tmp_file'
+ posterior_sampler_options.save_tmp_file = options_list{i,2};
+ case 'number_of_particles'
+ posterior_sampler_options.particles = options_list{i,2};
+ otherwise
+ warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
+ end
+ end
+ end
+
+ options_.mode_compute = 0;
+ options_.cova_compute = 0;
+ options_.mh_replic = 0;
+ options_.mh_posterior_mode_estimation = true;
+
+ otherwise
+ error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
+ end
- otherwise
- error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
- end
+ return
- return
-end
+end % if init
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index 6e2bd2bee0d2b881644be389f367d46fb99acd8e..12c0444b5aa1c1b18960e95752a7ff28d83ce1ef 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -237,11 +237,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
optimizer_vec = [options_.mode_compute;num2cell(options_.additional_optimizer_steps)];
for optim_iter = 1:length(optimizer_vec)
current_optimizer = optimizer_vec{optim_iter};
- if current_optimizer~=11
- [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
- else
- [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
- end
+ [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
if isnumeric(current_optimizer)
@@ -271,7 +267,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
- end
+ end
hh = reshape(hh, nx, nx);
elseif isnumeric(current_optimizer) && isequal(current_optimizer,5)
% other numerical hessian options available with optimizer
@@ -398,7 +394,11 @@ end
% Run SMC sampler.
%
-if ishssmc(options_)
+if isonline(options_)
+ [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
+ options_.posterior_sampler_options.current_options = posterior_sampler_options;
+ online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_);
+elseif ishssmc(options_)
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
options_.posterior_sampler_options.current_options = posterior_sampler_options;
oo_.MarginalDensity.hssmc = hssmc(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
@@ -471,7 +471,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
oo_.lprob = trace_plot_dime(options_, M_);
end
% Estimation of the marginal density from the Mh draws:
- if ishssmc(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
+ if ishssmc(options_) || isonline(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
if ~issmc(options_)
[~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
end
diff --git a/matlab/estimation/get_posterior_folder_name.m b/matlab/estimation/get_posterior_folder_name.m
index 25b86f1c8b6375b691dfff125cbee03b0752aafd..fc0fd0b8e463d03b60347dc27fc6f5b1a95112d6 100644
--- a/matlab/estimation/get_posterior_folder_name.m
+++ b/matlab/estimation/get_posterior_folder_name.m
@@ -36,6 +36,8 @@ else
folder_name='hssmc';
elseif isdime(options_)
folder_name='dime';
+ elseif isonline(options_)
+ folder_name='online';
else
error('get_posterior_folder_name:: case should not happen. Please contact the developers')
end
diff --git a/matlab/estimation/online/isonline.m b/matlab/estimation/online/isonline.m
new file mode 100644
index 0000000000000000000000000000000000000000..b44670ea524b57d299f6d4349180e99c94a0fdf6
--- /dev/null
+++ b/matlab/estimation/online/isonline.m
@@ -0,0 +1,25 @@
+function bool = isonline(options_)
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+bool = false;
+if isfield(options_, 'posterior_sampler_options')
+ if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'online')
+ bool = true;
+ end
+end
diff --git a/matlab/nonlinear-filters/online_auxiliary_filter.m b/matlab/estimation/online/online_auxiliary_filter.m
similarity index 80%
rename from matlab/nonlinear-filters/online_auxiliary_filter.m
rename to matlab/estimation/online/online_auxiliary_filter.m
index e392474ea5681c7e6bc319aec860cfb0ef88d87c..ebe1394ff1d2e18c98bfd4e4e69aa3f5f848e7f9 100644
--- a/matlab/nonlinear-filters/online_auxiliary_filter.m
+++ b/matlab/estimation/online/online_auxiliary_filter.m
@@ -1,5 +1,5 @@
-function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
-% [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
+function online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
+% online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
%
% INPUTS
@@ -12,16 +12,7 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
% - bayestopt_ [struct] Prior definition.
% - oo_ [struct] Results.
%
-% OUTPUTS
-% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters).
-% - pmode [double] n×1 vector, mode of the particles at the end of the sample (for the parameters).
-% - pmedian [double] n×1 vector, median of the particles at the end of the sample (for the parameters).
-% - pstdev [double] n×1 vector, st. dev. of the particles at the end of the sample (for the parameters).
-% - p025 [double] n×1 vector, 2.5 percent of the particles are below p025(i) for i=1,…,n.
-% - p975 [double] n×1 vector, 97.5 percent of the particles are below p975(i) for i=1,…,n.
-% - covariance [double] n×n matrix, covariance of the particles at the end of the sample.
-
-% Copyright © 2013-2023 Dynare Team
+% Copyright © 2013-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -41,8 +32,11 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
% Set seed for randn().
options_ = set_dynare_seed_local_options(options_,'default');
pruning = options_.particle.pruning;
-second_resample = options_.particle.resampling.status.systematic;
variance_update = true;
+online_opt = options_.posterior_sampler_options.current_options;
+
+% Set location for the simulated particles.
+SimulationFolder = CheckPath('online', M_.dname);
bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
@@ -52,13 +46,13 @@ bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb an
order = options_.order;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
-number_of_particles = options_.particle.number_of_particles;
+number_of_particles = online_opt.particles;
number_of_parameters = size(xparam1,1);
Y = dataset_.data;
sample_size = size(Y,1);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
-liu_west_delta = options_.particle.liu_west_delta;
+liu_west_delta = online_opt.liu_west_delta;
% Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean;
@@ -99,7 +93,6 @@ 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);
-mode_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);
@@ -194,7 +187,6 @@ for t=1:sample_size
% Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails.
z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ;
ddl = 3 ;
- %tau_tilde(i) = weights(i).*(tpdf(z,3*ones(size(z)))+1e-99) ;
tau_tilde(i) = weights(i)*(exp(gammaln((ddl + 1) / 2) - gammaln(ddl/2))/(sqrt(ddl*pi)*(1 + (z^2)/ddl)^((ddl + 1)/2))+1e-99) ;
else
tau_tilde(i) = 0 ;
@@ -202,7 +194,7 @@ for t=1:sample_size
end
% particles selection
tau_tilde = tau_tilde/sum(tau_tilde);
- indx = resample(0, tau_tilde', options_.particle);
+ indx = kitagawa(tau_tilde');
StateVectors = StateVectors(:,indx);
xparam = fore_xparam(:,indx);
if pruning
@@ -214,7 +206,7 @@ for t=1:sample_size
for i=1:number_of_particles
info = 12042009;
counter=0;
- while info(1) && counter <options_.particle.liu_west_max_resampling_tries
+ while info(1) && counter <online_opt.liu_west_max_resampling_tries
counter=counter+1;
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
@@ -298,12 +290,12 @@ for t=1:sample_size
wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1)));
end
end
- if counter==options_.particle.liu_west_max_resampling_tries
- fprintf('\nLiu & West particle filter: I haven''t been able to solve the model in %u tries.\n',options_.particle.liu_west_max_resampling_tries)
- fprintf('Liu & West particle filter: The last error message was: %s\n',get_error_message(info))
- fprintf('Liu & West particle filter: You can try to increase liu_west_max_resampling_tries, but most\n')
- fprintf('Liu & West particle filter: likely there is an issue with the model.\n')
- error('Liu & West particle filter: unable to solve the model.')
+ if counter==online_opt.liu_west_max_resampling_tries
+ fprintf('\nLiu & West online filter: I haven''t been able to solve the model in %u tries.\n',online_opt.liu_west_max_resampling_tries)
+ fprintf('Liu & West online filter: The last error message was: %s\n',get_error_message(info))
+ fprintf('Liu & West online filter: You can try to increase liu_west_max_resampling_tries, but most\n')
+ fprintf('Liu & West online filter: likely there is an issue with the model.\n')
+ error('Liu & West online filter: unable to solve the model.')
end
end
end
@@ -313,10 +305,8 @@ for t=1:sample_size
variance_update = false;
end
% final resampling (not advised)
- if second_resample
- [~, idmode] = max(weights);
- mode_xparam(:,t) = xparam(:,idmode);
- indx = resample(0, weights,options_.particle);
+ if online_opt.second_resampling
+ indx = kitagawa(weights);
StateVectors = StateVectors(:,indx) ;
if pruning
StateVectors_ = StateVectors_(:,indx);
@@ -333,10 +323,11 @@ for t=1:sample_size
median_xparam(i,t) = temp(0.5*number_of_particles);
ub95_xparam(i,t) = temp(0.975*number_of_particles);
end
- end
- if second_resample
- [~, idmode] = max(weights);
- mode_xparam(:,t) = xparam(:,idmode);
+ if t==sample_size
+ param = xparam ;
+ save(sprintf('%s%sparameters_particles_final.mat', SimulationFolder, filesep()), 'param');
+ end
+ else
mean_xparam(:,t) = xparam*(weights');
mat_var_cov = bsxfun(@minus, xparam,mean_xparam(:,t));
mat_var_cov = mat_var_cov*(bsxfun(@times, mat_var_cov, weights)');
@@ -362,22 +353,21 @@ for t=1:sample_size
end
end
end
+ if t==sample_size
+ param = xparam(:,kitagawa(weights));
+ save(sprintf('%s%sparameters_particles_final.mat', SimulationFolder, filesep()), 'param');
+ end
end
+ save(sprintf('%s%sparameters_particles-%u.mat', SimulationFolder, filesep(), t), 'xparam');
+
str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)');
for l=1:size(xparam,1)
str = sprintf('%s\n %5.4f \t\t %7.5f \t\t %5.4f', str, lb95_xparam(l,t), mean_xparam(l,t), ub95_xparam(l,t));
end
disp(str)
disp('')
-end
-pmean = xparam(:,sample_size);
-pmode = mode_xparam(:,sample_size);
-pstdev = std_xparam(:,sample_size) ;
-p025 = lb95_xparam(:,sample_size) ;
-p975 = ub95_xparam(:,sample_size) ;
-pmedian = median_xparam(:,sample_size) ;
-covariance = mat_var_cov;
+end
%% Plot parameters trajectory
TeX = options_.TeX;
@@ -395,7 +385,7 @@ end
for plt = 1:nbplt
hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Trajectories');
- for k=1:length(pmean)
+ for k=1:number_of_parameters
subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
% Draw the surface for an interval containing 95% of the particles.
@@ -426,38 +416,3 @@ for plt = 1:nbplt
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
- hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
- for k=1:length(pmean)
- subplot(nr,nc,k)
- [name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
- optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function);
- [density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ...
- number_of_particles, optimal_bandwidth, kernel_function);
- plot(density(:,1), density(:,2));
- hold on
- if TeX
- title(texname,'interpreter','latex')
- else
- title(name,'interpreter','none')
- end
- hold off
- axis tight
- drawnow
- end
- dyn_saveas(hh_fig,[ M_.fname '_param_density' int2str(plt) ],options_.nodisplay,options_.graph_format);
- if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
- % TeX eps loader file
- fprintf(fidTeX, '\\begin{figure}[H]\n');
- fprintf(fidTeX,'\\centering \n');
- fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/nc,1),M_.fname,int2str(plt));
- fprintf(fidTeX,'\\caption{Parameter densities based on the Liu/West particle filter.}');
- fprintf(fidTeX,'\\label{Fig:ParameterDensities:%s}\n',int2str(plt));
- fprintf(fidTeX,'\\end{figure}\n');
- fprintf(fidTeX,' \n');
- end
-end
diff --git a/matlab/estimation/set_number_of_subdraws.m b/matlab/estimation/set_number_of_subdraws.m
index c6dcf6911d635903353e46e964907fdcc218f070..0241b32050e50376e00b3f4570f7c1cbb9ec13aa 100644
--- a/matlab/estimation/set_number_of_subdraws.m
+++ b/matlab/estimation/set_number_of_subdraws.m
@@ -57,6 +57,11 @@ else
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
NumberOfDraws = size(posterior.particles,2);
NumberOfDrawsPerChain=NumberOfDraws;
+ elseif isonline(options_)
+ % Load draws from the posterior distribution
+ posterior = load(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
+ NumberOfDraws = size(posterior.param,2);
+ NumberOfDrawsPerChain=NumberOfDraws;
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
diff --git a/matlab/issmc.m b/matlab/issmc.m
index 0f80730a55635be528e9fc6cf648921ff57b5957..45080a164049e1e9f671ae4dd8842e1d4a7b5934 100644
--- a/matlab/issmc.m
+++ b/matlab/issmc.m
@@ -17,4 +17,4 @@ function bool = issmc(options_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-bool = ishssmc(options_) || isdsmh(options_) || isdime(options_);
+bool = ishssmc(options_) || isonline(options_) || isdsmh(options_) || isdime(options_);
diff --git a/matlab/nonlinear-filters/solve_model_for_online_filter.m b/matlab/nonlinear-filters/solve_model_for_online_filter.m
index ca8d9fe08b305fca4640005be3ffea249384723d..6305f214e613ac8c5ee63a055e9bc9297ed12dfe 100644
--- a/matlab/nonlinear-filters/solve_model_for_online_filter.m
+++ b/matlab/nonlinear-filters/solve_model_for_online_filter.m
@@ -127,7 +127,6 @@ M_.H = H;
warning('off', 'MATLAB:nearlySingularMatrix')
[oo_.dr, info, M_.params] = ...
compute_decision_rules(M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
-% resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
warning('on', 'MATLAB:nearlySingularMatrix')
if info(1)~=0
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 990f9b146e08310cd2dc72077bb60a7800bcd32b..0452ca29559a56cbc2e7f68fbf4a4eadc0b51359 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -496,9 +496,11 @@ switch minimizer_algorithm
[LB, UB]=set_bounds_to_finite_values(bounds, options_.huge_number);
[opt_par_values, fval, exitflag] = simpsa(func2str(objective_function),start_par_value,LB,UB,simpsaOptions,varargin{:});
case 11
- options_.cova_compute = 0;
- subvarargin = [varargin(1), varargin(3:6), varargin(8)];
- opt_par_values = online_auxiliary_filter(start_par_value, subvarargin{:});
+ % waiting for validation
+ % options_.cova_compute = 0;
+ % subvarargin = [varargin(1), varargin(3:6), varargin(8)];
+ % opt_par_values = online_auxiliary_filter(start_par_value, subvarargin{:});
+ warning('Online particle filter is no more available with mode_compute=11; use posterior_sampler with option online')
case 12
if isoctave
error('Option mode_compute=12 is not available under Octave')
diff --git a/tests/particle/dsge_base2.mod b/tests/particle/dsge_base2.mod
index f07148e9f890408f5c946169220bcc8beed630ff..473e9366e750f8426650e30a4ae1544dad7ffa68 100644
--- a/tests/particle/dsge_base2.mod
+++ b/tests/particle/dsge_base2.mod
@@ -8,7 +8,7 @@
@#endif
@#ifndef ALGO_SIR
- @#define ALGO_SIR = 1
+ @#define ALGO_SIR = 0
@#endif
@#ifndef ALGO_APF
@@ -24,7 +24,7 @@
@#endif
@#ifndef ALGO_ONLINE
- @#define ALGO_ONLINE = 0
+ @#define ALGO_ONLINE = 1
@#endif
@#ifndef MCMC
@@ -150,9 +150,12 @@ estimation(order=3,nograph,filter_algorithm=gf,proposal_approximation=montecarlo
@#endif
@#if ALGO_ONLINE
- estimation(order=1,nograph,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
- estimation(order=2,nograph,number_of_particles=10000,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
- estimation(order=3,nograph,number_of_particles=10000,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
+% estimation(order=1,nograph,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
+% estimation(order=2,nograph,number_of_particles=10000,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
+% estimation(order=3,nograph,number_of_particles=10000,mode_compute=11,mh_replic=0,particle_filter_options=('liu_west_delta',0.9));
+ estimation(order=1,posterior_sampling_method='online',posterior_sampler_options=('particles',1000));
+ estimation(order=2,posterior_sampling_method='online',posterior_sampler_options=('particles',1000));
+ estimation(order=3,posterior_sampling_method='online',filter_algorithm=nlkf,proposal_approximation=montecarlo,number_of_particles=500,posterior_sampler_options=('particles',500));
@#endif
@#if MCMC
@@ -187,9 +190,7 @@ estimation(order=3,filter_algorithm=nlkf,number_of_particles=10000,proposal_appr
@#endif
@#if SMC
-% estimation(order=1,nograph,posterior_sampling_method='hssmc');
-% estimation(order=2,nograph,posterior_sampling_method='hssmc',filter_algorithm=nlkf,proposal_approximation=montecarlo);
-% estimation(order=2,nograph,posterior_sampling_method='hssmc',number_of_particles=10000);
-% estimation(order=3,nograph,posterior_sampling_method='hssmc',number_of_particles=10000,filter_algorithm=apf,resampling=none);%,posterior_sampler_options=('particles',1000));
- estimation(order=3,nograph,posterior_sampling_method='hssmc',filter_algorithm=nlkf,proposal_approximation=montecarlo);
+ estimation(order=1,posterior_sampling_method='hssmc',posterior_sampler_options=('particles',1000));
+ estimation(order=2,posterior_sampling_method='hssmc',posterior_sampler_options=('particles',1000));
+ estimation(order=3,posterior_sampling_method='hssmc',filter_algorithm=nlkf,proposal_approximation=montecarlo,number_of_particles=500,posterior_sampler_options=('particles',500));
@#endif