function [Xdata,Ydata] = GraphBVARPriorDensity(ParamStr,DSGEModel,CurrINI)
% GraphBVARPriorDensity: Plots prior densities for a subset of the parameters of the
%                        Bayesian VAR.
%
% USAGE:
%
%       [Xdata,Ydata] = GraphBVARPriorDensity(ParamStr,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.
%
%                  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.
%
% OPTIONAL OUTPUT: Xdata (matrix) with dimension NumParamxNumPoints. It contains the x-axis data for all
%                            the densities.
%
%                  Ydata (matrix) with dimension NumParamxNumPoints. It contains the y-axis data for all
%                            the densities.
%
%
%                       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 9, 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: Finished the function.
%
% * 08-02-2007: Added a second "drawnow" call at the end of the figure loop.
%
% * 09-11-2007: Updated the documentation.
%
% * 21-11-2007: Updated the documentation.
%
% * 02-01-2008: Removed code using "DSGEModel.InitializeKalman".
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 23-05-2008: Updated the documentation.
%
% * 19-12-2008: Updated the documentation.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy command.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% collect the handle structure from the YADA gui
%
maingui = findobj('Type','figure','Tag','YADA');
controls = get(maingui,'UserData');
%
% number of grid points
%
NumGridPointsStr = get(controls.optimize.numbergridpoints,'String');
NumGridPoints = str2num(StringTrim(NumGridPointsStr(DSGEModel.NumberOfGridPoints,:)));
NumGridPointsHalf = NumGridPoints/2;
%
% number of endogenous variables and the number of lags
%
p = length(DSGEModel.BVARY);
k = DSGEModel.BVARLags;
q = length(DSGEModel.BVARX);
if strcmp(lower(ParamStr),'psi')==0;
   %
   % Determine the sample
   %
   [FirstPeriod,LastPeriod] = CreateSubSample(DSGEModel);
   %
   % Determine the first observation to use when computing log-likelihood with the Kalman filter
   %
   KalmanFirstObservation = DSGEModel.KalmanFirstObservation;
   %
   % Check if we should adjust this based on k, the lag order
   %
   if FirstPeriod+KalmanFirstObservation-1<k+1;
      KalmanFirstObservation = k+2-FirstPeriod;
   end;
   %
   % Setup all the basic data matrices
   %
   x = DSGEModel.Y(DSGEModel.BVARY,FirstPeriod+KalmanFirstObservation-1:LastPeriod);
   T = size(x,2);
   X = zeros(p*k,T);
   for i=1:k;
      X(1+(p*(i-1)):p*i,:) = DSGEModel.Y(DSGEModel.BVARY,FirstPeriod+KalmanFirstObservation-1-i:LastPeriod-i);
   end;
   if isempty(DSGEModel.BVARX)==0;
      d = DSGEModel.X(DSGEModel.BVARX,FirstPeriod+KalmanFirstObservation-1:LastPeriod);
      dLag = zeros(q*k,T);
      for i=1:k;
         dLag(1+(q*(i-1)):q*i,:) = -DSGEModel.X(DSGEModel.BVARX,FirstPeriod+KalmanFirstObservation-1-i:LastPeriod-i);
      end;
   else;
      d = [];
      dLag = [];
   end;
   %
   % Compute the Maximum likelihood estimate of the residual covariance matrix of the VAR
   %
   Z = [d;X];
   BML = (x*Z')/(Z*Z');
   Epsilon = x-(BML*Z);
   OmegaML = (1/T)*Epsilon*Epsilon';
end;
%
if strcmp(lower(ParamStr),'psi')==1;
   %
   % Collect all the data for the stadey state prior
   %
   [stat,msg] = dos(['copy /Y "' DSGEModel.SteadyStatePriorFile '" "' pwd '\tmp"']);
   %
   % collect the steady state prior data
   %
   eval(['PriorMatrix = ' GetName(DSGEModel.SteadyStatePriorFile) '(DSGEModel.BVARY,DSGEModel.BVARX);']);
   %
   thetaPsi = PriorMatrix(:,1);
   OmegaPsi = diag(PriorMatrix(:,2).*PriorMatrix(:,2));
   if det(OmegaPsi)<1.0e-10;
      if nargout==2;
         Xdata = [];
         Ydata = [];
         return;
      end;
      txt = ['The determinant of the prior covariance matrix is: ' num2str(det(OmegaPsi)) '. This is too small. YADA has therefore aborted posterior mode estimation of the Bayesian VAR.'];
      About(txt,'information','YADA - BVAR Model Problem',200,500,CurrINI);
      return;
   end;
   %
   % compute prior parameter values
   %
   if DSGEModel.PriorKernel==0;
      XValues = zeros(length(thetaPsi),NumGridPoints+1);
      YValues = XValues;
      for i=1:length(thetaPsi);
         lb = thetaPsi(i)-(DSGEModel.GridWidth*sqrt(OmegaPsi(i,i)));
         dev = (thetaPsi(i)-lb)/NumGridPointsHalf;
         XValues(i,:) = (lb:dev:lb+(NumGridPoints*dev));
         YValues(i,:) = exp(logNormalPDF(XValues(i,:),thetaPsi(i),sqrt(OmegaPsi(i,i))));
      end;
   else;
      XValues = zeros(length(thetaPsi),512);
      YValues = XValues;
      for i=1:length(thetaPsi);
         z = NormalRndFcn(thetaPsi(i),sqrt(OmegaPsi(i,i)),5000);
         %
         % compute the kernel density values
         %
         [XValue,YValue] = KernelDensityPlot(z,DSGEModel.KernelDensityEstimator);
         %
         XValues(i,:) = XValue';
         YValues(i,:) = YValue';
      end;
   end;
elseif strcmp(lower(ParamStr),'pi')==1;
   %
   % the Pi parameters
   %
   %
   % Collect the hyperparameters
   %
   PriorLevelMeanStr = get(controls.bvar.priorlevelmean,'String');
   PriorLevelMean = str2num(StringTrim(PriorLevelMeanStr(DSGEModel.PriorLevelMeanValue,:)));
   %
   PriorDiffMeanStr = get(controls.bvar.priordiffmean,'String');
   PriorDiffMean = str2num(StringTrim(PriorDiffMeanStr(DSGEModel.PriorDiffMeanValue,:)));
   %
   thetaPi = zeros(p*p*k,1);
   DiffVar = 1-DSGEModel.levels;
   DiffVar = DiffVar(DSGEModel.BVARY);
   LevVar = DSGEModel.levels(DSGEModel.BVARY);
   Pi1 = diag((PriorLevelMean*LevVar)+(PriorDiffMean*DiffVar));
   thetaPi(1:p*p,1) = vec(Pi1);
   %
   OverallTightnessStr = get(controls.bvar.overalltightness,'String');
   OverallTightness = str2num(StringTrim(OverallTightnessStr(DSGEModel.OverallTightnessValue,:)));
   %
   CrossEqTightnessStr = get(controls.bvar.crosseqtightness,'String');
   CrossEqTightness = str2num(StringTrim(CrossEqTightnessStr(DSGEModel.CrossEqTightnessValue,:)));
   %
   HarmonicLagDecayStr = get(controls.bvar.harmoniclagdecay,'String');
   HarmonicLagDecay = str2num(StringTrim(HarmonicLagDecayStr(DSGEModel.HarmonicLagDecayValue,:)));
   if DSGEModel.PriorType==1;
      %
      % we have a Minnesota-style prior
      %
      OmegaPi = MinnesotaPrior(OverallTightness,CrossEqTightness,HarmonicLagDecay,diag(OmegaML),k);
      %
      % we now have the prior mean and the prior covariances.
      %
      if DSGEModel.PriorKernel==0;
         XValues = zeros(length(thetaPi),NumGridPoints+1);
         YValues = XValues;
         for i=1:length(thetaPi);
            lb = thetaPi(i)-(DSGEModel.GridWidth*sqrt(OmegaPi(i,i)));
            dev = (thetaPi(i)-lb)/NumGridPointsHalf;
            XValues(i,:) = (lb:dev:lb+(NumGridPoints*dev));
            YValues(i,:) = exp(logNormalPDF(XValues(i,:),thetaPi(i),sqrt(OmegaPi(i,i))));
         end;
      else;
         XValues = zeros(length(thetaPi),512);
         YValues = XValues;
         for i=1:length(thetaPi);
            z = NormalRndFcn(thetaPi(i),sqrt(OmegaPi(i,i)),5000);
            %
            % compute the kernel density values
            %
            [XValue,YValue] = KernelDensityPlot(z,DSGEModel.KernelDensityEstimator);
            %
            XValues(i,:) = XValue';
            YValues(i,:) = YValue';
         end;
      end;
   else;
      %
      % we have a normal conditional on Omega prior with an informative Omega prior.
      % The prior is therefore matrix t.
      %
      OmegaPi = NormalConditionalPrior(OverallTightness,HarmonicLagDecay,p,k);
      qDF = p+get(controls.bvar.wishartdf,'Value')+1;
      if DSGEModel.WishartType==1;
         %
         % We use the ML estimate of Omega for the A matrix
         %
         Amatrix = OmegaML;
      else;
         %
         % collect the variance tightness hyperparameter
         %
         VarianceTightnessStr = get(controls.bvar.variancetightness,'String');
         VarianceTightness = str2num(StringTrim(VarianceTightnessStr(DSGEModel.VarianceTightnessValue,:)));
         Amatrix = VarianceTightness*eye(p);
      end;
      %
      if DSGEModel.PriorKernel==0;
         XValues = zeros(length(thetaPi),NumGridPoints+1);
         YValues = XValues;
         num = 1;
         for j=1:p*k;
            q_jj = OmegaPi(j,j);
            for i=1:p;
               p_ii = 1/Amatrix(i,i);
               lb = thetaPi(num)-(DSGEModel.GridWidth*sqrt(q_jj/p_ii));
               dev = (thetaPi(num)-lb)/NumGridPointsHalf;
               XValues(num,:) = (lb:dev:lb+(NumGridPoints*dev));
               YValues(num,:) = exp(logStudentTPDF(XValues(num,:),thetaPi(num),p_ii,q_jj,qDF-p));
               num = num+1;
            end;
         end;
      else;
         XValues = zeros(length(thetaPi),512);
         YValues = XValues;
         num = 1;
         for j=1:p*k;
            q_jj = OmegaPi(j,j);
            for i=1:p;
               p_ii = 1/Amatrix(i,i);
               z = StudentTRndFcn(thetaPi(num),p_ii,q_jj,qDF-p,5000);
               %
               % compute the kernel density values
               %
%               [XValue,YValue] = KernelDensityPlot(z,DSGEModel.KernelDensityEstimator);
               [XValue,YValue] = YADAMHDensEst(z,2^9,0,'gaussian',[]);
               %
               XValues(num,:) = XValue';
               YValues(num,:) = YValue';
               num = num+1;
            end;
         end;
      end;
   end;
else;
   %
   % we're dealing with the Omega parameters
   %
   qDF = p+get(controls.bvar.wishartdf,'Value')+1;
   if DSGEModel.WishartType==1;
      %
      % We use the ML estimate of Omega for the A matrix
      %
      Amatrix = OmegaML;
   else;
      %
      % collect the variance tightness hyperparameter
      %
      VarianceTightnessStr = get(controls.bvar.variancetightness,'String');
      VarianceTightness = str2num(StringTrim(VarianceTightnessStr(DSGEModel.VarianceTightnessValue,:)));
      Amatrix = VarianceTightness*eye(p);
   end;
   %
   % we always draw from an inverse Wishart distribution
   %
   OmegaDraws = zeros(p*(p+1)/2,5000);
   num = 0;
   num = 0;
   for i=1:5000;
      OmegaDraws(:,i) = vech(InvWishartRndFcn(Amatrix,qDF));
      num = num+1;
      if num==100;
         %
         % make sure that the dialog is refreshed
         %
         drawnow;
         num = 0;
      end;
   end;
   %
   % compute the kernel density estimates
   %
   for i=1:p*(p+1)/2;
%      [XValue,YValue] = KernelDensityPlot(OmegaDraws(i,:),'kepan');
      [XValue,YValue] = YADAMHDensEst(OmegaDraws(i,:),2^9,0,'gaussian',[]);
%      [YValue,XValue] = YADAksdensity(OmegaDraws(i,:));
      if i==1;
         XValues = zeros(p*(p+1)/2,length(XValue));
         YValues = XValues;
      end;
      XValues(i,:) = XValue';
      YValues(i,:) = YValue';
   end;
end;
%
% check if we should just exit
%
if nargout==2;
   Xdata = XValues;
   Ydata = YValues;
   return;
end;
%
% determine the number of parameters
%
NumParam = size(XValues,1);
%
% Determine the variable names
%
VariableNamesBVAR = DSGEModel.VariableNames(DSGEModel.BVARY,:);
if isempty(DSGEModel.BVARX)==0;
   XVariableNamesBVAR = DSGEModel.XVariableNames(DSGEModel.BVARX,:);
else;
   XVariableNamesBVAR = '';
end;
%
% 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 = '';

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',['Prior Densities of ' ParamStr]);
   else;
      set(FigHandles{i},'Name',['Prior 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,:)) ')'];
               RestrictX = 0;
            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') ')'];
               RestrictX = 0;
               ParamTestValue = 1.0e-03;
            else;
               if IndexOne==p;
                  IndexOne = IndexTwo+1;
                  IndexTwo = IndexTwo+1;
               else;
                  IndexOne = IndexOne+1;
               end;
               FigParmStr = ['Omega(' StringTrim(VariableNamesBVAR(IndexOne,:)) ',' StringTrim(VariableNamesBVAR(IndexTwo,:)) ')'];
               if IndexOne~=IndexTwo;
                  RestrictX = 0;
               else;
                  RestrictX = 1;
               end;
               ParamTestValue = 1.0e-06;
            end;
            %
            % plot the posterior density
            %
            PriorHandle = plot(XValues(NumPlots,:),YValues(NumPlots,:),'-','Color',[0.75 0 0],'LineWidth',2);
            hold('on');
            SetAxesFonts(CurrINI,'Title',FigParmStr);
            axis('tight');
            x_lim = get(gca,'XLim');
            %
            % check if we should restrict the XLim in some way
            %
            if (strcmp(lower(ParamStr),'omega')==1)|((strcmp(lower(ParamStr),'pi')==1)&(DSGEModel.PriorKernel==1)&(DSGEModel.PriorType==2));
               maxY = max(YValues(NumPlots,:));
               Ycheck = (YValues(NumPlots,:)/maxY>ParamTestValue) .* (1:length(YValues(NumPlots,:)));
               Ycheck = Ycheck(Ycheck>0);
               x_lim = [XValues(NumPlots,min(Ycheck)) XValues(NumPlots,max(Ycheck))];
            end;
            if (RestrictX==1)&(x_lim(1)<0);
               x_lim(1) = 0;
            end;
            ylim = get(gca,'YLim');
            add_on = 0.1*(ylim(2)-ylim(1));
            ylim = [0 ylim(2)+add_on];
            %
            set(gca,'XLim',x_lim,'YLim',ylim);
            hold('off');
         end;
      end;
   end;
   %
   % save it to disk before
   %
   if DirExist(GraphDirectory)==1;
      if NumFigs>1;
         file = [GraphDirectory '\PriorDensity-BVAR-' ParamStr '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(i,'%0.0f') '.fig'];
      else;
         file = [GraphDirectory '\PriorDensity-BVAR-' ParamStr '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%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 GraphBVARPriorDensity.m
%
