diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m
index 9efabb4d89b00fac5a65608aee9e4e4436e09d74..3e187a19eae79f87a3548e6b0670d9c4050d4d06 100644
--- a/matlab/posterior_sampler.m
+++ b/matlab/posterior_sampler.m
@@ -61,7 +61,7 @@ objective_function_penalty_base = Inf;
vv = sampler_options.invhess;
% Initialization of the sampler
-[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
+[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m
index ec749de63833aa0ccc7cf7e02b2d8d423ee713be..c677e02cceb28627103039868b0ee1948ff569e0 100644
--- a/matlab/posterior_sampler_core.m
+++ b/matlab/posterior_sampler_core.m
@@ -110,6 +110,9 @@ end
sampler_options.xparam1 = xparam1;
if ~isempty(d),
sampler_options.proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
+ %store information for load_mh_file
+ record.ProposalCovariance=d;
+ record.ProposalScaleVec=bayestopt_.jscale;
end
block_iter=0;
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index 75df703c8698d5f740f67b49528f31ed658f98e3..3adca10a50105e8ae45dfdb6c0ab13bfe47b19c1 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -1,6 +1,6 @@
-function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ...
+function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
-% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ...
+% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% Metropolis-Hastings initialization.
%
@@ -33,6 +33,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa
% draws
% o MAX_nruns [scalar] maximum number of draws in each MH-file on the harddisk
% o d [double] (p*p) matrix, Cholesky decomposition of the posterior covariance matrix (at the mode).
+% o bayestopt_ [structure] estimation options structure
%
% SPECIAL REQUIREMENTS
% None.
@@ -224,6 +225,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
record.LastFileNumber = AnticipatedNumberOfFiles ;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
record.MCMCConcludedSuccessfully = 0;
+ record.MCMC_sampler=options_.posterior_sampler_options.posterior_sampling_method;
fprintf('Ok!\n');
id = write_mh_history_file(MetropolisFolder, ModelName, record);
disp(['Estimation::mcmc: Details about the MCMC are available in ' BaseName '_mh_history_' num2str(id) '.mat'])
@@ -283,6 +285,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
end
ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
+ [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_);
FirstBlock = 1;
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
fprintf('Estimation::mcmc: I am writing a new mh-history file... ');
@@ -363,7 +366,7 @@ elseif options_.mh_recover
FirstLine = ones(NumberOfBlocks,1);
LastFileFullIndicator=1;
end
-
+ [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_);
%% Now find out what exactly needs to be redone
% 1. Check if really something needs to be done
% How many mh files should we have ?
@@ -447,3 +450,21 @@ elseif options_.mh_recover
write_mh_history_file(MetropolisFolder, ModelName, record);
end
end
+
+function [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_)
+ if isfield(record,'ProposalCovariance') && isfield(record,'ProposalCovariance')
+ if isfield(record,'MCMC_sampler')
+ if ~strcmp(record.MCMC_sampler,options_.posterior_sampler_options.posterior_sampling_method)
+ error(fprintf('Estimation::mcmc: The posterior_sampling_method differs from the one of the original chain. Please reset it to %s',record.MCMC_sampler))
+ end
+ end
+ fprintf('Estimation::mcmc: Recovering the previous proposal density\n.')
+ d=record.ProposalCovariance;
+ bayestopt_.jscale=record.ProposalScaleVec;
+ else
+ if options_.mode_compute~=0
+ fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by mode_compute\n.');
+ elseif ~isempty(options_.mode_file)
+ fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by the mode_file\n.');
+ end
+ end
diff --git a/tests/estimation/fs2000.mod b/tests/estimation/fs2000.mod
index 2eaeff302b9b54846d4e339dbbb7f1d08466fe58..640f1961f96b49a1d5f6fc79ae497e3fd31ade69 100644
--- a/tests/estimation/fs2000.mod
+++ b/tests/estimation/fs2000.mod
@@ -82,4 +82,21 @@ varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
-estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=3000,mh_nblocks=2,mh_jscale=0.8,moments_varendo,selected_variables_only,contemporaneous_correlation,smoother,forecast=8) y m;
+estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=3000,mh_nblocks=1,mh_jscale=0.8,moments_varendo,selected_variables_only,contemporaneous_correlation,smoother,forecast=8) y m;
+%test load_mh_file option
+options_.smoother=0;
+options_.moments_varendo=0;
+options_.forecast=0;
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
+estimation(mode_compute=0,mode_file=fs2000_mode,order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=1500, mh_nblocks=1, mh_jscale=0.8);
+hh=eye(size(bayestopt_.name,1));
+save('fs2000_mode','hh','-append')
+estimation(mode_compute=0,mode_file=fs2000_mode,order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=1500, mh_nblocks=1, mh_jscale=10,load_mh_file);
+
+
+temp1=load([M_.dname '_mh1_blck1.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of unaffected chain are not the same')
+end