diff --git a/matlab/nonlinear-filters/src/DSMH_initialization.m b/matlab/nonlinear-filters/src/DSMH_initialization.m
deleted file mode 100644
index a9b9ee1bd..000000000
--- a/matlab/nonlinear-filters/src/DSMH_initialization.m
+++ /dev/null
@@ -1,118 +0,0 @@
-function [ ix2, temperedlogpost, loglik, ModelName, MetropolisFolder, npar, NumberOfParticles, bayestopt_] = ...
- DSMH_initialization(TargetFun, xparam1, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
-% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfParticles, bayestopt_] = ...
-% DSMH_initialization(TargetFun, xparam1, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
-% Dynamic Striated Metropolis-Hastings initialization.
-%
-% INPUTS
-% o TargetFun [char] string specifying the name of the objective
-% function (tempered posterior kernel and likelihood).
-% 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
-% o options_ options structure
-% o M_ model structure
-% o estim_params_ estimated parameters structure
-% o bayestopt_ estimation options structure
-% o oo_ outputs structure
-%
-% OUTPUTS
-% o ix2 [double] (NumberOfParticles*npar) vector of starting points for different chains
-% o ilogpo2 [double] (NumberOfParticles*1) vector of initial posterior values for different chains
-% o iloglik2 [double] (NumberOfParticles*1) vector of initial likelihood values for different chains
-% o ModelName [string] name of the mod-file
-% o MetropolisFolder [string] path to the Metropolis subfolder
-% o npar [scalar] number of parameters estimated
-% o NumberOfParticles [scalar] Number of particles requested for the parameters distributions
-% o bayestopt_ [structure] estimation options structure
-%
-% SPECIAL REQUIREMENTS
-% None.
-
-% Copyright (C) 2006-2017 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 .
-
-%Initialize outputs
-ix2 = [];
-ilogpo2 = [];
-iloglik2 = [];
-ModelName = [];
-MetropolisFolder = [];
-npar = [];
-NumberOfParticles = [];
-
-ModelName = M_.fname;
-if ~isempty(M_.bvar)
- ModelName = [ModelName '_bvar'];
-end
-
-MetropolisFolder = CheckPath('dsmh',M_.dname);
-BaseName = [MetropolisFolder filesep ModelName];
-
-NumberOfParticles = options_.dsmh.number_of_particles; %Number of particles for the parameters
-npar = length(xparam1);
-
-% Here we start a new DS Metropolis-Hastings, previous draws are discarded.
-disp('Estimation::dsmh: Initialization...')
-% Delete old dsmh files if any...
-files = dir([BaseName '_dsmh*_blck*.mat']);
-%if length(files)
-% delete([BaseName '_dsmh*_blck*.mat']);
-% disp('Estimation::smc: Old dsmh-files successfully erased!')
-%end
-% Delete old log file.
-file = dir([ MetropolisFolder '/dsmh.log']);
-%if length(file)
-% delete([ MetropolisFolder '/dsmh.log']);
-% disp('Estimation::dsmh: Old dsmh.log file successfully erased!')
-% disp('Estimation::dsmh: Creation of a new dsmh.log file.')
-%end
-fidlog = fopen([MetropolisFolder '/dsmh.log'],'w');
-fprintf(fidlog,'%% DSMH log file (Dynare).\n');
-fprintf(fidlog,['%% ' datestr(now,0) '.\n']);
-fprintf(fidlog,' \n\n');
-fprintf(fidlog,'%% Session 1.\n');
-fprintf(fidlog,' \n');
-prior_draw(bayestopt_,options_.prior_trunc);
-% Find initial values for the NumberOfParticles chains...
-set_dynare_seed('default');
-fprintf(fidlog,[' Initial values of the parameters:\n']);
-disp('Estimation::dsmh: Searching for initial values...');
-ix2 = zeros(npar,NumberOfParticles);
-temperedlogpost = zeros(NumberOfParticles,1);
-loglik = zeros(NumberOfParticles,1);
-%stderr = sqrt(bsxfun(@power,mh_bounds.ub-mh_bounds.lb,2)/12)/10;
-for j=1:NumberOfParticles
- validate = 0;
- while validate == 0
- candidate = prior_draw()';
-% candidate = xparam1(:) + 0.001*randn(npar,1);%bsxfun(@times,stderr,randn(npar,1)) ;
- if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
- ix2(:,j) = candidate ;
- [temperedlogpost(j),loglik(j)] = tempered_likelihood(TargetFun,candidate,0.0,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
- if isfinite(loglik(j)) % if returned log-density is Inf or Nan (penalized value)
- validate = 1;
- end
- end
- end
-end
-fprintf(fidlog,' \n');
-disp('Estimation::dsmh: Initial values found!')
-skipline()
-
-
diff --git a/matlab/nonlinear-filters/src/DSMH_sampler.m b/matlab/nonlinear-filters/src/DSMH_sampler.m
deleted file mode 100644
index 9c522f2a1..000000000
--- a/matlab/nonlinear-filters/src/DSMH_sampler.m
+++ /dev/null
@@ -1,301 +0,0 @@
-function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
-% function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
-% Dynamic Striated Metropolis-Hastings algorithm.
-%
-% 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
-% o options_ options structure
-% o M_ model structure
-% o estim_params_ estimated parameters structure
-% o bayestopt_ estimation options structure
-% o oo_ outputs structure
-%
-% 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 (C) 2006-2021 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 .
-
-
-lambda = exp(bsxfun(@minus,options_.dsmh.H,1:1:options_.dsmh.H)/(options_.dsmh.H-1)*log(options_.dsmh.lambda1));
-c = 0.055 ;
-
-% Step 0: Initialization of the sampler
-[ param, tlogpost_iminus1, loglik, npar, ~, bayestopt_] = ...
- SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,options_.dsmh.number_of_particles);
-
-ESS = zeros(options_.dsmh.H,1) ;
-zhat = 1 ;
-
-% The DSMH starts here
-for i=2:options_.dsmh.H
- disp('');
- disp('Tempered iteration');
- disp(i) ;
- % Step 1: sort the densities and compute IS weigths
- [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param) ;
- [tlogpost_i,weights,zhat,ESS,mu,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_);
- % Step 3: Metropolis step
- [param,tlogpost_iminus1,loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-end
-
-weights = exp(loglik*(lambda(end)-lambda(end-1)));
-weights = weights/sum(weights);
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.dsmh.number_of_particles);
-distrib_param = param(:,indx_resmpl);
-
-mean_xparam = mean(distrib_param,2);
-%mat_var_cov = bsxfun(@minus,distrib_param,mean_xparam) ;
-%mat_var_cov = (mat_var_cov*mat_var_cov')/(options_.HSsmc.nparticles-1) ;
-%std_xparam = sqrt(diag(mat_var_cov)) ;
-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_.HSsmc.nparticles) ;
- ub95_xparam(i) = temp(0.975*options_.HSsmc.nparticles) ;
-end
-
-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_);
- 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
-
-[nbplt,nr,nc,lr,lc,nstar] = pltorg(npar);
-
-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 = 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_);
- optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.dsmh.number_of_particles,bandwidth,kernel_function);
- [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
- options_.dsmh.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,[ 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
-
-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,mu,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 ;
-weights = weights/sum(weights);
-ESS(i) = 1/sum(weights.^2);
-% estimates of mean and variance
-mu = param*weights;
-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) ;
-npar = size(param,1);
-lower_prob = (.5*(options_.dsmh.alpha0+options_.dsmh.alpha1))^5;
-upper_prob = (.5*(options_.dsmh.alpha0+options_.dsmh.alpha1))^(1/5);
-stop=0 ;
-while stop==0
- acpt = 0.0;
- indx_resmpl = smc_resampling(weights,rand(1,1),options_.dsmh.G);
- param0 = param(:,indx_resmpl);
- tlogpost0 = tlogpost_i(indx_resmpl);
- for j=1:options_.dsmh.G
- for l=1:options_.dsmh.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_);
- if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
- validate = 1;
- if rand(1,1)= 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_);
- if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
- validate = 1;
- if u2.
-
-
-% Create the tempering schedule
-phi = bsxfun(@power,(bsxfun(@minus,1:1:options_.HSsmc.nphi,1)/(options_.HSsmc.nphi-1)),options_.HSsmc.lambda) ;
-% tuning for MH algorithms matrices
-zhat = 0 ; % normalization constant
-csim = zeros(options_.HSsmc.nphi,1) ; % scale parameter
-ESSsim = zeros(options_.HSsmc.nphi,1) ; % ESS
-acptsim = zeros(options_.HSsmc.nphi,1) ; % average acceptance rate
-% Step 0: Initialization of the sampler
-[ param, tlogpost_i, loglik, npar, ~, bayestopt_] = ...
- SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,options_.HSsmc.nparticles);
-weights = ones(options_.HSsmc.nparticles,1)/options_.HSsmc.nparticles ;
-% The Herbst and Schorfheide sampler starts here
-for i=2:options_.HSsmc.nphi
- % (a) Correction
- % incremental weights
- incwt = exp((phi(i)-phi(i-1))*loglik) ;
- % update weights
- weights = bsxfun(@times,weights,incwt) ;
- sum_weights = sum(weights) ;
- zhat = zhat + log(sum_weights) ;
- % normalize weights
- weights = weights/sum_weights ;
- % (b) Selection
- ESSsim(i) = 1/sum(weights.^2) ;
- if (ESSsim(i) < options_.HSsmc.nparticles/2)
- indx_resmpl = smc_resampling(weights,rand(1,1),options_.HSsmc.nparticles) ;
- param = param(:,indx_resmpl) ;
- loglik = loglik(indx_resmpl) ;
- tlogpost_i = tlogpost_i(indx_resmpl) ;
- weights = ones(options_.HSsmc.nparticles,1)/options_.HSsmc.nparticles ;
- end
- % (c) Mutation
- options_.HSsmc.c = options_.HSsmc.c*modified_logit(0.95,0.1,16.0,options_.HSsmc.acpt-options_.HSsmc.trgt) ;
- % Calculate estimates of mean and variance
- mu = param*weights ;
- z = bsxfun(@minus,param,mu) ;
- R = z*(bsxfun(@times,z',weights)) ;
- Rchol = chol(R)' ;
- % Mutation
- if options_.HSsmc.option_mutation==1
- [param,tlogpost_i,loglik,options_.HSsmc.acpt] = mutation_RW(TargetFun,param,tlogpost_i,loglik,phi,i,options_.HSsmc.c*Rchol,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
- elseif options_.HSsmc.option_mutation==2
- inv_R = inv(options_.HSsmc.c^2*R) ;
- Rdiagchol = sqrt(diag(R)) ;
- [param,tlogpost_i,loglik,options_.HSsmc.acpt] = mutation_Mixture(TargetFun,param,tlogpost_i,loglik,phi,i,options_.HSsmc.c*Rchol,options_.HSsmc.c*Rdiagchol,inv_R,mu,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
- end
- acptsim(i) = options_.HSsmc.acpt ;
- csim(i) = options_.HSsmc.c ;
- % print information
- fprintf(' Iteration = %5.0f / %5.0f \n', i, options_.HSsmc.nphi);
- fprintf(' phi = %5.4f \n', phi(i));
- fprintf(' Neff = %5.4f \n', ESSsim(i));
- fprintf(' %accept. = %5.4f \n', acptsim(i));
-end
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.HSsmc.nparticles);
-distrib_param = param(:,indx_resmpl);
-fprintf(' Log_lik = %5.4f \n', zhat);
-
-mean_xparam = mean(distrib_param,2);
-%mat_var_cov = bsxfun(@minus,distrib_param,mean_xparam) ;
-%mat_var_cov = (mat_var_cov*mat_var_cov')/(options_.HSsmc.nparticles-1) ;
-%std_xparam = sqrt(diag(mat_var_cov)) ;
-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_.HSsmc.nparticles) ;
- ub95_xparam(i) = temp(0.975*options_.HSsmc.nparticles) ;
-end
-
-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_);
- 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
-
-[nbplt,nr,nc,lr,lc,nstar] = pltorg(npar);
-
-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 = 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_);
- optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.HSsmc.nparticles,bandwidth,kernel_function);
- [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
- options_.HSsmc.nparticles,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,[ 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 Herbst/Schorfheide sampler.}');
- fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
- fprintf(fidTeX,'\\end{figure}\n');
- fprintf(fidTeX,' \n');
-end
-%end
-
-function [param,tlogpost_i,loglik,acpt] = mutation_RW(TargetFun,param,tlogpost_i,loglik,phi,i,cRchol,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-acpt = 0.0 ;
-tlogpost_i = tlogpost_i + (phi(i)-phi(i-1))*loglik ;
-for j=1:options_.HSsmc.nparticles
- validate= 0;
- while validate==0
- candidate = param(:,j) + cRchol*randn(size(param,1),1) ;
- if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
- [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
- if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
- validate = 1;
- if rand(1,1)= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
- [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
- if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
- validate = 1;
- if rand(1,1).
-
-%Initialize outputs
-ix2 = [];
-ilogpo2 = [];
-iloglik2 = [];
-ModelName = [];
-MetropolisFolder = [];
-npar = [];
-
-ModelName = M_.fname;
-if ~isempty(M_.bvar)
- ModelName = [ModelName '_bvar'];
-end
-
-MetropolisFolder = CheckPath('dsmh',M_.dname);
-BaseName = [MetropolisFolder filesep ModelName];
-
-npar = length(xparam1);
-
-% Here we start a new DS Metropolis-Hastings, previous draws are discarded.
-disp('Estimation:: Initialization...')
-% Delete old dsmh files if any...
-files = dir([BaseName '_dsmh*_blck*.mat']);
-%if length(files)
-% delete([BaseName '_dsmh*_blck*.mat']);
-% disp('Estimation::smc: Old dsmh-files successfully erased!')
-%end
-% Delete old log file.
-file = dir([ MetropolisFolder '/dsmh.log']);
-%if length(file)
-% delete([ MetropolisFolder '/dsmh.log']);
-% disp('Estimation::dsmh: Old dsmh.log file successfully erased!')
-% disp('Estimation::dsmh: Creation of a new dsmh.log file.')
-%end
-fidlog = fopen([MetropolisFolder '/dsmh.log'],'w');
-fprintf(fidlog,'%% DSMH log file (Dynare).\n');
-fprintf(fidlog,['%% ' datestr(now,0) '.\n']);
-fprintf(fidlog,' \n\n');
-fprintf(fidlog,'%% Session 1.\n');
-fprintf(fidlog,' \n');
-prior_draw(bayestopt_,options_.prior_trunc);
-% Find initial values for the NumberOfParticles chains...
-set_dynare_seed('default');
-fprintf(fidlog,[' Initial values of the parameters:\n']);
-disp('Estimation::dsmh: Searching for initial values...');
-ix2 = zeros(npar,NumberOfParticles);
-temperedlogpost = zeros(NumberOfParticles,1);
-loglik = zeros(NumberOfParticles,1);
-%stderr = sqrt(bsxfun(@power,mh_bounds.ub-mh_bounds.lb,2)/12)/10;
-for j=1:NumberOfParticles
- validate = 0;
- while validate == 0
- candidate = prior_draw()';
-% candidate = xparam1(:) + 0.001*randn(npar,1);%bsxfun(@times,stderr,randn(npar,1)) ;
- if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
- ix2(:,j) = candidate ;
- [temperedlogpost(j),loglik(j)] = tempered_likelihood(TargetFun,candidate,0.0,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
- if isfinite(loglik(j)) % if returned log-density is Inf or Nan (penalized value)
- validate = 1;
- end
- end
- end
-end
-fprintf(fidlog,' \n');
-disp('Estimation:: Initial values found!')
-skipline()
-
-
diff --git a/matlab/nonlinear-filters/src/smc_resampling.m b/matlab/nonlinear-filters/src/smc_resampling.m
deleted file mode 100644
index bac7eeac0..000000000
--- a/matlab/nonlinear-filters/src/smc_resampling.m
+++ /dev/null
@@ -1,11 +0,0 @@
-function indx = smc_resampling(weights,noise,number)
- indx = zeros(number,1);
- cumweights = cumsum(weights);
- randvec = (transpose(1:number)-1+noise(:))/number;
- j = 1;
- for i=1:number
- while (randvec(i)>cumweights(j))
- j = j+1;
- end
- indx(i) = j;
- end
diff --git a/matlab/nonlinear-filters/src/tempered_likelihood.m b/matlab/nonlinear-filters/src/tempered_likelihood.m
deleted file mode 100644
index 4af11080a..000000000
--- a/matlab/nonlinear-filters/src/tempered_likelihood.m
+++ /dev/null
@@ -1,5 +0,0 @@
- function [tlogpostkern,loglik] = tempered_likelihood(TargetFun,xparam1,lambda,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_)
- logpostkern = -feval(TargetFun,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
- logprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
- loglik = logpostkern-logprior ;
- tlogpostkern = lambda*loglik + logprior;