function [DoneCalc,PredEventData,ShockCov,ParameterCov] = BVARPredictionPaths(PsiPostSample,PiPostSample,OmegaPostSample,DSGEModel,CurrINI,DPred,h,NumPaths,AnnualizeData)
% BVARPredictionPaths: Calculates a certain number of possible paths that the endogenous variables
%                      can take out-of-sample for a given sample of draws from the posterior of the
%                      parameters in the Bayesian VAR model.
%
% USAGE:
%
%       [DoneCalc,PredEventData,ShockCov,ParameterCov] = BVARPredictionPaths(PsiPostSample,PiPostSample,OmegaPostSample,DSGEModel,...
%                                                           CurrINI,DPred,h,NumPaths,AnnualizeData)
%
% REQUIRED INPUT:  PsiPostSample (matrix) of dimension NumDraws x pq with posterior draws of the steady state parameters.
%
%                  PiPostSample (matrix) of dimension NumDraws x p^2 k with posterior draws of the parameters on lagged
%                               endogenous variables.
%
%                  OmegaPostSample (matrix) of dimension NumDraws x p(p+1)/2 with posterior draws of the covariance matrix
%                               of the VAR residuals.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure= with initialization information.
%
%                  DPred (matrix) of dimension qxh with data for the exogenous variables over the prediction sample.
%
%                  h (integer) with the length of the prediction sample.
%
%                  NumPaths (integer) with the number of prediction paths to compute
%
%                  AnnualizeData (boolean) that is 1 if the prediction paths should be annualized and 0 otherwise.
%
% REQUIRED OUTPUT: DoneCalc (boolean) that is 1 if the calculations were completed and 0 otherwise.
%
%                  PredEventData (n x 3) matrix where the first column holds the number of times
%                               the PredPaths data satisfies the event, the second hold the number
%                               of times the PredPaths data is greater than the event, while the
%                               third column holds the number of times the event has been tested.
%
%                  ShockCov (p x p x h) matrix with the estimated shock uncertainty.
%
%                  ParameterCov (p x p x h) matrix with the estimated parameter uncertainty.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: May 10, 2007.
%                        This version: December 19, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:
%
% * 14-5-2007:  Fixed a bug when constructing the "Xinit" matrix. Moreover, the vectors "DSGEModel.annual" and
%               "DSGEModel.annualscale" didn't properly take the dimension of the BVAR model into account.
%
% * 15-5-2007:  Added an extra "drawnow" after shutting down the progress and the wait dialog.
%
% * 6-6-2007:   Changed "predictive paths" to "prediction paths" in the title for the progress dialog.
%
% * 6-7-2007:   Added an output-directory creation check.
%
% * 25-10-2007: The prediction event calculations now take the "AdjustPredictionPaths" boolean variable into
%               account.
%
% * 13-11-2007: Updated the documentation.
%
% * 2-1-2008:   Removed code using "DSGEModel.InitializeKalman".
%
% * 10-1-2008:  Changed the last drawnow call to pause(0.02).
%
% * 23-5-2008:  Updated the documentation.
%
% * 19-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize
%
DoneCalc = 0;
PredEventData = [];
NumDraws = size(PiPostSample,1);
%
% set up the data matrices
%
p = length(DSGEModel.BVARY);
k = DSGEModel.BVARLags;
q = length(DSGEModel.BVARX);
%
% create the output directory
%
stat = MakeDir(DSGEModel.OutputDirectory,'predictions');
if stat~=1;
   ShockCov = [];
   ParameterCov = [];
   txt = ['YADA was for some reason unable to create the directory "' DSGEModel.OutputDirectory '\predictions". The computation of predictive distributions has therefore been aborted.'];
   About(txt,'information','YADA - Directory Creation Problem',180,500,CurrINI);
   return;
end;
%
% 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 the initial value for Y at period T
%
T = LastPeriod-(FirstPeriod+KalmanFirstObservation-1)+1;
X_T = zeros(p*k,1);
for i=1:k;
   X_T(1+(p*(i-1)):p*i,1) = DSGEModel.Y(DSGEModel.BVARY,LastPeriod-i+1);
end;
%
% and now for the exogenous variables
%
if isempty(DSGEModel.BVARX)==0;
   d = DPred;
   dLags = zeros(q*k,1);
   for i=1:k;
      dLags(1+(q*(i-1)):q*i,:) = DSGEModel.X(DSGEModel.BVARX,LastPeriod-i+1);
   end;
else;
   d = [];
end;
%
% setup the J-matrix
%
if k>1;
   J = [eye(p) zeros(p,p*(k-1))];
else;
   J = eye(p);
end;
%
% setup an initialization data for the endogenous variables
%
Xinit = [];
if AnnualizeData==1;
   if strcmp(DSGEModel.DataFrequency,'q')==1;
      Xinit = DSGEModel.Y(DSGEModel.BVARY,LastPeriod-2:LastPeriod);
   elseif strcmp(DSGEModel.DataFrequency,'m')==1;
      Xinit = DSGEModel.Y(DSGEModel.BVARY,LastPeriod-10:LastPeriod);
   end;
   ShockCov = [];
   ParameterCov = [];
   TypeStr = 'Annualized';
else;
   ShockCov = zeros(p,p,h);
   ParameterCov = ShockCov;
   MeanPredictY = zeros(p,h);
   TypeStr = 'Original';
end;
if DSGEModel.RunPredictionEvent==1;
   PredEventData = zeros(p,3);
end;
%
% compute the p^2 x p(p+1)/2 duplication matrix
%
Dp = DuplicationMatrix(p);
%
% check if we should setup a wait handle
%
if DSGEModel.ShowProgress==1;
   ProgressStructure.title = ['Progress for ' num2str(NumDraws*NumPaths,'%0.0f') ' Prediction Paths'];
   ProgressStructure.facecolor = CurrINI.progress_facecolor;
   ProgressStructure.startfacecolor = CurrINI.progress_startfacecolor;
   ProgressStructure.edgecolor = CurrINI.progress_edgecolor;
   ProgressStructure.bgcolor = CurrINI.progress_bgcolor;
   ProgressStructure.stop = 0;
   ProgressStructure.clock = DSGEModel.ShowProgressClock;
   ProgressStructure.label = 'Path Number:';
   %
   ProgressStructure.name = 'BVAR Predictive Distributions';
   ProgressStructure.CurrINI = CurrINI;
   WaitHandle = ProgressDLG(0,ProgressStructure);
   set(WaitHandle,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
else;
   txt = ['Please wait while YADA computes predictive distributions for with the Bayesian VAR for the ' lower(TypeStr) ' data. Computations started at: ' StringTrim(datestr(now,14))];
   WaitHandle = WaitDLG(txt,'information',['BVAR Predictive Distributions - ' TypeStr ' Data'],500,150,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
end;
%
% begin the loop
%
abort = '';
StopEstimation = 0;
ParamIndex = 1;
for ParamIndex=1:NumDraws;
   %
   % setup the parameter matrices
   %
   if isempty(DSGEModel.BVARX)==0;
      Psi = reshape(PsiPostSample(ParamIndex,:)',p,q);
      Y_T = X_T-(kron(eye(k),Psi)*dLags);
   else;
      Y_T = X_T;
   end;
   Pi = reshape(PiPostSample(ParamIndex,:)',p,p*k);
   Omega = reshape(Dp*(OmegaPostSample(ParamIndex,:)'),p,p);
   OmegaStMat = chol(Omega)';
   %
   % rewrite the VAR on VAR form
   %
   if k>1;
      B = [Pi;eye(p*(k-1)) zeros(p*(k-1),p)];
   else;
      B = Pi;
   end;
   %
   % compute the shock and parameter uncertainty
   %
   if AnnualizeData==0;
      %
      % initialize output
      %
      PredictY = zeros(p,h);
      Bi = eye(p*k);
      BjSum = zeros(p,p);
      for i=1:h;
         Bj = J*Bi*J';
         BjSum = BjSum+(Bj*Omega*Bj');
         ShockCov(:,:,i) = ShockCov(:,:,i)+((1/NumDraws)*BjSum);
         %
         % update Fi
         %
         Bi = B*Bi;
         if isempty(DSGEModel.BVARX)==0;
            PredictY(:,i) = (Psi*d(:,i))+(J*Bi*Y_T);
         else;
            PredictY(:,i) = J*Bi*Y_T;
         end;
         ParameterCov(:,:,i) = ParameterCov(:,:,i)+((1/NumDraws)*PredictY(:,i)*PredictY(:,i)');
      end;
      MeanPredictY = MeanPredictY+((1/NumDraws)*PredictY);
   end;
   %
   if DSGEModel.AdjustPredictionPaths==1;
      %
      % compute the mean predictions at the given parameter value
      %
      PredMean = zeros(p,h);
      %
      % initialize the Y-vector
      %
      Y = Y_T;
      for i=1:h;
         %
         % update the Y-vector
         %
         Y = (B*Y);
         %
         % calculate the PredMean
         %
         if isempty(d)==0;
            PredMean(:,i) = (Psi*d(:,i))+(J*Y);
         else;
            PredMean(:,i) = J*Y;
         end;
      end;
      if AnnualizeData==1;
         PredMean = Annualize([Xinit PredMean],DSGEModel.annual(DSGEModel.BVARY),DSGEModel.annualscale(DSGEModel.BVARY));
      end;
   end;
   %
   % we are now ready to compute possible paths at the selected parameter value
   %
   PredPaths = zeros(p,h,NumPaths);
   for it=1:NumPaths;
      CurrPath = it+(NumPaths*(ParamIndex-1));
      %
      % check if we should update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(WaitHandle,'UserData');
         if strcmp(abort,'cancel')==1;
            StopEstimation = 1;
            break;
         else;
            ProgressDLG([CurrPath/(NumDraws*NumPaths) CurrPath]);
         end;
      end;
      %
      % initialize the Y-vector
      %
      Y = Y_T;
      %
      % draw the residuals
      %
      EpsilonMat = OmegaStMat*randn(p,h);
      X = zeros(p,h);
      %
      % create the path
      %
      i=1;
      for i=1:h;
         %
         % update the Y-vector
         %
         Y = (B*Y)+(J'*EpsilonMat(:,i));
         %
         % calculate the X-vector
         %
         if isempty(d)==0;
            X(:,i) = (Psi*d(:,i))+(J*Y);
         else;
            X(:,i) = J*Y;
         end;
      end;
      if AnnualizeData==0;
         PredPaths(:,:,it) = X;
      else;
         PredPaths(:,:,it) = Annualize([Xinit X],DSGEModel.annual(DSGEModel.BVARY),DSGEModel.annualscale(DSGEModel.BVARY));
      end;
      drawnow;
   end;
   if StopEstimation==1;
      break;
   end;
   %
   if DSGEModel.AdjustPredictionPaths==1;
      %
      % we need to fix the paths such that the sample mean is equal to the population mean
      %
      MeanPrediction = zeros(p,h);
      for i=1:NumPaths;
         MeanPrediction = MeanPrediction+((1/NumPaths)*PredPaths(:,:,i));
      end;
      AdjustMean = PredMean-MeanPrediction;
      for i=1:NumPaths;
         PredPaths(:,:,i) = PredPaths(:,:,i)+AdjustMean;
      end;
   end;
   %
   % save the calculations to disk
   %
   if AnnualizeData==0;
      file = [DSGEModel.OutputDirectory '\predictions\BVAR-Prediction-Original-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(h,'%0.0f') '-' num2str(ParamIndex,'%0.0f') '.' num2str(NumDraws,'%0.0f') '.' num2str(NumPaths,'%0.0f') '.mat'];;
   else;
      file = [DSGEModel.OutputDirectory '\predictions\BVAR-Prediction-Annual-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(h,'%0.0f') '-' num2str(ParamIndex,'%0.0f') '.' num2str(NumDraws,'%0.0f') '.' num2str(NumPaths,'%0.0f') '.mat'];;
   end;
   save(file,'PredPaths');
   if DSGEModel.RunPredictionEvent==1;
      %
      % compute the prediction events results
      %
      PredEventData = PredEventData+CalculatePredictionEvents(PredPaths,DSGEModel.PredictionEvent(DSGEModel.BVARY,:));
   end;
   drawnow;
end;
%
% close the wait dialog
%
if DSGEModel.ShowProgress==1;
   if ishandle(WaitHandle)==1;
      set(WaitHandle,'UserData','done');
      close(WaitHandle);
      drawnow;
      pause(0.02);
   end;
else;
   set(WaitControls.text,'UserData','done');
   delete(WaitHandle);
   drawnow;
   pause(0.02)
end;
if StopEstimation==1;
   return;
end;
%
% make sure the parameter uncertainty is corrected for the mean
%
DoneCalc = 1;
if isfield(DSGEModel,'PredictionEvent')==1;
   PredictionEvent = DSGEModel.PredictionEvent(DSGEModel.BVARY,:);
else;
   PredictionEvent = [Inf*ones(p,1) zeros(p,1) ones(p,1)];
end;
if AnnualizeData==0;
   for i=1:h;
      ParameterCov(:,:,i) = ParameterCov(:,:,i)-(MeanPredictY(:,i)*MeanPredictY(:,i)');
   end;
   save(file,'PredEventData','PredictionEvent','ShockCov','ParameterCov','-APPEND');
else;
   save(file,'PredEventData','PredictionEvent','-APPEND');
end;

%
% end of BVARPredictionPaths.m
%
