function [Xdata,Ydata,BoundsMatrix] = GraphPriorDensity(theta,thetaDist,LowerBound,PriorDist,ParameterNames,controls,DSGEModel,CurrINI)
% GraphPriorDensity: Creates plots of the prior densities for all original parameters. When the two output
%                    arguments are listed, then the function simply resturns matrices with x-axis and y-axis
%                    values for the densities.
%
% USAGE:
%
%       [Xdata,Ydata,BoundsMatrix] = GraphPriorDensity(theta,thetaDist,LowerBound,PriorDist,ParameterNames,controls,DSGEModel,CurrINI)
%
% REQUIRED INPUT:  theta (vector), with initial values of the parameters.
%
%                  thetaDist (vector), with integers 1 if gamma, 2 if beta, 3 if normal, 4 if invgamma, 5 if
%                            truncnormal, and 6 if uniform.
%
%                  LowerBound (vector) with lower bound values.
%
%                  PriorDist (structure), with fields "beta", "gamma", "normal", "invgamma", "truncnormal", and "uniform" where
%                            the values are matrices with the parameter values for the distributions. In the case of
%                            the normal, beta, and gamma these parameter values are the mean and the standard deviation.
%                            For the invgamma it is the s and q parameters, where q is a positive integer (degrees of
%                            freedom). Finally, for the left truncated normal there's a third parameter, the lower bound.
%
%                  ParameterNames (structure), with fields "calibrated", "beta", "gamma", "normal", "invgamma", "truncnormal",
%                            "uniform", "all" and "estimated". Each field returns a string matrix with the parameter names.
%
%                  controls (structure), whose fields provide handles to all the UI controls in YADA's main dialog.
%
%                  DSGEModel (structure) whose fields contain all required information about the model,
%                            such as input files, variable and sample data.
%
%                  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.
%
%                  BoundsMatrix (matrix) with dimension NumParamx2. It contains information about possible
%                            upper and lower bounds of the parameters.
%
%
%                       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: August 28, 2006.
%                        This version: January 29, 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:
%
% * 29-8-2006: Added the call to the SetToolbarImages function. Kernel density estimates added, but they are not yet
%              enabled. The random number generators appear to work well through! But this method is slower and also
%              tends to include too much of the tails.
%
% * 31-8-2006: Made use of the initialization options "PriorKernel".
%
% * 1-9-2006:  Added a CloseRequestFcn for the figure.
%
% * 22-9-2006: Fixed a bug with "scrsz" variable being called rather than "CurrINI.scrsz".
%
% * 7-12-2006:  Took the change in CurrINI.NumberOfGridPoints into DSGEModel.NumberOfGridPoints into account. This also means
%               that the function requires the input DSGEModel. Moreover, the CurrINI.GridWidth is now DSGEModel.GridWidth,
%               while CurrINI.PriorKernel is DSGEModel.PriorKernel and CurrINI.KernelDensityEstimator is
%               DSGEModel.KernelDensityEstimator.
%
% * 15-12-2006: Added the optional output.
%
% * 6-2-2007:   Updated the documentation.
%
% * 8-2-2007:   Added a second "drawnow" call at the end of the figure loop.
%
% * 9-11-2007:  Updated the documentation.
%
% * 10-1-2008:  Changed the last drawnow call to pause(0.02).
%
% * 23-5-2008:  Updated the documentation.
%
% * 9-7-2008:   Changed the setup up the grid for the left truncated normal density.
%
% * 11-7-2008:  Fixed a bug when setting up the grid for the truncated normal density.
%
% * 21-7-2008:  Added the "LowerBound" input variable. The gamma and inverse gamma priors now take
%               these bounds into account.
%
% * 22-7-2008:  Made sure that the grids for the gamma and inverse gamma take the lower bound into
%               account. Also, the the lower bound is subtracted from the "s" parameter of the
%               the inverse gamma.
%
% * 23-7-2008:  Took into account that the beta distribution now has general lower and upper bounds
%               rather then the 0-1 bounds for the standardized beta distribution.
%
% * 19-12-2008: Updated the documentation.
%
% * 29-01-2009: Fixed a bug when calling "UniformRndFcn" using too few input variables and from the
%               wrong distribution.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% First calculate the values to plot
%
NumGridPointsStr = get(controls.optimize.numbergridpoints,'String');
NumGridPoints = str2num(StringTrim(NumGridPointsStr(DSGEModel.NumberOfGridPoints,:)));
NumGridPointsHalf = NumGridPoints/2;
if DSGEModel.PriorKernel==0;
   thetaValues = zeros(length(theta),NumGridPoints+1);
   DensityValue = thetaValues;
else;
   KernelDensityX = zeros(length(theta),512);
   KernelDensityY = KernelDensityX;
end;
%
n_beta = 0;
n_gamma = 0;
n_normal = 0;
n_invgamma = 0;
n_truncnormal = 0;
n_uniform = 0;
DistNames = '';
BoundMatrix = zeros(length(theta),2);
%
% lower bounds for the gamma and inverse gamma distributions
%
for i=1:length(theta);
   %
   % Determine the distribution
   %
   if thetaDist(i)==1;
      %
      % gamma dist
      %
      n_gamma = n_gamma+1;
      %
      % determine the grid
      %
      p_std = PriorDist.gamma(n_gamma,2);
      lb = PriorDist.gamma(n_gamma,1)-(DSGEModel.GridWidth*p_std);
      dev = (PriorDist.gamma(n_gamma,1)-lb)/NumGridPointsHalf;
      if lb<=LowerBound(i);
         lb = 1.0e-004+LowerBound(i);
      end;
      if DSGEModel.PriorKernel==0;
         thetaVals = (lb:dev:lb+(NumGridPoints*dev));
         logDensityVals = logGammaPDF(thetaVals'-LowerBound(i),PriorDist.gamma_ab(n_gamma,1),PriorDist.gamma_ab(n_gamma,2));
      else;
         z = LowerBound(i)+GammaRndFcn(PriorDist.gamma_ab(n_gamma,1),PriorDist.gamma_ab(n_gamma,2),5000);
      end;
      BoundMatrix(i,:) = [LowerBound(i) Inf];
      DistNames = strvcat(DistNames,'gamma');
   elseif thetaDist(i)==2;
      %
      % beta dist
      %
      n_beta = n_beta+1;
      if size(PriorDist.beta,2)==4;
         %
         % get the lower and upper bounds (if available from PriorDist.)
         %
         c = PriorDist.beta(n_beta,3);
         d = PriorDist.beta(n_beta,4);
      else;
         c = 0;
         d = 1;
      end;
      p_std = PriorDist.beta(n_beta,2);
      lb = PriorDist.beta(n_beta,1)-(DSGEModel.GridWidth*p_std);
      dev = (PriorDist.beta(n_beta,1)-lb)/NumGridPointsHalf;
      if lb<=c;
         lb = c+1.0e-004;
      end;
      if DSGEModel.PriorKernel==0;
         thetaVals = (lb:dev:lb+(NumGridPoints*dev));
         thetaVals(thetaVals>=1) = d-1.0e-004;
         logDensityVals = logBetaPDF(thetaVals',PriorDist.beta_ab(n_beta,1),PriorDist.beta_ab(n_beta,2),c,d);
      else;
         z = BetaRndFcn(PriorDist.beta_ab(n_beta,1),PriorDist.beta_ab(n_beta,2),c,d,5000);
      end;
      BoundMatrix(i,:) = [c d];
      DistNames = strvcat(DistNames,'beta');
   elseif thetaDist(i)==3;
      %
      % normal
      %
      n_normal = n_normal+1;
      p_std = PriorDist.normal(n_normal,2);
      lb = PriorDist.normal(n_normal,1)-(DSGEModel.GridWidth*p_std);
      dev = (PriorDist.normal(n_normal,1)-lb)/NumGridPointsHalf;
      thetaVals = (lb:dev:lb+(NumGridPoints*dev));
      logDensityVals = logNormalPDF(thetaVals',PriorDist.normal(n_normal,1),PriorDist.normal(n_normal,2));
      if DSGEModel.PriorKernel==1;
         z = NormalRndFcn(PriorDist.normal(n_normal,1),PriorDist.normal(n_normal,2),5000);
      end;
      BoundMatrix(i,:) = [-Inf Inf];
      DistNames = strvcat(DistNames,'normal');
   elseif thetaDist(i)==4;
      %
      % inverse gamma
      %
      n_invgamma = n_invgamma+1;
      p_par1 = PriorDist.invgamma(n_invgamma,1);
      p_par2 = PriorDist.invgamma(n_invgamma,2);
      if p_par2>2;
         p_mean = gammaln((p_par2-1)/2)-gammaln(p_par2/2);
         p_mean = exp(p_mean)*sqrt(p_par2/2)*p_par1;
         p_std = sqrt(((p_par2*p_par1*p_par1)/(p_par2-2))-(p_mean^2));
      else;
         p_mean = gammaln(1)-gammaln(3/2);
         p_mean = exp(p_mean)*sqrt(3/2)*p_par1;
         p_std = sqrt(((3*p_par1*p_par1))-(p_mean^2));
      end;
      lb = p_par1-(DSGEModel.GridWidth*p_std);
      dev = (p_par1-lb)/NumGridPointsHalf;
      if lb<=LowerBound(i);
         lb = 1.0e-004+LowerBound(i);
      end;
      if DSGEModel.PriorKernel==0;
         thetaVals = (lb:dev:lb+(NumGridPoints*dev));
         logDensityVals = logInverseGammaPDF(thetaVals'-LowerBound(i),p_par1-LowerBound(i),p_par2);
      else;
         z = LowerBound(i)+InvGammaRndFcn(p_par1-LowerBound(i),p_par2,5000);
      end;
      BoundMatrix(i,:) = [LowerBound(i) Inf];
      DistNames = strvcat(DistNames,'inv gamma');
   elseif thetaDist(i)==5;
      %
      % left truncated normal
      %
      n_truncnormal = n_truncnormal+1;
      p_std = PriorDist.truncnormal(n_truncnormal,2);
      lb = theta(i)-(DSGEModel.GridWidth*p_std);
      dev = abs(PriorDist.truncnormal(n_truncnormal,1)-lb)/NumGridPointsHalf;
%      if lb<=PriorDist.truncnormal(n_truncnormal,3);
%         lb = PriorDist.truncnormal(n_truncnormal,3)+1.0e-008;
%      end;
      if DSGEModel.PriorKernel==0;
%         thetaVals = (lb:dev:lb+(NumGridPoints*dev));
         thetaVals = (max(lb,PriorDist.truncnormal(n_truncnormal,3)+1.0e-008):dev:max(lb,PriorDist.truncnormal(n_truncnormal,3)+1.0e-008)+(NumGridPoints*dev));
%         if length(thetaVals)>size(thetaValues,2);
%            thetaVals = thetaVals(1:size(thetaValues,2));
%         elseif length(thetaVals)<size(thetaValues,2);
%            thetaVals = [thetaVals (max(thetaVals)+dev:dev:max(thetaVals)+(dev*(size(thetaValues,2)-length(thetaVals))))];
%         end;
         logDensityVals = logLTNormalPDF(thetaVals',PriorDist.truncnormal(n_truncnormal,1),PriorDist.truncnormal(n_truncnormal,2),PriorDist.truncnormal(n_truncnormal,3));
      else;
         z = LTNormalRndFcn(PriorDist.truncnormal(n_truncnormal,1),PriorDist.truncnormal(n_truncnormal,2),PriorDist.truncnormal(n_truncnormal,3),5000);
      end;
      BoundMatrix(i,:) = [PriorDist.truncnormal(n_truncnormal,3) Inf];
      DistNames = strvcat(DistNames,'trunc normal');
   elseif thetaDist(i)==6;
      %
      % uniform
      %
      n_uniform = n_uniform+1;
      p_par1 = PriorDist.uniform(n_uniform,1);
      p_par2 = PriorDist.uniform(n_uniform,2);
      p_std = sqrt((p_par2-p_par1)/12);
      mid_point = (p_par2+p_par1)/2;
      lb = mid_point-(DSGEModel.GridWidth*p_std);
      dev = (mid_point-lb)/NumGridPointsHalf;
      if lb<p_par1;
         lb = p_par1;
      end;
      if DSGEModel.PriorKernel==0;
         thetaVals = (lb:dev:lb+(NumGridPoints*dev));
         thetaVals(thetaVals>p_par2) = p_par2;
         logDensityVals = logUniformPDF(p_par_1,p_par2)*ones(length(thetaVals),1);
      else;
         z = UniformRndFcn(PriorDist.uniform(n_uniform,1),PriorDist.uniform(n_uniform,2),5000);
      end;
      BoundMatrix(i,:) = [p_par_1 p_par2];
      DistNames = strvcat(DistNames,'uniform');
   end;
   if DSGEModel.PriorKernel==0;
      thetaValues(i,:) = thetaVals;
      DensityValue(i,:) = exp(logDensityVals)';
   else;
      %
      % Compute the Kernel density estimate of the density function
      %
      [XValues,YValues] = KernelDensityPlot(z,DSGEModel.KernelDensityEstimator);
      %
      % make sure the height is equal for the two estimates when using the normal distribution
      %
      if thetaDist(i)==3;
         YValues = (max(exp(logDensityVals))/max(YValues))*YValues;
      end;
      KernelDensityX(i,:) = XValues';
      KernelDensityY(i,:) = YValues';
   end;
end;
%
% check if we should just return Xdata and Ydata
%
if nargout==3;
   if DSGEModel.PriorKernel==0;
      Xdata = thetaValues;
      Ydata = DensityValue;
   else;
      Xdata = KernelDensityX;
      Ydata = KernelDensityY;
   end;
   BoundsMatrix = BoundMatrix;
   return;
end;
%
% Compute the number of figures (max 16 plots per graph)
%
NumFigs = ceil(length(theta)/16);
%
if NumFigs==1;
   [NumRows,NumCols] = VillaniConstructOptimalSubplot(length(theta));
else;
   NumRows = 4;
   NumCols = 4;
end;
NumPlots = 0;
pos = [(CurrINI.scrsz(3)-1000)/2 (CurrINI.scrsz(4)-680)/2 1000 650];
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 theta');
   else;
      set(FigHandles{i},'Name',['Prior Densities of theta - Figure ' num2str(i)]);
   end;
   SetFigureProperties(FigHandles{i},CurrINI);
   k = 0;
   for r=1:NumRows;
      for c=1:NumCols;
         NumPlots = NumPlots+1;
         if NumPlots<=length(theta);
            k = k+1;
            subplot(NumRows,NumCols,k);
            %
            % plot the density
            %
            if DSGEModel.PriorKernel==0;
               dh = plot(thetaValues(NumPlots,:),DensityValue(NumPlots,:),'-','Color',[0.75 0 0],'LineWidth',2);
            else;
               dh = plot(KernelDensityX(NumPlots,:),KernelDensityY(NumPlots,:),'-','Color',[0.75 0 0],'LineWidth',2);
            end;
            hold('on');
            SetAxesFonts(CurrINI,'Title',[StringTrim(ParameterNames.estimated(NumPlots,:)) ' - ' StringTrim(DistNames(NumPlots,:)) ' prior']);
            %
            % plot the initial value
            %
            y_lim = get(gca,'Ylim');
            y_lim(1) = 0;
            ivh = plot(theta(NumPlots)*ones(1,2),y_lim,'-','Color',[0 0 0]);
            kids = [dh ivh];
            set(gca,'Children',kids);
            axis('tight');
            x_lim = get(gca,'XLim');
            if DSGEModel.PriorKernel==1;
               if strcmp(StringTrim(DistNames(NumPlots,:)),'gamma')==1;
                  x_lim(1) = BoundMatrix(NumPlots,1);
               elseif strcmp(StringTrim(DistNames(NumPlots,:)),'beta')==1;
                  x_lim = BoundMatrix(NumPlots,:);
               elseif strcmp(StringTrim(DistNames(NumPlots,:)),'inv gamma')==1;
                  x_lim(1) = BoundMatrix(NumPlots,1);
               elseif strcmp(StringTrim(DistNames(NumPlots,:)),'trunc normal')==1;
                  x_lim(1) = BoundMatrix(NumPlots,1);
               elseif strcmp(StringTrim(DistNames(NumPlots,:)),'uniform')==1;
                  x_lim = BoundMatrix(NumPlots,:);
               end;
            end;
            set(gca,'YLim',y_lim,'XLim',x_lim);
            hold('off');
         end;
      end;
   end;
   if DSGEModel.PriorKernel==0;
      legHandle = legend([dh ivh],strvcat('Grid Density Estimate','Initial Value'));
   else;
      legHandle = legend([dh ivh],strvcat('Kernel Density Estimate','Initial Value'));
   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);
   %
   % 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;
   pause(0.02);
end;

%
% end of GraphPriorDensity.m
%
