chaostest.m

toolbox for chaos system

function [H, pValue, Lambda, Orders, CI] = chaostest(X, ActiveFN, maxL, maxM, maxQ, alpha, StartV)
%CHAOSTEST performs the test for Chaos to test the positivity of the
%   dominant Lyapunov Exponent LAMBDA and Local Lyapunov Exponents.
%
%   The test hypothesis are:
%   Null hypothesis: LAMBDA >= 0 which indicates the presence of chaos.
%   Alternative hypothesis: LAMBDA < 0 indicates no chaos.
%   This is a one tailed test.
%
%   [H, pValue, LAMBDA, Orders, CI] = ...
%           CHAOSTEST(Series, ActiveFN, maxL, maxM, maxQ, ALPHA, STARTV)
%
% Inputs:
%   Series - a vector of observation to test.
%
% Optional inputs:
%   ActiveFN - String containing the activation function to use in the
%     neural net estimation. ActiveFN can be the 'LOGISTIC' function
%     f(u) = 1 / (1 + exp(- u)), domain = [0, 1], or 'TANH' function
%     f(u) = tanh(u), domain = [-1, 1], or 'FUNFIT' function
%     f(u) = u * (1 + |u/2|) / (2 + |u| + u^2 / 2), domain = [-1, 1].
%     Default = 'LOGISTIC'.
%
%   maxL, maxM, maxQ - the maximum orders that the chaos function defined
%     in CHAOSFN can take. Increasing the model's orders can slow down
%     calculations. Default = [5, 6, 5].
%
%   ALPHA - The significance level for the test (default = 0.05)
%
%   STARTV - String specifying the method to carry out the test, by varying
%     the triplet (L, m, q) {'VARY' or anything else} or by fixing them {'FIX'}.
%     Default = {'VARY'}.
%
% Outputs:
%   H = 0 => Do not reject the null hypothesis of Chaos at significance level ALPHA.
%   H = 1 => Reject the null hypothesis of Chaos at significance level ALPHA.
%
%   pValue - is the p-value, or the probability of observing the given
%     result by chance given that the null hypothesis is true. Small values
%     of pValue cast doubt on the validity of the null hypothesis of Chaos.
%
%  LAMBDA - The dominant Lyapunov Exponent.
%     If LAMBDA is positive, this indicates the presence of Chaos.
%     If LAMBDA is negative, this indicates the absence of Chaos.
%
%   Orders - gives the triplet (L, m, q) that minimizes the Schwartz
%     Information Criterion to obtain the best coefficients to compute LAMBDA.
%
%   CI - Confidence interval for LAMBDA at level ALPHA.
%
%   The algorithm uses the Jacobian method in contrast to the direct method
%   as descibed by Ellner et. al (1991).

%
% References: 
%             Barnett W. and He Y. (2001), "Unsolved Econometric Problem in
%               Nonlinearity, Chaos, and Bufircation", Working paper.
%             Barnett W., Gallant R., Hinich M., Jensen M., Jungeilges J.,
%               and Kalpan D. (1995), "Robustness of nonlinearity and chaos
%               tests to measurement error, inference method, and sample
%               size", Journal of Economic Behavior and Organization 27,
%               301-320.
%             Ellner S., Gallant R., McCaffrey D. and Nychka D. (1991),
%               "Convergence rates and data requirements for Jacobian-based
%               estimates of Lyapunov exponents from data", Physics Letters
%               A 153. 357-363.
%             Nychka D., Ellner S., Gallant R. and McCaffrey D. (1992),
%               "Finding Chaos in noisy systems", Journal of the Royal
%               Statistical Society B 54, 2: 399-426.
%             Nychka D., Ellner S., Balley B. (1997), "Chaos with
%               confidence: asymptotics and application of local Lyapunov
%               exponents", American Mathematical Society.
%             Whang Y. and Linton O. (1999), "The asymptotic distribution
%               of nonparametric estimates of the Lyapunov exponent for
%               stochastic time series", Journal of Econometrics 91 p 1-42.
%             Shintani M. and Linton O. (2002), "Nonparametric neural
%               network estimation of Lyapunov exponents and direct test
%               for chaos", STICERD/LSE Econometrics Discussion Paper
%               EM/02/434.
%             Andrews D. (1991), "Heteroskedasticity and autocorrelation
%               consistent covariance matrix estimation", Econometrica
%               59-3, 817-858.
%             Abarbanel H.D.I., Brown R. and Kennel M.B. (1991), "Lyapunov
%               exponents in chaotic systems: their importance and their
%               evaluation using observed data", International Journal of
%               Modern Physics B 5.
%

% Set initial conditions.

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

   if numel(X) == length(X)   % Check for a vector.
      X  =  X(:);                  % Convert to a column vector.
   else
      error(' Observation ''Series'' must be a vector.');
   end

   % Remove any NaN (i.e. missing values) observation from 'Series'. 
   
    X(isnan(X)) =   [];
         
else
    error(' Must enter at least observation vector ''Series''.');
end

%
% Specify the activation function to use in CHAOSFN and ensure it to be a
% string. Set default if necessary.
%

if nargin >= 2 && ~isempty(ActiveFN)
    if ~ischar(ActiveFN)
        error(' Activation function ''ActiveFN'' must be a string.')
    end

    % Specify the activation function to use: ActiveFN = {'logistic', 'tanh', 'funfit'}.
    
    if ~any(strcmpi(ActiveFN , {'logistic' , 'tanh' , 'funfit'}))
        error(' Activaton function ''ActiveFN'' must be LOGISTIC, TANH or FUNFIT');
    end
    
else
    ActiveFN    =   'logistic';
end

%
% Ensure the maximum orders maxL, maxM and maxQ are positive integers, and
% set default if necessary.
%

if (nargin >= 3) && ~isempty(maxL)
    if numel(maxL) > 1
      error(' Maximum lag delay ''maxL'' must be a scalar.');
    end
    if (maxL - round(maxL) ~= 0) || (maxL <= 0)
        error(' Maximum lag delay ''maxL'' must be a positive integer.');
    end
else
    maxL = 5;
end

if (nargin >= 4) && ~isempty(maxM)
    if numel(maxM) > 1
      error(' Maximum series lag ''maxM'' must be a scalar.');
    end
    if (maxM - round(maxM) ~= 0) || (maxM <= 0)
        error(' Maximum series lag ''maxM'' must be a positive integer.');
    end
else
    maxM = 6;
end

if (nargin >= 5) && ~isempty(maxQ)
    if numel(maxQ) > 1
      error(' Maximum neural-net hidden layer ''maxQ'' must be a scalar.');
    end
    if (maxQ - round(maxQ) ~= 0) || (maxQ <= 0)
        error(' Maximum neural-net hidden layer ''maxQ'' must be a positive integer.');
    end
else
    maxQ = 5;
end

%
% Ensure the significance level, ALPHA, is a 
% scalar, and set default if necessary.
%

if (nargin >= 6) && ~isempty(alpha)
   if numel(alpha) > 1
      error(' Significance level ''Alpha'' must be a scalar.');
   end
   if (alpha <= 0 || alpha >= 1)
      error(' Significance level ''Alpha'' must be between 0 and 1.'); 
   end
else
   alpha  =  0.05;
end

%
% Check the method to carry out the test, by varying (L, m, q) or by fixing
% them, and set default if necessary. The method must be a string.
%

if nargin >= 7 && ~isempty(StartV)
    if ~ischar(StartV)
        error(' Regressions type ''StartV'' must be a string.')
    end    
else
    StartV  =   'VARY';
end

if strcmpi(StartV, 'FIX')
    StartL  =   maxL;
    StartM  =   maxM;
    StartQ  =   maxQ;
else
    StartL  =   1;
    StartM  =   1;
    StartQ  =   1;
end

% Initialize the Lyapunov Exponent to a small number.
Lambda =   - Inf;

%
% Create the structure OPTIONS to use for nonlinear least square function
% LSQNONLIN (included in the optimization toolbox).
%

Options =   optimset('lsqnonlin');
Options =   optimset(Options, 'Display', 'off');

% Use the user supplied analytical Jacobian in CHAOSFN.
Options =   optimset(Options, 'Jacobian', 'on');

% Initialize the starting value for the coefficient THETA.
StartValue  =   0.5;

for L = StartL:maxL
    for m = StartM:maxM
        for q = StartQ:maxQ