Merge pull request #1038 from JohannesPfeifer/kernel_density

Allow computing posterior kernel density
time-shift
Stéphane Adjemian 2015-10-12 23:23:07 +02:00
commit d3eca9338b
15 changed files with 227 additions and 149 deletions

View File

@ -5436,6 +5436,9 @@ controlled by the @code{ar} option.
@xref{contemporaneous_correlation}. Results are stored in @code{oo_.PosteriorTheoreticalMoments}.
Note that the @code{nocorr}-option has no effect.
@item no_posterior_kernel_density
Shuts off the computation of the kernel density estimator for the posterior objects (@pxref{density-field}).
@item conditional_variance_decomposition = @var{INTEGER}
See below.
@ -5814,8 +5817,10 @@ Variance of the posterior distribution
Deciles of the distribution.
@item density
Non parametric estimate of the posterior density. First and second
columns are respectively abscissa and ordinate coordinates.
@anchor{density-field}
Non parametric estimate of the posterior density following the approach outlined in @cite{Skoeld and Roberts (2003)}. First and second
columns are respectively abscissa and ordinate coordinates.
@end table
@item ESTIMATED_OBJECT
@ -13124,6 +13129,10 @@ Sims, Christopher A., Daniel F. Waggoner and Tao Zha (2008): ``Methods for
inference in large multiple-equation Markov-switching models,''
@i{Journal of Econometrics}, 146, 255--274
@item
Skoeld, Martin and Gareth O. Roberts (2003): ``Density Estimation for the
Metropolis-Hastings Algorithm,'' @i{Scandinavian Journal of Statistics}, 30, 699--718
@item
Smets, Frank and Rafael Wouters (2003): ``An Estimated Dynamic
Stochastic General Equilibrium Model of the Euro Area,'' @i{Journal of

View File

@ -389,6 +389,6 @@ oo.posterior_density.(type).(name) = density;
function [post_mean,hpd_interval,post_var] = Extractoo(oo,name,type)
hpd_interval = zeros(2,1);
post_mean = oo.posterior_mean.(type).(name);
hpd_interval(1) = oo.posterior_hpdinf.(type).(name);
hpd_interval(1) = oo.posterior_hpdinf.(type).(name);
hpd_interval(2) = oo.posterior_hpdsup.(type).(name);
post_var = oo.posterior_variance.(type).(name);

View File

@ -304,12 +304,12 @@ clear STOCK_IRF_DSGE;
for i = irf_shocks_indx
for j = 1:nvar
name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.deciles.' name ' = DistribIRF(:,:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
oo_.PosteriorIRF.dsge.Mean.(name) = MeanIRF(:,j,i);
oo_.PosteriorIRF.dsge.Median.(name) = MedianIRF(:,j,i);
oo_.PosteriorIRF.dsge.Var.(name) = VarIRF(:,j,i);
oo_.PosteriorIRF.dsge.deciles.(name) = DistribIRF(:,:,j,i);
oo_.PosteriorIRF.dsge.HPDinf.(name) = HPDIRF(:,1,j,i);
oo_.PosteriorIRF.dsge.HPDsup.(name) = HPDIRF(:,2,j,i);
end
end
@ -341,12 +341,12 @@ if MAX_nirfs_dsgevar
for i = irf_shocks_indx
for j = 1:nvar
name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
eval(['oo_.PosteriorIRF.bvardsge.Mean.' name ' = MeanIRFdsgevar(:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.Median.' name ' = MedianIRFdsgevar(:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.Var.' name ' = VarIRFdsgevar(:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.deciles.' name ' = DistribIRFdsgevar(:,:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.HPDinf.' name ' = HPDIRFdsgevar(:,1,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.HPDsup.' name ' = HPDIRFdsgevar(:,2,j,i);']);
oo_.PosteriorIRF.bvardsge.Mean.(name) = MeanIRFdsgevar(:,j,i);
oo_.PosteriorIRF.bvardsge.Median.(name) = MedianIRFdsgevar(:,j,i);
oo_.PosteriorIRF.bvardsge.Var.(name) = VarIRFdsgevar(:,j,i);
oo_.PosteriorIRF.bvardsge.deciles.(name) = DistribIRFdsgevar(:,:,j,i);
oo_.PosteriorIRF.bvardsge.HPDinf.(name) = HPDIRFdsgevar(:,1,j,i);
oo_.PosteriorIRF.bvardsge.HPDsup.(name) = HPDIRFdsgevar(:,2,j,i);
end
end
end

View File

@ -1,5 +1,5 @@
function oo_ = ...
conditional_variance_decomposition_mc_analysis(NumberOfSimulations, type, dname, fname, Steps, exonames, exo, var_list, endogenous_variable_index, mh_conf_sig, oo_)
conditional_variance_decomposition_mc_analysis(NumberOfSimulations, type, dname, fname, Steps, exonames, exo, var_list, endogenous_variable_index, mh_conf_sig, oo_,options_)
% This function analyses the (posterior or prior) distribution of the
% endogenous variables' conditional variance decomposition.
%
@ -62,11 +62,11 @@ end
name = [ var_list(endogenous_variable_index,:) '.' exo ];
if isfield(oo_, [ TYPE 'TheoreticalMoments' ])
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;'])
temporary_structure = oo_.([TYPE 'TheoreticalMoments']);
if isfield(temporary_structure,'dsge')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;'])
temporary_structure = oo_.([TYPE 'TheoreticalMoments']).dsge;
if isfield(temporary_structure,'ConditionalVarianceDecomposition')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean;'])
temporary_structure = oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Mean;
if isfield(temporary_structure,name)
if sum(Steps-temporary_structure.(name)(1,:)) == 0
% Nothing (new) to do here...
@ -91,25 +91,34 @@ p_mean = NaN(1,length(Steps));
p_median = NaN(1,length(Steps));
p_variance = NaN(1,length(Steps));
p_deciles = NaN(9,length(Steps));
p_density = NaN(2^9,2,length(Steps));
if options_.estimation.moments_posterior_density.indicator
p_density = NaN(2^9,2,length(Steps));
end
p_hpdinf = NaN(1,length(Steps));
p_hpdsup = NaN(1,length(Steps));
for i=1:length(Steps)
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ...
posterior_moments(tmp(:,i),1,mh_conf_sig);
if options_.estimation.moments_posterior_density.indicator
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ...
posterior_moments(tmp(:,i),1,mh_conf_sig);
p_density(:,:,i) = pp_density;
else
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ...
posterior_moments(tmp(:,i),0,mh_conf_sig);
end
p_mean(i) = pp_mean;
p_median(i) = pp_median;
p_variance(i) = pp_var;
p_deciles(:,i) = pp_deciles;
p_hpdinf(i) = hpd_interval(1);
p_hpdsup(i) = hpd_interval(2);
p_density(:,:,i) = pp_density;
end
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Steps = Steps;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean.' name ' = p_mean;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Median.' name ' = p_median;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Variance.' name ' = p_variance;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.HPDinf.' name ' = p_hpdinf;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.HPDsup.' name ' = p_hpdsup;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.deciles.' name ' = p_deciles;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.density.' name ' = p_density;']);
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Steps = Steps;
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Mean.(name) = p_mean;
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Median.(name) = p_median;
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Variance.(name) = p_variance;
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.HPDinf.(name) = p_hpdinf;
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.HPDsup.(name) = p_hpdsup;
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.deciles.(name) = p_deciles;
if options_.estimation.moments_posterior_density.indicator
oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.density.(name) = p_density;
end

View File

@ -44,15 +44,15 @@ else
end
if isfield(oo_,[TYPE 'TheoreticalMoments'])
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
if isfield(temporary_structure,'dsge')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge;
if isfield(temporary_structure,'correlation')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.correlation.Mean;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge.correlation.Mean;
if isfield(temporary_structure,deblank(var1))
eval(['temporary_structure_1 = oo_.' TYPE 'TheoreticalMoments.dsge.correlation.Mean.' var1 ';'])
temporary_structure_1 = oo_.([TYPE, 'TheoreticalMoments']).dsge.correlation.Mean.(var1);
if isfield(temporary_structure_1,deblank(var2))
eval(['temporary_structure_2 = temporary_structure_1.' var2 ';'])
temporary_structure_2 = temporary_structure_1.(var2);
l1 = length(temporary_structure_2);
if l1<nar
% INITIALIZATION:
@ -67,19 +67,19 @@ if isfield(oo_,[TYPE 'TheoreticalMoments'])
end
end
else
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_);
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_);
end
else
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_);
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_);
end
else
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_);
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_);
end
else
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_);
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_);
end
else
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_);
oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_);
end
ListOfFiles = dir([ PATH fname '_' TYPE 'Correlations*.mat']);
i1 = 1; tmp = zeros(SampleSize,1);
@ -91,12 +91,17 @@ for file = 1:length(ListOfFiles)
end
name = [ var1 '.' var2 ];
if ~isconst(tmp)
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp,0,mh_conf_sig);
end
if isfield(oo_,[ TYPE 'TheoreticalMoments'])
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
if isfield(temporary_structure,'dsge')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge;
if isfield(temporary_structure,'correlation')
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'Mean',nar,p_mean);
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'Median',nar,p_median);
@ -104,15 +109,17 @@ if ~isconst(tmp)
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDinf',nar,hpd_interval(1));
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDsup',nar,hpd_interval(2));
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'deciles',nar,p_deciles);
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,density);
if options_.estimation.moments_posterior_density.indicator
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,density);
end
end
end
end
else
if isfield(oo_,'PosteriorTheoreticalMoments')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
if isfield(temporary_structure,'dsge')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge;
if isfield(temporary_structure,'correlation')
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'Mean',nar,NaN);
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'Median',nar,NaN);
@ -120,33 +127,37 @@ else
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDinf',nar,NaN);
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDsup',nar,NaN);
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'deciles',nar,NaN);
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,NaN);
if options_.estimation.moments_posterior_density.indicator
oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,NaN);
end
end
end
end
end
function oo_ = initialize_output_structure(var1,var2,nar,type,oo_)
name = [ var1 '.' var2 ];
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.Mean.' name ' = NaN(' int2str(nar) ',1);']);
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.Median.' name ' = NaN(' int2str(nar) ',1);']);
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.Variance.' name ' = NaN(' int2str(nar) ',1);']);
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.HPDinf.' name ' = NaN(' int2str(nar) ',1);']);
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.HPDsup.' name ' = NaN(' int2str(nar) ',1);']);
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.deciles.' name ' = cell(' int2str(nar) ',1);']);
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.density.' name ' = cell(' int2str(nar) ',1);']);
function oo_ = initialize_output_structure(var1,var2,nar,type,oo_,options_)
oo_.([type, 'TheoreticalMoments']).dsge.correlation.Mean.(var1).(var2) = NaN(nar,1);
oo_.([type, 'TheoreticalMoments']).dsge.correlation.Median.(var1).(var2) = NaN(nar,1);
oo_.([type, 'TheoreticalMoments']).dsge.correlation.Variance.(var1).(var2) = NaN(nar,1);
oo_.([type, 'TheoreticalMoments']).dsge.correlation.HPDinf.(var1).(var2) = NaN(nar,1);
oo_.([type, 'TheoreticalMoments']).dsge.correlation.HPDsup.(var1).(var2) = NaN(nar,1);
oo_.([type, 'TheoreticalMoments']).dsge.correlation.deciles.(var1).(var2) = cell(nar,1);
if options_.estimation.moments_posterior_density.indicator
oo_.([type, 'TheoreticalMoments']).dsge.correlation.density.(var1).(var2) = cell(nar,1);
end
for i=1:nar
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.density.' name '(' int2str(i) ',1) = {NaN};']);
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.deciles.' name '(' int2str(i) ',1) = {NaN};']);
if options_.estimation.moments_posterior_density.indicator
oo_.([type, 'TheoreticalMoments']).dsge.correlation.density.(var1).(var2)(i,1) = {NaN};
end
oo_.([type, 'TheoreticalMoments']).dsge.correlation.deciles.(var1).(var2)(i,1) = {NaN};
end
function oo_ = fill_output_structure(var1,var2,type,oo_,moment,lag,result)
name = [ var1 '.' var2 ];
switch moment
case {'Mean','Median','Variance','HPDinf','HPDsup'}
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.' moment '.' name '(' int2str(lag) ',1) = result;']);
oo_.([type, 'TheoreticalMoments']).dsge.correlation.(moment).(var1).(var2)(lag,1) = result;
case {'deciles','density'}
eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.' moment '.' name '(' int2str(lag) ',1) = {result};']);
oo_.([type, 'TheoreticalMoments']).dsge.correlation.(moment).(var1).(var2)(lag,1) = {result};
otherwise
disp('fill_output_structure:: Unknown field!')
end

View File

@ -61,20 +61,20 @@ else
end
if isfield(oo_,[ TYPE 'TheoreticalMoments'])
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
if isfield(temporary_structure,'dsge')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge;
if isfield(temporary_structure,'covariance')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Mean;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Mean;
if isfield(temporary_structure,var1)
eval(['temporary_structure_1 = oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Mean.' var1 ';'])
temporary_structure_1 = oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Mean.(var1);
if isfield(temporary_structure_1,var2)
% Nothing to do (the covariance matrix is symmetric!).
return
end
else
if isfield(temporary_structure,var2)
eval(['temporary_structure_2 = oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Mean.' var2 ';'])
temporary_structure_2 = oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Mean.(var2);
if isfield(temporary_structure_2,var1)
% Nothing to do (the covariance matrix is symmetric!).
return
@ -104,35 +104,47 @@ for file = 1:length(ListOfFiles)
end
i1 = i2+1;
end
name = [var1 '.' var2];
if ~isconst(tmp)
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Mean.' name ' = p_mean;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Median.' name ' = p_median;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Variance.' name ' = p_var;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDinf.' name ' = hpd_interval(1);']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDsup.' name ' = hpd_interval(2);']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.deciles.' name ' = p_deciles;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = density;']);
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.density.(var1).(var2) = density;
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp,0,mh_conf_sig);
end
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Mean.(var1).(var2) = p_mean;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Median.(var1).(var2) = p_median;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Variance.(var1).(var2) = p_var;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.HPDinf.(var1).(var2) = hpd_interval(1);
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.HPDsup.(var1).(var2) = hpd_interval(2);
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.deciles.(var1).(var2) = p_deciles;
else
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Mean.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Median.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Variance.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDinf.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDsup.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.deciles.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = NaN;']);
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Mean.(var1).(var2) = NaN;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Median.(var1).(var2) = NaN;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Variance.(var1).(var2) = NaN;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.HPDinf.(var1).(var2) = NaN;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.HPDsup.(var1).(var2) = NaN;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.deciles.(var1).(var2) = NaN;
if options_.estimation.moments_posterior_density.indicator
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.density.(var1).(var2) = NaN;
end
end
if options_.contemporaneous_correlation
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp_corr_mat,1,mh_conf_sig);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Mean.' name ' = p_mean;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Median.' name ' = p_median;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Variance.' name ' = p_var;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.HPDinf.' name ' = hpd_interval(1);']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.HPDsup.' name ' = hpd_interval(2);']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.deciles.' name ' = p_deciles;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.density.' name ' = density;']);
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp_corr_mat,1,mh_conf_sig);
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.density.(var1).(var2) = density;
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp_corr_mat,0,mh_conf_sig);
end
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.Mean.(var1).(var2) = p_mean;
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.Median.(var1).(var2) = p_median;
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.Variance.(var1).(var2) = p_var;
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.HPDinf.(var1).(var2) = hpd_interval(1);
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.HPDsup.(var1).(var2) = hpd_interval(2);
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.deciles.(var1).(var2) = p_deciles;
end

View File

@ -487,6 +487,10 @@ options_.selected_variables_only = 0;
options_.contemporaneous_correlation = 0;
options_.initialize_estimated_parameters_with_the_prior_mode = 0;
options_.estimation_dll = 0;
options_.estimation.moments_posterior_density.indicator = 1;
options_.estimation.moments_posterior_density.gridpoints = 2^9;
options_.estimation.moments_posterior_density.bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
options_.estimation.moments_posterior_density.kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
% Misc
options_.conf_sig = 0.6;
oo_.exo_simul = [];

View File

@ -79,6 +79,9 @@ Median = zeros(n2,nvar);
Var = zeros(n2,nvar);
Distrib = zeros(9,n2,nvar);
HPD = zeros(2,n2,nvar);
if options_.estimation.moments_posterior_density.indicator
Density = zeros(options_.estimation.moments_posterior_density.gridpoints,2,n2,nvar);
end
fprintf(['Estimation::mcmc: ' tit1 '\n']);
stock1 = zeros(n1,n2,B);
k = 0;
@ -112,18 +115,33 @@ if filter_step_ahead_indicator
Var_filter_step_ahead = zeros(filter_steps,nvar,n2);
Distrib_filter_step_ahead = zeros(9,filter_steps,nvar,n2);
HPD_filter_step_ahead = zeros(2,filter_steps,nvar,n2);
if options_.estimation.moments_posterior_density.indicator
Density_filter_step_ahead = zeros(options_.estimation.moments_posterior_density.gridpoints,2,filter_steps,nvar,n2);
end
end
tmp =zeros(B,1);
for i = 1:nvar
for j = 1:n2
[Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig);
if options_.estimation.moments_posterior_density.indicator
[Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i),Density(:,:,j,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
else
[Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig);
end
if filter_step_ahead_indicator
for K_step = 1:length(options_.filter_step_ahead)
[Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ...
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig);
end
if options_.estimation.moments_posterior_density.indicator
for K_step = 1:length(options_.filter_step_ahead)
[Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j),Density_filter_step_ahead(:,:,K_step,i,j) ] = ...
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
end
else
for K_step = 1:length(options_.filter_step_ahead)
[Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ...
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig);
end
end
end
end
end
@ -142,21 +160,27 @@ end
for i = 1:nvar
name = deblank(names1(SelecVariables(i),:));
eval(['oo_.' name3 '.Mean.' name ' = Mean(:,i);']);
eval(['oo_.' name3 '.Median.' name ' = Median(:,i);']);
eval(['oo_.' name3 '.Var.' name ' = Var(:,i);']);
eval(['oo_.' name3 '.deciles.' name ' = Distrib(:,:,i);']);
eval(['oo_.' name3 '.HPDinf.' name ' = HPD(1,:,i)'';']);
eval(['oo_.' name3 '.HPDsup.' name ' = HPD(2,:,i)'';']);
oo_.(name3).Mean.(name) = Mean(:,i);
oo_.(name3).Median.(name) = Median(:,i);
oo_.(name3).Var.(name) = Var(:,i);
oo_.(name3).deciles.(name) = Distrib(:,:,i);
oo_.(name3).HPDinf.(name) = HPD(1,:,i)';
oo_.(name3).HPDsup.(name) = HPD(2,:,i)';
if options_.estimation.moments_posterior_density.indicator
oo_.(name3).density.(name) = Density(:,:,:,i);
end
if filter_step_ahead_indicator
for K_step = 1:length(options_.filter_step_ahead)
name4=['Filtered_Variables_',num2str(K_step),'_step_ahead'];
eval(['oo_.' name4 '.Mean.' name ' = squeeze(Mean_filter_step_ahead(K_step,i,:));']);
eval(['oo_.' name4 '.Median.' name ' = squeeze(Median_filter_step_ahead(K_step,i,:));']);
eval(['oo_.' name4 '.Var.' name ' = squeeze(Var_filter_step_ahead(K_step,i,:));']);
eval(['oo_.' name4 '.deciles.' name ' = squeeze(Distrib_filter_step_ahead(:,K_step,i,:));']);
eval(['oo_.' name4 '.HPDinf.' name ' = squeeze(HPD_filter_step_ahead(1,K_step,i,:));']);
eval(['oo_.' name4 '.HPDsup.' name ' = squeeze(HPD_filter_step_ahead(2,K_step,i,:));']);
oo_.(name4).Mean.(name) = squeeze(Mean_filter_step_ahead(K_step,i,:));
oo_.(name4).Median.(name) = squeeze(Median_filter_step_ahead(K_step,i,:));
oo_.(name4).Var.(name) = squeeze(Var_filter_step_ahead(K_step,i,:));
oo_.(name4).deciles.(name) = squeeze(Distrib_filter_step_ahead(:,K_step,i,:));
oo_.(name4).HPDinf.(name) = squeeze(HPD_filter_step_ahead(1,K_step,i,:));
oo_.(name4).HPDsup.(name) = squeeze(HPD_filter_step_ahead(2,K_step,i,:));
if options_.estimation.moments_posterior_density.indicator
oo_.(name4).density.(name) = squeeze(Density_filter_step_ahead(:,:,K_step,i,:));
end
end
end
end
@ -267,10 +291,4 @@ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fclose(fidTeX);
end
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']);
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']);

View File

@ -62,7 +62,7 @@ switch type
dsge_simulated_theoretical_variance_decomposition(SampleSize,M_,options_,oo_,'posterior');
end
oo_ = variance_decomposition_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_);
M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
case 'correlation'
if nargin==narg1
[nvar,vartan,NumberOfFiles] = ...
@ -76,7 +76,7 @@ switch type
dsge_simulated_theoretical_conditional_variance_decomposition(SampleSize,arg3,M_,options_,oo_,'posterior');
end
oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_);
arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
otherwise
disp('Not yet implemented')
end

View File

@ -1,4 +1,4 @@
function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info,mh_conf_sig)
function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info,mh_conf_sig,kernel_options)
% Computes posterior mean, median, variance, HPD interval, deciles, and density from posterior draws.
%
% INPUTS
@ -38,6 +38,16 @@ function [post_mean, post_median, post_var, hpd_interval, post_deciles, density]
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<4
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
else
number_of_grid_points = kernel_options.gridpoints;
bandwidth = kernel_options.bandwidth;
kernel_function = kernel_options.kernel_function;
end
xx = xx(:);
xx = sort(xx);
@ -76,10 +86,7 @@ else
end
density = [];
if info
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
if post_var > 1e-12
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
optimal_bandwidth = mh_optimal_bandwidth(xx,number_of_draws,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,...
number_of_draws,optimal_bandwidth,kernel_function);

View File

@ -63,7 +63,7 @@ switch type
dsge_simulated_theoretical_variance_decomposition(SampleSize,M_,options_,oo_,'prior');
end
oo_ = variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_);
M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
case 'correlation'
if nargin==narg1
[nvar,vartan,NumberOfFiles] = ...
@ -77,7 +77,7 @@ switch type
dsge_simulated_theoretical_conditional_variance_decomposition(SampleSize,arg3,M_,options_,oo_,'prior');
end
oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_);
arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
otherwise
disp('Not yet implemented')
end

View File

@ -1,14 +1,14 @@
function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_)
function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_,options_)
% function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_)
% This function analyses the (posterior or prior) distribution of the
% endogenous variables' variance decomposition.
%
%
% INPUTS
% NumberOfSimulations [integer] scalar, number of simulations.
% type [string] 'prior' or 'posterior'
% dname [string] directory name where to save
% fname [string] name of the mod-file
% exonames [string] (n_exo*char_length) character array with names of exogenous variables
% exonames [string] (n_exo*char_length) character array with names of exogenous variables
% exo [string] name of current exogenous
% variable
% vartan [string] (n_endo*char_length) character array with name
@ -18,6 +18,7 @@ function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname
% mh_conf_sig [double] 2 by 1 vector with upper
% and lower bound of HPD intervals
% oo_ [structure] Dynare structure where the results are saved.
% options_ [structure] Dynare options structure
%
% OUTPUTS
% oo_ [structure] Dynare structure where the results are saved.
@ -62,11 +63,11 @@ end
name = [ var '.' exo ];
if isfield(oo_, [ TYPE 'TheoreticalMoments'])
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
if isfield(temporary_structure,'dsge')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge;
if isfield(temporary_structure,'VarianceDecomposition')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.Mean;'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.Mean;
if isfield(temporary_structure,name)
% Nothing to do.
return
@ -93,7 +94,7 @@ if t3<1.0e-12
end
if abs(t1-1)<1.0e-12
t1 = 1;
end
end
p_mean = t1;
p_median = t1;
p_var = 0;
@ -101,13 +102,20 @@ if t3<1.0e-12
p_deciles = NaN(9,1);
density = NaN;
else
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp,0,mh_conf_sig);
end
end
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.Mean.' name ' = p_mean;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.Median.' name ' = p_median;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.Variance.' name ' = p_var;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.HPDinf.' name ' = hpd_interval(1);']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.HPDsup.' name ' = hpd_interval(2);']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.deciles.' name ' = p_deciles;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.density.' name ' = density;']);
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.Mean.(var).(exo) = p_mean;
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.Median.(var).(exo) = p_median;
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.Variance.(var).(exo) = p_var;
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.HPDinf.(var).(exo) = hpd_interval(1);
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.HPDsup.(var).(exo) = hpd_interval(2);
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.deciles.(var).(exo) = p_deciles;
if options_.estimation.moments_posterior_density.indicator
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.density.(var).(exo) = density;
end

View File

@ -114,7 +114,7 @@ class ParsingDriver;
%token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS
%token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE
%token PARALLEL_LOCAL_FILES PARAMETERS PARAMETER_SET PARTIAL_INFORMATION PERFECT_FORESIGHT PERIODS PERIOD PLANNER_OBJECTIVE PLOT_CONDITIONAL_FORECAST PLOT_PRIORS PREFILTER PRESAMPLE
%token PERFECT_FORESIGHT_SETUP PERFECT_FORESIGHT_SOLVER POSTERIOR_KERNEL_DENSITY
%token PERFECT_FORESIGHT_SETUP PERFECT_FORESIGHT_SOLVER NO_POSTERIOR_KERNEL_DENSITY
%token PRINT PRIOR_MC PRIOR_TRUNC PRIOR_MODE PRIOR_MEAN POSTERIOR_MODE POSTERIOR_MEAN POSTERIOR_MEDIAN PRUNING
%token <string_val> QUOTED_STRING
%token QZ_CRITERIUM QZ_ZERO_THRESHOLD FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT TRUNCATE
@ -1739,7 +1739,7 @@ estimation_options : o_datafile
| o_silent_optimizer
| o_proposal_distribution
| o_student_degrees_of_freedom
| o_posterior_kernel_density
| o_no_posterior_kernel_density
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@ -2692,8 +2692,8 @@ o_mh_jscale : MH_JSCALE EQUAL non_negative_number { driver.option_num("mh_jscale
o_optim : OPTIM EQUAL '(' optim_options ')';
o_tarb_optim : TARB_OPTIM EQUAL '(' tarb_optim_options ')';
o_proposal_distribution : PROPOSAL_DISTRIBUTION EQUAL symbol { driver.option_str("proposal_distribution", $3); };
o_posterior_kernel_density : POSTERIOR_KERNEL_DENSITY
{ driver.option_num("posterior_kernel_density.indicator", "1"); }
o_no_posterior_kernel_density : NO_POSTERIOR_KERNEL_DENSITY
{ driver.option_num("moments_posterior_density.indicator", "0"); }
;
o_student_degrees_of_freedom : STUDENT_DEGREES_OF_FREEDOM EQUAL INT_NUMBER { driver.option_num("student_degrees_of_freedom", $3); };
o_mh_init_scale : MH_INIT_SCALE EQUAL non_negative_number { driver.option_num("mh_init_scale", $3); };

View File

@ -377,7 +377,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>montecarlo {return token::MONTECARLO;}
<DYNARE_STATEMENT>distribution_approximation {return token::DISTRIBUTION_APPROXIMATION;}
<DYNARE_STATEMENT>proposal_distribution {return token::PROPOSAL_DISTRIBUTION;}
<DYNARE_STATEMENT>posterior_kernel_density {return token::POSTERIOR_KERNEL_DENSITY;}
<DYNARE_STATEMENT>no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;}
<DYNARE_STATEMENT>student_degrees_of_freedom {return token::STUDENT_DEGREES_OF_FREEDOM;}
<DYNARE_STATEMENT>alpha {

View File

@ -82,4 +82,4 @@ varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=3000,mh_nblocks=2,mh_jscale=0.8,moments_varendo,selected_variables_only) y m;
estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=3000,mh_nblocks=2,mh_jscale=0.8,moments_varendo,selected_variables_only,contemporaneous_correlation,smoother,forecast=8) y m;