Merge remote-tracking branch 'johannes/geweke_convergence' into geweke
commit
941cda7655
|
@ -4212,12 +4212,19 @@ graphs of smoothed shocks, smoothed observation errors, smoothed and historical
|
||||||
|
|
||||||
@algorithmshead
|
@algorithmshead
|
||||||
|
|
||||||
The Monte Carlo Markov Chain (MCMC) univariate diagnostics are generated
|
The Monte Carlo Markov Chain (MCMC) diagnostics are generated
|
||||||
by the estimation command if @ref{mh_nblocks} is larger than 1, if
|
by the estimation command if @ref{mh_replic} is larger than 2000 and if
|
||||||
@ref{mh_replic} is larger than 2000, and if option @ref{nodiagnostic} is
|
option @ref{nodiagnostic} is not used. If @ref{mh_nblocks} is equal to one,
|
||||||
not used. As described in section 3 of @cite{Brooks and Gelman (1998)}
|
the convergence diagnostics of @cite{Geweke (1992,1999)} is computed. It uses a
|
||||||
the convergence diagnostics are based on comparing pooled and within
|
chi square test to compare the means of the first and last draws specified in
|
||||||
MCMC moments (Dynare displays the second and third order moments, and
|
@ref{geweke_interval} (@pxref{geweke_interval}) after discarding the burnin of @ref{mh_drop}. The test is
|
||||||
|
computed using variance estimates under the assumption of no serial correlation
|
||||||
|
as well as using tapering windows specified in @ref{taper_steps} (@pxref{taper_steps}).
|
||||||
|
If @ref{mh_nblocks} is larger than 1, the convergence diagnostics of
|
||||||
|
@cite{Brooks and Gelman (1998)} are used instead.
|
||||||
|
As described in section 3 of @cite{Brooks and Gelman (1998)} the univariate
|
||||||
|
convergence diagnostics are based on comparing pooled and within MCMC moments
|
||||||
|
(Dynare displays the second and third order moments, and
|
||||||
the length of the Highest Probability Density interval covering 80% of
|
the length of the Highest Probability Density interval covering 80% of
|
||||||
the posterior distribution). Due to computational reasons, the
|
the posterior distribution). Due to computational reasons, the
|
||||||
multivariate convergence diagnostic does not follow @cite{Brooks and
|
multivariate convergence diagnostic does not follow @cite{Brooks and
|
||||||
|
@ -4358,8 +4365,8 @@ the total number of Metropolis draws available. Default:
|
||||||
@code{2}
|
@code{2}
|
||||||
|
|
||||||
@item mh_drop = @var{DOUBLE}
|
@item mh_drop = @var{DOUBLE}
|
||||||
The fraction of initially generated parameter vectors to be dropped
|
@anchor{mh_drop}
|
||||||
before using posterior simulations. Default: @code{0.5}
|
The fraction of initially generated parameter vectors to be dropped as a burnin before using posterior simulations. Default: @code{0.5}
|
||||||
|
|
||||||
@item mh_jscale = @var{DOUBLE}
|
@item mh_jscale = @var{DOUBLE}
|
||||||
The scale to be used for the jumping distribution in
|
The scale to be used for the jumping distribution in
|
||||||
|
@ -4758,6 +4765,18 @@ Value used to test if a generalized eigenvalue is 0/0 in the generalized
|
||||||
Schur decomposition (in which case the model does not admit a unique
|
Schur decomposition (in which case the model does not admit a unique
|
||||||
solution). Default: @code{1e-6}.
|
solution). Default: @code{1e-6}.
|
||||||
|
|
||||||
|
@item taper_steps = [@var{INTEGER1} @var{INTEGER2} @dots{}]
|
||||||
|
@anchor{taper_steps}
|
||||||
|
Percent tapering used for the spectral window in the @cite{Geweke (1992,1999)}
|
||||||
|
convergence diagnostics (requires @ref{mh_nblocks}=1). The tapering is used to
|
||||||
|
take the serial correlation of the posterior draws into account. Default: @code{[4 8 15]}.
|
||||||
|
|
||||||
|
@item geweke_interval = [@var{double} @var{double}]
|
||||||
|
@anchor{geweke_interval}
|
||||||
|
Percentage of MCMC draws at the beginning and end of the MCMC chain taken
|
||||||
|
to compute the @cite{Geweke (1992,1999)} convergence diagnostics (requires @ref{mh_nblocks}=1)
|
||||||
|
after discarding the first @ref{mh_drop} percent of draws as a burnin. Default: @code{[0.2 0.5]}.
|
||||||
|
|
||||||
@end table
|
@end table
|
||||||
|
|
||||||
@customhead{Note}
|
@customhead{Note}
|
||||||
|
@ -5048,6 +5067,56 @@ Upper/lower bound of the 90\% HPD interval taking into account both parameter an
|
||||||
|
|
||||||
@end defvr
|
@end defvr
|
||||||
|
|
||||||
|
@defvr {MATLAB/Octave variable} oo_.convergence.geweke
|
||||||
|
@anchor{convergence.geweke}
|
||||||
|
Variable set by the convergence diagnostics of the @code{estimation} command when used with @ref{mh_nblocks}=1 option (@pxref{mh_nblocks}).
|
||||||
|
|
||||||
|
Fields are of the form:
|
||||||
|
@example
|
||||||
|
@code{oo_.convergence.geweke.@var{VARIABLE_NAME}.@var{DIAGNOSTIC_OBJECT}}
|
||||||
|
@end example
|
||||||
|
where @var{DIAGNOSTIC_OBJECT} is one of the following:
|
||||||
|
|
||||||
|
@table @code
|
||||||
|
|
||||||
|
@item posteriormean
|
||||||
|
Mean of the posterior parameter distribution
|
||||||
|
|
||||||
|
@item posteriorstd
|
||||||
|
Standard deviation of the posterior parameter distribution
|
||||||
|
|
||||||
|
@item nse_iid
|
||||||
|
Numerical standard error (NSE) under the assumption of iid draws
|
||||||
|
|
||||||
|
@item rne_iid
|
||||||
|
Relative numerical efficiency (RNE) under the assumption of iid draws
|
||||||
|
|
||||||
|
@item nse_x
|
||||||
|
Numerical standard error (NSE) when using an x% taper
|
||||||
|
|
||||||
|
@item rne_x
|
||||||
|
Relative numerical efficiency (RNE) when using an x% taper
|
||||||
|
|
||||||
|
@item pooled_mean
|
||||||
|
Mean of the parameter when pooling the beginning and end parts of the chain
|
||||||
|
specified in @ref{geweke_interval} and weighting them with their relative precision.
|
||||||
|
It is a vector containing the results under the iid assumption followed by the ones
|
||||||
|
using the @ref{taper_steps} (@pxref{taper_steps}).
|
||||||
|
|
||||||
|
@item pooled_nse
|
||||||
|
NSE of the parameter when pooling the beginning and end parts of the chain and weighting them with their relative precision. See @code{pooled_mean}
|
||||||
|
|
||||||
|
@item prob_chi2_test
|
||||||
|
p-value of a chi squared test for equality of means in the beginning and the end
|
||||||
|
of the MCMC chain. See @code{pooled_mean}. A value above 0.05 indicates that
|
||||||
|
the null hypothesis of equal means and thus convergence cannot be rejected
|
||||||
|
at the 5 percent level. Differing values along the @ref{taper_steps} signal
|
||||||
|
the presence of significant autocorrelation in draws. In this case, the
|
||||||
|
estimates using a higher tapering are usually more reliable.
|
||||||
|
|
||||||
|
@end table
|
||||||
|
@end defvr
|
||||||
|
|
||||||
@deffn Command model_comparison @var{FILENAME}[(@var{DOUBLE})]@dots{};
|
@deffn Command model_comparison @var{FILENAME}[(@var{DOUBLE})]@dots{};
|
||||||
@deffnx Command model_comparison (marginal_density = laplace | modifiedharmonicmean) @var{FILENAME}[(@var{DOUBLE})]@dots{};
|
@deffnx Command model_comparison (marginal_density = laplace | modifiedharmonicmean) @var{FILENAME}[(@var{DOUBLE})]@dots{};
|
||||||
|
|
||||||
|
@ -8688,6 +8757,16 @@ Fernández-Villaverde, Jesús and Juan Rubio-Ramírez (2005): ``Estimating
|
||||||
Dynamic Equilibrium Economies: Linear versus Nonlinear Likelihood,''
|
Dynamic Equilibrium Economies: Linear versus Nonlinear Likelihood,''
|
||||||
@i{Journal of Applied Econometrics}, 20, 891--910
|
@i{Journal of Applied Econometrics}, 20, 891--910
|
||||||
|
|
||||||
|
@item
|
||||||
|
Geweke, John (1992): ``Evaluating the accuracy of sampling-based approaches
|
||||||
|
to the calculation of posterior moments'', in J.O. Berger, J.M. Bernardo,
|
||||||
|
A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of the Fourth Valencia
|
||||||
|
International Meeting on Bayesian Statistics, pp. 169--194, Oxford University Press
|
||||||
|
|
||||||
|
@item
|
||||||
|
Geweke, John (1999): ``Using simulation methods for Bayesian econometric models:
|
||||||
|
Inference, development and communication,'' @i{Econometric Reviews}, 18(1), 1--73
|
||||||
|
|
||||||
@item
|
@item
|
||||||
Ireland, Peter (2004): ``A Method for Taking Models to the Data,''
|
Ireland, Peter (2004): ``A Method for Taking Models to the Data,''
|
||||||
@i{Journal of Economic Dynamics and Control}, 28, 1205--26
|
@i{Journal of Economic Dynamics and Control}, 28, 1205--26
|
||||||
|
|
|
@ -1,4 +1,4 @@
|
||||||
function McMCDiagnostics(options_, estim_params_, M_)
|
function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
|
||||||
% function McMCDiagnostics
|
% function McMCDiagnostics
|
||||||
% Computes convergence tests
|
% Computes convergence tests
|
||||||
%
|
%
|
||||||
|
@ -8,7 +8,7 @@ function McMCDiagnostics(options_, estim_params_, M_)
|
||||||
% M_ [structure]
|
% M_ [structure]
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% none
|
% oo_ [structure]
|
||||||
%
|
%
|
||||||
% SPECIAL REQUIREMENTS
|
% SPECIAL REQUIREMENTS
|
||||||
% none
|
% none
|
||||||
|
@ -38,10 +38,6 @@ MhDirectoryName = CheckPath('metropolis',M_.dname);
|
||||||
|
|
||||||
TeX = options_.TeX;
|
TeX = options_.TeX;
|
||||||
nblck = options_.mh_nblck;
|
nblck = options_.mh_nblck;
|
||||||
% Brooks and Gelman tests need more than one block
|
|
||||||
if nblck == 1
|
|
||||||
return;
|
|
||||||
end
|
|
||||||
npar = estim_params_.nvx;
|
npar = estim_params_.nvx;
|
||||||
npar = npar + estim_params_.nvn;
|
npar = npar + estim_params_.nvn;
|
||||||
npar = npar + estim_params_.ncx;
|
npar = npar + estim_params_.ncx;
|
||||||
|
@ -56,6 +52,57 @@ NumberOfMcFilesPerBlock = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blc
|
||||||
% check if all previous files are there for block 1
|
% check if all previous files are there for block 1
|
||||||
check_presence_consecutive_MC_files(MhDirectoryName,M_.fname,1)
|
check_presence_consecutive_MC_files(MhDirectoryName,M_.fname,1)
|
||||||
|
|
||||||
|
if nblck == 1 % Brooks and Gelman tests need more than one block
|
||||||
|
convergence_diagnostics_geweke=zeros(npar,4+2*length(options_.convergence.geweke.taper_steps));
|
||||||
|
if any(options_.convergence.geweke.geweke_interval<0) || any(options_.convergence.geweke.geweke_interval>1) || length(any(options_.convergence.geweke.geweke_interval<0))~=2 ...
|
||||||
|
|| (options_.convergence.geweke.geweke_interval(2)-options_.convergence.geweke.geweke_interval(1)<0)
|
||||||
|
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
|
||||||
|
first_obs_begin_sample = max(1,ceil(options_.mh_drop*options_.mh_replic));
|
||||||
|
last_obs_begin_sample = first_obs_begin_sample+round(options_.convergence.geweke.geweke_interval(1)*options_.mh_replic*options_.mh_drop);
|
||||||
|
first_obs_end_sample = first_obs_begin_sample+round(options_.convergence.geweke.geweke_interval(2)*options_.mh_replic*options_.mh_drop);
|
||||||
|
param_name=[];
|
||||||
|
for jj=1:npar
|
||||||
|
param_name = strvcat(param_name,get_the_name(jj,options_.TeX,M_,estim_params_,options_));
|
||||||
|
end
|
||||||
|
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,options_.mh_replic);
|
||||||
|
fprintf('p-values are for Chi2-test for equality of means.\n');
|
||||||
|
Geweke_header={'Parameter', 'Post. Mean', 'Post. Std', 'p-val No Taper'};
|
||||||
|
print_string=['%',num2str(size(param_name,2)+3),'s \t %12.3f \t %12.3f \t %12.3f'];
|
||||||
|
print_string_header=['%',num2str(size(param_name,2)+3),'s \t %12s \t %12s \t %12s'];
|
||||||
|
for ii=1:length(options_.convergence.geweke.taper_steps)
|
||||||
|
Geweke_header=[Geweke_header, ['p-val ' num2str(options_.convergence.geweke.taper_steps(ii)),'% Taper']];
|
||||||
|
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,Geweke_header{1,:});
|
||||||
|
for jj=1:npar
|
||||||
|
startline=0;
|
||||||
|
for n = 1:NumberOfMcFilesPerBlock
|
||||||
|
load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck1.mat'],'x2');
|
||||||
|
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;
|
||||||
|
|
||||||
|
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_);
|
||||||
|
|
||||||
|
results_struct = geweke_chi2_test(results_vec1,results_vec2,results_struct,options_);
|
||||||
|
eval(['oo_.convergence.geweke.',param_name(jj,:),'=results_struct;'])
|
||||||
|
fprintf(print_string,param_name(jj,:),results_struct.posteriormean,results_struct.posteriorstd,results_struct.prob_chi2_test)
|
||||||
|
end
|
||||||
|
skipline(2);
|
||||||
|
return;
|
||||||
|
end
|
||||||
|
|
||||||
for blck = 2:nblck
|
for blck = 2:nblck
|
||||||
tmp = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blck' int2str(blck) '.mat']),1);
|
tmp = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blck' int2str(blck) '.mat']),1);
|
||||||
if tmp~=NumberOfMcFilesPerBlock
|
if tmp~=NumberOfMcFilesPerBlock
|
||||||
|
|
|
@ -611,8 +611,8 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
|
||||||
CutSample(M_, options_, estim_params_);
|
CutSample(M_, options_, estim_params_);
|
||||||
return
|
return
|
||||||
else
|
else
|
||||||
if ~options_.nodiagnostic && options_.mh_replic > 2000 && options_.mh_nblck > 1
|
if ~options_.nodiagnostic && options_.mh_replic > 2000
|
||||||
McMCDiagnostics(options_, estim_params_, M_);
|
oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
|
||||||
end
|
end
|
||||||
%% Here i discard first half of the draws:
|
%% Here i discard first half of the draws:
|
||||||
CutSample(M_, options_, estim_params_);
|
CutSample(M_, options_, estim_params_);
|
||||||
|
|
|
@ -0,0 +1,74 @@
|
||||||
|
function results_struct = geweke_chi2_test(results1,results2,results_struct,options)
|
||||||
|
% results_struct = geweke_chi2_test(results1,results2,results_struct,options)
|
||||||
|
% PURPOSE: computes Geweke's chi-squared test for two sets of MCMC sample draws
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% results1 [1 by (4+n_taper*2) vector] vector with post. mean,
|
||||||
|
% std, NSE_iid, RNE_iid, and tapered NSE and RNE
|
||||||
|
% for chain part 1
|
||||||
|
% results2 [1 by (4+n_taper*2) vector] vector with post. mean,
|
||||||
|
% std, NSE_iid, RNE_iid, and tapered NSE and RNE
|
||||||
|
% for chain part 2
|
||||||
|
% results_struct [structure] results structure generated by geweke_moments
|
||||||
|
% Dynareoptions [structure]
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% results_struct [structure] containing the following fields:
|
||||||
|
% pooled_mean Pooled mean of the chain parts, weighted
|
||||||
|
% with precision
|
||||||
|
% rpooled_nse Pooled NSE of the chain parts, weighted
|
||||||
|
% with precision
|
||||||
|
% prob_chi2_test p-value of Chi2-test for equality of
|
||||||
|
% means in both chain parts
|
||||||
|
% -----------------------------------------------------------------
|
||||||
|
|
||||||
|
%
|
||||||
|
% SPECIAL REQUIREMENTS
|
||||||
|
% None.
|
||||||
|
|
||||||
|
% Copyright (C) 2013 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/>.
|
||||||
|
%
|
||||||
|
% ------------------------------------------------
|
||||||
|
% REFERENCES: Geweke (1992), `Evaluating the accuracy of sampling-based
|
||||||
|
% approaches to the calculation of posterior moments', in J.O. Berger,
|
||||||
|
% J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of
|
||||||
|
% the Fourth Valencia International Meeting on Bayesian Statistics,
|
||||||
|
% pp. 169-194, Oxford University Press
|
||||||
|
% Geweke (1999): `Using simulation methods for Bayesian econometric models:
|
||||||
|
% Inference, development and communication', Econometric Reviews, 18(1),
|
||||||
|
% 1-73
|
||||||
|
|
||||||
|
% written by: Johannes Pfeifer,
|
||||||
|
% based on code by James P. LeSage, who in turn
|
||||||
|
% drew on MATLAB programs written by Siddartha Chib
|
||||||
|
|
||||||
|
for k=1:length(options.convergence.geweke.taper_steps)+1;
|
||||||
|
NSE=[results1(:,3+(k-1)*2) results2(:,3+(k-1)*2)];
|
||||||
|
means=[results1(:,1) results2(:,1)];
|
||||||
|
diff_Means=means(:,1)-means(:,2);
|
||||||
|
sum_of_weights=sum(1./(NSE.^2),2);
|
||||||
|
pooled_mean=sum(means./(NSE.^2),2)./sum_of_weights;
|
||||||
|
pooled_NSE=1./sqrt(sum_of_weights);
|
||||||
|
|
||||||
|
test_stat=diff_Means.^2./sum(NSE.^2,2);
|
||||||
|
p = 1-chi2cdf(test_stat,1);
|
||||||
|
results_struct.pooled_mean(:,k) = pooled_mean;
|
||||||
|
results_struct.pooled_nse(:,k) = pooled_NSE;
|
||||||
|
results_struct.prob_chi2_test(:,k) = p;
|
||||||
|
end;
|
||||||
|
|
|
@ -0,0 +1,109 @@
|
||||||
|
function [results_vec, results_struct] = geweke_moments(draws,Dynareoptions)
|
||||||
|
%[results_vec, results_struct] = geweke_moments(draws,Dynareoptions)
|
||||||
|
% PURPOSE: computes Gewke's convergence diagnostics NSE and RNE
|
||||||
|
% (numerical std error and relative numerical efficiencies)
|
||||||
|
|
||||||
|
% INPUTS
|
||||||
|
% draws [ndraws by 1 vector]
|
||||||
|
% Dynareoptions [structure]
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% results_vec
|
||||||
|
% results_struct [structure] containing the following fields:
|
||||||
|
% posteriormean= posterior parameter mean
|
||||||
|
% posteriorstd = posterior standard deviation
|
||||||
|
% nse_iid = nse assuming no serial correlation for variable i
|
||||||
|
% rne_iid = rne assuming no serial correlation for variable i
|
||||||
|
% nse_x = nse using x% autocovariance tapered estimate
|
||||||
|
% rne_x = rne using x% autocovariance taper
|
||||||
|
% -----------------------------------------------------------------
|
||||||
|
|
||||||
|
%
|
||||||
|
% SPECIAL REQUIREMENTS
|
||||||
|
% None.
|
||||||
|
|
||||||
|
% Copyright (C) 2013 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/>.
|
||||||
|
|
||||||
|
% REFERENCES: Geweke (1992), `Evaluating the accuracy of sampling-based
|
||||||
|
% approaches to the calculation of posterior moments', in J.O. Berger,
|
||||||
|
% J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of
|
||||||
|
% the Fourth Valencia International Meeting on Bayesian Statistics,
|
||||||
|
% pp. 169-194, Oxford University Press
|
||||||
|
% Geweke (1999): `Using simulation methods for Bayesian econometric models:
|
||||||
|
% Inference, development and communication', Econometric Reviews, 18(1),
|
||||||
|
% 1-73
|
||||||
|
% -----------------------------------------------------------------
|
||||||
|
|
||||||
|
% written by: Johannes Pfeifer,
|
||||||
|
% based on code by James P. LeSage, who in turn
|
||||||
|
% drew on MATLAB programs written by Siddartha Chib
|
||||||
|
|
||||||
|
|
||||||
|
ndraw = size(draws,1);
|
||||||
|
n_groups=100;
|
||||||
|
taper_steps=Dynareoptions.convergence.geweke.taper_steps;
|
||||||
|
results_vec=zeros(1,4+2*length(taper_steps));
|
||||||
|
|
||||||
|
ns = floor(ndraw/n_groups); %step_size
|
||||||
|
n_draws_used = ns*n_groups; %effective number of draws used after rounding down
|
||||||
|
|
||||||
|
window_means= zeros(n_groups,1);
|
||||||
|
window_uncentered_variances= zeros(n_groups,1);
|
||||||
|
for ig=1:n_groups;
|
||||||
|
window_means(ig,1)=sum(draws((ig-1)*ns+1:ig*ns,1))/ns;
|
||||||
|
window_uncentered_variances(ig,1)=sum(draws((ig-1)*ns+1:ig*ns,1).^2)/ns;
|
||||||
|
end; %for ig
|
||||||
|
total_mean=mean(window_means);
|
||||||
|
total_variance=mean(window_uncentered_variances)-total_mean^2;
|
||||||
|
|
||||||
|
% save posterior means and std deviations to results_struct structure
|
||||||
|
results_vec(1,1)=total_mean;
|
||||||
|
results_vec(1,2)=sqrt(total_variance);
|
||||||
|
results_struct.posteriormean = total_mean;
|
||||||
|
results_struct.posteriorstd = results_vec(1,2);
|
||||||
|
|
||||||
|
% numerical standard error assuming no serial correlation
|
||||||
|
NSE=std(draws(1:n_draws_used,1),1)/sqrt(n_draws_used);
|
||||||
|
% save to results_struct structure
|
||||||
|
results_vec(1,3)=NSE;
|
||||||
|
results_vec(1,4)=total_variance/(n_draws_used*NSE^2);
|
||||||
|
results_struct.nse_iid = NSE;
|
||||||
|
results_struct.rne_iid = results_vec(1,4);
|
||||||
|
|
||||||
|
%get autocovariance of grouped means
|
||||||
|
centered_window_means=window_means-total_mean;
|
||||||
|
autocov_grouped_means=zeros(n_groups,1);
|
||||||
|
for lag=0:n_groups-1;
|
||||||
|
autocov_grouped_means(lag+1)=centered_window_means(lag+1:n_groups,1)'*centered_window_means(1:n_groups-lag,1)/100;
|
||||||
|
end;
|
||||||
|
|
||||||
|
% numerical standard error with tapered autocovariance functions
|
||||||
|
for taper_index=1:length(taper_steps)
|
||||||
|
taper=taper_steps(taper_index);
|
||||||
|
taper_lags=(1:taper-1)';
|
||||||
|
taper_lag_weight=1-taper_lags/taper;
|
||||||
|
tapered_sum_of_covariances=autocov_grouped_means(1)+sum(2*taper_lag_weight.*autocov_grouped_means(1+taper_lags));
|
||||||
|
NSE_taper=sqrt(tapered_sum_of_covariances/n_groups);
|
||||||
|
% save results_struct in structure
|
||||||
|
results_vec(1,4+taper_index*2-1)=NSE_taper;
|
||||||
|
results_vec(1,4+taper_index*2)=total_variance/(n_draws_used*NSE_taper^2);
|
||||||
|
|
||||||
|
eval(['results_struct.nse_taper_',num2str(taper),'= NSE_taper;']);
|
||||||
|
eval(['results_struct.rne_taper_',num2str(taper),'= total_variance/(n_draws_used*NSE_taper^2);']);
|
||||||
|
end; % end of for mm loop
|
||||||
|
|
|
@ -577,6 +577,10 @@ options_.osr.verbose=2;
|
||||||
% use GPU
|
% use GPU
|
||||||
options_.gpu = 0;
|
options_.gpu = 0;
|
||||||
|
|
||||||
|
%Geweke convergence diagnostics
|
||||||
|
options_.convergence.geweke.taper_steps=[4 8 15];
|
||||||
|
options_.convergence.geweke.geweke_interval=[0.2 0.5];
|
||||||
|
|
||||||
% initialize persistent variables in priordens()
|
% initialize persistent variables in priordens()
|
||||||
priordens([],[],[],[],[],[],1);
|
priordens([],[],[],[],[],[],1);
|
||||||
% initialize persistent variables in dyn_first_order_solver()
|
% initialize persistent variables in dyn_first_order_solver()
|
||||||
|
|
Loading…
Reference in New Issue