Added function to compute of inefficiency factors, used in McMCDiagnostics.m

Requires missing autocorr (available in some MATLAB toolboxes).
Requires factor out CutSample in dynare_estimation_1.m

matlab/McMCDiagnostics.m

@@ -69,6 +69,46 @@ if issue_an_error_message
     error('Estimation::mcmc::diagnostics: I cannot proceed because some MCMC files are missing. Check your MCMC files...')
 end
 
+% 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);
+
+param_name=[];
+for jj=1:npar   
+    param_name = strvcat(param_name,get_the_name(jj,options_.TeX,M_,estim_params_,options_));
+end
+fprintf('\nMCMC Inefficiency factors per block.\n');
+IFAC_header={'Parameter'};
+print_string=['%',num2str(size(param_name,2)+3),'s'];
+print_string_header=['%',num2str(size(param_name,2)+3),'s'];
+for j=1:nblck,
+    IFAC_header=[IFAC_header, {['block ' int2str(j)]}];
+    print_string=[print_string,' \t %12.3f'];
+    print_string_header=[print_string_header,' \t %12s'];
+end
+print_string=[print_string,'\n'];
+print_string_header=[print_string_header,'\n'];
+fprintf(print_string_header,IFAC_header{1,:});
+
+for jj=1:npar   
+    Draws = GetAllPosteriorDraws(jj,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+    Draws = reshape(Draws,[NumberOfDraws nblck]);
+    Nc = min(1000, NumberOfDraws/2);
+    for ll=1:nblck
+        Ifac(ll,jj) = mcmc_ifac(Draws(:,ll), Nc);
+    end
+    tmp = num2cell(Ifac(:,jj));
+    fprintf(print_string,param_name(jj,:),tmp{:})
+end
+skipline()
+record.InefficiencyFactorsPerBlock = Ifac;
+update_last_mh_history_file(MetropolisFolder, ModelName, record);
+
+
 PastDraws = sum(record.MhDraws,1);
 LastFileNumber = PastDraws(2);
 LastLineNumber = record.MhDraws(end,3);
@@ -436,4 +476,5 @@ if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
     fprintf(fidTeX,'\n');
     fprintf(fidTeX,'% End Of TeX file.');
     fclose(fidTeX);
-end
+end
+

matlab/dynare_estimation_1.m

@@ -427,15 +427,15 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
         end
         options_.analytic_derivation = ana_deriv_old;
     end
+    %% Here I discard first mh_drop percent of the draws:
+    CutSample(M_, options_, estim_params_);
     if options_.mh_posterior_mode_estimation
-        CutSample(M_, options_, estim_params_);
         return
     else
         if ~options_.nodiagnostic && options_.mh_replic>0
             oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
         end
-        %% Here I discard first mh_drop percent of the draws:
-        CutSample(M_, options_, estim_params_);
+
         %% Estimation of the marginal density from the Mh draws:
         if options_.mh_replic
             [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);

matlab/mcmc_ifac.m

@@ -0,0 +1,44 @@
+function Ifac = mcmc_ifac(X, Nc)
+% function Ifac = mcmc_ifac(X, Nc)
+% Compute inefficiency factor of a MCMC sample X
+%
+% INPUTS
+%   X:       time series
+%   Nc:      # of lags
+%
+% OUTPUTS
+%   Ifac:       inefficiency factor of MCMC sample
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2015 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/>.
+
+Nc = min(Nc, length(X)/2);
+AcorrXSIM = autocorr(X(:), Nc);
+%
+%Calculate the Parzen Weight
+Parzen=zeros(Nc+1,1);
+for i=1: Nc/2+1
+    Parzen(i)=1 - 6*(i/Nc)^2+ 6*(i/Nc)^3;
+end
+for i=(Nc/2)+1: Nc+1
+    Parzen(i)=2 * (1-(i/Nc))^3;
+end
+Parzen=Parzen';
+Ifac= 1+2*sum(Parzen(:).* AcorrXSIM);

matlab/missing/stats/autocorr.m

@@ -0,0 +1,40 @@
+function acf = autocorr(y, ar)
+% function acf = autocorr(y, ar)
+% autocorrelation function of y
+%
+% INPUTS
+%   y:       time series
+%   ar:      # of lags
+%
+% OUTPUTS
+%   acf:       autocorrelation for lags 1 to ar
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2015 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/>.
+
+
+y=y(:);
+acf = NaN(ar+1,1);
+acf(1)=1;
+m = mean(y);
+sd = std(y,1);
+for i=1:ar,
+    acf(i+1) = (y(i+1:end)-m)'*(y(1:end-i)-m)./((size(y,1))*sd^2);
+end