Merge pull request #911 from JohannesPfeifer/MH_cosmetics
Cosmetic Fixes to Metropolis-Hastings routinestime-shift
commit
ddc7464b9b
|
@ -1,7 +1,7 @@
|
||||||
function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
|
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_)
|
metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||||
%function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] =
|
%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_)
|
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||||
% Metropolis-Hastings initialization.
|
% Metropolis-Hastings initialization.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
|
@ -11,6 +11,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck,
|
||||||
% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode).
|
% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode).
|
||||||
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
|
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
|
||||||
% o dataset_ data structure
|
% o dataset_ data structure
|
||||||
|
% o dataset_info dataset info structure
|
||||||
% o options_ options structure
|
% o options_ options structure
|
||||||
% o M_ model structure
|
% o M_ model structure
|
||||||
% o estim_params_ estimated parameters structure
|
% o estim_params_ estimated parameters structure
|
||||||
|
@ -18,12 +19,25 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck,
|
||||||
% o oo_ outputs structure
|
% o oo_ outputs structure
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% None
|
% o ix2 [double] (nblck*npar) vector of starting points for different chains
|
||||||
|
% o ilogpo2 [double] (nblck*1) vector of initial posterior values for different chains
|
||||||
|
% o ModelName [string] name of the mod-file
|
||||||
|
% o MetropolisFolder [string] path to the Metropolis subfolder
|
||||||
|
% o fblck [scalar] number of the first MH chain to be run (not equal to 1 in case of recovery)
|
||||||
|
% o fline [double] (nblck*1) vector of first draw in each chain (not equal to 1 in case of recovery)
|
||||||
|
% o npar [scalar] number of parameters estimated
|
||||||
|
% o nblck [scalar] Number of MCM chains requested
|
||||||
|
% o nruns [double] (nblck*1) number of draws in each chain
|
||||||
|
% o NewFile [scalar] (nblck*1) vector storing the number
|
||||||
|
% of the first MH-file to created for each chain when saving additional
|
||||||
|
% 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).
|
||||||
%
|
%
|
||||||
% SPECIAL REQUIREMENTS
|
% SPECIAL REQUIREMENTS
|
||||||
% None.
|
% None.
|
||||||
|
|
||||||
% Copyright (C) 2006-2013 Dynare Team
|
% Copyright (C) 2006-2015 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -40,6 +54,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck,
|
||||||
% You should have received a copy of the GNU General Public License
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
%Initialize outputs
|
||||||
ix2 = [];
|
ix2 = [];
|
||||||
ilogpo2 = [];
|
ilogpo2 = [];
|
||||||
ModelName = [];
|
ModelName = [];
|
||||||
|
@ -68,7 +83,7 @@ MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
|
||||||
d = chol(vv);
|
d = chol(vv);
|
||||||
|
|
||||||
if ~options_.load_mh_file && ~options_.mh_recover
|
if ~options_.load_mh_file && ~options_.mh_recover
|
||||||
% Here we start a new metropolis-hastings, previous draws are discarded.
|
% Here we start a new Metropolis-Hastings, previous draws are discarded.
|
||||||
if nblck > 1
|
if nblck > 1
|
||||||
disp('Estimation::mcmc: Multiple chains mode.')
|
disp('Estimation::mcmc: Multiple chains mode.')
|
||||||
else
|
else
|
||||||
|
@ -80,7 +95,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
|
||||||
delete([BaseName '_mh*_blck*.mat']);
|
delete([BaseName '_mh*_blck*.mat']);
|
||||||
disp('Estimation::mcmc: Old mh-files successfully erased!')
|
disp('Estimation::mcmc: Old mh-files successfully erased!')
|
||||||
end
|
end
|
||||||
% Delete old metropolis log file.
|
% Delete old Metropolis log file.
|
||||||
file = dir([ MetropolisFolder '/metropolis.log']);
|
file = dir([ MetropolisFolder '/metropolis.log']);
|
||||||
if length(file)
|
if length(file)
|
||||||
delete([ MetropolisFolder '/metropolis.log']);
|
delete([ MetropolisFolder '/metropolis.log']);
|
||||||
|
@ -128,11 +143,11 @@ if ~options_.load_mh_file && ~options_.mh_recover
|
||||||
if init_iter > 100 && validate == 0
|
if init_iter > 100 && validate == 0
|
||||||
disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.'])
|
disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.'])
|
||||||
if options_.nointeractive
|
if options_.nointeractive
|
||||||
disp(['Estimation::mcmc: I reduce mh_init_scale by ten percent:'])
|
disp(['Estimation::mcmc: I reduce mh_init_scale by 10 percent:'])
|
||||||
options_.mh_init_scale = .9*options_.mh_init_scale;
|
options_.mh_init_scale = .9*options_.mh_init_scale;
|
||||||
disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.',options_.mh_init_scale))
|
disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.',options_.mh_init_scale))
|
||||||
else
|
else
|
||||||
disp(['Estimation::mcmc: You should Reduce mh_init_scale...'])
|
disp(['Estimation::mcmc: You should reduce mh_init_scale...'])
|
||||||
disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale))
|
disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale))
|
||||||
options_.mh_init_scale = input('Estimation::mcmc: Enter a new value... ');
|
options_.mh_init_scale = input('Estimation::mcmc: Enter a new value... ');
|
||||||
end
|
end
|
||||||
|
@ -217,7 +232,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
|
||||||
fclose(fidlog);
|
fclose(fidlog);
|
||||||
elseif options_.load_mh_file && ~options_.mh_recover
|
elseif options_.load_mh_file && ~options_.mh_recover
|
||||||
% Here we consider previous mh files (previous mh did not crash).
|
% Here we consider previous mh files (previous mh did not crash).
|
||||||
disp('Estimation::mcmc: I am loading past metropolis-hastings simulations...')
|
disp('Estimation::mcmc: I am loading past Metropolis-Hastings simulations...')
|
||||||
load_last_mh_history_file(MetropolisFolder, ModelName);
|
load_last_mh_history_file(MetropolisFolder, ModelName);
|
||||||
mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']);
|
mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']);
|
||||||
if ~length(mh_files)
|
if ~length(mh_files)
|
||||||
|
@ -238,7 +253,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
|
||||||
nblck = past_number_of_blocks;
|
nblck = past_number_of_blocks;
|
||||||
options_.mh_nblck = nblck;
|
options_.mh_nblck = nblck;
|
||||||
end
|
end
|
||||||
% I read the last line of the last mh-file for initialization of the new metropolis-hastings simulations:
|
% I read the last line of the last mh-file for initialization of the new Metropolis-Hastings simulations:
|
||||||
LastFileNumber = record.LastFileNumber;
|
LastFileNumber = record.LastFileNumber;
|
||||||
LastLineNumber = record.LastLineNumber;
|
LastLineNumber = record.LastLineNumber;
|
||||||
if LastLineNumber < MAX_nruns
|
if LastLineNumber < MAX_nruns
|
||||||
|
@ -252,7 +267,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
|
||||||
ix2 = record.LastParameters;
|
ix2 = record.LastParameters;
|
||||||
fblck = 1;
|
fblck = 1;
|
||||||
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
|
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
|
||||||
fprintf('Estimation::mcmc: I am writting a new mh-history file... ');
|
fprintf('Estimation::mcmc: I am writing a new mh-history file... ');
|
||||||
record.MhDraws = [record.MhDraws;zeros(1,3)];
|
record.MhDraws = [record.MhDraws;zeros(1,3)];
|
||||||
NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber;
|
NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber;
|
||||||
NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile;
|
NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile;
|
||||||
|
@ -318,7 +333,7 @@ elseif options_.mh_recover
|
||||||
disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!')
|
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(' You have to edit the mod file and remove the mh_recover option')
|
||||||
disp(' in the estimation command')
|
disp(' in the estimation command')
|
||||||
error()
|
error('Estimation::mcmc: mh_recover option not required!')
|
||||||
end
|
end
|
||||||
% I count the number of saved mh files per block.
|
% I count the number of saved mh files per block.
|
||||||
NumberOfMhFilesPerBlock = zeros(nblck,1);
|
NumberOfMhFilesPerBlock = zeros(nblck,1);
|
||||||
|
@ -339,7 +354,7 @@ elseif options_.mh_recover
|
||||||
end
|
end
|
||||||
% How many mh-files are saved in this block?
|
% How many mh-files are saved in this block?
|
||||||
NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck);
|
NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck);
|
||||||
% Correct the number of saved mh files if the crashed metropolis was not the first session (so
|
% 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
|
% that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain
|
||||||
% of the current session).
|
% of the current session).
|
||||||
if OldMh
|
if OldMh
|
||||||
|
|
|
@ -3,28 +3,38 @@ function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock,
|
||||||
% This is the most important function for the management of DYNARE parallel
|
% This is the most important function for the management of DYNARE parallel
|
||||||
% computing.
|
% computing.
|
||||||
% It is the top-level function called on the master computer when parallelizing a task.
|
% It is the top-level function called on the master computer when parallelizing a task.
|
||||||
|
%
|
||||||
% This function has two main computational strategies for managing the matlab worker (slave process).
|
% This function has two main computational strategies for managing the
|
||||||
|
% matlab worker (slave process):
|
||||||
|
%
|
||||||
% 0 Simple Close/Open Stategy:
|
% 0 Simple Close/Open Stategy:
|
||||||
% In this case the new Matlab instances (slave process) are open when
|
% In this case the new Matlab instances (slave process) are open when
|
||||||
% necessary and then closed. This can happen many times during the
|
% necessary and then closed. This can happen many times during the
|
||||||
% simulation of a model.
|
% simulation of a model.
|
||||||
|
%
|
||||||
% 1 Always Open Strategy:
|
% 1 Always Open Strategy:
|
||||||
% In this case we have a more sophisticated management of slave processes,
|
% In this case we have a more sophisticated management of slave processes,
|
||||||
% which are no longer closed at the end of each job. The slave processes
|
% which are no longer closed at the end of each job. The slave processes
|
||||||
% wait for a new job (if it exists). If a slave does not receive a new job after a
|
% wait for a new job (if it exists). If a slave does not receive a new job after a
|
||||||
% fixed time it is destroyed. This solution removes the computational
|
% fixed time it is destroyed. This solution removes the computational
|
||||||
% time necessary to Open/Close new Matlab instances.
|
% time necessary to Open/Close new Matlab instances.
|
||||||
|
%
|
||||||
% The first (point 0) is the default Strategy
|
% The first (point 0) is the default Strategy
|
||||||
% i.e.(Parallel_info.leaveSlaveOpen=0). This value can be changed by the
|
% i.e.(Parallel_info.leaveSlaveOpen=0). This value can be changed by the
|
||||||
% user in xxx.mod file or it is changed by the programmer if it is necessary to
|
% user in xxx.mod file or it is changed by the programmer if it is necessary to
|
||||||
% reduce the overall computational time. See for example the
|
% reduce the overall computational time. See for example the
|
||||||
% prior_posterior_statistics.m.
|
% prior_posterior_statistics.m.
|
||||||
|
%
|
||||||
% The number of parallelized threads will be equal to (nBlock-fBlock+1).
|
% The number of parallelized threads will be equal to (nBlock-fBlock+1).
|
||||||
%
|
%
|
||||||
|
% Treatment of global variables:
|
||||||
|
% Global variables used within the called function (e.g.
|
||||||
|
% objective_function_penalty_base) are wrapped and passed by storing their
|
||||||
|
% values at the start of the parallel computation in a file via
|
||||||
|
% storeGlobalVars.m. This file is then loaded in the separate,
|
||||||
|
% independent slave Matlab sessions. By keeping them separate, no
|
||||||
|
% interaction via global variables can take place.
|
||||||
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% o Parallel [struct vector] copy of options_.parallel
|
% o Parallel [struct vector] copy of options_.parallel
|
||||||
% o fBlock [int] index number of the first thread
|
% o fBlock [int] index number of the first thread
|
||||||
|
@ -53,7 +63,7 @@ function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock,
|
||||||
% the number of CPUs declared in "Parallel", if
|
% the number of CPUs declared in "Parallel", if
|
||||||
% the number of required threads is lower)
|
% the number of required threads is lower)
|
||||||
|
|
||||||
% Copyright (C) 2009-2013 Dynare Team
|
% Copyright (C) 2009-2015 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -126,7 +136,7 @@ if isHybridMatlabOctave || isoctave
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
if exist('fGlobalVar') && ~isempty(fGlobalVar),
|
if exist('fGlobalVar','var') && ~isempty(fGlobalVar),
|
||||||
fInputNames = fieldnames(fGlobalVar);
|
fInputNames = fieldnames(fGlobalVar);
|
||||||
for j=1:length(fInputNames),
|
for j=1:length(fInputNames),
|
||||||
TargetVar = fGlobalVar.(fInputNames{j});
|
TargetVar = fGlobalVar.(fInputNames{j});
|
||||||
|
@ -161,7 +171,7 @@ switch Strategy
|
||||||
save([fname,'_input.mat'],'fInputVar','Parallel','-append')
|
save([fname,'_input.mat'],'fInputVar','Parallel','-append')
|
||||||
|
|
||||||
case 1
|
case 1
|
||||||
if exist('fGlobalVar'),
|
if exist('fGlobalVar','var'),
|
||||||
save(['temp_input.mat'],'fInputVar','fGlobalVar')
|
save(['temp_input.mat'],'fInputVar','fGlobalVar')
|
||||||
else
|
else
|
||||||
save(['temp_input.mat'],'fInputVar')
|
save(['temp_input.mat'],'fInputVar')
|
||||||
|
@ -540,7 +550,6 @@ else
|
||||||
'Renderer','Painters', ...
|
'Renderer','Painters', ...
|
||||||
'Resize','off');
|
'Resize','off');
|
||||||
|
|
||||||
vspace = 0.1;
|
|
||||||
ncol = ceil(totCPU/10);
|
ncol = ceil(totCPU/10);
|
||||||
hspace = 0.9/ncol;
|
hspace = 0.9/ncol;
|
||||||
hstatus(1) = axes('position',[0.05/ncol 0.92 0.9/ncol 0.03], ...
|
hstatus(1) = axes('position',[0.05/ncol 0.92 0.9/ncol 0.03], ...
|
||||||
|
@ -882,8 +891,4 @@ switch Strategy
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -14,11 +14,11 @@ function prior_posterior_statistics(type,dataset,dataset_info)
|
||||||
%
|
%
|
||||||
% SPECIAL REQUIREMENTS
|
% SPECIAL REQUIREMENTS
|
||||||
% none
|
% none
|
||||||
|
%
|
||||||
% PARALLEL CONTEXT
|
% PARALLEL CONTEXT
|
||||||
% See the comments random_walk_metropolis_hastings.m funtion.
|
% See the comments in the random_walk_metropolis_hastings.m funtion.
|
||||||
|
%
|
||||||
|
%
|
||||||
% Copyright (C) 2005-2013 Dynare Team
|
% Copyright (C) 2005-2013 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
|
|
|
@ -1,14 +1,17 @@
|
||||||
function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||||
%function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
|
% function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||||
% Random walk Metropolis-Hastings algorithm.
|
% Random Walk Metropolis-Hastings algorithm.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% o TargetFun [char] string specifying the name of the objective
|
% o TargetFun [char] string specifying the name of the objective
|
||||||
% function (posterior kernel).
|
% function (posterior kernel).
|
||||||
|
% o ProposalFun [char] string specifying the name of the proposal
|
||||||
|
% density
|
||||||
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
|
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
|
||||||
% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode).
|
% o vv [double] (p*p) matrix, posterior covariance matrix (at the mode).
|
||||||
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
|
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
|
||||||
% o dataset_ data structure
|
% o dataset_ data structure
|
||||||
|
% o dataset_info dataset info structure
|
||||||
% o options_ options structure
|
% o options_ options structure
|
||||||
% o M_ model structure
|
% o M_ model structure
|
||||||
% o estim_params_ estimated parameters structure
|
% o estim_params_ estimated parameters structure
|
||||||
|
@ -16,27 +19,27 @@ function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou
|
||||||
% o oo_ outputs structure
|
% o oo_ outputs structure
|
||||||
%
|
%
|
||||||
% ALGORITHM
|
% ALGORITHM
|
||||||
% Metropolis-Hastings.
|
% Random-Walk Metropolis-Hastings.
|
||||||
%
|
%
|
||||||
% SPECIAL REQUIREMENTS
|
% SPECIAL REQUIREMENTS
|
||||||
% None.
|
% None.
|
||||||
%
|
%
|
||||||
% PARALLEL CONTEXT
|
% PARALLEL CONTEXT
|
||||||
% The most computationally intensive part of this function may be executed
|
% The most computationally intensive part of this function may be executed
|
||||||
% in parallel. The code sutable to be executed in
|
% in parallel. The code suitable to be executed in
|
||||||
% parallel on multi core or cluster machine (in general a 'for' cycle),
|
% parallel on multi core or cluster machine (in general a 'for' cycle)
|
||||||
% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
|
% has been removed from this function and been placed in the random_walk_metropolis_hastings_core.m funtion.
|
||||||
% Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in
|
%
|
||||||
% parallel and called name_function.m and name_function_core.m.
|
% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
|
||||||
% In addition in parallel package we have second set of functions used
|
% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
|
||||||
% to manage the parallel computation.
|
% to manage the parallel computations.
|
||||||
%
|
%
|
||||||
% This function was the first function to be parallelized, later other
|
% This function was the first function to be parallelized. Later, other
|
||||||
% functions have been parallelized using the same methodology.
|
% functions have been parallelized using the same methodology.
|
||||||
% Then the comments write here can be used for all the other pairs of
|
% Then the comments write here can be used for all the other pairs of
|
||||||
% parallel functions and also for management funtions.
|
% parallel functions and also for management functions.
|
||||||
|
|
||||||
% Copyright (C) 2006-2013 Dynare Team
|
% Copyright (C) 2006-2015 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -54,7 +57,7 @@ function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou
|
||||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
|
||||||
% In Metropolis, we set penalty to Inf to as to reject all parameter sets triggering error in target density computation
|
% In Metropolis, we set penalty to Inf so as to reject all parameter sets triggering an error during target density computation
|
||||||
global objective_function_penalty_base
|
global objective_function_penalty_base
|
||||||
objective_function_penalty_base = Inf;
|
objective_function_penalty_base = Inf;
|
||||||
|
|
||||||
|
@ -73,16 +76,16 @@ load_last_mh_history_file(MetropolisFolder, ModelName);
|
||||||
% First run in serial mode, and then comment the follow line.
|
% First run in serial mode, and then comment the follow line.
|
||||||
% save('recordSerial.mat','-struct', 'record');
|
% save('recordSerial.mat','-struct', 'record');
|
||||||
|
|
||||||
% For parllel runs after serial runs with the abobe line active.
|
% For parallel runs after serial runs with the abobe line active.
|
||||||
% TempRecord=load('recordSerial.mat');
|
% TempRecord=load('recordSerial.mat');
|
||||||
% record.Seeds=TempRecord.Seeds;
|
% record.Seeds=TempRecord.Seeds;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
% Snapshot of the current state of computing. It necessary for the parallel
|
% Snapshot of the current state of computing. It necessary for the parallel
|
||||||
% execution (i.e. to execute in a corretct way portion of code remotely or
|
% execution (i.e. to execute in a corretct way a portion of code remotely or
|
||||||
% on many core). The mandatory variables for local/remote parallel
|
% on many cores). The mandatory variables for local/remote parallel
|
||||||
% computing are stored in localVars struct.
|
% computing are stored in the localVars struct.
|
||||||
|
|
||||||
localVars = struct('TargetFun', TargetFun, ...
|
localVars = struct('TargetFun', TargetFun, ...
|
||||||
'ProposalFun', ProposalFun, ...
|
'ProposalFun', ProposalFun, ...
|
||||||
|
@ -110,9 +113,8 @@ localVars = struct('TargetFun', TargetFun, ...
|
||||||
'varargin',[]);
|
'varargin',[]);
|
||||||
|
|
||||||
|
|
||||||
% The user don't want to use parallel computing, or want to compute a
|
% User doesn't want to use parallel computing, or wants to compute a
|
||||||
% single chain. In this cases Random walk Metropolis-Hastings algorithm is
|
% single chain compute Random walk Metropolis-Hastings algorithm sequentially.
|
||||||
% computed sequentially.
|
|
||||||
|
|
||||||
if isnumeric(options_.parallel) || (nblck-fblck)==0,
|
if isnumeric(options_.parallel) || (nblck-fblck)==0,
|
||||||
fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0);
|
fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0);
|
||||||
|
@ -152,6 +154,7 @@ NewFile = fout(1).NewFile;
|
||||||
|
|
||||||
update_last_mh_history_file(MetropolisFolder, ModelName, record);
|
update_last_mh_history_file(MetropolisFolder, ModelName, record);
|
||||||
|
|
||||||
|
% Provide diagnostic output
|
||||||
skipline()
|
skipline()
|
||||||
disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.'])
|
disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.'])
|
||||||
disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.'])
|
disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.'])
|
||||||
|
|
|
@ -1,25 +1,25 @@
|
||||||
function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
|
function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
|
||||||
% PARALLEL CONTEXT
|
% function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
|
||||||
% This function contain the most computationally intensive portion of code in
|
% Contains the most computationally intensive portion of code in
|
||||||
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for'
|
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in that 'for'
|
||||||
% cycle and are completely independent than suitable to be executed in parallel way.
|
% cycle are completely independent to be suitable for parallel execution.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% o myimput [struc] The mandatory variables for local/remote
|
% o myimput [struc] The mandatory variables for local/remote
|
||||||
% parallel computing obtained from random_walk_metropolis_hastings.m
|
% parallel computing obtained from random_walk_metropolis_hastings.m
|
||||||
% function.
|
% function.
|
||||||
% o fblck and nblck [integer] The Metropolis-Hastings chains.
|
% o fblck and nblck [integer] The Metropolis-Hastings chains.
|
||||||
% o whoiam [integer] In concurrent programming a modality to refer to the differents thread running in parallel is needed.
|
% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed.
|
||||||
% The integer whoaim is the integer that
|
% The integer whoaim is the integer that
|
||||||
% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
|
% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
|
||||||
% cluster.
|
% cluster.
|
||||||
% o ThisMatlab [integer] Allows us to distinguish between the
|
% o ThisMatlab [integer] Allows us to distinguish between the
|
||||||
% 'main' matlab, the slave matlab worker, local matlab, remote matlab,
|
% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab,
|
||||||
% ... Then it is the index number of this slave machine in the cluster.
|
% ... Then it is the index number of this slave machine in the cluster.
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% o myoutput [struc]
|
% o myoutput [struc]
|
||||||
% If executed without parallel is the original output of 'for b =
|
% If executed without parallel, this is the original output of 'for b =
|
||||||
% fblck:nblck' otherwise a portion of it computed on a specific core or
|
% fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or
|
||||||
% remote machine. In this case:
|
% remote machine. In this case:
|
||||||
% record;
|
% record;
|
||||||
% irun;
|
% irun;
|
||||||
|
@ -31,23 +31,12 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,wh
|
||||||
%
|
%
|
||||||
% SPECIAL REQUIREMENTS.
|
% SPECIAL REQUIREMENTS.
|
||||||
% None.
|
% None.
|
||||||
|
%
|
||||||
% PARALLEL CONTEXT
|
% PARALLEL CONTEXT
|
||||||
% The most computationally intensive part of this function may be executed
|
% See the comments in the random_walk_metropolis_hastings.m funtion.
|
||||||
% in parallel. The code sutable to be executed in parallel on multi core or cluster machine,
|
|
||||||
% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
|
|
||||||
% Then the DYNARE parallel package contain a set of pairs matlab functios that can be executed in
|
|
||||||
% parallel and called name_function.m and name_function_core.m.
|
|
||||||
% In addition in the parallel package we have second set of functions used
|
|
||||||
% to manage the parallel computation.
|
|
||||||
%
|
|
||||||
% This function was the first function to be parallelized, later other
|
|
||||||
% functions have been parallelized using the same methodology.
|
|
||||||
% Then the comments write here can be used for all the other pairs of
|
|
||||||
% parallel functions and also for management funtions.
|
|
||||||
|
|
||||||
|
|
||||||
% Copyright (C) 2006-2013 Dynare Team
|
% Copyright (C) 2006-2015 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -74,11 +63,9 @@ end
|
||||||
TargetFun=myinputs.TargetFun;
|
TargetFun=myinputs.TargetFun;
|
||||||
ProposalFun=myinputs.ProposalFun;
|
ProposalFun=myinputs.ProposalFun;
|
||||||
xparam1=myinputs.xparam1;
|
xparam1=myinputs.xparam1;
|
||||||
vv=myinputs.vv;
|
|
||||||
mh_bounds=myinputs.mh_bounds;
|
mh_bounds=myinputs.mh_bounds;
|
||||||
ix2=myinputs.ix2;
|
last_draw=myinputs.ix2;
|
||||||
ilogpo2=myinputs.ilogpo2;
|
last_posterior=myinputs.ilogpo2;
|
||||||
ModelName=myinputs.ModelName;
|
|
||||||
fline=myinputs.fline;
|
fline=myinputs.fline;
|
||||||
npar=myinputs.npar;
|
npar=myinputs.npar;
|
||||||
nruns=myinputs.nruns;
|
nruns=myinputs.nruns;
|
||||||
|
@ -94,11 +81,9 @@ estim_params_ = myinputs.estim_params_;
|
||||||
options_ = myinputs.options_;
|
options_ = myinputs.options_;
|
||||||
M_ = myinputs.M_;
|
M_ = myinputs.M_;
|
||||||
oo_ = myinputs.oo_;
|
oo_ = myinputs.oo_;
|
||||||
varargin=myinputs.varargin;
|
|
||||||
|
|
||||||
% Necessary only for remote computing!
|
% Necessary only for remote computing!
|
||||||
if whoiam
|
if whoiam
|
||||||
Parallel=myinputs.Parallel;
|
|
||||||
% initialize persistent variables in priordens()
|
% initialize persistent variables in priordens()
|
||||||
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
|
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
|
||||||
end
|
end
|
||||||
|
@ -116,53 +101,51 @@ elseif strcmpi(ProposalFun,'rand_multivariate_student')
|
||||||
end
|
end
|
||||||
|
|
||||||
%
|
%
|
||||||
% NOW i run the (nblck-fblck+1) metropolis-hastings chains
|
% Now I run the (nblck-fblck+1) Metropolis-Hastings chains
|
||||||
%
|
%
|
||||||
|
|
||||||
proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
|
proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
|
||||||
|
|
||||||
jloop=0;
|
block_iter=0;
|
||||||
|
|
||||||
JSUM = 0;
|
for curr_block = fblck:nblck,
|
||||||
for b = fblck:nblck,
|
block_iter=block_iter+1;
|
||||||
jloop=jloop+1;
|
|
||||||
try
|
try
|
||||||
% This will not work if the master uses a random generator not
|
% This will not work if the master uses a random number generator not
|
||||||
% available in the slave (different Matlab version or
|
% available in the slave (different Matlab version or
|
||||||
% Matlab/Octave cluster). Therefor the trap.
|
% Matlab/Octave cluster). Therefore the trap.
|
||||||
%
|
%
|
||||||
% This set the random generator type (the seed is useless but
|
% Set the random number generator type (the seed is useless but needed by the function)
|
||||||
% needed by the function)
|
|
||||||
set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
|
set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
|
||||||
% This set the state
|
% Set the state of the RNG
|
||||||
set_dynare_random_generator_state(record.InitialSeeds(b).Unifor, record.InitialSeeds(b).Normal);
|
set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);
|
||||||
catch
|
catch
|
||||||
% If the state set by master is incompatible with the slave, we
|
% If the state set by master is incompatible with the slave, we only reseed
|
||||||
% only reseed
|
set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
|
||||||
set_dynare_seed(options_.DynareRandomStreams.seed+b);
|
|
||||||
end
|
end
|
||||||
if (options_.load_mh_file~=0) && (fline(b)>1) && OpenOldFile(b)
|
if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
|
||||||
load([BaseName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'])
|
load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
|
||||||
x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
|
x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
|
||||||
logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
|
logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
|
||||||
OpenOldFile(b) = 0;
|
OpenOldFile(curr_block) = 0;
|
||||||
else
|
else
|
||||||
x2 = zeros(InitSizeArray(b),npar);
|
x2 = zeros(InitSizeArray(curr_block),npar);
|
||||||
logpo2 = zeros(InitSizeArray(b),1);
|
logpo2 = zeros(InitSizeArray(curr_block),1);
|
||||||
end
|
end
|
||||||
|
%Prepare waiting bars
|
||||||
if whoiam
|
if whoiam
|
||||||
prc0=(b-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
|
prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
|
||||||
hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
|
hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
|
||||||
else
|
else
|
||||||
hh = dyn_waitbar(0,['Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
|
hh = dyn_waitbar(0,['Metropolis-Hastings (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
|
||||||
set(hh,'Name','Metropolis-Hastings');
|
set(hh,'Name','Metropolis-Hastings');
|
||||||
end
|
end
|
||||||
isux = 0;
|
accepted_draws_this_chain = 0;
|
||||||
jsux = 0;
|
accepted_draws_this_file = 0;
|
||||||
irun = fline(b);
|
draw_index_current_file = fline(curr_block); %get location of first draw in current block
|
||||||
j = 1;
|
draw_iter = 1;
|
||||||
while j <= nruns(b)
|
while draw_iter <= nruns(curr_block)
|
||||||
par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n);
|
par = feval(ProposalFun, last_draw(curr_block,:), proposal_covariance_Cholesky_decomposition, n);
|
||||||
if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub )
|
if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub )
|
||||||
try
|
try
|
||||||
logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
|
logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
|
||||||
|
@ -172,29 +155,29 @@ for b = fblck:nblck,
|
||||||
else
|
else
|
||||||
logpost = -inf;
|
logpost = -inf;
|
||||||
end
|
end
|
||||||
if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b))
|
if (logpost > -inf) && (log(rand) < logpost-last_posterior(curr_block))
|
||||||
x2(irun,:) = par;
|
x2(draw_index_current_file,:) = par;
|
||||||
ix2(b,:) = par;
|
last_draw(curr_block,:) = par;
|
||||||
logpo2(irun) = logpost;
|
logpo2(draw_index_current_file) = logpost;
|
||||||
ilogpo2(b) = logpost;
|
last_posterior(curr_block) = logpost;
|
||||||
isux = isux + 1;
|
accepted_draws_this_chain = accepted_draws_this_chain + 1;
|
||||||
jsux = jsux + 1;
|
accepted_draws_this_file = accepted_draws_this_file + 1;
|
||||||
else
|
else
|
||||||
x2(irun,:) = ix2(b,:);
|
x2(draw_index_current_file,:) = last_draw(curr_block,:);
|
||||||
logpo2(irun) = ilogpo2(b);
|
logpo2(draw_index_current_file) = last_posterior(curr_block);
|
||||||
end
|
end
|
||||||
prtfrc = j/nruns(b);
|
prtfrc = draw_iter/nruns(curr_block);
|
||||||
if (mod(j, 3)==0 && ~whoiam) || (mod(j,50)==0 && whoiam)
|
if (mod(draw_iter, 3)==0 && ~whoiam) || (mod(draw_iter,50)==0 && whoiam)
|
||||||
dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', isux/j)]);
|
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
|
end
|
||||||
if (irun == InitSizeArray(b)) || (j == nruns(b)) % Now I save the simulations
|
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(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
|
save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2');
|
||||||
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
|
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
|
||||||
fprintf(fidlog,['\n']);
|
fprintf(fidlog,['\n']);
|
||||||
fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']);
|
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']);
|
||||||
fprintf(fidlog,' \n');
|
fprintf(fidlog,' \n');
|
||||||
fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']);
|
fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']);
|
||||||
fprintf(fidlog,[' Acceptance ratio......: ' num2str(jsux/length(logpo2)) '\n']);
|
fprintf(fidlog,[' Acceptance ratio......: ' num2str(accepted_draws_this_file/length(logpo2)) '\n']);
|
||||||
fprintf(fidlog,[' Posterior mean........:\n']);
|
fprintf(fidlog,[' Posterior mean........:\n']);
|
||||||
for i=1:length(x2(1,:))
|
for i=1:length(x2(1,:))
|
||||||
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
|
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
|
||||||
|
@ -212,32 +195,32 @@ for b = fblck:nblck,
|
||||||
fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']);
|
fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']);
|
||||||
fprintf(fidlog,' \n');
|
fprintf(fidlog,' \n');
|
||||||
fclose(fidlog);
|
fclose(fidlog);
|
||||||
jsux = 0;
|
accepted_draws_this_file = 0;
|
||||||
if j == nruns(b) % I record the last draw...
|
if draw_iter == nruns(curr_block) % I record the last draw...
|
||||||
record.LastParameters(b,:) = x2(end,:);
|
record.LastParameters(curr_block,:) = x2(end,:);
|
||||||
record.LastLogPost(b) = logpo2(end);
|
record.LastLogPost(curr_block) = logpo2(end);
|
||||||
end
|
end
|
||||||
% size of next file in chain b
|
% size of next file in chain curr_block
|
||||||
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
|
InitSizeArray(curr_block) = min(nruns(curr_block)-draw_iter,MAX_nruns);
|
||||||
% initialization of next file if necessary
|
% initialization of next file if necessary
|
||||||
if InitSizeArray(b)
|
if InitSizeArray(curr_block)
|
||||||
x2 = zeros(InitSizeArray(b),npar);
|
x2 = zeros(InitSizeArray(curr_block),npar);
|
||||||
logpo2 = zeros(InitSizeArray(b),1);
|
logpo2 = zeros(InitSizeArray(curr_block),1);
|
||||||
NewFile(b) = NewFile(b) + 1;
|
NewFile(curr_block) = NewFile(curr_block) + 1;
|
||||||
irun = 0;
|
draw_index_current_file = 0;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
j=j+1;
|
draw_iter=draw_iter+1;
|
||||||
irun = irun + 1;
|
draw_index_current_file = draw_index_current_file + 1;
|
||||||
end% End of the simulations for one mh-block.
|
end% End of the simulations for one mh-block.
|
||||||
record.AcceptanceRatio(b) = isux/j;
|
record.AcceptanceRatio(curr_block) = accepted_draws_this_chain/draw_iter;
|
||||||
dyn_waitbar_close(hh);
|
dyn_waitbar_close(hh);
|
||||||
[record.LastSeeds(b).Unifor, record.LastSeeds(b).Normal] = get_dynare_random_generator_state();
|
[record.LastSeeds(curr_block).Unifor, record.LastSeeds(curr_block).Normal] = get_dynare_random_generator_state();
|
||||||
OutputFileName(jloop,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']};
|
OutputFileName(block_iter,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(curr_block) '.mat']};
|
||||||
end% End of the loop over the mh-blocks.
|
end% End of the loop over the mh-blocks.
|
||||||
|
|
||||||
|
|
||||||
myoutput.record = record;
|
myoutput.record = record;
|
||||||
myoutput.irun = irun;
|
myoutput.irun = draw_index_current_file;
|
||||||
myoutput.NewFile = NewFile;
|
myoutput.NewFile = NewFile;
|
||||||
myoutput.OutputFileName = OutputFileName;
|
myoutput.OutputFileName = OutputFileName;
|
|
@ -6,6 +6,7 @@ MODFILES = \
|
||||||
estimation/fs2000_MCMC_jumping_covariance.mod \
|
estimation/fs2000_MCMC_jumping_covariance.mod \
|
||||||
estimation/fs2000_initialize_from_calib.mod \
|
estimation/fs2000_initialize_from_calib.mod \
|
||||||
estimation/fs2000_calibrated_covariance.mod \
|
estimation/fs2000_calibrated_covariance.mod \
|
||||||
|
estimation/MH_recover/fs2000_recover.mod \
|
||||||
gsa/ls2003.mod \
|
gsa/ls2003.mod \
|
||||||
gsa/ls2003a.mod \
|
gsa/ls2003a.mod \
|
||||||
ramst.mod \
|
ramst.mod \
|
||||||
|
|
|
@ -0,0 +1,139 @@
|
||||||
|
/*
|
||||||
|
* 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_.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'])
|
||||||
|
|
||||||
|
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']);
|
||||||
|
|
||||||
|
if max(max(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']);
|
||||||
|
|
||||||
|
if max(max(temp1.x2-temp2.x2))>1e-10
|
||||||
|
error('Draws of affected chain''s unaffected files are not the same')
|
||||||
|
end
|
Loading…
Reference in New Issue