function [H, pValue, Lambda, Orders, CI, allLambda] = chaostest(X, ActiveFN, maxOrders, alpha, flag)
%CHAOSTEST performs the test for Chaos to test the positivity of the
% dominant Lyapunov Exponent LAMBDA.
%
% 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, AllLambda] = ...
% CHAOSTEST(Series, ActiveFN, maxOrders, ALPHA, FLAG)
%
% 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 'SIGMOID' function
% f(u) = u * (1 + |u/2|) / (2 + |u| + u^2 / 2), domain = [-1, 1].
% Default = 'TANH'.
%
% maxOrders - the maximum orders that the chaos function defined
% in CHAOSFN can take. This must be a vector containing 3 elements.
% maxOrders = [L, m, q].
% Increasing the model's orders can slow down calculations.
% Default = [5, 6, 5].
%
% ALPHA - The significance level for the test (default = 0.05)
%
% FLAG - 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 maximizes the Lyapunov
% exponent computed from all L*m*q estimations.
%
% CI - Confidence interval for LAMBDA at level ALPHA.
%
% allLambda - all the computed Lyapunov exponents for all L-by-m-by-q
% estimations.
%
% The algorithm uses the Jacobian method in contrast to the direct
% method, it needs the optimazition and the statistics toolboxes.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 17 March 2009 by Ahmed BenSa�da %
% Department of Finance, IHEC Sousse - Tunisia %
% Email: ahmedbensaida@yahoo.com %
% $ Revision 5.1 $ Date: 2 April 2016 $ %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% References:
% BenSa�da, A. and Litimi, H. (2013), "High level chaos in the
% exchange and index markets", Chaos, Solitons & Fractals 54:
% 90-95.
% BenSa�da, A. (2014), "Noisy chaos in intraday financial
% data: Evidence from the American index", Applied
% Mathematics and Computation 226: 258-265.
% BenSa�da, A. (2015), "A practical test for noisy chaotic
% dynamics", SoftwareX 3-4:1-5.
%
% Check if the Optimization and Statistics toolboxes are installed.
% 1 if installed, 0 if not.
existToolbox = license('test','optimization_toolbox') && ...
license('test','statistics_toolbox');
if ~existToolbox
error('The ''CHAOSTEST'' function needs both Optimization and Statistics toolboxes.')
end
% Set initial conditions.
if (nargin >= 1) && (~isempty(X))
if istable(X) % Convert Table to an array.
X = table2array(X);
end
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 , {'tanh' , 'logistic' , 'sigmoid'}))
error('Activaton function ''ActiveFN'' must be TANH, LOGISTIC, or SIGMOID');
end
else
ActiveFN = 'tanh';
end
%
% Ensure the maximum orders maxL, maxM and maxQ are positive integers, and
% set default if necessary.
%
if (nargin >= 3) && ~isempty(maxOrders)
if numel(maxOrders) ~= 3 || any((maxOrders - round(maxOrders) ~= 0)) || ...
any((maxOrders <= 0))
error(' Maximum orders ''maxOrders'' must a vector of 3 positive integers.');
end
else
maxOrders = [5, 6, 5];
end
maxL = maxOrders(1);
maxM = maxOrders(2);
maxQ = maxOrders(3);
% Check for a minimum size of vector X.
if length(X) <= maxL * maxM
error(' Observation ''Series'' must have at least %d obeservations.',maxL*maxM+1)
end
%
% Ensure the significance level, ALPHA, is a
% scalar, and set default if necessary.
%
if (nargin >= 4) && ~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 >= 5 && ~isempty(flag)
if ~ischar(flag)
error(' Regressions type ''FLAG'' must be a string.')
end
else
flag = 'VARY';
end
if strcmpi(flag, '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;
Orders = [maxL, maxM, maxQ];
% Initialized all value of Lambdas.
allLambda = -Inf(maxL, maxM, maxQ);
%
% 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
A = StartValue * ones(q, m); % Contains the coefficients gamma(i,j).
B = StartValue * ones(q, 1); % Contains the coefficients gamma(0,j).
C = StartValue * ones(1, q); % Contains the coefficients beta(j).
D = StartValue; % Constant.
%E = StartValue * ones(1, m); % Linear coefficients for lagged X.
Theta0 = [reshape(A', 1, q*m), B', C, D];
%Theta0 = [reshape(A', 1, q*m), B', C, D, E];
Theta = lsqnonlin(@chaosfn, Theta0, [], [], Options, X, L, m, q, ActiveFN);
T = length(X) - m*L; % New sample size.
%
% Use an inner function JACOBMATX to compute the Jacobian needed
% to compute the Lyapunov Exponent LAMBDA. This jacobian is relative
% to X and not to parameter THETA. N.B: the jacobian gi
