From 6390830b4d924df9dba14b3620a77a10409a1880 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Sun, 17 Mar 2013 22:49:28 +0100
Subject: [PATCH] Store MCMC information recorded in record in oo_
Closes issue #315 (https://github.com/DynareTeam/dynare/issues/315)
---
doc/dynare.texi | 20 ++++++++++++++++++-
matlab/adaptive_metropolis_hastings.m | 4 ++--
matlab/dynare_estimation_1.m | 12 +++++------
matlab/independent_metropolis_hastings.m | 6 +++---
matlab/independent_metropolis_hastings_core.m | 2 +-
matlab/metropolis_hastings_initialization.m | 6 +++---
matlab/random_walk_metropolis_hastings.m | 2 +-
matlab/random_walk_metropolis_hastings_core.m | 2 +-
8 files changed, 36 insertions(+), 18 deletions(-)
diff --git a/doc/dynare.texi b/doc/dynare.texi
index 0f01b3fa1..0d30190ba 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -4461,6 +4461,9 @@ Median of the posterior distribution
@item Std
Standard deviation of the posterior distribution
+@item Variance
+Variance of the posterior distribution
+
@item deciles
Deciles of the distribution.
@@ -4649,6 +4652,21 @@ oo_.posterior_mean.shocks_std.ex
oo_.posterior_hpdsup.measurement_errors_corr.gdp_conso
@end example
+@defvr {MATLAB/Octave variable} oo_.MC_record.Seeds
+Variable set by the @code{estimation} command. Stores seeds used in MCMC chains
+@end defvr
+
+@defvr {MATLAB/Octave variable} oo_.MC_record.AcceptationRates
+Variable set by the @code{estimation} command. Stores acceptation rates of the MCMC chains
+@end defvr
+
+@defvr {MATLAB/Octave variable} oo_.MC_record.LastParameters
+Variable set by the @code{estimation} command. Stores parameter vector of final MCMC chain draw
+@end defvr
+
+@defvr {MATLAB/Octave variable} oo_.MC_record.LastLogPost
+Variable set by the @code{estimation} command. Stores log-posterior of final MCMC chain draw
+@end defvr
@deffn Command model_comparison @var{FILENAME}[(@var{DOUBLE})]@dots{};
@deffnx Command model_comparison (marginal_density = laplace | modifiedharmonicmean) @var{FILENAME}[(@var{DOUBLE})]@dots{};
@@ -5316,7 +5334,7 @@ The discussion of the methodologies and their application is described in
With respect to the previous version of the toolbox, in order to work
properly, the GSA toolbox no longer requires that the Dynare
-estimation environment is setup.
+estimation environment is set up.
Sensitivity analysis results are saved locally in @code{<mod_file>/GSA},
where @code{<mod_file>.mod} is the name of the DYNARE model file.
diff --git a/matlab/adaptive_metropolis_hastings.m b/matlab/adaptive_metropolis_hastings.m
index ea71612f0..2b0cb451a 100644
--- a/matlab/adaptive_metropolis_hastings.m
+++ b/matlab/adaptive_metropolis_hastings.m
@@ -1,4 +1,4 @@
-function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
+function record=adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
%function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
% Random walk Metropolis-Hastings algorithm.
%
@@ -11,7 +11,7 @@ function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds
% o varargin list of argument following mh_bounds
%
% OUTPUTS
-% None
+% o record [struct] structure describing the iterations
%
% ALGORITHM
% Metropolis-Hastings.
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 45e47397f..f4a6dcada 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -156,7 +156,7 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.m
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
if options_.filter_covariance
- oo_.Smoother.variance = P;
+ oo_.Smoother.Variance = P;
end
i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0
@@ -482,12 +482,12 @@ if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation
end
oo_.posterior.optimization.mode = xparam1;
-oo_.posterior.optimization.variance = [];
+oo_.posterior.optimization.Variance = [];
if ~options_.mh_posterior_mode_estimation
if options_.cova_compute
invhess = inv(hh);
stdh = sqrt(diag(invhess));
- oo_.posterior.optimization.variance = invhess;
+ oo_.posterior.optimization.Variance = invhess;
end
else
variances = bayestopt_.p2.*bayestopt_.p2;
@@ -921,7 +921,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 0;
if options_.cova_compute
- feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ oo_.MC_record=feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
else
error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.')
end
@@ -943,7 +943,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if ~options_.nograph
oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
end
- [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.variance] ...
+ [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
= GetPosteriorMeanVariance(M_,options_.mh_drop);
else
load([M_.fname '_results'],'oo_');
@@ -974,7 +974,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
- oo_.Smoother.variance = P;
+ oo_.Smoother.Variance = P;
i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0
oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
diff --git a/matlab/independent_metropolis_hastings.m b/matlab/independent_metropolis_hastings.m
index bf198383f..8e24ad178 100644
--- a/matlab/independent_metropolis_hastings.m
+++ b/matlab/independent_metropolis_hastings.m
@@ -1,4 +1,4 @@
-function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
+function record=independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
% Independent Metropolis-Hastings algorithm.
%
@@ -11,7 +11,7 @@ function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou
% o varargin list of argument following mh_bounds
%
% OUTPUTS
-% None
+% o record [struct] structure describing the iterations
%
% ALGORITHM
% Metropolis-Hastings.
@@ -103,7 +103,7 @@ else
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
- record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));
+ record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)));
record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:);
record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)));
record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
diff --git a/matlab/independent_metropolis_hastings_core.m b/matlab/independent_metropolis_hastings_core.m
index 04b2683aa..efb50965e 100644
--- a/matlab/independent_metropolis_hastings_core.m
+++ b/matlab/independent_metropolis_hastings_core.m
@@ -229,7 +229,7 @@ for b = fblck:nblck,
jsux = 0;
if j == nruns(b) % I record the last draw...
record.LastParameters(b,:) = x2(end,:);
- record.LastLogLiK(b) = logpo2(end);
+ record.LastLogPost(b) = logpo2(end);
end
% size of next file in chain b
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m
index 06dc8a4ca..e56ddcf5f 100644
--- a/matlab/metropolis_hastings_initialization.m
+++ b/matlab/metropolis_hastings_initialization.m
@@ -194,7 +194,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
record.InitialParameters = ix2;
record.InitialLogLiK = ilogpo2;
record.LastParameters = zeros(nblck,npar);
- record.LastLogLiK = zeros(nblck,1);
+ record.LastLogPost = zeros(nblck,1);
record.LastFileNumber = AnticipatedNumberOfFiles ;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
@@ -253,7 +253,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
else
NewFile = ones(nblck,1)*(LastFileNumber+1);
end
- ilogpo2 = record.LastLogLiK;
+ ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
fblck = 1;
fline = ones(nblck,1)*(LastLineNumber+1);
@@ -292,7 +292,7 @@ elseif options_.mh_recover
end
%% Default initialization:
if OldMh
- ilogpo2 = record.LastLogLiK;
+ ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
else
ilogpo2 = record.InitialLogLiK;
diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m
index 453e4d615..1c75e2102 100644
--- a/matlab/random_walk_metropolis_hastings.m
+++ b/matlab/random_walk_metropolis_hastings.m
@@ -148,7 +148,7 @@ else
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
- record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));
+ record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)));
record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:);
record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)));
record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m
index 69062f66c..a92c3c21c 100644
--- a/matlab/random_walk_metropolis_hastings_core.m
+++ b/matlab/random_walk_metropolis_hastings_core.m
@@ -247,7 +247,7 @@ for b = fblck:nblck,
jsux = 0;
if j == nruns(b) % I record the last draw...
record.LastParameters(b,:) = x2(end,:);
- record.LastLogLiK(b) = logpo2(end);
+ record.LastLogPost(b) = logpo2(end);
end
% size of next file in chain b
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
--
GitLab