function PrintBVARPosteriorSampling(DSGEModel,CurrINI,LogMargs,PsiParData,PiParData,OmegaParData,PsiMode,PiMode,OmegaMode,NumPostDraws,NumBurnin,OverallTightness,CrossEqTightness,HarmonicLagDecay,Amatrix,qDF,NewStartYear,NewStartPeriod,GibbsDrawTime)
% PrintBVARPosteriorSampling: Prints the results from the posterior sampler of the Bayesian VAR to file.
%
% USAGE:
%
%       PrintBVARPosteriorSampling(DSGEModel,CurrINI,LogMargs,PsiParData,PiParData,OmegaParData,...
%                                 PsiMode,PiMode,OmegaMode,NumPostDraws,NumBurnin,OverallTightness,...
%                                 CrossEqTightness,HarmonicLagDecay,Amatrix,qDF,NewStartYear,...
%                                 NewStartPeriod,GibbsDrawTime)
%
% REQUIRED INPUT: DSGEModel (structure) whose fields contain all required information
%                           about the model, such as input files, variable and sample
%                           data. Actual observations are not stored in this file.
%
%                 CurrINI (structure), with initialization information.
%
%                 LogMargs (matrix) with sample size in the first column, estimates of the
%                           log marginal likelihood in the second column, and numerical
%                           standard errors in the third column.
%
%                 PsiParData (structure) with the fields "SampleMean", "SampleStd", "SampleStdErr",
%                           "SampleMedian", and "Quantiles". These results concern the Psi parameters.
%
%                 PiParData (structure) with the same fields as PsiParData, but with results for
%                           the Pi parameters.
%
%                 OmegaParData (structure) with the same fields as PsiParData, but with results for
%                           the Omega parameters.
%
%                 PsiMode (matrix) of dimension pxq with the posterior mode estimate of the steady
%                           state parameters,
%
%                 PiMode (matrix) of dimension pxpk with the posterior mode estimate of the parameters
%                           on lagged endogenous variables.
%
%                 OmegaMode (matrix) with the posterior mode estimate of the residual covariance matrix.
%
%                 NumPostDraws (integer) with the number of posterior draws.
%
%                 NumBurnin (integer) with the number of burn-in draws for the posterior sampler.
%
%                 OverallTightness (scalar) with the overall tightness of the prior on Pi.
%
%                 CrossEqTightness (scalar) with the cross-equation tightness of the prior on Pi.
%
%                 HarmonicLagDecay (scalar) with the harmonic lag decay of the prior on Pi.
%
%                 Amatrix (pxp) matrix with the prior data on Omega.
%
%                 qDF (integer) with the number of degrees of freedom of the prior on Omega.
%
%                 NewStartYear (string vector) with the start year for the estimation.
%
%                 NewStartPeriod (string vector) with the start period for the estimation.
%
%                 GibbsDrawTime (time) with the length it took the posterior sampler to finish from where
%                           it started.
%
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2009 European Central Bank.
%
%                       First version: February 1, 2007.
%                        This version: January 21, 2009.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 2-2-2007:   Continued work on the function at the sounds of Pallas.
%
% * 5-2-2007:   Fixed the sample size presentation for the marginal likelihood estimates.
%
% * 13-6-2007:  Made sure that the strings "head_prior_mean_str" and "b_prt_str" are properly
%               initialized when printing Omega estimation results.
%
% * 13-11-2007: Updated the documentation.
%
% * 23-5-2008:  Updated the documentation.
%
% * 21-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

TextFile = [DSGEModel.OutputDirectory '\BVAR-GibbsSummary-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NumPostDraws,'%0.0f') '.txt'];
%
fid = fopen(TextFile,'wt');
%
% Setup the names of the variables used by the BVAR estimation.
%
VariableNamesBVAR = DSGEModel.VariableNames(DSGEModel.BVARY,:);
if isempty(DSGEModel.BVARX)==0;
   XVariableNamesBVAR = DSGEModel.XVariableNames(DSGEModel.BVARX,:);
else;
   XVariableNamesBVAR = '';
end;
%
% some useful data parameters.
%
p = length(DSGEModel.BVARY);
k = DSGEModel.BVARLags;
q = length(DSGEModel.BVARX);
%
% the header information
%
fprintf(fid,'********************************************************************************\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'*                     B A Y E S I A N   V A R   M O D E L                      *\n');
fprintf(fid,'*           P O S T E R I O R   S A M P L I N G   W I T H   Y A D A            *\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'********************************************************************************\n\n');
%
% Model info
%
fprintf(fid,'Name of model:                    %s\n',DSGEModel.NameOfModel);
fprintf(fid,'Data construction file:           %s\n',GetFilename(DSGEModel.DataConstructionFile));
fprintf(fid,'Steady state prior file:          %s\n\n',GetFilename(DSGEModel.SteadyStatePriorFile));
%
fprintf(fid,'********************************************************************************\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'*              S A M P L E   &   D A T A   I N F O R M A T I O N               *\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'********************************************************************************\n\n');
%
% Sample information
%
fprintf(fid,'Estimation sample:                %s:%s -- %s:%s\n',NewStartYear,NewStartPeriod,DSGEModel.SubEndYear,DSGEModel.SubEndPeriod);
fprintf(fid,'Number of lags of the VAR model:  %0.0f\n',k);
fprintf(fid,'Number of endogenous variables:   %0.0f\n',p);
fprintf(fid,'Number of exogenous variables:    %0.0f\n\n',q);
%
% variable names
%
VarNames = StrWrap(vecstr(VariableNamesBVAR,' '),60);
fprintf(fid,'Names of endogenous variables:    %s\n',VarNames(1,:));
for i=2:size(VarNames,1);
   if strcmp(StringTrim(VarNames(i,:)),'')==0;
      fprintf(fid,'                                  %s\n',VarNames(i,:));
   end;
end;
if q>0;
   XVarNames = StrWrap(vecstr(XVariableNamesBVAR,' '),60);
   fprintf(fid,'Names of exogenous variables:     %s\n',XVarNames(1,:));
   for i=2:size(XVarNames,1);
      if strcmp(StringTrim(XVarNames(i,:)),'')==0;
         fprintf(fid,'                                  %s\n',XVarNames(i,:));
      end;
   end;
end;
%
% General information about the prior
%
fprintf(fid,'\n********************************************************************************\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'*                      P R I O R   I N F O R M A T I O N                       *\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'********************************************************************************\n\n');
if q>0;
   fprintf(fid,'Psi:                              Normal distribution\n');
end;
if DSGEModel.PriorType==1;
   %
   % we have a Minnesota-style prior on Pi
   %
   fprintf(fid,'Pi:                               Minnesota-style\n');
   fprintf(fid,'Overall tightness:                %0.6f\n',OverallTightness);
   fprintf(fid,'Cross-equation tightness:         %0.6f\n',CrossEqTightness);
   fprintf(fid,'Harmonic lag decay:               %0.6f\n',HarmonicLagDecay);
elseif DSGEModel.PriorType==2;
   %
   % we have a normal conditional on redidual covariance matrix prior on Pi
   %
   fprintf(fid,'Pi:                               Normal conditional on Omega\n');
   fprintf(fid,'Overall tightness:                %0.6f\n',OverallTightness);
   fprintf(fid,'Harmonic lag decay:               %0.6f\n',HarmonicLagDecay);
else;
   %
   % we have a diffuse prior on Pi
   %
   fprintf(fid,'Pi:                               Diffuse\n');
end;
if DSGEModel.OmegaPriorType==1;
   %
   % We have a diffuse prior on Omega
   %
   fprintf(fid,'Omega:                            Diffuse\n');
else;
   %
   % we have an inverse Wishart prior on Omega
   %
   fprintf(fid,'Omega:                            Inverse Wishart\n');
   if DSGEModel.WishartType==1;
      %
      % We use the ML estimate of Omega for the A matrix
      %
      fprintf(fid,'A matrix:                         ML estimate of Omega\n');
   else;
      %
      % we have a common variance tightness hyperparameters
      %
      fprintf(fid,'Variance tightness:               %0.6f\n',Amatrix(1,1));
   end;
   fprintf(fid,'Degrees of freedom:               %0.0f\n',qDF);
end;
%
% posterior sampler information
%
fprintf(fid,'\n********************************************************************************\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'*            P O S T E R I O R   D R A W S   I N F O R M A T I O N             *\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'********************************************************************************\n\n');
%
fprintf(fid,'Number of posterior draws:             %0.0f\n',NumPostDraws);
fprintf(fid,'Number of burn-in draws:               %0.0f\n',NumBurnin);
fprintf(fid,'Number of parameters to estimate:      %0.0f\n',p*(q+(p*k)+((p+1)/2)));
fprintf(fid,'Posterior sampling time:               %s (hrs:min:sec)\n\n',SecondsToHoursStr(GibbsDrawTime));
%
%
%
prt_val = ['%' num2str(9+CurrINI.decimals,'%0.0f') '.' num2str(CurrINI.decimals,'%0.0f') 'f'];
%
% write results for log marginal likelihoods
%
fprintf(fid,'********************************************************************************\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'*         M A R G I N A L   L I K E L I H O O D   E S T I M A T I O N          *\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'********************************************************************************\n\n');
%
fprintf(fid,'======================================================================\n');
fprintf(fid,' Chib''s BMLI          sample size    point estimate        std. error\n');
for i=1:size(LogMargs,1);
   fprintf(fid,'----------------------------------------------------------------------\n');
   fprintf(fid,['                          %7.0f   ' prt_val '   ' prt_val '\n'],LogMargs(i,1),LogMargs(i,2),LogMargs(i,3));
end;
fprintf(fid,'======================================================================\n\n');
%
% deal with the parameters
%
prt_val = ['%' num2str(6+CurrINI.decimals,'%0.0f') '.' num2str(CurrINI.decimals,'%0.0f') 'f'];
fprintf(fid,'********************************************************************************\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'*                             P A R A M E T E R S                              *\n');
fprintf(fid,'*     P O I N T   E S T I M A T E S   F R O M   T H E   P O S T E R I O R      *\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'********************************************************************************\n\n');
%
% The steady state
%
if q>0;
   %
   % get the posterior median and mean estimates
   %
   PsiMean = reshape(PsiParData.SampleMean,p,q);
   PsiMedian = reshape(PsiParData.SampleMedian,p,q);
   fprintf(fid,'STEADY STATE PARAMETERS - Psi\n\n');
   head_str = '  Eq\Var  ';
   add_head_str = '';
   if size(XVariableNamesBVAR,2)<10;
      add_head_str = SpaceStr(10-size(XVariableNamesBVAR,2));
   end;
   head_prior_mean_str = '          ';
   head_prior_var_str = '          ';
   add_name_str = '';
   if size(VariableNamesBVAR,2)<10;
      add_name_str = SpaceStr(10-size(VariableNamesBVAR,2));
   end;
   prt_str = '%s%s';
   b_prt_str = '          ';
   for i=1:q;
      head_str = [head_str sprintf('       %s%s',add_head_str,XVariableNamesBVAR(i,:))];
      head_prior_mean_str = [head_prior_mean_str '           median'];
      head_prior_var_str = [head_prior_var_str '             mode'];
      prt_str = [prt_str '     ' prt_val];
      b_prt_str = [b_prt_str '     ' prt_val];
   end;
   TopStr = MultiCharStr('=',length(head_str)+2);
   MidStr = MultiCharStr('-',length(head_str)+2);
   fprintf(fid,'%s\n',TopStr);
   fprintf(fid,'%s\n',head_str);
   fprintf(fid,'%s\n',strrep(head_prior_var_str,'mode','mean'));
   fprintf(fid,'%s\n',head_prior_mean_str);
   fprintf(fid,'%s\n',head_prior_var_str);
   for i=1:p;
      fprintf(fid,'%s\n',MidStr);
      fprintf(fid,[prt_str '\n'],add_name_str,VariableNamesBVAR(i,:),PsiMean(i,:)');
      fprintf(fid,[b_prt_str '\n'],PsiMedian(i,:)');
      fprintf(fid,[b_prt_str '\n'],PsiMode(i,:)');
   end;
   fprintf(fid,'%s\n\n',TopStr);
end;
%
% The parameters on lagged values
%
fprintf(fid,'PARAMETERS ON LAGGED ENDOGENOUS VARIABLES - Pi\n\n');
%
% get the posterior mean and median estimates
%
PiMean = reshape(PiParData.SampleMean,p,p*k);
PiMedian = reshape(PiParData.SampleMedian,p,p*k);
for l=1:k;
   PilMean = PiMean(:,1+p*(l-1):p*l);
   PilMedian = PiMedian(:,1+p*(l-1):p*l);
   PilMode = PiMode(:,1+p*(l-1):p*l);
   if l==1;
      add_name_str = '';
      if size(VariableNamesBVAR,2)<10;
         add_name_str = SpaceStr(10-size(VariableNamesBVAR,2));
      end;
      head_str = '  Eq\Var  ';
      prt_str = '%s%s';
      head_prior_mean_str = '          ';
      head_prior_var_str = '          ';
      b_prt_str = '          ';
      for i=1:p;
         head_str = [head_str sprintf('       %s%s',add_name_str,VariableNamesBVAR(i,:))];
         prt_str = [prt_str '     ' prt_val];
         head_prior_mean_str = [head_prior_mean_str '           median'];
         head_prior_var_str = [head_prior_var_str '             mode'];
         b_prt_str = [b_prt_str '     ' prt_val];
      end;
      TopStr = MultiCharStr('=',length(head_str)+2);
      MidStr = MultiCharStr('-',length(head_str)+2);
   end;
   %
   % ready to print
   %
   fprintf(fid,'%s\n',TopStr);
   fprintf(fid,'Lag: %2.0f\n\n',l);
   fprintf(fid,'%s\n',head_str);
   if l==1;
      fprintf(fid,'%s\n',strrep(head_prior_var_str,'mode','mean'));
      fprintf(fid,'%s\n',head_prior_mean_str);
      fprintf(fid,'%s\n',head_prior_var_str);
   end;
   for i=1:p;
      fprintf(fid,'%s\n',MidStr);
      fprintf(fid,[prt_str '\n'],add_name_str,VariableNamesBVAR(i,:),PilMean(i,:)');
      fprintf(fid,[b_prt_str '\n'],PilMedian(i,:)');
      fprintf(fid,[b_prt_str '\n'],PilMode(i,:)');
   end;
   fprintf(fid,'%s\n\n',TopStr);
end;
%
% covariance parameters
%
fprintf(fid,'RESIDUAL COVARIANCE PARAMETERS - Omega\n\n');
%
% get the posterior mean and median
%
Dp = DuplicationMatrix(p);
OmegaMean = reshape(Dp*OmegaParData.SampleMean',p,p);
OmegaMedian = reshape(Dp*OmegaParData.SampleMedian',p,p);
%
head_str = '  Eq\Res  ';
prt_str = '%s%s';
head_prior_mean_str = '          ';
b_prt_str = '          ';
for i=1:p;
   head_str = [head_str sprintf('  %sRes(%s)',add_name_str,VariableNamesBVAR(i,:))];
   prt_str = [prt_str '     ' prt_val];
   head_prior_mean_str = [head_prior_mean_str '             mean'];
   b_prt_str = [b_prt_str '     ' prt_val];
end;
%
fprintf(fid,'%s\n',TopStr);
fprintf(fid,'%s\n',head_str);
fprintf(fid,'%s\n',head_prior_mean_str);
fprintf(fid,'%s\n',strrep(head_prior_mean_str,'  mean','median'));
fprintf(fid,'%s\n',strrep(head_prior_mean_str,'mean','mode'));
for i=1:p;
   fprintf(fid,'%s\n',MidStr);
   fprintf(fid,[prt_str '\n'],add_name_str,VariableNamesBVAR(i,:),OmegaMean(i,:)');
   fprintf(fid,[b_prt_str '\n'],OmegaMedian(i,:)');
   fprintf(fid,[b_prt_str '\n'],OmegaMode(i,:)');
end;
fprintf(fid,'%s\n\n',TopStr);
%
% Some information about the filename and its location
%
fprintf(fid,'Directory for file: %s\n',GetPath(TextFile));
fprintf(fid,'Name of file:       %s\n',GetFilename(TextFile));
fprintf(fid,'Output created on:  %s\n\n',datestr(now,0));
fclose(fid);

%
% end of PrintBVARPosteriorSampling.m
%
