diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 6d6919166d40b8ec9a69c727b90545e1e4820f29..9288951466ddf311b10dd224e80fe014cf7a1fa4 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -515,6 +515,7 @@ options_.posterior_sampler_options.dsmh.tau = 10 ;
 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.systematic_resampling = false;
 
 options_.trace_plot_ma = 200;
 options_.mh_autocorrelation_function_size = 30;
diff --git a/matlab/estimation/check_posterior_sampler_options.m b/matlab/estimation/check_posterior_sampler_options.m
index 12d570057df990f62498a7ef23f4e07e4bdb95a3..ccf43fed140479358fa53dcc5b49ad21091b2a9a 100644
--- a/matlab/estimation/check_posterior_sampler_options.m
+++ b/matlab/estimation/check_posterior_sampler_options.m
@@ -458,7 +458,9 @@ if init
                       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};
+                          posterior_sampler_options.liu_west_max_resampling_tries = options_list{i,2};
+                      case 'systematic_resampling'
+                          posterior_sampler_options.systematic_resampling = options_list{i,2};
                       otherwise
                           warning(['online: Unknown option (' options_list{i,1} ')!'])
                   end
diff --git a/matlab/estimation/online/online_auxiliary_filter.m b/matlab/estimation/online/online_auxiliary_filter.m
index d1fce18e659e0fc459dd1c8e24b111a0b536a8a3..fbd693638d218533c2baea19796e5d0dacfc47de 100644
--- a/matlab/estimation/online/online_auxiliary_filter.m
+++ b/matlab/estimation/online/online_auxiliary_filter.m
@@ -40,7 +40,6 @@ function online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_,
 % 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;
 
@@ -102,7 +101,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);
@@ -315,9 +313,7 @@ for t=1:sample_size
         variance_update = false;
     end
     % final resampling (not advised)
-    if second_resample
-%        [~, idmode] = max(weights);
-%        mode_xparam(:,t) = xparam(:,idmode);
+    if online_opt.systematic_resampling
         indx = kitagawa(weights);
         StateVectors = StateVectors(:,indx) ;
         if pruning
@@ -340,8 +336,6 @@ for t=1:sample_size
             save(sprintf('%s%sparameters_particles_final.mat', SimulationFolder, filesep()), 'param');
         end 
     else
-%        [~, idmode] = max(weights);
-%        mode_xparam(:,t) = xparam(:,idmode);
         mean_xparam(:,t) = xparam*(weights');
         mat_var_cov = bsxfun(@minus, xparam,mean_xparam(:,t));
         mat_var_cov = mat_var_cov*(bsxfun(@times, mat_var_cov, weights)');