function [innovations , sigma , LLF] = garchinfer(spec , y , X)
%GARCHINFER Infer GARCH innovations process from an observed return series.
%   Given a conditional mean specification of ARMAX form and conditional
%   variance specification of GARCH form, infer the innovations and conditional
%   standard deviations from an observed univariate return series. Since 
%   GARCHINFER is a wrapper around the appropriate log-likelihood objective 
%   function, the log-likelihood value is also computed for convenience.
%
%   [Innovations , Sigma , LogLikelihood] = garchinfer(Spec , Series)
%   [Innovations , Sigma , LogLikelihood] = garchinfer(Spec , Series, X)
%
%   Optional Input: X
%
% Inputs:
%   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.
%
%   Series - Matrix of observations of the underlying univariate return series 
%     of interest. Series is the response variable representing the time 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.
%
% Optional Input:
%   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 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.
%
% Outputs:
%   Innovations - Innovations matrix inferred from the input Series matrix. The
%     size of Innovations is the same as the size of Series.
%
%   Sigma - Conditional standard deviation matrix corresponding to Innovations.
%     The size of Sigma is the same as the size of Series.
%
%   LogLikelihood - Vector of log-likelihood objective function values for each 
%     realization of Series. The length of LogLikelihood is the same as the 
%     number of columns in Series.
%
% See also GARCHSET, GARCHSIM, GARCHPRED, GARCHFIT, GARCHLLFN.

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

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

%
% 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). During estimation, the innovations
% process e(t) must be inferred from the conditional mean specification, which 
% may include a regression component if desired. In contrast to simulation, 
% estimation of the innovations process e(t) is NOT independent of X. 
%

if (nargin >= 3) & ~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

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

%
% Format the parameter estimation vector.
%
% NOTES:
% 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)
%
% The input coefficient vector 'Parameters' is formatted exactly as the 
% coefficients would be read from the recursive difference equations when 
% solving for the current values of the y(t) and h(t) time series. 
%
% For example, consider the following equations for the conditional mean 
% and variance of an ARMAX(R=2,M=2,nX=1)/GARCH(P=2,Q=2) composite model:
%
%   y(t) =  1.3 + 0.5y(t-1) - 0.8y(t-2) + e(t) - 0.6e(t-1) + 0.08e(t-2) + 1.2X(t)
%   h(t) =  0.5 + 0.2h(t-1) + 0.1h(t-2) + 0.3e(t-1)^2  +  0.4e(t-2)^2
%
% The above equations are examples of the following general 
% ARMAX(R,M,nX)/GARCH(P,Q) form:
%
%   y(t) =  C + AR(1)y(t-1) + ... + AR(R)y(t-R) + e(t) 
%             + MA(1)e(t-1) + ... + MA(M)e(t-M) + B(1)X(t,1) + ... + B(nX)X(t,nX)
%
%   h(t) =  K + GARCH(1)h(t-1)   + ... + GARCH(P)h(t-P) 
%             +  ARCH(1)e(t-1)^2 + ... +  ARCH(Q)e(t-Q)^2
%
% For the example listed above, the coefficient vector 'Parameters' would be
%
%   Parameters = [ C     AR(1:R)     MA(1:M)  B(1:nX)  K   GARCH(1:P)   ARCH(1:Q)]'
%              = [1.3   0.5 -0.8   -0.6 0.08    1.2   0.5    0.2 0.1     0.3 0.4 ]'
%
% Notice that the coefficient of e(t) in the conditional mean equation is
% defined to be 1, and is NOT included in 'Parameters' vector because it
% is not estimated.
%

Parameters = [ C ; AR(:) ; MA(:) ; regress(:) ; K ; GARCH(:) ; ARCH(:)];

%
% Get the probability distribution of the innovations process e(t) 
% and call the appropriate log-likelihood objective function.
%

distribution  =  garchget(spec , 'Distribution');
distribution  =  distribution(~isspace(distribution));

switch upper(distribution)
   case 'GAUSSIAN'
      [LLF,G,H,innovations,sigma] = garchllfn(Parameters , y , R , M , P , Q , X);
   otherwise
      error(' Distribution of innovations must be ''Gaussian''.')
end

%
% Negate objective function value to compensate for the minimization function FMINCON.
%

LLF  =  -LLF;

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

if rowY
   innovations  =  innovations(:).';
   sigma        =  sigma(:).';
end