matlab/estimation/GetAllPosteriorDraws.m

@@ -51,6 +51,18 @@ if issmc(options_)
         else
             draws = posterior.tlogpostkernel;
         end
+    elseif isdsmh(options_)
+        % Load draws from the posterior distribution
+        posterior = load(sprintf('%s/dsmh/parameters_particles_final.mat', dname));
+        if column>0 || strcmp(column,'all')
+            if strcmp(column,'all')
+                draws = transpose(posterior.particles);
+            else
+                draws = transpose(posterior.particles(column,:));
+            end
+        else
+            draws = posterior.tlogpostkernel;
+        end
     elseif isonline(options_)
         % Load draws from the posterior distribution
         posterior = load(sprintf('%s/online/parameters_particles_final.mat', dname));

matlab/estimation/GetPosteriorMeanVariance.m

@@ -36,6 +36,13 @@ 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 isdsmh(options_)
+        % Load draws from the posterior distribution
+        posterior = load(sprintf('%s/dsmh/parameters_particles_final.mat', M_.dname));
+        % Compute the posterior mean
+        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

matlab/estimation/GetPosteriorParametersStatistics.m

@@ -58,6 +58,9 @@ skipline()
 if ishssmc(options_) 
     dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
     hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
+elseif isdsmh(options_) 
+    dprintf('Log data density is %f.', oo_.MarginalDensity.dsmh);
+    hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
 elseif isonline(options_)
     hpd_draws = num_draws;
 elseif isdime(options_)
@@ -89,7 +92,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_) || isonline(options_)
+        if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_) || isdsmh(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 isdsmh(options_)
+        [MetropolisFolder, info] = CheckPath('dsmh',M_.dname);
     elseif isonline(options_)
         [MetropolisFolder, info] = CheckPath('online',M_.dname);
     elseif isdime(options_)
@@ -78,6 +80,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 isdsmh(options_)
+            % Load draws from the posterior distribution
+            pfiles = dir(sprintf('%s/dsmh/particles-*.mat', M_.dname));
+            mhdate = get_date_of_a_file(sprintf('%s/dsmh/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));

matlab/estimation/check_posterior_sampler_options.m

@@ -482,34 +482,37 @@ if init
               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'
+                      case 'H'
+                          posterior_sampler_options.H = options_list{i,2};
+                      case 'N'
+                          posterior_sampler_options.N = options_list{i,2};
+                      case 'G'
+                          posterior_sampler_options.G = options_list{i,2};
+                      case 'K'
+                          posterior_sampler_options.K = options_list{i,2};
+                      case 'lambda1'
+                          posterior_sampler_options.lambda1 = options_list{i,2};
+                      case 'tau'
+                          posterior_sampler_options.tau = options_list{i,2};
+                      case 'particles'
                           posterior_sampler_options.particles = options_list{i,2};
+                      case 'alpha0'
+                          posterior_sampler_options.alpha0 = options_list{i,2};
+                      case 'alpha1'
+                          posterior_sampler_options.alpha1 = options_list{i,2};
                       otherwise
-                          warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
-                  end
+                          warning(['dsmh: Unknown option (' options_list{i,1} ')!'])
                   end
+              end
+          end
+
+          if posterior_sampler_options.particles<posterior_sampler_options.N*posterior_sampler_options.G
+              error('check_posterior_sampler_options:: DSMH requires particles to be at least than N*G = %u ',posterior_sampler_options.N*posterior_sampler_options.G);             
           end
-          
           options_.mode_compute = 0;
           options_.cova_compute = 0;
           options_.mh_replic = 0;
-          options_.mh_posterior_mode_estimation = true;
+          options_.mh_posterior_mode_estimation = false;
               
       otherwise
           error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);

matlab/estimation/dynare_estimation_1.m

@@ -405,20 +405,18 @@ end
 % Run SMC sampler.
 %
 
-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_);
-elseif isdime(options_)
+if issmc(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;
-    dime(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
-elseif isdsmh(options_)
-    dsmh(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
+    if isonline(options_)
+        online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_);
+    elseif ishssmc(options_)
+        oo_.MarginalDensity.hssmc = hssmc(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+    elseif isdime(options_)
+        dime(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+    elseif isdsmh(options_)
+        oo_.MarginalDensity.dsmh = dsmh(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+    end 
 end
 
 %
@@ -482,7 +480,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_) || isonline(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
+        if isdsmh(options_) || 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

@@ -34,6 +34,8 @@ if ~issmc(options_)
 else
     if ishssmc(options_)
         folder_name='hssmc';
+    elseif isdsmh(options_)
+        folder_name='dsmh';
     elseif isdime(options_)
         folder_name='dime';
     elseif isonline(options_)

matlab/estimation/online/online_auxiliary_filter.m

@@ -5,7 +5,6 @@ function online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_,
 % INPUTS
 % - xparam1                  [double]    n×1 vector, Initial condition for the estimated parameters.
 % - dataset_                 [dseries]   Sample used for estimation.
-% - dataset_info             [struct]    Description of the sample.
 % - options_                 [struct]    Option values.
 % - M_                       [struct]    Description of the model.
 % - estim_params_            [struct]    Description of the estimated parameters.
@@ -373,7 +372,7 @@ nc = floor(sqrt(number_of_parameters));
 nbplt = 1 ;
 
 if TeX
-    fidTeX = fopen([M_.fname '_param_traj.tex'],'w');
+    fidTeX = fopen([M_.dname filesep 'online' filesep M_.fname '_param_traj.tex'],'w');
     fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
     fprintf(fidTeX,' \n');
@@ -399,12 +398,12 @@ for plt = 1:nbplt
         axis tight
         drawnow
     end
-    dyn_saveas(hh_fig, [M_.fname '_param_traj' int2str(plt)], options_.nodisplay, options_.graph_format);
+    dyn_saveas(hh_fig, [M_.dname filesep 'online' filesep M_.fname '_param_traj' int2str(plt)], options_.nodisplay, options_.graph_format);
     if TeX
         % TeX eps loader file
         fprintf(fidTeX,'\\begin{figure}[H]\n');
         fprintf(fidTeX,'\\centering \n');
-        fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',M_.fname,int2str(plt));
+        fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s/online/%s_ParamTraj%s}\n',M_.dname,M_.fname,int2str(plt));
         fprintf(fidTeX,'\\caption{Parameters trajectories.}');
         fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt));
         fprintf(fidTeX,'\\end{figure}\n');

matlab/estimation/resampler/kitagawa.m

-function indices = kitagawa(weights, noise)
-
+function indices = kitagawa(weights, noise, m)
+% function indices = kitagawa(weights, noise, m)
 % Return indices for resampling.
 %
 % INPUTS
-% - weights   [double]    n×1 vector of partcles' weights.
+% - weights   [double]    n×1 vector of particles' weights.
 % - noise     [double]    scalar, uniform random deviates in [0,1]
+% - m         [integer]   scalar, number of particles to resample 
 %
 % OUTPUTS
-% - indices   [integer]   n×1 vector of indices in [1:n]
+% - indices   [integer]   m×1 vector of indices in [1:n]
 
-% Copyright © 2022-2023 Dynare Team
+% Copyright © 2022-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -26,7 +27,11 @@ function indices = kitagawa(weights, noise)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-n= length(weights);
+if nargin<3
+    n = length(weights);
+else 
+    n = m ;
+end 
 
 if nargin<2, noise = rand; 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 isdsmh(options_)
+        % Load draws from the posterior distribution
+        posterior = load(sprintf('%s/dsmh/parameters_particles_final.mat', M_.dname));
+        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));

matlab/estimation/smc/dsmh.m

-function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
-
+function mdd = dsmh(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
+% function mdd = dsmh(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
 % Dynamic Striated Metropolis-Hastings algorithm.
+% based on Waggoner/Wu/Zha (2016): "Striated Metropolis–Hastings sampler for high-dimensional models,"
+% Journal of Econometrics, 192(2): 406-420, https://doi.org/10.1016/j.jeconom.2016.02.007
 %
 % INPUTS
 %   o TargetFun  [char]     string specifying the name of the objective
 %                           function (posterior kernel).
-%   o xparam1    [double]   (p*1) vector of parameters to be estimated (initial values).
 %   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters.
 %   o dataset_              data structure
 %   o dataset_info          dataset info structure
@@ -17,23 +18,8 @@ function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M
 %
 % SPECIAL REQUIREMENTS
 %   None.
-%
-% PARALLEL CONTEXT
-% The most computationally intensive part of this function may be executed
-% in parallel. The code suitable to be executed in
-% parallel on multi core or cluster machine (in general a 'for' cycle)
-% has been removed from this function and been placed in the posterior_sampler_core.m funtion.
-%
-% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
-% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
-% to manage the parallel computations.
-%
-% This function was the first function to be parallelized. Later, other
-% functions have been parallelized using the same methodology.
-% Then the comments write here can be used for all the other pairs of
-% parallel functions and also for management functions.
 
-% Copyright © 2022-2023 Dynare Team
+% Copyright © 2022-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -50,123 +36,71 @@ function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-opts = options_.posterior_sampler_options.dsmh;
+opts = options_.posterior_sampler_options.current_options;
+
+% Set location for the simulated particles.
+SimulationFolder = CheckPath('dsmh', M_.dname);
+%delete old stale files before creating new ones
+delete_stale_file(sprintf('%s%sparticles-*.mat', SimulationFolder,filesep))
+
+% Define prior distribution
+Prior = dprior(bayestopt_, options_.prior_trunc);
 
-lambda = exp(bsxfun(@minus,options_.posterior_sampler_options.dsmh.H,1:1:options_.posterior_sampler_options.dsmh.H)/(options_.posterior_sampler_options.dsmh.H-1)*log(options_.posterior_sampler_options.dsmh.lambda1));
+% Set function handle for the objective
+eval(sprintf('%s = @(x) %s(x, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, mh_bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);', 'funobj', func2str(TargetFun)));
+
+lambda = exp(bsxfun(@minus,opts.H,1:1:opts.H)/(opts.H-1)*log(opts.lambda1));
 c = 0.055 ;
-MM = int64(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/10) ;
+MM = int64(opts.N*opts.G/10) ;
 
 % Step 0: Initialization of the sampler
-[param, tlogpost_iminus1, loglik, bayestopt_] = ...
-    smc_samplers_initialization(TargetFun, 'dsmh', opts.particles, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+[param, tlogpost_iminus1, loglik] = ...
+    smc_samplers_initialization(funobj, 'dsmh', opts.particles, Prior, SimulationFolder, opts.H) ;
 
-ESS = zeros(options_.posterior_sampler_options.dsmh.H,1) ;
-zhat = 1 ;
+ESS = zeros(opts.H,1) ;
+zhat = 0 ;
 
 % The DSMH starts here
-for i=2:options_.posterior_sampler_options.dsmh.H
-    disp('');
-    disp('Tempered iteration');
-    disp(i) ;
+dprintf('#Iter.       lambda         ESS                 c     Accept. rate       scale      resample    seconds')
+for i=2:opts.H
+    t0 = tic;
     % Step 1: sort the densities and compute IS weigths
-    [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param) ;
-    [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS) ;
+    [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1, loglik,param) ;
+    [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param, tlogpost_iminus1, loglik, lambda, i, zhat, ESS) ;
     % Step 2: tune c_i
-    c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+    [c,acpt] = tune_c(funobj, param, tlogpost_i, lambda, i, c, Omegachol, weights, mh_bounds, opts, Prior) ;
     % Step 3: Metropolis step
-    [param,tlogpost_iminus1,loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+    [param,tlogpost_iminus1,loglik] = mutation_DSMH(funobj, param, tlogpost_i, tlogpost_iminus1, loglik, lambda, i, c, MM, Omegachol, weights, mh_bounds, opts, Prior) ;
+    tt = toc(t0) ;
+    dprintf('%3u          %5.4f     %9.5E         %5.4f        %5.4f        %+5.4f        %3s       %5.2f', i, lambda(i), ESS(i), c, acpt, zhat, 'no', tt)
+    %save(sprintf('%s%sparticles-%u-%u.mat', SimulationFolder, filesep(), i, opts.H), 'param', 'tlogpost', 'loglik')
 end
 
+tlogpost = tlogpost_iminus1 + loglik*(lambda(end)-lambda(end-1));
 weights = exp(loglik*(lambda(end)-lambda(end-1)));
 weights = weights/sum(weights);
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.particles);
-distrib_param = param(:,indx_resmpl);
-
-mean_xparam = mean(distrib_param,2);
-npar  = length(xparam1);
-lb95_xparam = zeros(npar,1) ;
-ub95_xparam = zeros(npar,1) ;
-for i=1:npar
-    temp = sortrows(distrib_param(i,:)') ;
-    lb95_xparam(i) = temp(0.025*options_.posterior_sampler_options.dsmh.particles) ;
-    ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.particles) ;
-end
+iresample = kitagawa(weights);
+particles = param(:,iresample);
+tlogpostkernel = tlogpost(iresample);
+loglikelihood = loglik(iresample);
+save(sprintf('%s%sparameters_particles_final.mat', SimulationFolder, filesep()), 'particles', 'tlogpostkernel', 'loglikelihood')
 
-TeX = options_.TeX;
-
-str = sprintf(' Param. \t Lower Bound (95%%) \t Mean \t Upper Bound (95%%)');
-for l=1:npar
-    name = get_the_name(l,TeX,M_,estim_params_,options_.varobs);
-    str = sprintf('%s\n %s \t\t %5.4f \t\t %7.5f \t\t %5.4f', str, name, lb95_xparam(l), mean_xparam(l), ub95_xparam(l));
-end
-disp(str)
-disp('')
-
-%% Plot parameters densities
-
-if TeX
-    fidTeX = fopen([M_.fname '_param_density.tex'],'w');
-    fprintf(fidTeX,'%% TeX eps-loader file generated by DSMH.m (Dynare).\n');
-    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
-    fprintf(fidTeX,' \n');
-end
-
-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.
-
-plt = 1 ;
-%for plt = 1:nbplt,
-if TeX
-    NAMES = [];
-    TeXNAMES = [];
-end
-hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
-for k=1:npar %min(nstar,npar-(plt-1)*nstar)
-    subplot(ceil(sqrt(npar)),floor(sqrt(npar)),k)
-    %kk = (plt-1)*nstar+k;
-    [name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
-    optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.dsmh.particles,bandwidth,kernel_function);
-    [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
-        options_.posterior_sampler_options.dsmh.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/floor(sqrt(npar)),1),M_.fname,int2str(plt));
-    fprintf(fidTeX,'\\caption{Parameter densities based on the Dynamic Striated Metropolis-Hastings algorithm.}');
-    fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
-    fprintf(fidTeX,'\\end{figure}\n');
-    fprintf(fidTeX,' \n');
-end
-%end
+mdd = zhat;
 
-function [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param)
+function [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1, loglik,param)
 [~,indx_ord] = sortrows(tlogpost_iminus1);
 tlogpost_iminus1 = tlogpost_iminus1(indx_ord);
 param = param(:,indx_ord);
 loglik = loglik(indx_ord);
 
-function [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS)
+function [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param, tlogpost_iminus1, loglik, lambda, i, zhat, ESS)
 if i==1
     tlogpost_i = tlogpost_iminus1 + loglik*lambda(i);
 else
     tlogpost_i = tlogpost_iminus1 + loglik*(lambda(i)-lambda(i-1));
 end
 weights = exp(tlogpost_i-tlogpost_iminus1);
-zhat = (mean(weights))*zhat ;
+zhat = zhat + log(mean(weights)) ;
 weights = weights/sum(weights);
 ESS(i) = 1/sum(weights.^2);
 % estimates of mean and variance
@@ -175,30 +109,30 @@ z = bsxfun(@minus,param,mu);
 Omega = z*diag(weights)*z';
 Omegachol = chol(Omega)';
 
-function c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-disp('tuning c_i...');
-disp('Initial value =');
-disp(c) ;
+function [c,acpt] = tune_c(funobj, param, tlogpost_i, lambda, i, c, Omegachol, weights, mh_bounds, opts, Prior)
+%        disp('tuning c_i...');
+%        disp('Initial value =');
+%        disp(c) ;
 npar = size(param,1);
-lower_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^5;
-upper_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^(1/5);
+lower_prob = (.5*(opts.alpha0+opts.alpha1))^5;
+upper_prob = (.5*(opts.alpha0+opts.alpha1))^(1/5);
 stop=0 ;
 while stop==0
     acpt = 0.0;
-    indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
+    indx_resmpl = kitagawa(weights,rand(1,1),opts.G);
     param0 = param(:,indx_resmpl);
     tlogpost0 = tlogpost_i(indx_resmpl);
-    for j=1:options_.posterior_sampler_options.dsmh.G
-        for l=1:options_.posterior_sampler_options.dsmh.K
+    for j=1:opts.G
+        for l=1:opts.K
             validate = 0;
             while validate == 0
                 candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
                 if all(candidate >= mh_bounds.lb) && all(candidate <= mh_bounds.ub)
-                    [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+                    [tlogpostx,loglikx] = tempered_likelihood(funobj, candidate, lambda(i), Prior);
                     if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
                         validate = 1;
                         if rand(1,1)<exp(tlogpostx-tlogpost0(j)) % accept
-                            acpt = acpt + 1/(options_.posterior_sampler_options.dsmh.G*options_.posterior_sampler_options.dsmh.K);
+                            acpt = acpt + 1/(opts.G*opts.K);
                             param0(:,j)= candidate;
                             tlogpost0(j) = tlogpostx;
                         end
@@ -207,29 +141,29 @@ while stop==0
             end
         end
     end
-    disp('Acceptation rate =') ;
-    disp(acpt) ;
-    if options_.posterior_sampler_options.dsmh.alpha0<=acpt && acpt<=options_.posterior_sampler_options.dsmh.alpha1
-        disp('done!');
+    %           disp('Acceptation rate =') ;
+    %           disp(acpt) ;
+    if opts.alpha0<=acpt && acpt<=opts.alpha1
+        %                disp('done!');
         stop=1;
     else
         if acpt<lower_prob
             c = c/5;
         elseif lower_prob<=acpt && acpt<=upper_prob
-            c = c*log(.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))/log(acpt);
+            c = c*log(.5*(opts.alpha0+opts.alpha1))/log(acpt);
         else
             c = 5*c;
         end
-        disp('Trying with c= ') ;
-        disp(c)
+        %                disp('Trying with c= ') ;
+        %                disp(c)
     end
 end
 
-function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-indx_levels = (1:1:MM-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM;
+function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(funobj, param, tlogpost_i, tlogpost_iminus1, loglik, lambda, i, c, MM, Omegachol, weights, mh_bounds, opts, Prior)
+indx_levels = (1:1:MM-1)*opts.N*opts.G/MM;
 npar = size(param,1) ;
-p = 1/(10*options_.posterior_sampler_options.dsmh.tau);
-disp('Metropolis step...');
+p = 1/(10*opts.tau);
+%        disp('Metropolis step...');
 % build the dynamic grid of levels
 levels = [0.0;tlogpost_iminus1(indx_levels)];
 % initialize the outputs
@@ -237,14 +171,14 @@ out_param = param;
 out_tlogpost_iminus1 = tlogpost_i;
 out_loglik = loglik;
 % resample and initialize the starting groups
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
+indx_resmpl = kitagawa(weights,rand(1,1),opts.G);
 param0 = param(:,indx_resmpl);
 tlogpost_iminus10 = tlogpost_iminus1(indx_resmpl);
 tlogpost_i0 = tlogpost_i(indx_resmpl);
 loglik0 = loglik(indx_resmpl);
 % Start the Metropolis
-for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.tau
-    for j=1:options_.posterior_sampler_options.dsmh.G
+for l=1:opts.N*opts.tau
+    for j=1:opts.G
         u1 = rand(1,1);
         u2 = rand(1,1);
         if u1<p
@@ -255,7 +189,7 @@ for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_opt
                     break
                 end
             end
-            indx = floor( (k-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM+1 + u2*(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM-1) );
+            indx = floor( (k-1)*opts.N*opts.G/MM+1 + u2*(opts.N*opts.G/MM-1) );
             if i==1
                 alp = (loglik(indx)-loglik0(j))*lambda(i);
             else
@@ -272,7 +206,7 @@ for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_opt
             while validate==0
                 candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
                 if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
-                    [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
+                    [tlogpostx, loglikx] = tempered_likelihood(funobj, candidate, lambda(i), Prior);
                     if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
                         validate = 1;
                         if u2<exp(tlogpostx-tlogpost_i0(j)) % accept
@@ -290,8 +224,8 @@ for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_opt
             end
         end
     end
-    if mod(l,options_.posterior_sampler_options.dsmh.tau)==0
-        rang = (l/options_.posterior_sampler_options.dsmh.tau-1)*options_.posterior_sampler_options.dsmh.G+1:l*options_.posterior_sampler_options.dsmh.G/options_.posterior_sampler_options.dsmh.tau;
+    if mod(l,opts.tau)==0
+        rang = (l/opts.tau-1)*opts.G+1:l*opts.G/opts.tau;
         out_param(:,rang) = param0;
         out_tlogpost_iminus1(rang) = tlogpost_i0;
         out_loglik(rang) = loglik0;

matlab/estimation/smc/hssmc.m

@@ -4,7 +4,6 @@ function mdd = hssmc(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_,
 %
 % INPUTS
 % - TargetFun        [char]     string specifying the name of the objective function (posterior kernel).
-% - xparam1          [double]   p×1 vector of parameters to be estimated (initial values).
 % - mh_bounds        [double]   p×2 matrix defining lower and upper bounds for the parameters.
 % - dataset_         [dseries]  sample
 % - dataset_info     [struct]   informations about the dataset

meson.build

@@ -1768,11 +1768,14 @@ mod_and_m_tests = [
   { 'test' : [ 'ecb/pooled_fgls/test_param_names.mod' ] },
 
   # Particle files
-  { 'test' : [ 'particle/dsge_base2.mod' ],
+  { 'test' : [ 'particle/dsge_base2.mod',
+               'particle/dsge_base_Liu_West.mod',
+               'particle/dsge_base_HSSMC_DSMH.mod',],
     'extra' : [ 'particle/risky.m',
                 'particle/extreme.m',
                 'particle/benchmark.m',
-                'particle/mysample.m'] },
+                'particle/mysample.m',
+                'particle/dsge_base.inc'] },
   { 'test' : [ 'particle/first_spec.mod',
                'particle/first_spec_hssmc.mod',
                'particle/local_state_space_iteration_k_test.mod',

tests/particle/dsge_base.inc

+var k A c l i y;
+varexo e_a;
+
+parameters alp bet tet tau delt rho ;
+alp = 0.4;
+bet = 0.99;
+tet = 0.357 ;
+tau =  50 ;
+delt = 0.02;
+rho = 0.95;
+
+model;
+	c = ((1 - alp)*tet/(1-tet))*A*(1-l)*((k(-1)/l)^alp) ;
+	y = A*(k(-1)^alp)*(l^(1-alp)) ;
+	i = y-c ;
+	k = (1-delt)*k(-1) + i ;
+	log(A) = rho*log(A(-1)) + e_a ;
+	(((c^(tet))*((1-l)^(1-tet)))^(1-tau))/c - bet*((((c(+1)^(tet))*((1-l(+1))^(1-tet)))^(1-tau))/c(+1))*(1 -delt+alp*(A(1)*(k^alp)*(l(1)^(1-alp)))/k)=0 ;
+end;
+
+steady_state_model;
+	k = -(alp-1)*(alp^(1/(1-alp)))*(bet^(1/(1-alp)))*((bet*(delt-1)+1)^(alp/(alp-1)))*tet/(-alp*delt*bet+delt*bet+alp*tet*bet-bet-alp*tet+1);
+	l = (alp-1)*(bet*(delt-1)+1)*tet/(alp*tet+bet*((alp-1)*delt-alp*tet+1)-1) ;
+	y = (k^alp)*(l^(1-alp)) ;
+	i = delt*k ;
+	c = y - i ;
+	A = 1;
+end;
+
+shocks;
+    var e_a; stderr 0.035;
+end;
+
+steady;
+
+estimated_params;
+    alp, uniform_pdf,,, 0.0001, 0.99;
+    bet, uniform_pdf,,, 0.0001, 0.99999;
+    tet, uniform_pdf,,, 0.0001, .999;
+    tau, uniform_pdf,,, 0.0001, 100;
+    delt, uniform_pdf,,, 0.0001, 0.05;
+    rho, uniform_pdf,,, 0.0001, 0.9999;
+    stderr e_a, uniform_pdf,,, 0.00001, 0.1;
+    stderr y, uniform_pdf,,, 0.00001, 0.1;
+    stderr l, uniform_pdf,,, 0.00001, 0.1;
+    stderr i, uniform_pdf,,, 0.00001, 0.1;
+end;
+
+varobs y l i ;
+
+%data(file='./risky.m');
+%stoch_simul(order=5,periods=1000);
+%datatomfile('mysample')
+%return;
+
+data(file='./mysample.m',first_obs=801Y,nobs=50); %no measurement errors added in the simulated data

tests/particle/dsge_base2.mod

@@ -8,7 +8,7 @@
 @#endif
 
 @#ifndef ALGO_SIR
-    @#define ALGO_SIR = 0
+    @#define ALGO_SIR = 1
 @#endif
 
 @#ifndef ALGO_APF
@@ -24,66 +24,25 @@
 @#endif
 
 @#ifndef ALGO_ONLINE
-    @#define ALGO_ONLINE = 1
+    @#define ALGO_ONLINE = 0
 @#endif
 
 @#ifndef MCMC
     @#define MCMC = 0
 @#endif
 
-@#ifndef SMC
-    @#define SMC = 0
+@#ifndef HSSMC
+    @#define HSSMC = 0
 @#endif
 
-%%
-
-var k A c l i y;
-varexo e_a;
-
-parameters alp bet tet tau delt rho ;
-alp = 0.4;
-bet = 0.99;
-tet = 0.357 ;
-tau =  50 ;
-delt = 0.02;
-rho = 0.95;
-
-model;
-	c = ((1 - alp)*tet/(1-tet))*A*(1-l)*((k(-1)/l)^alp) ;
-	y = A*(k(-1)^alp)*(l^(1-alp)) ;
-	i = y-c ;
-	k = (1-delt)*k(-1) + i ;
-	log(A) = rho*log(A(-1)) + e_a ;
-	(((c^(tet))*((1-l)^(1-tet)))^(1-tau))/c - bet*((((c(+1)^(tet))*((1-l(+1))^(1-tet)))^(1-tau))/c(+1))*(1 -delt+alp*(A(1)*(k^alp)*(l(1)^(1-alp)))/k)=0 ;
-end;
+@#ifndef DSMH
+    @#define DSMH = 0
+@#endif
 
-steady_state_model;
-	k = -(alp-1)*(alp^(1/(1-alp)))*(bet^(1/(1-alp)))*((bet*(delt-1)+1)^(alp/(alp-1)))*tet/(-alp*delt*bet+delt*bet+alp*tet*bet-bet-alp*tet+1);
-	l = (alp-1)*(bet*(delt-1)+1)*tet/(alp*tet+bet*((alp-1)*delt-alp*tet+1)-1) ;
-	y = (k^alp)*(l^(1-alp)) ;
-	i = delt*k ;
-	c = y - i ;
-	A = 1;
-end;
+%%
 
-shocks;
-    var e_a; stderr 0.035;
-end;
+@#include "dsge_base.inc"
 
-steady;
-
-estimated_params;
-    alp, uniform_pdf,,, 0.0001, 0.99;
-    bet, uniform_pdf,,, 0.0001, 0.99999;
-    tet, uniform_pdf,,, 0.0001, .999;
-    tau, uniform_pdf,,, 0.0001, 100;
-    delt, uniform_pdf,,, 0.0001, 0.05;
-    rho, uniform_pdf,,, 0.0001, 0.9999;
-    stderr e_a, uniform_pdf,,, 0.00001, 0.1;
-    stderr y, uniform_pdf,,, 0.00001, 0.1;
-    stderr l, uniform_pdf,,, 0.00001, 0.1;
-    stderr i, uniform_pdf,,, 0.00001, 0.1;
-end;
 
 estimated_params_init;
 	alp, 0.4;
@@ -98,14 +57,6 @@ estimated_params_init;
 	stderr i, .001;
 end;
 
-varobs y l i ;
-
-%data(file='./risky.m');
-%stoch_simul(order=5,periods=1000);
-%datatomfile('mysample')
-%return;
-
-data(file='./mysample.m',first_obs=801Y,nobs=50); %no measurement errors added in the simulated data
 
 @#if LINEAR_KALMAN
 	estimation(nograph,order=1,mode_compute=8,silent_optimizer,mh_replic=0,additional_optimizer_steps=[8 4],mode_check);
@@ -188,9 +139,3 @@ estimation(order=3,nograph,filter_algorithm=gf,proposal_approximation=montecarlo
   end;  
 estimation(order=3,filter_algorithm=nlkf,number_of_particles=10000,proposal_approximation=montecarlo,resampling=none,silent_optimizer,mode_compute=0,cova_compute=0,MCMC_jumping_covariance=prior_variance,mh_init_scale_factor=0.01);
 @#endif
-
-@#if SMC
-  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

tests/particle/dsge_base_HSSMC_DSMH.mod

+@#ifndef HSSMC
+    @#define HSSMC = 1
+@#endif
+
+@#ifndef DSMH
+    @#define DSMH = 1
+@#endif
+
+@#include "dsge_base.inc"
+
+    
+estimated_params_init;
+	alp, 0.4;
+	bet, 0.99;
+	tet, 0.357;
+	tau, 50;
+	delt, 0.02;
+	rho, 0.95;
+	stderr e_a, .035;
+	stderr y, .001;
+	stderr l, .001;
+	stderr i, .001;
+end;
+
+@#if HSSMC
+  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
+
+@#if DSMH
+  estimation(order=1,posterior_sampling_method='dsmh',posterior_sampler_options=('particles',200));
+%  estimation(order=2,posterior_sampling_method='dsmh',posterior_sampler_options=('particles',1000));
+%  estimation(order=3,posterior_sampling_method='dsmh',filter_algorithm=nlkf,proposal_approximation=montecarlo,number_of_particles=500,posterior_sampler_options=('particles',500));
+@#endif
\ No newline at end of file

tests/particle/dsge_base_Liu_West.mod

+@#ifndef ALGO_ONLINE
+    @#define ALGO_ONLINE = 1
+@#endif
+
+@#include "dsge_base.inc"
+
+estimated_params_init;
+	alp, 0.4;
+	bet, 0.99;
+	tet, 0.357;
+	tau, 50;
+	delt, 0.02;
+	rho, 0.95;
+	stderr e_a, .035;
+	stderr y, .001;
+	stderr l, .001;
+	stderr i, .001;
+end;
+options_.TeX=true;
+
+@#if ALGO_ONLINE
+  // estimation(order=1,posterior_sampling_method='online',posterior_sampler_options=('particles',100));
+  // estimation(order=2,posterior_sampling_method='online',posterior_sampler_options=('particles',100));
+  estimation(order=3,posterior_sampling_method='online',filter_algorithm=nlkf,proposal_approximation=montecarlo,posterior_sampler_options=('particles',100));
+@#endif
+
+collect_latex_files;
+[status, cmdout]=system(['pdflatex -halt-on-error -interaction=nonstopmode ' M_.fname '_TeX_binder.tex']);
+if status
+    cmdout
+    error('TeX-File did not compile.')
+end