Simplex projection and delay embedding

EDM

time series

nonlinear dynamics

ecology tutorial

Reconstruct a population’s attractor from a single time series, forecast it with simplex projection in base R, and separate deterministic dynamics from red noise.

Author

Tidy Ecology

Published

2026-08-22

An annual count series wobbles up and down for thirty years. The autocorrelation function says the wobble is not white noise, and an ARIMA model fits it well enough. Nothing in that workflow can tell you whether the wobble comes from a low-dimensional deterministic system that you could, in principle, forecast, or from a linear process pushed around by weather.

Empirical dynamic modelling asks the question the other way round. Instead of fitting an equation, it treats the series as a shadow of a trajectory moving on an attractor, and reconstructs that attractor from the series itself. This post builds the two pieces you need: the delay embedding, and simplex projection, the nearest-neighbour forecast that runs on it. Both are twenty lines of base R.

The payoff is a test that no linear model can perform: does knowing where the system is now, in reconstructed state space, tell you more about the next step than the best linear predictor can? For the series below the answer is yes, and the same machinery says no for a linear series with an identical autocorrelation function.

The idea: one variable carries the whole attractor

A population does not live in one dimension. Abundance next year depends on abundance now, on the predator, on the nutrient pool, on whatever else is coupled to it. That full state is a point moving on an attractor in some multi-dimensional space, and you measured one coordinate of it.

Takens showed that this is enough. Under mild conditions, the vector of lagged values

\[\mathbf{v}_t = \big(x_t,\; x_{t-\tau},\; x_{t-2\tau},\; \ldots,\; x_{t-(E-1)\tau}\big)\]

traces out a shape that is a faithful copy of the true attractor: distinct states stay distinct, and nearby points stay nearby. The unmeasured variables leave their fingerprint in the history of the one you did measure. E is the embedding dimension and tau the lag; for annual ecological series tau = 1 is almost always the right choice, so the only real knob is E.

The practical consequence is the one that matters here. If the reconstruction is faithful, then points that are close in delay coordinates are close in the true state, so they should do the same thing next. That turns forecasting into a lookup: find the neighbours, see where they went, average.

The engine

Two functions. make_block stacks the delay vectors; simplex finds the E+1 nearest neighbours of each point, weights them by distance, and predicts tp steps ahead. The prediction for each point leaves that point out, so nothing is fitted to itself.

make_block <- function(x, E, tau = 1) {
  n <- length(x)
  tt <- seq((E - 1) * tau + 1, n)
  M <- matrix(NA_real_, nrow = length(tt), ncol = E)
  for (j in seq_len(E)) M[, j] <- x[tt - (j - 1) * tau]
  list(M = M, t = tt)
}

simplex <- function(x, E, tau = 1, tp = 1, theiler = 0, target = NULL) {
  if (is.null(target)) target <- x
  bl <- make_block(x, E, tau)
  M <- bl$M; tt <- bl$t
  ok <- (tt + tp) >= 1 & (tt + tp) <= length(target)
  M <- M[ok, , drop = FALSE]; tt <- tt[ok]
  y <- target[tt + tp]
  D <- as.matrix(dist(M))
  pred <- rep(NA_real_, length(tt))
  for (i in seq_along(tt)) {
    d <- D[i, ]
    d[abs(tt - tt[i]) <= theiler] <- Inf   # drop self and temporal neighbours
    o <- order(d)[seq_len(E + 1)]
    dn <- d[o]
    if (!all(is.finite(dn))) next
    w <- if (dn[1] == 0) as.numeric(dn == 0) else exp(-dn / dn[1])
    w <- pmax(w, 1e-8)
    pred[i] <- sum(w * y[o]) / sum(w)
  }
  keep <- !is.na(pred)
  list(rho = suppressWarnings(cor(pred[keep], y[keep])),
       pred = pred, obs = y, t = tt)
}

The weights exp(-d / d1) scale each neighbour by its distance relative to the closest one, so an identical past state dominates and a far-off state barely contributes. The simplex of E+1 points is the smallest set that can enclose a point in E dimensions, which is where the method gets its name.

theiler deserves a note now and a warning later. Setting it to zero drops only the point being predicted. Setting it higher drops every library point within that many time steps, which is how you stop the method from “forecasting” a point by interpolating between its immediate neighbours in time.

A population that is deterministic and unpredictable

The Ricker map is the standard single-species chaotic population model: density dependence with overcompensation, no external forcing, no process noise. At an intrinsic growth rate of 3.0 it is chaotic. Every value follows exactly from the one before, and yet the series looks like noise.

We measure it with 10 per cent lognormal observation error, which is optimistic for a real count survey and still enough to matter.

gen_ricker <- function(n, burn = 500, r = 3.0, x0 = 0.6) {
  N <- n + burn; x <- numeric(N); x[1] <- x0
  for (t in 1:(N - 1)) x[t + 1] <- x[t] * exp(r * (1 - x[t]))
  x[(burn + 1):N]
}

set.seed(4264)
n <- 300
z_true <- gen_ricker(n)
z <- z_true * exp(rnorm(n, 0, 0.10))     # 10% observation error

acf1 <- acf(z, plot = FALSE)$acf[2]

The lag-1 autocorrelation is -0.71. The comparison series is the honest linear counterpart: fit an autoregressive model to the chaotic series, then simulate from it. It has the same second-order structure by construction, and no determinism whatsoever.

ar_fit <- ar(z, order.max = 10, method = "mle")
set.seed(9264)
z_lin <- as.numeric(arima.sim(n = n, model = list(ar = ar_fit$ar),
                              sd = sqrt(ar_fit$var.pred)))
acf1_lin <- acf(z_lin, plot = FALSE)$acf[2]

An AR(6) was selected, and the simulated series has lag-1 autocorrelation -0.75 against the original’s -0.71. To the autocorrelation function these two series are twins.

The picture that gives the game away

Plot each value against the one before it. This is the two-dimensional delay embedding, and for a one-dimensional map it is the map.

Two scatter plots. The left shows values against their previous value for a chaotic population, forming a clean single-humped curve. The right shows the same plot for a linear autoregressive series, forming a broad shapeless cloud with a weak negative slope. — Figure 1: Lag plots of two series with the same autocorrelation function. The chaotic population traces its own deterministic map; the linear series with matching autocorrelation is a cloud.

Both series have the same lag-1 autocorrelation. One of them is a curve with a little observation scatter around it, the other is a cloud. Every summary that stops at second-order structure treats them as equivalent; the lag plot does not, and neither does simplex projection.

For a real system with more than one dimension the picture is not this generous, which is exactly why the forecast skill, not the eyeball, has to do the work.

Choosing the embedding dimension

Run simplex at E = 1 to 8 and keep the E that forecasts best. For a one-dimensional map like the Ricker the true answer is E = 1.

Es <- 1:8
rho_E <- sapply(Es, function(E) simplex(z, E = E, tp = 1)$rho)
rho_E_lin <- sapply(Es, function(E) simplex(z_lin, E = E, tp = 1)$rho)
E_best <- Es[which.max(rho_E)]
E_best_lin <- Es[which.max(rho_E_lin)]

Two line plots of forecast correlation against embedding dimension from one to eight. The left panel for the chaotic series is high near 0.96 and declines slowly. The right panel for the linear series rises from about 0.62 to a peak near dimension four at 0.76 and then falls. — Figure 2: One-step forecast skill against embedding dimension. The chaotic series peaks at its true dimension of one; the linear series also has a best E, which measures nothing.

The chaotic series peaks at E = 2 with rho = 0.962, which is the right answer. The linear series peaks at E = 4, and that number means nothing at all: the series has no attractor and no dimension. A peak in this curve is not evidence of determinism, it is just the point where the bias and variance of a nearest-neighbour average happen to balance. The peak location is the first thing people over-read in EDM, and the first honest limit of this post.

The test: does the nonlinear forecast beat the linear one?

The comparison that carries information is not simplex against nothing, it is simplex against the best linear model you can build from the same data. If a linear autoregressive model does as well, the reconstruction has bought you nothing.

We give the linear model an advantage on purpose. Its coefficients are fitted on the whole series and it predicts in sample, while simplex leaves each point out. If simplex still wins, the result is not an artefact of a generous comparison.

lin_rho <- function(x, p, tp) {
  n <- length(x)
  X <- embed(x, p + 1)
  cf <- coef(lm(X[, 1] ~ X[, -1, drop = FALSE]))
  pr <- rep(NA_real_, n)
  for (t in p:(n - tp)) {
    st <- x[t:(t - p + 1)]
    for (k in 1:tp) {                       # iterate the linear map forward
      nx <- cf[1] + sum(cf[-1] * st)
      st <- c(nx, st[-p])
    }
    pr[t + tp] <- st[1]
  }
  ok <- !is.na(pr)
  cor(pr[ok], x[ok])
}

tps <- 1:8
sx_chaos <- sapply(tps, function(tp) simplex(z, E = E_best, tp = tp)$rho)
ln_chaos <- sapply(tps, function(tp) lin_rho(z, ar_fit$order, tp))

ar_lin <- ar(z_lin, order.max = 10, method = "mle")
sx_lin <- sapply(tps, function(tp) simplex(z_lin, E = E_best_lin, tp = tp)$rho)
ln_lin <- sapply(tps, function(tp) lin_rho(z_lin, ar_lin$order, tp))

Two line plots of forecast correlation against horizon one to eight. The left panel for the chaotic series shows simplex starting near 0.96 above the linear model near 0.83, both declining, with simplex above throughout. The right panel for the linear series shows the linear model above simplex at every horizon. — Figure 3: Forecast skill against horizon for the nonlinear and linear predictors. On the chaotic series simplex wins and then decays as trajectories diverge; on the linear series the linear model wins at every horizon.

On the chaotic series simplex reaches rho = 0.962 one step ahead against the linear model’s 0.816, and the gap widens at two steps (0.917 against 0.449). On the linear series the ordering flips: simplex gets 0.691 and the linear model 0.820. The nonlinear machinery loses, as it should, because there is nothing nonlinear to find.

Two separate things are visible in the chaotic panel. The gap between the curves is the determinism: state-space structure that a linear model cannot represent. The decline of the simplex curve, from 0.962 at one step to 0.223 at eight, is chaos itself: nearby trajectories separate exponentially, so knowing the state now buys less and less about the state later. A deterministic system can be perfectly knowable and still unforecastable beyond a few steps. That decay rate is a rough read on the Lyapunov exponent, and it is the thing that makes “deterministic” and “predictable” different words.

How reliable is the verdict?

The gap is a difference of two correlations, not a test statistic. It has no null distribution attached, so the sensible question is how often it points the wrong way. Sixty replicate series each of the linear process, and thirty of the chaotic one:

gain <- function(x, Emax = 6) {
  rr <- sapply(1:Emax, function(E) simplex(x, E = E, tp = 1)$rho)
  p <- max(1, ar(x, order.max = 6, method = "mle")$order)
  X <- embed(x, p + 1)
  max(rr) - cor(fitted(lm(X[, 1] ~ X[, -1, drop = FALSE])), X[, 1])
}

set.seed(3264)
g_null <- replicate(60, gain(as.numeric(arima.sim(n = 200, model = list(ar = ar_fit$ar),
                                                  sd = sqrt(ar_fit$var.pred)))))
g_chaos <- replicate(30, gain(gen_ricker(200, x0 = runif(1, 0.3, 0.9)) *
                              exp(rnorm(200, 0, 0.10))))
fp <- mean(g_null > 0)
tp_rate <- mean(g_chaos > 0)

The gain is positive in 0% of the linear replicates and 100% of the chaotic ones. The direction is trustworthy here, but notice what is doing the work: simplex is handicapped by leave-one-out while the linear model is not, so the comparison leans against a false positive by construction. That is a design choice, not a calibrated size. If you want a p-value, you need a null that preserves the linear structure and destroys the determinism, which is what the phase-randomised surrogate does, and what the checking post uses.

Where this breaks

The best E moves with the noise. Repeat the sweep at three observation-error levels:

set.seed(5264)
noise_levels <- c(0, 0.05, 0.10, 0.20)
E_by_noise <- sapply(noise_levels, function(s) {
  zz <- z_true * exp(rnorm(n, 0, s))
  which.max(sapply(1:8, function(E) simplex(zz, E = E, tp = 1)$rho))
})
rho_flat <- diff(range(rho_E[1:4]))

Left panel: a step plot of selected embedding dimension against observation error showing the choice drifting upward as noise increases. Right panel: forecast skill against embedding dimension for four noise levels, showing nearly flat curves whose maxima are barely distinguishable. — Figure 4: Selected embedding dimension against observation error, and the flatness of the skill curve that makes the selection fragile.

With no observation error the sweep returns E = 1, the truth. Add error and the choice climbs to 1, 2, 2 at 5, 10, 20 per cent. The true dimension never changed.

Look closely at the 10 per cent case and it gets worse. The main analysis above, on a series with exactly that much error, selected E = 2; this sweep, at the same error level but a different random draw, selects E = 2. Same system, same noise level, different answer. Across E = 1 to 4 the skill on the original series spans only 0.010, so the argmax is picked out of a nearly flat curve by whichever way the noise happened to fall. Report the curve, not just its maximum, and check that your conclusions survive the neighbouring values of E.

Temporal neighbours are not state-space neighbours. With a smooth, densely sampled series the closest point in delay space is often the one immediately before or after in time. Predicting from those is interpolation dressed up as forecasting, and the skill it reports is not out-of-sample in any useful sense. The theiler argument exists for this; set it to something like the decorrelation time and see whether the skill survives. It matters most for high-frequency data (loggers, camera traps, daily counts) and least for annual series, where consecutive years are already far apart on the attractor.

Determinism is not the same as low dimension, and neither is the same as knowing the mechanism. A good simplex forecast says the system has state-dependent structure that a linear model misses. It does not say which variables make up the state, whether the dynamics are stationary, or what causes what. The S-map turns the first question into a quantitative test, and convergent cross mapping takes on the last one.

Where to go next

The reconstruction is the foundation for everything else in this cluster. Once you trust that M_x is a faithful copy of the attractor, two things follow. First, the local structure can be turned into a smoothly tunable model that measures how nonlinear the dynamics are, which is the S-map. Second, if two species are dynamically coupled, each one’s reconstruction carries the other’s signature, which is what makes causal inference from time series possible at all.

Before applying any of it to a real series, be clear about what you have. Thirty annual counts are not going to reconstruct anything: the library is too small for nearest neighbours to be near. The methods here need hundreds of points, or a great deal of luck.

References

Takens 1981 Lecture Notes in Mathematics 898:366-381 (10.1007/BFb0091924)

Sugihara & May 1990 Nature 344(6268):734-741 (10.1038/344734a0)

Sugihara, May, Ye, Hsieh, Deyle, Fogarty & Munch 2012 Science 338(6106):496-500 (10.1126/science.1227079)

Hsieh, Glaser, Lucas & Sugihara 2005 Nature 435(7040):336-340 (10.1038/nature03553)

Deyle & Sugihara 2011 PLoS ONE 6(3):e18295 (10.1371/journal.pone.0018295)

Chang, Ushio & Hsieh 2017 Ecological Research 32(6):785-796 (10.1007/s11284-017-1469-9)

Theiler 1986 Physical Review A 34(3):2427-2432 (10.1103/PhysRevA.34.2427)

--- title: "Simplex projection and delay embedding" description: "Reconstruct a population's attractor from a single time series, forecast it with simplex projection in base R, and separate deterministic dynamics from red noise." date: "2026-08-22 09:00" categories: [R, EDM, time series, nonlinear dynamics, ecology tutorial] image: thumbnail.png image-alt: "Lag plot of a chaotic population series showing a clean deterministic curve beside a structureless cloud from linear noise." --- An annual count series wobbles up and down for thirty years. The autocorrelation function says the wobble is not white noise, and an ARIMA model fits it well enough. Nothing in that workflow can tell you whether the wobble comes from a low-dimensional deterministic system that you could, in principle, forecast, or from a linear process pushed around by weather. Empirical dynamic modelling asks the question the other way round. Instead of fitting an equation, it treats the series as a shadow of a trajectory moving on an attractor, and reconstructs that attractor from the series itself. This post builds the two pieces you need: the delay embedding, and simplex projection, the nearest-neighbour forecast that runs on it. Both are twenty lines of base R. The payoff is a test that no linear model can perform: does knowing where the system is *now*, in reconstructed state space, tell you more about the next step than the best linear predictor can? For the series below the answer is yes, and the same machinery says no for a linear series with an identical autocorrelation function. ## The idea: one variable carries the whole attractor A population does not live in one dimension. Abundance next year depends on abundance now, on the predator, on the nutrient pool, on whatever else is coupled to it. That full state is a point moving on an attractor in some multi-dimensional space, and you measured one coordinate of it. Takens showed that this is enough. Under mild conditions, the vector of lagged values $$\mathbf{v}_t = \big(x_t,\; x_{t-\tau},\; x_{t-2\tau},\; \ldots,\; x_{t-(E-1)\tau}\big)$$ traces out a shape that is a faithful copy of the true attractor: distinct states stay distinct, and nearby points stay nearby. The unmeasured variables leave their fingerprint in the history of the one you did measure. `E` is the embedding dimension and `tau` the lag; for annual ecological series `tau = 1` is almost always the right choice, so the only real knob is `E`. The practical consequence is the one that matters here. If the reconstruction is faithful, then points that are close in delay coordinates are close in the true state, so they should do the same thing next. That turns forecasting into a lookup: find the neighbours, see where they went, average. ## The engine Two functions. `make_block` stacks the delay vectors; `simplex` finds the `E+1` nearest neighbours of each point, weights them by distance, and predicts `tp` steps ahead. The prediction for each point leaves that point out, so nothing is fitted to itself. ```{r} #| label: edm-core make_block <- function(x, E, tau = 1) { n <- length(x) tt <- seq((E - 1) * tau + 1, n) M <- matrix(NA_real_, nrow = length(tt), ncol = E) for (j in seq_len(E)) M[, j] <- x[tt - (j - 1) * tau] list(M = M, t = tt) } simplex <- function(x, E, tau = 1, tp = 1, theiler = 0, target = NULL) { if (is.null(target)) target <- x bl <- make_block(x, E, tau) M <- bl$M; tt <- bl$t ok <- (tt + tp) >= 1 & (tt + tp) <= length(target) M <- M[ok, , drop = FALSE]; tt <- tt[ok] y <- target[tt + tp] D <- as.matrix(dist(M)) pred <- rep(NA_real_, length(tt)) for (i in seq_along(tt)) { d <- D[i, ] d[abs(tt - tt[i]) <= theiler] <- Inf # drop self and temporal neighbours o <- order(d)[seq_len(E + 1)] dn <- d[o] if (!all(is.finite(dn))) next w <- if (dn[1] == 0) as.numeric(dn == 0) else exp(-dn / dn[1]) w <- pmax(w, 1e-8) pred[i] <- sum(w * y[o]) / sum(w) } keep <- !is.na(pred) list(rho = suppressWarnings(cor(pred[keep], y[keep])), pred = pred, obs = y, t = tt) } ``` The weights `exp(-d / d1)` scale each neighbour by its distance relative to the closest one, so an identical past state dominates and a far-off state barely contributes. The simplex of `E+1` points is the smallest set that can enclose a point in `E` dimensions, which is where the method gets its name. `theiler` deserves a note now and a warning later. Setting it to zero drops only the point being predicted. Setting it higher drops every library point within that many time steps, which is how you stop the method from "forecasting" a point by interpolating between its immediate neighbours in time. ## A population that is deterministic and unpredictable The Ricker map is the standard single-species chaotic population model: density dependence with overcompensation, no external forcing, no process noise. At an intrinsic growth rate of 3.0 it is chaotic. Every value follows exactly from the one before, and yet the series looks like noise. We measure it with 10 per cent lognormal observation error, which is optimistic for a real count survey and still enough to matter. ```{r} #| label: data gen_ricker <- function(n, burn = 500, r = 3.0, x0 = 0.6) { N <- n + burn; x <- numeric(N); x[1] <- x0 for (t in 1:(N - 1)) x[t + 1] <- x[t] * exp(r * (1 - x[t])) x[(burn + 1):N] } set.seed(4264) n <- 300 z_true <- gen_ricker(n) z <- z_true * exp(rnorm(n, 0, 0.10)) # 10% observation error acf1 <- acf(z, plot = FALSE)$acf[2] ``` The lag-1 autocorrelation is `r sprintf("%.2f", acf1)`. The comparison series is the honest linear counterpart: fit an autoregressive model to the chaotic series, then simulate from it. It has the same second-order structure by construction, and no determinism whatsoever. ```{r} #| label: surrogate ar_fit <- ar(z, order.max = 10, method = "mle") set.seed(9264) z_lin <- as.numeric(arima.sim(n = n, model = list(ar = ar_fit$ar), sd = sqrt(ar_fit$var.pred))) acf1_lin <- acf(z_lin, plot = FALSE)$acf[2] ``` An AR(`r ar_fit$order`) was selected, and the simulated series has lag-1 autocorrelation `r sprintf("%.2f", acf1_lin)` against the original's `r sprintf("%.2f", acf1)`. To the autocorrelation function these two series are twins. ```{r} #| label: setup #| echo: false library(ggplot2) library(grid) te_ink <- "#16241d"; te_forest <- "#275139"; te_gold <- "#c9b458" te_rust <- "#b5534e"; te_paper <- "#f5f4ee"; te_line <- "#dad9ca" te_muted <- "#5d6b61"; te_green <- "#2f8f63" theme_te <- function(base_size = 11) { theme_minimal(base_size = base_size) + theme( plot.title = element_text(colour = te_ink, face = "bold", size = rel(1.0)), plot.subtitle = element_text(colour = te_muted, size = rel(0.85)), axis.title = element_text(colour = te_ink, size = rel(0.9)), axis.text = element_text(colour = te_muted, size = rel(0.8)), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), panel.background = element_rect(fill = te_paper, colour = NA), plot.background = element_rect(fill = te_paper, colour = NA), legend.position = "bottom", legend.title = element_text(colour = te_ink, size = rel(0.8)), legend.text = element_text(colour = te_muted, size = rel(0.75)), strip.text = element_text(colour = te_ink, size = rel(0.85)) ) } two_panel <- function(pA, pB) { grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(pA, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) } ``` ## The picture that gives the game away Plot each value against the one before it. This is the two-dimensional delay embedding, and for a one-dimensional map it *is* the map. ```{r} #| label: fig-lagplot #| fig-width: 8 #| fig-height: 4 #| fig-cap: "Lag plots of two series with the same autocorrelation function. The chaotic population traces its own deterministic map; the linear series with matching autocorrelation is a cloud." #| fig-alt: "Two scatter plots. The left shows values against their previous value for a chaotic population, forming a clean single-humped curve. The right shows the same plot for a linear autoregressive series, forming a broad shapeless cloud with a weak negative slope." #| echo: false lagdf <- function(v, lab) data.frame(prev = v[-length(v)], now = v[-1], panel = lab) dA <- lagdf(z, "chaotic population") dB <- lagdf(z_lin, "linear series, same ACF") pA <- ggplot(dA, aes(prev, now)) + geom_point(colour = te_forest, alpha = 0.7, size = 1.1) + labs(title = "Chaotic population", subtitle = "deterministic skeleton visible", x = expression(x[t-1]), y = expression(x[t])) + theme_te() pB <- ggplot(dB, aes(prev, now)) + geom_point(colour = te_rust, alpha = 0.7, size = 1.1) + labs(title = "Linear series, same ACF", subtitle = "no skeleton", x = expression(x[t-1]), y = expression(x[t])) + theme_te() two_panel(pA, pB) ``` Both series have the same lag-1 autocorrelation. One of them is a curve with a little observation scatter around it, the other is a cloud. Every summary that stops at second-order structure treats them as equivalent; the lag plot does not, and neither does simplex projection. For a real system with more than one dimension the picture is not this generous, which is exactly why the forecast skill, not the eyeball, has to do the work. ## Choosing the embedding dimension Run simplex at `E = 1` to `8` and keep the `E` that forecasts best. For a one-dimensional map like the Ricker the true answer is `E = 1`. ```{r} #| label: e-sweep Es <- 1:8 rho_E <- sapply(Es, function(E) simplex(z, E = E, tp = 1)$rho) rho_E_lin <- sapply(Es, function(E) simplex(z_lin, E = E, tp = 1)$rho) E_best <- Es[which.max(rho_E)] E_best_lin <- Es[which.max(rho_E_lin)] ``` ```{r} #| label: fig-esweep #| fig-width: 8 #| fig-height: 4 #| fig-cap: "One-step forecast skill against embedding dimension. The chaotic series peaks at its true dimension of one; the linear series also has a best E, which measures nothing." #| fig-alt: "Two line plots of forecast correlation against embedding dimension from one to eight. The left panel for the chaotic series is high near 0.96 and declines slowly. The right panel for the linear series rises from about 0.62 to a peak near dimension four at 0.76 and then falls." #| echo: false eA <- data.frame(E = Es, rho = rho_E) eB <- data.frame(E = Es, rho = rho_E_lin) pA <- ggplot(eA, aes(E, rho)) + geom_line(colour = te_forest, linewidth = 0.7) + geom_point(colour = te_forest, size = 2) + geom_point(data = eA[eA$E == E_best, ], colour = te_gold, size = 3.5) + scale_x_continuous(breaks = Es) + labs(title = "Chaotic population", subtitle = paste0("best E = ", E_best, " (true dimension 1)"), x = "embedding dimension E", y = expression(rho)) + theme_te() pB <- ggplot(eB, aes(E, rho)) + geom_line(colour = te_rust, linewidth = 0.7) + geom_point(colour = te_rust, size = 2) + geom_point(data = eB[eB$E == E_best_lin, ], colour = te_gold, size = 3.5) + scale_x_continuous(breaks = Es) + labs(title = "Linear series, same ACF", subtitle = paste0("best E = ", E_best_lin, " (no true dimension)"), x = "embedding dimension E", y = expression(rho)) + theme_te() two_panel(pA, pB) ``` The chaotic series peaks at `E = `r E_best`` with `rho = `r sprintf("%.3f", max(rho_E))``, which is the right answer. The linear series peaks at `E = `r E_best_lin``, and that number means nothing at all: the series has no attractor and no dimension. A peak in this curve is not evidence of determinism, it is just the point where the bias and variance of a nearest-neighbour average happen to balance. The peak location is the first thing people over-read in EDM, and the first honest limit of this post. ## The test: does the nonlinear forecast beat the linear one? The comparison that carries information is not simplex against nothing, it is simplex against the best linear model you can build from the same data. If a linear autoregressive model does as well, the reconstruction has bought you nothing. We give the linear model an advantage on purpose. Its coefficients are fitted on the whole series and it predicts in sample, while simplex leaves each point out. If simplex still wins, the result is not an artefact of a generous comparison. ```{r} #| label: horizon lin_rho <- function(x, p, tp) { n <- length(x) X <- embed(x, p + 1) cf <- coef(lm(X[, 1] ~ X[, -1, drop = FALSE])) pr <- rep(NA_real_, n) for (t in p:(n - tp)) { st <- x[t:(t - p + 1)] for (k in 1:tp) { # iterate the linear map forward nx <- cf[1] + sum(cf[-1] * st) st <- c(nx, st[-p]) } pr[t + tp] <- st[1] } ok <- !is.na(pr) cor(pr[ok], x[ok]) } tps <- 1:8 sx_chaos <- sapply(tps, function(tp) simplex(z, E = E_best, tp = tp)$rho) ln_chaos <- sapply(tps, function(tp) lin_rho(z, ar_fit$order, tp)) ar_lin <- ar(z_lin, order.max = 10, method = "mle") sx_lin <- sapply(tps, function(tp) simplex(z_lin, E = E_best_lin, tp = tp)$rho) ln_lin <- sapply(tps, function(tp) lin_rho(z_lin, ar_lin$order, tp)) ``` ```{r} #| label: fig-horizon #| fig-width: 8 #| fig-height: 4 #| fig-cap: "Forecast skill against horizon for the nonlinear and linear predictors. On the chaotic series simplex wins and then decays as trajectories diverge; on the linear series the linear model wins at every horizon." #| fig-alt: "Two line plots of forecast correlation against horizon one to eight. The left panel for the chaotic series shows simplex starting near 0.96 above the linear model near 0.83, both declining, with simplex above throughout. The right panel for the linear series shows the linear model above simplex at every horizon." #| echo: false hA <- rbind(data.frame(tp = tps, rho = sx_chaos, method = "simplex"), data.frame(tp = tps, rho = ln_chaos, method = "linear AR")) hB <- rbind(data.frame(tp = tps, rho = sx_lin, method = "simplex"), data.frame(tp = tps, rho = ln_lin, method = "linear AR")) pal <- c("simplex" = te_forest, "linear AR" = te_rust) pA <- ggplot(hA, aes(tp, rho, colour = method)) + geom_line(linewidth = 0.7) + geom_point(size = 1.8) + scale_colour_manual(values = pal, name = NULL) + scale_x_continuous(breaks = tps) + labs(title = "Chaotic population", subtitle = "determinism pays, then decays", x = "forecast horizon (steps)", y = expression(rho)) + theme_te() pB <- ggplot(hB, aes(tp, rho, colour = method)) + geom_line(linewidth = 0.7) + geom_point(size = 1.8) + scale_colour_manual(values = pal, name = NULL) + scale_x_continuous(breaks = tps) + labs(title = "Linear series, same ACF", subtitle = "nothing to gain", x = "forecast horizon (steps)", y = expression(rho)) + theme_te() two_panel(pA, pB) ``` On the chaotic series simplex reaches `rho = `r sprintf("%.3f", sx_chaos[1])`` one step ahead against the linear model's `r sprintf("%.3f", ln_chaos[1])`, and the gap widens at two steps (`r sprintf("%.3f", sx_chaos[2])` against `r sprintf("%.3f", ln_chaos[2])`). On the linear series the ordering flips: simplex gets `r sprintf("%.3f", sx_lin[1])` and the linear model `r sprintf("%.3f", ln_lin[1])`. The nonlinear machinery loses, as it should, because there is nothing nonlinear to find. Two separate things are visible in the chaotic panel. The *gap* between the curves is the determinism: state-space structure that a linear model cannot represent. The *decline* of the simplex curve, from `r sprintf("%.3f", sx_chaos[1])` at one step to `r sprintf("%.3f", sx_chaos[8])` at eight, is chaos itself: nearby trajectories separate exponentially, so knowing the state now buys less and less about the state later. A deterministic system can be perfectly knowable and still unforecastable beyond a few steps. That decay rate is a rough read on the Lyapunov exponent, and it is the thing that makes "deterministic" and "predictable" different words. ## How reliable is the verdict? The gap is a difference of two correlations, not a test statistic. It has no null distribution attached, so the sensible question is how often it points the wrong way. Sixty replicate series each of the linear process, and thirty of the chaotic one: ```{r} #| label: mc-null gain <- function(x, Emax = 6) { rr <- sapply(1:Emax, function(E) simplex(x, E = E, tp = 1)$rho) p <- max(1, ar(x, order.max = 6, method = "mle")$order) X <- embed(x, p + 1) max(rr) - cor(fitted(lm(X[, 1] ~ X[, -1, drop = FALSE])), X[, 1]) } set.seed(3264) g_null <- replicate(60, gain(as.numeric(arima.sim(n = 200, model = list(ar = ar_fit$ar), sd = sqrt(ar_fit$var.pred))))) g_chaos <- replicate(30, gain(gen_ricker(200, x0 = runif(1, 0.3, 0.9)) * exp(rnorm(200, 0, 0.10)))) fp <- mean(g_null > 0) tp_rate <- mean(g_chaos > 0) ``` The gain is positive in `r sprintf("%.0f%%", 100 * fp)` of the linear replicates and `r sprintf("%.0f%%", 100 * tp_rate)` of the chaotic ones. The direction is trustworthy here, but notice what is doing the work: simplex is handicapped by leave-one-out while the linear model is not, so the comparison leans against a false positive by construction. That is a design choice, not a calibrated size. If you want a p-value, you need a null that preserves the linear structure and destroys the determinism, which is what the phase-randomised surrogate does, and what the [checking post](../checking-an-empirical-dynamic-model/) uses. ## Where this breaks **The best E moves with the noise.** Repeat the sweep at three observation-error levels: ```{r} #| label: noise-e set.seed(5264) noise_levels <- c(0, 0.05, 0.10, 0.20) E_by_noise <- sapply(noise_levels, function(s) { zz <- z_true * exp(rnorm(n, 0, s)) which.max(sapply(1:8, function(E) simplex(zz, E = E, tp = 1)$rho)) }) rho_flat <- diff(range(rho_E[1:4])) ``` ```{r} #| label: fig-noise #| fig-width: 8 #| fig-height: 4 #| fig-cap: "Selected embedding dimension against observation error, and the flatness of the skill curve that makes the selection fragile." #| fig-alt: "Left panel: a step plot of selected embedding dimension against observation error showing the choice drifting upward as noise increases. Right panel: forecast skill against embedding dimension for four noise levels, showing nearly flat curves whose maxima are barely distinguishable." #| echo: false nA <- data.frame(noise = noise_levels, E = E_by_noise) pA <- ggplot(nA, aes(noise, E)) + geom_step(colour = te_forest, linewidth = 0.8) + geom_point(colour = te_gold, size = 3) + scale_y_continuous(breaks = 1:8, limits = c(0.5, 8.5)) + labs(title = "Selected E drifts with noise", subtitle = "true dimension is 1 throughout", x = "observation error (lognormal sd)", y = "selected E") + theme_te() set.seed(5264) curves <- do.call(rbind, lapply(noise_levels, function(s) { zz <- z_true * exp(rnorm(n, 0, s)) data.frame(E = Es, rho = sapply(Es, function(E) simplex(zz, E = E, tp = 1)$rho), noise = factor(s)) })) pB <- ggplot(curves, aes(E, rho, colour = noise)) + geom_line(linewidth = 0.6) + geom_point(size = 1.4) + scale_colour_manual(values = c("0" = te_ink, "0.05" = te_forest, "0.1" = te_green, "0.2" = te_gold), name = "obs sd") + scale_x_continuous(breaks = Es) + labs(title = "The curves are flat", subtitle = "the maximum is a weak signal", x = "embedding dimension E", y = expression(rho)) + theme_te() two_panel(pA, pB) ``` With no observation error the sweep returns `E = `r E_by_noise[1]``, the truth. Add error and the choice climbs to `r paste(E_by_noise[-1], collapse = ", ")` at `r paste(noise_levels[-1] * 100, collapse = ", ")` per cent. The true dimension never changed. Look closely at the 10 per cent case and it gets worse. The main analysis above, on a series with exactly that much error, selected `E = `r E_best``; this sweep, at the same error level but a different random draw, selects `E = `r E_by_noise[3]``. Same system, same noise level, different answer. Across `E = 1` to `4` the skill on the original series spans only `r sprintf("%.3f", rho_flat)`, so the argmax is picked out of a nearly flat curve by whichever way the noise happened to fall. Report the curve, not just its maximum, and check that your conclusions survive the neighbouring values of `E`. **Temporal neighbours are not state-space neighbours.** With a smooth, densely sampled series the closest point in delay space is often the one immediately before or after in time. Predicting from those is interpolation dressed up as forecasting, and the skill it reports is not out-of-sample in any useful sense. The `theiler` argument exists for this; set it to something like the decorrelation time and see whether the skill survives. It matters most for high-frequency data (loggers, camera traps, daily counts) and least for annual series, where consecutive years are already far apart on the attractor. **Determinism is not the same as low dimension, and neither is the same as knowing the mechanism.** A good simplex forecast says the system has state-dependent structure that a linear model misses. It does not say which variables make up the state, whether the dynamics are stationary, or what causes what. The [S-map](../s-map-and-state-dependence/) turns the first question into a quantitative test, and [convergent cross mapping](../convergent-cross-mapping/) takes on the last one. ## Where to go next The reconstruction is the foundation for everything else in this cluster. Once you trust that `M_x` is a faithful copy of the attractor, two things follow. First, the local structure can be turned into a smoothly tunable model that measures *how* nonlinear the dynamics are, which is the S-map. Second, if two species are dynamically coupled, each one's reconstruction carries the other's signature, which is what makes causal inference from time series possible at all. Before applying any of it to a real series, be clear about what you have. Thirty annual counts are not going to reconstruct anything: the library is too small for nearest neighbours to be near. The methods here need hundreds of points, or a great deal of luck. ## References Takens 1981 Lecture Notes in Mathematics 898:366-381 (10.1007/BFb0091924) Sugihara & May 1990 Nature 344(6268):734-741 (10.1038/344734a0) Sugihara, May, Ye, Hsieh, Deyle, Fogarty & Munch 2012 Science 338(6106):496-500 (10.1126/science.1227079) Hsieh, Glaser, Lucas & Sugihara 2005 Nature 435(7040):336-340 (10.1038/nature03553) Deyle & Sugihara 2011 PLoS ONE 6(3):e18295 (10.1371/journal.pone.0018295) Chang, Ushio & Hsieh 2017 Ecological Research 32(6):785-796 (10.1007/s11284-017-1469-9) Theiler 1986 Physical Review A 34(3):2427-2432 (10.1103/PhysRevA.34.2427) ## Related tutorials - [The S-map and state-dependent dynamics](../s-map-and-state-dependence/) - [Convergent cross mapping in ecology](../convergent-cross-mapping/) - [Checking an empirical dynamic model](../checking-an-empirical-dynamic-model/) - [Autocorrelation in time series](../autocorrelation-in-time-series/)