function [h , y , yTotalRMSE, yRMSE] = garchpred(spec , y , horizon , X , XF)
%GARCHPRED Univariate GARCH process forecasting.
%   Given specifications for the conditional mean and variance of an observed 
%   univariate return series, forecast the conditional mean and standard 
%   deviation of the return series NUMPERIODS into the future. Additionally, 
%   volatility forecasts of asset returns over multi-period holding intervals, 
%   and the standard errors of conditional mean forecasts are also computed. 
%   The conditional mean and variance may be of general ARMAX and GARCH form, 
%   respectively.
%   
%   [SigmaForecast, MeanForecast, ...
%    SigmaTotal   , MeanRMSE    ] = garchpred(Spec , Series)
%
%   [SigmaForecast, MeanForecast, ...
%    SigmaTotal   , MeanRMSE]     = garchpred(Spec, Series, NumPeriods, X, XF)
%
%   Optional Inputs: NumPeriods, X, XF
%
% Inputs:
%   Spec - Structure specification for the conditional mean and variance models.
%     Spec is a MATLAB structure with fields generated by calling the function 
%     GARCHSET, or the output of the estimation function GARCHFIT. For details, 
%     type "help garchset" or "help garchfit".
%
%   Series - Matrix of observations of the underlying univariate return series 
%     of interest. Series is the response variable representing the return 
%     series fit to conditional mean and variance specifications. Each column 
%     of Series in an independent realization (i.e., path). The last row of 
%     Series holds the most recent observation of each realization. Series is 
%     assumed to be a stationary stochastic process, and the ARMA component of 
%     the conditional mean model (if any) is assumed to stationary and 
%     invertible.
%
% Optional Inputs:
%   NumPeriods - Positive, scalar integer representing the forecast horizon of
%     interest, expressed in periods compatible with the sampling frequency 
%     of the input Series. If empty or missing, the default is 1.
%
%   X - Time series regression matrix of observed explanatory data. Typically, 
%     X is a matrix of observed asset returns (e.g., the return series of an 
%     equity index), and represents the 'past history' of the explanatory data. 
%     Each column of X is an individual explanatory variable in the regression 
%     component of the conditional mean. In each column of X, the first row 
%     contains the oldest observation and the last row the most recent. If X is 
%     specified, the most recent number of valid (non-NaN) observations in each 
%     column of X must equal or exceed the most recent number of valid 
%     observations in Series. When the number of valid observations in each 
%     column of X exceeds that of Series, only the most recent observations of 
%     X are used. If empty or missing, the conditional mean will have no 
%     regression component. 
%
%   XF - Time series matrix of forecasted explanatory data. XF represents the 
%     evolution of the same explanatory data found in X projected into the 
%     'future'. Thus, XF and X must have the same number of columns. In each 
%     column of XF, the first row contains the 1-period-ahead forecast, the 
%     second row the 2-period-ahead forecast, and so on. If XF is specified, 
%     the number of rows (forecasts) in each column (time series) of XF must 
%     equal or exceed the forecast horizon NUMPERIODS. When the number of 
%     forecasts in XF exceeds NUMPERIODS, only the first NUMPERIODS forecasts 
%     are used. If empty or missing, the conditional mean forecast will have 
%     no regression component.
%
% Outputs:
%   SigmaForecast - Matrix of minimum mean square error (MSE) forecasts of 
%     the conditional standard deviation of Series on a per period basis. 
%     SigmaForecast will have NUMPERIODS rows and the same number of columns as
%     Series. The first row contains the volatility forecast in the 1st period 
%     for each realization of Series, the second row contains the volatility 
%     forecast in the 2nd period, and so on. Thus, if a forecast horizon 
%     greater than one is specified (NUMPERIODS > 1), the per-period forecasts 
%     of all intermediate horizons are returned as well; in this case, the last
%     row will contain the forecast at the specified horizon.
%
%   MeanForecast - Matrix of minimum mean square error (MSE) forecasts of the
%     conditional mean of Series on a per period basis. MeanForecast will be 
%     the same size as SigmaForecast. The first row contains the forecast in 
%     the 1st period for each realization of Series, the second row contains 
%     the forecast in the 2nd period, and so on.
%
%   SigmaTotal - Matrix of volatility forecasts of Series over multi-period 
%     holding intervals. The first row contains the standard deviation of returns
%     expected for assets held for 1 period for each realization of Series, the 
%     second row contains the standard deviation of returns expected for assets 
%     held for 2 periods, and so on. Thus, the last row will contain the 
%     volatility forecast of the cumulative return that would be obtained if an 
%     asset was held for the entire NUMPERIODS forecast horizon. These forecasts
%     are correct for continuously-compounded returns, and approximate for 
%     periodically compounded returns. SigmaTotal is the same size as 
%     SigmaForecast provided the conditional mean is modeled as a stationary/
%     invertible ARMA process; for conditional mean models with regression 
%     components (i.e., X and/or XF are specified), an empty matrix, [], is 
%     returned.
%
%   MeanRMSE - Matrix of root mean square errors (RMSE) associated with 
%     MeanForecast. That is, MeanRMSE is the conditional standard deviation 
%     of the forecast errors (i.e., the standard error of the forecast) of the
%     corresponding MeanForecast matrix. MeanRMSE is the same size as 
%     MeanForecast and is organized in exactly the same manner, provided the
%     conditional mean is modeled as a stationary/invertible ARMA process; for 
%     conditional mean models with regression components (i.e., X and/or XF 
%     are specified), an empty matrix, [], is returned.
%
% Notes:
%   (1) Since a complete conditional mean specification is required to correctly
%     infer the innovations process which drives the forecasts, the regression 
%     matrix of observed returns (X), is typically the same X (if any) used for 
%     simulation (via GARCHSIM) and/or estimation (via GARCHFIT). XF, however, 
%     is just the forecast of X, and is ONLY needed for forecasting the 
%     conditional mean (MeanForecast). Thus, if only the conditional variance 
%     forecast (SigmaForecast) is needed, XF is unnecessary. Furthermore, 
%     although X may be specified without XF, the converse is NOT true. XF may 
%     be specified ONLY when X is specified.
%   (2) This function calls the function GARCHINFER to access the past history 
%     of innovations and conditional standard deviations inferred from Series. 
%     If the innovations and conditional standard deviations are needed, call 
%     GARCHINFER directly.
%
% See also GARCHSET, GARCHSIM, GARCHINFER, GARCHFIT.

% Copyright 1999-2002 The MathWorks, Inc.   
% $Revision: 1.10 $   $ Date: 1998/01/30 13:45:34 $

%
% References:
%   Baillie, R.T., Bollerslev, T. (1992), "Prediction in Dynamic Models with 
%     Time-Dependent Conditional Variances", Journal of Econometrics, 
%     vol. 52, pp. 91-113.
%   Bollerslev, T. (1986), "Generalized Autoregressive Conditional 
%     Heteroskedasticity", Journal of Econometrics, vol. 31, pp. 307-327.
%   Box, G.E.P., Jenkins, G.M., Reinsel, G.C., "Time Series Analysis: 
%     Forecasting and Control", 3rd edition, Prentice Hall, 1994.
%   Engle, Robert (1982), "Autoregressive Conditional Heteroskedasticity 
%     with Estimates of the Variance of United Kingdom Inflation", 
%     Econometrica, vol. 50, pp. 987-1007.
%   Hamilton, J.D., "Time Series Analysis", Princeton University Press, 1994.
%

%
% Test for stationarity & invertibility of the ARMA model. This should be 
% done here because the input SPEC structure may be the output of the
% estimation function GARCHFIT which does NOT guarantee stationarity or
% invertibility of the ARMA model. 
%

AR  =  garchget(spec , 'AR');      % Conditional mean AR coefficients.
MA  =  garchget(spec , 'MA');      % Conditional mean MA coefficients.

if any(abs(roots([1 ; -AR(:)])) >= 1)
   error(' Auto-Regressive Polynomial must be Stationary.')
end

if any(abs(roots([1 ;  MA(:)])) >= 1)
   error(' Moving-Average Polynomial must be Invertible.')
end

%
% Check & scrub the observed return series matrix y(t).
%

if (nargin < 2)

   error(' Observed return series ''Series'' must be specified.');

else

   rowY  =  logical(0);     

   if prod(size(y)) == length(y)   % Check for a vector (single return series).
      rowY  =  size(y,1) == 1;     % Flag a row vector for outputs.
      y     =  y(:);               % Convert to a column vector.
   end

%
%  The following code segment assumes that missing observations are indicated
%  by the presence of NaN's. Any initial rows with NaN's are removed, and 
%  processing proceeds with the remaining block of contiguous non-NaN rows. 
%  Put another way, NaN's are allowed, but they MUST appear as a contiguous 
%  sequence in the initial rows of the y(t) matrix. Since the log-likelihood 
%  functions are designed to process y(t) as a matrix (instead of individual 
%  vectors!), any initial rows with NaN's are stripped. Thus, realizations 
%  with no missing observations will lose data if other realizations have 
%  missing values.
%
   i1  =  find(isnan(y));
   i2  =  find(isnan(diff([y ; zeros(1,size(y,2))]) .* y));

   if (length(i1) ~= length(i2)) | any(i1 - i2)
      error(' Only initial observations in ''Series'' may be missing (NaN''s).')
   end

   if any(sum(isnan(y)) == size(y,1))
      error(' A realization of ''Series'' is completely missing (all NaN''s).')
   end

   firstValidRow  =  max(sum(isnan(y))) + 1;
   y              =  y(firstValidRow:end , :);

end

%
% Check input parameters and set defaults.
%

if (nargin >= 3) & ~isempty(horizon)
   if prod(size(horizon)) > 1
      error(' Forecast horizon ''NumPeriods'' must be a scalar.');
   end
   if (round(horizon) ~= horizon) | (horizon <= 0)
      error(' Forecast horizon ''NumPeriods'' must be a positive integer.');
   end
else
   horizon  =  1;   % Set default.
end

%
% Scrub the regression matrix, and ensure the observed return series matrix y(t) 
% and the regression matrix X(t) have the same number of valid (i.e., non-NaN)
% rows (i.e., impose time index compatibility).
%

if (nargin >= 4) & ~isempty(X)

   if prod(size(X)) == length(X)   % Check for a vector.
      X  =  X(:);                  % Convert to a column vector.
   end
%
%  Retain the last contiguous block of non-NaN (i.e, non-missing valued) observations only. 
%
   if any(isnan(X(:)))
      X  =  X((max(find(isnan(sum(X,2)))) + 1):end , :);
   end

   if size(X,1) < size(y,1)
      error(' Regression matrix ''X'' has insufficient number of observations.');
   else
      X  =  X(size(X,1) - (size(y,1) - 1):end , :);    % Retain only the most recent samples.
   end
%
%  Ensure number of regression coefficients match number of regressors.
%
   regress =  garchget(spec , 'Regress'); % Conditional mean regression coefficients.

   if size(X,2) ~= length(regress)
      error(' Number of ''Regress'' coefficients unequal to number of regressors in ''X''.');
   end

else

   X        =  [];   % Ensure X exists.
   regress  =  [];

end

%
% Scrub the regression forecast matrix XF(t) and ensure sufficient forecasts exist 
% to (at least) cover the forecast horizon. As opposed to X(t), which may contain 
% some missing data (i.e., NaN's), XF(t) is created by the user and should contain
% no missing data (i.e., XF(t) should be NaN-free). Also, XF(t) and X(t) MUST have 
% the number of columns (i.e., explanatory variables).
%

if (nargin >= 5) & ~isempty(XF)

%
%  Check XF(t) versus X(t).
%
   if isempty(X) & ~isempty(XF)
      error(' Forecast matrix ''XF'' cannot be specified without ''X''.');
   end

   if size(X,2) ~= size(XF,2)
      error(' Matrices ''XF'' and ''X'' must have the same number of columns.');
   end