Projection on the softest eigenmode of \(A\) (the one that is unstable when \(\delta<0\)) and on the stiffest one. Over a short window the soft mode is indistinguishable from a random walk.
Many independent systems and trajectories at the current settings. Press Run ensemble.
Mass to the right of zero is a verdict of stability (teal); the dashed rule marks the true \(\delta\).
Runs where \(\hat\Lambda\) has a negative real eigenvalue have no real logarithm and are dropped.
Vary the number of samples \(T\) at fixed \(N\), \(\Delta t\) and \(\delta\). Press Sweep horizon.
Mean \(\lVert A-\hat A\rVert_F\) over up to 40 realisations per point, log–log. The dotted guide has slope \(-1/2\), anchored at the longest horizon.
Share of realisations whose estimate looks stable, over a grid of distances \(\delta\) and observation horizons \(T\Delta t\), at the current \(N\), \(\Delta t\) and coupling ensemble. Left of the vertical rule the true system is unstable, so teal there means the estimator is fooled. Press Phase diagram.
Dashed curves: \(|\delta|\,T\Delta t=1.78\), the scalar unit-root bias scale. Cells needing more than 40 000 samples are skipped.
Fit a multivariate time series: one column per variable, one row per time step, comma, semicolon, tab or space separated, with an optional header row. The data are demeaned, then fitted by the same VAR(1) maximum likelihood and, for comparison, by the Lyapunov inverse method.
Parametric bootstrap with the fitted noise covariance and the same length and \(\Delta t\). Null: the fitted model with its softest mode moved exactly onto the stability boundary. Fitted: the fitted model itself. The rule marks your estimate.
The model
We use the Majumdar–Bialek parametrisation of a random linear system,
$$dX=-AX\,dt+dW,\qquad A=I-\alpha J,$$with \(W\) a standard Brownian motion (noise covariance \(I\,dt\)). In the symmetric case \(J\) is a Gaussian matrix with off-diagonal variance \(1/(2N)\) and diagonal variance \(1/N\), so its spectrum fills a semicircle of radius \(\sqrt2\) and the transition sits at \(\alpha_c=1/\sqrt2\). Rather than fixing \(\alpha\), each sample of \(J\) gets the \(\alpha\) that puts the smallest real part of the spectrum of \(A\) exactly at a chosen \(\delta\):
$$\alpha=\frac{1-\delta}{\max\operatorname{Re}\lambda(J)}\quad\Longrightarrow\quad\min\operatorname{Re}\lambda(A)=\delta .$$Negative \(\delta\) means one mode grows as \(e^{|\delta|t}\): the true system is unstable. The elliptic options draw pairs \((J_{ij},J_{ji})\) with correlation τ (same off-diagonal variance, diagonal variance \((1+\tau)/(2N)\)); τ = 1 recovers the symmetric model, τ = 0 is the real Ginibre ensemble. The same tuning of \(\alpha\) applies, now to the rightmost eigenvalue of \(J\), which may be complex.
Exact simulation
Trajectories start at rest and are sampled every \(\Delta t\) as a VAR(1),
$$X_{n+1}=e^{-A\Delta t}X_n+\varepsilon_n,\qquad \varepsilon_n\sim\mathcal N\!\Big(0,\ \textstyle\int_0^{\Delta t}e^{-As}e^{-A^\top s}\,ds\Big).$$For symmetric \(A\) this is done in its eigenbasis, mode by mode, with variance \((1-e^{-2a_k\Delta t})/(2a_k)\). For non-symmetric \(A\) the propagator and the noise covariance come from one matrix exponential of a \(2N\times2N\) block matrix (Van Loan's method).
The estimator
The exact maximum-likelihood estimate from sufficient statistics is
$$\hat\Lambda=\Big(\sum_n X_{n+1}X_n^\top\Big)\Big(\sum_n X_nX_n^\top\Big)^{-1},\qquad \hat A=-\frac{\log\hat\Lambda}{\Delta t}.$$The estimated system looks stable when every eigenvalue of \(\hat\Lambda\) lies inside the unit circle, which is the same as every eigenvalue of \(\hat A\) having positive real part. Rates are read from the eigenvalues \(\mu\) of \(\hat\Lambda\) as \(-\log\mu/\Delta t\).
Why it happens
Least squares on an autoregression is biased toward stationarity. For a scalar unit root started at rest, the mean bias of the AR coefficient is about \(-1.78/T\), so the implied rate is pushed up by roughly \(1.78/(T\Delta t)\), and the Dickey–Fuller distribution puts about two-thirds of its mass on the stable side even exactly at criticality. A slightly unstable mode keeps looking stable until the observation horizon \(T\Delta t\) is a few times \(1/|\delta|\). What matters is elapsed physical time, not the number of samples: shrinking \(\Delta t\) at fixed \(T\) makes things worse, not better.
In \(N\) dimensions the finite-\(q\) noise in \(\hat\Lambda\) adds a second, random-matrix effect: the whole estimated spectrum spreads out, which can push the soft eigenvalue either way, and inflates the Frobenius error \(\lVert A-\hat A\rVert_F\), which grows roughly like \(\sqrt{N/(q\Delta t)}\) because \(N^2\) entries are fitted from \(T\Delta t=qN\Delta t\) of physical time. Near the transition the Fisher information along the soft direction scales with the variance accumulated by that mode, so the uncertainty on \(\delta\) is set by the horizon, which is exactly what the sweep shows.
The practical message: a fitted drift matrix with a small positive gap is not evidence that a system is stable. If the gap is comparable to \(1/(T\Delta t)\), the data cannot tell stable from unstable.
The Lyapunov (equilibrium) inverse method
The OUinference code fits coupled heterogeneous OU processes \(\dot x_i=-\sum_jJ_{ij}x_j+\sqrt{T_i}\,\eta_i\), \(\langle\eta_i(t)\eta_j(t')\rangle=2\delta_{ij}\delta(t-t')\), with symmetric \(J\), from the equilibrium covariance \(C\) (Barucca, Phys. Rev. E 90, 062129, 2014). The covariance satisfies the Lyapunov equation \(\{C,J\}=2\tilde T\) with \(\tilde T_{ij}=T_i\delta_{ij}\), which inverts in the eigenbasis \(\{c_a,u^a\}\) of \(C\) as
$$J=\sum_{a,b}\frac{2\,(u^{a\top}\tilde T\,u^b)}{c_a+c_b}\,u^a u^{b\top},\qquad T_i=\frac{\langle(\Delta x_i)^2\rangle}{2\,\Delta t}.$$This estimator assumes the series is stationary, and that assumption decides the verdict: if \(Jv=\lambda v\), then \(v^\top(CJ+JC)v=2\lambda\,v^\top Cv=2\,v^\top\tilde Tv>0\), so with a positive-definite \(C\) every eigenvalue of \(J\) is positive. The equilibrium method can only ever return a stable system, which is why this page shows it next to the maximum-likelihood estimate.
Checking your own data
The stability check in Your data is a simple parametric bootstrap in the spirit of a Dickey–Fuller test. It fits \(\hat\Lambda\) and the residual covariance, builds a null model by moving the softest eigenvalue of \(\hat\Lambda\) exactly onto the unit circle (a rank-one change that leaves the other modes untouched when that eigenvalue is real; otherwise the whole matrix is rescaled to spectral radius one), simulates paths of the same length from the null and from the fitted model, and refits each one. The reported share is the fraction of null fits whose smallest rate is at least as large as yours: a small share means your positive gap is unlikely to be produced by the finite-horizon bias alone. It assumes a linear Gaussian VAR(1) with constant coefficients, so treat it as a sanity check rather than a formal test.
Notes
- Errors for symmetric systems are computed in the eigenbasis, where the Frobenius norm is unchanged.
- The matrix logarithm uses inverse scaling and squaring. When \(\hat\Lambda\) has a negative real eigenvalue no real logarithm exists and that run's error is dropped.
- Code: github.com/paolobarucca/OUinference (MATLAB, Lyapunov inverse method for coupled OU processes).