Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
Loading items
Show changes
Showing
with 879 additions and 156 deletions
function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMatlab)
% PARALLEL CONTEXT
% Core functionality for prior_posterior.m function, which can be parallelized.
% See also the comment in posterior_sampler_core.m funtion.
% See also the comment in posterior_sampler_core.m function.
%
% INPUTS
% See the comment in posterior_sampler_core.m funtion.
% See the comment in posterior_sampler_core.m function.
%
% OUTPUTS
% o myoutput [struc]
......@@ -76,7 +76,7 @@ naK=myinputs.naK;
horizon=myinputs.horizon;
iendo=myinputs.iendo;
IdObs=myinputs.IdObs; %index of observables
if horizon
if horizon && ~options_.occbin.smoother.status
i_last_obs=myinputs.i_last_obs;
MAX_nforc1=myinputs.MAX_nforc1;
MAX_nforc2=myinputs.MAX_nforc2;
......@@ -119,7 +119,8 @@ end
% DirectoryName = myinputs.DirectoryName;
if strcmpi(type,'posterior')
DirectoryName = CheckPath('metropolis',M_.dname);
folder_name=get_posterior_folder_name(options_);
DirectoryName = CheckPath(folder_name,M_.dname);
elseif strcmpi(type,'gsa')
if options_.opt_gsa.pprior
DirectoryName = CheckPath(['gsa',filesep,'prior'],M_.dname);
......@@ -167,7 +168,7 @@ if run_smoother
stock_smoothed_constant=NaN(endo_nbr,gend,MAX_n_smoothed_constant);
stock_smoothed_trend=NaN(endo_nbr,gend,MAX_n_smoothed_trend);
stock_trend_coeff = zeros(endo_nbr,MAX_n_trend_coeff);
if horizon
if horizon && ~options_.occbin.smoother.status
stock_forcst_mean= NaN(endo_nbr,horizon,MAX_nforc1);
stock_forcst_point = NaN(endo_nbr,horizon,MAX_nforc2);
if ~isequal(M_.H,0)
......@@ -229,12 +230,29 @@ for b=fpar:B
opts_local.occbin.simul.waitbar=0;
opts_local.occbin.smoother.waitbar = 0;
opts_local.occbin.smoother.linear_smoother=false; % speed-up
if options_.occbin.smoother.inversion_filter
dataset_.data=Y';
[~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, alphahat, etahat] = ...
occbin.IVF_posterior(deep,dataset_,[],options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if info(1)
message=get_error_message(info,opts_local);
fprintf('\nprior_posterior_statistics: IVF smoother failed for one of the draws:\n%s\n',message)
else
alphatilde = alphahat*nan;
SteadyState=oo_.dr.ys;
trend_coeff = zeros(length(options_.varobs_id),1);
trend_addition=zeros(options_.number_of_observed_variables,gend);
end
%epsilonhat not available as no measurement error allowed
else
opts_local.verbosity=0;
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_] = ...
occbin.DSGE_smoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);
if oo_.occbin.smoother.error_flag(1)
message=get_error_message(oo_.occbin.smoother.error_flag,opts_local);
fprintf('\nprior_posterior_statistics: One of the draws failed with the error:\n%s\n',message)
end
end
else
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,oo_,bayestopt_] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);
......@@ -310,7 +328,7 @@ for b=fpar:B
end
end
end
if horizon
if horizon && ~options_.occbin.smoother.status
yyyy = alphahat(iendo,i_last_obs);
yf = simulate_posterior_forecasts(yyyy,dr,horizon,false,M_.Sigma_e,1);
if options_.prefilter
......@@ -352,7 +370,7 @@ for b=fpar:B
stock_forcst_point(:,:,irun(7)) = yf1(maxlag+1:end,:)';
if ~isequal(M_.H,0)
ME_shocks=zeros(length(varobs),horizon);
i_exo_var = setdiff([1:length(varobs)],find(diag(M_.H) == 0));
i_exo_var = setdiff(1:length(varobs),find(diag(M_.H) == 0));
nxs = length(i_exo_var);
chol_H = chol(M_.H(i_exo_var,i_exo_var));
if ~isempty(M_.H)
......@@ -368,7 +386,7 @@ for b=fpar:B
stock_smoothed_uncert(dr.order_var,dr.order_var,:,irun(13)) = state_uncertainty;
end
else
[~,~,SteadyState,info] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~,~,SteadyState] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
stock_param(irun(5),:) = deep;
stock_logpo(irun(5),1) = logpo;
......@@ -518,9 +536,10 @@ for b=fpar:B
end
irun(irun_index) = 1;
end
if mod(b-fpar+1, 5)==0
dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
end
end
myoutput.ifil=ifil;
if RemoteFlag==1
......@@ -544,7 +563,7 @@ dyn_waitbar_close(h);
function yf=simulate_posterior_forecasts(y0,dr,horizon,stochastic_indicator,Sigma_e,n)
% function yf=forcst2(y0,horizon,dr,n)
% function yf=simulate_posterior_forecasts(y0,horizon,dr,n)
%
% computes forecasts based on first order model solution, given shocks
% drawn from the shock distribution, but not including measurement error
......
......@@ -6,6 +6,8 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
% M_ [structure] Model description.
% bayestopt_ [structure] Prior distribution description.
% options_ [structure] Global options of Dynare.
% oo_ [structure] Dynare structure where the results are saved.
% estim_params_ [structure] characterizing parameters to be estimated
%
% OUTPUTS:
% results [structure] Various statistics.
......@@ -13,7 +15,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
% SPECIAL REQUIREMENTS
% none
% Copyright © 2009-2023 Dynare Team
% Copyright © 2009-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -33,14 +35,11 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
% Initialization.
Prior = dprior(bayestopt_, options_.prior_trunc);
PriorDirectoryName = CheckPath('prior/draws',M_.dname);
work = ~drsave;
iteration = 0;
loop_indx = 0;
file_indx = [];
count_bk_indeterminacy = 0;
count_bk_unstability = 0;
count_bk_singularity = 0;
count_static_var_def = 0;
count_no_steadystate = 0;
count_steadystate_file_exit = 0;
count_dll_problem = 0;
......@@ -83,7 +82,7 @@ sampled_prior_covariance = zeros(NumberOfParameters,NumberOfParameters);
file_line_number = 0;
file_indx_number = 0;
oo_.dr=set_state_space(oo_.dr,M_,options_);
oo_.dr=set_state_space(oo_.dr,M_);
hh_fig = dyn_waitbar(0,'Please wait. Prior sampler...');
set(hh_fig,'Name','Prior sampler.');
......@@ -98,7 +97,7 @@ while iteration < NumberOfSimulations
M_ = set_all_parameters(params, estim_params_, M_);
[T, R, ~, INFO, oo_.dr,M_.params] = dynare_resolve(M_, options_, oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state, 'restrict');
if ~INFO(1)
INFO=endogenous_prior_restrictions(T,R,M_,options_,oo_);
INFO=endogenous_prior_restrictions(T,R,M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
file_line_number = file_line_number + 1;
iteration = iteration + 1;
......@@ -163,7 +162,7 @@ end
dyn_waitbar_close(hh_fig);
% Get informations about BK conditions and other things...
% Get information about BK conditions and other things...
results.bk.indeterminacy_share = count_bk_indeterminacy/loop_indx;
results.bk.unstability_share = count_bk_unstability/loop_indx;
results.bk.singularity_share = count_bk_singularity/loop_indx;
......
function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4, initialization)
function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4)
% Computes a prior density for the structural parameters of DSGE models
%
% INPUTS
......@@ -8,14 +8,13 @@ function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape,
% p7: [double] vector with n elements, second parameter of the prior distribution (bayestopt_.p7).
% p3: [double] vector with n elements, lower bounds of the untruncated standard or generalized distribution
% p4: [double] vector with n elements, upper bound of the untruncated standard or generalized distribution
% initialization [integer] if 1: initialize persistent variables
%
% OUTPUTS
% logged_prior_density [double] scalar, log of the prior density evaluated at x.
% info [double] error code for index of Inf-prior parameter
%
% Copyright © 2003-2023 Dynare Team
% Copyright © 2003-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -32,12 +31,8 @@ function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
persistent id1 id2 id3 id4 id5 id6 id8
persistent tt1 tt2 tt3 tt4 tt5 tt6 tt8
info=0;
if nargin > 6 && initialization
% Beta indices.
tt1 = true;
id1 = find(pshape==1);
......@@ -80,7 +75,6 @@ if nargin > 6 && initialization
if isempty(id8)
tt8 = false;
end
end
logged_prior_density = 0.0;
dlprior = zeros(1,length(x));
......@@ -95,9 +89,9 @@ if tt1
return
end
if nargout == 2
[tmp, dlprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1));
[~, dlprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1));
elseif nargout == 3
[tmp, dlprior(id1), d2lprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1));
[~, dlprior(id1), d2lprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1));
end
end
......@@ -110,18 +104,18 @@ if tt2
return
end
if nargout == 2
[tmp, dlprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2));
[~, dlprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2));
elseif nargout == 3
[tmp, dlprior(id2), d2lprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2));
[~, dlprior(id2), d2lprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2));
end
end
if tt3
logged_prior_density = logged_prior_density + sum(lpdfnorm(x(id3),p6(id3),p7(id3))) ;
if nargout == 2
[tmp, dlprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3));
[~, dlprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3));
elseif nargout == 3
[tmp, dlprior(id3), d2lprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3));
[~, dlprior(id3), d2lprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3));
end
end
......@@ -134,9 +128,9 @@ if tt4
return
end
if nargout == 2
[tmp, dlprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4));
[~, dlprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4));
elseif nargout == 3
[tmp, dlprior(id4), d2lprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4));
[~, dlprior(id4), d2lprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4));
end
end
......@@ -166,9 +160,9 @@ if tt6
return
end
if nargout == 2
[tmp, dlprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6));
[~, dlprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6));
elseif nargout == 3
[tmp, dlprior(id6), d2lprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6));
[~, dlprior(id6), d2lprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6));
end
end
......@@ -181,9 +175,9 @@ if tt8
return
end
if nargout==2
[tmp, dlprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8));
[~, dlprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8));
elseif nargout==3
[tmp, dlprior(id8), d2lprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8));
[~, dlprior(id8), d2lprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8));
end
end
......@@ -257,8 +251,8 @@ end
% Call the tested routine
try
% Initialization of priordens.
lpdstar = priordens(p5, p0, p6, p7, p3, p4, true);
% Prior (logged) density at the mode
lpdstar = priordens(p5, p0, p6, p7, p3, p4);
% Do simulations in a loop and estimate recursively the mean and teh variance.
LPD = NaN(10000,1);
for i = 1:10000
......
......@@ -6,10 +6,10 @@ function [Q,R] = qr2(varargin)
% & Schorfheides's identification scheme.
%
% INPUTS
% See matlab's documentation for QR decomposition.
% See MATLAB's documentation for QR decomposition.
%
% OUTPUTS
% See matlab's documentation for QR decomposition.
% See MATLAB's documentation for QR decomposition.
%
% ALGORITHM
% None.
......
......@@ -36,7 +36,7 @@ function [mu,sigma,offset] = recursive_moments(m0,s0,data,offset)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[T,n] = size(data);
[T,~] = size(data);
for t = 1:T
tt = t+offset;
......
function y = myprctilecol(x,p)
function indices = kitagawa(weights, noise)
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Return indices for resampling.
%
% INPUTS
% - weights [double] n×1 vector of partcles' weights.
% - noise [double] scalar, uniform random deviates in [0,1]
%
% OUTPUTS
% - indices [integer] n×1 vector of indices in [1:n]
% Copyright © 2012 European Commission
% Copyright © 2012-2017 Dynare Team
% Copyright © 2022-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -22,22 +26,20 @@ function y = myprctilecol(x,p)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
xx = sort(x);
[m,n] = size(x);
n= length(weights);
if nargin<2, noise = rand; end
indices = NaN(n, 1);
if m==1 | n==1
m = max(m,n);
if m == 1
y = x*ones(length(p),1);
return
cweights = cumsum(weights);
wweights = (transpose(0:n-1)+noise)*(1.0/n);
j = 1;
for i=1:n
while wweights(i)>cweights(j)
j = j+1;
end
n = 1;
q = 100*(0.5:m - 0.5)./m;
xx = [min(x); xx(:); max(x)];
else
q = 100*(0.5:m - 0.5)./m;
xx = [min(x); xx; max(x)];
indices(i) = j;
end
q = [0 q 100];
y = interp1(q,xx,p);
\ No newline at end of file
......@@ -2,12 +2,11 @@ function SampleAddress = selec_posterior_draws(M_,options_,dr,endo_steady_state,
% Selects a sample of draws from the posterior distribution and if nargin>1
% saves the draws in _pdraws mat files (metropolis folder). If drsize>0
% the dr structure, associated to the parameters, is also saved in _pdraws.
% This routine is more efficient than metropolis_draw.m because here an
% _mh file cannot be opened twice.
% This routine assures an _mh file cannot be opened twice.
%
% INPUTS
% o M_ [structure] Matlab's structure describing the model
% o options_ [structure] Matlab's structure describing the current options
% o M_ [structure] MATLAB's structure describing the model
% o options_ [structure] MATLAB's structure describing the current options
% o dr [structure] Reduced form model.
% o endo_steady_state [vector] steady state value for endogenous variables
% o exo_steady_state [vector] steady state value for exogenous variables
......@@ -26,7 +25,7 @@ function SampleAddress = selec_posterior_draws(M_,options_,dr,endo_steady_state,
% None.
%
% Copyright © 2006-2022 Dynare Team
% Copyright © 2006-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -55,7 +54,7 @@ switch nargin
case 8
info = 0;
case 9
MAX_mega_bytes = 10;% Should be an option...
MAX_mega_bytes = 100;% Should be an option...
if drsize>0
info=2;
else
......@@ -63,105 +62,58 @@ switch nargin
end
drawsize = drsize+npar*8/1048576;
otherwise
error(['selec_posterior_draws:: Unexpected number of input arguments!'])
error('selec_posterior_draws:: Unexpected number of input arguments!')
end
if ~issmc(options_)
MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName];
% Get informations about the mcmc:
record=load_last_mh_history_file(MetropolisFolder, ModelName);
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
LastMhFile = TotalNumberOfMhFiles;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
mh_nblck = options_.mh_nblck;
% Randomly select draws in the posterior distribution:
SampleAddress = zeros(SampleSize,4);
for i = 1:SampleSize
ChainNumber = ceil(rand*mh_nblck);
DrawNumber = ceil(rand*NumberOfDraws);
SampleAddress(i,1) = DrawNumber;
SampleAddress(i,2) = ChainNumber;
if DrawNumber <= MAX_nruns-FirstLine+1
MhFileNumber = FirstMhFile;
MhLineNumber = FirstLine+DrawNumber-1;
else
DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
if ishssmc(options_)
[MetropolisFolder] = CheckPath('hssmc',M_.dname);
elseif isdime(options_)
[MetropolisFolder] = CheckPath('dime',M_.dname);
else
error('check_posterior_analysis_data:: case should not happen. Please contact the developers')
end
SampleAddress(i,3) = MhFileNumber;
SampleAddress(i,4) = MhLineNumber;
end
SampleAddress = sortrows(SampleAddress,[3 2]);
ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName];
parameter_draws=get_posterior_subsample(M_,options_,SampleSize);
% Selected draws in the posterior distribution, and if drsize>0
% reduced form solutions, are saved on disk.
if info
%delete old stale files before creating new ones
delete_stale_file([BaseName '_posterior_draws*.mat'])
if SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file.
pdraws = cell(SampleSize,info);
old_mhfile = 0;
old_mhblck = 0;
for i = 1:SampleSize
mhfile = SampleAddress(i,3);
mhblck = SampleAddress(i,2);
if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
end
pdraws(i,1) = {x2(SampleAddress(i,4),:)};
if info==2
M_ = set_parameters_locally(M_,pdraws{i,1});
[dr,~,M_.params] =compute_decision_rules(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
pdraws(i,2) = { dr };
end
old_mhfile = mhfile;
old_mhblck = mhblck;
end
clear('x2')
save([BaseName '_posterior_draws1.mat'],'pdraws','estim_params_')
else% The posterior draws are saved in xx files.
NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
NumberOfLinesLastFile = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
linee = 0;
fnum = 1;
if fnum < NumberOfFiles
pdraws = cell(NumberOfDrawsPerFile,info);
old_mhfile = 0;
old_mhblck = 0;
else
pdraws = cell(NumberOfLinesLastFile,info);
end
for i=1:SampleSize
linee = linee+1;
mhfile = SampleAddress(i,3);
mhblck = SampleAddress(i,2);
if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
end
pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
pdraws(i,1) = {parameter_draws(i,:)};
if info==2
M_ = set_parameters_locally(M_,pdraws{i,1});
[dr,~,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
M_ = set_parameters_locally(M_,pdraws{linee,1});
[dr,~,M_.params] = compute_decision_rules(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
pdraws(linee,2) = { dr };
end
old_mhfile = mhfile;
old_mhblck = mhblck;
if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile
if (fnum < NumberOfFiles && linee == NumberOfDrawsPerFile) || (fnum <= NumberOfFiles && i==SampleSize)
linee = 0;
save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
fnum = fnum+1;
if fnum < NumberOfFiles
pdraws = cell(NumberOfDrawsPerFile,info);
else
pdraws = cell(NumberOfLines,info);
pdraws = cell(NumberOfLinesLastFile,info);
end
end
end
save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
end
end
\ No newline at end of file
function [sub_draws, error_flag, NumberOfDrawsPerChain]=set_number_of_subdraws(M_,options_)
% function [sub_draws, error_flag, NumberOfDrawsPerChain]=set_number_of_subdraws(M_,options_)
% Set option field sub_draws based on size of available sample
% Inputs:
% M_ [structure] Dynare model structure
% options_ [structure] Dynare options structure
%
% Outputs:
% sub_draws [scalar] number of sub-draws
% error_flag [boolean] error indicate of insufficient draws are
% available
% NumberOfDrawsPerChain [integer] number of available posterior draws
% per chain
% Copyright © 2024 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/>.
error_flag=false;
if ~issmc(options_)
record=load_last_mh_history_file([M_.dname filesep 'metropolis'], M_.fname);
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDrawsPerChain = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
NumberOfDraws=NumberOfDrawsPerChain*record.Nblck;
if isempty(options_.sub_draws)
sub_draws = min(options_.posterior_max_subsample_draws, ceil(NumberOfDraws));
else
if options_.sub_draws>NumberOfDraws*record.Nblck
skipline()
disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*options_.mh_nblck) ')!'])
disp('You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).')
skipline()
error_flag = true;
sub_draws=options_.sub_draws; %pass through
else
sub_draws=options_.sub_draws; %pass through
end
end
else
if ishssmc(options_)
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
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));
NumberOfDraws = size(posterior.param,2);
NumberOfDrawsPerChain=NumberOfDraws;
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
chains = reshape(posterior.chains(end-tune:end,:,:), [], size(posterior.chains, 3));
NumberOfDraws=size(chains,1);
NumberOfDrawsPerChain=NumberOfDraws;
else
error('set_number_of_subdraws:: case should not happen. Please contact the developers')
end
if isempty(options_.sub_draws)
sub_draws = min(options_.posterior_max_subsample_draws, NumberOfDraws);
else
if options_.sub_draws>NumberOfDraws
skipline()
disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws) ')!'])
disp('You can either change the value of sub_draws or run another MCMC.')
skipline()
error_flag = true;
sub_draws=options_.sub_draws; %pass through
else
sub_draws=options_.sub_draws; %pass through
end
end
end
\ No newline at end of file
......@@ -133,12 +133,12 @@ if ncn
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrn(:,10)]; %take generalized distribution into account
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrn(:,11)];
baseid = length(bayestopt_.name);
bayestopt_.name = [bayestopt_.name; cell(ncn, 1)];;
bayestopt_.name = [bayestopt_.name; cell(ncn, 1)];
for i=1:ncn
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
bayestopt_.name(baseid+i) = {sprintf('corr %s, %s', M_.endo_names{k1}, M_.endo_names{k2})};
% find correspondence to varobs to construct H in set_all_paramters
% find correspondence to varobs to construct H in set_all_parameters
obsi1 = strmatch(M_.endo_names{k1}, options_.varobs, 'exact');
obsi2 = strmatch(M_.endo_names{k2}, options_.varobs, 'exact');
% save correspondence
......@@ -163,8 +163,8 @@ bayestopt_.p7 = bayestopt_.p6 ;
%% check for point priors and disallow them as they do not work with MCMC
if any(bayestopt_.p2 ==0)
error(sprintf(['Error in prior for %s: you cannot use a point prior in estimation. Either increase the prior standard deviation',...
' or fix the parameter completely.'], bayestopt_.name{bayestopt_.p2 ==0}))
error('Error in prior for %s: you cannot use a point prior in estimation. Either increase the prior standard deviation',...
' or fix the parameter completely.', bayestopt_.name{bayestopt_.p2 ==0})
end
% generalized location parameters by default for beta distribution
......@@ -285,7 +285,7 @@ end
CheckPath('prior',M_.dname);
% I save the prior definition if the prior has changed.
if exist([ M_.dname '/prior/definition.mat'])
if exist([ M_.dname '/prior/definition.mat'],'file')
old = load([M_.dname '/prior/definition.mat'],'bayestopt_');
prior_has_changed = 0;
if length(bayestopt_.p1)==length(old.bayestopt_.p1)
......@@ -314,7 +314,3 @@ if exist([ M_.dname '/prior/definition.mat'])
else
save([M_.dname '/prior/definition.mat'],'bayestopt_');
end
% initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
bayestopt_.p3,bayestopt_.p4,1);
\ No newline at end of file
......@@ -5,7 +5,7 @@ function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetapri
% W: optimal value in the range (3,10)*std(x)
% - see C.Planas and A.Rossi (2014)
% objective_function(theta,varargin): -log of any unnormalized pdf
% with varargin (optional) a vector of auxiliaty parameters
% with varargin (optional) a vector of auxiliary parameters
% to be passed to f( ).
% ----------------------------------------------------------
%
......@@ -55,7 +55,6 @@ npar = length(theta);
W1 = sampler_options.W1;
neval = zeros(npar,1);
% % % fname0=fname;
fname = [ int2str(sampler_options.curr_block)];
Prior = dprior(varargin{6},varargin{3}.prior_trunc);
......
function dime(objective_function, init_x, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, dr , steady_state, exo_steady_state, exo_det_steady_state)
%dime(objective_function, init_x, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, dr , steady_state, exo_steady_state, exo_det_steady_state)
% Differential-Independence Mixture Ensemble ("DIME") MCMC sampling
% as proposed in "Ensemble MCMC Sampling for Robust Bayesian Inference"
% (Gregor Boehl, 2022, SSRN No. 4250395):
%
% https://gregorboehl.com/live/dime_mcmc_boehl.pdf
%
% INPUTS
% - objective_function [char] string specifying the name of the objective function (posterior kernel).
% - init_x [double] p×1 vector of parameters to be estimated (initial values, not used).
% - mh_bounds [double] p×2 matrix defining lower and upper bounds for the parameters.
% - dataset_ [dseries] sample
% - dataset_info [struct] information about the dataset
% - options_ [struct] Dynare's options
% - M_ [struct] model description
% - estim_params_ [struct] estimated parameters
% - bayestopt_ [struct] describing the priors
% - steady_state [double] steady state of endogenous variables
% - exo_steady_state [double] steady state of exogenous variables
% - exo_det_steady_state [double] steady state of exogenous deterministic variables
%
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2024 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/>.
% initialize
ndim = length(init_x);
prop_cov = eye(ndim);
prop_mean = zeros(1,ndim);
naccepted = 1;
cumlweight = -inf;
fixPSD = 1e-16*prop_cov;
% get options & assign problem specific default value
opts = options_.posterior_sampler_options.current_options;
if isfield(opts,'gamma')
g0 = opts.gamma;
else
g0 = 2.38 / sqrt(2 * ndim);
end
nchain = opts.nchain;
dft = opts.df_proposal_dist;
tune = opts.tune;
if tune > opts.niter
error('dime: number of ensemble iterations to keep (tune=%d) exceeds number of iterations (niter=%d)!', tune, opts.niter)
end
% temporary workaround to deal with parallelization
if opts.parallel
% check if parallel is possible
installed_toolboxes = ver;
toolbox_installed = any(strcmp("Parallel Computing Toolbox", {installed_toolboxes.Name}));
if ~toolbox_installed || feature('numcores') == 1
error('dime: parallel processing option is chosen but Parallel Computing Toolbox is not installed or machine has only one core')
end
end
% Set location for the simulated particles.
SimulationFolder = CheckPath('dime', M_.dname);
% Define prior distribution
Prior = dprior(bayestopt_, options_.prior_trunc);
bounds = check_bounds(Prior, mh_bounds);
% Set function handle for the objective
eval(sprintf('%s = @(x) %s(x, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, mh_bounds, dr , steady_state, exo_steady_state, exo_det_steady_state, []);', 'funobj', func2str(objective_function)));
% Initialization of the sampler (draws from the prior distribution with finite logged likelihood)
t0 = tic;
[x, ~, lprob] = ...
smc_samplers_initialization(funobj, 'dime', nchain, Prior, SimulationFolder, opts.niter);
x = ptransform(x', bounds, true);
disp_verbose(sprintf('Estimation:dime: log-posterior standard deviation (post.std) of a multivariate normal would be %.2f.\n', sqrt(0.5*ndim)), options_.verbosity);
% preallocate
chains = zeros(opts.niter, nchain, ndim);
lprobs = zeros(opts.niter, nchain);
isplit = fix(nchain/2) + 1;
for iter = 1:opts.niter
for complementary_ensemble = [false true]
% define current ensemble
if complementary_ensemble
idcur = 1:isplit;
idref = (isplit+1):nchain;
else
idcur = (isplit+1):nchain;
idref = 1:isplit;
end
cursize = length(idcur);
refsize = nchain - cursize;
% get differential evolution proposal
% draw the indices of the complementary chains
i1 = (0:cursize-1) + randi([1 cursize-1], 1, cursize);
i2 = (0:cursize-1) + randi([1 cursize-2], 1, cursize);
i2(i2 >= i1) = i2(i2 >= i1) + 1;
% add small noise and calculate proposal
f = opts.sigma*randn(cursize, 1);
q = x(idcur,:) + g0 .* (x(idref(mod(i1, refsize) + 1),:) - x(idref(mod(i2, refsize) + 1),:)) + f;
factors = zeros(cursize,1);
% get AIMH proposal
xchnge = rand(cursize,1) <= opts.aimh_prob;
% draw alternative candidates and calculate their proposal density
prop_cov_chol = chol(prop_cov*(dft - 2)/dft + fixPSD);
xcand = rand_multivariate_student(prop_mean, prop_cov_chol, dft, sum(xchnge));
lprop_old = log(multivariate_student_pdf(x(idcur(xchnge),:), prop_mean, prop_cov_chol, dft));
lprop_new = log(multivariate_student_pdf(xcand, prop_mean, prop_cov_chol, dft));
% update proposals and factors
q(xchnge,:) = xcand;
factors(xchnge) = lprop_old - lprop_new;
% Metropolis-Hastings
newlprob = log_prob_fun(funobj, Prior, bounds, opts.parallel, q);
lnpdiff = factors + newlprob - lprob(idcur);
accepted = lnpdiff > log(rand(cursize,1));
naccepted = naccepted + sum(accepted);
% update chains
x(idcur(accepted),:) = q(accepted,:);
lprob(idcur(accepted)) = newlprob(accepted);
end
% store
chains(iter,:,:) = ptransform(x, bounds, false);
lprobs(iter,:) = lprob;
% log weight of current ensemble
lweight = logsumexp(lprob) + log(naccepted) - log(nchain);
% calculate stats for current ensemble
ncov = cov(x);
nmean = mean(x);
% update AIMH proposal distribution
newcumlweight = logsumexp([cumlweight lweight]);
prop_cov = exp(cumlweight - newcumlweight) .* prop_cov + exp(lweight - newcumlweight) .* ncov;
prop_mean = exp(cumlweight - newcumlweight) .* prop_mean + exp(lweight - newcumlweight) .* nmean;
cumlweight = newcumlweight + log(opts.rho);
% be informative
tt = toc(t0);
if iter == 1 || ~mod(iter,15)
disp_verbose(' #iter. post.mode post.std %accept lapsed', options_.verbosity)
disp_verbose(' ------ --------- -------- ------- ------', options_.verbosity)
end
counter = sprintf('%u/%u', iter, opts.niter);
disp_verbose(sprintf('%11s %10.3f %1.2e %5.2f%% %5.2fs', counter, max(lprob), std(lprob), 100*naccepted/nchain, tt), options_.verbosity)
naccepted = 0;
end
save(sprintf('%s%schains.mat', SimulationFolder, filesep()), 'chains', 'lprobs', 'tune')
end
function lprobs = log_prob_fun(loglikefun, Prior, bounds, runs_in_parallel, x)
dim = length(x);
x = ptransform(x, bounds, false);
lprobs = zeros(dim,1);
if runs_in_parallel
parfor i=1:dim
par = x(i,:)';
loglikelihood = -loglikefun(par);
% temporary workaround for bug #1930
logprior = Prior.density(par);
lprobs(i) = loglikelihood + logprior;
end
else
for i=1:dim
lprobs(i) = -loglikefun(x(i,:)');
end
end
end
function nbounds = check_bounds(Prior, bounds)
% ensure that prior is continuous at bounds
nbounds = bounds;
prior_at_mean = Prior.density(Prior.mean);
dim = length(Prior.mean);
for j = 1:dim
if ~isinf(bounds.ub(j))
x = Prior.mean;
x(j) = bounds.ub(j);
if Prior.density(x - 1e-16) - prior_at_mean > log(1e-16)
nbounds.ub(j) = nbounds.ub(j) + 1e-2;
end
x = Prior.mean;
x(j) = bounds.lb(j);
if Prior.density(x + 1e-16) - prior_at_mean > log(1e-16)
nbounds.lb(j) = nbounds.lb(j) - 1e-2;
end
end
end
end
function x = ptransform(x, bounds, con2unc)
% parameter transformation
lb = bounds.lb;
ub = bounds.ub;
dim = size(x,2);
for j = 1:dim
% one-sided
if ~isinf(lb(j)) && isinf(ub(j))
if con2unc
x(:,j) = log(x(:,j) - lb(j));
else
x(:,j) = lb(j) + exp(x(:,j));
end
% two-sided
elseif ~isinf(lb(j)) && ~isinf(ub(j))
if con2unc
y = (x(:,j) - lb(j))/(ub(j) - lb(j));
x(:,j) = log(y ./ (1 - y));
else
x(:,j) = lb(j) + (ub(j) - lb(j)) ./ (1 + exp(-x(:,j)));
end
end
end
end
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_)
function dsmh(objective_function, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
% dsmh(objective_function, 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
% o objective_function [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.
......@@ -22,9 +22,9 @@ function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_
% 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.
% has been removed from this function and been placed in the posterior_sampler_core.m function.
%
% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
% 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.
%
......@@ -33,7 +33,7 @@ function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_
% Then the comments write here can be used for all the other pairs of
% parallel functions and also for management functions.
% Copyright © 2006-2023 Dynare Team
% Copyright © 2022-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -50,6 +50,7 @@ function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_
% 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;
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));
c = 0.055 ;
......@@ -57,7 +58,7 @@ MM = int64(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_
% Step 0: Initialization of the sampler
[param, tlogpost_iminus1, loglik, bayestopt_] = ...
SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,options_.posterior_sampler_options.dsmh.nparticles);
smc_samplers_initialization(objective_function, 'dsmh', opts.particles, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
ESS = zeros(options_.posterior_sampler_options.dsmh.H,1) ;
zhat = 1 ;
......@@ -67,47 +68,42 @@ for i=2:options_.posterior_sampler_options.dsmh.H
disp('');
disp('Tempered iteration');
disp(i) ;
% Step 1: sort the densities and compute IS weigths
% Step 1: sort the densities and compute IS weights
[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 = tune_c(objective_function,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,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
[param,tlogpost_iminus1,loglik] = mutation_DSMH(objective_function,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,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_.posterior_sampler_options.dsmh.nparticles);
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);
%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_.posterior_sampler_options.dsmh.nparticles) ;
ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.nparticles) ;
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
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_);
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(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');
......@@ -120,19 +116,14 @@ 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_);
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.dsmh.nparticles,bandwidth,kernel_function);
[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.nparticles,optimal_bandwidth,kernel_function);
options_.posterior_sampler_options.dsmh.particles,optimal_bandwidth,kernel_function);
plot(density(:,1),density(:,2));
hold on
if TeX
......@@ -179,7 +170,7 @@ 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_)
function c = tune_c(objective_function,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) ;
......@@ -198,7 +189,7 @@ while stop==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(objective_function,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)<exp(tlogpostx-tlogpost0(j)) % accept
......@@ -229,7 +220,7 @@ while stop==0
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_)
function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(objective_function,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;
npar = size(param,1) ;
p = 1/(10*options_.posterior_sampler_options.dsmh.tau);
......@@ -276,7 +267,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(objective_function,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<exp(tlogpostx-tlogpost_i0(j)) % accept
......
function mdd = hssmc(objective_function, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
% mdd = hssmc(objective_function, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
% Sequential Monte-Carlo sampler, Herbst and Schorfheide (JAE, 2014).
%
% INPUTS
% - objective_function [char] string specifying the name of the objective function (posterior kernel).
% - mh_bounds [double] p×2 matrix defining lower and upper bounds for the parameters.
% - dataset_ [dseries] sample
% - dataset_info [struct] Information about the dataset
% - options_ [struct] Dynare options
% - M_ [struct] model description
% - estim_params_ [struct] estimated parameters
% - bayestopt_ [struct] estimated parameters
% - oo_ [struct] outputs
%
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2022-2023 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/>.
smcopt = options_.posterior_sampler_options.current_options;
% Set location for the simulated particles.
SimulationFolder = CheckPath('hssmc', 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);
% 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(objective_function)));
mlogit = @(x) .95 + .1/(1 + exp(-16*x)); % Update of the scale parameter
% Create the tempering schedule
phi = ((0:smcopt.steps-1)/(smcopt.steps-1)).^smcopt.lambda;
% Initialize the estimate of the marginal density of the data
mdd = .0;
% tuning for MH algorithms matrices
scl = zeros(smcopt.steps, 1); % scale parameter
ESS = zeros(smcopt.steps, 1); % ESS
acpt = zeros(smcopt.steps, 1); % average acceptance rate
% Initialization of the sampler (draws from the prior distribution with finite logged likelihood)
t0 = tic;
[particles, tlogpostkernel, loglikelihood] = ...
smc_samplers_initialization(funobj, 'hssmc', smcopt.particles, Prior, SimulationFolder, smcopt.steps);
tt = toc(t0);
dprintf('#Iter. lambda ESS Acceptance rate scale resample seconds')
dprintf('%3u %5.4f %9.5E %5.4f %4.3f %3s %5.2f', 1, 0, 0, 0, 0, 'no', tt)
weights = ones(smcopt.particles, 1)/smcopt.particles;
resampled_particle_swarm = false;
for i=2:smcopt.steps % Loop over the weight on the likelihood (phi)
weights = weights.*exp((phi(i)-phi(i-1))*loglikelihood);
sweight = sum(weights);
weights = weights/sweight;
mdd = mdd + log(sweight);
ESS(i) = 1.0/sum(weights.^2);
if (2*ESS(i) < smcopt.particles) % Resampling
resampled_particle_swarm = true;
iresample = kitagawa(weights);
particles = particles(:,iresample);
loglikelihood = loglikelihood(iresample);
tlogpostkernel = tlogpostkernel(iresample);
weights = ones(smcopt.particles, 1)/smcopt.particles;
end
smcopt.scale = smcopt.scale*mlogit(smcopt.acpt-smcopt.target); % Adjust the scale parameter
scl(i) = smcopt.scale; % Scale parameter (for the jumping distribution in MH mutation step).
mu = particles*weights; % Weighted average of the particles.
z = particles-mu;
R = z*(z'.*weights); % Weighted covariance matrix of the particles.
t0 = tic;
acpt_ = false(smcopt.particles, 1);
tlogpostkernel = tlogpostkernel + (phi(i)-phi(i-1))*loglikelihood;
[acpt_, particles, loglikelihood, tlogpostkernel] = ...
randomwalk(funobj, chol(R, 'lower'), mu, scl(i), phi(i), acpt_, Prior, particles, loglikelihood, tlogpostkernel);
smcopt.acpt = sum(acpt_)/smcopt.particles; % Acceptance rate.
tt = toc(t0);
acpt(i) = smcopt.acpt;
if resampled_particle_swarm
dprintf('%3u %5.4f %9.5E %5.4f %4.3f %3s %5.2f', i, phi(i), ESS(i), acpt(i), scl(i), 'yes', tt)
else
dprintf('%3u %5.4f %9.5E %5.4f %4.3f %3s %5.2f', i, phi(i), ESS(i), acpt(i), scl(i), 'no', tt)
end
if i==smcopt.steps && ~resampled_particle_swarm
iresample = kitagawa(weights);
particles = particles(:,iresample);
end
if resampled_particle_swarm
save(sprintf('%s%sparticles-%u-%u.mat', SimulationFolder, filesep(), i, smcopt.steps), 'particles', 'tlogpostkernel', 'loglikelihood')
resampled_particle_swarm = false;
else
save(sprintf('%s%sparticles-%u-%u.mat', SimulationFolder, filesep(), i, smcopt.steps), 'particles', 'tlogpostkernel', 'loglikelihood', 'weights')
end
end
end
function [accept, particles, loglikelihood, tlogpostkernel] = randomwalk(funobj, RR, mu, scale, phi, accept, Prior, particles, loglikelihood, tlogpostkernel)
parfor j=1:size(particles, 2)
notvalid= true;
while notvalid
candidate = particles(:,j) + scale*(RR*randn(size(mu)));
if Prior.admissible(candidate)
[tlogpost, loglik] = tempered_likelihood(funobj, candidate, phi, Prior);
if isfinite(loglik)
notvalid = false;
if rand<exp(tlogpost-tlogpostkernel(j))
accept(j) = true ;
particles(:,j) = candidate;
loglikelihood(j) = loglik;
tlogpostkernel(j) = tlogpost;
end
end
end
end
end
end
function bool = isdime(options_)
% Copyright © 2024 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/>.
bool = false;
if isfield(options_, 'posterior_sampler_options')
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'dime_mcmc')
bool = true;
end
end
function bool = isdsmh(options_)
% Copyright © 2023 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/>.
bool = false;
if isfield(options_, 'posterior_sampler_options')
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'dsmh')
bool = true;
end
end
function bool = ishssmc(options_)
% Copyright © 2023 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/>.
bool = false;
if isfield(options_, 'posterior_sampler_options')
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'hssmc')
bool = true;
end
end
function [lse,sm] = logsumexp(x)
%LOGSUMEXP Log-sum-exp function.
% lse = LOGSUMEXP(x) returns the log-sum-exp function evaluated at
% the vector x, defined by lse = log(sum(exp(x)).
% [lse,sm] = LOGSUMEXP(x) also returns the softmax function evaluated
% at x, defined by sm = exp(x)/sum(exp(x)).
% The functions are computed in a way that avoids overflow and
% optimizes numerical stability.
%
% Reference:
% P. Blanchard, D. J. Higham, and N. J. Higham.
% Accurately computing the log-sum-exp and softmax functions.
% IMA J. Numer. Anal., Advance access, 2020.
%
% Taken from https://de.mathworks.com/matlabcentral/fileexchange/84892-logsumexp-softmax
%
% Copyright (c) 2020, Nicholas J. Higham
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% * Redistributions of source code must retain the above copyright notice, this
% list of conditions and the following disclaimer.
%
% * Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
if ~isvector(x), error('Input x must be a vector.'), end
n = length(x);
s = 0; e = zeros(n,1);
[xmax,k] = max(x); a = xmax;
s = 0;
for i = 1:n
e(i) = exp(x(i)-xmax);
if i ~= k
s = s + e(i);
end
end
lse = a + log1p(s);
if nargout > 1
sm = e/(1+s);
end