function GraphBVARPosteriorDensity(ParamStr,PostSample,NumBurnin,DSGEModel,CurrINI)
% GraphBVARPosteriorDensity: Creates plots of the posterior densities for a certain subset of the
%                            parameters of the Bayesian VAR.
%
% USAGE:
%
%       GraphBVARPosteriorDensity(ParamStr,PostSample,NumBurnin,DSGEModel,CurrINI)
%
% REQUIRED INPUT: ParamStr (string vector) taking on the values "Psi", "Pi", or "Omega". These string values
%                          correspond to the steady state parameters, the parameters on lagged endogenous
%                          variables, and the parameters of the residual covariance matrix.
%
%                 PostSample (matrix) with all the posterior draws of the specified groups of parameters.
%
%                 NumBurnin (integer), with the number of burn-in draws.
%
%                 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.
%
%
%                       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 6, 2007.
%                        This version: March 2, 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:
%
% * 07-02-2007: Added the call to the GraphBVARPriorDensity function.
%
% * 08-02-2007: Updated the documentation and added a second "drawnow" call at the end of the figure loop.
%
% * 18-04-2007: Changed the plot-style for the prior from dotted to dashed.
%
% * 09-11-2007: Updated the documentation.
%
% * 21-11-2007: Updated the documentation.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 23-05-2008: Updated the documentation.
%
% * 19-12-2008: Updated the documentation.
%
% * 02-03-2009: The function now check which gaussian kernel density estimation function
%               the user has when "DSGEModel.PosteriorDensityValue" is 1.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

[NumIter,NumParam] = size(PostSample);
%
% Try to load the posterior mode estimates
%
ModeFile = [DSGEModel.OutputDirectory '\mode\PosteriorMode-BVAR-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '.mat'];
ErrorStr = '';
try;
   ModeData = load(ModeFile);
catch;
   ErrorStr = ['Unable to load the file "' ModeFile '". Message caught is: ' lasterr];
end;
if isempty(ErrorStr)==0;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% data loaded fine
%
PsiMode = ModeData.Psi;
PiMode = vec(ModeData.Pi);
OmegaMode = ModeData.Omega;
%
% first we display a wait dialog
%
if DSGEModel.PosteriorDensityValue>2;
   WaittxtAdd = [' Estimation started at: ' StringTrim(datestr(now,14)) '. This may take a very long time...'];
   ShowCancel = 1;
else;
   WaittxtAdd = [' Estimation started at: ' StringTrim(datestr(now,14)) '.'];
   ShowCancel = 0;
end;
txt = ['Please wait while YADA computes kernel density estimates of the marginal posterior densities of the ' ParamStr ' parameters.' WaittxtAdd ' Current parameter: - 1 out of ' num2str(NumParam,'%0.0f') '.'];
WaitHandle = WaitDLG(txt,'information','Kernel Density Estimation - Posterior',500,150,CurrINI,ShowCancel);
WaitControls = get(WaitHandle,'UserData');
drawnow;
pause(0.02);
%
% now we try to construct the x-axis and y-axis data for the posterior densities
%
StopEstimation = 0;
for i=1:NumParam;
   txt = ['Please wait while YADA computes kernel density estimates of the marginal posterior densities of the ' ParamStr ' parameters.' WaittxtAdd ' Current parameter: ' num2str(i,'%0.0f') ' out of ' num2str(NumParam,'%0.0f') '.'];
   set(WaitControls.text,'String',txt);
   drawnow;
   if DSGEModel.PosteriorDensityValue==1;
      if exist('YADAksdensity')==2;
         [f,x] = YADAksdensity(PostSample(NumBurnin+1:NumIter,i),'npoints',256);
      else;
         if exist('ksdensity')==2;
            [f,x] = ksdensity(PostSample(NumBurnin+1:NumIter,i),'npoints',256);
         else;
            [x,f] = KernelDensityPlot(PostSample(NumBurnin+1:NumIter,i),'knorm');
         end;
      end;
   elseif DSGEModel.PosteriorDensityValue==2;
      [x,f] = YADAMHDensEst(PostSample(NumBurnin+1:NumIter,i),2^7,0,'gaussian',[]);
   elseif DSGEModel.PosteriorDensityValue==3;
      [x,f] = YADAMHDensEst(PostSample(NumBurnin+1:NumIter,i),2^7,-1,'gaussian',WaitControls.cancel);
   elseif DSGEModel.PosteriorDensityValue==4;
      [x,f] = YADAMHDensEst(PostSample(NumBurnin+1:NumIter,i),2^7,-2,'gaussian',WaitControls.cancel);
   end;
   if isempty(x)==1;
      StopEstimation = 1;
      break;
   end;
   if i==1;
      Xposterior = zeros(NumParam,length(x));
      Yposterior = zeros(NumParam,length(x));
   end;
   if size(x,1)>size(x,2);
      Xposterior(i,:) = x';
   else;
      Xposterior(i,:) = x;
   end;
   if size(f,1)>size(f,2);
      Yposterior(i,:) = f';
   else;
      Yposterior(i,:) = f;
   end;
end;
%
% close the wait dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);
%
if StopEstimation==1;
   return;
end;
%
if (strcmp(lower(ParamStr),'psi')==1)|((strcmp(lower(ParamStr),'pi')==1)&(DSGEModel.PriorType==1))|((strcmp(lower(ParamStr),'pi')==1)&(DSGEModel.PriorType==2)&(DSGEModel.OmegaPriorType==2))|((strcmp(lower(ParamStr),'omega')==1)&(DSGEModel.OmegaPriorType==2));
   [Xprior,Yprior] = GraphBVARPriorDensity(ParamStr,DSGEModel,CurrINI);
else;
   Xprior = [];
   Yprior = [];
end;
%
% Determine the variable names
%
VariableNamesBVAR = DSGEModel.VariableNames(DSGEModel.BVARY,:);
if isempty(DSGEModel.BVARX)==0;
   XVariableNamesBVAR = DSGEModel.XVariableNames(DSGEModel.BVARX,:);
else;
   XVariableNamesBVAR = '';
end;
p = length(DSGEModel.BVARY);
%
% Directory where the plots are stored
%
GraphDirectory = [DSGEModel.OutputDirectory '\graphs'];
%
% Compute the number of figures (max 16 plots per graph)
%
NumFigs = ceil(NumParam/16);
%
if NumFigs==1;
   [NumRows,NumCols] = VillaniConstructOptimalSubplot(NumParam);
else;
   NumRows = 4;
   NumCols = 4;
end;
NumPlots = 0;
pos = [(CurrINI.scrsz(3)-1000)/2 (CurrINI.scrsz(4)-680)/2 1000 650];
%
% now we can start the loops
%
IndexOne = 0;
IndexTwo = 1;
IndexThree = 1;
FigParmStr = '';
ModeEstimate = 0;
for i=1:NumFigs;
   if i>1;
      pos(1) = pos(1)+22;
      pos(2) = pos(2)-22;
      if pos(1)+1000>CurrINI.scrsz(3);
         pos(1) = 0;
      end;
      if pos(2)<=32;
         pos(2) = CurrINI.scrsz(4)-700;
      end;
   end;
   FigHandles{i} = figure('Visible','off','Position',pos);
   if NumFigs==1;
      set(FigHandles{i},'Name',['Posterior Densities of ' ParamStr]);
   else;
      set(FigHandles{i},'Name',['Posterior Densities of ' ParamStr ' - Figure ' num2str(i,'%0.0f')]);
   end;
   SetFigureProperties(FigHandles{i},CurrINI);
   %
   k = 0;
   for r=1:NumRows;
      for c=1:NumCols;
         NumPlots = NumPlots+1;
         if NumPlots<=NumParam;
            k = k+1;
            subplot(NumRows,NumCols,k);
            if strcmp(lower(ParamStr),'psi')==1;
               %
               % determine the variables
               %
               if IndexOne==p;
                  IndexOne = 1;
                  IndexTwo = IndexTwo+1;
               else;
                  IndexOne = IndexOne+1;
               end;
               FigParmStr = ['Psi(' StringTrim(VariableNamesBVAR(IndexOne,:)) ',' StringTrim(XVariableNamesBVAR(IndexTwo,:)) ')'];
               ModeEstimate = PsiMode(IndexOne,IndexTwo);
            elseif strcmp(lower(ParamStr),'pi')==1;
               if IndexOne==p;
                  IndexOne = 1;
                  IndexTwo = IndexTwo+1;
                  %
                  % if we should add a lag or not
                  %
                  if IndexTwo>p;
                     IndexTwo = 1;
                     IndexThree = IndexThree+1;
                  end;
               else;
                  IndexOne = IndexOne+1;
               end;
               FigParmStr = ['Pi(' StringTrim(VariableNamesBVAR(IndexOne,:)) ',' StringTrim(VariableNamesBVAR(IndexTwo,:)) ';Lag=' num2str(IndexThree,'%0.0f') ')'];
               ModeEstimate = PiMode(NumPlots,1);
            else;
               if IndexOne==p;
                  IndexOne = IndexTwo+1;
                  IndexTwo = IndexTwo+1;
               else;
                  IndexOne = IndexOne+1;
               end;
               FigParmStr = ['Omega(' StringTrim(VariableNamesBVAR(IndexOne,:)) ',' StringTrim(VariableNamesBVAR(IndexTwo,:)) ')'];
               ModeEstimate = OmegaMode(IndexOne,IndexTwo);
            end;
            %
            % plot the posterior density
            %
            PosteriorHandle = plot(Xposterior(NumPlots,:),Yposterior(NumPlots,:),'-','Color',[0 0 0.75],'LineWidth',2);
            hold('on');
            SetAxesFonts(CurrINI,'Title',FigParmStr);
            axis('tight');
            xlim = [min(PostSample(NumBurnin+1:NumIter,NumPlots)) max(PostSample(NumBurnin+1:NumIter,NumPlots))];
            add_x = 0.1*(xlim(2)-xlim(1));
            x_lim = [xlim(1)-add_x xlim(2)+add_x];
            set(gca,'XLim',x_lim);
            ylim = get(gca,'YLim');
            add_on = 0.1*(ylim(2)-ylim(1));
            ylim = [0 ylim(2)+add_on];
            %
            % check if we should plot a prior
            %
            if isempty(Xprior)==0;
               PriorHandle = plot(Xprior(NumPlots,:),Yprior(NumPlots,:),'--','Color',[0.75 0 0],'LineWidth',2);
               hold('on');
               %
               % check 
               %
               if (strcmp(lower(ParamStr),'omega')==1)&(IndexOne==IndexTwo);
                  x_lim(1) = 0;
               end;
               if ylim(2)<max(Yprior(NumPlots,:));
                  ylim(2) = max(Yprior(NumPlots,:))+add_on;
               end;
            end;
            %
            % plot the posterior mode estimate
            %
            ModeHandle = plot(ModeEstimate*ones(1,2),ylim,'-','Color',[0 0 0],'LineWidth',1);
            %
            set(gca,'XLim',x_lim,'YLim',ylim);
            hold('off');
         end;
      end;
   end;
   %
   % create a legend
   %
   if isempty(Xprior)==0;
      legHandle = legend([PriorHandle PosteriorHandle ModeHandle],strvcat('Prior','Posterior','Estimated Mode'));
   else;
      legHandle = legend([PosteriorHandle ModeHandle],strvcat('Posterior','Estimated Mode'));
   end;
   %
   SetAxesFonts(CurrINI,legHandle);
   %
   % transform units temporarily into pixels
   %
   OldAxesUnits = get(gca,'Units');
   OldLegUnits = get(legHandle,'Units');
   set(gca,'Units','pixels');
   set(legHandle,'Units','pixels');
   %
   % check if we can move legend
   %
   LegPos = get(legHandle,'Position');
   AxPos = get(gca,'Position');
   if LegPos(1)<AxPos(1)+AxPos(3);
      if AxPos(1)+AxPos(3)+10+LegPos(3)<1000;
         LegPos(1) = AxPos(1)+AxPos(3)+10;
         set(legHandle,'Position',LegPos);
      end;
   end;
   set(gca,'Units',OldAxesUnits);
   set(legHandle,'Units',OldLegUnits);
   %
   % save it to disk before
   %
   if DirExist(GraphDirectory)==1;
      if NumFigs>1;
         file = [GraphDirectory '\PosteriorDensity-BVAR-' ParamStr '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NumIter,'%0.0f') '-' num2str(i,'%0.0f') '.fig'];
      else;
         file = [GraphDirectory '\PosteriorDensity-BVAR-' ParamStr '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NumIter,'%0.0f') '.fig'];
      end;
      saveas(FigHandles{i},file);
   end;
   %
   % Fix the toolbar cdata
   %
   SetToolbarImages(FigHandles{i},CurrINI.images);
   %
   % Display the current figure
   %
   set(FigHandles{i},'Visible','on','CloseRequestFcn','delete(gcf); drawnow; pause(0.02);');
   drawnow;
   drawnow;
end;

%
% end of GraphBVARPosteriorDensity.m
%
