Nonlinear Timeseries Analysis

Numerical Lyapunov Exponent

Given any timeseries, one can first embed it using delay coordinates, and then calculate a maximum Lyapunov exponent for it. This is done with

numericallyapunov(R::Dataset, ks;  refstates, w, distance, ntype)

Return E = [E(k) for k ∈ ks], where E(k) is the average logarithmic distance between states of a neighborhood that are evolved in time for k steps (k must be integer). Typically R is the result of delay coordinates of a single timeseries.

Keyword Arguments

  • refstates = 1:(length(R) - ks[end]) : Vector of indices that notes which states of the reconstruction should be used as "reference states", which means that the algorithm is applied for all state indices contained in refstates.
  • w::Int = 1 : The Theiler window.
  • ntype = NeighborNumber(1) : The neighborhood type. Either NeighborNumber or WithinRange. See Neighborhoods for more info.
  • distance::Metric = Cityblock() : The distance function used in the logarithmic distance of nearby states. The allowed distances are Cityblock() and Euclidean(). See below for more info. The metric for finding neighbors is always the Euclidean one.


If the dataset exhibits exponential divergence of nearby states, then it should hold

\[E(k) \approx \lambda\cdot k \cdot \Delta t + E(0)\]

for a well defined region in the k axis, where $\lambda$ is the approximated maximum Lyapunov exponent. $\Delta t$ is the time between samples in the original timeseries. You can use linear_region with arguments (ks .* Δt, E) to identify the slope (= $\lambda$) immediatelly, assuming you have choosen sufficiently good ks such that the linear scaling region is bigger than the saturated region.

The algorithm used in this function is due to Parlitz[Skokos2016], which itself expands upon Kantz [Kantz1994]. In sort, for each reference state a neighborhood is evaluated. Then, for each point in this neighborhood, the logarithmic distance between reference state and neighborhood state(s) is calculated as the "time" index k increases. The average of the above over all neighborhood states over all reference states is the returned result.

If the Metric is Euclidean() then use the Euclidean distance of the full D-dimensional points (distance $d_E$ in ref.[Skokos2016]). If however the Metric is Cityblock(), calculate the absolute distance of only the first elements of the m+k and n+k points of R (distance $d_F$ in ref.[Skokos2016], useful when R comes from delay embedding).

The function numericallyapunov has a total of 4 different approaches for the algorithmic process, by combining 2 types of distances with 2 types of neighborhoods.

Example of Numerical Lyapunov computation

using DynamicalSystems, PyPlot

ds = Systems.henon()
data = trajectory(ds, 100000)
x = data[:, 1] #fake measurements for the win!

ks = 1:20
ℜ = 1:10000
fig = figure(figsize=(10,6))

for (i, di) in enumerate([Euclidean(), Cityblock()])
    subplot(1, 2, i)
    ntype = NeighborNumber(2)
    title("Distance: $(di)", size = 18)
    for D in [2, 4, 7]
        R = embed(x, D, 1)
        E = numericallyapunov(R, ks;
        refstates = ℜ, distance = di, ntype = ntype)
        Δt = 1
        λ = linear_region(ks.*Δt, E)[2]
        # gives the linear slope, i.e. the Lyapunov exponent
        plot(ks .- 1, E .- E[1], label = "D=$D, λ=$(round(λ, digits = 3))")

Bad Time-axis (ks) length

Large `ks`

This simply cannot be stressed enough! It is just too easy to overshoot the range at which the exponential expansion region is valid!

Let's revisit the example of the previous section:

ds = Systems.henon()
data = trajectory(ds, 100000)
x = data[:, 1]

The timeseries of such length could be considered big. A time length of 100 seems very small. Yet it turns out it is way too big! The following

ks = 1:100
R = embed(x, 2, 1)
E = numericallyapunov(R, ks, ntype = NeighborNumber(2))
fig = figure()
plot(ks .- 1, E .- E[1])
title("Lyappunov: $(linear_region(ks, E)[2])")
┌ Warning: Found linear region spans less than a 3rd of the available x-axis and might imply inaccurate slope or insufficient data. Recommended: plot `x` vs `y`.
└ @ ChaosTools ~/.julia/packages/ChaosTools/Ax8S6/src/dimensions/linear_regions.jl:131

Notice that even though this value for the Lyapunov exponent is correct, it happened to be correct simply due to the jitter of the saturated region. Since the saturated region is much bigger than the linear scaling region, if it wasn't that jittery the function linear_region would not give the scaling of the linear region, but instead a slope near 0! (or if you were to give bigger tolerance as a keyword argument)

Case of a Continuous system

The process for continuous systems works identically with discrete, but one must be a bit more thoughtful when choosing parameters. The following example helps the users get familiar with the process:

using DynamicalSystems, PyPlot

ds = Systems.lorenz()
# create a timeseries of 1 dimension
dt = 0.05
x = trajectory(ds, 1000.0; dt = dt)[:, 1]
20001-element Array{Float64,1}:

We know that we have to use much bigger ks than 1:20, because this is a continuous case! (See reference given in numericallyapunovspectrum)

ks1 = 0:200

and in fact it is even better to not increment the ks one by one but instead do

ks2 = 0:4:200

Now we plot some example computations

ntype = NeighborNumber(5) #5 nearest neighbors of each state

for d in [4, 8], τ in [7, 15]
    r = embed(x, d, τ)

    # E1 = numericallyapunov(r, ks1; ntype)
    # λ1 = linear_region(ks1 .* dt, E1)[2]
    # plot(ks1,E1.-E1[1], label = "dense, d=$(d), τ=$(τ), λ=$(round(λ1, 3))")

    E2 = numericallyapunov(r, ks2; ntype)
    λ2 = linear_region(ks2 .* dt, E2)[2]
    plot(ks2,E2.-E2[1], label = "d=$(d), τ=$(τ), λ=$(round(λ2, digits = 3))")

xlabel("k (0.05×t)")
ylabel("E - E(0)")
title("Continuous Reconstruction Lyapunov")
┌ Warning: Found linear region spans less than a 3rd of the available x-axis and might imply inaccurate slope or insufficient data. Recommended: plot `x` vs `y`.
└ @ ChaosTools ~/.julia/packages/ChaosTools/Ax8S6/src/dimensions/linear_regions.jl:131
┌ Warning: Found linear region spans less than a 3rd of the available x-axis and might imply inaccurate slope or insufficient data. Recommended: plot `x` vs `y`.
└ @ ChaosTools ~/.julia/packages/ChaosTools/Ax8S6/src/dimensions/linear_regions.jl:131
┌ Warning: Found linear region spans less than a 3rd of the available x-axis and might imply inaccurate slope or insufficient data. Recommended: plot `x` vs `y`.
└ @ ChaosTools ~/.julia/packages/ChaosTools/Ax8S6/src/dimensions/linear_regions.jl:131
┌ Warning: Found linear region spans less than a 3rd of the available x-axis and might imply inaccurate slope or insufficient data. Recommended: plot `x` vs `y`.
└ @ ChaosTools ~/.julia/packages/ChaosTools/Ax8S6/src/dimensions/linear_regions.jl:131

As you can see, using τ = 15 is not a great choice! The estimates with τ = 7 though are very good (the actual value is around λ ≈ 0.89...).

Broomhead-King Coordinates

broomhead_king(s::AbstractVector, d::Int) -> U, S, Vtr

Return the Broomhead-King coordinates of a timeseries s by performing svd on high-dimensional embedding if s with dimension d with minimum delay.


Broomhead and King coordinates is an approach proposed in [Broomhead1987] that applies the Karhunen–Loève theorem to delay coordinates embedding with smallest possible delay.

The function performs singular value decomposition on the d-dimensional matrix $X$ of $s$,

\[X = \frac{1}{\sqrt{N}}\left( \begin{array}{cccc} x_1 & x_2 & \ldots & x_d \\ x_2 & x_3 & \ldots & x_{d+1}\\ \vdots & \vdots & \vdots & \vdots \\ x_{N-d+1} & x_{N-d+2} &\ldots & x_N \end{array} \right) = U\cdot S \cdot V^{tr}.\]

where $x := s - \bar{s}$. The columns of $U$ can then be used as a new coordinate system, and by considering the values of the singular values $S$ you can decide how many columns of $U$ are "important". See the documentation page for example application.

This alternative/improvement of the traditional delay coordinates can be a very powerful tool. An example where it shines is noisy data where there is the effect of superficial dimensions due to noise.

Take the following example where we produce noisy data from a system and then use Broomhead-King coordinates as an alternative to "vanilla" delay coordinates:

using DynamicalSystems, PyPlot

ds = Systems.gissinger()
data = trajectory(ds, 1000.0, dt = 0.05)
x = data[:, 1]

L = length(x)
s = x .+ 0.5rand(L) #add noise

U, S = broomhead_king(s, 40)
"19962×40 Array{Float64,2}"

Now let's simply compare the above result with the one you get from doing a "standard" call to embed:

fig=figure(figsize= (10,6))
plot(U[:, 1], U[:, 2])
title("Broomhead-King of s")

R = embed(s, 2, 30)
plot(columns(R)...; color = "C3")
title("2D embedding of s")

we have used the same system as in the Delay Coordinates Embedding example, and picked the optimal delay time of τ = 30 (for same dt = 0.05). Regardless, the vanilla delay coordinates is much worse than the Broomhead-King coordinates.

Nearest Neighbor Prediction

Nearest neighbor timeseries prediction is a method commonly listed under nonlinear timeseries analysis. This is not part of DynamicalSystems.jl, because in JuliaDynamics we have a dedicated package for this, TimeseriesPrediction.jl.

  • Skokos2016Skokos, C. H. et al., Chaos Detection and Predictability - Chapter 1 (section 1.3.2), Lecture Notes in Physics 915, Springer (2016)
  • Kantz1994Kantz, H., Phys. Lett. A 185, pp 77–87 (1994)
  • Broomhead1987D. S. Broomhead, R. Jones and G. P. King, J. Phys. A 20, 9, pp L563 (1987)