Description of McMCDiagnostics

0001 function McMCDiagnostics(options_, estim_params_, M_)
0002 % function McMCDiagnostics
0003 % Computes convergence tests
0004 %
0005 % INPUTS
0006 %   options_         [structure]
0007 %   estim_params_    [structure]
0008 %   M_               [structure]
0009 %
0010 % OUTPUTS
0011 %   none
0012 %
0013 % SPECIAL REQUIREMENTS
0014 %   none
0015 %
0016 % PARALLEL CONTEXT
0017 % See the comment in random_walk_metropolis_hastings.m funtion.
0018 
0019 % Copyright (C) 2005-2011 Dynare Team
0020 %
0021 % This file is part of Dynare.
0022 %
0023 % Dynare is free software: you can redistribute it and/or modify
0024 % it under the terms of the GNU General Public License as published by
0025 % the Free Software Foundation, either version 3 of the License, or
0026 % (at your option) any later version.
0027 %
0028 % Dynare is distributed in the hope that it will be useful,
0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0031 % GNU General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License
0034 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 DirectoryName = CheckPath('Output',M_.dname);
0037 MhDirectoryName = CheckPath('metropolis',M_.dname);
0038 
0039 TeX = options_.TeX;
0040 nblck = options_.mh_nblck;
0041 % Brooks and Gelman tests need more than one block
0042 if nblck == 1
0043     return;
0044 end
0045 npar = estim_params_.nvx;
0046 npar = npar + estim_params_.nvn;
0047 npar = npar + estim_params_.ncx;
0048 npar = npar + estim_params_.ncn;
0049 npar = npar + estim_params_.np ;
0050 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
0051 
0052 load([MhDirectoryName '/'  M_.fname '_mh_history.mat'])
0053 
0054 NumberOfMcFilesPerBlock = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blck1.mat']),1);
0055 for blck = 2:nblck
0056     tmp = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blck' int2str(blck) '.mat']),1);
0057     if tmp~=NumberOfMcFilesPerBlock
0058         disp(['McMCDiagnostics:: The number of mh files in chain ' int2str(blck) ' is ' int2str(tmp) ' while'])
0059         disp(['                  the number of mh files in chain 1 is ' int2str(mcfiles) '!'])
0060         error('The number of mh files has to be constant across chains!')
0061     end
0062 end
0063 
0064 PastDraws = sum(record.MhDraws,1);
0065 LastFileNumber = PastDraws(2);
0066 LastLineNumber = record.MhDraws(end,3);
0067 NumberOfDraws  = PastDraws(1);
0068 
0069 Origin = 1000;
0070 StepSize = ceil((NumberOfDraws-Origin)/100);% So that the computational time does not
0071 ALPHA = 0.2;                                % increase too much with the number of simulations.
0072 time = 1:NumberOfDraws;
0073 xx = Origin:StepSize:NumberOfDraws;
0074 NumberOfLines = length(xx);
0075 tmp = zeros(NumberOfDraws*nblck,3);
0076 UDIAG = zeros(NumberOfLines,6,npar);
0077 
0078 if NumberOfDraws < Origin
0079     disp('MCMC Diagnostics :: The number of simulations is to small to compute the MCMC convergence diagnostics.')
0080     return
0081 end
0082 
0083 if TeX
0084     fidTeX = fopen([DirectoryName '/' M_.fname '_UnivariateDiagnostics.TeX'],'w');
0085     fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n');
0086     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
0087     fprintf(fidTeX,' \n');
0088 end
0089 
0090 disp('MCMC Diagnostics: Univariate convergence diagnostic, Brooks and Gelman (1998):')
0091 
0092 % The mandatory variables for local/remote parallel
0093 % computing are stored in localVars struct.
0094 
0095 localVars.MhDirectoryName = MhDirectoryName;
0096 localVars.nblck = nblck;
0097 localVars.NumberOfMcFilesPerBlock = NumberOfMcFilesPerBlock;
0098 localVars.Origin = Origin;
0099 localVars.StepSize = StepSize;
0100 localVars.mh_drop = options_.mh_drop;
0101 localVars.NumberOfDraws = NumberOfDraws;
0102 localVars.NumberOfLines = NumberOfLines;
0103 localVars.time = time;
0104 localVars.M_ = M_;
0105 
0106 
0107 % Like sequential execution!
0108 if isnumeric(options_.parallel),
0109     fout = McMCDiagnostics_core(localVars,1,npar,0);
0110     UDIAG = fout.UDIAG;
0111     clear fout
0112     % Parallel execution!
0113 else
0114     ModelName = M_.fname;
0115     if ~isempty(M_.bvar)
0116         ModelName = [M_.fname '_bvar'];
0117     end
0118     NamFileInput={[M_.dname '/metropolis/'],[ModelName '_mh*_blck*.mat']};
0119     
0120     [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,NamFileInput,'McMCDiagnostics_core', localVars, [], options_.parallel_info);
0121     UDIAG = fout(1).UDIAG;
0122     for j=2:totCPU,
0123         UDIAG = cat(3,UDIAG ,fout(j).UDIAG);
0124     end
0125 end
0126 
0127 UDIAG(:,[2 4 6],:) = UDIAG(:,[2 4 6],:)/nblck;
0128 disp(' ')
0129 clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp;    
0130 pages = floor(npar/3);
0131 k = 0;  
0132 for i = 1:pages
0133     h=dyn_figure(options_,'Name','MCMC univariate diagnostic (Brooks and Gelman,1998)');
0134     boxplot = 1;
0135     for j = 1:3 % Loop over parameters
0136         k = k+1;
0137         [nam,namtex] = get_the_name(k,TeX,M_,estim_params_,options_);
0138         for crit = 1:3% Loop over criteria
0139             if crit == 1
0140                 plt1 = UDIAG(:,1,k);
0141                 plt2 = UDIAG(:,2,k);
0142                 namnam  = [nam , ' (Interval)']; 
0143             elseif crit == 2
0144                 plt1 = UDIAG(:,3,k);
0145                 plt2 = UDIAG(:,4,k);
0146                 namnam  = [nam , ' (m2)'];
0147             elseif crit == 3    
0148                 plt1 = UDIAG(:,5,k);
0149                 plt2 = UDIAG(:,6,k);
0150                 namnam  = [nam , ' (m3)'];
0151             end
0152             if TeX
0153                 if j==1
0154                     NAMES = deblank(namnam);
0155                     TEXNAMES = deblank(namtex);
0156                 else
0157                     NAMES = char(NAMES,deblank(namnam));
0158                     TEXNAMES = char(TEXNAMES,deblank(namtex));
0159                 end
0160             end
0161             subplot(3,3,boxplot);
0162             plot(xx,plt1,'-b');     % Pooled
0163             hold on;
0164             plot(xx,plt2,'-r');     % Within (mean)
0165             hold off;
0166             xlim([xx(1) xx(NumberOfLines)])
0167             title(namnam,'Interpreter','none')
0168             boxplot = boxplot + 1;
0169         end
0170     end
0171     dyn_saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(i)],options_);
0172     if TeX
0173         fprintf(fidTeX,'\\begin{figure}[H]\n');
0174         for jj = 1:size(NAMES,1)
0175             fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
0176         end    
0177         fprintf(fidTeX,'\\centering \n');
0178         fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_udiag%s}\n',M_.fname,int2str(i));
0179         fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n');
0180         fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n');
0181         fprintf(fidTeX,'the eighty percent interval, the second and third moments.}');
0182         fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(i));
0183         fprintf(fidTeX,'\\end{figure}\n');
0184         fprintf(fidTeX,'\n');
0185     end
0186 end
0187 reste = npar-k;
0188 if reste
0189     if reste == 1
0190         nr = 3;
0191         nc = 1;
0192     elseif reste == 2;
0193         nr = 2;
0194         nc = 3;
0195     end
0196     h = dyn_figure(options_,'Name','MCMC univariate diagnostic (Brooks and Gelman, 1998)');
0197     boxplot = 1;
0198     for j = 1:reste
0199         k = k+1;
0200         [nam,namtex] = get_the_name(k,TeX,M_,estim_params_,options_);
0201         for crit = 1:3
0202             if crit == 1
0203                 plt1 = UDIAG(:,1,k);
0204                 plt2 = UDIAG(:,2,k);
0205                 namnam  = [nam , ' (Interval)']; 
0206             elseif crit == 2
0207                 plt1 = UDIAG(:,3,k);
0208                 plt2 = UDIAG(:,4,k);
0209                 namnam  = [nam , ' (m2)'];
0210             elseif crit == 3    
0211                 plt1 = UDIAG(:,5,k);
0212                 plt2 = UDIAG(:,6,k);
0213                 namnam  = [nam , ' (m3)'];
0214             end
0215             if TeX
0216                 if j==1
0217                     NAMES = deblank(namnam);
0218                     TEXNAMES = deblank(namtex);
0219                 else
0220                     NAMES = char(NAMES,deblank(namnam));
0221                     TEXNAMES = char(TEXNAMES,deblank(namtex));
0222                 end
0223             end
0224             subplot(nr,nc,boxplot);
0225             plot(xx,plt1,'-b');                                 % Pooled
0226             hold on;
0227             plot(xx,plt2,'-r');                                 % Within (mean)
0228             hold off;
0229             xlim([xx(1) xx(NumberOfLines)]);
0230             title(namnam,'Interpreter','none');
0231             boxplot = boxplot + 1;
0232         end
0233     end
0234     dyn_saveas(h,[ DirectoryName '/' M_.fname '_udiag' int2str(pages+1)],options_);
0235     if TeX
0236         fprintf(fidTeX,'\\begin{figure}[H]\n');
0237         for jj = 1:size(NAMES,1);
0238             fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
0239         end    
0240         fprintf(fidTeX,'\\centering \n');
0241         fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_udiag%s}\n',M_.fname,int2str(pages+1));
0242         if reste == 2
0243             fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n');
0244             fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n');
0245             fprintf(fidTeX,'the eighty percent interval, the second and third moments.}');
0246         elseif reste == 1
0247             fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n');
0248             fprintf(fidTeX,'The first, second and third rows are respectively the criteria based on\n');
0249             fprintf(fidTeX,'the eighty percent interval, the second and third moments.}');
0250         end
0251         fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(pages+1));
0252         fprintf(fidTeX,'\\end{figure}\n');
0253         fprintf(fidTeX,'\n');
0254         fprintf(fidTeX,'% End Of TeX file.');
0255         fclose(fidTeX);
0256     end
0257 end % if reste > 0
0258 clear UDIAG;
0259 %%
0260 %% Multivariate diagnostic.
0261 %%
0262 if TeX
0263     fidTeX = fopen([DirectoryName '/' M_.fname '_MultivariateDiagnostics.TeX'],'w');
0264     fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n');
0265     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
0266     fprintf(fidTeX,' \n');
0267 end
0268 tmp = zeros(NumberOfDraws*nblck,3);
0269 MDIAG = zeros(NumberOfLines,6);
0270 for b = 1:nblck
0271     startline = 0;
0272     for n = 1:NumberOfMcFilesPerBlock
0273         %load([MhDirectoryName '/' mcfiles(n,1,b).name],'logpo2');
0274         load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) '.mat'],'logpo2');
0275         nlogpo2 = size(logpo2,1);
0276         tmp((b-1)*NumberOfDraws+startline+(1:nlogpo2),1) = logpo2;
0277         startline = startline+nlogpo2;
0278     end
0279 % $$$   %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'logpo2');
0280 % $$$   load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'logpo2');
0281 % $$$   tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+ MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = logpo2;
0282 end
0283 clear logpo2;
0284 tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
0285 tmp(:,3) = kron(ones(nblck,1),time'); 
0286 tmp = sortrows(tmp,1);
0287 ligne   = 0;
0288 for iter  = Origin:StepSize:NumberOfDraws
0289     ligne = ligne+1;
0290     linea = ceil(options_.mh_drop*iter);
0291     n     = iter-linea+1;
0292     cinf  = round(n*ALPHA/2);
0293     csup  = round(n*(1-ALPHA/2));
0294     CINF  = round(nblck*n*ALPHA/2);
0295     CSUP  = round(nblck*n*(1-ALPHA/2));
0296     temp  = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2);
0297     MDIAG(ligne,1) = temp(CSUP,1)-temp(CINF,1);
0298     moyenne = mean(temp(:,1));%% Pooled mean.
0299     MDIAG(ligne,3) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1);
0300     MDIAG(ligne,5) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1);
0301     for i=1:nblck
0302         pmet = temp(find(temp(:,2)==i));
0303         MDIAG(ligne,2) = MDIAG(ligne,2) + pmet(csup,1)-pmet(cinf,1);
0304         moyenne = mean(pmet,1); %% Within mean.
0305         MDIAG(ligne,4) = MDIAG(ligne,4) + sum((pmet(:,1)-moyenne).^2)/(n-1);
0306         MDIAG(ligne,6) = MDIAG(ligne,6) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1);
0307     end
0308 end
0309 MDIAG(:,[2 4 6],:) = MDIAG(:,[2 4 6],:)/nblck;  
0310 
0311 h = dyn_figure(options_,'Name','Multivariate diagnostic');
0312 boxplot = 1;
0313 for crit = 1:3
0314     if crit == 1
0315         plt1 = MDIAG(:,1);
0316         plt2 = MDIAG(:,2);
0317         namnam  = 'Interval'; 
0318     elseif crit == 2
0319         plt1 = MDIAG(:,3);
0320         plt2 = MDIAG(:,4);
0321         namnam  = 'm2';
0322     elseif crit == 3    
0323         plt1 = MDIAG(:,5);
0324         plt2 = MDIAG(:,6);
0325         namnam  = 'm3';
0326     end
0327     if TeX
0328         if crit == 1
0329             NAMES = deblank(namnam);
0330         else
0331             NAMES = char(NAMES,deblank(namnam));
0332         end
0333     end
0334     subplot(3,1,boxplot);
0335     plot(xx,plt1,'-b');  % Pooled
0336     hold on
0337     plot(xx,plt2,'-r');  % Within (mean)
0338     hold off
0339     xlim([xx(1) xx(NumberOfLines)])
0340     title(namnam,'Interpreter','none');
0341     boxplot = boxplot + 1;
0342 end
0343 dyn_saveas(h,[ DirectoryName '/' M_.fname '_mdiag'],options_);
0344 
0345 if TeX
0346     fprintf(fidTeX,'\\begin{figure}[H]\n');
0347     for jj = 1:3
0348         fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),' ');
0349     end    
0350     fprintf(fidTeX,'\\centering \n');
0351     fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_mdiag}\n',M_.fname);
0352     fprintf(fidTeX,'\\caption{Multivariate convergence diagnostics for the Metropolis-Hastings.\n');
0353     fprintf(fidTeX,'The first, second and third rows are respectively the criteria based on\n');
0354     fprintf(fidTeX,'the eighty percent interval, the second and third moments. The different \n');
0355     fprintf(fidTeX,'parameters are aggregated using the posterior kernel.}');
0356     fprintf(fidTeX,'\\label{Fig:MultivariateDiagnostics}\n');
0357     fprintf(fidTeX,'\\end{figure}\n');
0358     fprintf(fidTeX,'\n');
0359     fprintf(fidTeX,'% End Of TeX file.');
0360     fclose(fidTeX);
0361 end
McMCDiagnostics

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
