Delay Coordinates
A timeseries recorded in some manner from a dynamical system can be used to gain information about the dynamics of the entire phase-space of the system. This can be done by reconstructing a new phase-space from the timeseries. One method that can do this is what is known as delay coordinates embedding or delay coordinates reconstruction.
This is done through the reconstruct interface:
#DynamicalSystemsBase.reconstruct — Function.
reconstruct(s, D, τ)
Reconstruct s using the delay coordinates embedding with D temporal neighbors and delay τ and return the result as a Dataset.
Description
Single Timeseries
If τ is an integer, then the n-th entry of the embedded space is
If instead τ is a vector of integers, so that length(τ) == D, then the n-th entry is
The reconstructed dataset can have same invariant quantities (like e.g. lyapunov exponents) with the original system that the timeseries were recorded from, for proper D and τ [1, 2]. The case of different delay times allows reconstructing systems with many time scales, see [3].
Notice - The dimension of the returned dataset is D+1!
Multiple Timeseries
To make a reconstruction out of a multiple timeseries (i.e. trajectory) the number of timeseries must be known by type, so s can be either:
s::AbstractDataset{B}s::SizedAray{A, B}
If the trajectory is for example (x, y) and τ is integer, then the n-th entry of the embedded space is
If τ is an AbstractMatrix{Int}, so that size(τ) == (D, B), then we have
Notice - The dimension of the returned dataset is (D+1)*B!
References
[1] : F. Takens, Detecting Strange Attractors in Turbulence — Dynamical Systems and Turbulence, Lecture Notes in Mathematics 366, Springer (1981)
[2] : T. Sauer et al., J. Stat. Phys. 65, pp 579 (1991)
[3] : K. Judd & A. Mees, Physica D 120, pp 273 (1998)
Here are some examples of reconstructing a 3D continuous chaotic system:
using DynamicalSystems, PyPlot ds = Systems.gissinger(ones(3)) data = trajectory(ds, 1000.0, dt = 0.05) xyz = columns(data) figure(figsize = (12,10)) k = 1 for i in 1:3 for τ in [5, 30, 100] R = reconstruct(xyz[i], 1, τ) ax = subplot(3,3,k) plot(R[:, 1], R[:, 2], color = "C$(k-1)", lw = 0.8) title("var = $i, τ = $τ") global k+=1 end end tight_layout() suptitle("2D reconstructed space") subplots_adjust(top=0.9)

τ and dt
Keep in mind that whether a value of τ is "reasonable" for continuous systems depends on dt. In the above example the value τ=30 is good, only for the case of using dt = 0.05. For shorter/longer dt one has to adjust properly τ so that their product τ*dt is the same.
You can also reconstruct multidimensional timeseries. For this to be possible, the number of timeseries must be known by Type:
using StaticArrays: Size a = rand(1000, 3) # my trajectory A = Size(1000, 3)(a) # create array with the size as Type information R = reconstruct(A, 2, 2) #aaaall good
9-dimensional Dataset{Float64} with 996 points
 0.0333477  0.967267   0.162712     …  0.194024   0.169191   0.252666
 0.630255   0.222639   0.755131        0.94983    0.455341   0.346721
 0.521731   0.709986   0.639572        0.569019   0.154882   0.747215
 0.332819   0.842421   0.717376        0.409041   0.749972   0.311119
 0.194024   0.169191   0.252666        0.234317   0.24932    0.938725
 0.94983    0.455341   0.346721     …  0.705545   0.595288   0.440244
 0.569019   0.154882   0.747215        0.874132   0.5615     0.165738
 0.409041   0.749972   0.311119        0.215505   0.942835   0.905429
 0.234317   0.24932    0.938725        0.465969   0.0103991  0.422649
 0.705545   0.595288   0.440244        0.452805   0.725377   0.0549314
 ⋮                                  ⋱
 0.352406   0.302666   0.000731456     0.379557   0.766445   0.6546
 0.370318   0.445704   0.176596        0.76893    0.477902   0.903401
 0.194646   0.639107   0.192121        0.852681   0.265334   0.349134
 0.328198   0.330089   0.605354        0.762692   0.0287925  0.699436
 0.379557   0.766445   0.6546       …  0.931412   0.702562   0.146645
 0.76893    0.477902   0.903401        0.852133   0.324826   0.0784649
 0.852681   0.265334   0.349134        0.641115   0.688045   0.331545
 0.762692   0.0287925  0.699436        0.855515   0.271941   0.0331421
 0.931412   0.702562   0.146645        0.0861718  0.380648   0.398455
ds = Systems.towel(); tr = trajectory(ds, 10000) R = reconstruct(tr, 2, 2) # Dataset size is also known by Type!
9-dimensional Dataset{Float64} with 9997 points
 0.085     -0.121       0.075     …  0.837347   0.0372633   0.555269
 0.285813  -0.0675286   0.238038     0.51969    0.0616256   0.940906
 0.76827   -0.038933    0.672094     0.966676  -0.00171595  0.2225
 0.681871   0.0508933   0.825263     0.112748   0.0674955   0.653573
 0.837347   0.0372633   0.555269     0.386547  -0.0886542   0.869349
 0.51969    0.0616256   0.940906  …  0.910741  -0.0316828   0.411607
 0.966676  -0.00171595  0.2225       0.306095   0.0689305   0.909129
 0.112748   0.0674955   0.653573     0.824263  -0.056185    0.326064
 0.386547  -0.0886542   0.869349     0.545332   0.0508239   0.819404
 0.910741  -0.0316828   0.411607     0.954994   0.00453815  0.569534
 ⋮                                ⋱
 0.914702  -0.0315439   0.294266     0.90246    0.0242141   0.539502
 0.289932   0.0641239   0.778698     0.335976   0.0735803   0.943945
 0.793854  -0.0552801   0.664223     0.86657   -0.0497658   0.214728
 0.62671    0.0557527   0.832001     0.430816   0.0535742   0.62743
 0.90246    0.0242141   0.539502  …  0.936955  -0.0200112   0.894333
 0.335976   0.0735803   0.943945     0.237481   0.0983265   0.353212
 0.86657   -0.0497658   0.214728     0.681538  -0.0476555   0.883219
 0.430816   0.0535742   0.62743      0.836353   0.0363264   0.380351
 0.936955  -0.0200112   0.894333     0.515471   0.0534613   0.898152