From 0aee5cbf5aa8d0285d3f60d942dda219f86dda59 Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@jrc.ec.europa.eu>
Date: Fri, 25 Sep 2015 18:23:21 +0200
Subject: [PATCH] Fill in posterior_mode with info from posterior samples.
---
matlab/GetPosteriorParametersStatistics.m | 2 +-
matlab/dynare_estimation_1.m | 2 +-
matlab/marginal_density.m | 90 ++++++++++++++++++++++-
3 files changed, 90 insertions(+), 4 deletions(-)
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index a8ee9a58bb..820b6e53f0 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -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;
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 79087715a0..2f2e6d4688 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -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
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index f24c344f45..32880f34ae 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -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));
\ No newline at end of file
+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
\ No newline at end of file
--
GitLab