function [PredPaths,PredData,PredEventData,YObsEventData,status,kalmanstatus] = DSGEPredictionPathsTheta(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,X,h,NumPaths,AnnualizeData,TransData,NameStr)
% DSGEPredictionPathsTheta: Calculates a certain number of possible paths that the endogenous variables
%                           can take out-of-sample for a given value of the parameters in the DSGE model.
%
% USAGE:
%
%       [PredPaths,PredData,PredEventData,YObsEventData,status,kalmanstatus] = DSGEPredictionPathsTheta(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,...
%                                                                  X,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 kxMaxPH with the observations of the exogenous variables for
%                                 the prediction paths.
%
%                  h (integer) with the prediction horizon.
%
%                  NumPaths (integer)the number of paths to draw.
%
%                  AnnualizeData (boolean) that is 1 if the prediction paths should be annualized and 0 otherwise.
%
%                  TransData (boolean) that is 1 if the prediction paths should be transformed and 0 otherwise.
%
%                  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.
%                                 etaShocks  - (qxh) matrix with the values of the shocks at the mean
%                                              predictions.
%                                 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 holds 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.
%
% 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.
%                                      mocde = 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: February 16, 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:
%
% * 19-02-2007: Finished the function.
%
% * 21-02-2007: Added the outputs "StateCov", "ShockCov", "MeasureCov", and "PredictY".
%
% * 22-02-2007: Updated the documentation.
%
% * 26-02-2007: Removed the outputs "StateCov", "ShockCov", "MeasureCov", "PredictY", and "WaitHandle". The
%               inputs "WaitHandle" was removed as were inputs measuring the current parameter draw (always 1 before)
%               and the total number of draws (always 1 before) since they were made redundant.
%
% * 03-05-2007: Added the output "PredEventData".
%
% * 15-05-2007: Added an extra "drawnow" after shutting down the progress dialog.
%
% * 06-06-2007: Changed "predictive paths" to "prediction paths" in the title for the progress dialog.
%
% * 21-08-2007: Added the output "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.
%
% * 18-12-2007: Added the input "h". This makes it possible to allow for X being empty.
%
% * 21-12-2007: Added the code for dealing with unit-root state variables when running the Kalman filter.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 16-04-2008: The code now checks if the model has a unique convergent solution also at the posterior mode
%               estimates.
%
% * 02-05-2008: Added the "TransData" input variable.
%
% * 05-05-2008: Made sure that "PredData.PredMean" is tarnsformed when "TransData" is 1.
%
% * 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 "etaShocks" and "xiStates" 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.
%
% * 09-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 15-07-2008: Fixed an initialization bug for wTh.
%
% * 15-01-2009: Added the fields "ShockNames", "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 = [];
%
% set up the ModelParameters structure based on theta and thetaPositions
%
if nargout==4;
   %
   % 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>4;
      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==4;
   [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>4;
   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);
%
% 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;
%
% 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);
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 and forecasts (T+h = ' int2str(DSGEModel.LastPeriod+h) '). YADA has therefore aborted from the unconditional predictions.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
% 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>4;
   kalmanstatus = stat;
   if kalmanstatus~=0;
      return;
   end;
end;
%
% Now we need to calculate KsiTT and PTT for just the last period
%
[r,T] = size(Ksitt1);
q = size(B0,2);
if length(size(H))==2;
   inv_msey = inv((H'*Ptt1(:,:,T)*H)+R);
   KsiTT = Ksitt1(:,T)+(Ptt1(:,:,T)*H*inv_msey*(DSGEModel.Y(:,DSGEModel.LastPeriod)-Yhat(:,T)));
   PTT = Ptt1(:,:,T)-(Ptt1(:,:,T)*H*inv_msey*H'*Ptt1(:,:,T));
else;
   inv_msey = inv((H(:,:,DSGEModel.LastPeriod)'*Ptt1(:,:,T)*H(:,:,DSGEModel.LastPeriod))+R);
   KsiTT = Ksitt1(:,T)+(Ptt1(:,:,T)*H(:,:,DSGEModel.LastPeriod)*inv_msey*(DSGEModel.Y(:,DSGEModel.LastPeriod)-Yhat(:,T)));
   PTT = Ptt1(:,:,T)-(Ptt1(:,:,T)*H(:,:,DSGEModel.LastPeriod)*inv_msey*H(:,:,DSGEModel.LastPeriod)'*Ptt1(:,:,T));
end;
%
% use the Schur decomposition for the reduced rank matrix PTT instead of
% the Choleski decomposition
%
[S,L] = SchurDecomposition(PTT);
%
% check if we have measurement errors and if so compute a Schur decomposition of the R matrix
%
if rank(R)>0;
   [Rm,Lm] = SchurDecomposition(R);
   Rmeasure = Rm*sqrt(Lm);
   NumMeasure = size(Lm,1);
else;
   Rmeasure = [];
   wTh = zeros(DSGEModel.n,h);
   NumMeasure = 0;
end;
%
% set up the output
%
k = DSGEModel.k;
PredPaths = zeros(DSGEModel.n,h,NumPaths);
%
% 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;
elseif TransData==1;
   Yinit = DSGEModel.Y(:,1:DSGEModel.LastPeriod);
end;
%
% add info to the PredData structure
%
PredData.KernelX = [];
PredData.KernelY = [];
PredData.X = X;
PredData.A = A;
PredData.H = H;
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);
KsiT = KsiTT;
PredData.etaShocks = zeros(q,h);
PredData.xiStates = zeros(r,h);
if length(size(H))==2;
   for i=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;
else;
   for i=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)'*KsiT);
      else;
         PredData.PredMean(:,i) = H(:,:,DSGEModel.LastPeriod)'*KsiT;
      end;
   end;
end;
%
% annualize the mean predictions when needed
%
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 = ['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(KsiTT,PTT) distribution.
      %
      KsiT = KsiTT+(S*sqrt(L)*randn(size(L,2),1));
      PredData.HistState(:,it) = KsiT;
      %
      % Draw a sequence of state shocks
      %
      etaTh = randn(q,h);
      PredData.Shocks(:,:,it) = etaTh;
      vTh = B0*etaTh;
      PredData.etaMean = PredData.etaMean+((1/NumPaths)*etaTh);
      %
      % calculate a path for the states
      %
      KsiTh = zeros(r,h);
      KsiTh(:,1) = (F*KsiT)+vTh(:,1);
      for i=2:h;
         KsiTh(:,i) = (F*KsiTh(:,i-1))+vTh(:,i);
      end;
      %
      % Draw a sequence of measurement errors
      %
      if isempty(Rmeasure)==0;
         wTh = Rmeasure*randn(NumMeasure,h);
         PredData.wMean = PredData.wMean+((1/NumPaths)*wTh);
         PredData.MeasError(:,:,it) = wTh;
      end;
      %
      % add everything
      %
      if AnnualizeData==0;
         if TransData==0;
            if k>0;
               PredPaths(:,:,it) = (A'*X)+(H'*KsiTh)+wTh;
            else;
               PredPaths(:,:,it) = (H'*KsiTh)+wTh;
            end;
         else;
            if k>0;
               YTmp = [Yinit (A'*X)+(H'*KsiTh)+wTh];
            else;
               YTmp = [Yinit (H'*KsiTh)+wTh];
            end;
            %
            % transform the data
            %
            [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
            PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
         end;
      else;
         if k>0;
            PredPaths(:,:,it) = Annualize([Yinit (A'*X)+(H'*KsiTh)+wTh],DSGEModel.annual,DSGEModel.annualscale);
         else;
            PredPaths(:,:,it) = Annualize([Yinit (H'*KsiTh)+wTh],DSGEModel.annual,DSGEModel.annualscale);
         end;
      end;
      drawnow;
   end;
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(KsiTT,PTT) distribution.
      %
      KsiT = KsiTT+(S*sqrt(L)*randn(size(L,2),1));
      PredData.HistState(:,it) = KsiT;
      %
      % Draw a sequence of state shocks
      %
      etaTh = randn(q,h);
      PredData.Shocks(:,:,it) = etaTh;
      vTh = B0*etaTh;
      PredData.etaMean = PredData.etaMean+((1/NumPaths)*etaTh);
      %
      % calculate a path for the states
      %
      KsiTh = zeros(r,h);
      KsiTh(:,1) = (F*KsiT)+vTh(:,1);
      for i=2:h;
         KsiTh(:,i) = (F*KsiTh(:,i-1))+vTh(:,i);
      end;
      %
      % Draw a sequence of measurement errors
      %
      if isempty(Rmeasure)==0;
         wTh = Rmeasure*randn(NumMeasure,h);
         PredData.wMean = PredData.wMean+((1/NumPaths)*wTh);
         PredData.MeasError(:,:,it) = wTh;
      end;
      %
      % add everything
      %
      if AnnualizeData==0;
         if TransData==0;
            if k>0;
               for i=1:h;
                  PredPaths(:,i,it) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            else;
               for i=1:h;
                  PredPaths(:,i,it) = (H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            end;
         else;
            YTmp = zeros(size(Yinit,1),h);
            if k>0;
               for i=1:h;
                  YTmp(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            else;
               for i=1:h;
                  YTmp(:,i) = (H(:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
               end;
            end;
            YTmp = [Yinit YTmp];
            %
            % transform the data
            %
            [YTStatus,YData,ErrorStr,HeaderStr] = CalcTransformData(DSGEModel.VariableNames,YTmp',DSGEModel.YTransformation);
            PredPaths(:,:,it) = YData(Tt-h+1:Tt,:)';
         end;
      else;
         YTmp = zeros(size(Yinit,1),h);
         if k>0;
            for i=1:h;
               YTmp(:,i) = (A'*X(:,i))+(H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
            end;
         else;
            for i=1:h;
               YTmp(:,i) = (H(:,:,DSGEModel.LastPeriod+i)'*KsiTh(:,i))+wTh(:,i);
            end;
         end;
         PredPaths(:,:,it) = Annualize([Yinit YTmp],DSGEModel.annual,DSGEModel.annualscale);
      end;
      drawnow;
   end;
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;
         YPEPredPaths(:,:,i) = PredPaths(:,:,i)+AdjustMean;
      end;
      YObsEventData = CalculatePredictionEvents(YPEPredPaths,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 DSGEPredictionPathsTheta.m
%
