From 8f7356463442f0d882930d6f27c142cf4c9dac7b Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 19 Jul 2023 08:46:21 +0200 Subject: [PATCH 1/5] bug fix with non-zero lb bound of invgamma distribution --- matlab/distributions/inverse_gamma_specification.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m index e5136381a..e5ca2d3ba 100644 --- a/matlab/distributions/inverse_gamma_specification.m +++ b/matlab/distributions/inverse_gamma_specification.m @@ -85,7 +85,7 @@ s = []; nu = []; sigma = sqrt(sigma2); -mu2 = mu*mu; +mu2 = (mu-lb)*(mu-lb); if type == 2 % Inverse Gamma 2 nu = 2*(2+mu2/sigma2); @@ -132,10 +132,10 @@ elseif type == 1 % Inverse Gamma 1 end s = (sigma2+mu2)*(nu-2); if check_solution_flag - if abs(log(mu)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7 + if abs(log(mu-lb)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!'); end - if abs(sigma-sqrt(s/(nu-2)-mu*mu))>1e-7 + if abs(sigma-sqrt(s/(nu-2)-(mu-lb)*(mu-lb)))>1e-7 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!'); end end From de152a3de351d99d6731db42d2443e9c9702a668 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 19 Jul 2023 08:51:52 +0200 Subject: [PATCH 2/5] bug fix: indexing must also contain smpl+1 (needed for 1 step ahead forecast in last period when filtering). --- matlab/dynare_estimation_init.m | 8 ++++---- matlab/kalman/get_Qvec_heteroskedastic_filter.m | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index 4932e3042..01db632a8 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -500,18 +500,18 @@ if options_.heteroskedastic_filter 'for heteroskedastic_filter']) end - M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs); - M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs); + M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs+1); + M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs+1); for k=1:length(M_.heteroskedastic_shocks.Qvalue_orig) v = M_.heteroskedastic_shocks.Qvalue_orig(k); - temp_periods=v.periods(v.periods=options_.first_obs); M_.heteroskedastic_shocks.Qvalue(v.exo_id, temp_periods-(options_.first_obs-1)) = v.value^2; end for k=1:length(M_.heteroskedastic_shocks.Qscale_orig) v = M_.heteroskedastic_shocks.Qscale_orig(k); - temp_periods=v.periods(v.periods=options_.first_obs); M_.heteroskedastic_shocks.Qscale(v.exo_id, temp_periods-(options_.first_obs-1)) = v.scale^2; end diff --git a/matlab/kalman/get_Qvec_heteroskedastic_filter.m b/matlab/kalman/get_Qvec_heteroskedastic_filter.m index 7946b48d5..ae4e68868 100644 --- a/matlab/kalman/get_Qvec_heteroskedastic_filter.m +++ b/matlab/kalman/get_Qvec_heteroskedastic_filter.m @@ -27,7 +27,7 @@ function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,M_) isqdiag = all(all(abs(Q-diag(diag(Q)))<1e-14)); % ie, the covariance matrix is diagonal... Qvec=repmat(Q,[1 1 smpl+1]); -for k=1:smpl +for k=1:smpl+1 inx = ~isnan(M_.heteroskedastic_shocks.Qvalue(:,k)); if any(inx) if isqdiag From aad5c36081187fca1fabe221075ab591242ede89 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 19 Jul 2023 08:55:13 +0200 Subject: [PATCH 3/5] bug fix: with option mh_initialize_from_previous_mcmc, we need also to check if some prior changed, which may lead last draw in previous mcmc falling outside new prior bounds. --- matlab/posterior_sampler_initialization.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index a557cc652..f18ee5558 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -171,6 +171,12 @@ if ~options_.load_mh_file && ~options_.mh_recover ilogpo2 = zeros(NumberOfBlocks,1); ilogpo2(:,1) = record0.LastLogPost; ix2(:,IA) = record0.LastParameters(:,IB); + for j=1:NumberOfBlocks + if not(all(ix2(j,:)' >= mh_bounds.lb) && all(ix2(j,:)' <= mh_bounds.ub)) + new_estimated_parameters = logical(new_estimated_parameters + (ix2(j,:)' < mh_bounds.lb)); + new_estimated_parameters = logical(new_estimated_parameters + (ix2(j,:)' > mh_bounds.ub)); + end + end else new_estimated_parameters = true(1,npar); end From 53b57da8badce859de03c930f1cd5c5563f796da Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 19 Jul 2023 08:57:14 +0200 Subject: [PATCH 4/5] fix computation of initial prc0 under mh_recover (to avoid 0% being always displayed when recovery starts) --- matlab/posterior_sampler_core.m | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m index 551ba4daf..4b4ba5964 100644 --- a/matlab/posterior_sampler_core.m +++ b/matlab/posterior_sampler_core.m @@ -172,18 +172,6 @@ for curr_block = fblck:nblck logpo2 = zeros(InitSizeArray(curr_block),1); end end - %Prepare waiting bars - if whoiam - refresh_rate = sampler_options.parallel_bar_refresh_rate; - bar_title = sampler_options.parallel_bar_title; - prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode); - hh_fig = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); - else - refresh_rate = sampler_options.serial_bar_refresh_rate; - bar_title = sampler_options.serial_bar_title; - hh_fig = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); - set(hh_fig,'Name',bar_title); - end if mh_recover_flag==0 accepted_draws_this_chain = 0; accepted_draws_this_file = 0; @@ -192,6 +180,18 @@ for curr_block = fblck:nblck draw_iter = 1; draw_index_current_file = fline(curr_block); %get location of first draw in current block end + %Prepare waiting bars + if whoiam + refresh_rate = sampler_options.parallel_bar_refresh_rate; + bar_title = sampler_options.parallel_bar_title; + prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode)+(draw_iter-1)/nruns(curr_block); + hh_fig = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); + else + refresh_rate = sampler_options.serial_bar_refresh_rate; + bar_title = sampler_options.serial_bar_title; + hh_fig = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); + set(hh_fig,'Name',bar_title); + end sampler_options.curr_block = curr_block; while draw_iter <= nruns(curr_block) From f102a992aa91bdd050039b2833abbf9d8246f287 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 10 Nov 2023 17:51:41 +0100 Subject: [PATCH 5/5] fixed for the case when mcmc is incomplete WITHIN a block file (useful for expensive models and expensive methods like slice or TaRB) --- matlab/posterior_sampler_core.m | 47 ++++++++++++----------- matlab/posterior_sampler_initialization.m | 12 +++++- 2 files changed, 34 insertions(+), 25 deletions(-) diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m index 4b4ba5964..fb4bf8f99 100644 --- a/matlab/posterior_sampler_core.m +++ b/matlab/posterior_sampler_core.m @@ -141,31 +141,32 @@ for curr_block = fblck:nblck options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block); end mh_recover_flag=0; - if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood - load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat']) - x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)]; - logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)]; + if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2 && OpenOldFile(curr_block) + % this should be done whatever value of load_mh_file + load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']); + draw_iter = size(neval_this_chain,2)+1; + draw_index_current_file = draw_iter+fline(curr_block)-1; + feval_this_chain = sum(sum(neval_this_chain)); + feval_this_file = sum(sum(neval_this_chain)); + if feval_this_chain>draw_index_current_file-fline(curr_block) + % non Metropolis type of sampler + accepted_draws_this_chain = draw_index_current_file-fline(curr_block); + accepted_draws_this_file = draw_index_current_file-fline(curr_block); + else + accepted_draws_this_chain = 0; + accepted_draws_this_file = 0; + end + mh_recover_flag=1; + set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal); + last_draw(curr_block,:)=x2(draw_index_current_file-1,:); + last_posterior(curr_block)=logpo2(draw_index_current_file-1); OpenOldFile(curr_block) = 0; else - if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2 - load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']); - draw_iter = size(neval_this_chain,2)+1; - draw_index_current_file = draw_iter; - feval_this_chain = sum(sum(neval_this_chain)); - feval_this_file = sum(sum(neval_this_chain)); - if feval_this_chain>draw_iter-1 - % non Metropolis type of sampler - accepted_draws_this_chain = draw_iter-1; - accepted_draws_this_file = draw_iter-1; - else - accepted_draws_this_chain = 0; - accepted_draws_this_file = 0; - end - mh_recover_flag=1; - set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal); - last_draw(curr_block,:)=x2(draw_iter-1,:); - last_posterior(curr_block)=logpo2(draw_iter-1); - + if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood + load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat']) + x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)]; + logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)]; + OpenOldFile(curr_block) = 0; else x2 = zeros(InitSizeArray(curr_block),npar); diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index f18ee5558..d80644b1b 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -465,7 +465,7 @@ elseif options_.mh_recover AllMhFiles = dir([BaseName '_mh*_blck*.mat']); TotalNumberOfMhFiles = length(AllMhFiles)-length(dir([BaseName '_mh_tmp*_blck*.mat'])); % Quit if no crashed mcmc chain can be found as there are as many files as expected - if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) + if (ExpectedNumberOfMhFilesPerBlock>LastFileNumberInThePreviousMh) && (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) if isnumeric(options_.parallel) fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString); fprintf(' You have to edit the mod file and remove the mh_recover option\n'); @@ -476,15 +476,23 @@ elseif options_.mh_recover % 2. Something needs to be done; find out what % Count the number of saved mh files per block. NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1); + is_chain_complete = true(NumberOfBlocks,1); for b = 1:NumberOfBlocks BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']); NumberOfMhFilesPerBlock(b) = length(BlckMhFiles)-length(dir([BaseName '_mh_tmp*_blck' int2str(b) '.mat'])); + if (ExpectedNumberOfMhFilesPerBlock==LastFileNumberInThePreviousMh) + % here we need to check for the number of lines + tmpdata = load([BaseName '_mh' int2str(LastFileNumberInThePreviousMh) '_blck' int2str(b) '.mat'],'logpo2'); + if length(tmpdata.logpo2) == LastLineNumberInThePreviousMh + is_chain_complete(b) = false; + end + end end % Find FirstBlock (First block), an integer targeting the crashed mcmc chain. FirstBlock = 1; %initialize FBlock = zeros(NumberOfBlocks,1); while FirstBlock <= NumberOfBlocks - if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock + if (NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock) || not(is_chain_complete(FirstBlock)) fprintf('%s: Chain %u is not complete!\n', dispString,FirstBlock); FBlock(FirstBlock)=1; else