Delay coordinates embedding
A timeseries recorded in some manner from a dynamical system can be used to gain information about the dynamics of the entire state space of the system. This can be done by constructing a new state space from the timeseries. One method that can do this is what is known as delay coordinates embedding or delay coordinates reconstruction.
The main functions to use for embedding some input data are embed
or genembed
. Both functions return a StateSpaceSet
.
Timeseries embedding
DelayEmbeddings.embed
— Functionembed(s, d, τ [, h])
Embed s
using delay coordinates with embedding dimension d
and delay time τ
and return the result as a StateSpaceSet
. Optionally use weight h
, see below.
Here τ > 0
, use genembed
for a generalized version.
Description
If τ
is an integer, then the $n$-th entry of the embedded space is
\[(s(n), s(n+\tau), s(n+2\tau), \dots, s(n+(d-1)\tau))\]
If instead τ
is a vector of integers, so that length(τ) == d-1
, then the $n$-th entry is
\[(s(n), s(n+\tau[1]), s(n+\tau[2]), \dots, s(n+\tau[d-1]))\]
The resulting set can have same invariant quantities (like e.g. Lyapunov exponents) with the original system that the timeseries were recorded from, for proper d
and τ
. This is known as the Takens embedding theorem [Takens1981] [Sauer1991]. The case of different delay times allows embedding systems with many time scales, see[Judd1998].
If provided, h
can be weights to multiply the entries of the embedded space. If h isa Real
then the embedding is
\[(s(n), h \cdot s(n+\tau), h^2 \cdot s(n+2\tau), \dots,h^{d-1} \cdot s(n+γ\tau))\]
Otherwise h
can be a vector of length d-1
, which the decides the weights of each entry directly.
References
[Takens1981] : F. Takens, Detecting Strange Attractors in Turbulence — Dynamical Systems and Turbulence, Lecture Notes in Mathematics 366, Springer (1981)
[Sauer1991] : T. Sauer et al., J. Stat. Phys. 65, pp 579 (1991)
If the data values are very strongly discretized (e.g., integers or floating-point numbers with very small bits), this can result to distances between points in the embedded space being 0. This is problematic for several library functions. Best practice here is to add noise to your original timeseries before embedding, e.g., s = s .+ 1e-15randn(length(s))
.
Here are some examples of embedding a 3D continuous chaotic system:
using DelayEmbeddings
x = cos.(0:0.1:1)
11-element Vector{Float64}:
1.0
0.9950041652780258
0.9800665778412416
0.955336489125606
0.9210609940028851
0.8775825618903728
0.8253356149096783
0.7648421872844885
0.6967067093471654
0.6216099682706644
0.5403023058681398
embed(x, 3, 1)
3-dimensional StateSpaceSet{Float64} with 9 points
1.0 0.995004 0.980067
0.995004 0.980067 0.955336
0.980067 0.955336 0.921061
0.955336 0.921061 0.877583
0.921061 0.877583 0.825336
0.877583 0.825336 0.764842
0.825336 0.764842 0.696707
0.764842 0.696707 0.62161
0.696707 0.62161 0.540302
Keep in mind that whether a value of τ
is "reasonable" for continuous time systems depends on the sampling time Δt
.
Embedding Structs
The high level function embed
utilizes a low-level interface for creating embedded vectors on-the-fly. The high level interface simply loops over the low level interface.
DelayEmbeddings.DelayEmbedding
— TypeDelayEmbedding(γ, τ, h = nothing) → `embedding`
Return a delay coordinates embedding structure to be used as a function-like-object, given a timeseries and some index. Calling
embedding(s, n)
will create the n
-th delay vector of the embedded space, which has γ
temporal neighbors with delay(s) τ
. γ
is the embedding dimension minus 1, τ
is the delay time(s) while h
are extra weights, as in embed
for more.
Be very careful when choosing n
, because @inbounds
is used internally. Use τrange
!
DelayEmbeddings.τrange
— Functionτrange(s, de::AbstractEmbedding)
Return the range r
of valid indices n
to create delay vectors out of s
using de
.
Generalized embeddings
DelayEmbeddings.genembed
— Functiongenembed(s, τs, js = ones(...); ws = nothing) → ssset
Create a generalized embedding of s
which can be a timeseries or arbitrary StateSpaceSet
, and return the result as a new StateSpaceSet
.
The generalized embedding works as follows:
τs
denotes what delay times will be used for each of the entries of the delay vector. It is recommended thatτs[1] = 0
.τs
is allowed to have negative entries as well.js
denotes which of the timeseries contained ins
will be used for the entries of the delay vector.js
can contain duplicate indices.ws
are optional weights that weight each embedded entry (the i-th entry of the delay vector is weighted byws[i]
). If provided, it is recommended thatws[1] == 1
.
τs, js, ws
are tuples (or vectors) of length D
, which also coincides with the embedding dimension. For example, imagine input trajectory $s = [x, y, z]$ where $x, y, z$ are timeseries (the columns of the StateSpaceSet
). If js = (1, 3, 2)
and τs = (0, 2, -7)
the created delay vector at each step $n$ will be
\[(x(n), z(n+2), y(n-7))\]
Using ws = (1, 0.5, 0.25)
as well would create
\[(x(n), \frac{1}{2} z(n+2), \frac{1}{4} y(n-7))\]
js
can be skipped, defaulting to index 1 (first timeseries) for all delay entries, while it has no effect if s
is a timeseries instead of a StateSpaceSet
.
See also embed
. Internally uses GeneralizedEmbedding
.
DelayEmbeddings.GeneralizedEmbedding
— TypeGeneralizedEmbedding(τs, js = ones(length(τs)), ws = nothing) -> `embedding`
Return a delay coordinates embedding structure to be used as a function. Given a timeseries or trajectory (i.e. StateSpaceSet
) s
and calling
embedding(s, n)
will create the delay vector of the n
-th point of s
in the embedded space using generalized embedding (see genembed
).
js
is ignored for timeseries input s
(since all entries of js
must be 1
in this case) and in addition js
defaults to (1, ..., 1)
for all τ
.
Be very careful when choosing n
, because @inbounds
is used internally. Use τrange
!
StateSpaceSet reference
StateSpaceSets.StateSpaceSet
— TypeStateSpaceSet{D, T} <: AbstractStateSpaceSet{D,T}
A dedicated interface for sets in a state space. It is an ordered container of equally-sized points of length D
. Each point is represented by SVector{D, T}
. The data are a standard Julia Vector{SVector}
, and can be obtained with vec(ssset::StateSpaceSet)
. Typically the order of points in the set is the time direction, but it doesn't have to be.
When indexed with 1 index, StateSpaceSet
is like a vector of points. When indexed with 2 indices it behaves like a matrix that has each of the columns be the timeseries of each of the variables. When iterated over, it iterates over its contained points. See description of indexing below for more.
StateSpaceSet
also supports almost all sensible vector operations like append!, push!, hcat, eachrow
, among others.
Description of indexing
In the following let i, j
be integers, typeof(X) <: AbstractStateSpaceSet
and v1, v2
be <: AbstractVector{Int}
(v1, v2
could also be ranges, and for performance benefits make v2
an SVector{Int}
).
X[i] == X[i, :]
gives thei
th point (returns anSVector
)X[v1] == X[v1, :]
, returns aStateSpaceSet
with the points in those indices.X[:, j]
gives thej
th variable timeseries (or collection), asVector
X[v1, v2], X[:, v2]
returns aStateSpaceSet
with the appropriate entries (first indices being "time"/point index, while second being variables)X[i, j]
value of thej
th variable, at thei
th timepoint
Use Matrix(ssset)
or StateSpaceSet(matrix)
to convert. It is assumed that each column of the matrix
is one variable. If you have various timeseries vectors x, y, z, ...
pass them like StateSpaceSet(x, y, z, ...)
. You can use columns(dataset)
to obtain the reverse, i.e. all columns of the dataset in a tuple.
- Judd1998K. Judd & A. Mees, Physica D 120, pp 273 (1998)
- Farmer1988Farmer & Sidorowich, Exploiting Chaos to Predict the Future and Reduce Noise"