Merge pull request #1015 from JohannesPfeifer/model_comparison

Bugfixes and improvements for model_comparison

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{};
 

matlab/model_comparison.m

@@ -1,18 +1,26 @@
-function PosteriorOddsTable = model_comparison(ModelNames,ModelPriors,oo,options_,fname)
+function oo = model_comparison(ModelNames,ModelPriors,oo,options_,fname)
-% Bayesian model comparison. This function computes Odds ratios and
+% function oo = model_comparison(ModelNames,ModelPriors,oo,options_,fname)
-% estimate a posterior density over a colletion of models.
+% 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...')
+        if ~all(ModelPriors==1)
-        disp('model_comparison:: The distribution is automatically rescaled!')
+            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')
+if isfield(options_,'mc_marginal_density')
-    type = options_.model_comparison_approximation;
+    type = options_.mc_marginal_density;
-    if strcmp(type,'Laplace')
+    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);
 
 

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) :

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
 

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);