Fill in posterior_mode with info from posterior samples.

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;

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

matlab/marginal_density.m

-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
@@ -121,3 +124,86 @@ while check_coverage
 end
 
 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