%
%  Ensure no missing values (i.e., NaN's) exist in XF(t).
%
   if any(isnan(XF(:)))
      error(' Regression forecast matrix ''XF'' has missing data (i.e., NaN''s).');
   end

   if size(XF,1) < horizon
      error(' Regression matrix ''XF'' has insufficient number of forecasts.');
   else
      XF  =  XF(1:horizon , :);      % Retain only the initial forecasts.
   end

else

   XF  =  [];   % Ensure XF exists.

end

%
% Extract model orders for the conditional mean & variance processes.
%

R  =  garchget(spec , 'R');       % Conditional mean AR order.
M  =  garchget(spec , 'M');       % Conditional mean MA order.
P  =  garchget(spec , 'P');       % Conditional variance order for lagged variances.
Q  =  garchget(spec , 'Q');       % Conditional variance order for lagged squared residuals.

%
% Ensure all coefficients of the conditional variance
% model exist and have proper dimensions.
%

K      =  garchget(spec , 'K');       % Conditional variance constant.
GARCH  =  garchget(spec , 'GARCH');   % Conditional variance coefficients for lagged variances.
ARCH   =  garchget(spec , 'ARCH');    % Conditional variance coefficients for lagged squared residuals.

if isempty(K)
   error(' Conditional variance constant ''K'' must be specified.');
end
if isempty(GARCH) & (P ~= 0)
   error(' ''GARCH'' coefficients of lagged variances must be specified.');
end
if isempty(ARCH)  & (Q ~= 0)
   error(' ''ARCH'' coefficients of lagged squared residuals must be specified.');
end

%
% Let y(t) = return series of interest (assumed stationary)
%     e(t) = innovations of the model noise process (assumed invertible)
%     h(t) = conditional variance of the innovations process e(t)
%
% ARMA(R,M)/GARCH(P,Q) processing requires pre-sample values for conditioning. 
%
% We require R pre-sample lags of y(t), max(M,Q) pre-sample lags of e(t), and 
% P pre-sample lags of h(t). To be safe, work with max([R M P Q]) pre-sample lags
% for all processes in the ARMA(R,M)/GARCH(P,Q) model.
%

preSamples  =  max([R M P Q]);   % # of pre-sample lags for 'jump-starting' the process.

%
% Set the storage size of y(t), h(t), and e(t) required for forecasting.
%

T  =  horizon  +  preSamples;

%
% Given specifications for the conditional mean & variance models, infer the 
% innovations, e(t), and conditional variance, h(t), time series from the
% observed return series y(t). The code segment below represents a re-construction 
% of the 'past' e(t) and h(t) processes, which are required for forecasting into 
% the 'future'. 
%
% Note that GARCHINFER actually returns the conditional standard deviation, 
% rather than the conditional variance. The output h(t) is converted later 
% on instead of immediately for performance purposes only.
%

[e , h] =  garchinfer(spec , y , X);  % At this point, h(t) = standard deviation.

%
% Re-format the y(t), h(t), and e(t) processes for the purpose of forecasting. 
% For forecasting, only the most recent R lags of y(t), max(M,Q) lags of e(t), 
% and P lags of h(t) are required. Thus, truncate y(t), h(t), and e(t), retaining
% only the most recent 'preSamples' values, then append 'horizon' zeros to 
% accommodate the user-specified forecast horizon. Since forecasting requires 
% squared innovations e(t)^2, work with them directly instead of the e(t) process.
%
% Note that, upon input, y(t) is a univariate return series observed from the 'past'.
% Upon output, y(t) represents the MSE forecast of the conditional mean of the 
% return series projected into the 'future'.
%

y  =  y((end - preSamples + 1):end , :);
e  =  e((end - preSamples + 1):end , :);
h  =  h((end - preSamples + 1):end , :).^2;  % Convert STDs to variances.

e2 =  e.^2;    % Innovations squared.

y (T,:) =  0;  % Allocate sufficient storage for forecasting.
e (T,:) =  0;
h (T,:) =  0;
e2(T,:) =  0;

%
% Construct the conditional variance parameter vector of length (1 + P + Q).
%

varianceCoefficients  =  [K ; GARCH(:) ; ARCH(:)]';  % GARCH(P,Q) parameters.

%
% Apply iterative expectations one forecast step at a time. Note that the second 
% line inside the FOR loop applies the identity that the conditional expectation
% of the squared innovation at time 't' is also the conditional expectation
% of the variance forecast at time 't' for any period in the future.
%
% Note that the forecasts constructed below make no assumptions about the 
% distribution of the innovations process e(t) other than that the conditional 
% distribution is symmetric with zero-mean (see Baillie & Bollerslev, 1992).
%

nPaths  =  size(y,2);    % # of realizations (i.e., sample paths).

for t = (preSamples + 1):T
    h (t,:) =  varianceCoefficients * [ones(1,nPaths) ; h(t-(1:P),:) ; e2(t-(1:Q),:)];
    e2(t,:) =  h(t,:);
end

%
% Now strip off the first 'preSamples' values of the conditional variance h(t). 
% After the first 'preSamples' values have been excluded, h(j) = j-period-ahead 
% forecast of the conditional variance.
%

h  =  h((preSamples + 1):end , :);

%
% Forecast the conditional mean of y(t) only if it's requested.
%

if nargout >= 2

%
%  Ensure all coefficients needed for the conditional 
%  mean exist and have proper dimensions.
%

   C  =  garchget(spec , 'C');       % Conditional mean constant.

   if isempty(C)
      error(' Conditional mean constant ''C'' must be specified.');
   end
   if isempty(AR) & (R ~= 0)
      error(' Auto-regressive ''AR'' coefficients must be specified.');
   end
   if isempty(MA) & (M ~= 0)
      error(' Moving-average ''MA'' coefficients must be specified.');
   end

%
%  Construct ARMA portion of the conditional mean 
%  parameter vector of length (1 + R + (1 + M)).
%

   armaCoefficients  =  [C ; AR(:) ; [1 ; MA(:)]]';

   if isempty(XF)

      for t = (preSamples + 1):T
          armaData =  [ones(1,nPaths) ; y(t-(1:R),:) ; e(t-(0:M),:)];
          y(t,:)   =  armaCoefficients * armaData;
      end

   else
%
%     Pre-pend 'preSamples' rows to XF for initialization.
%
      XF  =  [zeros(preSamples,size(XF,2)) ; XF];

      for t = (preSamples + 1):T
          armaData =  [ones(1,nPaths) ; y(t-(1:R),:) ; e(t-(0:M),:)];
          y(t,:)   =  armaCoefficients * armaData  +  regress * XF(t,:).';
      end

   end

   y  =  y((preSamples + 1):end , :);

end

%
% Compute the following for stationary/invertible ARMA conditional 
% mean models ONLY (i.e., no regression component):
%
%   (1) The forecast of the standard deviation of the returns that
%       would be obtained if an asset were held for multiple periods.
%       In other words, for each period 1, 2, ... horizon, this is 
%       the volatility forecast of the effective, or cumulative,
%       return that would be expected if the asset was held for 1
%       period then sold, for 2 periods then sold, and so on. For
%       conditional mean models without an ARMAX component such that
% 
%           y(t) = C + e(t)
%
%       the result is the conditional analog of the 'square-root-of-time-rule'
%       for scaling standard deviations over multi-period holding intervals.
%
%       These volatility forecasts are correct for continuously-compounded
%       returns, but only approximations for periodically-compounded returns.
%
%   (2) The root mean square error (RMSE) of the forecast of the 
%       conditional mean (i.e., the standard error of the forecast
%       of future values of y(t)) in each period. This is the square
%       root of equation (19) of Baillie & Bollerslev for each 
%       period 1, 2, ... horizon.
%

yRMSE      =  [];
yTotalRMSE =  [];

if (nargout >= 3) & isempty(X)

   InfiniteMA  =  1;
%
%  Compute the coefficients of the (truncated) infinite-order MA 
%  representation associated with the finite-order ARMA(R,M) model.
%
   if horizon > 1
      InfiniteMA  =  [InfiniteMA ; garchma(AR(:) , MA(:) , horizon - 1)];
   end

%
%  Compute the forecast of the standard deviation of the returns that
%  would be obtained if an asset were held for multiple periods. These 
%  are the volatility forecasts of returns over a multi-period holding
%  interval.
%

   yTotalRMSE  =  sqrt(toeplitz(cumsum(InfiniteMA(:)).^2 , [1  zeros(1 , horizon - 1)]) * h);

%
%  Compute the root mean square error (RMSE) of the forecast of the 
%  conditional mean in each period (i.e., on a per-period basis). This
%  is the square root of equation (19) of Baillie & Bollerslev for each 
%  period 1, 2, ... horizon.
%

   if nargout >= 4
      yRMSE  =  sqrt(toeplitz(InfiniteMA(:).^2 , [1  zeros(1 , horizon - 1)]) * h);
   end

end

%
% Since h(t) is the conditional variance, convert it to a standard deviation.
%

h  =  sqrt(h);

%
% Re-format outputs for compatibility with the return series input. 
% When return series input is specified as a single row vector, then 
% pass the outputs as row vectors. 
%

if rowY
   h          =  h(:).'; 
   y          =  y(:).';
   yRMSE      =  yRMSE(:).';
   yTotalRMSE =  yTotalRMSE(:).';
end