diff --git a/doc/dynare.texi b/doc/dynare.texi index 79635bee5..69294112a 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5840,12 +5840,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 @@ -6180,6 +6182,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 @@ -6191,6 +6207,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 43cf0cc8f..4e4c904c3 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. % @@ -39,21 +47,28 @@ 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 % 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. @@ -98,22 +113,29 @@ lmpd = log(ModelPriors)+MarginalLogDensity; [maxval,k] = max(lmpd); 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', ... 'Posterior Model Probability'); + 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={'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 +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 72a26f350..46ec69cc2 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) : diff --git a/tests/Makefile.am b/tests/Makefile.am index 60c601faf..1516d2641 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 \ @@ -298,6 +299,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);