Empirical dynamic modeling

Empirical dynamic modeling (EDM) is a framework for analysis and prediction of nonlinear dynamical systems. Applications include population dynamics,¹²³⁴⁵⁶ ecosystem service,⁷ medicine,⁸ neuroscience,⁹¹⁰¹¹ dynamical systems,¹²¹³¹⁴ geophysics,¹⁵¹⁶¹⁷ and human-computer interaction.¹⁸ EDM was originally developed by Robert May and George Sugihara. It can be considered a methodology for data modeling, predictive analytics, dynamical system analysis, machine learning and time series analysis.

Mathematical models have tremendous power to describe observations of real-world systems. They are routinely used to test hypothesis, explain mechanisms and predict future outcomes. However, real-world systems are often nonlinear and multidimensional, in some instances rendering explicit equation-based modeling problematic. Empirical models, which infer patterns and associations from the data instead of using hypothesized equations, represent a natural and flexible framework for modeling complex dynamics.

Donald DeAngelis and Simeon Yurek illustrated that canonical statistical models are ill-posed when applied to nonlinear dynamical systems.¹⁹ A hallmark of nonlinear dynamics is state-dependence: system states are related to previous states governing transition from one state to another. EDM operates in this space, the multidimensional state-space of system dynamics rather than on one-dimensional observational time series. EDM does not presume relationships among states, for example, a functional dependence, but projects future states from localised, neighboring states. EDM is thus a state-space, nearest-neighbors paradigm where system dynamics are inferred from states derived from observational time series. This provides a model-free representation of the system naturally encompassing nonlinear dynamics.

A cornerstone of EDM is recognition that time series observed from a dynamical system can be transformed into higher-dimensional state-spaces by time-delay embedding with Takens's theorem. The state-space models are evaluated based on in-sample fidelity to observations, conventionally with Pearson correlation between predictions and observations.

Primary EDM algorithms include Simplex projection,²⁰ Sequential locally weighted global linear maps (S-Map) projection,²¹ Multivariate embedding in Simplex or S-Map,¹ Convergent cross mapping (CCM),²² and Multiview Embeding,²³ described below.

Nomenclature

Parameter	Description
${\displaystyle E}$	embedding dimension
${\displaystyle k}$	number of nearest neighbors
${\displaystyle T_{p}}$	prediction interval
${\displaystyle X\in \mathbb {R} }$	observed time series
${\displaystyle y\in \mathbb {R} ^{E}}$	vector of lagged observations
${\displaystyle \theta \geq 0}$	S-Map localization
${\displaystyle X_{t}^{E}=(X_{t},X_{t-1},\dots ,X_{t-E+1})\in \mathbb {R} ^{E}}$	lagged embedding vectors
${\displaystyle \|v\|}$	norm of v
${\displaystyle N=\{N_{1},\dots ,N_{k}\}}$	list of nearest neighbors

Nearest neighbors are found according to: ${\displaystyle {\text{NN}}(y,X,k)=\|X_{N_{i}}^{E}-y\|\leq \|X_{N_{j}}^{E}-y\|{\text{ if }}1\leq i\leq j\leq k}$

Simplex projection^20242526 is a nearest neighbor projection. It locates the ${\displaystyle k}$ nearest neighbors to the location in the state-space from which a prediction is desired. To minimize the number of free parameters ${\displaystyle k}$ is typically set to ${\displaystyle E+1}$ defining an ${\displaystyle E+1}$ dimensional simplex in the state-space. The prediction is computed as the average of the weighted phase-space simplex projected ${\displaystyle Tp}$ points ahead. Each neighbor is weighted proportional to their distance to the projection origin vector in the state-space.

1. Find ${\displaystyle k}$ nearest neighbor: ${\displaystyle N_{k}\gets {\text{NN}}(y,X,k)}$
2. Define the distance scale: ${\displaystyle d\gets \|X_{N_{1}}^{E}-y\|}$
3. Compute weights: For{${\displaystyle i=1,\dots ,k}$} : ${\displaystyle w_{i}\gets \exp(-\|X_{N_{i}}^{E}-y\|/d)}$
4. Average of state-space simplex: ${\displaystyle {\hat {y}}\gets \sum _{i=1}^{k}\left(w_{i}X_{N_{i}+T_{p}}\right)/\sum _{i=1}^{k}w_{i}}$

S-Map²¹ extends the state-space prediction in Simplex from an average of the ${\displaystyle E+1}$ nearest neighbors to a linear regression fit to all neighbors, but localised with an exponential decay kernel. The exponential localisation function is ${\displaystyle F(\theta )={\text{exp}}(-\theta d/D)}$, where ${\displaystyle d}$ is the neighbor distance and ${\displaystyle D}$ the mean distance. In this way, depending on the value of ${\displaystyle \theta }$, neighbors close to the prediction origin point have a higher weight than those further from it, such that a local linear approximation to the nonlinear system is reasonable. This localisation ability allows one to identify an optimal local scale, in-effect quantifying the degree of state dependence, and hence nonlinearity of the system.²⁷

Another feature of S-Map is that for a properly fit model, the regression coefficients between variables have been shown to approximate the gradient (directional derivative) of variables along the manifold.²⁸ These Jacobians represent the time-varying interaction strengths between system variables.

1. Find ${\displaystyle k}$ nearest neighbor: ${\displaystyle N\gets {\text{NN}}(y,X,k)}$
2. Sum of distances: ${\displaystyle D\gets {\frac {1}{k}}\sum _{i=1}^{k}\|X_{N_{i}}^{E}-y\|}$
3. Compute weights: For{${\displaystyle i=1,\dots ,k}$} : ${\displaystyle w_{i}\gets \exp(-\theta \|X_{N_{i}}^{E}-y\|/D)}$
4. Reweighting matrix: ${\displaystyle W\gets {\text{diag}}(w_{i})}$
5. Design matrix: ${\displaystyle A\gets {\begin{bmatrix}1&X_{N_{1}}&X_{N_{1}-1}&\dots &X_{N_{1}-E+1}\\1&X_{N_{2}}&X_{N_{2}-1}&\dots &X_{N_{2}-E+1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&X_{N_{k}}&X_{N_{k}-1}&\dots &X_{N_{k}-E+1}\end{bmatrix}}}$
6. Weighted design matrix: ${\displaystyle A\gets WA}$
7. Response vector at ${\displaystyle Tp}$: ${\displaystyle b\gets {\begin{bmatrix}X_{N_{1}+T_{p}}\\X_{N_{2}+T_{p}}\\\vdots \\X_{N_{k}+T_{p}}\end{bmatrix}}}$
8. Weighted response vector: ${\displaystyle b\gets Wb}$
9. Least squares solution (SVD): ${\displaystyle {\hat {c}}\gets {\text{argmin}}_{c}\|Ac-b\|_{2}^{2}}$
10. Local linear model ${\displaystyle {\hat {c}}}$ is prediction: ${\displaystyle {\hat {y}}\gets {\hat {c}}_{0}+\sum _{i=1}^{E}{\hat {c}}_{i}y_{i}}$

Multivariate Embedding¹¹²²⁹ recognizes that time-delay embeddings are not the only valid state-space construction. In Simplex and S-Map one can generate a state-space from observational vectors, or time-delay embeddings of a single observational time series, or both.

Convergent cross mapping (CCM)²² leverages a corollary to the Generalized Takens Theorem¹² that it should be possible to cross predict or cross map between variables observed from the same system. Suppose that in some dynamical system involving variables ${\displaystyle X}$ and ${\displaystyle Y}$, ${\displaystyle X}$ causes ${\displaystyle Y}$. Since ${\displaystyle X}$ and ${\displaystyle Y}$ belong to the same dynamical system, their reconstructions (via embeddings) ${\displaystyle M_{x}}$, and ${\displaystyle M_{y}}$, also map to the same system.

The causal variable ${\displaystyle X}$ leaves a signature on the affected variable ${\displaystyle Y}$, and consequently, the reconstructed states based on ${\displaystyle Y}$ can be used to cross predict values of ${\displaystyle X}$. CCM leverages this property to infer causality by predicting ${\displaystyle X}$ using the ${\displaystyle M_{y}}$ library of points (or vice versa for the other direction of causality), while assessing improvements in cross map predictability as larger and larger random samplings of ${\displaystyle M_{y}}$ are used. If the prediction skill of ${\displaystyle X}$ increases and saturates as the entire ${\displaystyle M_{y}}$ is used, this provides evidence that ${\displaystyle X}$ is casually influencing ${\displaystyle Y}$.

Multiview Embedding²³ is a Dimensionality reduction technique where a large number of state-space time series vectors are combitorially assessed towards maximal model predictability.

Extensions to EDM techniques include:

• Generalized Theorems for Nonlinear State Space Reconstruction¹²
• Extended Convergent Cross Mapping¹³
• Dynamic stability⁴
• S-Map regularization³⁰
• Visual analytics with EDM³¹
• Convergent Cross Sorting³²
• Expert system with EDM hybrid³³
• Sliding windows based on the extended convergent cross-mapping³⁴
• Empirical Mode Modeling¹⁷
• Accounting for missing data and variable step sizes³⁵
• Accounting for observation noise³⁶
• Hierarchical Bayesian EDM via Gaussian processes³⁷
• Intelligent and Adaptive Control³⁸
• Optimal control via Empirical dynamic programming³⁹
• Multiview distance regularised S-map⁴⁰