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

ChaosTools.numericallyapunov — Function

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.

Description

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))")
        legend()
        tight_layout()
    end
end
tight_layout(pad=0.3); fig

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]
length(x)

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])")
fig.tight_layout(pad=0.3); fig

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 Vector{Float64}:
   0.0
   4.285178117517708
   8.924780522479637
  15.012203311102235
  20.05533894475613
  18.062350952804728
   9.898343637398332
   2.199113375749754
  -2.6729722259323863
  -5.33812377718313
   ⋮
  -3.810825689423999
  -4.900396720248087
  -6.241066898514638
  -7.940150591858717
  -9.900960829002475
 -11.678879507784437
 -12.495362895268917
 -11.76970885072464
  -9.786325967697445

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

0:200

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

ks2 = 0:4:200

0:4:200

Now we plot some example computations

fig = figure()
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))")
end

legend()
xlabel("k (0.05×t)")
ylabel("E - E(0)")
title("Continuous Reconstruction Lyapunov")
fig.tight_layout(pad=0.3); fig

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

ChaosTools.broomhead_king — Function

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.

Description

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)
summary(U)

"19962×40 Matrix{Float64}"

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))
subplot(1,2,1)
plot(U[:, 1], U[:, 2])
title("Broomhead-King of s")

subplot(1,2,2)
R = embed(s, 2, 30)
plot(columns(R)...; color = "C3")
title("2D embedding of s")
fig.tight_layout(pad=0.3)

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.

DyCA - Dynamical Component Analysis

ChaosTools.dyca — Function

dyca(data, eig_thresold) -> eigenvalues, proj_mat, projected_data

Compute the Dynamical Component analysis (DyCA) of the given data ^[Uhl2018] used for dimensionality reduction.

Return the eigenvalues, projection matrix, and reduced-dimension data (which are just data*proj_mat).

Description

Dynamical Component Analysis (DyCA) is a method to detect projection vectors to reduce the dimensionality of multi-variate, high-dimensional deterministic datasets. Unlike methods like PCA or ICA that make a stochasticity assumption, DyCA relies on a determinacy assumption on the time-series and is based on the solution of a generalized eigenvalue problem. After choosing an appropriate eigenvalue threshold and solving the eigenvalue problem, the obtained eigenvectors are used to project the high-dimensional dataset onto a lower dimension. The obtained eigenvalues measure the quality of the assumption of linear determinism for the investigated data. Furthermore, the number of the generalized eigenvalues with a value of approximately 1.0 are a measure of the number of linear equations contained in the dataset. This property is useful in detecting regions with highly deterministic parts in the time-series and also as a preprocessing step for reservoir computing of high-dimensional spatio-temporal data.

The generalised eigenvalue problem we solve is:

\[C_1 C_0^{-1} C_1^{\top} \bar{u} = \lambda C_2 \bar{u} \]

where $C_0$ is the correlation matrix of the data with itself, $C_1$ the correlation matrix of the data with its derivative, and $C_2$ the correlation matrix of the derivative of the data with itself. The eigenvectors $\bar{u}$ with eigenvalues approximately 1 and their $C_1^{-1} C_2 u$ counterpart, form the space where the data is projected onto.

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)
Uhl2018B Seifert, K Korn, S Hartmann, C Uhl, Dynamical Component Analysis (DYCA): Dimensionality Reduction for High-Dimensional Deterministic Time-Series, 10.1109/mlsp.2018.8517024, 2018 IEEE 28th International Workshop on Machine Learning for Signal Processing (MLSP)