function [DoneCalc,PredEventData,StateCov,ShockCov,MeasureCov,ParameterCov] = DSGEPredictionPaths(thetaPostSample,thetaMode,thetaPositions,ModelParameters,DSGEModel,CurrINI,X,h,NumPaths,AnnualizeData,TransData,CurrChain)
% DSGEPredictionPaths: 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 DSGE model.
%
% USAGE:
%
%       [DoneCalc,PredEventData,StateCov,ShockCov,MeasureCov,ParameterCov] = DSGEPredictionPaths(thetaPostSample,thetaMode,thetaPositions,ModelParameters,...
%                                                                                  DSGEModel,CurrINI,X,h,NumPaths,AnnualizeData,TransData,CurrChain)
%
% REQUIRED INPUT:  thetaPostSample (matrix) of dimension NumDrawsxNumParam with draws from the posterior distribution
%                                  of the original DSGE model parameters.
%
%                  thetaMode (vector) with the posterior mode value of the original DSGE model parameters.
%
%                  thetaPositions (structure) with fields parameter and distribution. This function
%                                  needs the former field to know which parameter name corresponds to
%                                  position i in theta.
%
%                  ModelParameters (structure) whose fields have names given by the parameter names in the
%                                  DSGE model.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure), with initialization information.
%
%                  X (matrix) of dimension kxMaxPH with the observations of the exogenous variables for
%                                  the prediction paths.
%
%                  h (integer) with the prediction horizon.
%
%                  NumPaths (integer) the number of paths to draw.
%
%                  AnnualizeData (boolean) that is 1 if the prediction paths should be annualized and 0 otherwise.
%
%                  TransData (boolean) that is 1 if the prediction paths should be transformed and 0 otherwise.
%
%                  CurrChain (integer) with the number for the current Markov chain.
%
% 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.
%
%                  StateCov (n x n x h) matrix with the estimated state uncertainty.
%
%                  ShockCov (n x n x h) matrix with the estimated shock uncertainty.
%
%                  MeasureCov (n x n) matrix the the estimated measurement error uncertainty.
%
%                  ParameterCov (n x n 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-2009 European Central Bank.
%
%                       First version: February 22, 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:
%
% * 23-02-2007: Finished the function.
%
% * 26-02-2007: The names of the inputs "PredictedX" and "NumPredPaths" were changed to "X" and "NumPaths".
%
% * 27-02-2007: Simplified the check if user cancelled part.
%
% * 28-02-2007: Added a wait dialog when the progress dialog is not displayed.
%
% * 09-03-2007: Updated the documentation.
%
% * 04-05-2007: Added the prediction event calculations.
%
% * 15-05-2007: Added an extra "drawnow" after shutting down the progress and the wait dialog.
%
% * 06-06-2007: Changed "predictive paths" to "prediction paths" in the title for the progress dialog.
%
% * 25-10-2007: The code now takes the boolean variable "AdjustPredictionPaths" into account. Hence, for each
%               selected parameter vector, the 3D matrix with prediction paths is adjusted such that its
%               sample mean equals its population mean.
%
% * 13-11-2007: Updated the documentation.
%
% * 03-12-2007: Made sure the Kalman filter function includes the correct X matrix.
%
% * 18-12-2007: Added the "h" input to allow for X being empty.
%
% * 21-12-2007: Added the code for dealing with unit-root state variables when running the Kalman filter.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 23-05-2008: Updated the documentation.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 09-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 15-07-2008: Fixed an initialization bug for wTh.
%
% * 22-01-2009: Updated the documentation.
%
% * 13-02-2009: The mcode values 7 and 8 are now properly mentioned when encountered.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize the output
%
[NumDraws NumParam] = size(thetaPostSample);
k = DSGEModel.k;
DoneCalc = 0;
PredEventData = [];
if (AnnualizeData==1)|(TransData==1);
   StateCov = [];
   ShockCov = [];
   MeasureCov = [];
   ParameterCov = [];
   if TransData==0;
      TypeStr = 'Annualized';
   else;
      TypeStr = 'Transformed';
   end;
else;
   StateCov = zeros(DSGEModel.n,DSGEModel.n,h);
   ShockCov = StateCov;
   MeasureCov = zeros(DSGEModel.n,DSGEModel.n);
   ParameterCov = StateCov;
   MeanPredictY = zeros(DSGEModel.n,h);
   TypeStr = 'Original';
end;
%
% copy files to the tmp folder
%
[stat,msg] = dos(['copy /Y "' GetPath(DSGEModel.AIMDataFile) 'compute_aim_matrices.m" "' pwd '\tmp"']);
[stat,msg] = dos(['copy /Y "' DSGEModel.MeasurementEquationFile '" "' pwd '\tmp"']);
if FileExist(DSGEModel.UpdateParameterFile)==1;
   [stat,msg] = dos(['copy /Y "' DSGEModel.UpdateParameterFile '" "' pwd '\tmp"']);
end;
%
% load the AIMData file
%
ErrorStr = '';
try;
   AIMData = load(DSGEModel.AIMDataFile);
catch;
   ErrorStr = ['Unable to load the file "' DSGEModel.AIMDataFile '". Message caught is: ' lasterr];
end;
if isempty(ErrorStr)==0;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% create the output directory
%
stat = MakeDir(DSGEModel.OutputDirectory,'predictions');
if stat~=1;
   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;
%
% set the InitialStateVector value
%
if DSGEModel.UseOwnInitialState==1;
   if length(DSGEModel.InitialStateValues)==size(AIMData.endog,1);
      DSGEModel.InitialStateVector = DSGEModel.InitialStateValues(DSGEModel.StateVariablePositions);
   else;
      DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
   end;
else;
   DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
end;
%
% Get the sub-sample
%
[DSGEModel.FirstPeriod,DSGEModel.LastPeriod] = CreateSubSample(DSGEModel);
%
% try to solve the model for the posterior mode estimates
%
ErrorStr = '';
try;
   [A,H,R,F,B0,mcode] = SolveDSGEModel(DSGEModel,thetaMode,thetaPositions,ModelParameters,AIMData);
catch;
   ErrorStr = ['YADA caught an error when trying to solve the DSGE model. Message caught is:' lasterr];
end;
if isempty(ErrorStr)==0;
   About(ErrorStr,'information','YADA - Error Message',200,500,CurrINI);
   return;
end;
if mcode~=1;
   if mcode==2;
      mcodeStr = 'Roots are not correctly computed by real_schur.'
   elseif mcode==3;
      mcodeStr = 'Too many big roots.';
   elseif mcode==35;
      mcodeStr = 'Too many big roots, and q(:,right) is singular.';
   elseif mcode==4;
      mcodeStr = 'Too few big roots.';
   elseif mcode==45;
      mcodeStr = 'Too few big roots, and q(:,right) is singular.';
   elseif mcode==5;
      mcodeStr = 'q(:,right) is singular.';
   elseif mcode==61;
      mcodeStr = 'Too many exact shiftrights.';
   elseif mcode==62;
      mcodeStr = 'Too many numeric shiftrights.';
   elseif mcode==7;
      mcodeStr = 'Infinite or NaN values detected.';
   elseif mcode==8;
      mcodeStr = 'The function "compute_aim_matrices" returns complex numbers.';
   else;
      mcodeStr = 'Return code not properly specified.';
   end;
   txt = ['The AiM solver provided the return code: ' num2str(mcode) ', i.e., "' mcodeStr '"'];
   About(txt,'information','YADA - AiM Solver Error',200,500,CurrINI);
   return;
end;
if length(size(H))==3;
   if size(H,3)<DSGEModel.LastPeriod+h;
      %
      % too few data points in the 3D H matrix.
      %
      ErrorStr = ['The number of time periods for the time-varying measurement matrix H (' int2str(size(H,3)) ') is lower than the number of observations plus forecasts (T+h = ' int2str(DSGEModel.LastPeriod+h) '). YADA has therefore aborted from the unconditional predictions.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
% set up the data matrices
%
Y = DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod);
if isempty(DSGEModel.X)==0;
   Xdata = DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod);
else;
   Xdata = [];
end;
if strcmp(DSGEModel.DataFrequency,'q')==1;
   Yinit = DSGEModel.Y(:,DSGEModel.LastPeriod-2:DSGEModel.LastPeriod);
elseif strcmp(DSGEModel.DataFrequency,'m')==1;
   Yinit = DSGEModel.Y(:,DSGEModel.LastPeriod-10:DSGEModel.LastPeriod);
else;
   Yinit = [];
end;
if DSGEModel.RunPredictionEvent==1;
   PredEventData = zeros(DSGEModel.n,3);
end;
%
if TransData==1;
   Yinitial = DSGEModel.Y(:,1:DSGEModel.LastPeriod);
   %
   % try to covert: h
   %
   YTmp = [Yinitial DSGEModel.Y(:,DSGEModel.LastPeriod)*ones(1,h)];
   [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
   Tt = size(YData,1);
end;
%
% 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 = '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 ' lower(TypeStr) ' data. Computations started at: ' StringTrim(datestr(now,14))];
   WaitHandle = WaitDLG(txt,'information',['Predictive Distributions - ' TypeStr ' Data'],500,150,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
end;
%
% begin the loop
%
abort = '';
StopEstimation = 0;
if length(size(H))==2;
   for ParamIndex=1:NumDraws;
      %
      % Solve the model for the current parameter vector
      %
      theta = thetaPostSample(ParamIndex,:)';
      [A,H,R,F,B0] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
      %
      % We can now run the Kalman filter
      %
      if isempty(DSGEModel.UnitRootStates)==1;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilter(Y,Xdata,A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      else;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilter(Y,Xdata,A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      end;
      %
      % Now we need to calculate KsiTT and PTT for just the last period
      %
      [r,T] = size(Ksitt1);
      q = size(B0,2);
      inv_msey = inv((H'*Ptt1(:,:,T)*H)+R);
      KsiTT = Ksitt1(:,T)+(Ptt1(:,:,T)*H*inv_msey*(DSGEModel.Y(:,DSGEModel.LastPeriod)-Yhat(:,T)));
      PTT = Ptt1(:,:,T)-(Ptt1(:,:,T)*H*inv_msey*H'*Ptt1(:,:,T));
      %
      % use the Schur decomposition for the reduced rank matrix PTT instead of
      % the Choleski decomposition
      %
      [S,L] = SchurDecomposition(PTT);
      %
      % check if we have measurement errors and if so compute a Schur decomposition of the R matrix
      %
      if rank(R)>0;
         [Rm,Lm] = SchurDecomposition(R);
         Rmeasure = Rm*sqrt(Lm);
         NumMeasure = size(Lm,1);
      else;
         Rmeasure = [];
         wTh = 0;
         NumMeasure = 0;
      end;
      if (AnnualizeData==0)&(TransData==0);
         MeasureCov = MeasureCov+((1/NumDraws)*R);
         %
         % initialize output
         %
         PredictY = zeros(DSGEModel.n,h);
         Fi = eye(r);
         FjSum = zeros(r,r);
         Q = B0*B0';
         for i=1:h;
            FjSum = FjSum+(Fi*Q*Fi');
            ShockCov(:,:,i) = ShockCov(:,:,i)+((1/NumDraws)*H'*FjSum*H);
            %
            % update Fi
            %
            Fi = F*Fi;
            StateCov(:,:,i) = StateCov(:,:,i)+((1/NumDraws)*H'*Fi*PTT*Fi'*H);
            if k>0;
               PredictY(:,i) = (A'*X(:,i))+(H'*Fi*KsiTT);
            else;
               PredictY(:,i) = H'*Fi*KsiTT;
            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 prediction at the given parameter values
         %
         PredMean = zeros(DSGEModel.n,h);
         KsiTm = KsiTT;
         for i=1:h;
            %
            % update the states
            %
            KsiTm = F*KsiTm;
            %
            % compute the conditional mean predictions
            %
            if k>0;
               PredMean(:,i) = (A'*X(:,i))+(H'*KsiTm);
            else;
               PredMean(:,i) = H'*KsiTm;
            end;
         end;
         %
         % annualize the mean predictions when needed
         %
         if AnnualizeData==1;
            PredMean = Annualize([Yinit PredMean],DSGEModel.annual,DSGEModel.annualscale);
         elseif TransData==1;
            YTmp = [Yinitial PredMean];
            [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
            PredMean = YData(Tt-h+1:Tt,:)';
         end;
      end;
      %
      % we are now ready to compute possible paths at theta
      %
      PredPaths = zeros(DSGEModel.n,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;
         %
         % Draw a state vector for time T using the N(KsiTT,PTT) distribution.
         %
         KsiT = KsiTT+(S*sqrt(L)*randn(size(L,2),1));
         %
         % Draw a sequence of state shocks
         %
         vTh = B0*randn(q,h);
         %
         % calculate a path for the states
         %
         KsiTh = zeros(r,h);
         KsiTh(:,1) = (F*KsiT)+vTh(:,1);
         for i=2:h;
            KsiTh(:,i) = (F*KsiTh(:,i-1))+vTh(:,i);
         end;
         %
         % Draw a sequence of measurement errors
         %
         if isempty(Rmeasure)==0;
            wTh = Rmeasure*randn(NumMeasure,h);
         end;
         %
         % add everything
         %
         if AnnualizeData==0;
            if k>0;
               PredPaths(:,:,it) = (A'*X)+(H'*KsiTh)+wTh;
            else;
               PredPaths(:,:,it) = (H'*KsiTh)+wTh;
            end;
            if TransData==1;
               YTmp = [Yinitial PredPaths(:,:,it)];
               [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
               PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
            end;
         else;
            if k>0;
               PredPaths(:,:,it) = Annualize([Yinit (A'*X)+(H'*KsiTh)+wTh],DSGEModel.annual,DSGEModel.annualscale);
            else;
               PredPaths(:,:,it) = Annualize([Yinit (H'*KsiTh)+wTh],DSGEModel.annual,DSGEModel.annualscale);
            end;
         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(DSGEModel.n,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;
         if TransData==0;
            file = [DSGEModel.OutputDirectory '\predictions\DSGE-Prediction-Original-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' num2str(h,'%0.0f') '-' num2str(ParamIndex,'%0.0f') '.' num2str(NumDraws,'%0.0f') '.' num2str(NumPaths,'%0.0f') '.mat'];;
         else;
            file = [DSGEModel.OutputDirectory '\predictions\DSGE-Prediction-Transformed-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' num2str(h,'%0.0f') '-' num2str(ParamIndex,'%0.0f') '.' num2str(NumDraws,'%0.0f') '.' num2str(NumPaths,'%0.0f') '.mat'];;
         end;
      else;
         file = [DSGEModel.OutputDirectory '\predictions\DSGE-Prediction-Annual-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%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);
      end;
   end;
else;
   for ParamIndex=1:NumDraws;
      %
      % Solve the model for the current parameter vector
      %
      theta = thetaPostSample(ParamIndex,:)';
      [A,H,R,F,B0] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
      %
      % We can now run the Kalman filter
      %
      if isempty(DSGEModel.UnitRootStates)==1;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(Y,Xdata,A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      else;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilterHt(Y,Xdata,A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      end;
      %
      % Now we need to calculate KsiTT and PTT for just the last period
      %
      [r,T] = size(Ksitt1);
      q = size(B0,2);
      inv_msey = inv((H(:,:,DSGEModel.LastPeriod)'*Ptt1(:,:,T)*H(:,:,DSGEModel.LastPeriod))+R);
      KsiTT = Ksitt1(:,T)+(Ptt1(:,:,T)*H(:,:,DSGEModel.LastPeriod)*inv_msey*(DSGEModel.Y(:,DSGEModel.LastPeriod)-Yhat(:,T)));
      PTT = Ptt1(:,:,T)-(Ptt1(:,:,T)*H(:,:,DSGEModel.LastPeriod)*inv_msey*H(:,:,DSGEModel.LastPeriod)'*Ptt1(:,:,T));
      %
      % use the Schur decomposition for the reduced rank matrix PTT instead of
      % the Choleski decomposition
      %
      [S,L] = SchurDecomposition(PTT);
      %
      % check if we have measurement errors and if so compute a Schur decomposition of the R matrix
      %
      if rank(R)>0;
         [Rm,Lm] = SchurDecomposition(R);
         Rmeasure = Rm*sqrt(Lm);
         NumMeasure = size(Lm,1);
      else;
         Rmeasure = [];
         wTh = zeros(DSGEModel.n,h);
         NumMeasure = 0;
      end;
      if (AnnualizeData==0)&(TransData==0);
         MeasureCov = MeasureCov+((1/NumDraws)*R);
         %
         % initialize output
         %
         PredictY = zeros(DSGEModel.n,h);
         Fi = eye(r);
         FjSum = zeros(r,r);
         Q = B0*B0';
         for i=1:h;
            FjSum = FjSum+(Fi*Q*Fi');
            ShockCov(:,:,i) = ShockCov(:,:,i)+((1/NumDraws)*H(:,:,DSGEModel.LastPeriod+i)'*FjSum*H(:,:,DSGEModel.LastPeriod+i));
            %
            % update Fi
            %
            Fi = F*Fi;
            StateCov(:,:,i) = StateCov(:,:,i)+((1/NumDraws)*H(:,:,DSGEModel.LastPeriod+i)'*Fi*PTT*Fi'*H(:,:,DSGEModel.LastPeriod+i));
            if k>0;
               PredictY(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*Fi*KsiTT);
            else;
               PredictY(:,i) = H(:,:,DSGEModel.LastPeriod+i)'*Fi*KsiTT;
            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 prediction at the given parameter values
         %
         PredMean = zeros(DSGEModel.n,h);
         KsiTm = KsiTT;
         for i=1:h;
            %
            % update the states
            %
            KsiTm = F*KsiTm;
            %
            % compute the conditional mean predictions
            %
            if k>0;
               PredMean(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTm);
            else;
               PredMean(:,i) = H(:,:,DSGEModel.LastPeriod+i)'*KsiTm;
            end;
         end;
         %
         % annualize the mean predictions when needed
         %
         if AnnualizeData==1;
            PredMean = Annualize([Yinit PredMean],DSGEModel.annual,DSGEModel.annualscale);
         elseif TransData==1;
            YTmp = [Yinitial PredMean];
            [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
            PredMean = YData(Tt-h+1:Tt,:)';
         end;
      end;
      %
      % we are now ready to compute possible paths at theta
      %
      PredPaths = zeros(DSGEModel.n,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;
         %
         % Draw a state vector for time T using the N(KsiTT,PTT) distribution.
         %
         KsiT = KsiTT+(S*sqrt(L)*randn(size(L,2),1));
         %
         % Draw a sequence of state shocks
         %
         vTh = B0*randn(q,h);
         %
         % calculate a path for the states
         %
         KsiTh = zeros(r,h);
         KsiTh(:,1) = (F*KsiT)+vTh(:,1);
         for i=2:h;
            KsiTh(:,i) = (F*KsiTh(:,i-1))+vTh(:,i);
         end;
         %
         % Draw a sequence of measurement errors
         %
         if isempty(Rmeasure)==0;
            wTh = Rmeasure*randn(NumMeasure,h);
         end;
         %
         % add everything
         %
         if AnnualizeData==0;
            if k>0;
               for i=1:h;
                  PredPaths(:,i,it) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            else;
               for i=1:h;
                  PredPaths(:,i,it) = (H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            end;
            if TransData==1;
               YTmp = [Yinitial PredPaths(:,:,it)];
               [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
               PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
            end;
         else;
            YTmp = zeros(DSGEModel.n,h);
            if k>0;
               for i=1:h;
                  YTmp(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            else;
               for i=1:h;
                  YTmp(:,i) = (H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            end;
            PredPaths(:,:,it) = Annualize([Yinit YTmp],DSGEModel.annual,DSGEModel.annualscale);
         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(DSGEModel.n,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;
         if TransData==0;
            file = [DSGEModel.OutputDirectory '\predictions\DSGE-Prediction-Original-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' num2str(h,'%0.0f') '-' num2str(ParamIndex,'%0.0f') '.' num2str(NumDraws,'%0.0f') '.' num2str(NumPaths,'%0.0f') '.mat'];;
         else;
            file = [DSGEModel.OutputDirectory '\predictions\DSGE-Prediction-Transformed-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' num2str(h,'%0.0f') '-' num2str(ParamIndex,'%0.0f') '.' num2str(NumDraws,'%0.0f') '.' num2str(NumPaths,'%0.0f') '.mat'];;
         end;
      else;
         file = [DSGEModel.OutputDirectory '\predictions\DSGE-Prediction-Annual-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%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);
      end;
   end;
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;
else;
   PredictionEvent = [Inf*ones(DSGEModel.n,1) zeros(DSGEModel.n,1) ones(DSGEModel.n,1)];
end;
if (AnnualizeData==0)&(TransData==0);
   for i=1:h;
      ParameterCov(:,:,i) = ParameterCov(:,:,i)-(MeanPredictY(:,i)*MeanPredictY(:,i)');
   end;
   save(file,'PredEventData','PredictionEvent','StateCov','ShockCov','MeasureCov','ParameterCov','-APPEND');
else;
   save(file,'PredEventData','PredictionEvent','-APPEND');
end;

%
% end of DSGEPredictionPaths.m
%
