User Guide#

Empirical Dynamic Modeling#

Empirical dynamic modeling (EDM) is a framework for analysis and prediction of nonlinear dynamical systems with over 4,000 citations to its core algorithms: convergent cross mapping (CCM), simplex, and sequential locally weighted global linear maps (s-map). EDM continues to evolve with new applications and algorithmic extensions as documented in the EDM Wikipedia article. An in-depth introduction is provided in EDM docs.

`sciedm`#

sciedm provides a scikit-learn compliant implementation of core EDM algorithms along with utilities to discover appropriate embedding dimension and scale of state-dependence (nonlinearity). sciedm classes are based on the sklearn.base.BaseEstimator with Regressor and Transformer mixin’s defining the core sciedm classes:

Class	Function
`Simplex (RegressorMixin, BaseEstimator)`	simplex
`SMap (RegressorMixin, BaseEstimator)`	s-map
`CCM (TransformerMixin, BaseEstimator)`	CCM
`CCM_Matrix (TransformerMixin, BaseEstimator)`	CCM tensor
`EmbedDimension (TransformerMixin, BaseEstimator)`	embedding dimension
`PredictNonlinear (TransformerMixin, BaseEstimator)`	nonlinear dependence

sciedm leverages pandas DataFrame for I/O. Input data are expected to be a DataFrame with named columns including the target variable. The first column is expected to be a series of time or date values or strings, however this requirement can be removed with the noTime=True flag.

Example Data#

We use two data sets to illustrate sciedm. First a 5-dimensional coupled system generated from the Lorenz’96 model.

Lorenz'96 5-D — First 3 components of a Lorenz’96 5-D system.#

Second, cumulative daily water flow into Everglades National Park through the S12C, S12D and S333 spillways.

Daily cumulative flow into ENP. — Daily cumulative flow into Everglades National Park through S12C,D S333.#

Predictors#

sciedm predictors include Simplex and SMap.

Simplex#

Simplex is a nearest neighbor projection in the embedding from a convex simplex query state to a future (or past) state. It is the core of cross mapping, convergent cross mapping, and dimension estimation, inherently provides out of sample (train:test) validation with the lib and pred parameters, automatically excludes state vectors within the temporal prediction horizon and provides a temporal exclusionRadius parameter to exclude serially correlated states to focus on nonlinear dynamics.

Here we use Simplex to predict variable V3 from a 4-D multivariate embedding of [V1,V2,V4,V5] of the 5-D Lorenz’96 system. To specify a multivariate embedding instead of a time-delay embedding we set embedded=True which sets the embedding dimension to E=4. The default prediction horizon is Tp=1 points ahead.

>>> from sciedm import Simplex
>>> df = read_csv("../sciedm/data/Lorenz5D.csv")
>>> lib, pred = [1,500], [801,900] # out of sample library : prediction sets
>>> columns, target = ['V1','V2','V4','V5'], 'V3'
>>> smpx = Simplex(columns=columns, target=target, lib=lib, pred=pred, embedded=True)
>>> smpx.fit(df)
>>> rho = smpx.score(df, df[target])

>>> from sciedm.aux_func import PlotObsPred
>>> ax = PlotObsPred(smpx.Projection_,
...      title=f"Simplex: {smpx.columns} : {smpx.target} rho={rho:.2f}")

Simplex projection of Lorenz'96 5-D — Simplex prediction at `Tp=1` of Lorenz’96 5-D variable `V3` from a multivariate embedding of `['V1','V2','V4','V5']`#

Even though the embedding does not contain V3 the shared dynamical information in the other variables enables good prediction of V3 through simplex cross mapping.

SMap#

Sequential locally weighted global linear maps (s-map) can be viewed as a generalized forerunner of locally linear embedding. Instead of a fixed number of embedding neighbors an exponential localization kernel selects neighbors based on distance scale in the embedding allowing one to identify an optimal scale along system trajectories reflecting state dependence (nonlinearity). The scale is determined by a weight function: \(F(\theta )={\text{exp}}(-\theta d/D)\) where θ is the localization parameter, d a specific neighbor distance, and D the mean distance to all neighbors. At θ=0 all neighbors are equally weighted corresponding to the global linear map, if predictability is maximal at θ=0 there is not evidence for state dependence (nonlinearity).

When an s-map model is tuned to the scale appropriate for the dynamics, it is known s-map coefficients correspond to the time and state dependent derivatives between variables. As such, s-map facilitates prediction and quantification of both intervariable dependencies (Jacobians) and the scale of nonlinearity of multivariate dynamical systems.

Here we demonstrate SMap in multivariate time series prediction and variable interaction. The nearest neighbor localization parameter theta=8 weights local neighbors at small scale (distances) in the embedding.

>>> from sciedm import SMap
>>> smap = SMap(columns=columns, target=target, theta=8.,
...             embedded=True, lib=lib, pred=pred)
>>> smap.fit(df)
>>> rho = smap.score(df, df[target])

>>> from sciedm.aux_func import PlotObsPred, PlotCoeff
>>> title = f"SMap: {smap.columns} : {smap.target} rho={rho:.2f}"
>>> ax = PlotObsPred(smap.Projection_, title=title)
>>> ax = PlotCoeff(smap.Coefficients_, title=title)

SMap projection of Lorenz'96 5-D — SMap prediction at `Tp=1` of Lorenz’96 5-D variable `V3` from a multivariate embedding of `['V1','V2','V4','V5']` and intervariable derivatives (coefficients).#

Transformers#

Convergent cross mapping, embedding dimension estimation and assessment of nonlinear state dependence are implemented as Transformer classes.

CCM#

Convergent cross mapping (CCM) identifies whether two time series belong to the same dynamical system and are therefore causally related. The CCM Wikipedia page provides a detailed description.

Although we know the variables of the Lorenz’96 system are causally related, we use them as a demonstration.

>>> from sciedm import CCM
>>> df = read_csv("../sciedm/data/Lorenz5D.csv")
>>> libSizes = [20,50,100,200,500,900,1000]
>>> ccm = CCM(columns='V1', target='V5', E=5, libSizes=libSizes)
>>> ccm_V1_V5 = ccm.fit_transform(X=df)

>>> from sciedm.aux_func import PlotCCM
>>> ax = PlotCCM(ccm.libMeans_, title=f"E={ccm.E} {ccm.columns} : {ccm.target}")

CCM Lorenz'96 V1:V5 — CCM between Lorenz’96 V1 and V5.#

As expected both mappings V1:V5 (V5 drives V1) and V5:V1 (V1 drives V5) exhibit convergence and high predictability.

CCM_Matrix#

Given a set of N timeseries observations of M variables in an NxM DataFrame, CCM_Matrix computes the MxMxL CCM tensor of each variable against all others where L is the number of library sizes in libSizes.

Here we build a toy Lorenz’96 matrix as a demonstration.

>>> from sciedm import CCM_Matrix, PlotMatrix
>>> from pandas import concat, read_csv
>>> df_ = read_csv("../sciedm/data/Lorenz5D.csv")
>>> df20 = concat([df_.iloc[:,1:],df_.iloc[:,1:],df_.iloc[:,1:],df_.iloc[:,1:]], axis=1)
>>> df = concat( [df_, df20.sample(frac=1, axis=1)], axis=1 )
>>> df.columns = ['Time'] + [f"V{x}" for x in range(25)]
>>> libSizes = [50,100,500,900,1000]
>>> cmat = CCM_Matrix(E=5, libSizes=libSizes)
>>> tensor,columns = cmat.fit_transform(X=df)

>>> PlotMatrix(tensor[:,:,2], columns, title="Cross Map L=500 Lorenz'96")
>>> PlotMatrix(cmat.slope_, columns, title="Cross Map Slope Lorenz'96")
>>> ccm_L500 = tensor[:,:,2] # MxM cross map at L=500
>>> ccm_L500[ cmat.slope_ < 0.2 ] = np.nan # Mask library slope < 0.2
>>> PlotMatrix(ccm_L500, columns, title="CCM Converged Lorenz'96")

CCM_Matrix Lorenz'96 — Cross map matrix of Lorenz’96 5D expanded to 25 variables at library size 500. Slope of cross map vs. library size. CCM matrix (converged cross mapping).#

Embedding Dimension#

EDM methods are predicated on an E (or higher) dimensional embedding from which predictions are made and variables characterized. Robust estimation of embedding dimension remains an open problem (Tan 2023) with EDM taking a data-driven and practical approach: optimization of system predictability as a function of embedding dimension. Here we estimate embedding dimension of the Everglades flow data with exclusionRadius=3 to avoid serial correlation at time steps less than 3 days.

>>> from sciedm import EmbedDimension
>>> df = read_csv("../sciedm/data/S12CD-S333-SumFlow_1980-2005.csv")
>>> edim = EmbedDimension(columns='SumFlow', target='SumFlow', exclusionRadius=3)
>>> EDim = edim.fit_transform(df)

>>> from sciedm.aux_func import PlotEmbedDimension
>>> title = f"{edim.columns} Tp={edim.Tp} exclusionRadius={edim.exclusionRadius}"
>>> ax = PlotEmbedDimension(edim.E_rho_, title=title)

Embedding dimension of Everglades flow — Simplex prediction correlation of Everglades flow as a function of embedding dimension.#

The result suggests an embedding dimension of E=4 is a reasonable choice for simplex and s-map application.

Predict Nonlinear#

Class PredictNonlinear evaluates SMap predictions across a range of theta values.

>>> from pandas import read_csv
>>> from sciedm import PredictNonlinear
>>> df = read_csv("../sciedm/data/S12CD-S333-SumFlow_1980-2005.csv")
>>> pnl = PredictNonlinear(columns='SumFlow', target='SumFlow', E=4, Tp=3)
>>> theta_rho = pnl.fit_transform(df)

>>> from sciedm.aux_func import PlotPredictNonlinear
>>> title=f"{pnl.columns} : {pnl.target}  E={pnl.E}  Tp={pnl.Tp}"
>>> PlotPredictNonlinear(pnl.theta_rho_, title=title)

S-map theta dependence of Everglades flow — SMap prediction correlation of Everglades flow as a function of SMap localization parameter `theta`.#

Here we find evidence of nonlinear state dependence with peak predictability at theta = 3 at E = 4 Tp = 3.