function [DoneCalc,PredEventData,ShockMean,ShockNames,StateCov,ShockCov,MeasureCov,ParameterCov] = DSGECondPredictionPaths(thetaPostSample,thetaMode,thetaPositions,ModelParameters,DSGEModel,CurrINI,X,Z,U,h,NumPaths,AnnualizeData,TransData,CurrChain)
% DSGECondPredictionPaths: 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 and a given path for some conditioning variables.
%
% USAGE:
%
%       [DoneCalc,PredEventData,ShockMean,ShockNames,StateCov,ShockCov,MeasureCov,ParameterCov] = DSGECondPredictionPaths(thetaPostSample,thetaMode,thetaPositions,ModelParameters,...
%                                                                                      DSGEModel,CurrINI,X,Z,U,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) of dimension NumParam with the posterior mode estimate 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 kxh with the observations of the exogenous variables for
%                                  the prediction paths.
%
%                  Z (matrix) of dimension mxg with the out-of-sample data we shall condition the
%                                  predictions on.
%
%                  U (matrix) of dimension mxg with initial values for the mapping between the observed
%                                  data and the conditioning data.
%
%                  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.
%
%                  ShockMean (q x h) matrix with the average of the shocks under conditioning when
%                                  the measurement error is set to zero and the state vector at T
%                                  is set to KsiTT.
%
%                  ShockNames (q x qstr) string matrix with the names of the economic shocks.
%
%                  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: March 9, 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:
%
% * 12-03-2007: Finished the function.
%
% * 14-03-2007: Made sure that the K1 and K2 matrices took the field "ConditionalVariables" of the structure "DSGEModel"
%               into account.
%
% * 26-03-2007: Updated the documentation. The modesty statistics now use KsiT = KsiTT so that the state is treated as
%               known for the conditional and unconditional forecasts for these statistics.
%
% * 27-03-2007: Made sure the Leeper-Zha statistic uses the correct state vector.
%
% * 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.
%
% * 13-07-2007: Fixed a bug when computing "MeanPredictY" where shocks for i=g+1,...,h were not taklen into
%               account.
%
% * 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.
%
% * 17-12-2007: Added the input "h". This makes it possible for X to be an empty matrix.
%
% * 18-12-2007: Added the outputs "ShockMean" for the mean value of the shocks under conditioning, and "ShockNames"
%               with the names for the shocks.
%
% * 20-12-2007: Fixed an annualization bug when the theoretical mean is used for correcting the sampled
%               mean. The non-annualized mean was erroneously used.
%
% * 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).
%
% * 24-01-2008: The input U is now a matrix. Smooth estimation of the state and state covariance using
%               the conditioning assumptions has been included.
%
% * 05-05-2008: Added the input variable "TransData".
%
% * 23-05-2008: Updated the documentation.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 07-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 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;
g = size(Z,2);
DoneCalc = 0;
PredEventData = [];
ShockMean = [];
ShockNames = DSGEModel.StateShockNames;
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 conditional 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
   %
   YTmp = [Yinitial DSGEModel.Y(:,DSGEModel.LastPeriod)*ones(1,h)];
   [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
   Tt = size(YData,1);
end;
%
q = length(DSGEModel.StateShockPositions);
%
% create the matrices for selecting shocks
%
N = eye(q);
M = N(:,DSGEModel.ConditionalShocks);
M_orth = N(:,setdiff((1:q),DSGEModel.ConditionalShocks));
m = size(Z,1);
ShockMean = zeros(q,h);
%
% create the mapping matrices for conditioning variables
%
K1 = DSGEModel.K1(:,DSGEModel.ConditionalVariables);
K2 = DSGEModel.K2(:,DSGEModel.ConditionalVariables);
%
% test if we can invert the model for the standardized shocks at posterior mode
%
if length(size(H))==2;
   etaMat = K1'*H'*B0*M;
   if rank(etaMat)<m;
      txt = ['The shocks you have selected for conditioning cannot be used due to a singularity in the solution algorithm. Estimation of the conditional predictive distribution has therfore been aborted.'];
      About(txt,'error','Error - Conditional Shocks',150,500,CurrINI);
      PredPaths = [];
      return;
   end;
   etaMat = etaMat\eye(m);
else;
   BadEtaMat = 0;
   for i=1:g;
      etaMat = K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*M;
      if rank(etaMat)<m;
         txt = ['The shocks you have selected for conditioning cannot be used due to a singularity in the solution algorithm. Estimation of the conditional predictive distribution has therfore been aborted.'];
         About(txt,'error','Error - Conditional Shocks',150,500,CurrINI);
         PredPaths = [];
         BadEtaMat = 1;
         break;
      end;
   end;
   if BadEtaMat==1;
      return;
   end;
end;
%
% check if we should setup a wait handle
%
if DSGEModel.ShowProgress==1;
   ProgressStructure.title = ['Progress for ' num2str(NumDraws*NumPaths,'%0.0f') ' Conditional 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 = 'Conditional Predictive Distributions';
   ProgressStructure.CurrINI = CurrINI;
   WaitHandle = ProgressDLG(0,ProgressStructure);
   set(WaitHandle,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
else;
   txt = ['Please wait while YADA computes conditional predictive distributions for ' lower(TypeStr) ' data. Computations started at: ' StringTrim(datestr(now,14))];
   WaitHandle = WaitDLG(txt,'information',['Conditional Predictive Distributions - ' TypeStr ' Data'],500,150,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
end;
%
% begin the loop
%
abort = '';
StopEstimation = 0;
ParamIndex = 1;
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);
      %
      % check if we can use these values.
      %
      etaMat = K1'*H'*B0*M;
      if rank(etaMat)==m;
         etaMat = etaMat\eye(m);
         %
         % 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 smooth values for the states for g periods
         %
         if ParamIndex==1;
            [r,T] = size(Ksitt1);
            q = size(B0,2);
         end;
         [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmoother(Y(:,T-g+1:T),Yhat(:,T-g+1:T),Ksitt1(:,T-g+1:T),Ptt1(:,:,T-g+1:T),H,F,R);
         %
         % calculate measurement errors for the sample T-g+1 until T
         %
         if DSGEModel.KsiUseCondData==1;
            if rank(R)>0;
               if k>0;
                  wtT = Y(:,T-g+1:T)-(A'*Xdata(:,T-g+1:T))-(H'*KsitT);
               else;
                  wtT = Y(:,T-g+1:T)-(H'*KsitT);
               end;
            else;
               wtT = zeros(DSGEModel.n,g);
            end;
            %
            % Smooth estimation of the states at time T using the conditioning assumptions.
            %
            if k>0;
               if g>1;
                  Xtra = Xdata(:,T-g+2:T);
               else;
                  Xtra = [];
               end;
               [KsiTTz,PTTz] = CondPredictionKalmanSmoother(KsitT,PtT,wtT,Z,U,[Xtra X(:,1:g)],K1,K2,A,H,R,F,B0);
            else;
               [KsiTTz,PTTz] = CondPredictionKalmanSmoother(KsitT,PtT,wtT,Z,U,[],K1,K2,A,H,R,F,B0);
            end;
         else;
            KsiTTz = KsitT(:,g);
            PTTz = PtT(:,:,g);
         end;
         %
         KsiTT = KsitT(:,g);
         %
         % use the Schur decomposition for the reduced rank matrix PTT instead of
         % the Choleski decomposition
         %
         [S,L] = SchurDecomposition(PTTz);
         %
         % 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;
            if TransData==0;
               MeasureCov = MeasureCov+((1/NumDraws)*R);
            end;
            %
            % initialize output
            %
            PredictY = zeros(DSGEModel.n,h);
            PredictYC = PredictY;
            Fi = eye(r);
            FjSum = zeros(r,r);
            Q = B0*B0';
            QM_orth = B0*M_orth*M_orth'*B0';
            QM = B0*M*M'*B0';
            %
            % initialize the state uncertainty for modesty statistics
            %
            PTTg = zeros(r,r);
            PTTgLZ = PTTg;
            for i=1:g;
               FjSum = FjSum+(Fi*QM_orth*Fi');
               %
               % update Fi
               %
               Fi = F*Fi;
               if TransData==0;
                  ShockCov(:,:,i) = ShockCov(:,:,i)+((1/NumDraws)*H'*FjSum*H);
                  StateCov(:,:,i) = StateCov(:,:,i)+((1/NumDraws)*H'*Fi*PTTz*Fi'*H);
               end;
               if k>0;
                  PredictY(:,i) = (A'*X(:,i))+(H'*Fi*KsiTT);
                  PredictYC(:,i) = (A'*X(:,i))+(H'*Fi*KsiTTz);
               else;
                  PredictY(:,i) = H'*Fi*KsiTT;
                  PredictYC(:,i) = H'*Fi*KsiTTz;
               end;
               PTTg = (F*PTTg*F')+Q;
               PTTgLZ = (F*PTTgLZ*F')+QM;
            end;
            for i=g+1:h;
               FjSum = FjSum+(Fi*Q*Fi');
               %
               % update Fi
               %
               Fi = F*Fi;
               if TransData==0;
                  ShockCov(:,:,i) = ShockCov(:,:,i)+((1/NumDraws)*H'*FjSum*H);
                  StateCov(:,:,i) = StateCov(:,:,i)+((1/NumDraws)*H'*Fi*PTTz*Fi'*H);
               end;
            end;
         end;
         %
         % compute the mean prediction at the given parameter values
         %
         PredMean = zeros(DSGEModel.n,h);
         KsiTm = KsiTTz;
         etaShocksm = zeros(q,h);
         K2summ = zeros(m,1);
         for i=1:g;
            %
            % compute the shock values that need to be manipulated
            %
            if k>0;
               etaShocksm(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*H'*F*KsiTm) - K2summ -U(:,g));
            else;
               etaShocksm(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*H'*F*KsiTm) - K2summ -U(:,g));
            end;
            %
            % update the states
            %
            KsiTm = (F*KsiTm)+(B0*etaShocksm(:,i));
            %
            % compute the conditional mean predictions
            %
            if k>0;
               PredMean(:,i) = (A'*X(:,i))+(H'*KsiTm);
            else;
               PredMean(:,i) = H'*KsiTm;
            end;
            %
            % update the K2sum measure
            %
            K2summ = K2summ+(K2'*PredMean(:,i));
         end;
         for i=g+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;
         ShockMean = ShockMean+((1/NumDraws)*etaShocksm);
         %
         % we are now ready to compute possible paths at theta
         %
         PredPaths = zeros(DSGEModel.n,h,NumPaths);
         LZPredPaths = zeros(DSGEModel.n,g);
         if AnnualizeData==0;
            MultiModestyStat = zeros(NumPaths,2);
            UniModestyStat = zeros(NumPaths,DSGEModel.n);
            UniModestyStatLZ = zeros(1,DSGEModel.n);
            if TransData==1;
               YTmp = [Yinitial PredMean];
               [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
               PredMean = YData(Tt-h+1:Tt,:)';
            end;
         else;
            PredMean = Annualize([Yinit PredMean],DSGEModel.annual,DSGEModel.annualscale);
            MultiModestyStat = [];
            UniModestyStat = [];
            UniModestyStatLZ = [];
         end;
         %
         % start the loop over the paths
         %
         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;
            if AnnualizeData==0;
               PredictYShocks = zeros(r,h);
            end;
            %
            % Draw a state vector for time T using the N(KsiTTz,PTTz) distribution.
            %
            KsiTinit = KsiTTz+(S*sqrt(L)*randn(size(L,2),1));
            %
            % Draw a sequence of economic shocks and set the conditioning shocks to zero
            %
            etaTh = randn(q,h);
            etaTmp = etaTh(:,1:g);
            etaCTmp = etaTh(:,1:g);
            etaTh(DSGEModel.ConditionalShocks,1:g) = zeros(length(DSGEModel.ConditionalShocks),g);
            etaLZ = zeros(q,h);
            %
            % Draw a sequence of measurement errors
            %
            if isempty(Rmeasure)==0;
               wTh = Rmeasure*randn(NumMeasure,h);
            end;
            %
            % initialize the K2 part
            %
            K2sum = zeros(m,1);
            K2sumLZ = zeros(m,1);
            K2sumCTmp = zeros(m,1);
            KsiT = KsiTinit;
   %         KsiTLZ = KsiTinit;
            KsiTLZ = KsiTT;
            KsiCTmp = KsiTT;
   %         KsiTmp = KsiTinit;
            KsiTmp = KsiTT;
            %
            % start the loop that calculates shocks, states, and, observed variables
            %
            for i=1:g;
               if k>0;
                  etaTh(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*wTh(:,i)) - (K1'*H'*F*KsiT) - (K1'*H'*B0*etaTh(:,i)) - K2sum -U(:,g));
                  etaCTmp(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*wTh(:,i)) - (K1'*H'*F*KsiCTmp) - (K1'*H'*B0*etaTh(:,i)) - K2sumCTmp -U(:,g));
               else;
                  etaTh(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*wTh(:,i)) - (K1'*H'*F*KsiT) - (K1'*H'*B0*etaTh(:,i)) - K2sum -U(:,g));
                  etaCTmp(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*wTh(:,i)) - (K1'*H'*F*KsiCTmp) - (K1'*H'*B0*etaTh(:,i)) - K2sumCTmp -U(:,g));
               end;
               %
               % recompute the state
               %
               KsiT = (F*KsiT)+(B0*etaTh(:,i));
               KsiCTmp = (F*KsiT)+(B0*etaCTmp(:,i));
               KsiTmp = (F*KsiTmp)+(B0*etaTmp(:,i));
               %
               % compute the observed variables
               %
               if k>0;
                  PredPaths(:,i,it) = (A'*X(:,i))+(H'*KsiT)+wTh(:,i);
               else;
                  PredPaths(:,i,it) = (H'*KsiT)+wTh(:,i);
               end;
               %
               % update the sum of K2 times lagged values of the observed variables
               %
               K2sum = K2sum+(K2'*PredPaths(:,i,it));
               if k>0;
                  K2sumCTmp = K2sumCTmp+(K2'*((A'*X(:,i))+(H'*KsiCTmp)+wTh(:,i)));
               else;
                  K2sumCTmp = K2sumCTmp+(K2'*((H'*KsiCTmp)+wTh(:,i)));
               end;
               %
               % deal with Leeper-Zha matters
               %
               if it==1;
                  if k>0;
                     etaLZ(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*H'*F*KsiTLZ) - K2sumLZ -U(:,g));
                  else;
                     etaLZ(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*H'*F*KsiTLZ) - K2sumLZ -U(:,g));
                  end;
                  KsiTLZ = (F*KsiTLZ)+(B0*etaLZ(:,i));
                  if k>0;
                     LZPredPaths(:,i) = (A'*X(:,i))+(H'*KsiTLZ);
                  else;
                     LZPredPaths(:,i) = H'*KsiTLZ;
                  end;
                  K2sumLZ = K2sumLZ+(K2'*LZPredPaths(:,i));
               end;
               %
               % compute the average of the shocks that we shall add to PredictY(:,i)
               %
               if AnnualizeData==0;
                  if i==1;
                     PredictYShocks(:,i) = PredictYShocks(:,i)+(B0*M*etaTh(DSGEModel.ConditionalShocks,i));
                  else;
                     PredictYShocks(:,i) = PredictYShocks(:,i)+(F*PredictYShocks(:,i-1))+((B0*M*etaTh(DSGEModel.ConditionalShocks,i)));
                  end;
                  %
                  % compute modesty statistics
                  %
                  if i==g;
                     if k>0;
                        PredPathCTmp = (A'*X(:,g))+(H'*KsiCTmp)+wTh(:,g);
                     else;
                        PredPathCTmp = (H'*KsiCTmp)+wTh(:,g);
                     end;
                     PredDiff = PredPathCTmp-PredictY(:,g);
                     PredCov = (H'*PTTg*H)+R;
                     MultiModestyStat(it,1) = PredDiff'*inv(PredCov)*PredDiff;
                     UniModestyStat(it,:) = (PredDiff ./ sqrt(diag(PredCov)))';
                     %
                     % comparison-path for multivariate statistic
                     %
                     if k>0;
                        PredPathTmp = (A'*X(:,g))+(H'*KsiTmp)+wTh(:,g);
                     else;
                        PredPathTmp = (H'*KsiTmp)+wTh(:,g);
                     end;
                     PredDiffTmp = PredPathTmp-PredictY(:,g);
                     MultiModestyStat(it,2) = PredDiffTmp'*inv(PredCov)*PredDiffTmp;
                     %
                     % Leeper-Zha
                     %
                     PredDiff = LZPredPaths(:,g)-PredictY(:,g);
                     if it==1;
                        UniModestyStatLZ(1,:) = (PredDiff ./ sqrt(diag(H'*PTTgLZ*H)))';
                     end;
                  end;
               end;
            end;
            %
            % if h>g we need to run the loop
            %
            for i=g+1:h;
               KsiT = (F*KsiT)+(B0*etaTh(:,i));
               if k>0;
                  PredPaths(:,i,it) = (A'*X(:,i))+(H'*KsiT)+wTh(:,i);
               else;
                  PredPaths(:,i,it) = (H'*KsiT)+wTh(:,i);
               end;
               if AnnualizeData==0;
                  PredictYShocks(:,i) = F*PredictYShocks(:,i-1);
               end;
            end;
            %
            if AnnualizeData==0;
               if TransData==1;
                  YTmp = [Yinitial PredPaths(:,:,it)];
                  [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
                  PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
               else;
                  PredictYit = PredictYC+(H'*PredictYShocks);
                  for i=1:h;
                     ParameterCov(:,:,i) = ParameterCov(:,:,i)+((1/(NumDraws*NumPaths))*PredictYit(:,i)*PredictYit(:,i)');
                  end;
                  MeanPredictY = MeanPredictY+((1/(NumDraws*NumPaths))*PredictYit);
               end;
            else;
               %
               % we should annualize the conditional prediction path
               %
               PredPaths(:,:,it) = Annualize([Yinit PredPaths(:,:,it)],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-CondPrediction-Original-Direct-' 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-CondPrediction-Transformed-Direct-' 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-CondPrediction-Annual-Direct-' 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','ShockMean','ShockNames','MultiModestyStat','UniModestyStat','UniModestyStatLZ');
         %
         if DSGEModel.RunPredictionEvent==1;
            %
            % compute the prediction events results
            %
            PredEventData = PredEventData+CalculatePredictionEvents(PredPaths,DSGEModel.PredictionEvent);
         end;
      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);
      %
      % check if we can use these values.
      %
      BadRank = 0;
      for i=1:g;
         etaMat = K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*M;
         if rank(etaMat)<m;
            BadRank = 1;
            break;
         end;
      end;
      if BadRank==0;
         %
         % 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 smooth values for the states for g periods
         %
         if ParamIndex==1;
            [r,T] = size(Ksitt1);
            q = size(B0,2);
         end;
         [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(Y(:,T-g+1:T),Yhat(:,T-g+1:T),Ksitt1(:,T-g+1:T),Ptt1(:,:,T-g+1:T),H(:,:,DSGEModel.LastPeriod-g+1:DSGEModel.LastPeriod),F,R);
         %
         % calculate measurement errors for the sample T-g+1 until T
         %
         if DSGEModel.KsiUseCondData==1;
            wtT = zeros(DSGEModel.n,g);
            if rank(R)>0;
               if k>0;
                  for i=1:g;
                     wtT(:,i) = Y(:,i)-(A'*Xdata(:,i))-(H(:,:,DSGEModel.LastPeriod-g+i)'*KsitT(:,i));
                  end;
               else;
                  for i=1:g;
                     wtT(:,i) = Y(:,i)-(H(:,:,DSGEModel.LastPeriod-g+i)'*KsitT(:,i));
                  end;
               end;
            end;
            %
            % Smooth estimation of the states at time T using the conditioning assumptions.
            %
            if k>0;
               if g>1;
                  Xtra = Xdata(:,T-g+2:T);
               else;
                  Xtra = [];
               end;
               [KsiTTz,PTTz] = CondPredictionKalmanSmootherHt(KsitT,PtT,wtT,Z,U,[Xtra X(:,1:g)],K1,K2,A,H(:,:,DSGEModel.LastPeriod-g+2:DSGEModel.LastPeriod+g),R,F,B0);
            else;
               [KsiTTz,PTTz] = CondPredictionKalmanSmootherHt(KsitT,PtT,wtT,Z,U,[],K1,K2,A,H(:,:,DSGEModel.LastPeriod-g+2:DSGEModel.LastPeriod+g),R,F,B0);
            end;
         else;
            KsiTTz = KsitT(:,g);
            PTTz = PtT(:,:,g);
         end;
         %
         KsiTT = KsitT(:,g);
         %
         % use the Schur decomposition for the reduced rank matrix PTT instead of
         % the Choleski decomposition
         %
         [S,L] = SchurDecomposition(PTTz);
         %
         % 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;
            if TransData==0;
               MeasureCov = MeasureCov+((1/NumDraws)*R);
            end;
            %
            % initialize output
            %
            PredictY = zeros(DSGEModel.n,h);
            PredictYC = PredictY;
            Fi = eye(r);
            FjSum = zeros(r,r);
            Q = B0*B0';
            QM_orth = B0*M_orth*M_orth'*B0';
            QM = B0*M*M'*B0';
            %
            % initialize the state uncertainty for modesty statistics
            %
            PTTg = zeros(r,r);
            PTTgLZ = PTTg;
            for i=1:g;
               FjSum = FjSum+(Fi*QM_orth*Fi');
               %
               % update Fi
               %
               Fi = F*Fi;
               if TransData==0;
                  ShockCov(:,:,i) = ShockCov(:,:,i)+((1/NumDraws)*H(:,:,DSGEModel.LastPeriod+i)'*FjSum*H(:,:,DSGEModel.LastPeriod+i));
                  StateCov(:,:,i) = StateCov(:,:,i)+((1/NumDraws)*H(:,:,DSGEModel.LastPeriod+i)'*Fi*PTTz*Fi'*H(:,:,DSGEModel.LastPeriod+i));
               end;
               if k>0;
                  PredictY(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*Fi*KsiTT);
                  PredictYC(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*Fi*KsiTTz);
               else;
                  PredictY(:,i) = H(:,:,DSGEModel.LastPeriod+i)'*Fi*KsiTT;
                  PredictYC(:,i) = H(:,:,DSGEModel.LastPeriod+i)'*Fi*KsiTTz;
               end;
               PTTg = (F*PTTg*F')+Q;
               PTTgLZ = (F*PTTgLZ*F')+QM;
            end;
            for i=g+1:h;
               FjSum = FjSum+(Fi*Q*Fi');
               %
               % update Fi
               %
               Fi = F*Fi;
               if TransData==0;
                  ShockCov(:,:,i) = ShockCov(:,:,i)+((1/NumDraws)*H(:,:,DSGEModel.LastPeriod+i)'*FjSum*H(:,:,DSGEModel.LastPeriod+i));
                  StateCov(:,:,i) = StateCov(:,:,i)+((1/NumDraws)*H(:,:,DSGEModel.LastPeriod+i)'*Fi*PTTz*Fi'*H(:,:,DSGEModel.LastPeriod+i));
               end;
            end;
         end;
         %
         % compute the mean prediction at the given parameter values
         %
         PredMean = zeros(DSGEModel.n,h);
         KsiTm = KsiTTz;
         etaShocksm = zeros(q,h);
         K2summ = zeros(m,1);
         for i=1:g;
            %
            % compute the shock values that need to be manipulated
            %
            etaMat = (K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*M)\eye(m);
            if k>0;
               etaShocksm(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiTm) - K2summ -U(:,g));
            else;
               etaShocksm(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiTm) - K2summ -U(:,g));
            end;
            %
            % update the states
            %
            KsiTm = (F*KsiTm)+(B0*etaShocksm(:,i));
            %
            % 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;
            %
            % update the K2sum measure
            %
            K2summ = K2summ+(K2'*PredMean(:,i));
         end;
         for i=g+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;
         ShockMean = ShockMean+((1/NumDraws)*etaShocksm);
         %
         % we are now ready to compute possible paths at theta
         %
         PredPaths = zeros(DSGEModel.n,h,NumPaths);
         LZPredPaths = zeros(DSGEModel.n,g);
         if AnnualizeData==0;
            MultiModestyStat = zeros(NumPaths,2);
            UniModestyStat = zeros(NumPaths,DSGEModel.n);
            UniModestyStatLZ = zeros(1,DSGEModel.n);
            if TransData==1;
               YTmp = [Yinitial PredMean];
               [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
               PredMean = YData(Tt-h+1:Tt,:)';
            end;
         else;
            PredMean = Annualize([Yinit PredMean],DSGEModel.annual,DSGEModel.annualscale);
            MultiModestyStat = [];
            UniModestyStat = [];
            UniModestyStatLZ = [];
         end;
         %
         % start the loop over the paths
         %
         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;
            if AnnualizeData==0;
               PredictYShocks = zeros(r,h);
            end;
            %
            % Draw a state vector for time T using the N(KsiTTz,PTTz) distribution.
            %
            KsiTinit = KsiTTz+(S*sqrt(L)*randn(size(L,2),1));
            %
            % Draw a sequence of economic shocks and set the conditioning shocks to zero
            %
            etaTh = randn(q,h);
            etaTmp = etaTh(:,1:g);
            etaCTmp = etaTh(:,1:g);
            etaTh(DSGEModel.ConditionalShocks,1:g) = zeros(length(DSGEModel.ConditionalShocks),g);
            etaLZ = zeros(q,h);
            %
            % Draw a sequence of measurement errors
            %
            if isempty(Rmeasure)==0;
               wTh = Rmeasure*randn(NumMeasure,h);
            end;
            %
            % initialize the K2 part
            %
            K2sum = zeros(m,1);
            K2sumLZ = zeros(m,1);
            K2sumCTmp = zeros(m,1);
            KsiT = KsiTinit;
   %         KsiTLZ = KsiTinit;
            KsiTLZ = KsiTT;
            KsiCTmp = KsiTT;
   %         KsiTmp = KsiTinit;
            KsiTmp = KsiTT;
            %
            % start the loop that calculates shocks, states, and, observed variables
            %
            for i=1:g;
               etaMat = (K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*M)\eye(m);
               if k>0;
                  etaTh(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*wTh(:,i)) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiT) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*etaTh(:,i)) - K2sum -U(:,g));
                  etaCTmp(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*wTh(:,i)) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiCTmp) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*etaTh(:,i)) - K2sumCTmp -U(:,g));
               else;
                  etaTh(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*wTh(:,i)) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiT) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*etaTh(:,i)) - K2sum -U(:,g));
                  etaCTmp(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*wTh(:,i)) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiCTmp) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*B0*etaTh(:,i)) - K2sumCTmp -U(:,g));
               end;
               %
               % recompute the state
               %
               KsiT = (F*KsiT)+(B0*etaTh(:,i));
               KsiCTmp = (F*KsiT)+(B0*etaCTmp(:,i));
               KsiTmp = (F*KsiTmp)+(B0*etaTmp(:,i));
               %
               % compute the observed variables
               %
               if k>0;
                  PredPaths(:,i,it) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiT)+wTh(:,i);
               else;
                  PredPaths(:,i,it) = (H(:,:,DSGEModel.LastPeriod+i)'*KsiT)+wTh(:,i);
               end;
               %
               % update the sum of K2 times lagged values of the observed variables
               %
               K2sum = K2sum+(K2'*PredPaths(:,i,it));
               if k>0;
                  K2sumCTmp = K2sumCTmp+(K2'*((A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiCTmp)+wTh(:,i)));
               else;
                  K2sumCTmp = K2sumCTmp+(K2'*((H(:,:,DSGEModel.LastPeriod+i)'*KsiCTmp)+wTh(:,i)));
               end;
               %
               % deal with Leeper-Zha matters
               %
               if it==1;
                  if k>0;
                     etaLZ(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiTLZ) - K2sumLZ -U(:,g));
                  else;
                     etaLZ(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiTLZ) - K2sumLZ -U(:,g));
                  end;
                  KsiTLZ = (F*KsiTLZ)+(B0*etaLZ(:,i));
                  if k>0;
                     LZPredPaths(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTLZ);
                  else;
                     LZPredPaths(:,i) = H(:,:,DSGEModel.LastPeriod+i)'*KsiTLZ;
                  end;
                  K2sumLZ = K2sumLZ+(K2'*LZPredPaths(:,i));
               end;
               %
               % compute the average of the shocks that we shall add to PredictY(:,i)
               %
               if AnnualizeData==0;
                  if i==1;
                     PredictYShocks(:,i) = PredictYShocks(:,i)+(B0*M*etaTh(DSGEModel.ConditionalShocks,i));
                  else;
                     PredictYShocks(:,i) = PredictYShocks(:,i)+(F*PredictYShocks(:,i-1))+((B0*M*etaTh(DSGEModel.ConditionalShocks,i)));
                  end;
                  %
                  % compute modesty statistics
                  %
                  if i==g;
                     if k>0;
                        PredPathCTmp = (A'*X(:,g))+(H(:,:,DSGEModel.LastPeriod+g)'*KsiCTmp)+wTh(:,g);
                     else;
                        PredPathCTmp = (H(:,:,DSGEModel.LastPeriod+g)'*KsiCTmp)+wTh(:,g);
                     end;
                     PredDiff = PredPathCTmp-PredictY(:,g);
                     PredCov = (H(:,:,DSGEModel.LastPeriod+g)'*PTTg*H(:,:,DSGEModel.LastPeriod+g))+R;
                     MultiModestyStat(it,1) = PredDiff'*inv(PredCov)*PredDiff;
                     UniModestyStat(it,:) = (PredDiff ./ sqrt(diag(PredCov)))';
                     %
                     % comparison-path for multivariate statistic
                     %
                     if k>0;
                        PredPathTmp = (A'*X(:,g))+(H(:,:,DSGEModel.LastPeriod+g)'*KsiTmp)+wTh(:,g);
                     else;
                        PredPathTmp = (H(:,:,DSGEModel.LastPeriod+g)'*KsiTmp)+wTh(:,g);
                     end;
                     PredDiffTmp = PredPathTmp-PredictY(:,g);
                     MultiModestyStat(it,2) = PredDiffTmp'*inv(PredCov)*PredDiffTmp;
                     %
                     % Leeper-Zha
                     %
                     PredDiff = LZPredPaths(:,g)-PredictY(:,g);
                     if it==1;
                        UniModestyStatLZ(1,:) = (PredDiff ./ sqrt(diag(H(:,:,DSGEModel.LastPeriod+g)'*PTTgLZ*H(:,:,DSGEModel.LastPeriod+g))))';
                     end;
                  end;
               end;
            end;
            %
            % if h>g we need to run the loop
            %
            for i=g+1:h;
               KsiT = (F*KsiT)+(B0*etaTh(:,i));
               if k>0;
                  PredPaths(:,i,it) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiT)+wTh(:,i);
               else;
                  PredPaths(:,i,it) = (H(:,:,DSGEModel.LastPeriod+i)'*KsiT)+wTh(:,i);
               end;
               if AnnualizeData==0;
                  PredictYShocks(:,i) = F*PredictYShocks(:,i-1);
               end;
            end;
            %
            if AnnualizeData==0;
               if TransData==1;
                  YTmp = [Yinitial PredPaths(:,:,it)];
                  [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
                  PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
               else;
                  PredictYit = zeros(size(PredictYC));
                  for i=1:h;
                     PredictYit(:,i) = PredictYC(:,i)+(H(:,:,DSGEModel.LastPeriod+i)'*PredictYShocks(:,i));
                  end;
                  for i=1:h;
                     ParameterCov(:,:,i) = ParameterCov(:,:,i)+((1/(NumDraws*NumPaths))*PredictYit(:,i)*PredictYit(:,i)');
                  end;
                  MeanPredictY = MeanPredictY+((1/(NumDraws*NumPaths))*PredictYit);
               end;
            else;
               %
               % we should annualize the conditional prediction path
               %
               PredPaths(:,:,it) = Annualize([Yinit PredPaths(:,:,it)],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-CondPrediction-Original-Direct-' 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-CondPrediction-Transformed-Direct-' 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-CondPrediction-Annual-Direct-' 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','ShockMean','ShockNames','MultiModestyStat','UniModestyStat','UniModestyStatLZ');
         %
         if DSGEModel.RunPredictionEvent==1;
            %
            % compute the prediction events results
            %
            PredEventData = PredEventData+CalculatePredictionEvents(PredPaths,DSGEModel.PredictionEvent);
         end;
      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 DSGECondPredictionPaths.m
%
