Manual Merge of branch 'mh_recover' (Johannes)
commit
a55f00073b
|
@ -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{}
|
||||
|
|
|
@ -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']);
|
||||
|
|
|
@ -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)
|
||||
|
|
|
@ -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
|
|
@ -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']);
|
||||
|
|
|
@ -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
|
||||
|
||||
|
||||
|
|
|
@ -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;
|
|
@ -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
|
|
@ -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
|
|
@ -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
|
|
@ -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
|
Loading…
Reference in New Issue