2013-09-16 17:19:56 +02:00
function oo_ = McMCDiagnostics ( options_, estim_params_, M_, oo_)
2008-06-09 12:40:48 +02:00
% function McMCDiagnostics
2017-05-16 15:10:20 +02:00
% Computes convergence tests
%
% INPUTS
2008-09-03 18:05:35 +02:00
% options_ [structure]
% estim_params_ [structure]
% M_ [structure]
%
2017-05-16 15:10:20 +02:00
% OUTPUTS
% oo_ [structure]
2008-01-21 18:29:24 +01:00
%
% SPECIAL REQUIREMENTS
% none
2011-02-04 17:27:33 +01:00
%
2010-05-31 11:46:38 +02:00
% PARALLEL CONTEXT
2016-05-19 17:13:39 +02:00
% See the comment in posterior_sampler.m funtion.
2010-05-31 11:46:38 +02:00
2017-10-10 10:05:59 +02:00
% Copyright (C) 2005-2018 Dynare Team
2008-08-01 14:40:33 +02:00
%
% 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/>.
2005-09-10 12:04:22 +02:00
2013-11-20 18:03:12 +01:00
OutputFolder = CheckPath ( ' Output' , M_ . dname ) ;
MetropolisFolder = CheckPath ( ' metropolis' , M_ . dname ) ;
ModelName = M_ . fname ;
2005-09-10 12:04:22 +02:00
TeX = options_ . TeX ;
nblck = options_ . mh_nblck ;
npar = estim_params_ . nvx ;
npar = npar + estim_params_ . nvn ;
npar = npar + estim_params_ . ncx ;
npar = npar + estim_params_ . ncn ;
npar = npar + estim_params_ . np ;
MAX_nruns = ceil ( options_ . MaxNumberOfBytes / ( npar + 2 ) / 8 ) ;
2013-11-20 18:03:12 +01:00
load_last_mh_history_file ( MetropolisFolder , ModelName ) ;
2013-11-22 11:28:24 +01:00
NumberOfMcFilesPerBlock = record . LastFileNumber ;
2005-09-10 12:04:22 +02:00
2013-11-22 11:28:24 +01:00
% Check that the declared number of blocks is consistent with informations saved in mh-history files.
if ~ isequal ( nblck , record . Nblck )
disp ( [ ' Estimation::mcmc::diagnostics: The number of MCMC chains you declared (' num2str ( nblck ) ' ) is inconsistent with the information available in the mh-history files (' num2str ( record . Nblck ) ' chains)!' ] )
disp ( [ ' I reset the number of MCMC chains to ' num2str ( record . Nblck ) ' .' ] )
nblck = record . Nblck ;
end
% check if all the mh files are available.
issue_an_error_message = 0 ;
for b = 1 : nblck
nfiles = length ( dir ( [ MetropolisFolder , filesep , ModelName ' _mh*_blck' num2str ( b ) ' .mat' ] ) ) ;
if ~ isequal ( NumberOfMcFilesPerBlock , nfiles )
issue_an_error_message = 1 ;
disp ( [ ' Estimation::mcmc::diagnostics: The number of MCMC files in chain ' num2str ( b ) ' is ' num2str ( nfiles ) ' while the mh-history files indicate that we should have ' num2str ( NumberOfMcFilesPerBlock ) ' MCMC files per chain!' ] )
end
end
if issue_an_error_message
error ( ' Estimation::mcmc::diagnostics: I cannot proceed because some MCMC files are missing. Check your MCMC files...' )
end
2013-05-18 19:00:11 +02:00
2015-09-25 12:39:05 +02:00
% compute inefficiency factor
FirstMhFile = record . KeepedDraws . FirstMhFile ;
FirstLine = record . KeepedDraws . FirstLine ;
TotalNumberOfMhFiles = sum ( record . MhDraws ( : , 2 ) ) ;
TotalNumberOfMhDraws = sum ( record . MhDraws ( : , 1 ) ) ;
FirstMhFile = record . KeepedDraws . FirstMhFile ;
NumberOfDraws = TotalNumberOfMhDraws - floor ( options_ . mh_drop * TotalNumberOfMhDraws ) ;
2017-10-10 10:05:59 +02:00
param_name = { } ;
param_name_tex = { } ;
for jj = 1 : npar
2016-06-07 16:12:10 +02:00
if options_ . TeX
2017-10-10 10:05:59 +02:00
[ par_name_temp , par_name_tex_temp ] = get_the_name ( jj , options_ . TeX , M_ , estim_params_ , options_ ) ;
param_name = vertcat ( param_name , par_name_temp ) ;
2016-06-07 16:12:10 +02:00
par_name_tex_temp = strrep ( par_name_tex_temp , ' $' , ' ' ) ;
2017-10-10 10:05:59 +02:00
param_name_tex = vertcat ( param_name_tex , par_name_tex_temp ) ;
2016-06-07 16:12:10 +02:00
else
2017-10-10 10:05:59 +02:00
par_name_temp = get_the_name ( jj , options_ . TeX , M_ , estim_params_ , options_ ) ;
param_name = vertcat ( param_name , par_name_temp ) ;
2016-06-07 16:12:10 +02:00
end
2017-10-10 10:05:59 +02:00
Draws = GetAllPosteriorDraws ( jj , FirstMhFile , FirstLine , TotalNumberOfMhFiles , NumberOfDraws ) ;
Draws = reshape ( Draws , [ NumberOfDraws nblck ] ) ;
2015-09-25 12:39:05 +02:00
Nc = min ( 1000 , NumberOfDraws / 2 ) ;
2017-10-10 10:05:59 +02:00
for ll = 1 : nblck
2015-09-25 12:39:05 +02:00
Ifac ( ll , jj ) = mcmc_ifac ( Draws ( : , ll ) , Nc ) ;
end
tmp = num2cell ( Ifac ( : , jj ) ) ;
end
2016-05-13 21:38:25 +02:00
my_title = ' MCMC Inefficiency factors per block' ;
2017-10-10 10:05:59 +02:00
IFAC_header = { ' Parameter' } ;
IFAC_header_tex = { ' Parameter' } ;
for j = 1 : nblck
IFAC_header = vertcat ( IFAC_header , [ ' Block ' int2str ( j ) ] ) ;
IFAC_header_tex = vertcat ( IFAC_header_tex , [ ' Block~' int2str ( j ) ] ) ;
2016-05-13 21:38:25 +02:00
end
2017-10-10 10:05:59 +02:00
lh = cellofchararraymaxlength ( param_name ) + 2 ;
dyntable ( options_ , my_title , IFAC_header , param_name , Ifac ' , lh , 12 , 3 ) ;
2015-09-25 12:39:05 +02:00
skipline ( )
2016-05-13 21:38:25 +02:00
if options_ . TeX
2017-10-10 10:05:59 +02:00
dyn_latex_table ( M_ , options_ , my_title , ' MCMC_inefficiency_factors' , IFAC_header_tex , param_name_tex , Ifac ' , lh , 12 , 3 ) ;
2016-05-13 21:38:25 +02:00
end
2015-09-25 12:39:05 +02:00
record . InefficiencyFactorsPerBlock = Ifac ;
update_last_mh_history_file ( MetropolisFolder , ModelName , record ) ;
2014-02-25 17:24:34 +01:00
PastDraws = sum ( record . MhDraws , 1 ) ;
LastFileNumber = PastDraws ( 2 ) ;
LastLineNumber = record . MhDraws ( end , 3 ) ;
NumberOfDraws = PastDraws ( 1 ) ;
2014-03-01 17:13:49 +01:00
if NumberOfDraws < = 2000
2016-11-22 10:47:09 +01:00
warning ( [ ' estimation:: MCMC convergence diagnostics are not computed because the total number of iterations is not bigger than 2000!' ] )
2014-02-25 17:24:34 +01:00
return
end
2013-05-18 19:00:11 +02:00
2013-11-22 11:28:24 +01:00
if nblck == 1 % Brooks and Gelman tests need more than one block
2013-09-16 17:19:56 +02:00
convergence_diagnostics_geweke = zeros ( npar , 4 + 2 * length ( options_ . convergence . geweke . taper_steps ) ) ;
2014-03-24 13:50:08 +01:00
if any ( options_ . convergence . geweke . geweke_interval < 0 ) || any ( options_ . convergence . geweke . geweke_interval > 1 ) || length ( options_ . convergence . geweke . geweke_interval ) ~= 2 ...
2017-05-16 15:10:20 +02:00
|| ( options_ . convergence . geweke . geweke_interval ( 2 ) - options_ . convergence . geweke . geweke_interval ( 1 ) < 0 )
2013-09-16 17:19:56 +02:00
fprintf ( ' \nCONVERGENCE DIAGNOSTICS: Invalid option for geweke_interval. Using the default of [0.2 0.5].\n' )
options_ . convergence . geweke . geweke_interval = [ 0.2 0.5 ] ;
end
2015-05-02 14:23:23 +02:00
first_obs_begin_sample = max ( 1 , ceil ( options_ . mh_drop * NumberOfDraws ) ) ;
last_obs_begin_sample = first_obs_begin_sample + round ( options_ . convergence . geweke . geweke_interval ( 1 ) * NumberOfDraws * ( 1 - options_ . mh_drop ) ) ;
first_obs_end_sample = first_obs_begin_sample + round ( options_ . convergence . geweke . geweke_interval ( 2 ) * NumberOfDraws * ( 1 - options_ . mh_drop ) ) ;
2017-10-10 10:05:59 +02:00
param_name = { } ;
2015-12-13 16:41:41 +01:00
if options_ . TeX
2017-10-10 10:05:59 +02:00
param_name_tex = { } ;
2015-12-13 16:41:41 +01:00
end
2017-05-16 15:10:20 +02:00
for jj = 1 : npar
2015-12-13 16:41:41 +01:00
if options_ . TeX
2017-10-10 10:05:59 +02:00
[ param_name_temp , param_name_tex_temp ] = get_the_name ( jj , options_ . TeX , M_ , estim_params_ , options_ ) ;
param_name_tex = vertcat ( param_name_tex , strrep ( param_name_tex_temp , ' $' , ' ' ) ) ;
param_name = vertcat ( param_name , param_name_temp ) ;
2015-12-13 16:41:41 +01:00
else
2017-10-10 10:05:59 +02:00
param_name_temp = get_the_name ( jj , options_ . TeX , M_ , estim_params_ , options_ ) ;
param_name = vertcat ( param_name , param_name_temp ) ;
2015-12-13 16:41:41 +01:00
end
2013-09-16 17:19:56 +02:00
end
2015-05-02 14:23:23 +02:00
fprintf ( ' \nGeweke (1992) Convergence Tests, based on means of draws %d to %d vs %d to %d.\n' , first_obs_begin_sample , last_obs_begin_sample , first_obs_end_sample , NumberOfDraws ) ;
2017-05-16 15:10:20 +02:00
fprintf ( ' p-values are for Chi2-test for equality of means.\n' ) ;
2017-10-10 10:05:59 +02:00
Geweke_header = { ' Parameter' ; ' Post. Mean' ; ' Post. Std' ; ' p-val No Taper' } ;
for ii = 1 : length ( options_ . convergence . geweke . taper_steps )
Geweke_header = vertcat ( Geweke_header , [ ' p-val ' num2str ( options_ . convergence . geweke . taper_steps ( ii ) ) , ' % Taper' ] ) ;
2013-09-16 17:19:56 +02:00
end
2015-12-13 16:41:41 +01:00
datamat = NaN ( npar , 3 + length ( options_ . convergence . geweke . taper_steps ) ) ;
2013-09-16 17:19:56 +02:00
for jj = 1 : npar
startline = 0 ;
for n = 1 : NumberOfMcFilesPerBlock
2013-11-20 18:03:12 +01:00
load ( [ MetropolisFolder ' /' ModelName ' _mh' , int2str ( n ) , ' _blck1.mat' ] , ' x2' ) ;
2013-09-16 17:19:56 +02:00
nx2 = size ( x2 , 1 ) ;
param_draws ( startline + ( 1 : nx2 ) , 1 ) = x2 ( : , jj ) ;
startline = startline + nx2 ;
end
[ results_vec , results_struct ] = geweke_moments ( param_draws , options_ ) ;
convergence_diagnostics_geweke ( jj , : ) = results_vec ;
2017-05-16 15:10:20 +02:00
2013-09-16 17:19:56 +02:00
param_draws1 = param_draws ( first_obs_begin_sample : last_obs_begin_sample , : ) ;
param_draws2 = param_draws ( first_obs_end_sample : end , : ) ;
[ results_vec1 ] = geweke_moments ( param_draws1 , options_ ) ;
[ results_vec2 ] = geweke_moments ( param_draws2 , options_ ) ;
2017-05-16 15:10:20 +02:00
2013-09-16 17:19:56 +02:00
results_struct = geweke_chi2_test ( results_vec1 , results_vec2 , results_struct , options_ ) ;
2017-10-10 10:05:59 +02:00
oo_ . convergence . geweke . ( param_name { jj } ) = results_struct ;
2015-12-13 16:41:41 +01:00
datamat ( jj , : ) = [ results_struct . posteriormean , results_struct . posteriorstd , results_struct . prob_chi2_test ] ;
end
lh = size ( param_name , 2 ) + 2 ;
2017-10-10 10:05:59 +02:00
dyntable ( options_ , ' ' , Geweke_header , param_name , datamat , lh , 12 , 3 ) ;
2015-12-13 16:41:41 +01:00
if options_ . TeX
2017-10-10 10:05:59 +02:00
Geweke_tex_header = { ' Parameter' ; ' Mean' ; ' Std' ; ' No\ Taper' } ;
additional_header = { [ ' & \multicolumn{2}{c}{Posterior} & \multicolumn{' , num2str ( 1 + length ( options_ . convergence . geweke . taper_steps ) ) , ' }{c}{p-values} \\' ] ,
2019-12-20 16:28:06 +01:00
[ ' \cmidrule(r{.75em}){2-3} \cmidrule(r{.75em}){4-' , num2str ( 4 + length ( options_ . convergence . geweke . taper_steps ) ) , ' }' ] } ;
2015-12-13 16:41:41 +01:00
for ii = 1 : length ( options_ . convergence . geweke . taper_steps )
2017-10-10 10:05:59 +02:00
Geweke_tex_header = vertcat ( Geweke_tex_header , [ num2str ( options_ . convergence . geweke . taper_steps ( ii ) ) , ' \%%\ Taper' ] ) ;
2015-12-13 16:41:41 +01:00
end
2017-10-10 10:05:59 +02:00
headers = Geweke_tex_header ;
lh = cellofchararraymaxlength ( param_name_tex ) + 2 ;
2015-12-13 16:41:41 +01:00
my_title = sprintf ( ' Geweke (1992) Convergence Tests, based on means of draws %d to %d vs %d to %d. p-values are for $\\\\chi^2$-test for equality of means.' , first_obs_begin_sample , last_obs_begin_sample , first_obs_end_sample , NumberOfDraws ) ;
2017-10-10 10:05:59 +02:00
dyn_latex_table ( M_ , options_ , my_title , ' geweke' , headers , param_name_tex , datamat , lh , 12 , 4 , additional_header ) ;
2013-09-16 17:19:56 +02:00
end
skipline ( 2 ) ;
2017-05-16 15:10:20 +02:00
2016-11-22 10:47:09 +01:00
if options_ . convergence . rafterylewis . indicator
if any ( options_ . convergence . rafterylewis . qrs < 0 ) || any ( options_ . convergence . rafterylewis . qrs > 1 ) || length ( options_ . convergence . rafterylewis . qrs ) ~= 3 ...
2017-05-16 15:10:20 +02:00
|| ( options_ . convergence . rafterylewis . qrs ( 1 ) - options_ . convergence . rafterylewis . qrs ( 2 ) < = 0 )
2016-11-22 10:47:09 +01:00
fprintf ( ' \nCONVERGENCE DIAGNOSTICS: Invalid option for raftery_lewis_qrs. Using the default of [0.025 0.005 0.95].\n' )
options_ . convergence . rafterylewis . qrs = [ 0.025 0.005 0.95 ] ;
2017-05-16 15:10:20 +02:00
end
2016-11-22 10:47:09 +01:00
Raftery_Lewis_q = options_ . convergence . rafterylewis . qrs ( 1 ) ;
Raftery_Lewis_r = options_ . convergence . rafterylewis . qrs ( 2 ) ;
Raftery_Lewis_s = options_ . convergence . rafterylewis . qrs ( 3 ) ;
oo_ . Raftery_Lewis = raftery_lewis ( x2 , Raftery_Lewis_q , Raftery_Lewis_r , Raftery_Lewis_s ) ;
oo_ . Raftery_Lewis . parameter_names = param_name ;
2017-05-16 15:10:20 +02:00
my_title = sprintf ( ' Raftery/Lewis (1992) Convergence Diagnostics, based on quantile q=%4.3f with precision r=%4.3f with probability s=%4.3f.' , Raftery_Lewis_q , Raftery_Lewis_r , Raftery_Lewis_s ) ;
2017-10-10 10:05:59 +02:00
headers = { ' Variables' ; ' M (burn-in)' ; ' N (req. draws)' ; ' N+M (total draws)' ; ' k (thinning)' } ;
2016-11-22 10:47:09 +01:00
raftery_data_mat = [ oo_ . Raftery_Lewis . M_burn , oo_ . Raftery_Lewis . N_prec , oo_ . Raftery_Lewis . N_total , oo_ . Raftery_Lewis . k_thin ] ;
raftery_data_mat = [ raftery_data_mat ; max ( raftery_data_mat ) ] ;
2017-10-10 10:05:59 +02:00
labels_Raftery_Lewis = vertcat ( param_name , ' Maximum' ) ;
lh = cellofchararraymaxlength ( labels_Raftery_Lewis ) + 2 ;
dyntable ( options_ , my_title , headers , labels_Raftery_Lewis , raftery_data_mat , lh , 10 , 0 ) ;
2016-11-22 10:47:09 +01:00
if options_ . TeX
2017-10-10 10:05:59 +02:00
labels_Raftery_Lewis_tex = vertcat ( param_name_tex , ' Maximum' ) ;
lh = cellofchararraymaxlength ( labels_Raftery_Lewis_tex ) + 2 ;
dyn_latex_table ( M_ , options_ , my_title , ' raftery_lewis' , headers , labels_Raftery_Lewis_tex , raftery_data_mat , lh , 10 , 0 ) ;
2017-05-16 15:10:20 +02:00
end
2016-11-22 10:47:09 +01:00
end
2017-05-16 15:10:20 +02:00
2017-05-16 12:42:01 +02:00
return
2013-09-16 17:19:56 +02:00
end
2005-09-10 12:04:22 +02:00
Origin = 1000 ;
2017-05-16 15:10:20 +02:00
StepSize = ceil ( ( NumberOfDraws - Origin ) / 100 ) ; % So that the computational time does not
ALPHA = 0.2 ; % increase too much with the number of simulations.
2005-09-10 12:04:22 +02:00
time = 1 : NumberOfDraws ;
xx = Origin : StepSize : NumberOfDraws ;
NumberOfLines = length ( xx ) ;
tmp = zeros ( NumberOfDraws * nblck , 3 ) ;
UDIAG = zeros ( NumberOfLines , 6 , npar ) ;
2006-03-06 11:00:32 +01:00
2005-09-10 12:04:22 +02:00
if NumberOfDraws < Origin
2013-11-22 11:28:24 +01:00
disp ( ' Estimation::mcmc::diagnostics: The number of simulations is too small to compute the MCMC convergence diagnostics.' )
2005-09-10 12:04:22 +02:00
return
end
2006-03-06 11:00:32 +01:00
2015-05-12 08:53:57 +02:00
if TeX && any ( strcmp ( ' eps' , cellstr ( options_ . graph_format ) ) )
2016-06-14 09:48:35 +02:00
fidTeX = fopen ( [ OutputFolder ' /' ModelName ' _UnivariateDiagnostics.tex' ] , ' w' ) ;
2009-09-24 14:48:28 +02:00
fprintf ( fidTeX , ' %% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n' ) ;
fprintf ( fidTeX , [ ' %% ' datestr ( now , 0 ) ' \n' ] ) ;
fprintf ( fidTeX , ' \n' ) ;
2005-09-10 12:04:22 +02:00
end
2006-03-06 11:00:32 +01:00
2013-11-22 11:28:24 +01:00
disp ( ' Estimation::mcmc::diagnostics: Univariate convergence diagnostic, Brooks and Gelman (1998):' )
2009-05-27 18:30:30 +02:00
2010-05-31 11:46:38 +02:00
% The mandatory variables for local/remote parallel
% computing are stored in localVars struct.
2013-11-20 18:03:12 +01:00
localVars . MetropolisFolder = MetropolisFolder ;
2009-05-27 18:30:30 +02:00
localVars . nblck = nblck ;
localVars . NumberOfMcFilesPerBlock = NumberOfMcFilesPerBlock ;
localVars . Origin = Origin ;
localVars . StepSize = StepSize ;
2010-02-18 16:59:26 +01:00
localVars . mh_drop = options_ . mh_drop ;
2009-05-27 18:30:30 +02:00
localVars . NumberOfDraws = NumberOfDraws ;
localVars . NumberOfLines = NumberOfLines ;
localVars . time = time ;
localVars . M_ = M_ ;
2010-05-31 11:46:38 +02:00
% Like sequential execution!
2017-05-16 12:42:01 +02:00
if isnumeric ( options_ . parallel )
2009-05-27 18:30:30 +02:00
fout = McMCDiagnostics_core ( localVars , 1 , npar , 0 ) ;
UDIAG = fout . UDIAG ;
clear fout
2011-02-04 17:17:48 +01:00
% Parallel execution!
2009-05-27 18:30:30 +02:00
else
2013-11-20 18:03:12 +01:00
ModelName = ModelName ;
2010-04-14 17:20:29 +02:00
if ~ isempty ( M_ . bvar )
2013-11-20 18:03:12 +01:00
ModelName = [ ModelName ' _bvar' ] ;
2010-04-14 17:20:29 +02:00
end
NamFileInput = { [ M_ . dname ' /metropolis/' ] , [ ModelName ' _mh*_blck*.mat' ] } ;
2017-05-16 15:10:20 +02:00
2010-04-14 17:20:29 +02:00
[ fout , nBlockPerCPU , totCPU ] = masterParallel ( options_ . parallel , 1 , npar , NamFileInput , ' McMCDiagnostics_core' , localVars , [ ] , options_ . parallel_info ) ;
2009-06-09 12:03:29 +02:00
UDIAG = fout ( 1 ) . UDIAG ;
2017-05-16 12:42:01 +02:00
for j = 2 : totCPU
2009-05-27 18:30:30 +02:00
UDIAG = cat ( 3 , UDIAG , fout ( j ) . UDIAG ) ;
2005-09-10 12:04:22 +02:00
end
end
2009-05-27 18:30:30 +02:00
2005-09-10 12:04:22 +02:00
UDIAG ( : , [ 2 4 6 ] , : ) = UDIAG ( : , [ 2 4 6 ] , : ) / nblck ;
2013-07-10 17:12:34 +02:00
skipline ( )
2017-05-16 15:10:20 +02:00
clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp ;
2005-09-10 12:04:22 +02:00
pages = floor ( npar / 3 ) ;
2017-05-16 15:10:20 +02:00
k = 0 ;
2005-09-10 12:04:22 +02:00
for i = 1 : pages
2017-03-22 16:54:11 +01:00
h = dyn_figure ( options_ . nodisplay , ' Name' , ' MCMC univariate convergence diagnostic (Brooks and Gelman,1998)' ) ;
2009-09-24 14:48:28 +02:00
boxplot = 1 ;
for j = 1 : 3 % Loop over parameters
k = k + 1 ;
2011-09-14 19:03:56 +02:00
[ nam , namtex ] = get_the_name ( k , TeX , M_ , estim_params_ , options_ ) ;
2009-09-24 14:48:28 +02:00
for crit = 1 : 3 % Loop over criteria
if crit == 1
plt1 = UDIAG ( : , 1 , k ) ;
plt2 = UDIAG ( : , 2 , k ) ;
2017-05-16 15:10:20 +02:00
namnam = [ nam , ' (Interval)' ] ;
2009-09-24 14:48:28 +02:00
elseif crit == 2
plt1 = UDIAG ( : , 3 , k ) ;
plt2 = UDIAG ( : , 4 , k ) ;
namnam = [ nam , ' (m2)' ] ;
2017-05-16 15:10:20 +02:00
elseif crit == 3
2009-09-24 14:48:28 +02:00
plt1 = UDIAG ( : , 5 , k ) ;
plt2 = UDIAG ( : , 6 , k ) ;
namnam = [ nam , ' (m3)' ] ;
end
subplot ( 3 , 3 , boxplot ) ;
plot ( xx , plt1 , ' -b' ) ; % Pooled
hold on ;
plot ( xx , plt2 , ' -r' ) ; % Within (mean)
hold off ;
xlim ( [ xx ( 1 ) xx ( NumberOfLines ) ] )
2020-06-17 20:23:49 +02:00
if TeX
title ( namtex , ' interpreter' , ' latex' )
else
title ( namnam , ' Interpreter' , ' none' )
end
2009-09-24 14:48:28 +02:00
boxplot = boxplot + 1 ;
end
end
2017-03-23 17:59:05 +01:00
dyn_saveas ( h , [ OutputFolder ' /' ModelName ' _udiag' int2str ( i ) ] , options_ . nodisplay , options_ . graph_format ) ;
2015-05-12 08:53:57 +02:00
if TeX && any ( strcmp ( ' eps' , cellstr ( options_ . graph_format ) ) )
2009-09-24 14:48:28 +02:00
fprintf ( fidTeX , ' \\begin{figure}[H]\n' ) ;
fprintf ( fidTeX , ' \\centering \n' ) ;
2016-06-19 12:38:31 +02:00
fprintf ( fidTeX , ' \\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n' , options_ . figures . textwidth * min ( ( boxplot - 1 ) / 3 , 1 ) , [ OutputFolder ' /' ModelName ] , int2str ( i ) ) ;
2009-09-24 14:48:28 +02:00
fprintf ( fidTeX , ' \\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n' ) ;
fprintf ( fidTeX , ' The first, second and third columns are respectively the criteria based on\n' ) ;
fprintf ( fidTeX , ' the eighty percent interval, the second and third moments.}' ) ;
fprintf ( fidTeX , ' \\label{Fig:UnivariateDiagnostics:%s}\n' , int2str ( i ) ) ;
fprintf ( fidTeX , ' \\end{figure}\n' ) ;
fprintf ( fidTeX , ' \n' ) ;
end
2005-09-10 12:04:22 +02:00
end
reste = npar - k ;
if reste
2009-09-24 14:48:28 +02:00
if reste == 1
nr = 3 ;
nc = 1 ;
2017-05-16 12:42:01 +02:00
elseif reste == 2
2009-09-24 14:48:28 +02:00
nr = 2 ;
nc = 3 ;
2005-09-10 12:04:22 +02:00
end
2017-03-22 16:54:11 +01:00
h = dyn_figure ( options_ . nodisplay , ' Name' , ' MCMC univariate convergence diagnostic (Brooks and Gelman, 1998)' ) ;
2009-09-24 14:48:28 +02:00
boxplot = 1 ;
for j = 1 : reste
k = k + 1 ;
2011-09-14 19:03:56 +02:00
[ nam , namtex ] = get_the_name ( k , TeX , M_ , estim_params_ , options_ ) ;
2009-09-24 14:48:28 +02:00
for crit = 1 : 3
if crit == 1
plt1 = UDIAG ( : , 1 , k ) ;
plt2 = UDIAG ( : , 2 , k ) ;
2017-05-16 15:10:20 +02:00
namnam = [ nam , ' (Interval)' ] ;
2020-06-17 20:23:49 +02:00
if TeX
namnamtex = [ namtex , ' (Interval)' ] ;
end
2009-09-24 14:48:28 +02:00
elseif crit == 2
plt1 = UDIAG ( : , 3 , k ) ;
plt2 = UDIAG ( : , 4 , k ) ;
namnam = [ nam , ' (m2)' ] ;
2020-06-17 20:23:49 +02:00
if TeX
namnamtex = [ namtex , ' (m2)' ] ;
end
2017-05-16 15:10:20 +02:00
elseif crit == 3
2009-09-24 14:48:28 +02:00
plt1 = UDIAG ( : , 5 , k ) ;
plt2 = UDIAG ( : , 6 , k ) ;
namnam = [ nam , ' (m3)' ] ;
2020-06-17 20:23:49 +02:00
if TeX
namnamtex = [ namtex , ' (m3)' ] ;
2010-09-26 13:42:41 +02:00
end
2009-09-24 14:48:28 +02:00
end
subplot ( nr , nc , boxplot ) ;
2010-01-05 11:46:10 +01:00
plot ( xx , plt1 , ' -b' ) ; % Pooled
2009-09-24 14:48:28 +02:00
hold on ;
2010-01-05 11:46:10 +01:00
plot ( xx , plt2 , ' -r' ) ; % Within (mean)
2009-09-24 14:48:28 +02:00
hold off ;
xlim ( [ xx ( 1 ) xx ( NumberOfLines ) ] ) ;
2020-06-17 20:23:49 +02:00
if TeX
title ( namnamtex , ' Interpreter' , ' latex' ) ;
else
title ( namnam , ' Interpreter' , ' none' ) ;
end
2009-09-24 14:48:28 +02:00
boxplot = boxplot + 1 ;
end
end
2017-03-23 17:59:05 +01:00
dyn_saveas ( h , [ OutputFolder ' /' ModelName ' _udiag' int2str ( pages + 1 ) ] , options_ . nodisplay , options_ . graph_format ) ;
2015-05-12 08:53:57 +02:00
if TeX && any ( strcmp ( ' eps' , cellstr ( options_ . graph_format ) ) )
2009-09-24 14:48:28 +02:00
fprintf ( fidTeX , ' \\begin{figure}[H]\n' ) ;
fprintf ( fidTeX , ' \\centering \n' ) ;
2016-06-19 12:38:31 +02:00
fprintf ( fidTeX , ' \\includegraphics[width=%2.2f\\textwidth]{%s_udiag%s}\n' , options_ . figures . textwidth * min ( ( boxplot - 1 ) / nc , 1 ) , [ OutputFolder ' /' ModelName ] , int2str ( pages + 1 ) ) ;
2009-09-24 14:48:28 +02:00
if reste == 2
fprintf ( fidTeX , ' \\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n' ) ;
fprintf ( fidTeX , ' The first, second and third columns are respectively the criteria based on\n' ) ;
fprintf ( fidTeX , ' the eighty percent interval, the second and third moments.}' ) ;
elseif reste == 1
fprintf ( fidTeX , ' \\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n' ) ;
fprintf ( fidTeX , ' The first, second and third rows are respectively the criteria based on\n' ) ;
fprintf ( fidTeX , ' the eighty percent interval, the second and third moments.}' ) ;
end
fprintf ( fidTeX , ' \\label{Fig:UnivariateDiagnostics:%s}\n' , int2str ( pages + 1 ) ) ;
fprintf ( fidTeX , ' \\end{figure}\n' ) ;
fprintf ( fidTeX , ' \n' ) ;
fprintf ( fidTeX , ' % End Of TeX file.' ) ;
fclose ( fidTeX ) ;
2005-09-10 12:04:22 +02:00
end
end % if reste > 0
clear UDIAG ;
2013-11-20 18:03:12 +01:00
%
% Multivariate diagnostic.
%
2015-05-12 08:53:57 +02:00
if TeX && any ( strcmp ( ' eps' , cellstr ( options_ . graph_format ) ) )
2016-06-14 09:48:35 +02:00
fidTeX = fopen ( [ OutputFolder ' /' ModelName ' _MultivariateDiagnostics.tex' ] , ' w' ) ;
2009-09-24 14:48:28 +02:00
fprintf ( fidTeX , ' %% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n' ) ;
fprintf ( fidTeX , [ ' %% ' datestr ( now , 0 ) ' \n' ] ) ;
fprintf ( fidTeX , ' \n' ) ;
2005-09-10 12:04:22 +02:00
end
tmp = zeros ( NumberOfDraws * nblck , 3 ) ;
MDIAG = zeros ( NumberOfLines , 6 ) ;
for b = 1 : nblck
2009-09-24 14:48:28 +02:00
startline = 0 ;
for n = 1 : NumberOfMcFilesPerBlock
2013-11-20 18:03:12 +01:00
load ( [ MetropolisFolder ' /' ModelName ' _mh' , int2str ( n ) , ' _blck' int2str ( b ) ' .mat' ] , ' logpo2' ) ;
2009-09-24 14:48:28 +02:00
nlogpo2 = size ( logpo2 , 1 ) ;
tmp ( ( b - 1 ) * NumberOfDraws + startline + ( 1 : nlogpo2 ) , 1 ) = logpo2 ;
startline = startline + nlogpo2 ;
end
2005-09-10 12:04:22 +02:00
end
clear logpo2 ;
tmp ( : , 2 ) = kron ( transpose ( 1 : nblck ) , ones ( NumberOfDraws , 1 ) ) ;
2017-05-16 15:10:20 +02:00
tmp ( : , 3 ) = kron ( ones ( nblck , 1 ) , time ' ) ;
2005-09-10 12:04:22 +02:00
tmp = sortrows ( tmp , 1 ) ;
ligne = 0 ;
for iter = Origin : StepSize : NumberOfDraws
2009-09-24 14:48:28 +02:00
ligne = ligne + 1 ;
2010-04-14 17:20:29 +02:00
linea = ceil ( options_ . mh_drop * iter ) ;
2009-09-24 14:48:28 +02:00
n = iter - linea + 1 ;
cinf = round ( n * ALPHA / 2 ) ;
csup = round ( n * ( 1 - ALPHA / 2 ) ) ;
CINF = round ( nblck * n * ALPHA / 2 ) ;
CSUP = round ( nblck * n * ( 1 - ALPHA / 2 ) ) ;
temp = tmp ( find ( ( tmp ( : , 3 ) > = linea ) & ( tmp ( : , 3 ) < = iter ) ) , 1 : 2 ) ;
MDIAG ( ligne , 1 ) = temp ( CSUP , 1 ) - temp ( CINF , 1 ) ;
moyenne = mean ( temp ( : , 1 ) ) ; %% Pooled mean.
MDIAG ( ligne , 3 ) = sum ( ( temp ( : , 1 ) - moyenne ) .^ 2 ) / ( nblck * n - 1 ) ;
MDIAG ( ligne , 5 ) = sum ( abs ( temp ( : , 1 ) - moyenne ) .^ 3 ) / ( nblck * n - 1 ) ;
for i = 1 : nblck
pmet = temp ( find ( temp ( : , 2 ) == i ) ) ;
MDIAG ( ligne , 2 ) = MDIAG ( ligne , 2 ) + pmet ( csup , 1 ) - pmet ( cinf , 1 ) ;
2017-05-16 15:10:20 +02:00
moyenne = mean ( pmet , 1 ) ; %% Within mean.
2009-09-24 14:48:28 +02:00
MDIAG ( ligne , 4 ) = MDIAG ( ligne , 4 ) + sum ( ( pmet ( : , 1 ) - moyenne ) .^ 2 ) / ( n - 1 ) ;
MDIAG ( ligne , 6 ) = MDIAG ( ligne , 6 ) + sum ( abs ( pmet ( : , 1 ) - moyenne ) .^ 3 ) / ( n - 1 ) ;
end
2005-09-10 12:04:22 +02:00
end
2017-05-16 15:10:20 +02:00
MDIAG ( : , [ 2 4 6 ] , : ) = MDIAG ( : , [ 2 4 6 ] , : ) / nblck ;
2012-04-20 20:17:20 +02:00
2017-03-22 16:54:11 +01:00
h = dyn_figure ( options_ . nodisplay , ' Name' , ' Multivariate convergence diagnostic' ) ;
2005-09-10 12:04:22 +02:00
boxplot = 1 ;
for crit = 1 : 3
2009-09-24 14:48:28 +02:00
if crit == 1
plt1 = MDIAG ( : , 1 ) ;
plt2 = MDIAG ( : , 2 ) ;
2017-05-16 15:10:20 +02:00
namnam = ' Interval' ;
2009-09-24 14:48:28 +02:00
elseif crit == 2
plt1 = MDIAG ( : , 3 ) ;
plt2 = MDIAG ( : , 4 ) ;
namnam = ' m2' ;
2017-05-16 15:10:20 +02:00
elseif crit == 3
2009-09-24 14:48:28 +02:00
plt1 = MDIAG ( : , 5 ) ;
plt2 = MDIAG ( : , 6 ) ;
namnam = ' m3' ;
end
subplot ( 3 , 1 , boxplot ) ;
plot ( xx , plt1 , ' -b' ) ; % Pooled
hold on
plot ( xx , plt2 , ' -r' ) ; % Within (mean)
hold off
xlim ( [ xx ( 1 ) xx ( NumberOfLines ) ] )
title ( namnam , ' Interpreter' , ' none' ) ;
boxplot = boxplot + 1 ;
2005-09-10 12:04:22 +02:00
end
2017-03-23 17:59:05 +01:00
dyn_saveas ( h , [ OutputFolder ' /' ModelName ' _mdiag' ] , options_ . nodisplay , options_ . graph_format ) ;
2012-04-20 20:17:20 +02:00
2015-05-12 08:53:57 +02:00
if TeX && any ( strcmp ( ' eps' , cellstr ( options_ . graph_format ) ) )
2009-09-24 14:48:28 +02:00
fprintf ( fidTeX , ' \\begin{figure}[H]\n' ) ;
fprintf ( fidTeX , ' \\centering \n' ) ;
2016-06-18 14:25:00 +02:00
fprintf ( fidTeX , ' \\includegraphics[width=0.8\\textwidth]{%s_mdiag}\n' , [ OutputFolder ' /' ModelName ] ) ;
2009-09-24 14:48:28 +02:00
fprintf ( fidTeX , ' \\caption{Multivariate convergence diagnostics for the Metropolis-Hastings.\n' ) ;
fprintf ( fidTeX , ' The first, second and third rows are respectively the criteria based on\n' ) ;
fprintf ( fidTeX , ' the eighty percent interval, the second and third moments. The different \n' ) ;
fprintf ( fidTeX , ' parameters are aggregated using the posterior kernel.}' ) ;
fprintf ( fidTeX , ' \\label{Fig:MultivariateDiagnostics}\n' ) ;
fprintf ( fidTeX , ' \\end{figure}\n' ) ;
fprintf ( fidTeX , ' \n' ) ;
2016-11-21 19:25:54 +01:00
fprintf ( fidTeX , ' %% End Of TeX file.' ) ;
2009-09-24 14:48:28 +02:00
fclose ( fidTeX ) ;
2015-09-25 12:39:05 +02:00
end