diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m
index 1333f7cf874dbab4ed71024451a1f91a10fdf737..4291d0dfa921e41bfae447f0f7125c26a1fb779a 100644
--- a/matlab/McMCDiagnostics.m
+++ b/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
\ No newline at end of file
+end
+
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 4da9c12e575abbcc206265f93dd6cf2175efe0a6..2eaad390824df5b56f0ff9ea6c343ea531262fdc 100644
--- a/matlab/dynare_estimation_1.m
+++ b/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_);
diff --git a/matlab/mcmc_ifac.m b/matlab/mcmc_ifac.m
new file mode 100644
index 0000000000000000000000000000000000000000..bc7be30c43e2a0a9aa11e6234730b48b04985f5e
--- /dev/null
+++ b/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);
diff --git a/matlab/missing/stats/autocorr.m b/matlab/missing/stats/autocorr.m
new file mode 100644
index 0000000000000000000000000000000000000000..5fb53337c1388d81fda3bc2b1af44894b4d76dc5
--- /dev/null
+++ b/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