function [e , h , y] = garchsim(spec, nSamples, nPaths, seed, X)
%GARCHSIM Univariate GARCH process simulation.
%   Given specifications for the conditional mean and variance of a univariate 
%   time series, simulate sample paths for the return series, innovations, and 
%   conditional standard deviation processes. Each of NUMPATHS sample paths are 
%   sampled at NUMSAMPLES observations. The conditional mean may be of general 
%   ARMAX form and conditional variance of general GARCH form.
%
%   [Innovations, Sigma, Series] = garchsim(Spec)
%   [Innovations, Sigma, Series] = garchsim(Spec, NumSamples, NumPaths, Seed, X)
%
%   Optional Inputs: NumSamples, NumPaths, Seed, X
%
% Input:
%   Spec - Structure specification for the conditional mean and variance 
%     models. Spec is a MATLAB structure with fields generated by calling the 
%     function GARCHSET. Type "help garchset" for details.
%
% Optional Inputs:
%   NumSamples - Positive integer indicating the number of samples generated for 
%     each path of Innovations, Sigma, and Series outputs. If empty or missing, 
%     the default is 100.
%
%   NumPaths - Positive integer indicating the number of sample paths 
%     (realizations) generated for the Innovations, Sigma, and Series outputs.
%     If empty or missing, the default is 1.
%
%   Seed - Scalar random number generator seed. If empty or missing, default 
%     is 0 (the MATLAB initial state).
%
%   X - Time series regression matrix of explanatory variable(s). Typically, X 
%     is a regression matrix of asset returns (e.g., the return series of an 
%     equity index). Each column of X is an individual time series used as an 
%     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 empty or missing, the conditional mean 
%     will have no regression component. If specified, then at least the most 
%     recent NUMSAMPLES observations of each return series must be valid (i.e.,
%     non-NaN). When the number of valid observations in each series exceeds 
%     NUMSAMPLES, only the most recent NUMSAMPLES observations of X are used.
%
% Outputs:
%   Innovations - NUMSAMPLES by NUMPATHS matrix of innovations, representing a 
%     mean zero, discrete-time stochastic process. The Innovations time series 
%     follows the input conditional variance (GARCH) specification. Rows are 
%     sequential times samples, columns are independent realizations.
%
%   Sigma - NUMSAMPLES by NUMPATHS matrix of conditional standard deviations 
%     of the corresponding Innovations matrix. Innovations and Sigma are the 
%     same size, and form a matching pair of matrices. Rows are sequential 
%     times samples, columns are independent realizations.
%
%   Series - NUMSAMPLES by NUMPATHS matrix of the return series of interest. 
%     Series is the dependent stochastic process and follows the input 
%     conditional mean specification of general ARMAX form. Rows are 
%     sequential times samples, columns are independent realizations.
%
% See also GARCHSET, GARCHGET, GARCHPRED, GARCHFIT.

%
% References:
%   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.
%

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

%
% Check input parameters and set defaults.
%

if (nargin >= 2) & ~isempty(nSamples)
   if prod(size(nSamples)) > 1
      error(' Number of observations ''NumSamples'' must be a scalar.');
   end
   if (round(nSamples) ~= nSamples) | (nSamples <= 0)
      error(' Number of observations ''NumSamples'' must be a positive integer.');
   end
else
   nSamples  =  100;   % Set default.
end

if (nargin >= 3) & ~isempty(nPaths)
   if prod(size(nPaths)) > 1
      error(' Number of sample paths ''NumPaths'' must be a scalar.');
   end
   if (round(nPaths) ~= nPaths) | (nPaths <= 0)
      error(' Number of sample paths ''NumPaths'' must be a positive integer.');
   end
else
   nPaths  =  1;       % Set default.
end

if (nargin >= 4) & ~isempty(seed)
   if prod(size(seed)) > 1
      error(' Random number generator seed ''Seed'' must be a scalar.');
   end
else
   seed  =  0;         % Set default.
end

%
% Scrub the regression matrix and ensure sufficient observations exist. 
% Note that, in contrast to estimation, simulation of the innovations 
% process is independent of X. However, during estimation, the innovations
% process must be inferred from the conditional mean specification, which
% includes a regression component if desired.
%

if (nargin >= 5) & ~isempty(X) & (nargout >= 3)

   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) < nSamples
      error(' Regression matrix ''X'' has insufficient number of observations.');
   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

%
% Ensure all coefficients exist and have proper dimensions.
%

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.

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

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(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
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

%
% ARMA(R,M)/GARCH(P,Q) processing requires pre-sample values for conditioning. 
% 
% 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)
%
% 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, create max([R M P Q]) pre-sample lags
% for all processes in the ARMA(R,M)/GARCH(P,Q) model.
%