Recovering a system’s dynamics from its data

The starting point is a general multivariate Langevin equation for a vector of ℓ observables Xₜ = {Xₜⁱ}: dXₜⁱ = μⁱ(Xₜ)dt + σᵢⱼ(Xₜ)dWₜʲ, where μⁱ is the drift for observable i, σᵢⱼ is the diffusion matrix coupling observables i and j, and Wₜ is a vector of uncorrelated Wiener processes. This is a Markovian model: it carries no memory kernel, so the future depends only on the present state, not on the path taken to reach it.

The drift and diffusion terms are tied to the Kramers–Moyal (KM) coefficients of the process, D⁽¹⁾ᵢ(X) = limτ→0 ⟨Xᵢₜ₊τ − Xᵢₜ⟩/τ and D⁽²⁾ᵢⱼ(X) = limτ→0 ⟨(Xᵢₜ₊τ − Xᵢₜ)(Xⱼₜ₊τ − Xⱼₜ)⟩/2τ, which coincide with μⁱ and ½σᵢₖσⱼₖ respectively. Estimating these two coefficients from data is therefore equivalent to estimating the equation itself, without ever writing down what that equation should look like beforehand.

The paper truncates the Kramers–Moyal expansion at second order, which is the assumption of Gaussian noise increments. Pawula’s theorem says that if any coefficient beyond second order is non-zero, all higher ones must be too, or the resulting probability density can go negative, a mathematically valid but physically meaningless outcome. Truncating at second order is therefore not an approximation of convenience; it is the only truncation that keeps the reconstructed dynamics well-posed, and the paper checks its validity for each system rather than assuming it.

Computing the KM coefficients from a finite time series is the actual technical problem, and the paper solves it with kernel density estimation (KDE) rather than the more common histogram-binning approach, because binning degrades badly once the state space is multivariate.

For each observable, the empirical dataset of positions and short-time increments is treated as a cloud of samples, and a Gaussian kernel is placed on every sample and summed to approximate the underlying probability density, using Scott’s rule to set the kernel bandwidth from the sample size and the dimensionality of the estimation problem. The drift D⁽¹⁾ requires only a two-dimensional density (a value and its increment); the diffusion D⁽²⁾ requires a three-dimensional one, since it also depends on the correlated observable. The paper notes this stays tractable precisely because each drift term depends only on its own observable and each diffusion term on at most a pair, so the curse of dimensionality never touches the full ℓ-dimensional state at once.

Higher-order coefficients, third and fourth, are also estimated for both real-world systems as a check on the Gaussian assumption itself, rather than taken on faith.

The benchmark is a Brownian particle in a double-well potential, U(X) = ¼X⁴ − ½X², simulated for 2 × 10⁶ steps. The reconstructed drift comes out as a cubic function with roots at the three known equilibria, X = −1, 0, 1, and the reconstructed diffusion matches the true value across 96% of the domain, degrading only at the extremes where the simulated trajectory rarely visits, exactly where a data-driven method should be expected to struggle.

On the Spanish electricity day-ahead market (2004–2020, a 24-dimensional system, one dimension per hour), the reconstructed drift reveals a distinct equilibrium price for every hour rather than the single mean value a univariate model would assume: valley hours (03:00–06:00) settle between €30 and €35/MWh, mid-day hours around €46–48/MWh, and peak hours (20:00–22:00) as high as €55/MWh. The diffusion matrix, which a univariate model has no way to produce at all, shows the largest cross-hour coupling between 08:00–13:00 and 17:00–19:00, at €56–74/MWh² per day, meaning price movements in those hours tend to move together.

On daily EURUSD and GBPUSD exchange rates (December 2003 to March 2024), the reconstructed potentials show multiple minima rather than one, indicating metastability: the rate can settle in a region and stay there for years, as EURUSD did around 1.45 from 2007–2012, before a large enough fluctuation carries it over a saddle point to a different minimum. The effect is real but weaker than in the bistable benchmark: a soft metastability rather than two equally likely, sharply separated states.

For the bistable-potential benchmark, where ground truth is known exactly, the reconstructed drift and diffusion visually and numerically match the analytical forms used to generate the data, which is the paper’s check that the KDE machinery itself is sound before trusting it on real data where no ground truth exists.

For the two real-world systems, the check instead falls on the truncation assumption. Third- and fourth-order KM coefficients are estimated and found to be several orders of magnitude smaller than the first and second, for both the electricity and currency-exchange data, confirming after the fact that the second-order, Gaussian-noise truncation was the right modelling choice rather than merely a convenient one. One caveat is noted directly: the electricity market’s occasional large but infrequent price spikes are a form of non-Gaussian behaviour that a Gaussian-truncated model will not fully capture, and the paper flags large-deviation theory as the natural extension for that regime rather than treating it as already solved.

The electricity day-ahead market is not an incidental test case. It is the market our forecasting products operate in, and a method that recovers metastable regions without being told what the market is makes fewer assumptions than a price-equation model requiring specific domain knowledge.

This paper supplies the stationary picture; the companion work on forecasting couples the same Langevin description to a neural ODE to handle the non-stationary behaviour it cannot capture alone.