The question underlying the expression retro-modelling
is the following: is it possible to retrieve
the original formulation of a given dynamical system
from a set of observed time series generated by
this system (and without any extra information)?
To answer this question, the Rössler system^[O. Rössler,
An Equation for Continuous Chaos, *Physics Letters*,
**57A**(5), 1976, 397-398.] was extensively used
in the previous vignettes, in particular in vignettes
`3 Modelling`

and `6 Sensitivity`

, to show
that the global modelling is very efficient to retrieve
this system under various conditions
of noise, subsampling, short length time series, etc.
In the present vignette, other three-dimensional chaotic
systems are considered in order to check the
generality of the technique.
Seventeen other systems are considered in the present
vignette which equations can be loaded as follows:

data("allMod_nVar3_dMax2") names(allMod_nVar3_dMax2)

These systems were used to produce time series. These time series can also be directly loaded as follows:

data("TSallMod_nVar3_dMax2")

In the following experiments, these time series are used
as an input of the `GPoM`

package to check if
the sets of original equations can be retrieved.

The Nosé-Hoover system was discovered by
Shūichi Nosé^[S. Nosé, A molecular dynamics method
for simulations in the canonical ensemble,
*Molecular Physics*, **52** (2), 255-268, 1984.]
and William Hoover^[W. G. Hoover, Nonlinear conductivity
and entropy in the two-body Boltzmann gas,
*Journal of Statistical Physics*, **42**(3/4), 587-600, 1986.]
in 1986 (it was rediscovered by Sprott in 1994, known as Sprott-A).
This system reads

visuEq(K = allMod_nVar3_dMax2$NH86, substit = 1)

This system is particularly interesting because, despite their deterministic couplings, the three variables are characterized by quite different signals:

TS <- TSallMod_nVar3_dMax2$NH86$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-6.2,5.2), xlab = 't', ylab = '', main = 'Nosé-Hoover', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(10,-3.5, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

And the three variables are almost completely decorrelated.

# Compute the correlation cor(TS[,2:4])

Indeed, correlation does not exceed ~0.142 (the maximum correlation is observed between $x$ and $z$), and it is almost zero between $x$ and $y$ (~0.025) and between $y$ and $z$ (~0.016). Obviously, for such a type of dynamic, the correlation may be particulary misleading information for whom would like to detect the existence of a causal coupling.

The `GPoM`

package ^[S. Mangiarotti, 2015. Low dimensional chaotic models
for the plague epidemic in Bombay, *Chaos, Solitons & Fractals*,
**81**(A), 184-196.]$^,$ ^[S. Mangiarotti, M. Peyre & M. Huc, 2016.
A chaotic model for the epidemic of Ebola virus disease
in West Africa (2013–2016). *Chaos*, **26**, 113112.]
is used to investigate the relations existing between
the three time series $x(t)$, $y(t)$ and $z(t)$ shown above.
A drastic model selection is performed by the `gPoMo`

function when trying to obtain a global model from
a set of observational time series.
The numerical integration is then used to obtain
an optimal model among the few models preselected
by the function:

# modelling from time series starting outNH86 <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 3, nPmax = 5, method = 'rk4')

With the only exception of model #6, the preselected models are either diverging (models #1, #2, #4, #5, #7 and #10) or converging to a fixed point (#3, #8 and #9). Equations for model #6 effectively correspond to the original system:

visuEq(K = outNH86$models$model6, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out4model6, substit = 1, approx = 2)

Therefore, despite the absence of correlation observed
between the three variables, the proper formulation is directly
and completely retrieved.
Note that allowing models of larger size (e.g. with `nPmax = 6`

),
three other models of nontrivial solution can be obtained
(models #12, #14 and #15).
The model #6 should be preferred from these for its more
parsimonious formulation.

The results obtained for this system are especially interesting since providing an obvious evidence that correlation is not a good tool for the detection of linkages between dynamical variables. Contrarily, the retro-modelling technique enables not only to detect the couplings. Furthermore, it goes far beyound, by making it possible to retrieve the original and complete formulation of the coupling between the three variables.

The following differential equations was introduced
by Roberto Genesio and Alberto Tesi in 1992^[R. Genesio & A. Tesi,
Harmonic balance methods for the analysis of chaotic dynamics in nonlinear
systems, *Automatica*, **28**(3), 531-548, 1992.]

$d^3x/dt^3 + a . d^2x/dt^2 + b.dx/dt + x (1 + x) = 0$,

(for a = 0.446 and b = 1.1) which can be rewritten in the canonical form:

visuEq(K = allMod_nVar3_dMax2$GT92, substit = 1)

The following set of time series $x(t)$, $y(t)$ and $z(t)$ was produced from this system.

TS <- TSallMod_nVar3_dMax2$GT92$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-1.1, 0.8), xlab = 't', ylab = '', main= 'Genesio-Tesi (1992)', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(22,-0.5, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

And the global modelling technique was then applied to this set of time series:

# modelling from time series starting outGT92 <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

Performing the same analysis as for the system Nosé-Hoover (see upper in the text), the formulation of the original system can be directly retrieved (model #1):

visuEq(K = outGT92$models$model1, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out5model1, substit = 1, approx = 2)

Note that obtained models of smaller size have
period-1 cycles (these can be obtained by setting
`nPmin = 3`

). The complexity of the observed
signal is thus not reproduced by these models
and these can be rejected.
Allowing models of larger size (nPmax = 7), interesting
solutions will be otbained (in terms of phase portrait),
but their formulation will be less concise and the previous
formulation obtained with model #1 should thus be preferred.

Numerous systems were discovered (some rediscovered)
by Julian Clinton Sprott in the 1990s^[J. C. Sprott,
Some simple chaotic flows, *Physical Review E*, **50**(2),
647-650, 1994.]. Several of these are considered
in this section, namely systems
Sprott-F, H, K, O, P, G, M, Q and S.

The Sprott-F system reads:

visuEq(K = allMod_nVar3_dMax2$SprF, substit = 1)

A specific attention was paid uuper to the Nosé-Hoover system for its very low correlation. The Sprott-F system is interesting for the opposite reason: its variables are characterized by a high level of correlation (or anticorrelation) as shown by the following plots $x(t)$, $y(t)$ and $z(t)$:

TS <- TSallMod_nVar3_dMax2$SprF$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-3.6,4.2), xlab = 't', ylab = '', main = 'Sprott-F', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(0,4, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

Correlation is high between the three variables: ~0.68 between $x(t)$ and $y(t)$, ~0.74 between $x(t)$ and $-z(t)$, and ~0.77 between $y(t)$ and $-z(t)$:

# Compute the correlation cor(TS[,2:4])

The global modelling technique was applied to this data set
using the `gPoMo`

function:

# modelling from time series starting outSprF <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

Note that allowing a larger range of model size (`nPmin = 3`

and `nPmax = 7`

),
models solutions can be rejected either because they are diverging or
oversimplified (models #3-14, #16-24, #26-27 and #30-35), or because
they are less parsimonious (models #15, #25, #28 and #29).

The formulation of the original system is thus effectively retrieved:

visuEq(K = outSprF$models$model5, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out6model5, substit = 1, approx = 2)

Similarly, results are found to be straightforward for all the Sprott systems.

The Sprott-H system reads

visuEq(K = allMod_nVar3_dMax2$SprH, substit = 1)

Although less marked, the behavior of this system is more or less similar to Sprott-F in terms of correlation (or anticorrelation).

TS <- TSallMod_nVar3_dMax2$SprH$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-5,5), xlab = 't', ylab = '', main = 'Sprott-H', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(0, -2, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

Correlation is high between $x(t)$ and $z(t)$ (~0.77), and $x(t)$ and $-y(t)$ (~0.63), and lower between $y(t)$ and $z(t)$ (~0.38):

# Compute the correlation cor(TS[,2:4])

The global modelling technique was applied to this data set:

# modelling from time series starting outSprH <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

The formulation of the original system is directly retrieved:

visuEq(K = outSprH$models$model5, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out7model5, substit = 1, approx = 2)

No simpler formulation could enable to reproduce the observed dynamic.

The Sprott-K system reads

visuEq(K = allMod_nVar3_dMax2$SprK, substit = 1)

The following data set was used:

TS <- TSallMod_nVar3_dMax2$SprK$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-3,3.6), xlab = 't', ylab= '', main = 'Sprott-K', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(0,3.5, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

and the global modelling technique was applied.

# modelling from time series starting outSprK <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1),show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

The formulation of the original system is directly retrieved:

visuEq(K = outSprK$models$model5, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out8model5, substit = 1, approx = 2)

and no simpler formulation enabling to reproduce the observed dynamic was obtained.

The Sprott-O system reads

visuEq(K = allMod_nVar3_dMax2$SprO, substit = 1)

The following data set was used:

TS <- TSallMod_nVar3_dMax2$SprO$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-1.5,1), main = 'Sprott-O', xlab = 't', ylab='', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(0,-1, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied.

# modelling from time series starting outSprO <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

The formulation of the original was directly retrieved:

visuEq(K = outSprO$models$model2, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out9model2, substit = 1, approx = 2)

and no simpler formulation enabling to reproduce the observed dynamic was obtained.

The Sprott-P system reads

visuEq(K = allMod_nVar3_dMax2$SprP, substit = 1)

The following data set was used:

TS <- TSallMod_nVar3_dMax2$SprP$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-1.2,2.2), main='Sprott-P', xlab = 't', ylab = '', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(8,2, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied.

# modelling from time series starting outSprP<- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

The formulation of the original system is directly retrieved:

visuEq(K = outSprP$models$model5, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out10model5, substit = 1, approx = 2)

No simpler formulation enabling to reproduce the observed dynamic was obtained.

The Sprott-G system reads

visuEq(K = allMod_nVar3_dMax2$SprG, substit = 1)

The following data set was used:

TS <- TSallMod_nVar3_dMax2$SprG$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-2.5,2), xlab = 't', ylab = '', main='Sprott-G', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(0, 2 , c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was then applied.

# modelling from time series starting outSprG <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 3, nPmax = 6, method = 'rk4')

The formulation of the original system was directly retrieved:

visuEq(K = outSprG$models$model5, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out11model5, substit = 1, approx = 2)

No simpler formulation enabling to reproduce the observed dynamic was obtained.

The Sprott-M system reads

visuEq(K = allMod_nVar3_dMax2$SprM, substit = 1)

The following data set was used:

TS <- TSallMod_nVar3_dMax2$SprM$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-5.1,2.8), main='Sprott-M', xlab = 't', ylab = '', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(2, -3, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was then applied.

# modelling from time series starting outSprM <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

The formulation of the original system was directly retrieved:

visuEq(K = outSprM$models$model2, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out12model2, substit = 1, approx = 2)

No simpler formulation enabling to reproduce the observed dynamic was obtained.

The Sprott-Q system reads

visuEq(K = allMod_nVar3_dMax2$SprQ, substit = 1)

The following data set was used:

TS <- TSallMod_nVar3_dMax2$SprQ$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-7,8.8), main = 'Sprott-Q', xlab = 't', ylab = '', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(4, 7, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied.

# modelling from time series starting outSprQ <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

The formulation of the original is directly retrieved:

visuEq(K = outSprQ$models$model2, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out13model2, substit = 1, approx = 2)

No simpler formulation enabling to reproduce the observed dynamic was obtained.

The Sprott-S system reads

visuEq(K = allMod_nVar3_dMax2$SprS, substit = 1)

The following data set was used:

TS <- TSallMod_nVar3_dMax2$SprS$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-4,2), main = 'Sprott-S', xlab = 't', ylab = '', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(0,-2.5, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied.

# modelling from time series starting outSprS<- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 6, nPmax = 6, method = 'rk4')

The formulation of the original system is directly retrieved:

visuEq(K = outSprS$models$model5, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out14model5, substit = 1, approx = 2)

No simpler formulation enabling to reproduce the observed dynamic was obtained.

The Lorenz system^[E. N. Lorenz, Deterministic non-periodic flow,
*Journal of the Atmospheric Sciences*, **20**(2), 130–141 (1963).]
discovered by Edouard N. Lorenz in 1963 is now considered.
This system reads

visuEq(K = allMod_nVar3_dMax2$L63, substit = 1, approx = 4)

The following data set was considered for the modelling:

TS <- TSallMod_nVar3_dMax2$L63$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-25,45), xlab = 't', ylab = '', main = 'Lorenz 1963', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(8, -5, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied to this set of series with the hope to retrieve the original equations.

# modelling from time series starting outL63 <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 1500, nPmin = 6, nPmax = 7, method = 'rk4')

Surprisingly, the simplest model able to reproduce the original phase portrait has six terms, that is, one less term than the original model: the term $-y$ is missing in the second equation. This obtained morel reads

visuEq(K = outL63$models$model5, substit = 1, approx = 2)

#data("data_vignetteVII$out1model5") visuEq(K = data_vignetteVII$out1model5, substit = 1, approx = 2)

This term is not detected because it is of small influence on the dynamics (actually the dynamic of variable $y$ is more or less similar to the dynamic of $x$ -- See upper the observed time series -- and its coefficient is 28 times smaller). Three 7-parameter models were also obtained, all having a phase space very similar to the original one. The original model (#18) is included among them :

visuEq(K = outL63$models$model18, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out1model18, substit = 1, approx = 2)

Since model #5 is more parsomonious, it can be considered that six of the seven terms of the Lorenz system are thus clearly detected. A more in-deep investigation might enable the detection of the term of weaker influence, but this remains unsure. The phase portraits of four potentially valid models are very close to the original phase portraits, it may be quite quite challenging to retrieve the full formulation of the original equations.

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #5') lines(outL63$stockoutreg$model5[,1], outL63$stockoutreg$model5[,2], type = 'l', col = 'red')

library(float) plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #5') lines(dbl(data_vignetteVII$outL63_stockoutreg_model5[,1]), dbl(data_vignetteVII$outL63_stockoutreg_model5[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #15') lines(outL63$stockoutreg$model15[,1], outL63$stockoutreg$model15[,2], type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #15') lines(dbl(data_vignetteVII$outL63_stockoutreg_model15[,1]), dbl(data_vignetteVII$outL63_stockoutreg_model15[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #18') lines(dbl(outL63$stockoutreg$model18[,1]), dbl(outL63$stockoutreg$model18[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #18') lines(dbl(data_vignetteVII$outL63_stockoutreg_model18[,1]), dbl(data_vignetteVII$outL63_stockoutreg_model18[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #19') lines(dbl(outL63$stockoutreg$model19[,1]), dbl(outL63$stockoutreg$model19[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #19') lines(dbl(data_vignetteVII$outL63_stockoutreg_model19[,1]), dbl(data_vignetteVII$outL63_stockoutreg_model19[,2]), type = 'l', col = 'red')

A very refined validation technique based on
the transition matrices^[Mangiarotti S., Coudret R.,
Drapeau L. & Jarlan L., Global models’ formulation and
associated transition matrices for the Rössler
system, the electrodissolution of copper in phosphoric
acid and the cycles of rainfed wheat observed from satellite
in North Morocco. *Supplemental Material* to article
“Polynomial Model Search and Global Modelling – two new
algorithms for global modelling of chaos” by Mangiarotti S.,
Coudret R., Drapeau L. & Jarlan L., 2012, *Physical Review E*,
**86**(4), 046205.]
may be a powerful tool to reach such a level of precision
for the distinction.

The system discovered by Bill Burke and Robert
Show^[R. Shaw, Strange attractor, chaotic behavior
and information flow, *Zeitschrift für Naturforsch A*,
**36**, 80-112, 1981]
in 1981 is a Lorenz-like system. This system reads

visuEq(K = allMod_nVar3_dMax2$BS81, substit = 1, approx = 4)

The following set of variable was used to test the possibility to retrieve the original set of equations:

TS <- TSallMod_nVar3_dMax2$BS81$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-2,2.2), main = 'Burke-Shaw 1981', xlab = 't', ylab = '', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(9, 2.2, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied:

# modelling from time series starting outBS81 <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 1500, nPmin = 5, nPmax = 6, method = 'rk4')

Here also, the simplest model has one less term than the original model:

visuEq(K = outBS81$models$model3, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out2model3, substit = 1, approx = 2)

Three 6-parameter models were also obtained which all have a phase space very similar to the original one. The original model

visuEq(K = outBS81$models$model11, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out2model11, substit = 1, approx = 2)

is included among them.

A more in-deep investigation will be required to check if the term of weaker influence can be detected. At this stage, it can already be observed that the phase portrait of the model with the proper formulation (model #11 in red) seems closer to the original phase space (in black).

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #3') lines(outBS81$stockoutreg$model3[,1], outBS81$stockoutreg$model3[,2], type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #3') lines(dbl(data_vignetteVII$outBS81_stockoutreg_model3[,1]), dbl(data_vignetteVII$outBS81_stockoutreg_model3[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #9') lines(outBS81$stockoutreg$model9[,1], outBS81$stockoutreg$model9[,2], type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #9') lines(dbl(data_vignetteVII$outBS81_stockoutreg_model9[,1]), dbl(data_vignetteVII$outBS81_stockoutreg_model9[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #11') lines(outBS81$stockoutreg$model11[,1], outBS81$stockoutreg$model11[,2], type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #11') lines(dbl(data_vignetteVII$outBS81_stockoutreg_model11[,1]), dbl(data_vignetteVII$outBS81_stockoutreg_model11[,2]), type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #12') lines(outBS81$stockoutreg$model12[,1], outBS81$stockoutreg$model12[,2], type = 'l', col = 'red')

plot(TS[,2], TS[,3], type = 'l', xlab = 'x(t)', ylab = 'y(t)', main = 'model #12') lines(dbl(data_vignetteVII$outBS81_stockoutreg_model12[,1]), dbl(data_vignetteVII$outBS81_stockoutreg_model12[,2]), type = 'l', col = 'red')

The Lorenz-1984 system^[E. N. Lorenz, Irregularity: a fundamental property of the atmosphere, *Tellus A*, **36**, 98-110, 1984.] was discovered by E.N. Lorenz in 1984.
This system is much more complex than the systems previously
considered.

visuEq(K = allMod_nVar3_dMax2$L84, substit = 1)

It is also a quadratic system, but it includes eleven terms. The following data set was used for the analysis:

TS <- TSallMod_nVar3_dMax2$L84$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-2,2.2), main = 'Lorenz 1984', xlab = 't', ylab = '', col = 'blue') lines(TS[,1], TS[,3], type = 'l', col = 'red') lines(TS[,1], TS[,4], type = 'l', col = 'green') legend(16, -1, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied to this set of series.

# modelling from time series starting outL84 <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 5, nPmax = 11, method = 'rk4')

The simplest obtained model has only five terms

visuEq(K = outL84$models$model3, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out3model3, substit = 1, approx = 2)

However, its attractor is periodic and symmetric. Therefore, it is obviously not a precise approximation of the observed time series. The first model showing a high level of complexity is the model #9. This model has 6 terms

visuEq(K = outL84$models$model9, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out3model9, substit = 1, approx = 2)

Although several terms are missing (two in the first equation, three in the second one), it can be noticed that the detected terms are effectively present in the original system.

Many other models were obtained (models with 7 terms: #19 and #22; with 8 terms: #34, #40 and #49; with 9 terms: #55, #64, #76 and #77; with 10 terms: #83, #92-94, #97-99, #111 and #113; or with 11 terms: #119, #128, #129, #133, #134, #136, #139, #140-142, #147-149 and #155-158). Most of these models have strongly negative values along the $x$ abscissa and can thus be rejected. Only a few of them may be considered as potentially acceptable (models #9, #55, #77, #113, #119, #141 and #158). A visual inspection of the phase portraits obviously shows that the best agreement is obtained with the model #141.

visuEq(K = outL84$models$model141, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out3model141, substit = 1, approx = 2)

which formulation effectively corresponds to the original Lorenz-1984 system.

The Chlouverakis-Sprott system was adapted from the Lorenz-1984 system.
It was introduced by Konstantinos Chlouverakis and Julian Sprott
in 2004.^[K. E. Chlouverakis & J. C. Sprott, A comparison of correlation
and Lyapunov dimensions, *Physica D*, **200**, 156-164, 2004.]

visuEq(K = allMod_nVar3_dMax2$CS2004, substit = 1)

The following data set was used for the analysis:

TS <- TSallMod_nVar3_dMax2$CS2004$reconstr plot(TS[,1], TS[,2], type = 'l', xlab = 't', ylab ='', main='Chlouverakis-Sprott', col = 'blue') lines(TS[,1], TS[,3], type = 'l', xlab = 't', ylab = 'y(t)', col = 'red') lines(TS[,1], TS[,4], type = 'l', xlab = 't', ylab = 'z(t)', col = 'green') legend(40, -1.5, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied to this set of time series.

# modelling from time series starting outCS2004 <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 3, nPmax = 8, method = 'rk4')

Three types of models were obtained (type-1: #7, #13, #15, #23, #25, #38, #40, #43; type-2: #12, #22, #27, #37; and type-3: #26, #30, #41, #44-46, #49) among which only type-3 models seem compatible with the original phase portrait (the other types can only reproduce a small part of the original phase space). Some spurious effects being observed in models #30 and #49, only five models may be potential candidates: models #26, #41 and #44-46. The smallest model (#26) obtained with the third type has the following formulation

visuEq(K = outCS2004$models$model26, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out15model26, substit = 1, approx = 2)

It is almost identical to the original system since only one term is missing in the second equation.

The complete formulation is obtained in model #41. However, it is difficult to tell which model is the best when comparing models #26, #41, #44 and #45. These results suggest that the original system may be reductible to a 7-term sytem.

A toroidal and symmetrical system was introduced
by Dequan Li in 2007^[Dequan Li, A three-scroll chaotic
attractor, *Physics Letters A*, **372**, 387-393, 2007.]:

visuEq(K = allMod_nVar3_dMax2$Li2007, substit = 1, approx = 2)

The following data set was obtained from this system:

TS <- TSallMod_nVar3_dMax2$Li2007$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-160, 230), xlab = 't', ylab = '', main= ' Li 2007', col = 'blue') lines(TS[,1], TS[,3], type = 'l', xlab = 't', ylab = 'y(t)', col = 'red') lines(TS[,1], TS[,4], type = 'l', xlab = 't', ylab = 'z(t)', col = 'green') legend(0.3, -50, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was then applied to this set of time series.

# modelling from time series starting outLi2007<- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 9, nPmax = 9, method = 'rk4')

Only the model #13 is able to reproduce the original phase portrait.

visuEq(K = outLi2007$models$model13, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out16model13, substit = 1, approx = 2)

The same result is obtained when allowing models
of smaller sizes (e.g. `nPmin = 3`

).
The original formulation is thus appropriately retrieved.

The cord attractor^[C. Letellier & L. A. Aguirre, Required criteria
for recognizing new types of chaos : Application to the “cord” attractor,
*Physical Review E*, **85**, 036204, 2012.] was introduced by Luis Aguirre and
Christophe Letellier in 2012 by a modification of the Lorenz-1984 system.
The system reads:

visuEq(K = allMod_nVar3_dMax2$Cord2012, substit = 1)

The following data set was used for the analysis:

TS <- TSallMod_nVar3_dMax2$Cord2012$reconstr plot(TS[,1], TS[,2], type = 'l', ylim = c(-20,30), xlab = 't', ylab = '', main = 'Cord 2012', col = 'blue') lines(TS[,1], TS[,3], type = 'l', xlab = 't', ylab = 'y(t)', col = 'red') lines(TS[,1], TS[,4], type = 'l', xlab = 't', ylab = 'z(t)', col = 'green') legend(17, -8, c("x(t)", "y(t)", "z(t)"), col=c('blue', 'red', 'green'), lty=1, cex = 0.8)

The global modelling technique was applied to this set of time series.

# modelling from time series starting outCord2012 <- gPoMo(data = TS[,2:4], tin = TS[,1], dMax = 2, nS = c(1,1,1), show = 1, IstepMin = 10, IstepMax = 15000, nPmin = 7, nPmax = 11, method = 'rk4')

Several obtained models were able to reproduce, more or less, the original phase space are obtained. The smallest models (#9 and #10) have seven terms and are not chaotic. Their formulation is similar to the original system. All the terms detected in model #9 are effectively included in the original system (four terms are not detected).

visuEq(K = outCord2012$models$model9, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out17model9, substit = 1, approx = 2)

Only one spurious term (lower in magnitude compared to the others) is detected in model #10 (term $z^2$ in the first equations) and five terms are not detected:

visuEq(K = outCord2012$models$model10, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out17model10, substit = 1, approx = 2)

Among all the obtained models, the phase portrait of model #62 is the more similar to the original one.

visuOutGP(out17,selecmod = 62, prioMinMax = 'model')

In this formulation, eight of the eleven terms are correctly detected.
The same spurious term is wrongly detected (term `z^2`

in the first equations).

visuEq(K = outCord2012$models$model62, substit = 1, approx = 2)

visuEq(K = data_vignetteVII$out17model62, substit = 1, approx = 2)

The cord system cannot be retrieved completely in the present case, possibly because its formulation is too complex (too many terms) and not minimal. The most important terms can be detected.

In this vignette, it is shown that, in many cases, the original systems can be fully retrieved. It is in particular the case when the formulation of the original system is concise enough. But it is not always the case. When the systems formulation is more complex (more terms), an exhaustive detection becomes more difficult and only one part of the algebraic terms can be detected. In such conditions the global modelling technique can be used to obtain a reduced formulation of the dynamical systems.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.