function [PredPaths,PredData,PredEventData,YObsEventData,MultiModestyStat,UniModestyStat,UniModestyStatLZ,status,kalmanstatus] = DSGECondPredictionPathsTheta(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,X,Z,U,h,NumPaths,AnnualizeData,TransData,NameStr)
% DSGECondPredictionPathsTheta: Calculates a certain number of possible paths that the endogenous variablescan take out-of-sample
%                               conditional on the certain out-of-sample data for a given value of the parameters in the DSGE model.
%
% USAGE:
%
%       [PredPaths,PredData,PredEventData,YObsEventData,MultiModestyStat,UniModestyStat,UniModestyStatLZ,status,kalmanstatus] = DSGECondPredictionPathsTheta(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,...
%                                                                      X,Z,U,h,NumPaths,AnnualizeData,TransData,NameStr)
%
% REQUIRED INPUT:  theta (vector) with values 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 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.
%
%                  NameStr (string vector) that is sent to the progress dialog name after the string
%                                 "Predictive Distribution".
%
% REQUIRED OUTPUT: PredPaths (n x h x NumPaths) 3D matrix with the NumPaths possible paths for the n
%                                 endogenous variables.
%
%                  PredData (structure) with fields:
%                                 PredMean             - (nxh) matrix with the mean predictions of the observed
%                                                        variables.
%                                 ShockScenarioResults - (nxhxq_used) matrix with shock scenarios for the mean
%                                                        predictions of the observed variables. These are
%                                                        obtained by manipulating the shocks used to
%                                                        compute the mean predictions. q_used<=q sinze some
%                                                        of the q shocks can be zero always by assumption.
%                                 OrigPredMean         - (nxh) matrix with the non-annualized mean predictions
%                                                        of the observed variables.
%                                 MeanPredPath         - (nxh) matrix with the non-annualized average of the
%                                                        predictions paths for the observed variables.
%                                 etaShocks            - (qxh) matrix with the values of the shocks at the mean
%                                                        predictions.
%                                 KsitT                - (rxT) matrix with the smooth estimates of the states
%                                                        over the estimation sample.
%                                 xiStates             - (rxh) matrix with the values of the states at the mean
%                                                        predictions.
%                                 X                    - (kxh) matrix with the deterministic variables for the
%                                                        prediction horizon.
%                                 etaMean              - (qxh) matrix with the mean values for the economic shocks over
%                                                        the forecast period.
%                                 wMean                - (nxh) matrix with the mean values for the measurement errors over
%                                                        the forecast period.
%                                 A                    - (kxn) matrix with the parameters on the exogenous variables in
%                                                        the measurement equation.
%                                 H                    - (rxn) matrix with parameters on the state variabels in the
%                                                        measurement equation.
%                                 R                    - (nxn) covariance matrix for the measurement errors.
%                                 F                    - (rxr) matrix with state transitions parameters for the state
%                                                        equation.
%                                 B0                   - (rxq) matrix with parameters on the economic shocks in the
%                                                        state equation.
%                                 ShockNames           - (qxqstr) string matrix with the names of the economic shocks.
%                                 Shocks               - (qxhxNumPaths) matrix with the economic shocks for the paths.
%                                 MeasError            - (nxhxNumPaths) matrix with the measurement errors for the paths.
%                                 HistState            - (rxNumPaths) matrix with the time T period draw of the state
%                                                        vector for the paths.
%                                 KernelX    - (3D matrix) of dimension n x 2^8 x h with the x-axis values
%                                              from kernel density estimation of the marginal predictive
%                                              densities.
%                                 KernelY    - (3D matrix) of dimension n x 2^8 x h with the y-axis values
%                                              from kernel density estimation of the marginal predictive
%                                              densities.
%
%                  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.
%
%                  YObsEventData (n x 3) matrix where the first column hold the number of times
%                                 the PredPaths (adjusted such that the mean is equal to the realized
%                                 value) satisfies the event, the second holds the number of times
%                                 the adjusted paths data is greater than the event, while the third
%                                 column holds the number of tiems the event has been tested.
%
%                  MultiModestyStat (matrix) with NumPaths x 2 elements. The first column is the the multivariate
%                                 modesty statistic based on the conditioning consistent shocks and the second
%                                 is the value of this statistic for randomly drawn shocks.
%
%                  UniModestyStat (matrix) with NumPaths x n elements. Each column gives the univariate modesty
%                                 for a particular observed variable.
%
%                  UniModestyStatLZ (vector) with the values of the univariate modesty statistics suggested
%                                 by Leeper-Zha (other shocks are zero over examined prediction horizon) for
%                                 all observed variables.
%
% OPTIONAL OUTPUT: status (integer) Indicates whether the solution of the model is unique or not.
%                                 mcode is 1 when unique, 0 if not computed, and various other values when
%                                 there are problems with the solution. Specifically:
%
%                                      mcode = 2:      Roots not correctly computed by real_schur
%                                      mcode = 3:      Too many big roots
%                                      mcode = 35:     Too many big roots, and q(:,right) is singular
%                                      mcode = 4:      Too few big roots
%                                      mcode = 45:     Too few big roots, and q(:,right) is singular
%                                      mcode = 5:      q(:,right) is singular
%                                      mcode = 61:     Too many exact shiftrights
%                                      mcode = 62:     Too many numeric shiftrights
%                                      mcode = 7:      The a matrix has infinite or NaN values. Eigenvalues
%                                                      cannot be calculated.
%                                      mcode = 8:      The function "compute_aim_matrices" returns complex
%                                                      numbers.
%                                      else            Return code not properly specified
%
%                  kalmanstatus (boolean) indicating the covariance stationarity of the state variables. If all the
%                                 eigenvalues of F are inside the unit circle, then kalmanstatus=0 and the log-likelihood function
%                                 is calculated. If not, then status=1 and initP is set to 0.
%
%
%                       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 8, 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:
%
% * 09-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: Added the calculation of the multivariate and univariate ALLV modesty statistics for given parameter values.
%
% * 27-03-2007: Added the calculation of the univariate Leeper-Zha tyep modesty statistics for given parameter values.
%
% * 02-04-2007: Updated the documentation.
%
% * 19-04-2007: Fixed two bugs related to "K2sumCTmp" which (1) wasn't initialized properly and (2) not subtracted from etaCTmp.
%               This led to incorrect calculations of the modesty statistics.
%
% * 03-05-2007: Added the "PredEventData" output.
%
% * 15-05-2007: Added an extra "drawnow" after shutting down the progress dialog.
%
% * 05-06-2007: Made sure that the Leeper-Zha modesty statistic takes possible zero variances into account.
%
% * 06-06-2007: Changed "predictive paths" to "prediction paths" in the title for the progress dialog.
%
% * 20-08-2007: Added the output "PredData".
%
% * 21-08-2007: Added the fields "PredMean", "OrigPredMean", "MeanPredPath", and "etaShocks" to "PredData".
%
% * 28-08-2007: Added the field "ShockScenarioResults" to "PredData".
%
% * 25-10-2007: The prediction event calculations now take the "AdjustPredictionPaths" boolean variable into
%               account.
%
% * 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". The matrix X can now be empty.
%
% * 18-12-2007: Added the field "ShockNames" to "PredData".
%
% * 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).
%
% * 15-01-2008: Started work on smooth estimation of state variables using the conditioning assumptions.
%
% * 23-01-2008: Finished work on smooth estiamtion of the state variables using the conditioning assumptions.
%
% * 24-01-2008: Implemented the smooth estimates of the state variables based on the conditioning assumptions.
%
% * 16-04-2008: The code now checks if the model has a unique convergent solution also at the posterior mode
%               estimates.
%
% * 05-05-2008: Added the "TransData" input variable.
%
% * 16-05-2008: Took the "DSGEModel.AIMTolerance" input for AiMSolver into account.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: The mcode value 7 is now supported.
%
% * 01-07-2008: Added the fields "xiStates" and "KsitT" to PredData.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%               Added the try-catch code for AiMSolver.
%
% * 07-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 09-07-2008: Made sure that the matrix H for the PredData output has the correct dimension when
%               it is time-varying.
%
% * 15-07-2008: Fixed a bug when computing the population predictive mean for time-varying H.
%
% * 15-01-2009: Added the fields "Shocks", "MeasError", and "HistState" to the PredData output variable.
%               The output variable "YObsEventData" was also added along with the code that can compute
%               prediction events based on a predictive distribution whose mean is equal to the observed
%               data.
%
% * 30-01-2009: Made sure that status is only set to mcode when the latter exists.
%
% * 09-02-2009: Changed the input "length(DSGEModel.StateVariablePositions)" to "size(F,1)" when calling the
%               measurement equation function.
%
% * 13-02-2009: The mcode values 7 and 8 are now properly mentioned when encountered.
%
% * 02-03-2009: Added code for kernel density estimation of the marginal predictive distributions.
%
% * 03-03-2009: Made sure that the fields "KernelX" and "KernelY" of PredData are initialized.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize the output
%
PredPaths = [];
PredData = [];
PredEventData = [];
YObsEventData = [];
MultiModestyStat = [];
UniModestyStat = [];
UniModestyStatLZ = [];
%
% set up the ModelParameters structure based on theta and thetaPositions
%
if nargout==7;
   %
   % Make sure that the ModelParameter structure takes the correct values.
   % This is not necessary to do when initial values are given to the function.
   % In that case, more substantial checking is required!
   %
   ModelParameters = ThetaToModelParameters(ModelParameters,theta,thetaPositions);
   %
   % 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;
   %
   % Update any parameteters that need updating
   %
   if FileExist(DSGEModel.UpdateParameterFile)==1;
      eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
   end;
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;
   if nargout>7;
      status = 0;
      kalmanstatus = 0;
   end;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% See if we can solve the model
%
if nargout==7;
   [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
else;
   try;
      [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
   catch;
      if exist('mcode')>0;
         status = mcode;
      else;
         status = 0;
      end;
      ErrorStr = ['There was an error running AiM. YADA received the following message: ' lasterr];
      About(ErrorStr,'error','AiM Error',300,500,CurrINI);
      drawnow;
      return;
   end;
end;
if nargout>7;
   status = mcode;
   if status~=1;
      kalmanstatus = 0;
      return;
   end;
else;
   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;
end;
%
% Convert to the state equation matrices
%
[F,B0] = AiMtoStateSpace(BMatrix,SMatrix,DSGEModel.StateVariablePositions,DSGEModel.StateShockPositions);
%
% check B0 to make sure that the selected shocks for conditioning are valid
%
Btest = B0(:,DSGEModel.ConditionalShocks);
if rank(Btest)<length(DSGEModel.ConditionalShocks);
   txt = ['The shocks you have selected for conditioning have a singular covariance matrix. The number of shocks for conditioning is therefore too small relative to the number of conditioning variables.'];
   About(txt,'error','Error - Conditional Shocks',150,500,CurrINI);
   return;
end;
%
% obtain the measurement equation matrices
%
try;
   eval(['[A,H,R] = ' GetName(DSGEModel.MeasurementEquationFile) '(ModelParameters,AIMData.endog(DSGEModel.StateVariablePositions,:),DSGEModel.VariableNames,DSGEModel.XVariableNames,DSGEModel.n,size(F,1),DSGEModel.k);']);
catch;
   %
   % display an error message
   %
   ErrorStr = ['Error found in "' GetName(DSGEModel.MeasurementEquationFile) '.m". Message caught is: ' lasterr];
   About(ErrorStr,'error','Error - Measurement Equation',150,500,CurrINI);
   return;
end;
%
% Get the sub-sample
%
[DSGEModel.FirstPeriod,DSGEModel.LastPeriod] = CreateSubSample(DSGEModel);
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 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;
%
% We can now run the Kalman filter
%
if length(size(H))==2;
   if isempty(DSGEModel.UnitRootStates)==1;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),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] = KalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      end;
   else;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      else;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],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;
   end;
else;
   if isempty(DSGEModel.UnitRootStates)==1;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),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] = KalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      end;
   else;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),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));
      else;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],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;
   end;
end;
%
% we should adjust everything based on DSGEModel.KalmanFirstObservation
%
if nargout>7;
   kalmanstatus = stat;
   if kalmanstatus~=0;
      return;
   end;
end;
%
% get some parameter values
%
[r,T] = size(Ksitt1);
k = DSGEModel.k;
g = size(Z,2);
q = size(B0,2);
%
% create the mapping matrices for conditioning variables
%
K1 = DSGEModel.K1(:,DSGEModel.ConditionalVariables);
K2 = DSGEModel.K2(:,DSGEModel.ConditionalVariables);
%
% Now we need to calculate KsiTT and PTT for just the last period
%
%
% run the Kalman smoother for g periods
%
if length(size(H))==2;
   [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmoother(DSGEModel.Y(:,DSGEModel.LastPeriod-g+1:DSGEModel.LastPeriod),Yhat(:,T-g+1:T),Ksitt1(:,T-g+1:T),Ptt1(:,:,T-g+1:T),H,F,R);
else;
   [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(DSGEModel.Y(:,DSGEModel.LastPeriod-g+1:DSGEModel.LastPeriod),Yhat(:,T-g+1:T),Ksitt1(:,T-g+1:T),Ptt1(:,:,T-g+1:T),H(:,:,DSGEModel.LastPeriod-g+1:DSGEModel.LastPeriod),F,R);
end;
%
% calculate measurement errors for the sample T-g+1 until T
%
if DSGEModel.KsiUseCondData==1;
   if length(size(H))==2;
      if rank(R)>0;
         if k>0;
            wtT = DSGEModel.Y(:,DSGEModel.LastPeriod-g+1:DSGEModel.LastPeriod)-(A'*DSGEModel.X(:,DSGEModel.LastPeriod-g+1:DSGEModel.LastPeriod))-(H'*KsitT);
         else;
            wtT = DSGEModel.Y(:,DSGEModel.LastPeriod-g+1:DSGEModel.LastPeriod)-(H'*KsitT);
         end;
      else;
         wtT = zeros(DSGEModel.n,g);
      end;
   else;
      wtT = zeros(DSGEModel.n,g);
      if rank(R)>0;
         if k>0;
            for i=1:g;
               wtT(:,i) = DSGEModel.Y(:,DSGEModel.LastPeriod-g+i)-(A'*DSGEModel.X(:,DSGEModel.LastPeriod-g+i))-(H(:,:,DSGEModel.LastPeriod-g+i)'*KsitT(:,i));
            end;
         else;
            for i=1:g;
               wtT(:,i) = DSGEModel.Y(:,DSGEModel.LastPeriod-g+i)-(H(:,:,DSGEModel.LastPeriod-g+i)'*KsitT(:,i));
            end;
         end;
      end;
   end;
   %
   % Smooth estimation of the states at time T using the conditioning assumptions.
   %
   if length(size(H))==2;
      if k>0;
         if g>1;
            Xtra = DSGEModel.X(:,DSGEModel.LastPeriod-g+2:DSGEModel.LastPeriod);
         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;
      if k>0;
         if g>1;
            Xtra = DSGEModel.X(:,DSGEModel.LastPeriod-g+2:DSGEModel.LastPeriod);
         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;
   end;
else;
   KsiTTz = KsitT(:,g);
   PTTz = PtT(:,:,g);
end;
%
% Now we need to determine KsiTT for just the last period
%
KsiTT = KsitT(:,g);
%
% use the Schur decomposition for the reduced rank matrix PTTz instead of
% the Choleski decomposition
%
[S,L] = SchurDecomposition(PTTz);
%
% set up the output
%
k = DSGEModel.k;
g = size(Z,2);
PredPaths = zeros(DSGEModel.n,h,NumPaths);
LZPredPaths = zeros(DSGEModel.n,g);
%
% 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;
%
%
% compute the state and shock covariances at the given parameters
% as well as the predictions of the observed variables.
%
if (DSGEModel.RunPredictionEvent==1)&(size(DSGEModel.Y,2)>=DSGEModel.LastPeriod+h);
   YObserved = DSGEModel.Y(:,DSGEModel.LastPeriod+1:DSGEModel.LastPeriod+h);
else;
   YObserved = [];
end;
if AnnualizeData==1;
   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;
else;
   MultiModestyStat = zeros(NumPaths,2);
   UniModestyStat = zeros(NumPaths,DSGEModel.n);
   UniModestyStatLZ = zeros(1,DSGEModel.n);
   if TransData==1;
      Yinit = DSGEModel.Y(:,1:DSGEModel.LastPeriod);
   end;
end;
%
% create the matrices for selecting shocks
%
N = eye(q);
M = N(:,DSGEModel.ConditionalShocks);
m = size(Z,1);
%
% check if shocks can be used
%
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' lower(NameStr) ' 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' lower(NameStr) ' 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;
%
% add info to the PredData structure
%
PredData.KernelX = [];
PredData.KernelY = [];
PredData.X = X;
PredData.A = A;
if length(size(H))==2;
   PredData.H = H;
else;
   PredData.H = H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod+h);
end;
PredData.R = R;
PredData.B0 = B0;
PredData.F = F;
PredData.etaMean = zeros(q,h);
PredData.wMean = zeros(DSGEModel.n,h);
PredData.ShockNames = DSGEModel.StateShockNames;
PredData.Shocks = zeros(q,h,NumPaths);
PredData.MeasError = zeros(DSGEModel.n,h,NumPaths);
PredData.HistState = zeros(r,NumPaths);
%
% compute the mean prediction
%
PredData.PredMean = zeros(DSGEModel.n,h);
PredData.MeanPredPath = zeros(DSGEModel.n,h);
KsiT = KsiTTz;
PredData.etaShocks = zeros(q,h);
PredData.xiStates = zeros(r,h);
PredData.KsitT = KsitT;
K2sum = zeros(m,1);
if length(size(H))==2;
   for i=1:g;
      %
      % compute the shock values that need to be manipulated
      %
      if k>0;
         PredData.etaShocks(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*H'*F*KsiT) - K2sum -U(:,g));
      else;
         PredData.etaShocks(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*H'*F*KsiT) - K2sum -U(:,g));
      end;
      %
      % update the states
      %
      KsiT = (F*KsiT)+(B0*PredData.etaShocks(:,i));
      PredData.xiStates(:,i) = KsiT;
      %
      % compute the conditional mean predictions
      %
      if k>0;
         PredData.PredMean(:,i) = (A'*X(:,i))+(H'*KsiT);
      else;
         PredData.PredMean(:,i) = H'*KsiT;
      end;
      %
      % update the K2sum measure
      %
      K2sum = K2sum+(K2'*PredData.PredMean(:,i));
   end;
   for i=g+1:h;
      %
      % update the states
      %
      KsiT = (F*KsiT);
      PredData.xiStates(:,i) = KsiT;
      %
      % compute the conditional mean predictions
      %
      if k>0;
         PredData.PredMean(:,i) = (A'*X(:,i))+(H'*KsiT);
      else;
         PredData.PredMean(:,i) = H'*KsiT;
      end;
   end;
   [B0H,KeepShocks] = RemoveRedundantColumns(B0);
   q_used = size(B0H,2);
   PredData.ShockScenarioResults = zeros(DSGEModel.n,h,q_used);
   for j=1:q_used;
      %
      % set shock j of the shocks manipulated for the theoretical mean of the conditional
      % forecast
      %
      etaShocks = PredData.etaShocks(KeepShocks,:);
      etaShocks(j,:) = zeros(1,h);
      %
      % inistalize the state
      %
      KsiT = KsiTTz;
      for i=1:h;
         KsiT = (F*KsiT)+(B0H*etaShocks(:,i));
         if k>0;
            PredData.ShockScenarioResults(:,i,j) = (A'*X(:,i))+(H'*KsiT);
         else;
            PredData.ShockScenarioResults(:,i,j) = H'*KsiT;
         end;
      end;
   end;
else;
   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;
         PredData.etaShocks(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*A'*X(:,i)) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiT) - K2sum -U(:,g));
      else;
         PredData.etaShocks(DSGEModel.ConditionalShocks,i) = etaMat*(Z(:,i) - (K1'*H(:,:,DSGEModel.LastPeriod+i)'*F*KsiT) - K2sum -U(:,g));
      end;
      %
      % update the states
      %
      KsiT = (F*KsiT)+(B0*PredData.etaShocks(:,i));
      PredData.xiStates(:,i) = KsiT;
      %
      % compute the conditional mean predictions
      %
      if k>0;
         PredData.PredMean(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiT);
      else;
         PredData.PredMean(:,i) = H(:,:,DSGEModel.LastPeriod+i)'*KsiT;
      end;
      %
      % update the K2sum measure
      %
      K2sum = K2sum+(K2'*PredData.PredMean(:,i));
   end;
   for i=g+1:h;
      %
      % update the states
      %
      KsiT = (F*KsiT);
      PredData.xiStates(:,i) = KsiT;
      %
      % compute the conditional mean predictions
      %
      if k>0;
         PredData.PredMean(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiT);
      else;
         PredData.PredMean(:,i) = H(:,:,DSGEModel.LastPeriod+i)'*KsiT;
      end;
   end;
   [B0H,KeepShocks] = RemoveRedundantColumns(B0);
   q_used = size(B0H,2);
   PredData.ShockScenarioResults = zeros(DSGEModel.n,h,q_used);
   for j=1:q_used;
      %
      % set shock j of the shocks manipulated for the theoretical mean of the conditional
      % forecast
      %
      etaShocks = PredData.etaShocks(KeepShocks,:);
      etaShocks(j,:) = zeros(1,h);
      %
      % inistalize the state
      %
      KsiT = KsiTTz;
      for i=1:h;
         KsiT = (F*KsiT)+(B0H*etaShocks(:,i));
         if k>0;
            PredData.ShockScenarioResults(:,i,j) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiT);
         else;
            PredData.ShockScenarioResults(:,i,j) = H(:,:,DSGEModel.LastPeriod+i)'*KsiT;
         end;
      end;
   end;
end;
%
% annualize the mean predictions when needed
%
PredData.OrigPredMean = PredData.PredMean;
if AnnualizeData==1;
   PredData.PredMean = Annualize([Yinit PredData.PredMean],DSGEModel.annual,DSGEModel.annualscale);
   if isempty(YObserved)==0;
      YObserved = Annualize([Yinit YObserved],DSGEModel.annual,DSGEModel.annualscale);;
   end;
elseif TransData==1;
   YTmp = [Yinit PredData.PredMean];
   [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
   Tt = size(YData,1);
   PredData.PredMean = YData(Tt-h+1:Tt,:)';
   if isempty(YObserved)==0;
      YObserved = [Yinit YObserved];
      [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YObserved',DSGEModel.YTransformation);
      Tt = size(YData,1);
      YObserved = YData(Tt-h+1:Tt,:)';
   end;
end;
%
% check if we should setup a wait handle
%
if DSGEModel.ShowProgress==1;
   ProgressStructure.title = ['Progress for ' num2str(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 = ['Conditional Predictive Distribution' NameStr];
   ProgressStructure.CurrINI = CurrINI;
   WaitHandle = ProgressDLG(0,ProgressStructure);
   set(WaitHandle,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
end;
%
% start the iterations
%
it = 1;
abort = '';
if length(size(H))==2;
   for it=1:NumPaths;
      %
      % check if we should update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(WaitHandle,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            PredPaths = [];
            break;
         else;
            ProgressDLG([it/NumPaths it]);
         end;
      end;
      %
      % Draw a state vector for time T using the N(KsiTTz,PTTz) distribution.
      %
      KsiT = KsiTTz+(S*sqrt(L)*randn(size(L,2),1));
      PredData.HistState(:,it) = KsiT;
      KsiCTmp = KsiTT;
      KsiTmp = KsiTT;
      KsiTLZ = KsiTT;
      PTTg = zeros(r,r);
      PTTgLZ = PTTg;
      Fg = eye(r);
      Q = B0*B0';
      QM = B0*M*M'*B0';
      %
      % Draw a sequence of economic shocks and set the conditioning shocks to zero
      %
      etaTh = randn(q,h);
      etaCTmp = etaTh;
      etaTmp = etaTh;
      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);
         PredData.wMean = PredData.wMean+((1/NumPaths)*wTh);
         PredData.MeasError(:,:,it) = wTh;
      end;
      %
      % initialize the K2 part
      %
      K2sum = zeros(m,1);
      K2sumCTmp = K2sum;
      K2sumLZ = zeros(m,1);
      %
      % 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*KsiCTmp)+(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;
         Fg = Fg*F;
         PTTg = (F*PTTg*F')+Q;
         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;
               K2sumLZ = K2sumLZ+(K2'*((A'*X(:,i))+(H'*KsiTLZ)+wTh(:,i)));
            else;
               K2sumLZ = K2sumLZ+(K2'*((H'*KsiTLZ)+wTh(:,i)));
            end;
            PTTgLZ = (F*PTTgLZ*F')+QM;
         end;
      end;
      PredData.Shocks(:,:,it) = etaTh;
      PredData.etaMean = PredData.etaMean+((1/NumPaths)*etaTh);
      %
      % 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;
      end;
      PredData.MeanPredPath = PredData.MeanPredPath+((1/NumPaths)*PredPaths(:,:,it));
      if AnnualizeData==1;
         PredPaths(:,:,it) = Annualize([Yinit PredPaths(:,:,it)],DSGEModel.annual,DSGEModel.annualscale);
      else;
         if TransData==1;
            %
            % transform the data
            %
            YTmp = [Yinit PredPaths(:,:,it)];
            [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
            PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
         end;
         %
         % compute multivariate modesty statistic
         %
         if k>0;
            PredictCTmp = (A'*X(:,g))+(H'*KsiCTmp)+wTh(:,g);
            PredictY = (A'*X(:,g))+(H'*Fg*KsiTT);
         else;
            PredictCTmp = (H'*KsiCTmp)+wTh(:,g);
            PredictY = H'*Fg*KsiTT;
         end;
         PredCov = (H'*PTTg*H)+R;
         PredDiff = PredictCTmp-PredictY;
         MultiModestyStat(it,1) = PredDiff'*inv(PredCov)*PredDiff;
         UniModestyStat(it,:) = (PredDiff ./ sqrt(diag(PredCov)))';
         %
         % comparison-path for multivariate statistic
         %
         if k>0;
            PredictTmp = (A'*X(:,g))+(H'*KsiTmp)+wTh(:,g);
         else;
            PredictTmp = (H'*KsiTmp)+wTh(:,g);
         end;
         PredDiff = PredictTmp-PredictY;
         MultiModestyStat(it,2) = PredDiff'*inv(PredCov)*PredDiff;
         if it==1;
            if k>0;
               PredDiff = (A'*X(:,g))+(H'*KsiTLZ)+wTh(:,g)-PredictY;
            else;
               PredDiff = (H'*KsiTLZ)+wTh(:,g)-PredictY;
            end;
            CovPredDiff = diag(H'*PTTgLZ*H);
            for i=1:length(CovPredDiff);
               if CovPredDiff(i)==0;
                  PredDiff(i) = 0;
                  CovPredDiff(i) = 1;
               end;
            end;
            UniModestyStatLZ(1,:) = (PredDiff ./ sqrt(CovPredDiff))';
         end;
      end;
      drawnow;
   end;
else;
   for it=1:NumPaths;
      %
      % check if we should update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(WaitHandle,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            PredPaths = [];
            break;
         else;
            ProgressDLG([it/NumPaths it]);
         end;
      end;
      %
      % Draw a state vector for time T using the N(KsiTTz,PTTz) distribution.
      %
      KsiT = KsiTTz+(S*sqrt(L)*randn(size(L,2),1));
      PredData.HistState(:,it) = KsiT;
      KsiCTmp = KsiTT;
      KsiTmp = KsiTT;
      KsiTLZ = KsiTT;
      PTTg = zeros(r,r);
      PTTgLZ = PTTg;
      Fg = eye(r);
      Q = B0*B0';
      QM = B0*M*M'*B0';
      %
      % Draw a sequence of economic shocks and set the conditioning shocks to zero
      %
      etaTh = randn(q,h);
      etaCTmp = etaTh;
      etaTmp = etaTh;
      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);
         PredData.wMean = PredData.wMean+((1/NumPaths)*wTh);
         PredData.MeasError(:,:,it) = wTh;
      end;
      %
      % initialize the K2 part
      %
      K2sum = zeros(m,1);
      K2sumCTmp = K2sum;
      K2sumLZ = zeros(m,1);
      %
      % 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*KsiCTmp)+(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;
         Fg = Fg*F;
         PTTg = (F*PTTg*F')+Q;
         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;
               K2sumLZ = K2sumLZ+(K2'*((A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTLZ)+wTh(:,i)));
            else;
               K2sumLZ = K2sumLZ+(K2'*((H(:,:,DSGEModel.LastPeriod+i)'*KsiTLZ)+wTh(:,i)));
            end;
            PTTgLZ = (F*PTTgLZ*F')+QM;
         end;
      end;
      PredData.Shocks(:,:,it) = etaTh;
      PredData.etaMean = PredData.etaMean+((1/NumPaths)*etaTh);
      %
      % 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;
      end;
      PredData.MeanPredPath = PredData.MeanPredPath+((1/NumPaths)*PredPaths(:,:,it));
      if AnnualizeData==1;
         PredPaths(:,:,it) = Annualize([Yinit PredPaths(:,:,it)],DSGEModel.annual,DSGEModel.annualscale);
      else;
         if TransData==1;
            %
            % transform the data
            %
            YTmp = [Yinit PredPaths(:,:,it)];
            [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
            PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
         end;
         %
         % compute multivariate modesty statistic
         %
         if k>0;
            PredictCTmp = (A'*X(:,g))+(H(:,:,DSGEModel.LastPeriod+g)'*KsiCTmp)+wTh(:,g);
            PredictY = (A'*X(:,g))+(H(:,:,DSGEModel.LastPeriod+g)'*Fg*KsiTT);
         else;
            PredictCTmp = (H(:,:,DSGEModel.LastPeriod+g)'*KsiCTmp)+wTh(:,g);
            PredictY = H(:,:,DSGEModel.LastPeriod+g)'*Fg*KsiTT;
         end;
         PredCov = (H(:,:,DSGEModel.LastPeriod+g)'*PTTg*H(:,:,DSGEModel.LastPeriod+g))+R;
         PredDiff = PredictCTmp-PredictY;
         MultiModestyStat(it,1) = PredDiff'*inv(PredCov)*PredDiff;
         UniModestyStat(it,:) = (PredDiff ./ sqrt(diag(PredCov)))';
         %
         % comparison-path for multivariate statistic
         %
         if k>0;
            PredictTmp = (A'*X(:,g))+(H(:,:,DSGEModel.LastPeriod+g)'*KsiTmp)+wTh(:,g);
         else;
            PredictTmp = (H(:,:,DSGEModel.LastPeriod+g)'*KsiTmp)+wTh(:,g);
         end;
         PredDiff = PredictTmp-PredictY;
         MultiModestyStat(it,2) = PredDiff'*inv(PredCov)*PredDiff;
         if it==1;
            if k>0;
               PredDiff = (A'*X(:,g))+(H(:,:,DSGEModel.LastPeriod+g)'*KsiTLZ)+wTh(:,g)-PredictY;
            else;
               PredDiff = (H(:,:,DSGEModel.LastPeriod+g)'*KsiTLZ)+wTh(:,g)-PredictY;
            end;
            CovPredDiff = diag(H(:,:,DSGEModel.LastPeriod+g)'*PTTgLZ*H(:,:,DSGEModel.LastPeriod+g));
            for i=1:length(CovPredDiff);
               if CovPredDiff(i)==0;
                  PredDiff(i) = 0;
                  CovPredDiff(i) = 1;
               end;
            end;
            UniModestyStatLZ(1,:) = (PredDiff ./ sqrt(CovPredDiff))';
         end;
      end;
      drawnow;
   end;
end;
%
% close the wait handle
%
if DSGEModel.ShowProgress==1;
   set(WaitHandle,'UserData','done');
   close(WaitHandle);
   drawnow;
   pause(0.02);
end;
if (strcmp(abort,'cancel')==0)&(DSGEModel.RunPredictionEvent==1);
   MeanPrediction = zeros(DSGEModel.n,h);
   for i=1:NumPaths;
      MeanPrediction = MeanPrediction+((1/NumPaths)*PredPaths(:,:,i));
   end;
   if DSGEModel.AdjustPredictionPaths==1;
      %
      % we need to fix the paths such that the mean is correct
      %
      AdjustMean = PredData.PredMean-MeanPrediction;
      for i=1:NumPaths;
         PEPredPaths(:,:,i) = PredPaths(:,:,i)+AdjustMean;
      end;
   else;
      PEPredPaths = PredPaths;
   end;
   %
   % compute the prediction events results
   %
   PredEventData = CalculatePredictionEvents(PEPredPaths,DSGEModel.PredictionEvent);
   %
   if isempty(YObserved)==0;
      AdjustMean = YObserved-MeanPrediction;
      for i=1:NumPaths;
         PEPredPaths(:,:,i) = PredPaths(:,:,i)+AdjustMean;
      end;
      YObsEventData = CalculatePredictionEvents(PEPredPaths,DSGEModel.PredictionEvent);
   end;
end;
%
% check if we should compute kernel density estimates
%
if (strcmp(abort,'cancel')==0)&(NumPaths>=500);
   if DSGEModel.RunPredictionEvent==0;
      if DSGEModel.AdjustPredictionPaths==1;
         MeanPrediction = zeros(DSGEModel.n,h);
         for i=1:NumPaths;
            MeanPrediction = MeanPrediction+((1/NumPaths)*PredPaths(:,:,i));
         end;
         AdjustMean = PredData.PredMean-MeanPrediction;
         for i=1:NumPaths;
            PEPredPaths(:,:,i) = PredPaths(:,:,i)+AdjustMean;
         end;
      else;
         PEPredPaths = PredPaths;
      end;
   end;
   %
   [PredData.KernelX,PredData.KernelY] = PredictionKernelDensity(PEPredPaths);
end;

%
% end of DSGECondPredictionPathsTheta.m
%
