diff --git a/doc/dynare.texi b/doc/dynare.texi
index 94278ba689ebeda81060aac22a27dc096b1d1805..a0e987320a014331fb07b9623dc23ed2f48b64e1 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -4986,8 +4986,12 @@ Metropolis-Hastings chain. Default: 2*@code{mh_scale}
@item mh_recover
@anchor{mh_recover} Attempts to recover a Metropolis-Hastings
-simulation that crashed prematurely. Shouldn't be used together with
-@code{load_mh_file}
+simulation that crashed prematurely, starting with the last available saved
+@code{mh}-file. Shouldn't be used together with
+@code{load_mh_file} or a different @code{mh_replic} than in the crashed run.
+To assure a neat continuation of the chain with the same proposal density, you should
+provide the @code{mode_file} used in the previous
+run or the same user-defined @code{mcmc_jumping_covariance} when using this option.
@item mh_mode = @var{INTEGER}
@dots{}
diff --git a/matlab/TaRB_metropolis_hastings_core.m b/matlab/TaRB_metropolis_hastings_core.m
index 7259c8ff8def2e99da6a224c484f1a0c090d2f1b..f8e5d89151c0c9e0e6243aae55d735d942cd7232 100644
--- a/matlab/TaRB_metropolis_hastings_core.m
+++ b/matlab/TaRB_metropolis_hastings_core.m
@@ -236,7 +236,8 @@ for curr_chain = fblck:nblck,
dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_chain) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/blocked_draws_counter_this_chain)]);
end
if (draw_index_current_file == InitSizeArray(curr_chain)) || (draw_iter == nruns(curr_chain)) % Now I save the simulations, either because the current file is full or the chain is done
- save([BaseName '_mh' int2str(NewFile(curr_chain)) '_blck' int2str(curr_chain) '.mat'],'x2','logpo2');
+ [LastSeeds.(['file' int2str(NewFile(curr_chain))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_chain))]).Normal] = get_dynare_random_generator_state();
+ save([BaseName '_mh' int2str(NewFile(curr_chain)) '_blck' int2str(curr_chain) '.mat'],'x2','logpo2','LastSeeds');
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_chain)) 'Blck' int2str(curr_chain) ' (' datestr(now,0) ')\n']);
diff --git a/matlab/load_last_mh_history_file.m b/matlab/load_last_mh_history_file.m
index 3e54f9853e44802b97e78efc3981e7e88feb7ca8..7dfcbb0effa0cd3b02790ed93e5fce4069bbd07a 100644
--- a/matlab/load_last_mh_history_file.m
+++ b/matlab/load_last_mh_history_file.m
@@ -10,7 +10,7 @@ function info = load_last_mh_history_file(MetropolisFolder, ModelName)
% Notes: The record structure is written to the caller workspace via an
% assignin statement.
-% Copyright (C) 2013 Dynare Team
+% Copyright (C) 2013-16 Dynare Team
%
% This file is part of Dynare.
%
@@ -53,6 +53,7 @@ if format_mh_history_file %needed to preserve backward compatibility
record.AcceptanceRatio = record.AcceptationRates;
record.InitialSeeds = NaN; % This information is forever lost...
record.MCMCConcludedSuccessfully = NaN; % This information is forever lost...
+ record.MAX_nruns=NaN(size(record.MhDraws,1),1); % This information is forever lost...
record = rmfield(record,'LastLogLiK');
record = rmfield(record,'InitialLogLiK');
record = rmfield(record,'Seeds');
@@ -64,6 +65,9 @@ else
if ~isfield(record,'MCMCConcludedSuccessfully')
record.MCMCConcludedSuccessfully = NaN; % This information is forever lost...
end
+ if ~isfield(record,'MAX_nruns')
+ record.MAX_nruns=NaN(size(record.MhDraws,1),1); % This information is forever lost...
+ end
end
if isequal(nargout,0)
diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m
index f31111d43712e16ed6088383b290ab2271869fa9..663e97db57ebca9713a3ccc26e20491c269d6ec3 100644
--- a/matlab/metropolis_hastings_initialization.m
+++ b/matlab/metropolis_hastings_initialization.m
@@ -1,4 +1,4 @@
-function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
+function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ...
metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
%function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
@@ -59,10 +59,10 @@ ix2 = [];
ilogpo2 = [];
ModelName = [];
MetropolisFolder = [];
-fblck = [];
-fline = [];
+FirstBlock = [];
+FirstLine = [];
npar = [];
-nblck = [];
+NumberOfBlocks = [];
nruns = [];
NewFile = [];
MAX_nruns = [];
@@ -76,15 +76,15 @@ end
MetropolisFolder = CheckPath('metropolis',M_.dname);
BaseName = [MetropolisFolder filesep ModelName];
-nblck = options_.mh_nblck;
-nruns = ones(nblck,1)*options_.mh_replic;
+NumberOfBlocks = options_.mh_nblck;
+nruns = ones(NumberOfBlocks,1)*options_.mh_replic;
npar = length(xparam1);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
d = chol(vv);
if ~options_.load_mh_file && ~options_.mh_recover
% Here we start a new Metropolis-Hastings, previous draws are discarded.
- if nblck > 1
+ if NumberOfBlocks > 1
disp('Estimation::mcmc: Multiple chains mode.')
else
disp('Estimation::mcmc: One Chain mode.')
@@ -108,7 +108,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,' \n\n');
fprintf(fidlog,'%% Session 1.\n');
fprintf(fidlog,' \n');
- fprintf(fidlog,[' Number of blocks...............: ' int2str(nblck) '\n']);
+ fprintf(fidlog,[' Number of blocks...............: ' int2str(NumberOfBlocks) '\n']);
fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']);
fprintf(fidlog,' \n');
% Find initial values for the nblck chains...
@@ -116,9 +116,9 @@ if ~options_.load_mh_file && ~options_.mh_recover
set_dynare_seed('default');
fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::mcmc: Searching for initial values...')
- ix2 = zeros(nblck,npar);
- ilogpo2 = zeros(nblck,1);
- for j=1:nblck
+ ix2 = zeros(NumberOfBlocks,npar);
+ ilogpo2 = zeros(NumberOfBlocks,1);
+ for j=1:NumberOfBlocks
validate = 0;
init_iter = 0;
trial = 1;
@@ -184,22 +184,23 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,' \n');
end
fprintf(fidlog,' \n');
- fblck = 1;
- fline = ones(nblck,1);
- NewFile = ones(nblck,1);
+ FirstBlock = 1;
+ FirstLine = ones(NumberOfBlocks,1);
+ NewFile = ones(NumberOfBlocks,1);
% Delete the mh-history files
delete_mh_history_files(MetropolisFolder, ModelName);
% Create a new record structure
fprintf(['Estimation::mcmc: Write details about the MCMC... ']);
AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
- record.Nblck = nblck;
+ record.Nblck = NumberOfBlocks;
record.MhDraws = zeros(1,3);
record.MhDraws(1,1) = nruns(1);
record.MhDraws(1,2) = AnticipatedNumberOfFiles;
record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
- record.AcceptanceRatio = zeros(1,nblck);
- for j=1:nblck
+ record.MAX_nruns=MAX_nruns;
+ record.AcceptanceRatio = zeros(1,NumberOfBlocks);
+ for j=1:NumberOfBlocks
% we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
set_dynare_seed(options_.DynareRandomStreams.seed+j);
% record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
@@ -207,8 +208,8 @@ if ~options_.load_mh_file && ~options_.mh_recover
end
record.InitialParameters = ix2;
record.InitialLogPost = ilogpo2;
- record.LastParameters = zeros(nblck,npar);
- record.LastLogPost = zeros(nblck,1);
+ record.LastParameters = zeros(NumberOfBlocks,npar);
+ record.LastLogPost = zeros(NumberOfBlocks,1);
record.LastFileNumber = AnticipatedNumberOfFiles ;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
record.MCMCConcludedSuccessfully = 0;
@@ -220,7 +221,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,[' Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']);
fprintf(fidlog,[' Expected number of lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']);
fprintf(fidlog,['\n']);
- for j = 1:nblck,
+ for j = 1:NumberOfBlocks,
fprintf(fidlog,[' Initial state of the Gaussian random number generator for chain number ',int2str(j),':\n']);
for i=1:length(record.InitialSeeds(j).Normal)
fprintf(fidlog,[' ' num2str(record.InitialSeeds(j).Normal(i)') '\n']);
@@ -248,30 +249,30 @@ elseif options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,'\n');
fprintf(fidlog,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\n']);
fprintf(fidlog,' \n');
- fprintf(fidlog,[' Number of blocks...............: ' int2str(nblck) '\n']);
+ fprintf(fidlog,[' Number of blocks...............: ' int2str(NumberOfBlocks) '\n']);
fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']);
fprintf(fidlog,' \n');
past_number_of_blocks = record.Nblck;
- if past_number_of_blocks ~= nblck
+ if past_number_of_blocks ~= NumberOfBlocks
disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!')
- disp(['Estimation::mcmc: You declared ' int2str(nblck) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.'])
+ disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.'])
disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
- nblck = past_number_of_blocks;
- options_.mh_nblck = nblck;
+ NumberOfBlocks = past_number_of_blocks;
+ options_.mh_nblck = NumberOfBlocks;
end
% I read the last line of the last mh-file for initialization of the new Metropolis-Hastings simulations:
LastFileNumber = record.LastFileNumber;
LastLineNumber = record.LastLineNumber;
if LastLineNumber < MAX_nruns
- NewFile = ones(nblck,1)*LastFileNumber;
- fline = ones(nblck,1)*(LastLineNumber+1);
+ NewFile = ones(NumberOfBlocks,1)*LastFileNumber;
+ FirstLine = ones(NumberOfBlocks,1)*(LastLineNumber+1);
else
- NewFile = ones(nblck,1)*(LastFileNumber+1);
- fline = ones(nblck,1);
+ NewFile = ones(NumberOfBlocks,1)*(LastFileNumber+1);
+ FirstLine = ones(NumberOfBlocks,1);
end
ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
- fblck = 1;
+ FirstBlock = 1;
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
fprintf('Estimation::mcmc: I am writing a new mh-history file... ');
record.MhDraws = [record.MhDraws;zeros(1,3)];
@@ -284,6 +285,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
record.MhDraws(end,1) = nruns(1);
record.MhDraws(end,2) = AnticipatedNumberOfFiles;
record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
+ record.MAX_nruns = [record.MAX_nruns;MAX_nruns];
record.InitialSeeds = record.LastSeeds;
write_mh_history_file(MetropolisFolder, ModelName, record);
fprintf('Done.\n')
@@ -294,15 +296,32 @@ elseif options_.mh_recover
% The previous metropolis-hastings crashed before the end! I try to recover the saved draws...
disp('Estimation::mcmc: Recover mode!')
load_last_mh_history_file(MetropolisFolder, ModelName);
- nblck = record.Nblck;% Number of "parallel" mcmc chains.
- options_.mh_nblck = nblck;
+ NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains.
+ options_.mh_nblck = NumberOfBlocks;
+
+ %% check consistency of options
+ if record.MhDraws(end,1)~=options_.mh_replic
+ fprintf('\nEstimation::mcmc: You cannot specify a different mh_replic than in the chain you are trying to recover\n')
+ fprintf('Estimation::mcmc: I am resetting mh_replic to %u\n',record.MhDraws(end,1))
+ options_.mh_replic=record.MhDraws(end,1);
+ nruns = ones(NumberOfBlocks,1)*options_.mh_replic;
+ end
+
+ if ~isnan(record.MAX_nruns(end,1)) %field exists
+ if record.MAX_nruns(end,1)~=MAX_nruns
+ fprintf('\nEstimation::mcmc: You cannot specify a different MaxNumberOfBytes than in the chain you are trying to recover\n')
+ fprintf('Estimation::mcmc: I am resetting MAX_nruns to %u\n',record.MAX_nruns(end,1))
+ MAX_nruns=record.MAX_nruns(end,1);
+ end
+ end
+ %% do tentative initialization as if full last run of MCMC were to be redone
if size(record.MhDraws,1) == 1
- OldMh = 0;% The crashed metropolis was the first session.
+ OldMhExists = 0;% The crashed metropolis was the first session.
else
- OldMh = 1;% The crashed metropolis wasn't the first session.
+ OldMhExists = 1;% The crashed metropolis wasn't the first session.
end
% Default initialization:
- if OldMh
+ if OldMhExists
ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
else
@@ -310,67 +329,110 @@ elseif options_.mh_recover
ix2 = record.InitialParameters;
end
% Set NewFile, a nblck*1 vector of integers, and fline (first line), a nblck*1 vector of integers.
- if OldMh
+ % Relevant for starting at the correct file and potentially filling a file from the previous run that was not yet full
+ if OldMhExists
LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);% Number of lines in the last mh files of the previous session.
LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);% Number of mh files in the the previous sessions.
- if LastLineNumberInThePreviousMh < MAX_nruns% Test if the last mh files of the previous session are not complete (yes)
- NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh;
- fline = ones(nblck,1)*(LastLineNumberInThePreviousMh+1);
- else% The last mh files of the previous session are complete.
- NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1);
- fline = ones(nblck,1);
+ %Test if the last mh files of the previous session were not full yet
+ if LastLineNumberInThePreviousMh < MAX_nruns%not full
+ %store starting point if whole chain needs to be redone
+ NewFile = ones(NumberOfBlocks,1)*LastFileNumberInThePreviousMh;
+ FirstLine = ones(NumberOfBlocks,1)*(LastLineNumberInThePreviousMh+1);
+ LastFileFullIndicator=0;
+ else% The last mh files of the previous session were full
+ %store starting point if whole chain needs to be redone
+ NewFile = ones(NumberOfBlocks,1)*(LastFileNumberInThePreviousMh+1);
+ FirstLine = ones(NumberOfBlocks,1);
+ LastFileFullIndicator=1;
end
else
LastLineNumberInThePreviousMh = 0;
LastFileNumberInThePreviousMh = 0;
- NewFile = ones(nblck,1);
- fline = ones(nblck,1);
+ NewFile = ones(NumberOfBlocks,1);
+ FirstLine = ones(NumberOfBlocks,1);
+ LastFileFullIndicator=1;
end
- % Set fblck (First block), an integer targeting the crashed mcmc chain.
- fblck = 1;
+
+ %% 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 ?
ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
- ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck;
- % I count the total number of saved mh files...
+ ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*NumberOfBlocks;
+ % How many mh files do we actually have ?
AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
TotalNumberOfMhFiles = length(AllMhFiles);
- % And I quit if I can't find a crashed mcmc chain.
+ % Quit if no crashed mcmc chain can be found as there are as many files as expected
if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!')
disp(' You have to edit the mod file and remove the mh_recover option')
disp(' in the estimation command')
error('Estimation::mcmc: mh_recover option not required!')
end
- % I count the number of saved mh files per block.
- NumberOfMhFilesPerBlock = zeros(nblck,1);
- for b = 1:nblck
+ % 2. Something needs to be done; find out what
+ % Count the number of saved mh files per block.
+ NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
+ for b = 1:NumberOfBlocks
BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
NumberOfMhFilesPerBlock(b) = length(BlckMhFiles);
end
- % Is there a chain with less mh files than expected ?
- while fblck <= nblck
- if NumberOfMhFilesPerBlock(fblck) < ExpectedNumberOfMhFilesPerBlock
- disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is not complete!'])
+ % Find fblck (First block), an integer targeting the crashed mcmc chain.
+ FirstBlock = 1; %initialize
+ while FirstBlock <= NumberOfBlocks
+ if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock
+ disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is not complete!'])
break
% The mh_recover session will start from chain fblck.
else
- disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is complete!'])
+ disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is complete!'])
end
- fblck = fblck+1;
+ FirstBlock = FirstBlock+1;
end
+
+ %% 3. Overwrite default settings for
% How many mh-files are saved in this block?
- NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck);
- % Correct the number of saved mh files if the crashed Metropolis was not the first session (so
- % that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain
- % of the current session).
- if OldMh
- NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh;
+ NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(FirstBlock);
+ ExistingDrawsInLastMCFile=0; %initialize: no MCMC draws of current MCMC are in file from last run
+ % Check whether last present file is a file included in the last MCMC run
+ if ~LastFileFullIndicator
+ if NumberOfSavedMhFilesInTheCrashedBlck==NewFile(FirstBlock) %only that last file exists, but no files from current MCMC
+ loaded_results=load([BaseName '_mh' int2str(NewFile(FirstBlock)) '_blck' int2str(FirstBlock) '.mat']);
+ %check whether that last file was filled
+ if size(loaded_results.x2,1)==MAX_nruns %file is full
+ NewFile(FirstBlock)=NewFile(FirstBlock)+1; %set first file to be created to next one
+ FirstLine(FirstBlock) = 1; %use first line of next file
+ ExistingDrawsInLastMCFile=MAX_nruns-record.MhDraws(end-1,3);
+ else
+ ExistingDrawsInLastMCFile=0;
+ end
+ end
+ elseif LastFileFullIndicator
+ ExistingDrawsInLastMCFile=0;
+ if NumberOfSavedMhFilesInTheCrashedBlck==NewFile(FirstBlock) %only the last file exists, but no files from current MCMC
+ NewFile(FirstBlock)=NewFile(FirstBlock)+1; %set first file to be created to next one
+ end
end
- NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh;
+% % Correct the number of saved mh files if the crashed Metropolis was not the first session (so
+% % that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain
+% % of the current session).
+% if OldMhExists
+% NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh;
+% end
+% NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh;
+
% Correct initial conditions.
- if NumberOfSavedMhFiles
- load([BaseName '_mh' int2str(NumberOfSavedMhFiles) '_blck' int2str(fblck) '.mat']);
- ilogpo2(fblck) = logpo2(end);
- ix2(fblck,:) = x2(end,:);
+ if NumberOfSavedMhFilesInTheCrashedBlck<ExpectedNumberOfMhFilesPerBlock
+ loaded_results=load([BaseName '_mh' int2str(NumberOfSavedMhFilesInTheCrashedBlck) '_blck' int2str(FirstBlock) '.mat']);
+ ilogpo2(FirstBlock) = loaded_results.logpo2(end);
+ ix2(FirstBlock,:) = loaded_results.x2(end,:);
+ nruns(FirstBlock)=nruns(FirstBlock)-ExistingDrawsInLastMCFile-(NumberOfSavedMhFilesInTheCrashedBlck-LastFileNumberInThePreviousMh)*MAX_nruns;
+ %reset seed if possible
+ if isfield(loaded_results,'LastSeeds')
+ record.InitialSeeds(FirstBlock).Unifor=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Unifor;
+ record.InitialSeeds(FirstBlock).Normal=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Normal;
+ else
+ fprintf('Estimation::mcmc: You are trying to recover a chain generated with an older Dynare version.\n')
+ fprintf('Estimation::mcmc: I am using the default seeds to continue the chain.\n')
+ end
+ write_mh_history_file(MetropolisFolder, ModelName, record);
end
end
\ No newline at end of file
diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m
index 2c26de0e5c5cb53ce1ff5adb9041ca35390ecad4..83bbb48aa47615ad5596fb7b38526fba09dc8f94 100644
--- a/matlab/random_walk_metropolis_hastings_core.m
+++ b/matlab/random_walk_metropolis_hastings_core.m
@@ -109,6 +109,7 @@ proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
block_iter=0;
for curr_block = fblck:nblck,
+ LastSeeds=[];
block_iter=block_iter+1;
try
% This will not work if the master uses a random number generator not
@@ -171,7 +172,8 @@ for curr_block = fblck:nblck,
dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]);
end
if (draw_index_current_file == InitSizeArray(curr_block)) || (draw_iter == nruns(curr_block)) % Now I save the simulations, either because the current file is full or the chain is done
- save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2');
+ [LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal] = get_dynare_random_generator_state();
+ save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2','LastSeeds');
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 00f5879643364ce8e217f1d635d667de9541af5e..42f8c29efebcf5646331b23e189cbd385e899735 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -7,6 +7,9 @@ MODFILES = \
estimation/fs2000_model_comparison.mod \
estimation/fs2000_fast.mod \
estimation/MH_recover/fs2000_recover.mod \
+ estimation/MH_recover/fs2000_recover_2.mod \
+ estimation/MH_recover/fs2000_recover_3.mod \
+ estimation/MH_recover/fs2000_recover_tarb.mod \
estimation/t_proposal/fs2000_student.mod \
estimation/TaRB/fs2000_tarb.mod \
moments/example1_var_decomp.mod \
@@ -464,6 +467,7 @@ EXTRA_DIST = \
optimal_policy/Ramsey/oo_ramsey_policy_initval.mat \
optimizers/optimizer_function_wrapper.m \
optimizers/fs2000.common.inc \
+ estimation/MH_recover/fs2000.common.inc \
prior_posterior_function/posterior_function_demo.m
diff --git a/tests/estimation/MH_recover/fs2000.common.inc b/tests/estimation/MH_recover/fs2000.common.inc
new file mode 100644
index 0000000000000000000000000000000000000000..bbcc293f9ba9ecb9991b51c7f94893381a0481d9
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000.common.inc
@@ -0,0 +1,116 @@
+/*
+ * This file replicates the estimation of the cash in advance model described
+ * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
+ * Journal of Applied Econometrics, 15(6), 645-670.
+ *
+ * The data are in file "fsdat_simul.m", and have been artificially generated.
+ * They are therefore different from the original dataset used by Schorfheide.
+ *
+ * The equations are taken from J. Nason and T. Cogley (1994): "Testing the
+ * implications of long-run neutrality for monetary business cycle models",
+ * Journal of Applied Econometrics, 9, S37-S70.
+ * Note that there is an initial minus sign missing in equation (A1), p. S63.
+ *
+ * This implementation was written by Michel Juillard. Please note that the
+ * following copyright notice only applies to this Dynare implementation of the
+ * model.
+ */
+
+/*
+ * Copyright (C) 2004-2010 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 <http://www.gnu.org/licenses/>.
+ */
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+ W = l/n;
+ ist = y-c;
+ q = 1 - d;
+
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+options_.plot_priors=0;
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover.mod b/tests/estimation/MH_recover/fs2000_recover.mod
index 52c5664ac8e20fa9c233f21ec81f2692b48937c1..3f421a8ee2d0f8bccf0df4924065ed3f0d982eb9 100644
--- a/tests/estimation/MH_recover/fs2000_recover.mod
+++ b/tests/estimation/MH_recover/fs2000_recover.mod
@@ -1,139 +1,35 @@
-/*
- * This file replicates the estimation of the cash in advance model described
- * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
- * Journal of Applied Econometrics, 15(6), 645-670.
- *
- * The data are in file "fsdat_simul.m", and have been artificially generated.
- * They are therefore different from the original dataset used by Schorfheide.
- *
- * The equations are taken from J. Nason and T. Cogley (1994): "Testing the
- * implications of long-run neutrality for monetary business cycle models",
- * Journal of Applied Econometrics, 9, S37-S70.
- * Note that there is an initial minus sign missing in equation (A1), p. S63.
- *
- * This implementation was written by Michel Juillard. Please note that the
- * following copyright notice only applies to this Dynare implementation of the
- * model.
- */
+//Test mh_recover function for RW-MH
-/*
- * Copyright (C) 2004-2010 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 <http://www.gnu.org/licenses/>.
- */
+@#include "fs2000.common.inc"
-var m P c e W R k d n l gy_obs gp_obs y dA;
-varexo e_a e_m;
-
-parameters alp bet gam mst rho psi del;
-
-alp = 0.33;
-bet = 0.99;
-gam = 0.003;
-mst = 1.011;
-rho = 0.7;
-psi = 0.787;
-del = 0.02;
-
-model;
-dA = exp(gam+e_a);
-log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
--P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
-W = l/n;
--(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
-R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
-1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
-c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
-P*c = m;
-m-1+d = l;
-e = exp(e_a);
-y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
-gy_obs = dA*y/y(-1);
-gp_obs = (P/P(-1))*m(-1)/dA;
-end;
-
-shocks;
-var e_a; stderr 0.014;
-var e_m; stderr 0.005;
-end;
-
-steady_state_model;
- dA = exp(gam);
- gst = 1/dA;
- m = mst;
- khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
- xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
- nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
- n = xist/(nust+xist);
- P = xist + nust;
- k = khst*n;
-
- l = psi*mst*n/( (1-psi)*(1-n) );
- c = mst/P;
- d = l - mst + 1;
- y = k^alp*n^(1-alp)*gst^alp;
- R = mst/bet;
- W = l/n;
- ist = y-c;
- q = 1 - d;
-
- e = 1;
-
- gp_obs = m/dA;
- gy_obs = dA;
-end;
-
-steady;
-
-check;
-
-estimated_params;
-alp, beta_pdf, 0.356, 0.02;
-bet, beta_pdf, 0.993, 0.002;
-gam, normal_pdf, 0.0085, 0.003;
-mst, normal_pdf, 1.0002, 0.007;
-rho, beta_pdf, 0.129, 0.223;
-psi, beta_pdf, 0.65, 0.05;
-del, beta_pdf, 0.01, 0.005;
-stderr e_a, inv_gamma_pdf, 0.035449, inf;
-stderr e_m, inv_gamma_pdf, 0.008862, inf;
-end;
-
-varobs gp_obs gy_obs;
-
-options_.MaxNumberOfBytes=2000*11*8/4;
-estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8);
-copyfile([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh1_blck1.mat'],'fs2000_mh1_blck1.mat')
-copyfile([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh3_blck2.mat'],'fs2000_mh3_blck2.mat')
-delete([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh4_blck2.mat'])
+options_.MaxNumberOfBytes=1000*11*8/2;
+estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=1000, mh_nblocks=2, mh_jscale=0.8);
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'],[M_.dname '_mh2_blck2.mat'])
+delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'])
estimation(order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
%check first unaffected chain
-temp1=load('fs2000_mh1_blck1.mat');
-temp2=load([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh1_blck1.mat']);
+temp1=load([M_.dname '_mh1_blck1.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
-if max(max(temp1.x2-temp2.x2))>1e-10
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of unaffected chain are not the same')
end
%check second, affected chain with last unaffected file
-temp1=load('fs2000_mh3_blck2.mat');
-temp2=load([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh3_blck2.mat']);
+temp1=load([M_.dname '_mh1_blck1.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
-if max(max(temp1.x2-temp2.x2))>1e-10
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s unaffected files are not the same')
end
+
+%check second, affected chain with affected file
+temp1=load([M_.dname '_mh2_blck2.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of affected chain''s affected files are not the same')
+end
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover_2.mod b/tests/estimation/MH_recover/fs2000_recover_2.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1d503de7c6dd29d4b3084bd41ce9f0b30cf5bea2
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000_recover_2.mod
@@ -0,0 +1,49 @@
+//Test mh_recover function for RW-MH when load_mh_file was also used
+//Previous chain did not fill last mh_file so mh_recover needs to fill that file
+
+@#include "fs2000.common.inc"
+
+options_.MaxNumberOfBytes=2000*11*8/4;
+estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=999, mh_nblocks=2, mh_jscale=0.8);
+estimation(order=1,mode_compute=0,mode_file=fs2000_recover_2_mode, datafile='../fsdat_simul',nobs=192, loglinear, load_mh_file,mh_replic=1002, mh_nblocks=2, mh_jscale=0.8);
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'],[M_.dname '_mh3_blck2.mat'])
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'],[M_.dname '_mh4_blck2.mat'])
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh5_blck2.mat'],[M_.dname '_mh5_blck2.mat'])
+delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'])
+delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'])
+delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh5_blck2.mat'])
+
+estimation(order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_2_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
+
+%check first unaffected chain
+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
+
+%check second, affected chain with last unaffected file
+temp1=load([M_.dname '_mh3_blck2.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of affected chain''s unaffected files are not the same')
+end
+
+%check second, affected chain with affected file
+temp1=load([M_.dname '_mh4_blck2.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of affected chain''s affected files are not the same')
+end
+
+%check second, affected chain with affected file
+temp1=load([M_.dname '_mh5_blck2.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh5_blck2.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of affected chain''s affected files are not the same')
+end
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover_3.mod b/tests/estimation/MH_recover/fs2000_recover_3.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3f7b920ec3346ba2720578d81928047b56b868b4
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000_recover_3.mod
@@ -0,0 +1,39 @@
+//Test mh_recover function for RW-MH when load_mh_file was also used
+//Previous chain did fill last mh_file so mh_recover does not need to fill that file
+
+@#include "fs2000.common.inc"
+
+options_.MaxNumberOfBytes=2000*11*8/4;
+estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=1000, mh_nblocks=2, mh_jscale=0.8);
+estimation(order=1,mode_compute=0,mode_file=fs2000_recover_3_mode, datafile='../fsdat_simul',nobs=192, loglinear, load_mh_file,mh_replic=1000, mh_nblocks=2, mh_jscale=0.8);
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'],[M_.dname '_mh3_blck2.mat'])
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'],[M_.dname '_mh4_blck2.mat'])
+delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'])
+delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'])
+
+estimation(order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_3_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
+
+%check first unaffected chain
+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
+
+%check second, affected chain with last unaffected file
+temp1=load([M_.dname '_mh3_blck2.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of affected chain''s unaffected files are not the same')
+end
+
+%check second, affected chain with affected file
+temp1=load([M_.dname '_mh4_blck2.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of affected chain''s affected files are not the same')
+end
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover_tarb.mod b/tests/estimation/MH_recover/fs2000_recover_tarb.mod
new file mode 100644
index 0000000000000000000000000000000000000000..9c4fca92a4feacc1c1c0f94963cebb78aedd6e23
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000_recover_tarb.mod
@@ -0,0 +1,35 @@
+//Test mh_recover function for RW-MH when use_tarb was also used
+
+@#include "fs2000.common.inc"
+
+options_.MaxNumberOfBytes=10*11*8/2;
+estimation(use_tarb,order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=10, mh_nblocks=2, mh_jscale=0.8);
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
+copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'],[M_.dname '_mh2_blck2.mat'])
+delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'])
+
+estimation(use_tarb,order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
+
+%check first unaffected chain
+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
+
+%check second, affected chain with last unaffected 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 affected chain''s unaffected files are not the same')
+end
+
+%check second, affected chain with affected file
+temp1=load([M_.dname '_mh2_blck2.mat']);
+temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat']);
+
+if max(max(abs(temp1.x2-temp2.x2)))>1e-10
+ error('Draws of affected chain''s affected files are not the same')
+end
\ No newline at end of file