From 9651f62aa0feefcad7ad8e7fee0e4664798063f4 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 6 Aug 2015 08:43:44 +0200 Subject: [PATCH 1/5] Fix bug preventing selecting modified harmonic mean for model_comparison --- matlab/model_comparison.m | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/matlab/model_comparison.m b/matlab/model_comparison.m index 43cf0cc8f..800d544ea 100644 --- a/matlab/model_comparison.m +++ b/matlab/model_comparison.m @@ -1,18 +1,26 @@ -function PosteriorOddsTable = model_comparison(ModelNames,ModelPriors,oo,options_,fname) -% Bayesian model comparison. This function computes Odds ratios and -% estimate a posterior density over a colletion of models. +function oo = model_comparison(ModelNames,ModelPriors,oo,options_,fname) +% function oo = model_comparison(ModelNames,ModelPriors,oo,options_,fname) +% Conducts Bayesian model comparison. This function computes Odds ratios and +% estimates a posterior density over a collection of models. % % INPUTS % ModelNames [string] m*1 cell array of string. % ModelPriors [double] m*1 vector of prior probabilities +% oo [struct] Dynare results structure +% options_ [struct] Dynare options structure +% fname [string] name of the current mod-file % % OUTPUTS -% none +% oo [struct] Dynare results structure containing the +% results in a field PosteriorOddsTable +% +% ALGORITHM +% See e.g. Koop (2003): Bayesian Econometrics % % SPECIAL REQUIREMENTS % none -% Copyright (C) 2007-2011 Dynare Team +% Copyright (C) 2007-2015 Dynare Team % % This file is part of Dynare. % @@ -47,13 +55,18 @@ end % The marginal densities are based on Laplace approxiations (default) or % modified harmonic mean estimators. -if isfield(options_,'model_comparison_approximation') - type = options_.model_comparison_approximation; - if strcmp(type,'Laplace') +if isfield(options_,'mc_marginal_density') + type = options_.mc_marginal_density; + if strcmp(type,'laplace') || strcmp(type,'Laplace') type = 'LaplaceApproximation'; + title = 'Model Comparison (based on Laplace approximation)'; + elseif strcmp(type,'modifiedharmonicmean') || strcmp(type,'ModifiedHarmonicMean') + type = 'ModifiedHarmonicMean'; + title = 'Model Comparison (based on Modified Harmonic Mean Estimator)'; end else type = 'LaplaceApproximation'; + title = 'Model Comparison (based on Laplace approximation)'; end % Get the estimated logged marginal densities. @@ -100,7 +113,6 @@ elmpd = exp(lmpd-maxval); % Now I display the posterior probabilities. -title = 'Model Comparison'; headers = char('Model',ShortModelNames{:}); if prior_flag labels = char('Priors','Log Marginal Density','Bayes Ratio', ... From 680625f260e4b9c70b5357a4d036811c59fd0c3b Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 6 Aug 2015 08:44:47 +0200 Subject: [PATCH 2/5] Save results of model_comparison to oo_ --- matlab/model_comparison.m | 10 +++++++++- preprocessor/ComputingTasks.cc | 2 +- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/matlab/model_comparison.m b/matlab/model_comparison.m index 800d544ea..b8847f51b 100644 --- a/matlab/model_comparison.m +++ b/matlab/model_comparison.m @@ -111,21 +111,29 @@ lmpd = log(ModelPriors)+MarginalLogDensity; [maxval,k] = max(lmpd); elmpd = exp(lmpd-maxval); - % Now I display the posterior probabilities. headers = char('Model',ShortModelNames{:}); if prior_flag labels = char('Priors','Log Marginal Density','Bayes Ratio', ... 'Posterior Model Probability'); + field_labels={'Priors','Log_Marginal_Density','Bayes_Ratio', ... + 'Posterior_Model_Probability'}; values = [ModelPriors';MarginalLogDensity';exp(lmpd-lmpd(1))'; ... elmpd'/sum(elmpd)]; else labels = char('Priors','Log Marginal Density','Bayes Ratio','Posterior Odds Ratio', ... 'Posterior Model Probability'); + field_labels={'Priors','Log_Marginal_Density','Bayes_Ratio','Posterior_Odds_Ratio','Posterior_Model_Probability'}; values = [ModelPriors';MarginalLogDensity'; exp(MarginalLogDensity-MarginalLogDensity(1))'; ... exp(lmpd-lmpd(1))'; elmpd'/sum(elmpd)]; end +for model_iter=1:NumberOfModels + for var_iter=1:size(labels,1) + oo.Model_Comparison.(deblank(headers(1+model_iter,:))).(field_labels{var_iter})=values(var_iter,model_iter); + end +end + dyntable(title,headers,labels,values, 0, 15, 6); diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 28de2ee55..37a5cea1a 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -1171,7 +1171,7 @@ ModelComparisonStatement::writeOutput(ostream &output, const string &basename, b output << "ModelNames_ = { ModelNames_{:} '" << (*it).first << "'};" << endl; output << "ModelPriors_ = [ ModelPriors_ ; " << (*it).second << "];" << endl; } - output << "model_comparison(ModelNames_,ModelPriors_,oo_,options_,M_.fname);" << endl; + output << "oo_ = model_comparison(ModelNames_,ModelPriors_,oo_,options_,M_.fname);" << endl; } PlannerObjectiveStatement::PlannerObjectiveStatement(StaticModel *model_tree_arg) : From 35bd2ed5d07c54b68a3db302497189a4a1ca98fd Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 6 Aug 2015 08:59:48 +0200 Subject: [PATCH 3/5] Add unit test for model_comparison --- tests/Makefile.am | 4 + tests/estimation/fs2000_model_comparison.mod | 110 +++++++++++++++++++ 2 files changed, 114 insertions(+) create mode 100644 tests/estimation/fs2000_model_comparison.mod diff --git a/tests/Makefile.am b/tests/Makefile.am index a639bf838..b2244a596 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -3,6 +3,7 @@ MODFILES = \ estimation/fs2000_MCMC_jumping_covariance.mod \ estimation/fs2000_initialize_from_calib.mod \ estimation/fs2000_calibrated_covariance.mod \ + estimation/fs2000_model_comparison.mod \ estimation/MH_recover/fs2000_recover.mod \ estimation/t_proposal/fs2000_student.mod \ estimation/TaRB/fs2000_tarb.mod \ @@ -294,6 +295,9 @@ AIM/fs2000x10_L9_L_AIM.o.trs: AIM/fs2000x10_L9_L.o.trs AIM/ls2003_2L0L_AIM.o.trs: AIM/ls2003_2L0L.o.trs AIM/ls2003_2L2L_AIM.o.trs: AIM/ls2003_2L2L.o.trs +estimation/fs2000_model_comparison.m.trs: estimation/fs2000.m.trs estimation/fs2000_initialize_from_calib.m.trs estimation/fs2000_calibrated_covariance.m.trs +estimation/fs2000_model_comparison.o.trs: estimation/fs2000.o.trs estimation/fs2000_initialize_from_calib.o.trs estimation/fs2000_calibrated_covariance.o.trs + optimizers/fs2000_102.m.trs: estimation/fs2000.m.trs optimizers/fs2000_102.o.trs: estimation/fs2000.o.trs diff --git a/tests/estimation/fs2000_model_comparison.mod b/tests/estimation/fs2000_model_comparison.mod new file mode 100644 index 000000000..b350eb337 --- /dev/null +++ b/tests/estimation/fs2000_model_comparison.mod @@ -0,0 +1,110 @@ +// See fs2000.mod in the examples/ directory for details on the model + +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +steady_state_model; + dA = exp(gam); + gst = 1/dA; + m = mst; + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; +end; + + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +check; + +estimated_params; +% alp, beta_pdf, 0.356, 0.02; +bet, beta_pdf, 0.993, 0.002; +% gam, normal_pdf, 0.0085, 0.003; +% mst, normal_pdf, 1.0002, 0.007; +% rho, beta_pdf, 0.129, 0.223; +% psi, beta_pdf, 0.65, 0.05; +% del, beta_pdf, 0.01, 0.005; +% stderr e_a, inv_gamma_pdf, 0.035449, inf; +% stderr e_m, inv_gamma_pdf, 0.008862, inf; +end; + +varobs gp_obs gy_obs; + +options_.solve_tolf = 1e-12; + +estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=2000,mh_nblocks=1,mh_jscale=0.8); + + +model_comparison fs2000(0.5) fs2000_calibrated_covariance(0.5); +if isfield(oo_,'Model_Comparison') + error('Model Comparison should not work with ML') +end +model_comparison fs2000 fs2000_initialize_from_calib; +model_comparison fs2000(0.5) fs2000_initialize_from_calib(0.5); +model_comparison fs2000(0.25) fs2000_initialize_from_calib(0.25); +model_comparison (marginal_density=laplace) fs2000(0.5) fs2000_initialize_from_calib(0.5); +if oo_.Model_Comparison.fs2000.Posterior_Model_Probability < oo_.Model_Comparison.fs2000_initialize_from_calib.Posterior_Model_Probability && ... + oo_.Model_Comparison.fs2000.Log_Marginal_Density > oo_.Model_Comparison.fs2000_initialize_from_calib.Log_Marginal_Density + error('Model Comparison not correct') +end +oo_laplace=oo_; +model_comparison (marginal_density=modifiedharmonicmean) fs2000(0.5) fs2000_initialize_from_calib(0.75); +if abs(oo_laplace.Model_Comparison.fs2000.Log_Marginal_Density-oo_.Model_Comparison.fs2000.Log_Marginal_Density)>0.1 + error('Laplace and Harmonic Mean do not match') +end +model_comparison (marginal_density=modifiedharmonicmean) fs2000(0) fs2000_initialize_from_calib(1); +if oo_.Model_Comparison.fs2000.Posterior_Model_Probability~=0 && oo_.Model_Comparison.fs2000_initialize_from_calib.Log_Marginal_Density~=1 + error('Incorporation of prior wrong') +end + +model_comparison(marginal_density=laplace) fs2000(0.333) fs2000_initialize_from_calib(0.333) fs2000_model_comparison(0.333); From a072294cab3901beac15a582637db0e22cafc9fa Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 6 Aug 2015 10:36:45 +0200 Subject: [PATCH 4/5] Improve documentation for model_comparison.m --- doc/dynare.texi | 46 +++++++++++++++++++++++++++++++++++++-- matlab/model_comparison.m | 4 ++-- 2 files changed, 46 insertions(+), 4 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index f55599867..c0079db0d 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5814,12 +5814,14 @@ Standard deviation of structural shocks @defvr {MATLAB/Octave variable} oo_.MarginalDensity.LaplaceApproximation -Variable set by the @code{estimation} command. +Variable set by the @code{estimation} command. Stores the marginal data density +based on the Laplace Approximation. @end defvr @defvr {MATLAB/Octave variable} oo_.MarginalDensity.ModifiedHarmonicMean Variable set by the @code{estimation} command, if it is used with -@code{mh_replic > 0} or @code{load_mh_file} option. +@code{mh_replic > 0} or @code{load_mh_file} option. Stores the marginal data density +based on @cite{Geweke (1999)} Modified Harmonic Mean estimator. @end defvr @defvr {MATLAB/Octave variable} oo_.FilteredVariables @@ -6154,6 +6156,20 @@ Dynare: it is the user's responsibility to make sure that model comparison is based on proper priors. Note that, for obvious reasons, this is not an issue if the compared marginal densities are based on Laplace approximations. +@optionshead + +@table @code + +@item marginal_density = @var{ESTIMATOR} +Specifies the estimator for computing the marginal data density. @var{ESTIMATOR} can +take one of the following two values: +@code{laplace} for the Laplace estimator or @code{modifiedharmonicmean} for the +@cite{Geweke (1999)} Modified Harmonic Mean estimator. Default value: @code{laplace} +@end table + +@outputhead + +The results are stored in @code{oo_.Model_Comparison}, which is described below. @examplehead @@ -6165,6 +6181,32 @@ over @code{alt_model}. @end deffn +@defvr {MATLAB/Octave variable} oo_.Model_Comparison +Variable set by the @code{model_comparison} command. Fields are of the form: +@example +@code{oo_.Model_Comparison.@var{FILENAME}.@var{VARIABLE_NAME}} +@end example +where @var{FILENAME} is the file name of the model and @var{VARIABLE_NAME} is one of the following: + +@table @code + +@item Prior +(Normalized) prior density over the model + +@item Log_Marginal_Density +Logarithm of the marginal data density + +@item Bayes_Ratio +Ratio of the marginal data density of the model relative to the one of the first declared model + +@item Posterior_Model_Probability +Posterior probability of the respective model + +@end table + +@end defvr + + @deffn Command shock_decomposition [@var{VARIABLE_NAME}]@dots{}; @deffnx Command shock_decomposition (@var{OPTIONS}@dots{}) [@var{VARIABLE_NAME}]@dots{}; diff --git a/matlab/model_comparison.m b/matlab/model_comparison.m index b8847f51b..2f5157a71 100644 --- a/matlab/model_comparison.m +++ b/matlab/model_comparison.m @@ -116,14 +116,14 @@ headers = char('Model',ShortModelNames{:}); if prior_flag labels = char('Priors','Log Marginal Density','Bayes Ratio', ... 'Posterior Model Probability'); - field_labels={'Priors','Log_Marginal_Density','Bayes_Ratio', ... + field_labels={'Prior','Log_Marginal_Density','Bayes_Ratio', ... 'Posterior_Model_Probability'}; values = [ModelPriors';MarginalLogDensity';exp(lmpd-lmpd(1))'; ... elmpd'/sum(elmpd)]; else labels = char('Priors','Log Marginal Density','Bayes Ratio','Posterior Odds Ratio', ... 'Posterior Model Probability'); - field_labels={'Priors','Log_Marginal_Density','Bayes_Ratio','Posterior_Odds_Ratio','Posterior_Model_Probability'}; + field_labels={'Prior','Log_Marginal_Density','Bayes_Ratio','Posterior_Odds_Ratio','Posterior_Model_Probability'}; values = [ModelPriors';MarginalLogDensity'; exp(MarginalLogDensity-MarginalLogDensity(1))'; ... exp(lmpd-lmpd(1))'; elmpd'/sum(elmpd)]; end From 6e93c83b69f202bab4436cc53e9af4bd7245d42a Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 6 Aug 2015 11:01:39 +0200 Subject: [PATCH 5/5] model_comparison: Do not issue warning if all prior densities are 1 If the user specifies nothing, this is the default, which should not result in a warning --- matlab/model_comparison.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/matlab/model_comparison.m b/matlab/model_comparison.m index 2f5157a71..4e4c904c3 100644 --- a/matlab/model_comparison.m +++ b/matlab/model_comparison.m @@ -47,8 +47,10 @@ else % The prior density has to sum up to one. prior_flag = 1; improper = abs(sum(ModelPriors)-1)>1e-6; if improper - disp('model_comparison:: The user supplied prior distribution over models is improper...') - disp('model_comparison:: The distribution is automatically rescaled!') + if ~all(ModelPriors==1) + disp('model_comparison:: The user supplied prior distribution over models is improper...') + disp('model_comparison:: The distribution is automatically rescaled!') + end ModelPriors=ModelPriors/sum(ModelPriors); end end