model_comparison

 Bayesian model comparison. This function computes Odds ratios and
 estimate a posterior density over a colletion of models.

 INPUTS
    ModelNames       [string]     m*1 cell array of string.
    ModelPriors      [double]     m*1 vector of prior probabilities 

 OUTPUTS
    none

 SPECIAL REQUIREMENTS
    none

0001 function PosteriorOddsTable = model_comparison(ModelNames,ModelPriors,oo,options_,fname)
0002 % Bayesian model comparison. This function computes Odds ratios and
0003 % estimate a posterior density over a colletion of models.
0004 %
0005 % INPUTS
0006 %    ModelNames       [string]     m*1 cell array of string.
0007 %    ModelPriors      [double]     m*1 vector of prior probabilities
0008 %
0009 % OUTPUTS
0010 %    none
0011 %
0012 % SPECIAL REQUIREMENTS
0013 %    none
0014 
0015 % Copyright (C) 2007-2011 Dynare Team
0016 %
0017 % This file is part of Dynare.
0018 %
0019 % Dynare is free software: you can redistribute it and/or modify
0020 % it under the terms of the GNU General Public License as published by
0021 % the Free Software Foundation, either version 3 of the License, or
0022 % (at your option) any later version.
0023 %
0024 % Dynare is distributed in the hope that it will be useful,
0025 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0026 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0027 % GNU General Public License for more details.
0028 %
0029 % You should have received a copy of the GNU General Public License
0030 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0031 
0032 NumberOfModels = size(ModelNames,2);
0033 
0034 disp(' ')
0035 disp(' ')
0036 if isempty(ModelPriors)
0037     prior_flag = 0;% empty_prior=0
0038     ModelPriors = ones(NumberOfModels,1)/NumberOfModels;
0039 else % The prior density has to sum up to one.
0040     prior_flag = 1;
0041     improper = abs(sum(ModelPriors)-1)>1e-6;
0042     if improper
0043         disp('model_comparison:: The user supplied prior distribution over models is improper...')
0044         disp('model_comparison:: The distribution is automatically rescaled!')
0045         ModelPriors=ModelPriors/sum(ModelPriors);
0046     end
0047 end
0048 
0049 % The marginal densities are based on Laplace approxiations (default) or
0050 % modified harmonic mean estimators.
0051 if isfield(options_,'model_comparison_approximation')
0052     type = options_.model_comparison_approximation;
0053     if strcmp(type,'Laplace')
0054         type = 'LaplaceApproximation';
0055     end
0056 else
0057     type = 'LaplaceApproximation';
0058 end
0059 
0060 % Get the estimated logged marginal densities.
0061 MarginalLogDensity = zeros(NumberOfModels,1);
0062 ShortModelNames = get_short_names(ModelNames);
0063 iname = strmatch(fname,ShortModelNames,'exact');
0064 
0065 for i=1:NumberOfModels
0066     if i==iname
0067         mstruct.oo_ = oo;
0068     else
0069         if strcmpi(ModelNames{i}(end-3:end),'.mod') || strcmpi(ModelNames{i}(end-3:end),'.dyn')
0070             mstruct = load([ModelNames{i}(1:end-4) '_results.mat' ],'oo_');
0071         else
0072             mstruct = load([ModelNames{i} '_results.mat' ],'oo_');
0073         end
0074     end
0075     try
0076         eval(['MarginalLogDensity(i) = mstruct.oo_.MarginalDensity.' type ';']) 
0077     catch
0078         if strcmpi(type,'LaplaceApproximation')
0079             disp(['MODEL_COMPARISON: I cant''t find the Laplace approximation associated to model ' ModelNames{i}])
0080             return
0081         elseif strcmpi(type,'ModifiedHarmonicMean')
0082             disp(['MODEL_COMPARISON: I cant''t find the modified harmonic mean estimate associated to model ' ModelNames{i}])
0083             return
0084         end
0085     end
0086 end
0087 
0088 % In order to avoid overflow, we divide the numerator and the denominator
0089 % of the Posterior Odds Ratio by the largest Marginal Posterior Density
0090 lmpd = log(ModelPriors)+MarginalLogDensity;
0091 [maxval,k] = max(lmpd);
0092 elmpd = exp(lmpd-maxval);
0093 
0094 
0095 % Now I display the posterior probabilities.
0096 title = 'Model Comparison'; 
0097 headers = char('Model',ShortModelNames{:});
0098 if prior_flag
0099     labels = char('Priors','Log Marginal Density','Bayes Ratio', ...
0100                   'Posterior Model Probability');
0101     values = [ModelPriors';MarginalLogDensity';exp(lmpd-lmpd(1))'; ...
0102               elmpd'/sum(elmpd)];
0103 else
0104     labels = char('Priors','Log Marginal Density','Bayes Ratio','Posterior Odds Ratio', ...
0105                   'Posterior Model Probability');
0106     values = [ModelPriors';MarginalLogDensity'; exp(MarginalLogDensity-MarginalLogDensity(1))'; ...
0107               exp(lmpd-lmpd(1))'; elmpd'/sum(elmpd)];
0108 end
0109 
0110 dyntable(title,headers,labels,values, 0, 15, 6);
0111 
0112 
0113 function name = get_model_name_without_path(modelname)
0114 idx = strfind(modelname,'\');
0115 if isempty(idx)
0116     idx = strfind(modelname,'/');
0117 end
0118 if isempty(idx)
0119     name = modelname;
0120     return
0121 end
0122 name = modelname(idx(end)+1:end);
0123 
0124 
0125 function name = get_model_name_without_extension(modelname)
0126 idx = strfind(modelname,'.mod');
0127 if isempty(idx)
0128     idx = strfind(modelname,'.dyn')
0129 end
0130 if isempty(idx)
0131     name = modelname;
0132     return
0133 end
0134 name = modelname(1:end-4);
0135 
0136 
0137 function modellist = get_short_names(modelnames)
0138 n = length(modelnames);
0139 modellist = {};
0140 for i=1:n
0141     name = get_model_name_without_extension(modelnames{i});
0142     name = get_model_name_without_path(name);
0143     modellist = {modellist{:} name};
0144 end