Commit 6390830b authored by Johannes Pfeifer's avatar Johannes Pfeifer
Browse files

Store MCMC information recorded in record in oo_

Closes issue #315 (https://github.com/DynareTeam/dynare/issues/315)
parent f1406487
...@@ -4461,6 +4461,9 @@ Median of the posterior distribution ...@@ -4461,6 +4461,9 @@ Median of the posterior distribution
@item Std @item Std
Standard deviation of the posterior distribution Standard deviation of the posterior distribution
@item Variance
Variance of the posterior distribution
@item deciles @item deciles
Deciles of the distribution. Deciles of the distribution.
...@@ -4649,6 +4652,21 @@ oo_.posterior_mean.shocks_std.ex ...@@ -4649,6 +4652,21 @@ oo_.posterior_mean.shocks_std.ex
oo_.posterior_hpdsup.measurement_errors_corr.gdp_conso oo_.posterior_hpdsup.measurement_errors_corr.gdp_conso
@end example @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{}; @deffn Command model_comparison @var{FILENAME}[(@var{DOUBLE})]@dots{};
@deffnx Command model_comparison (marginal_density = laplace | modifiedharmonicmean) @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 ...@@ -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 With respect to the previous version of the toolbox, in order to work
properly, the GSA toolbox no longer requires that the Dynare 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}, Sensitivity analysis results are saved locally in @code{<mod_file>/GSA},
where @code{<mod_file>.mod} is the name of the DYNARE model file. where @code{<mod_file>.mod} is the name of the DYNARE model file.
......
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) %function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
% Random walk Metropolis-Hastings algorithm. % Random walk Metropolis-Hastings algorithm.
% %
...@@ -11,7 +11,7 @@ function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds ...@@ -11,7 +11,7 @@ function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds
% o varargin list of argument following mh_bounds % o varargin list of argument following mh_bounds
% %
% OUTPUTS % OUTPUTS
% None % o record [struct] structure describing the iterations
% %
% ALGORITHM % ALGORITHM
% Metropolis-Hastings. % Metropolis-Hastings.
......
...@@ -156,7 +156,7 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.m ...@@ -156,7 +156,7 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.m
oo_.Smoother.SteadyState = ys; oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff; oo_.Smoother.TrendCoeffs = trend_coeff;
if options_.filter_covariance if options_.filter_covariance
oo_.Smoother.variance = P; oo_.Smoother.Variance = P;
end end
i_endo = bayestopt_.smoother_saved_var_list; i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0 if options_.nk ~= 0
...@@ -482,12 +482,12 @@ if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation ...@@ -482,12 +482,12 @@ if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation
end end
oo_.posterior.optimization.mode = xparam1; oo_.posterior.optimization.mode = xparam1;
oo_.posterior.optimization.variance = []; oo_.posterior.optimization.Variance = [];
if ~options_.mh_posterior_mode_estimation if ~options_.mh_posterior_mode_estimation
if options_.cova_compute if options_.cova_compute
invhess = inv(hh); invhess = inv(hh);
stdh = sqrt(diag(invhess)); stdh = sqrt(diag(invhess));
oo_.posterior.optimization.variance = invhess; oo_.posterior.optimization.Variance = invhess;
end end
else else
variances = bayestopt_.p2.*bayestopt_.p2; variances = bayestopt_.p2.*bayestopt_.p2;
...@@ -921,7 +921,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... ...@@ -921,7 +921,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
ana_deriv = options_.analytic_derivation; ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 0; options_.analytic_derivation = 0;
if options_.cova_compute 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 else
error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.') error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.')
end end
...@@ -943,7 +943,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... ...@@ -943,7 +943,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if ~options_.nograph if ~options_.nograph
oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_); oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
end end
[oo_.posterior.metropolis.mean,oo_.posterior.metropolis.variance] ... [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
= GetPosteriorMeanVariance(M_,options_.mh_drop); = GetPosteriorMeanVariance(M_,options_.mh_drop);
else else
load([M_.fname '_results'],'oo_'); load([M_.fname '_results'],'oo_');
...@@ -974,7 +974,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha ...@@ -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); [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.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff; oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.variance = P; oo_.Smoother.Variance = P;
i_endo = bayestopt_.smoother_saved_var_list; i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0 if options_.nk ~= 0
oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ... oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
......
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. % Independent Metropolis-Hastings algorithm.
% %
...@@ -11,7 +11,7 @@ function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou ...@@ -11,7 +11,7 @@ function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou
% o varargin list of argument following mh_bounds % o varargin list of argument following mh_bounds
% %
% OUTPUTS % OUTPUTS
% None % o record [struct] structure describing the iterations
% %
% ALGORITHM % ALGORITHM
% Metropolis-Hastings. % Metropolis-Hastings.
...@@ -103,7 +103,7 @@ else ...@@ -103,7 +103,7 @@ else
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
for j=1:totCPU, for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1; 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.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.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))); record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
......
...@@ -229,7 +229,7 @@ for b = fblck:nblck, ...@@ -229,7 +229,7 @@ for b = fblck:nblck,
jsux = 0; jsux = 0;
if j == nruns(b) % I record the last draw... if j == nruns(b) % I record the last draw...
record.LastParameters(b,:) = x2(end,:); record.LastParameters(b,:) = x2(end,:);
record.LastLogLiK(b) = logpo2(end); record.LastLogPost(b) = logpo2(end);
end end
% size of next file in chain b % size of next file in chain b
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns); InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
......
...@@ -194,7 +194,7 @@ if ~options_.load_mh_file && ~options_.mh_recover ...@@ -194,7 +194,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
record.InitialParameters = ix2; record.InitialParameters = ix2;
record.InitialLogLiK = ilogpo2; record.InitialLogLiK = ilogpo2;
record.LastParameters = zeros(nblck,npar); record.LastParameters = zeros(nblck,npar);
record.LastLogLiK = zeros(nblck,1); record.LastLogPost = zeros(nblck,1);
record.LastFileNumber = AnticipatedNumberOfFiles ; record.LastFileNumber = AnticipatedNumberOfFiles ;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile; record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record'); save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
...@@ -253,7 +253,7 @@ elseif options_.load_mh_file && ~options_.mh_recover ...@@ -253,7 +253,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
else else
NewFile = ones(nblck,1)*(LastFileNumber+1); NewFile = ones(nblck,1)*(LastFileNumber+1);
end end
ilogpo2 = record.LastLogLiK; ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters; ix2 = record.LastParameters;
fblck = 1; fblck = 1;
fline = ones(nblck,1)*(LastLineNumber+1); fline = ones(nblck,1)*(LastLineNumber+1);
...@@ -292,7 +292,7 @@ elseif options_.mh_recover ...@@ -292,7 +292,7 @@ elseif options_.mh_recover
end end
%% Default initialization: %% Default initialization:
if OldMh if OldMh
ilogpo2 = record.LastLogLiK; ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters; ix2 = record.LastParameters;
else else
ilogpo2 = record.InitialLogLiK; ilogpo2 = record.InitialLogLiK;
......
...@@ -148,7 +148,7 @@ else ...@@ -148,7 +148,7 @@ else
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
for j=1:totCPU, for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1; 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.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.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))); record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
......
...@@ -247,7 +247,7 @@ for b = fblck:nblck, ...@@ -247,7 +247,7 @@ for b = fblck:nblck,
jsux = 0; jsux = 0;
if j == nruns(b) % I record the last draw... if j == nruns(b) % I record the last draw...
record.LastParameters(b,:) = x2(end,:); record.LastParameters(b,:) = x2(end,:);
record.LastLogLiK(b) = logpo2(end); record.LastLogPost(b) = logpo2(end);
end end
% size of next file in chain b % size of next file in chain b
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns); InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment