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

@@ -405,7 +405,46 @@ if init
           options_.mh_replic = 0;
           options_.mh_posterior_mode_estimation = false;
           
-      case 'dime_mcmc'
+      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 '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(['dsmh: 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);
@@ -472,44 +511,44 @@ if init
           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;
+      % 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;
               
       otherwise
           error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);

matlab/estimation/dynare_estimation_1.m

@@ -394,20 +394,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
 
 %
@@ -471,7 +469,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/resampler/kitagawa.m

-function indices = kitagawa(weights, noise)
+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/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

tests/particle/DSMH.mod

+/*
+ * Copyright © 2004-2025 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/>.
+ *   
+ */
+
+
+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;
+
+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;
+
+varobs y l i ;
+
+// simulate the dataset
+set_dynare_seed(123456);
+stoch_simul(order=5,periods=1000,noprint,nograph);
+ds = dseries(oo_.endo_simul','1Y',M_.endo_names);
+data(series=ds,first_obs=801Y,nobs=200);
+
+estimation(order=1,posterior_sampling_method='dsmh',posterior_sampler_options=('particles',10000));

tests/particle/dsge_base2.mod

@@ -24,15 +24,19 @@
 @#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
+
+@#ifndef DSMH
+    @#define DSMH = 1
 @#endif
 
 %%
@@ -189,8 +193,14 @@ estimation(order=3,nograph,filter_algorithm=gf,proposal_approximation=montecarlo
 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
+@#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));
+%  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',10000));
+%  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