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,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;
@@ -657,6 +659,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;
 

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

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;

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;

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

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

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

matlab/estimation/dynare_estimation_1.m

@@ -185,7 +185,7 @@ end
 % Run smoother if no estimation or mode-finding are requested
 %
 
-if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation && ~issmc(options_)
+if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation && ~issmc(options_) 
     if options_.order==1 && ~options_.particle.status
         if options_.smoother
             if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
@@ -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)
@@ -271,14 +267,24 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
                         penalized_objective_function = str2func('penalty_objective_function');
                         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
+                        if options_.one_step.outer_product==1                            
+                            crit = options_.newrat.tolerance.f/10;
+                            hess_info.gstep = options_.gstep;
+                            hess_info.htol = crit;
+                            hess_info.htol_rel = options_.newrat.tolerance.gstep_rel;
+                            hess_info.h1 = options_.gradient_epsilon*ones(nx,1);
+                            hess_info.robust = options_.newrat.robust;
+                            hh = mr_hessian(xparam1,objective_function,fval,0,crit,hess_info,[bounds.lb bounds.ub],bayestopt_.p2,0,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)
+                elseif (isnumeric(current_optimizer) && isequal(current_optimizer,5)) 
                     % other numerical hessian options available with optimizer
                     % 5 and dsge_likelihood
                     %
-                    % if newratflag == 0
+                    % if newratflag == 01
                     % compute outer product gradient of optimizer 5
                     %
                     % if newratflag == 2
@@ -313,6 +319,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
     else
         save([M_.dname filesep 'Output' filesep M_.fname '_mode.mat'],'xparam1','parameter_names','fval');
     end
+    options_.one_step.hessian = hh ;
+
 end
 
 if ~options_.mh_posterior_mode_estimation && options_.cova_compute && ~issmc(options_)
@@ -337,6 +345,8 @@ oo_.posterior.optimization.log_density=[];
 
 invhess = [];
 
+%%
+
 if ~issmc(options_)
     if ~options_.mh_posterior_mode_estimation
         oo_.posterior.optimization.mode = xparam1;
@@ -397,9 +407,12 @@ 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 +485,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

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

matlab/estimation/online/isonline.m

+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

matlab/nonlinear-filters/online_auxiliary_filter.m → matlab/estimation/online/online_auxiliary_filter.m

-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

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;

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_);

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

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