Lyapunov Exponent Python Implementation

Question

I have 10k data points like this:

0.010222
0.010345
0.010465
0.010611
0.010768
0.010890
0.011049
0.011206
0.011329
0.011465
0.011613
0.11763
0.011888
0.012015
0.012154
0.012282
0.012408
0.012524
....

I want to calculate Lyapunov exponent for that. This is what I've done so far:

lyapunovs = []
eps = 0.0001
for i in range(N):
    for j in range(i + 1, N):
        if np.abs(data[i] - data[j]) < eps:
            for k in range(1, min(N - i, N - j)):
                 d0 = np.abs(data[i] - data[j])
                 dn = np.abs(data[i + k] - data[j + k])
                 lyapunovs.append(math.log(dn) - math.log(d0))  # problem

My problem is that I don't know first Lyapunov exponent is average of all the lyapunovs when k = 1 or average of all the lyapunovs for the first time that data[i] - data[j] < eps?

Is this right implementation for Lyapunov exponent?

And this is the Numerical Calculation of Lyapunov Exponent

Answer 1

I see from the chosen loop structure in the question that a triangle of the Cartesian product of the points is being used. This might improve the estimate of the derivatives, which are susceptible to noise, but it is not part of the Lyapunov exponent explicitly. See this example of the calculations on a known function in the absence of measurement error. Feel free to look into that aspect more, but below I will assume the comparison of signal points adjacent in time.

Your original question uses NumPy, so I will also make use of it. One of the rules of thumb to using NumPy well is to avoid loops, although it is possible to vectorize functions that contain loops. With no explicit time measurements, and no repeated values, you could simply do:

import numpy as np
x = np.random.normal(0,1,size=10**4) # Mock signal data
np.mean(np.log(np.abs(np.diff(x))))

Or if the signal is paired with an array of timepoints, then the numerical derivative can involve time:

import numpy as np
x = np.random.normal(0,1,size=10**4) # Mock signal data
t = np.arange(10**4) # Mock time data
np.mean(np.log(np.abs(np.diff(x) / np.diff(t))))

However, in some datasets it is possible for adjacent values to repeat! This can occur when you've measured the signal only to a few decimal places, and it is a problem because it leads to np.log(0) (=-np.inf) which will blow up your calculation. A simple solution is to remove duplicated values, but this will only be suitable if duplicates are relatively rare and you have a large sample size. It is possible to estimate an upper bound on the estimate of the L-exponent by considering the precision of your measurements, but that is not the estimate of the L-exponent itself.