diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index d3773fc9377e6432e8e1689878908ca1e5743213..ff21c7cafad652352385a1cdebc83649a0cbf107 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -288,9 +288,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;
@@ -514,6 +511,10 @@ 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_.trace_plot_ma = 200;
options_.mh_autocorrelation_function_size = 30;
@@ -657,6 +658,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/estimation/GetAllPosteriorDraws.m b/matlab/estimation/GetAllPosteriorDraws.m
index 5c1847cf0283ecf87da59089dd78f843a7306ab4..5c9eb23265471e1c8651433ebe465ee0bcafa94a 100644
--- a/matlab/estimation/GetAllPosteriorDraws.m
+++ b/matlab/estimation/GetAllPosteriorDraws.m
@@ -33,6 +33,7 @@ function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstLine
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+<<<<<<< HEAD
if isnumeric(column) && column<0
error('GetAllPosteriorDraws:: column cannot be negative');
end
@@ -47,6 +48,48 @@ if issmc(options_)
draws = transpose(posterior.particles);
else
draws = transpose(posterior.particles(column,:));
+=======
+if ishssmc(options_)
+ % Load draws from the posterior distribution
+ pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
+ posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
+ if column==0
+ draws = posterior.tlogpostkernel;
+ else
+ draws = transpose(posterior.particles(column,:));
+ end
+elseif isdime(options_)
+ posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
+ tune = posterior.tune;
+ chains = posterior.chains(end-tune:end,:,:);
+ if column>0
+ chains = reshape(chains, [], size(chains, 3));
+ draws = chains(:,column);
+ else
+ draws = posterior.lprobs;
+ end
+ draws = transpose(posterior.particles(column,:));
+elseif isonline(options_)
+ posterior = load(sprintf('%s/online/parameters_particles_final.mat', dname));
+ draws = transpose(posterior.param(column,:));
+else
+ iline = FirstLine;
+ linee = 1;
+ DirectoryName = CheckPath('metropolis',dname);
+ if nblcks>1 && nargin<10
+ draws = zeros(NumberOfDraws*nblcks,1);
+ iline0=iline;
+ if column>0
+ for blck = 1:nblcks
+ iline=iline0;
+ for file = FirstMhFile:TotalNumberOfMhFile
+ load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
+ NumberOfLines = size(x2(iline:end,:),1);
+ draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
+ linee = linee+NumberOfLines;
+ iline = 1;
+ end
+>>>>>>> 377185759 (Make compatible the outputs of the online particle filter with the Dynare routines for presenting results + migrate the Online Particle Filter from nonlinear filters to posterior samplers. Now require to create an online option for posterior_sampling_method.)
end
else
draws = posterior.tlogpostkernel;
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..63a40adcb75b3b76976f59804a213181848470cd 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -55,7 +55,7 @@ skipline(2)
disp('ESTIMATION RESULTS')
skipline()
-if ishssmc(options_)
+if ishssmc(options_) || isonline(options_)
dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
elseif isdime(options_)
@@ -87,7 +87,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..12d570057df990f62498a7ef23f4e07e4bdb95a3 100644
--- a/matlab/estimation/check_posterior_sampler_options.m
+++ b/matlab/estimation/check_posterior_sampler_options.m
@@ -374,125 +374,148 @@ 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};
+ 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 a2e9ffd8c100a8ce8a4180369dde10be3eeba310..a554146b21090afe4efecf83081711800e25b29a 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -238,11 +238,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)
@@ -272,7 +268,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
@@ -399,7 +395,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_);
@@ -472,7 +472,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/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 81%
rename from matlab/nonlinear-filters/online_auxiliary_filter.m
rename to matlab/estimation/online/online_auxiliary_filter.m
index e392474ea5681c7e6bc319aec860cfb0ef88d87c..d1fce18e659e0fc459dd1c8e24b111a0b536a8a3 100644
--- a/matlab/nonlinear-filters/online_auxiliary_filter.m
+++ b/matlab/estimation/online/online_auxiliary_filter.m
@@ -1,4 +1,5 @@
-function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
+%function [xparam,mean_xparam,median_xparam,lb95_xparam,ub95_xparam] =
+function 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_)
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
%
@@ -13,14 +14,12 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
% - 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.
-
+% - xparam [double] n×M matrix, parameters particles at the end of the sample (posterior approx.).
+% - mean_xparam [double] n×T matrix, time series of the parameters particles mean.
+% - median_xparam [double] n×T matrix, time series of the parameters particles median.
+% - lb95_xparam [double] n×T matrix, time series of 2.5 percent of the particles are below p025(i) for i=1,…,n.
+% - ub95_xparam [double] n×T matrix, time series of 97.5 percent of the particles are below p975(i) for i=1,…,n.
+%
% Copyright © 2013-2023 Dynare Team
%
% This file is part of Dynare.
@@ -43,6 +42,10 @@ 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 +55,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 +102,7 @@ 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);
+%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 +197,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 +204,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 +216,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 +300,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
@@ -314,9 +316,9 @@ for t=1:sample_size
end
% final resampling (not advised)
if second_resample
- [~, idmode] = max(weights);
- mode_xparam(:,t) = xparam(:,idmode);
- indx = resample(0, weights,options_.particle);
+% [~, idmode] = max(weights);
+% mode_xparam(:,t) = xparam(:,idmode);
+ indx = kitagawa(weights);
StateVectors = StateVectors(:,indx) ;
if pruning
StateVectors_ = StateVectors_(:,indx);
@@ -333,10 +335,13 @@ 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
+% [~, idmode] = max(weights);
+% mode_xparam(:,t) = xparam(:,idmode);
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 +367,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 +399,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 +430,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..fd3bdec202e3a574318004eef8354e04fb273f52 100644
--- a/matlab/estimation/set_number_of_subdraws.m
+++ b/matlab/estimation/set_number_of_subdraws.m
@@ -57,6 +57,12 @@ 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
+ pfiles = dir(sprintf('%s/online/particles-*.mat', M_.dname));
+ posterior = load(sprintf('%s/online/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
+ 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