stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fill in posterior_mode with info from posterior samples.

time-shift
Marco Ratto 2015-09-25 18:23:21 +02:00 committed by Johannes Pfeifer
parent b830805cb6
commit 8bd963de64
3 changed files with 90 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/GetPosteriorParametersStatistics.m
View File

@ -72,7 +72,7 @@ skipline()
try
disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
catch
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
end
num_draws=NumberOfDraws*options_.mh_nblck;

2
matlab/dynare_estimation_1.m
View File

@ -443,7 +443,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
%% Estimation of the marginal density from the Mh draws:
if options_.mh_replic
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
% Store posterior statistics by parameter name
oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames);
if ~options_.nograph

90
matlab/marginal_density.m
View File

@ -1,4 +1,4 @@
function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_)
function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
% function marginal = marginal_density()
% Computes the marginal density
%
@ -58,6 +58,9 @@ lpost_mode = posterior_kernel_at_the_mode;
xparam1 = posterior_mean;
hh = inv(SIGMA);
fprintf(' Done!\n');
if ~isfield(oo_,'posterior_mode')
oo_=fill_mh_mode(posterior_mode',NaN(npar,1),M_,options_,estim_params_,bayestopt_,oo_,'posterior');
end
% save the posterior mean and the inverse of the covariance matrix
% (usefull if the user wants to perform some computations using
@ -120,4 +123,87 @@ while check_coverage
end
end
oo_.MarginalDensity.ModifiedHarmonicMean = mean(marginal(:,2));
oo_.MarginalDensity.ModifiedHarmonicMean = mean(marginal(:,2));
return
function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
%function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
%
% INPUTS
% o xparam1 [double] (p*1) vector of estimate parameters.
% o stdh [double] (p*1) vector of estimate parameters.
% o M_ Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o estim_params_ Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
% o options_ Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o bayestopt_ Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
% o oo_ Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%
% OUTPUTS
% o oo_ Matlab's structure gathering the results
%
% SPECIAL REQUIREMENTS
% None.
nvx = estim_params_.nvx; % Variance of the structural innovations (number of parameters).
nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters).
ncx = estim_params_.ncx; % Covariance of the structural innovations (number of parameters).
ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of parameters).
np = estim_params_.np ; % Number of deep parameters.
nx = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
if np
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
name = bayestopt_.name{ip};
eval(['oo_.' field_name '_mode.parameters.' name ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.parameters.' name ' = stdh(ip);']);
ip = ip+1;
end
end
if nvx
ip = 1;
for i=1:nvx
k = estim_params_.var_exo(i,1);
name = deblank(M_.exo_names(k,:));
eval(['oo_.' field_name '_mode.shocks_std.' name ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.shocks_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
if nvn
ip = nvx+1;
for i=1:nvn
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
eval(['oo_.' field_name '_mode.measurement_errors_std.' name ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.measurement_errors_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
if ncx
ip = nvx+nvn+1;
for i=1:ncx
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
eval(['oo_.' field_name '_mode.shocks_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.shocks_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
end
end
if ncn
ip = nvx+nvn+ncx+1;
for i=1:ncn
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
eval(['oo_.' field_name '_mode.measurement_errors_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.measurement_errors_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
end
end
return