diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m index 1333f7cf8..4291d0dfa 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 aa8aad754..4f408a4de 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 000000000..bc7be30c4 --- /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 . + +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 000000000..5fb53337c --- /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 . + + +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