Code covered by the BSD License  

Highlights from
Chaos test

Chaos test

by

 

13 Jan 2009 (Updated )

A test for chaotic dynamics of a noisy time series based on the Lyapunov exponent.

chaosfn(Theta, X, L, m, q, ActiveFN)
function [y, J] = chaosfn(Theta, X, L, m, q, ActiveFN)
%CHAOSFN is the objective neural net function used to run
%   test for chaos and to compute the Lyapunov exponent.
%   Y = CHAOSFN(THETA, X, L, m, q, ActiveFN) where X is a
%   time series and THETA is a vector of parameters.

A   =   reshape(Theta(1:q*m), m, q)'; % size(A) = [q, m].
B   =   Theta(q*m+1:q*m+q)'; % size(B) = [q, 1].
C   =   Theta(q*m+q+1:q*m+q+q); % size(C) = [1, q].
D   =   Theta(q*m+2*q+1); % size(D) = 1.
%E   =   Theta(q*m+2*q+2:q*m+2*q+1+m); % size(E) = [1, m].

XLAG   =  lagmatrix(X, L:L:L*m); % size(XLAG) = [T, m].
XLAG   =  XLAG(m*L+1:end, :); % Remove all NaN observation.
T      =  length(X) - m*L; % New sample size.

u   =   A * XLAG' + repmat(B,1,T); % size(u) = [q, T].

%
% Initialize the analytical Jacobian needed to speed up nonlinear least
% square optimization carried out by LSQNONLIN, this Jacobian is relative
% to parameter THETA: Jacobian = dF(X)/d(THETA). Not to confuse with the
% Jacobian needed to compute the Lyapunov exponent.
%

J   =   ones(length(Theta), T);

switch upper(ActiveFN)
    case 'LOGISTIC'
        % Use the logistic function as the sigmoid activation function.
        y   =   D + C * (1 ./ (1 + exp(- u))); % size(y) = [1, T].
        %y   =   E * XLAG' + D + C * (1 ./ (1 + exp(- u))); % size(y) = [1, T].        
        
        J(1:q*m,:)   =   repmat(C', m, T) .* (repmat(XLAG', q, 1) .* ...
                        repmat(exp(- u) ./ (1 + exp(- u)).^2, m, 1)); % A.

        J(q*m+1:q*m+q,:) =  repmat(C', 1, T) .* (exp(- u) ./ (1 + exp(- u)).^2); % B.
        J(q*m+q+1:q*m+q+q,:) =  1 ./ (1 + exp(- u)); % C.
        %J(q*m+2*q+2:q*m+2*q+1+m,:)  =   XLAG'; % E.
        
    case 'TANH'
        % Use a hyperbolic tangent as the activation function. This is a
        % two layered feed-forward neural network.
        y   =   D + C * tanh(u); % size(y) = [1, T].
        %y   =   E * XLAG' + D + C * tanh(u); % size(y) = [1, T].
        
        J(1:q*m,:)   =   repmat(C', m, T) .* (repmat(XLAG', q, 1) .* ...
                        repmat(sech(u).^2, m, 1)); % A.

        J(q*m+1:q*m+q,:) =  repmat(C', 1, T) .* (sech(u).^2); % B.
        J(q*m+q+1:q*m+q+q,:) =  tanh(u); % C.
        %J(q*m+2*q+2:q*m+2*q+1+m,:)  =   XLAG'; % E.
            
    case 'FUNFIT'
        % Use another type of the activation function.
        y   =   D + C * (u .* (1 + abs(u / 2)) ./ (2 + abs(u) + u.^2 / 2));
        %y   =   E * XLAG' + D + C * (u .* (1 + abs(u / 2)) ./ (2 + abs(u) + u.^2 / 2));
        
        J(1:q*m,:)   =   repmat(C', m, T) .* (repmat(XLAG', q, 1) .* ...
                repmat(8*(1 + abs(u)) ./ (3 + (1 + abs(u)).^2).^2, m, 1)); % A.
       
        J(q*m+1:q*m+q,:) =  repmat(C', 1, T) .* (8*(1 + abs(u)) ./ ...
                            (3 + (1 + abs(u)).^2).^2); % B.
        J(q*m+q+1:q*m+q+q,:) =  u .* (1 + abs(u/2)) ./ (2 + abs(u) + u.^2 / 2); % C.
        %J(q*m+2*q+2:q*m+2*q+1+m,:)  =   XLAG'; % E.
    
    otherwise
        error(' Unrecognized activation function!')
end

%
% The obtained Jacobian is for F(x), transform it to the Jacobian of
% (X - F(X)) so it can be used by LSQNONLIN to estimate THETA.
%

J   =   - J';

%
% CHAOSFN returns (X - F(X)), because the function LSQNONLIN
% minimizes the sum of squares of the objective function.
%

y   =   X(m*L+1:end)' - y;

%-------------------------------------------------------------------------%

function YLag = lagmatrix(Y,lags)
%LAGMATRIX Create matrix of lagged time series
%
% Description:
%
%   Create a matrix of lagged (time-shifted) series. Positive lags
%   correspond to delays; negative lags correspond to leads.
%
% Input Arguments:
%
%   Y - Time series data. Y may be a vector or a matrix. If Y is a vector,
%     it represents a single series. If Y is a numObs-by-numSeries matrix,
%     it represents numObs observations of numSeries series, with
%     observations across any row assumed to occur at the same time. The
%     last observation of any series is assumed to be the most recent.
%
%   lags - Vector of integer delays or leads, of length numLags, applied to
%     each series in Y. The first lag is applied to all series in Y, then
%     the second lag is applied to all series in Y, and so forth. To
%     include an unshifted copy of a series in the output, use a zero lag.
%
% Output Argument:
%
%   YLag - numObs-by-(numSeries*numLags) matrix of lagged versions of the
%     series in Y. Columns of YLag are, in order, all series in Y lagged by
%     the first lag in lags, all series in Y lagged by the second lag in
%     lags, and so forth. Unspecified observations (presample and
%     postsample data) are padded with NaN values.
%

% Check for a vector:

if numel(Y) == length(Y)
   Y = Y(:); % Ensure a column vector
end

% Ensure lags is a vector of integers:

lags = lags(:); % Ensure a column vector

% Cycle through the lags vector and shift the input time series. Positive 
% lags are delays, and can be processed by FILTER. Negative lags are leads,
% and series are flipped (reflected in time), run through FILTER, and then
% flipped again. Series with zero lags are simply copied.

numLags = length(lags); % Number of lags to apply to each time series

[numObs,numSeries] = size(Y);

YLag = NaN(ones(numObs,numSeries*numLags)); % Preallocate

for i = 1:numLags

    L       = lags(i);
    columns = (numSeries*(i-1)+1):i*numSeries; % Columns to fill, this lag

    if L > 0 % Time delays

       YLag((L + 1):end,columns) = Y(1:(end - L), :);

    elseif L < 0 % Time leads

       YLag(1:(end + L),columns) = Y((1 - L):end, :);

    else % No shifts

       YLag(:,columns) = Y;

    end

end

Contact us