Home > matlab > GetPosteriorParametersStatistics.m

GetPosteriorParametersStatistics

PURPOSE ^

This function prints and saves posterior estimates after the mcmc

SYNOPSIS ^

function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_)

DESCRIPTION ^

 This function prints and saves posterior estimates after the mcmc
 (+updates of oo_ & TeX output). 
 
 INPUTS 
   estim_params_    [structure] 
   M_               [structure]
   options_         [structure]
   bayestopt_       [structure]
   oo_              [structure]
  
 OUTPUTS 
   oo_              [structure]  

 SPECIAL REQUIREMENTS
   None.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_)
0002 % This function prints and saves posterior estimates after the mcmc
0003 % (+updates of oo_ & TeX output).
0004 %
0005 % INPUTS
0006 %   estim_params_    [structure]
0007 %   M_               [structure]
0008 %   options_         [structure]
0009 %   bayestopt_       [structure]
0010 %   oo_              [structure]
0011 %
0012 % OUTPUTS
0013 %   oo_              [structure]
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %   None.
0017 
0018 % Copyright (C) 2006-2011 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 %if ~options_.mh_replic && options_.load_mh_file
0036 %   load([M_.fname '_results.mat'],'oo_');
0037 %end
0038 
0039 TeX     = options_.TeX;
0040 nblck   = options_.mh_nblck;
0041 nvx     = estim_params_.nvx;
0042 nvn     = estim_params_.nvn;
0043 ncx     = estim_params_.ncx;
0044 ncn     = estim_params_.ncn;
0045 np      = estim_params_.np ;
0046 nx      = nvx+nvn+ncx+ncn+np;
0047 
0048 DirectoryName = CheckPath('metropolis',M_.dname);
0049 OutputDirectoryName = CheckPath('Output',M_.dname);
0050 
0051 load([ DirectoryName '/'  M_.fname '_mh_history'])
0052 FirstMhFile = record.KeepedDraws.FirstMhFile;
0053 FirstLine = record.KeepedDraws.FirstLine;
0054 TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
0055 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0056 FirstMhFile = record.KeepedDraws.FirstMhFile;
0057 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
0058 clear record;
0059 
0060 pnames=['     ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
0061 header_width = row_header_width(M_,estim_params_,bayestopt_);
0062 tit2 = sprintf('%-*s %10s %10s %16s %6s %10s\n',header_width+2,' ','prior mean','post. mean','conf. interval','prior','pstdev');
0063 pformat = '%-*s %10.3f %10.4f %10.4f %8.4f %6s %10.4f';
0064 
0065 disp(' ');disp(' ');disp('ESTIMATION RESULTS');disp(' ');
0066 try
0067     disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
0068 catch
0069     [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
0070     disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
0071 end
0072 if np
0073     type = 'parameters';
0074     if TeX
0075         fid = TeXBegin(OutputDirectoryName,M_.fname,1,type);
0076     end
0077     disp(' ')
0078     disp(type)
0079     disp(tit2)
0080     ip = nvx+nvn+ncx+ncn+1;
0081     for i=1:np
0082         if options_.mh_replic
0083             Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0084             [post_mean, post_median, post_var, hpd_interval, post_deciles, ...
0085              density] = posterior_moments(Draws,1,options_.mh_conf_sig);
0086             name = bayestopt_.name{ip};
0087             oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
0088         else
0089             try
0090                 name = bayestopt_.name{ip};
0091                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
0092             catch
0093                 Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0094                 [post_mean, post_median, post_var, hpd_interval, post_deciles, ...
0095                  density] = posterior_moments(Draws,1,options_.mh_conf_sig);
0096                 name = bayestopt_.name{ip};
0097                 oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);                
0098             end
0099         end
0100         disp(sprintf(pformat,header_width,name,bayestopt_.p1(ip),...
0101                      post_mean, ...
0102                      hpd_interval, ...
0103                      pnames(bayestopt_.pshape(ip)+1,:), ...
0104                      bayestopt_.p2(ip)));    
0105         if TeX
0106             TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
0107                     bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
0108         end
0109         ip = ip+1;
0110     end
0111     if TeX
0112         TeXEnd(fid,1,type);
0113     end
0114 end
0115 if nvx
0116     type = 'shocks_std';
0117     if TeX
0118         fid = TeXBegin(OutputDirectoryName,M_.fname,2,'standard deviation of structural shocks');
0119     end
0120     ip = 1;
0121     disp(' ')
0122     disp('standard deviation of shocks')
0123     disp(tit2)
0124     for i=1:nvx
0125         if options_.mh_replic
0126             Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0127             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0128                 posterior_moments(Draws,1,options_.mh_conf_sig);
0129             k = estim_params_.var_exo(i,1);
0130             name = deblank(M_.exo_names(k,:));
0131             oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
0132             M_.Sigma_e(k,k) = post_mean*post_mean;
0133         else
0134             try
0135                 k = estim_params_.var_exo(i,1);
0136                 name = deblank(M_.exo_names(k,:));
0137                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
0138             catch
0139                 Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0140                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0141                     posterior_moments(Draws,1,options_.mh_conf_sig);
0142                 k = estim_params_.var_exo(i,1);
0143                 name = deblank(M_.exo_names(k,:));
0144                 oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
0145                 M_.Sigma_e(k,k) = post_mean*post_mean;
0146             end
0147         end
0148         disp(sprintf(pformat,header_width,name,bayestopt_.p1(ip),post_mean,hpd_interval,...
0149                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
0150         if TeX
0151             TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
0152                     bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
0153         end
0154         ip = ip+1;
0155     end
0156     if TeX
0157         TeXEnd(fid,2,'standard deviation of structural shocks');        
0158     end
0159 end
0160 if nvn
0161     type = 'measurement_errors_std';
0162     if TeX
0163         fid = TeXBegin(OutputDirectoryName,M_.fname,3,'standard deviation of measurement errors')
0164     end
0165     disp(' ')
0166     disp('standard deviation of measurement errors')
0167     disp(tit2)
0168     ip = nvx+1;
0169     for i=1:nvn
0170         if options_.mh_replic
0171             Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0172             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0173                 posterior_moments(Draws,1,options_.mh_conf_sig);
0174             name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
0175             oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
0176         else
0177             try
0178                 name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
0179                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
0180             catch
0181                 Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0182                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0183                     posterior_moments(Draws,1,options_.mh_conf_sig);
0184                 name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
0185                 oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
0186             end
0187         end
0188         disp(sprintf(pformat,header_width,name,bayestopt_.p1(ip),post_mean,hpd_interval, ...
0189                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
0190         if TeX
0191             TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
0192                     bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
0193         end
0194         ip = ip+1;
0195     end
0196     if TeX
0197         TeXEnd(fid,3,'standard deviation of measurement errors');        
0198     end
0199 end
0200 if ncx
0201     type = 'shocks_corr';
0202     if TeX
0203         fid = TeXBegin(OutputDirectoryName,M_.fname,4,'correlation of structural shocks');
0204     end
0205     disp(' ')
0206     disp('correlation of shocks')
0207     disp(tit2)
0208     ip = nvx+nvn+1;
0209     for i=1:ncx
0210         if options_.mh_replic
0211             Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0212             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0213                 posterior_moments(Draws,1,options_.mh_conf_sig);
0214             k1 = estim_params_.corrx(i,1);
0215             k2 = estim_params_.corrx(i,2);
0216             name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
0217             NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
0218             oo_ = Filloo(oo_,NAME,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
0219             M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
0220             M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
0221         else
0222             try
0223                 k1 = estim_params_.corrx(i,1);
0224                 k2 = estim_params_.corrx(i,2);
0225                 name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
0226                 NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
0227                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,NAME,type);
0228             catch
0229                 Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0230                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0231                     posterior_moments(Draws,1,options_.mh_conf_sig);
0232                 k1 = estim_params_.corrx(i,1);
0233                 k2 = estim_params_.corrx(i,2);
0234                 name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
0235                 NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
0236                 oo_ = Filloo(oo_,NAME,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
0237                 M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
0238                 M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
0239             end
0240         end
0241         disp(sprintf(pformat, header_width,name,bayestopt_.p1(ip),post_mean,hpd_interval, ...
0242                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
0243         if TeX
0244             TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
0245                     bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
0246         end
0247         ip = ip+1;
0248     end
0249     if TeX
0250         TeXEnd(fid,4,'correlation of structural shocks');
0251     end
0252 end
0253 if ncn
0254     type = 'measurement_errors_corr';
0255     if TeX
0256         fid = TeXBegin(OutputDirectoryName,M_.fname,5,'correlation of measurement errors');
0257     end
0258     disp(' ')
0259     disp('correlation of measurement errors')
0260     disp(tit2)
0261     ip = nvx+nvn+ncx+1;
0262     for i=1:ncn
0263         if options_.mh_replic
0264             Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0265             [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0266                 posterior_moments(Draws,1,options_.mh_conf_sig);
0267             k1 = estim_params_.corrn(i,1);
0268             k2 = estim_params_.corrn(i,2);
0269             name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
0270             NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
0271             oo_ = Filloo(oo_,NAME,type,post_mean,hpd_interval,...
0272                          post_median,post_var,post_deciles,density);
0273         else
0274             try
0275                 k1 = estim_params_.corrn(i,1);
0276                 k2 = estim_params_.corrn(i,2);
0277                 name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
0278                 NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
0279                 [post_mean,hpd_interval,post_var] = Extractoo(oo_,NAME,type);
0280             catch
0281                 Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
0282                 [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
0283                     posterior_moments(Draws,1,options_.mh_conf_sig);
0284                 k1 = estim_params_.corrn(i,1);
0285                 k2 = estim_params_.corrn(i,2);
0286                 name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
0287                 NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
0288                 oo_ = Filloo(oo_,NAME,type,post_mean,hpd_interval,...
0289                              post_median,post_var,post_deciles,density);
0290             end
0291         end
0292         disp(sprintf(pformat, header_width,name,bayestopt_.p1(ip),post_mean,hpd_interval, ...
0293                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
0294         if TeX
0295             TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
0296                     bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);            
0297         end
0298         ip = ip+1;
0299     end
0300     if TeX
0301         TeXEnd(fid,5,'correlation of measurement errors');        
0302     end
0303 end
0304 
0305 
0306 %
0307 %% subfunctions:
0308 %
0309 function fid = TeXBegin(directory,fname,fnum,title)
0310 TeXfile = [directory '/' fname '_Posterior_Mean_' int2str(fnum) '.TeX'];
0311 fidTeX = fopen(TeXfile,'w');
0312 fprintf(fidTeX,'%% TeX-table generated by Dynare.\n');
0313 fprintf(fidTeX,['%% RESULTS FROM METROPOLIS HASTINGS (' title ')\n']);
0314 fprintf(fidTeX,['%% ' datestr(now,0)]);
0315 fprintf(fidTeX,' \n');
0316 fprintf(fidTeX,' \n');
0317 fprintf(fidTeX,'\\begin{table}\n');
0318 fprintf(fidTeX,'\\centering\n');
0319 fprintf(fidTeX,'\\begin{tabular}{l|lcccccc} \n');
0320 fprintf(fidTeX,'\\hline\\hline \\\\ \n');
0321 fprintf(fidTeX,['  & Prior distribution & Prior mean  & Prior ' ...
0322                 's.d. & Posterior mean & Posterior s.d.  & HPD inf & HPD sup\\\\ \n']);
0323 fprintf(fidTeX,'\\hline \\\\ \n');
0324 fid = fidTeX;
0325 
0326 
0327 function TeXCore(fid,name,shape,priormean,priorstd,postmean,poststd,hpd)
0328 fprintf(fid,['$%s$ & %s & %7.3f & %6.4f & %7.3f& %6.4f & %7.4f & %7.4f \\\\ \n'],... 
0329         name,...
0330         shape,...
0331         priormean,...
0332         priorstd,...
0333         postmean,...
0334         poststd,...
0335         hpd(1),...
0336         hpd(2));
0337 
0338 
0339 function TeXEnd(fid,fnum,title)
0340 fprintf(fid,'\\hline\\hline \n');
0341 fprintf(fid,'\\end{tabular}\n ');    
0342 fprintf(fid,['\\caption{Results from Metropolis-Hastings (' title ')}\n ']);
0343 fprintf(fid,['\\label{Table:MHPosterior:' int2str(fnum)  '}\n']);
0344 fprintf(fid,'\\end{table}\n');
0345 fprintf(fid,'%% End of TeX file.\n');
0346 fclose(fid);
0347 
0348 
0349 function oo = Filloo(oo,name,type,postmean,hpdinterval,postmedian,postvar,postdecile,density)
0350 eval(['oo.posterior_mean.' type '.' name ' = postmean;']);
0351 eval(['oo.posterior_hpdinf.' type '.' name ' = hpdinterval(1);']); 
0352 eval(['oo.posterior_hpdsup.' type '.' name ' = hpdinterval(2);']);      
0353 eval(['oo.posterior_median.' type '.' name ' = postmedian;']);
0354 eval(['oo.posterior_variance.' type '.' name ' = postvar;']);
0355 eval(['oo.posterior_deciles.' type '.' name ' = postdecile;']);
0356 eval(['oo.posterior_density.' type '.' name ' = density;']);
0357 
0358 function [post_mean,hpd_interval,post_var] = Extractoo(oo,name,type)
0359 hpd_interval = zeros(2,1);
0360 eval(['post_mean = oo.posterior_mean.' type '.' name ';']);
0361 eval(['hpd_interval(1) = oo.posterior_hpdinf.' type '.' name ';']); 
0362 eval(['hpd_interval(2) = oo.posterior_hpdsup.' type '.' name ';']);
0363 eval(['post_var = oo.posterior_variance.' type '.' name ';']);

Generated on Mon 21-May-2012 02:42:43 by m2html © 2005