%
% 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 would be formatted as
%
%   Coefficient = [ 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 the vector because it is not estimated.
%

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

      x0  =  [C ; AR(:) ; MA(:) ; regress(:) ; K ; GARCH(:) ; ARCH(:)];   % Initial guess.
%
%     Set lower bounds constraints. 
%
%     The parameters of the conditional mean, in theory, have no lower bounds. 
%     However, for robustness, set the mean constant (C) and the ARMA(R,M) 
%     parameters to -100 (i.e., an unrealistically large, yet totally finite 
%     number). The coefficients of the regression component of the mean are 
%     set to -infinity to reflect the fact that the origin and generating 
%     mechanism of a regression is completely unknown. This is consistent 
%     with the 'hands-off' approach regarding regression.
%
%     The conditional variance constant (K) must be positive, so it's set to 
%     a very small (yet positive) value; the conditional variance coefficients 
%     (GARCH(i), ARCH(i)) are constrained to be non-negative, so they're set 
%     to zero.
%
      infinity     =  inf;
      hundred      =  100;
      lowerBounds  =  [-hundred(ones(1+R+M,1)) ; -infinity(ones(size(X,2),1)) ; 1e-10 ; zeros(P+Q,1)];
      upperBounds  =  [];

%
%     Set linear inequality of the covariance-stationarity constraint of 
%     the conditional variance, Ax <= b. Since the conditional variance 
%     (GARCH(i), ARCH(i)) parameters are constrained to be non-negative, 
%     the covariance-stationarity constraint is just a summation constraint.
%     Also, adjust the summation constraint, b, to reflect a tolerance 
%     offset from a fully integrated conditional variance condition (i.e.,
%     an IGARCH process).
%

      if ((P + Q) > 0)
         A  =  [zeros(1,1+R+M+size(X,2)+1)  ones(1,P)  ones(1,Q)];
         b  =  1  -  2*optimget(spec.Optimization , 'TolCon', 1e-6);
      else
         A  =  [];
         b  =  [];
      end
%
%     Set any linear equality constraints.
%
      if any(Fix)
         i    =  find(Fix); 
         Aeq  =  [zeros(length(i) , 1 + R + M + size(X,2) + 1 + P + Q)];

         for j = 1:length(i)
             Aeq(j,i(j))  =  1;
         end

         beq  =  x0(logical(Fix));

      else
         Aeq  =  [];
         beq  =  [];
      end
%
%     Perform constrained non-linear optimization. 
%

      [coefficients , LLF    , ...
       exitFlag     , output , lambda] =  fmincon('garchllfn'  , x0 , A  , b       , Aeq , beq , ...
                                                   lowerBounds , upperBounds       , ...
                                                  'garchnlc'   , spec.Optimization , ...
                                                   y , R , M , P , Q , X);
%
%     Negate objective function value to compensate for FMINCON.
%
      LLF  =  -LLF;
%
%     Over-write all GARCH constraint-violating parameters that are less than 
%     zero. This will, occasionally, occur because FMINCON may violate constraints 
%     ever so slightly. Also, the constant term (K) of the conditional variance 
%     equation must be positive, so enforce this if necessary.
%
      varianceCoefficients =  coefficients((2 + R + M + size(X,2)):end);

      varianceCoefficients(varianceCoefficients <  0) =  0;              % All variance-related coefficients.
      varianceCoefficients(1) = max(varianceCoefficients(1) , realmin);  % Constant term 'K'.

      coefficients  =  [coefficients(1:(1 + R + M + size(X,2))) ; varianceCoefficients(:)];

      MLEparameters =  coefficients(:);  % Save MLE parameter values in vector form.
%
%     Extract the parameter estimates & equality constraints.
%
      C          =  coefficients(1);
      AR         =  coefficients(2:R+1);
      MA         =  coefficients(R+2:R+M+1);
      regress    =  coefficients(R+M+2:R+M+size(X,2)+1);

      K          =  coefficients(R+M+size(X,2)+2:R+M+size(X,2)+2);
      GARCH      =  coefficients(R+M+size(X,2)+3:R+M+size(X,2)+2+P);
      ARCH       =  coefficients(R+M+size(X,2)+3+P:R+M+size(X,2)+2+P+Q);

      FixC       =  Fix(1);
      FixAR      =  Fix(2:R+1);
      FixMA      =  Fix(R+2:R+M+1);
      FixRegress =  Fix(R+M+2:R+M+size(X,2)+1);

      FixK       =  Fix(R+M+size(X,2)+2:R+M+size(X,2)+2);
      FixGARCH   =  Fix(R+M+size(X,2)+3:R+M+size(X,2)+2+P);
      FixARCH    =  Fix(R+M+size(X,2)+3+P:R+M+size(X,2)+2+P+Q);

%
%     Test for stationarity & invertibility of the ARMA model (if any).
%
      summary.warning = 'No Warnings';

      if any(abs(roots([1 ; -AR(:)])) >= 1)
         summary.warning = 'ARMA Model is Not Stationary/Invertible';
         if DisplayFlag
            warning('Auto-Regressive Polynomial is Non-Stationary.')
         end
      end

      if any(abs(roots([1 ; MA(:)])) >= 1)
         summary.warning = 'ARMA Model is Not Stationary/Invertible';
         if DisplayFlag
            warning('Moving-Average Polynomial is Non-Invertible.')
         end
      end
%
%     Now pack the data into the output COEFFICIENTS structure. Note that 
%     the COEFFICIENTS output structure is of the same form as the SPEC 
%     input structure. This allows GARCHSET, GARCHGET, GARCHSIM, and 
%     GARCHPRED to accept either SPEC or COEFFICIENTS seamlessly. 
%
%     Strictly speaking, GARCHSET should be used to make the following
%     assignment, but will error-out if the stationarity/invertibility
%     constraints are violated. Simple assignment allows for graceful 
%     termination with a warning.
%
      coefficients            =  spec;
      coefficients.C          =  C;
      coefficients.AR         =  AR(:)';
      coefficients.MA         =  MA(:)';
      coefficients.Regress    =  regress(:)';
      coefficients.K          =  K;
      coefficients.GARCH      =  GARCH(:)';
      coefficients.ARCH       =  ARCH(:)';

      coefficients.FixC       =  FixC;
      coefficients.FixAR      =  FixAR(:)';
      coefficients.FixMA      =  FixMA(:)';
      coefficients.FixRegress =  FixRegress(:)';
      coefficients.FixK       =  FixK;
      coefficients.FixGARCH   =  FixGARCH(:)';
      coefficients.FixARCH    =  FixARCH(:)';
%
%     Update the 'comment' field to reflect a regression component (only if auto-generated).
%
      comment =  garchget(coefficients , 'comment');

      if length(findstr(comment, char(0))) == 2
         pOpen  =  findstr(comment, '(');
         pClose =  findstr(comment, ')');
         if ~isempty(pOpen) & ~isempty(pClose) 
            commas  =  findstr(comment(pOpen(1):pClose(1)) , ',');
            if length(commas) == 1
               coefficients.Comment =  [comment(1:pClose(1)-1) ',' num2str(size(X,2)) comment(pClose(1):end)];
            elseif length(commas) == 2
               coefficients.Comment =  [comment(1:(pOpen(1) + commas(2)-1)) num2str(size(X,2)) comment(pClose(1):end)];
            end
        end
      end

      if length(coefficients.AR)      == 0 , coefficients.AR          =  []; end   % Just for aesthetics.
      if length(coefficients.MA)      == 0 , coefficients.MA          =  []; end
      if length(coefficients.Regress) == 0 , coefficients.Regress     =  []; end
      if length(coefficients.GARCH)   == 0 , coefficients.GARCH       =  []; end
      if length(coefficients.ARCH)    == 0 , coefficients.ARCH        =  []; end

      if sum(coefficients.FixC)       == 0 , coefficients.FixC        =  []; end   % Just for aesthetics.
      if sum(coefficients.FixAR)      == 0 , coefficients.FixAR       =  []; end
      if sum(coefficients.FixMA)      == 0 , coefficients.FixMA       =  []; end
      if sum(coefficients.FixRegress) == 0 , coefficients.FixRegress  =  []; end
      if sum(coefficients.FixK)       == 0 , coefficients.FixK        =  []; end
      if sum(coefficients.FixGARCH)   == 0 , coefficients.FixGARCH    =  []; end
      if sum(coefficients.FixARCH)    == 0 , coefficients.FixARCH     =  []; end

%
%     Compute the variance-covariance matrix of the parameter 
%     estimates and extract the standard errors of the estimation.
%
      if (nargout >= 2) | (nargout == 0)

         covarianceMatrix =  varcov(MLEparameters , y , R , M , P , Q , X , Fix);
         standardErrors   =  sqrt(diag(covarianceMatrix))';

         errors.C         =  standardErrors(1);
         errors.AR        =  standardErrors(2:R+1);
         errors.MA        =  standardErrors(R+2:R+M+1);
         errors.Regress   =  standardErrors(R+M+2:R+M+size(X,2)+1);

         errors.K         =  standardErrors(R+M+size(X,2)+2:R+M+size(X,2)+2);
         errors.GARCH     =  standardErrors(R+M+size(X,2)+3:R+M+size(X,2)+2+P);
         errors.ARCH      =  standardErrors(R+M+size(X,2)+3+P:R+M+size(X,2)+2+P+Q);

         if length(errors.AR)      == 0 , errors.AR       =  []; end   % Just for aesthetics.
         if length(errors.MA)      == 0 , errors.MA       =  []; end
         if length(errors.Regress) == 0 , errors.Regress  =  []; end
         if length(errors.GARCH)   == 0 , errors.GARCH    =  []; end
         if length(errors.ARCH)    == 0 , errors.ARCH     =  []; end

      end
%
%     Return the innovations and conditional standard deviation vectors if requested.
%
      if (nargout >= 4) | (nargout == 0)
         [innovations , sigma]  =  garchinfer(coefficients , y , X);
      end
%
%     Return summary information if requested.
%
      if (nargout >= 6) | (nargout == 0)

         if exitFlag == 0
            summary.converge  =  'Maximum Function Evaluations or Iterations Reached';
         elseif exitFlag < 0
            summary.converge  =  'Function Did NOT Converge';
         elseif exitFlag > 0
            summary.converge  =  'Function Converged to a Solution';
         end
 
         summary.covMatrix      =  covarianceMatrix;
         summary.iterations     =  output.iterations;
         summary.functionCalls  =  output.funcCount;
%
%        Flag any boundary constraints enforced EXCEPT LINEAR EQUALITY CONSTRAINTS
%        specifically requested by the user. Whenever any constraints (excluding
%        linear equalities) are imposed, the log-likelihood function will probably NOT
%        be approximately quadratic at the solution. In this case, the standard errors
%        of the parameter estimates are unlikely to be accurate. However, if the ONLY
%        constraints are the linear equalities imposed by the user, then the resulting
%        log-likelihood value LLF may still be useful for post-fit assessment and
%        inference tests, such as likelihood ratio tests (LRT's).
%

         TolCon = optimget(spec.Optimization , 'TolCon', 1e-6);

         if (norm([lambda.lower(:) ; lambda.upper(:) ; lambda.ineqlin(:) ; lambda.ineqnonlin(:)] , 1) > TolCon)
            summary.constraints  =  'Boundary Constraints Active; Errors may be Inaccurate';
            if DisplayFlag
               warning('Boundary Constraints Active; Standard Errors may be Inaccurate.')
            end
         else
            summary.constraints  =  'No Boundary Constraints';
         end

      end

   otherwise

      error(' Distribution of innovations must be ''Gaussian''.')

end

%
% 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 & (nargout >= 4)
   innovations  =  innovations(:).';
   sigma        =  sigma(:).';
end

%
% Perform the default no-output action: 
%
%  (1) Print the parameter estimates to the screen, and 
%  (2) Display the estimated residuals, conditional standard
%      deviations, and input raw return series.
%

if nargout == 0

   garchdisp(coefficients , errors);
   garchplot(innovations  , sigma , y);

   disp(' ')
   fprintf('  Log Likelihood Value: %f\n\n' , LLF)

   clear coefficients     % Suppress unexpected printing to the command window.

end

%
%   * * * * *  Helper function for initial GARCH guesses.  * * * * *
%

function [K , GARCH , ARCH] = garch0(P , Q , unconditionalVariance)
%GARCH0 Initial GARCH process parameter estimates.
%   Given the orders of a GARCH(P,Q) model and an estimate of the unconditional 
%   variance of the innovations process, compute initial estimates for the 
%   (1 + P + Q) parameters of a GARCH(P,Q) conditional variance model. These
%   estimates serve as initial guesses for further refinement via maximum 
%   likelihood.
%
%   [K , GARCH , ARCH] = garch0(P , Q , Variance)
%
% Inputs:
%   P - Non-negative, scalar integer representing the number of lags of the
%     conditional variance included in the GARCH process.
%
%   Q - Non-negative, scalar integer representing the number of lags of the 
%     squared innovations included in the GARCH process.
%
%   Variance - Estimate of the unconditional variance of the innovations noise
%     process. Equivalently, it may also be viewed as the variance estimate of
%     the white noise innovations under the assumption of homoskedasticity.
%
% Outputs:
%   K - Conditional variance constant (scalar).
%
%   GARCH - P-element column vector of coefficients of lagged conditional 
%     variances.
%
%   ARCH - Q-element column vector of coefficients of lagged squared 
%     innovations.
%
% Note: 
%   This is an internal helper function for GARCHFIT. No error checking is 
%   performed.
%

%
% The following initial guesses are based on empirical observation of GARCH
% model parameters. This approach is rather ad hoc, but is very typical of
% GARCH models in financial time series. The most common (and very useful!) 
% GARCH model is the simple GARCH(1,1) model in which the coefficient of the
% lagged conditional variance (i.e., the 'GARCH' coefficient) is about 0.8 
% to 0.9, and the coefficient of lagged squared innovation (i.e., the 'ARCH' 
% coefficient) is about 0.05. Thus, a reasonable GARCH(1,1) assumption is:
%
%            h(t) = K + 0.85h(t-1) + 0.05e^2(t-1)
%
% In a GARCH(1,1) model, the unconditional variance of the innovations 
% process, V, is
%
%            V = K / (1 - (0.85 + 0.05))
% or,
%
%            K = V*(1 - (0.85 + 0.05))
%
% For higher-order GARCH(P,Q) models, this approach assumes the sum of the
% lagged conditional variance coefficient = 0.85, and the sum of the coefficients
% of lagged squared innovations = 0.05.
%

GARCH =  0.85;
GARCH =  GARCH(ones(P,1)) / max(P,1);
GARCH =  GARCH(:);

ARCH  =  0.05;
ARCH  =  ARCH(ones(Q,1)) / max(Q,1);