From 494248828735d762052877abc264d305146c366f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Karam=C3=A9?= Date: Wed, 9 Mar 2022 17:21:02 +0100 Subject: [PATCH] Routines moved out of the submodule. --- src/DSMH_initialization.m | 118 ------------ src/DSMH_sampler.m | 301 ------------------------------ src/Herbst_Schorfheide_sampler.m | 245 ------------------------ src/SMC_samplers_initialization.m | 117 ------------ src/smc_resampling.m | 11 -- src/tempered_likelihood.m | 5 - 6 files changed, 797 deletions(-) delete mode 100644 src/DSMH_initialization.m delete mode 100644 src/DSMH_sampler.m delete mode 100644 src/Herbst_Schorfheide_sampler.m delete mode 100644 src/SMC_samplers_initialization.m delete mode 100644 src/smc_resampling.m delete mode 100644 src/tempered_likelihood.m diff --git a/src/DSMH_initialization.m b/src/DSMH_initialization.m deleted file mode 100644 index a9b9ee1bd..000000000 --- a/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/src/DSMH_sampler.m b/src/DSMH_sampler.m deleted file mode 100644 index 9c522f2a1..000000000 --- a/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/src/smc_resampling.m b/src/smc_resampling.m deleted file mode 100644 index bac7eeac0..000000000 --- a/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/src/tempered_likelihood.m b/src/tempered_likelihood.m deleted file mode 100644 index 4af11080a..000000000 --- a/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;