Mean Reversion and Heavy Tails: Characterizing Time Series Data Using Ornstein-Uhlenbeck Processes and Machine Learning

1. Introduction

Mean reversion and heavy tails are two recurring features of real-world time series. The Ornstein–Uhlenbeck (OU) process captures mean reversion through a linear drift toward equilibrium [1]. In its Gaussian form, increments are light-tailed, and the second-order structure fully characterizes variability. Many applications, however, show departures from Gaussianity: returns in finance, intermittency in environmental series, and burst-like activity in solar and geophysical records. A tractable extension replaces Brownian motion by a symmetric $\alpha$-stable Lévy driver with stability index $\alpha \in (0,2]$, which preserves mean reversion while allowing heavy tails and occasional large jumps [2,3]. When $\alpha <2$, the variance is infinite, so methods relying on low-order moments can become unstable or biased.

Jointly estimating a mean-reversion rate $\theta$ and a tail index $\alpha$ from finite, noisy data is difficult. Likelihood-based inference for stable-driven OU models requires either long samples or numerical inversion of characteristic functions and careful treatment of normalization; performance deteriorates as $\alpha$ moves away from the Gaussian boundary [4]. Moment and cumulant approaches are not well-defined when moments diverge. Tail-quantile and peak-over-threshold methods focus on extremes but ignore the temporal structure that informs the mean reversion rate $\theta$. These gaps motivate an approach that uses short-window information beyond low-order moments that scales to diverse types of data/signals.

Why is this important? Many operational decisions depend on distinguishing transient noise from persistent regime changes. In finance, separating Gaussian from heavy-tailed periods and identifying changes in reversion speed affects risk limits, hedging, and stress testing [5,6,7,8,9,10]. [8] provide an early and systematic analysis of mean reversion in financial time series based on variance-ratio methods. Their results demonstrate that deviations from a random-walk benchmark can be detected in short windows and that mean-reverting dynamics appear intermittently across asset classes. This aligns with our objective of classifying local windows by reversion strength rather than assuming stable long-run parameters. Relatedly, [11] show a structural increase in complexity and a persistent decline in out-of-sample predictability for major equity indices. They show that increases in approximate and sample entropy correspond to weaker forecasting performance of a wide range of ML models. This supports our motivation to jointly screen windows for Gaussian vs. non-Gaussian features and for strong vs. weak reversion, since both characteristics are linked to shifts in local predictability and model stability. In environmental monitoring, heavy-tailed windows and weak short-horizon reversion indicate intermittent events and regime persistence, which influence planning and early warning [12,13,14]. In scientific data analysis, reliable identification of heavy tails and mean reversion supports model selection and hypothesis testing when moments are undefined or uninformative [2,15]. A window-based classifier that does not rely on strong distributional assumptions enables rolling-window-based analysis, cross-domain comparability, and robust summaries in settings where classical estimators fail, are expensive to deploy, or require long observation horizons, especially for live data stream analysis [16,17,18].

We address this by developing a supervised learning approach that maps window-level features to discrete classes for $(\theta ,\alpha )$. Each sample is a fixed-length segment of a dimensionless OU process; a vector of numeric descriptors summarizes ordinal patterns and short-horizon dynamics without assuming finite variance (e.g., permutation-based features [18]). We use CatBoost, a gradient-boosted decision-tree method designed for tabular inputs [19,20]. CatBoost’s ordered boosting and oblivious-tree structure provide stable and fast training and work well with categorical, as well as in our case, purely numeric features.

The synthetic-data, i.e. the employed Ornstein Uhlenbeck processes establish the ground truth and control the task’s granularity. We generate dimensionless OU series on a grid of $\alpha \in \{0.05,0.5,1.0,1.5,2.0\}$ and $\theta \in \{{10}^{-6},2,4,8,16,32\}$ using exact discrete transitions for both Gaussian and stable drivers. Base series have fixed length and are segmented into overlapping windows of sizes $w\in \{50,100,250\}$ with shift $w/2$. For every w, we compute the same feature set, train two independent CatBoost classifiers (one for $\alpha$, one for $\theta$), and evaluate on independent-seed test sets that share the parameter grids but not the noise realizations. We report per-class precision, recall, and $F_1$ score, overall accuracy, and weighted $F_1$ score, along with absolute and row-normalized confusion matrices. This approach quantifies (i) separability as a function of window length and (ii) where misclassifications occur in parameter space. To reflect a frequent screening use case, we also evaluate the same trained models on binary classifications as well, i.e., Gaussian vs. non-Gaussian ($\alpha$) and mean reversion v.s. no mean reversion ($\theta$)—without retraining. Finally we test, evaluate and interpret our approach on several real-life time series. Figure 1 depicts the main idea.

A central objective of this work is to demonstrate that short-window statistical diagnostics can recover meaningful regime information across real-world systems. We therefore analyze heterogeneous datasets—daily financial returns, daily sunspot numbers, and environmental measurements—not because they share a physical mechanism, but because they offer complementary testbeds with well-documented regime structure. First, these domains exhibit transitions between markedly different dynamical states (e.g., volatility shifts in financial markets, solar-cycle phases in sunspot activity), making them suitable for evaluating whether short-window classifiers can detect changes that traditionally require long-term aggregation [6,21]. Second, each dataset provides a distinct stress case for robustness: financial returns contain heavy tails and clustered volatility, sunspot numbers have quasi-periodic but irregular cycles, and NASPAW data include meteorological noise and intermittent extremes. Demonstrating consistent performance across such systems highlights that the proposed framework functions as a general-purpose analysis tool.

More broadly, the models introduced here should be understood as methodological tools rather than domain-specific predictors. Many classical approaches to estimating distributional properties or mean-reversion rates rely on long samples and strong parametric assumptions; this limits their usefulness when regime changes occur on short horizons or when assumptions (such as finite variance) do not hold. The presented machine-learning approach, by contrast, enables data-driven detection of local dynamical structure from relatively small windows, offering a scalable approach for identifying short-term regime changes in settings where traditional estimators fail or require prohibitively long data records.

This paper makes four contributions. First, it provides an exact, dimensionless simulation of stable-driven OU paths across a wide $(\theta ,\alpha )$ grid for supervised learning. Second, it shows that window-level features paired with boosted tree models recover both tail index and reversion rate with useful fidelity, with errors concentrated between neighboring classes and accuracy improving with window length. Third, it documents a compact binary screening mode that reuses the same models to flag Gaussian vs. non-Gaussian regimes and the presence vs. absence of short-horizon mean reversion. Fourth, it applies the method to financial, environmental, and solar time series to illustrate how regime mixes shift across assets, variables, and time intervals.

The remainder of the paper is organized as follows. Section 3 introduces the Ornstein–Uhlenbeck model with $\alpha$-stable increments, the exact discretization used for generating synthetic benchmarks, and the window-level feature set. Section 4 specifies the supervised learning setup, including the CatBoost classifiers, the windowing procedure, data splits, and the training protocol. Section 5 outlines the validation framework and reports the full multi-class results on synthetic data together with the corresponding confusion matrices; the complementary binary evaluations are provided in Appendix A. Section 6, Section 7 and Section 8 present the empirical applications to financial returns, daily sunspot numbers, and NASA POWER climate variables for Austria, illustrating how the method identifies short-window regime changes across heterogeneous real-world signals. Section 9 summarizes the main findings and discusses limitations and practical use cases. Section 10 conlcudes this article.

2. Related Work

A growing line of research connects time–series complexity with the achievable accuracy of machine and deep learning predictors. Survey evidence shows that entropy- and fractal-based descriptors (e.g., Hurst exponent, fractal dimension, Lyapunov spectrum) help quantify persistence, intermittency, and predictability, and that models which account for these properties tend to perform better on difficult, noisy signals [22]. Building on this, several studies operationalize the link by filtering or shaping model outputs according to signal complexity: ensemble LSTMs combined with complexity filters reduce error and limit overfitting without expensive hyperparameter tuning [23]; reconstructed phase spaces and phase-space–consistent filtering stabilize autoregressive neural forecasts across datasets [24]; and stochastic interpolation that preserves attractor geometry (via fractional Brownian bridges) improves data density while maintaining dynamical structure [25]. Complementary work shows that complexity-aware interpolation—especially when guided by fractal metrics—can lift LSTM performance on short or coarse data and better handle long-memory behavior [26].

Our contribution follows this arc but targets the interpretable regime classification that is useful beyond prediction alone. Instead of using complexity solely as auxiliary features, we estimate two window-local, physically meaningful quantities: tail heaviness ($\alpha$) and mean-reversion speed ($\theta$). Trained on OU processes with $\alpha$-stable increments, the resulting classifiers act as compact “regime sensors” that travel across domains. We validate them on synthetic data and then apply them to three heterogeneous real datasets (finance, sunspots, and environmental series), showing that the same complexity-informed framing supports both predictive tasks (where prior work emphasized LSTM ensembles and interpolation [23,24,25,26,27]) and diagnostic tasks (where our $\{\alpha ,\theta \}$ labels summarize distributional shape and recovery dynamics). In short, this extends the complexity–predictability link surveyed by Raubitzek and Neubauer [22] into a regime-mapping tool that is practical for out-of-sample monitoring in real systems.

3. Ornstein–Uhlenbeck Processes

The Ornstein–Uhlenbeck (OU) process was introduced by Uhlenbeck and Ornstein to model velocity relaxation in a gas [1]. It is now a standard mean-reverting stochastic differential equation (SDE) with broad use in physics, finance, and engineering [28,29,30]. In its basic form the process $\{X_t:t\ge 0\}$ satisfies

$$d{X}_t=-\theta (X_t-\mu )dt+\sigma d{W}_t,$$

(1)

where t denotes time, $\theta >0$ is the mean-reversion rate towards the level $\mu \in \mathbb{R}$, $\sigma >0$ is the diffusion (noise) scale, and ${\left\{W_t\right\}}_{t\ge 0}$ is a standard Brownian motion. We write $\mathcal{N}(m,v)$ for the normal distribution with mean m and variance v.

Equation (1) admits the explicit solution obtained by variation of constants,

$$X_t=\mu +(X_0-\mu )e^{-\theta t}+\sigma {\int }_0^t{e}^{-\theta (t-s)}d{W}_s,$$

(2)

where $X_0$ is the initial value at $t=0$. Conditioned on $X_0=x_0$, the process is Gaussian with

$$\mathbb{E}[X_t\mid X_0=x_0]=\mu +(x_0-\mu )e^{-\theta t},Var(X_t\mid X_0=x_0)={\displaystyle \frac{{\sigma }^2}{2\theta }}\left(1-e^{-2\theta t}\right).$$

(3)

As $t\to \infty$, the conditional mean converges to $\mu$ and the conditional variance converges to ${\sigma }^2/\left(2\theta \right)$, hence

$$X_t\Rightarrow \mathcal{N}\left(\mu ,{\displaystyle \frac{{\sigma }^2}{2\theta }}\right).$$

(4)

If the process is initialized from the stationary distribution (or after transients are negligible), the stationary autocovariance is

$$Cov(X_t,X_{t+\tau })={\displaystyle \frac{{\sigma }^2}{2\theta }}{e}^{-\theta \tau },\tau \ge 0,$$

(5)

so the characteristic relaxation time is $1/\theta$. This relaxation time quantifies the temporal scale over which correlations decay and the process loses memory of past states. Specifically, after a time interval of order $1/\theta$, the autocovariance has decreased by a factor of $e^{-1}$, indicating that deviations from the mean level $\mu$ are largely dissipated. Larger values of $\theta$ therefore correspond to faster mean reversion and shorter memory, while smaller values of $\theta$ imply slower relaxation and more persistent temporal dependence.

Over a fixed sampling interval $\Delta t>0$, the transition density is Gaussian. Sampling at times $t_n=n\Delta t$ and writing $X_n:=X_{t_n}$ gives the exact recursion

$$X_{n+1}=\mu +(X_n-\mu )e^{-\theta \Delta t}+\sigma \sqrt{{\displaystyle \frac{1-e^{-2\theta \Delta t}}{2\theta }}}{Z}_n,Z_n\sim \mathcal{N}(0,1)\mathrm{i}.\mathrm{i}.\mathrm{d}.$$

(6)

Equivalently, the conditional distribution is

$$X_{n+1}\mid X_n\sim \mathcal{N}\left(\mu +(X_n-\mu )e^{-\theta \Delta t},{\displaystyle \frac{{\sigma }^2}{2\theta }}\left(1-e^{-2\theta \Delta t}\right)\right).$$

(7)

This is the closed-form (non-Euler–Maruyama) update; therefore no additional $\sqrt{\Delta t}$ factor multiplies $Z_n$.

For parameter inference it is convenient to remove location and scale while preserving mean reversion. Define the dimensionless, centered, and scaled process

$$Y_t=\sqrt{{\displaystyle \frac{2\theta }{{\sigma }^2}}}(X_t-\mu ).$$

(8)

Since $Y_t$, Equation (8), is affine in $X_t$, Itô’s formula reduces to differentiation of the map $g\left(x\right)=\sqrt{2\theta /{\sigma }^2}(x-\mu )$ with $g^{\prime }\left(x\right)=\sqrt{2\theta /{\sigma }^2}$ and $g^{\prime \prime }\left(x\right)=0$:

$$\begin{array}{cc}\hfill d{Y}_t& =g^{\prime }\left(X_t\right)d{X}_t+\frac{1}{2}{g}^{\prime \prime }\left(X_t\right){\sigma }^{2}dt\hfill \\ \hfill & =\sqrt{{\displaystyle \frac{2\theta }{{\sigma }^2}}}\left(-\theta (X_t-\mu )dt+\sigma d{W}_t\right)\hfill \\ \hfill & =-\theta \sqrt{{\displaystyle \frac{2\theta }{{\sigma }^2}}}(X_t-\mu )dt+\sqrt{2\theta }d{W}_t\hfill \\ \hfill & =-\theta Y_{t}dt+\sqrt{2\theta }d{W}_t.\hfill \end{array}$$

(9)

Hence, in the stationary regime, $Y_t$ has mean 0 and variance 1. The exact discrete transition, Equation (9), over $\Delta t$ is

$$Y_{n+1}=e^{-\theta \Delta t}{Y}_n+\sqrt{1-e^{-2\theta \Delta t}}{Z}_n,Z_n\sim \mathcal{N}(0,1)\mathrm{i}.\mathrm{i}.\mathrm{d}.$$

(10)

Equation (9) keeps $\theta$ as the single free parameter in the Gaussian OU case.

3.1. Replacing the Wiener Driver by Symmetric $\alpha$-Stable Noise

Empirical series in several domains exhibit heavy tails and jumps that occur more often than Gaussian models predict. A natural extension replaces $d{W}_t$ in (9) by the increment of a symmetric ($\beta =0$) $\alpha$-stable Lévy process $L^{\left(\alpha \right)}$.

A random variable X is stable if sums of i.i.d. copies are equal in law to a location–scale transform of X. We write $X\sim S(\alpha ,\beta ,\gamma ,\delta )$ with stability $\alpha \in (0,2]$, skewness $\beta \in [-1,1]$, scale $\gamma >0$, and location $\delta \in \mathbb{R}$. For $\alpha \ne 1$, the characteristic function is

$${\phi }_X\left(t\right)=exp\left\{i\delta t-{\gamma }^{\alpha }{\left|t\right|}^{\alpha }\left[1-i\beta sgn\left(t\right)tan(\pi \alpha /2)\right]\right\}.$$

(11)

For $\alpha =1$ the characteristic function has a different (logarithmic) form; in this work we use the symmetric case $\beta =0$, for which the $\alpha =1$ characteristic function reduces to ${\phi }_X\left(t\right)=exp\{i\delta t-\gamma |t\left|\right\}$. Special cases: $\alpha =2$ gives $\mathcal{N}(\delta ,2{\gamma }^2)$; with $\alpha =1$ and $\beta =0$ one obtains the Cauchy law. Moments satisfy ${\mathbb{E}\left|X\right|}^p<\infty$ iff $p<\alpha$.

Let $\{L_t^{\left(\alpha \right)}:t\ge 0\}$ be a symmetric $\alpha$-stable Lévy motion with independent, stationary increments and

$$L_0^{\left(\alpha \right)}=0,L_t^{\left(\alpha \right)}-L_s^{\left(\alpha \right)}\sim S\left(\alpha ,0,{(t-s)}^{1/\alpha },0\right),0\le s<t,$$

i.e., the increment scale grows like ${(t-s)}^{1/\alpha }$ (unit scale at $t=1$).

We consider the Lévy-driven OU SDE

$$d{Y}_t=-\theta Y_{t}dt+{\sigma }_{\alpha }d{L}_t^{\left(\alpha \right)},\theta >0,0<\alpha \le 2.$$

(12)

The solution admits the moving-average representation

$$Y_t=e^{-\theta t}{Y}_0+{\sigma }_{\alpha }{\int }_0^t{e}^{-\theta (t-s)}d{L}_s^{\left(\alpha \right)}.$$

(13)

In the stationary regime (or for $t\to \infty$), the marginal distribution of $Y_t$ is again $\alpha$-stable [31]. To determine its scale, use the stability property of stochastic integrals with respect to $L^{\left(\alpha \right)}$: for deterministic $f\in L^{\alpha }\left({\mathbb{R}}_{+}\right)$,

$${\int }_0^{\infty }f\left(s\right)d{L}_s^{\left(\alpha \right)}\sim S\left(\alpha ,0,{\left({\int }_0^{\infty }{\left|f\left(s\right)\right|}^{\alpha }ds\right)}^{1/\alpha },0\right).$$

With $f\left(s\right)={\sigma }_{\alpha }{e}^{-\theta s}$, we obtain the stationary scale parameter

$${\gamma }_{\infty }={\sigma }_{\alpha }{\left({\int }_0^{\infty }{e}^{-\alpha \theta s}ds\right)}^{1/\alpha }={\sigma }_{\alpha }{\left({\displaystyle \frac{1}{\alpha \theta }}\right)}^{1/\alpha }.$$

(14)

To keep the dimensionless convention (unit stationary scale for all $\alpha$), we choose

$${\sigma }_{\alpha }={\left(\alpha \theta \right)}^{1/\alpha },$$

(15)

which yields ${\gamma }_{\infty }=1$ and hence the stationary marginal $S(\alpha ,0,1,0)$ for every $\alpha$.

Let $\Delta t>0$ and $t_n=n\Delta t$. From (13),

$$Y_{t_{n+1}}=e^{-\theta \Delta t}{Y}_{t_n}+{\sigma }_{\alpha }{\int }_{t_n}^{t_{n+1}}{e}^{-\theta (t_{n+1}-s)}d{L}_s^{\left(\alpha \right)}.$$

By stationarity and independence of Lévy increments, the integral term is independent of $Y_{t_n}$ and has the same law as ${\sigma }_{\alpha }{\int }_0^{\Delta t}{e}^{-\theta (\Delta t-u)}d{L}_u^{\left(\alpha \right)}$. Using the scaling rule for stable integrals,

$$\begin{array}{cc}\hfill {\sigma }_{\alpha }{\int }_0^{\Delta t}{e}^{-\theta (\Delta t-u)}d{L}_u^{\left(\alpha \right)}& \sim S\left(\alpha ,0,{\sigma }_{\alpha }{\left({\int }_0^{\Delta t}{e}^{-\alpha \theta (\Delta t-u)}du\right)}^{1/\alpha },0\right)\hfill \\ \hfill & =S\left(\alpha ,0,{\sigma }_{\alpha }{\left({\int }_0^{\Delta t}{e}^{-\alpha \theta v}dv\right)}^{1/\alpha },0\right),\hfill \end{array}$$

(16)

where $v=\Delta t-u$. Therefore the exact transition can be written as

$$Y_{n+1}=e^{-\theta \Delta t}{Y}_n+{\sigma }_{\alpha }{\left({\int }_0^{\Delta t}{e}^{-\alpha \theta v}dv\right)}^{1/\alpha }{S}_n,S_n\sim S(\alpha ,0,1,0)\mathrm{i}.\mathrm{i}.\mathrm{d}.$$

(17)

Evaluating the integral gives ${\int }_0^{\Delta t}{e}^{-\alpha \theta v}dv=(1-e^{-\alpha \theta \Delta t})/\left(\alpha \theta \right)$, so with the normalization (15) this simplifies to

$$Y_{n+1}=e^{-\theta \Delta t}{Y}_n+{\left(1-e^{-\alpha \theta \Delta t}\right)}^{1/\alpha }{S}_n,S_n\sim S(\alpha ,0,1,0)\mathrm{i}.\mathrm{i}.\mathrm{d}.$$

(18)

For $\alpha =2$, $S_n$ corresponds to a Gaussian draw and (18) reduces to the Gaussian transition (10).

In the implementation, we draw $S_n$ using scipy.stats.levy_stable, which follows Nolan’s $S^0$ parameterization. At $\alpha =2$, $\beta =0$, and scale$=1$, this sampler returns $\mathcal{N}(0,2)$ rather than $\mathcal{N}(0,1)$. Therefore, when $\alpha =2$ we rescale the draw by $1/\sqrt{2}$ to match $Z\sim \mathcal{N}(0,1)$ in (10). Alternatively, one may bypass the stable sampler at $\alpha =2$ and draw Z directly.

$\theta$

exponential pull toward zero; larger $\theta$ implies faster mean reversion.

$\alpha$

tail index; smaller $\alpha$ yields heavier tails and larger jumps; $\alpha =2$ is Gaussian.

${\sigma }_{\alpha }$

noise scale in (12); with ${\sigma }_{\alpha }={\left(\alpha \theta \right)}^{1/\alpha }$ the stationary marginal scale is 1 for all $\alpha$.

The combination of mean-reverting drift and heavy-tailed shocks provides a flexible baseline for time series where both features are present. See Samorodnitsky and Taqqu [2], Nolan [15], Marquardt [31] for detailed background.

3.2. Simulation Setup and Implementation

We use symmetric $\alpha$-stable noise ($\beta =0$) with $\alpha \in \{0.05,0.5,1.0,1.5,2.0\}$, spanning the range from Gaussian ($\alpha =2$) to extremely heavy-tailed regimes (small $\alpha$).

Stable variates are generated with

$$S_n\sim \mathtt{levy}\_\mathtt{stable}.\mathtt{rvs}(\alpha ,\beta =0,\mathrm{scale}=1,\mathrm{loc}=0),$$

which returns i.i.d. draws from $S(\alpha ,0,1,0)$. Inserting these into (18) yields the exact one-step update used for simulation. At $\alpha =2$, scale$=1$ corresponds to $\mathcal{N}(0,2)$ (i.e., norm with scale$=\sqrt{2}$), so we rescale by $1/\sqrt{2}$ to obtain $Z\sim \mathcal{N}(0,1)$ [32].

To generate time series for supervised learning, we discretize the dimensionless OU model on a uniform grid.

We use the exact transition (10):

$$Y_{n+1}=e^{-\theta \Delta t}{Y}_n+\sqrt{1-e^{-2\theta \Delta t}}{Z}_n,Z_n\sim \mathcal{N}(0,1)\mathrm{i}.\mathrm{i}.\mathrm{d}.$$

(19)

For $\alpha \in (0,2]$ we use the exact transition (18):

$$Y_{n+1}=e^{-\theta \Delta t}{Y}_n+{\left(1-e^{-\alpha \theta \Delta t}\right)}^{1/\alpha }{S}_n,S_n\sim S(\alpha ,0,1,0)\mathrm{i}.\mathrm{i}.\mathrm{d}.$$

(20)

This choice is consistent with the SDE normalization ${\sigma }_{\alpha }={\left(\alpha \theta \right)}^{1/\alpha }$ and keeps the stationary marginal scale equal to one for all $\alpha$.

The resulting dimensionless series provide a homogeneous input space across $(\theta ,\alpha )$ for feature extraction and learning.

These dimensionless Ornstein–Uhlenbeck processes are depicted in Figure 2 and Figure 3 for several parameter combinations $(\theta ,\alpha )$.

Stability index.

We explore the five values

$$\alpha \in \{0.05,0.5,1.0,1.5,2.0\},$$

each with skewness fixed at $\beta =0$ (symmetric law). The choice spans the full range from the Gaussian case ($\alpha =2$) to an extreme heavy-tail regime ($\alpha =0.05$). For reference:

(1)

$\alpha =2.0$: normal distribution; finite mean and variance.

(2)

$\alpha =1.5$: infinite variance but finite mean; moderate tails.

(3)

$\alpha =1.0$: Cauchy-type; mean and variance both infinite.

(4)

$\alpha =0.5$: very heavy tails; large jumps common.

(5)

$\alpha =0.05$: extreme heavy-tail regime with rare but dominating outliers.

Mean-reversion rate.

A separate grid

$$\theta \in \{{\theta }_1,\dots ,{\theta }_m\},$$

with ${\theta }_j>0$, controls how fast the process returns to equilibrium. Concrete values are specified in Section 5.

Noise scale.

We keep the stationary marginal scale equal to one by setting

$${\sigma }_{\alpha }={\left(\alpha \theta \right)}^{1/\alpha },$$

see equation (15). This normalization allows direct comparison of paths across different $\alpha$.

Stable increments are generated with scipy.stats.levy_stable.rvs$(\alpha ,\beta =0,\mathtt{scale}=1,\mathtt{loc}=0)$ [32,33]. The call returns i.i.d. variables $S_n\sim S(\alpha ,0,1,0)$ that feed the discrete update below. When $\alpha =2$, SciPy’s stable law corresponds to norm with scale$=\sqrt{2}$ [33]; therefore we rescale by $1/\sqrt{2}$ when we require $Z\sim \mathcal{N}(0,1)$.

Let $Y_n=Y_{t_n}$ with $t_n=n\Delta t$. The exact transition of equation (12) over one step $\Delta t$ is

$$Y_{n+1}=e^{-\theta \Delta t}{Y}_n+{\left(1-e^{-\alpha \theta \Delta t}\right)}^{1/\alpha }{S}_n,S_n\sim S(\alpha ,0,1,0)\mathrm{i}.\mathrm{i}.\mathrm{d}.$$

(21)

Derivation. Start from the moving-average solution $Y_{t_{n+1}}=e^{-\theta \Delta t}{Y}_{t_n}+{\sigma }_{\alpha }{\int }_{t_n}^{t_{n+1}}{e}^{-\theta (t_{n+1}-s)}d{L}_s^{\left(\alpha \right)}$. By independent increments, the integral term is independent of $Y_{t_n}$ and has the same law as ${\sigma }_{\alpha }{\int }_0^{\Delta t}{e}^{-\theta (\Delta t-u)}d{L}_u^{\left(\alpha \right)}$. For symmetric $\alpha$-stable Lévy motion, the scale of such an integral equals ${\sigma }_{\alpha }{\left({\int }_0^{\Delta t}{e}^{-\alpha \theta v}dv\right)}^{1/\alpha }={\sigma }_{\alpha }{\left({\displaystyle \frac{1-e^{-\alpha \theta \Delta t}}{\alpha \theta }}\right)}^{1/\alpha }$, and with ${\sigma }_{\alpha }={\left(\alpha \theta \right)}^{1/\alpha }$ this becomes ${(1-e^{-\alpha \theta \Delta t})}^{1/\alpha }$.

For $\alpha =2$ we obtain

$$Y_{n+1}=e^{-\theta \Delta t}{Y}_n+\sqrt{1-e^{-2\theta \Delta t}}{Z}_n,Z_n\sim \mathcal{N}(0,1)\mathrm{i}.\mathrm{i}.\mathrm{d}.,$$

which matches the exact Gaussian OU update.

$e^{-\theta \Delta t}$

deterministic decay over one step.

${(1-e^{-\alpha \theta \Delta t})}^{1/\alpha }$

innovation scale; small when $\Delta t\ll 1/\theta$.

$S_n$

independent symmetric $\alpha$-stable driver controlling tail heaviness via $\alpha$.

Equation (21) is used throughout the study to simulate training and test paths for every $(\theta ,\alpha )$ pair.

Algorithm implemented in scipy.stats.levy_stable.

SciPy draws random numbers from a Lévy–stable law with the rvs method of the class scipy.stats.levy_stable [33]. Internally this routine follows the Chambers–Mallows–Stuck (CMS) transformation [34].

Dimensionless scaling is critical for ML comparability because it removes the location and scale degrees of freedom and fixes a common stationary marginal scale across different $(\theta ,\alpha )$ settings. This ensures that learned features and decision boundaries reflect differences in mean reversion and tail behavior rather than trivial changes in units or amplitude. Using exact discrete transitions instead of Euler–Maruyama updates avoids discretization bias that is pronounced for short windows and coarse sampling, and it prevents artificial smoothing of jumps in the $\alpha$-stable case. This is important for the following because our classifier is trained on window-level statistics, where small systematic simulation errors can translate into consistent label noise and thus degrade generalization. This, as we will later see, will still be present and visible by confusions for non-extreme regimes where $\alpha$ and $\theta$ are not at their respective boundary (highest/lowest) values.

4. Machine Learning

We use CatBoost, a gradient-boosting decision tree method, as the primary classifier for all window-level experiments [19,20]. Gradient boosting builds an ensemble sequentially, each tree trained to reduce residual errors of the current model [35]. In contrast to bagging methods such as Random Forests [36], boosting focuses subsequent learners on remaining mistakes, which is effective for tabular feature sets with non-linear interactions.

CatBoost employs ordered boosting to limit target leakage and reduce overfitting, and uses symmetric (oblivious) trees that provide stable splits and fast inference [19,20]. Although CatBoost is often chosen for mixed numeric–categorical inputs, our features are entirely numeric (window-level summary and complexity metrics). The method remains suitable because it handles tabular numeric features with minimal preprocessing and performs reliably with default settings, as shown in [37,38,39].

Each sample is a fixed-length sliding window extracted from a simulated or empirical time series. The feature vector contains only numeric descriptors computed on that window. We train two independent multiclass models on the same feature representation: one predicts the stability index $\alpha$ and one predicts the mean-reversion rate $\theta$. The parameter grids for $\alpha$ and $\theta$ follow the simulation design defined earlier (Section 5). Training and test sets use the same grids but independent noise seeds, so test scores measure generalization to new noise realizations with the same statistics. Base series have length 2000 and are segmented into overlapping windows of length $w\in \{50,100,200,250\}$ with shift $w/2$.

Models use CatBoostClassifier with the library’s multiclass loss and otherwise default hyperparameters. For each window length w, we train a separate pair of models (one for $\alpha$, one for $\theta$).

For each w and each target ($\alpha$ or $\theta$), we construct a stratified split into 80% train and 20% test (stratification by the class label). From the training portion, we carve out an internal validation subset (10–20% of the training pool, stratified) used exclusively for early stopping and model selection.

We then perform a compact Bayesian hyperparameter search around the CatBoost GPU defaults using 3-fold cross-validation on the training set only. The search space centers on depth$\in [5,7]$, iterations$\in [900,1100]$, l2_leaf_reg$\in [2,5]$, and border_count$\in [120,136]$; the objective is cross-validated accuracy. This yields a tuned candidate model. Between the baseline and the tuned candidate, we select the final model using performance on the internal validation data (or equivalently by cross-validated performance within the training set). The held-out test split is used only once, for the final evaluation. The chosen model is saved per window length and target.

Evaluation is performed once, on the held-out test split. We report per-class precision, recall, and $F_1$-score, overall accuracy, and weighted $F_1$-score. We also produce confusion matrices in two forms: absolute counts and row-normalized (to show conditional accuracy by true class). These outputs are shown in Table 1, Table 2 and Table 3 and Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. For the coarse screening use case, we reuse the same trained models without retraining to form binary tasks (Gaussian vs. non-Gaussian for $\alpha$; mean reversion vs. no mean reversion for $\theta$). Reproducibility is ensured by fixed parameter grids defined in the simulation section, documented seeds for training and testing generation, identical windowing and feature extraction across runs, and a consistent train/validation/test protocol per window size and target.

5. Initial Model Validation

We estimate two window-level quantities: the stability index $\alpha$ (tail heaviness) and the mean-reversion rate $\theta$. Because both are local and depend on time scale, we evaluate models on sliding windows of fixed length $w\in \{50,100,200,250\}$. Larger w provides more information per decision; smaller w provides finer time localization. Windows come from the dimensionless OU series generated on the $(\alpha ,\theta )$ grid defined in the simulation section (Section 3.2). Training and test sets use the same parameter grids but independent noise seeds, so test scores reflect generalization to new realizations with the same statistics.

Each base series (length 2000) is split into overlapping windows with shift $w/2$. Every window inherits the pair $(\alpha ,\theta )$ used for its source series. We train two CatBoost classifiers, one for $\alpha$ and one for $\theta$, and evaluate them on a held-out test pool generated with different seeds. For every w, we report per-class precision, recall, and $F_1$-score, plus overall accuracy and weighted $F_1$-score. Classification reports are provided in Table 1, Table 2 and Table 3. To describe error structure, we use row-normalized confusion matrices (Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9).

Across window sizes, the $\alpha$ models reach mid–70% accuracy at $w=50$ and remain close to that level at larger w in the five-class setup. The extremes $\alpha =0.05$ and $\alpha =2.0$ are classified best; intermediate $\alpha$ values show most confusion with their neighbors. The $\theta$ models improve more with longer windows: accuracy rises from about 67% at $w=50$ to about 84% at $w=250$. The strongest reversion class $\theta =32$ has the highest per-class scores. Confusion is mostly between adjacent $\theta$ bins (e.g., 4 vs. 8, 8 vs. 16), indicating ordinal errors, i.e., confusing processes in similar in mean reversion, rather than arbitrary swaps.

The row-normalized confusion matrices show diagonal dominance at all window sizes and a tighter diagonal as w grows. For $\alpha$, longer windows increase the chance that rare, large deviations appear inside the window, which improves separation of very heavy-tailed and Gaussian regimes. For $\theta$, longer windows expose more of the pull-to-mean pattern or the lack of it, which improves separation of reversion rates. Figure 4 and Figure 5 illustrate the shortest-window setting, while Figure 8 and Figure 9 show the strongest separation at $w=250$.

We also run a binary evaluation without retraining, using the same models and windows. For $\alpha$, we group non-Gaussian states into a single “Lévy” class and test “Gaussian vs. Lévy”; for $\theta$, we test “no mean reversion ($\theta ={10}^{-6}$) vs. mean reversion (all other $\theta$).” Full binary reports are provided in Appendix A. In short, the Gaussian vs. Lévy task reaches about 94–96% accuracy even at $w=50$, and the no mean reversion vs. mean reversion task rises from about 89% to about 95% with window size. This indicates that the same trained models can be reused for coarse screening without additional fitting.

Overall, with independent-seed test sets and matched parameter grids, the models generalize consistently across window sizes. Multi-class errors are concentrated in neighboring labels, and the binary appendix confirms that the models provide reliable screening for Gaussian vs. non-Gaussian behavior and for mean reversion vs. no mean reversion.

5.1. Summary

The validation shows three main outcomes. First, the multiclass $\alpha$ classifier separates very heavy-tailed and Gaussian regimes reliably, with most errors confined to neighboring intermediate tail classes. Second, the multiclass $\theta$ classifier benefits strongly from longer windows, indicating that mean-reversion strength requires more temporal context than tail classification in this setup. Third, the same models achieve high accuracy on the two binary tasks, supporting their use as window-based detectors for Gaussian vs. non-Gaussian behavior and mean reversion vs. no mean reversion in analyses.

Table 1. Classification Reports for Sliding Window of 50 Points.

Table 1. Classification Reports for Sliding Window of 50 Points.

Table 2. Classification Reports for Sliding Window of 100 Points.

Table 2. Classification Reports for Sliding Window of 100 Points.

Table 3. Classification Reports for Sliding Window of 250 Points.

Table 3. Classification Reports for Sliding Window of 250 Points.

Figure 4. Row-normalized confusion matrix for $\alpha$ classification with a window size of 50.

Figure 4. Row-normalized confusion matrix for $\alpha$ classification with a window size of 50.

Figure 5. Row-normalized confusion matrix for $\theta$ classification with a window size of 50.

Figure 5. Row-normalized confusion matrix for $\theta$ classification with a window size of 50.

Figure 6. Row-normalized confusion matrix for $\alpha$ classification with a window size of 100.

Figure 6. Row-normalized confusion matrix for $\alpha$ classification with a window size of 100.

Figure 7. Row-normalized confusion matrix for $\theta$ classification with a window size of 100.

Figure 7. Row-normalized confusion matrix for $\theta$ classification with a window size of 100.

Figure 8. Row-normalized confusion matrix for $\alpha$ classification with a window size of 250.

Figure 8. Row-normalized confusion matrix for $\alpha$ classification with a window size of 250.

Figure 9. Row-normalized confusion matrix for $\theta$ classification with a window size of 250.

Figure 9. Row-normalized confusion matrix for $\theta$ classification with a window size of 250.

6. Daily Asset Data and Interval Comparison

Short-horizon distributional properties of financial returns are known to vary across market regimes and across assets. Heavy tails, volatility persistence, and state-dependent dependence structures are well-documented empirical features [6,40]. These characteristics can shift over time in response to changes in market microstructure, liquidity, or macroeconomic conditions [41]. Because our machine-learning models infer tail behavior ($\alpha$) and mean-reversion strength ($\theta$) from short windows, applying them to long daily time series allows us to examine whether the relative prevalence of distributional regimes has changed across decades. This situates the method within a broader empirical literature on evolving market structure and non-Gaussian return behavior.

We analyze daily end-of-day price series for four representative assets: Apple (AAPL), Microsoft (MSFT), the S&P 500 (GSPC), and the Dow Jones Industrial Average (DJI). The history is divided into two non-overlapping intervals:

$$\mathbf{Period}\mathbf{1}:1995--2004,\mathbf{Period}\mathbf{2}:2010--2024.$$

These intervals are selected to provide sufficiently long, contiguous samples for all assets and to contrast an earlier market environment with a more recent one. In practice, we chose them because they offer large and symmetric time spans with virtually no missing data for all four assets, i.e., primarily for reasons of data availability. All results are reported per asset and per interval.

To extract local distributional structure, each series is scanned using a sliding window of length 50 with step size 1. For every window, the trained classifiers assign one category for the stability index $\alpha$ and one for the mean-reversion rate $\theta$. This produces a categorical trace per asset. For each interval, we count how many windows fall into each discrete $\alpha$ or $\theta$ bin. Figure 10 and Figure 11 summarize the annualized distribution of category assignments, and Table 4 and Table 5 provide absolute bin counts for each asset and interval. The binary classification tables are collected in Appendix B.

Table 4 and Table 5 summarize the number of windows classified into each $\alpha$ or $\theta$ category. These counts provide a coarse view of how the empirical distribution of local regimes changes between intervals.

Interpretation. 

Two consistent shifts appear when comparing Period 1 and Period 2. First, both broad indices (DJI, GSPC) show a notable reduction in windows classified as Gaussian ($\alpha =2.0$), with a corresponding increase in $\alpha =1.5$. AAPL and MSFT exhibit only modest changes. This indicates that index-level returns display heavier effective tails in the post-2010 sample, aligning with documented findings that aggregate markets can experience structural shifts in volatility and tail behavior [6,40].

Second, windows classified into the strongest mean-reversion state ($\theta =32$) decline for all assets except MSFT. This shift is consistent with a regime in which large deviations occur more frequently and take longer to dissipate, reducing the prevalence of tight, high-frequency mean reversion. These patterns do not contradict informational efficiency: EMH does not require Gaussianity [42,43]. Rather, they indicate that distributional characteristics of returns can evolve even when markets remain informationally efficient.

7. Daily Sunspot Data and Interval Comparison

Short-term statistical properties of solar activity are known to vary across solar cycles, including changes in dispersion, intermittency, and persistence [21,44]. These variations reflect the underlying solar magnetic dynamo and its modulation across decades. Because our machine-learning models quantify local tail behavior ($\alpha$) and local mean-reversion strength ($\theta$), applying the method to long daily sunspot records allows us to examine whether the prevalence of these regimes shifts systematically across different phases of solar activity. This links the analysis to empirical studies documenting structural differences between strong and weak solar cycles and motivates the use of short-window statistical diagnostics as complementary indicators of cycle-dependent behavior.

The daily international sunspot number series, published by SILSO (Sunspot Index and Long-term Solar Observations), provides a standardized measure of solar magnetic activity corrected for observational biases [45,46]. The time series reflects the emergence, evolution, and decay of sunspot groups and exhibits the well-known $\sim 11$-year cycle. For this study, the record is divided into two non-overlapping intervals:

$$\mathbf{Period}\mathbf{1}:1985--2004,\mathbf{Period}\mathbf{2}:2005--2024.$$

These intervals were selected primarily based on data availability: they are long, contiguous segments with no missing observations and can be split into two approximately symmetric parts that permit a comparison between an earlier phase dominated by relatively strong solar cycles 22–23 and a later phase that includes the deep minimum between cycles 23 and 24, the comparatively weak solar cycle 24, and the early phase of solar cycle 25 [21,44,47]. All results are reported per interval. For clarity, the solar cycles discussed here span the following official SILSO/NOAA intervals: Cycle 22 (1986–1996), Cycle 23 (1996–2008), Cycle 24 (2008–2019), and Cycle 25 (2019–present).

To extract local statistical regimes, the daily series is scanned with a sliding window of length 50 and step size 1. Each window is assigned an $\alpha$ category (gaussian or levy) and a $\theta$ category (mean_rev or no_mean_rev) using the same trained models applied in the financial analysis. This produces two binary categorical sequences that summarize local tail behavior and local persistence. Annual counts of these binary categories are shown in Figure 12 for $\alpha$ and Figure 13 for $\theta$. Both use logarithmic scaling for the category counts and overlay the yearly average sunspot number for reference. Figure 12: The gaussian category peaks near solar maxima, indicating that high-activity phases are dominated by many moderate, approximately symmetric fluctuations. In contrast, levy windows cluster more strongly around solar minima, where long quiet intervals are punctuated by intermittent larger changes, yielding heavier tails. The transition from cycles 22–23 to the weaker cycle 24 and early cycle 25 is visible as a relative increase in levy windows during the later decades, consistent with more irregular and intermittent activity.

Figure 13: Years with high sunspot activity are dominated by mean_rev windows, indicating that local deviations tend to be rapidly corrected during strong cycles. Around extended minima and in the weaker recent cycles, the share of no_mean_rev windows increases, reflecting slower return toward equilibrium and more persistent deviations. The co-variation of these counts with the sunspot cycle suggests that local mean-reversion strength is itself a cycle-dependent property of the solar dynamo.

Interval-level totals are reported in Table 6 and Table 7. These aggregated counts complement the annual dynamics by showing how the relative frequencies of distributional and dynamical regimes change across the two multi-cycle periods.

The earlier interval (1985–2004) is dominated by windows classified as Gaussian and mean reversion. This suggests that local increments tend to be moderate, symmetric, and quickly corrected, consistent with the relatively strong and regular solar cycles characterizing this period. High-activity cycles naturally generate many small fluctuations driven by frequent sunspot emergence, producing distributions closer to Gaussian and dynamics consistent with strong mean reversion.

The later interval (2005–2024) shows a clear increase in levy-classified windows, rising from 2500 to 3575. This shift indicates heavier-tailed short-term changes, consistent with the long quiet stretches and intermittent bursts observed in weaker and more irregular cycles (especially the deep minimum between cycles 23 and 24 and the comparatively weak cycle 24). The $\theta$-classification aligns with this pattern: the number of no_mean_rev windows increases substantially, implying that deviations persist longer and the process exhibits weaker corrective dynamics. Such behavior reflects the extended minima, asymmetric cycle shapes, and lower amplitudes documented in recent solar activity  [21,44].

These findings form a coherent picture. Stronger cycles tend to generate frequent moderate fluctuations, leading to Gaussian-like and mean-reverting local dynamics. Weaker cycles produce intermittent, larger deviations that increase heavy-tailed behavior and decrease local mean reversion. Importantly, observing these differences is not only consistent with known solar-cycle variability but also valuable for future inference. If only short-window statistics ($\alpha$ and $\theta$) are available, the observed combination of heavier tails and weaker mean reversion provides evidence that the system is in a quieter and more irregular cycle, whereas predominantly Gaussian and mean-reverting behavior indicates a more active and regular phase. Since the intervals were intentionally chosen to contrast such regimes, the presence of clear statistical differences supports the usefulness of these local diagnostics for identifying the dynamical state of the solar cycle.

8. NASA POWER Daily Data for Austria

Daily reanalysis-based surface climate data provide a practical testbed for our regime-classification framework. They combine long, homogeneous time series with physically interpretable variables that are relevant for energy planning, agriculture, and impact studies [48,49,50]. Austria is a suitable case because it spans complex topography and continental-to-alpine climate gradients within a relatively small area, so that changes in short-horizon variability can be compared consistently across diverse conditions. By applying the $(\alpha ,\theta )$ classification to multi-decadal NASA POWER records, we assess whether the method recovers known contrasts between strongly forced and weather-modulated variables and whether it detects systematic shifts between two long, symmetric calendar intervals.

The Austrian surface-climate dataset analyzed here is derived from NASA’s POWER (Prediction Of Worldwide Energy Resources) project. POWER is an operational dissemination service hosted at NASA Langley Research Center that delivers solar and meteorological surface variables for energy, agrometeorology, and research applications, with consistent units, harmonized time standards, and long, gap-free coverage [48,50]. We retrieve daily records for a dense network of Austrian locations spanning all federal states—major urban centers (Vienna, Graz, Linz, Salzburg, Innsbruck, Klagenfurt, Bregenz), representative alpine valleys and high-elevation sites, and rural lowlands across Lower/Upper Austria, Styria, Carinthia, Tyrol, Vorarlberg, Salzburg, and Burgenland. Each site contributes a multi-decadal series from 1985-01-01 to 2024-12-31 for the variables used in this study: all-sky, clear-sky, and top-of-atmosphere shortwave irradiance; 10 m mean and daily-maximum wind speed; 2 m relative humidity; bias-corrected daily precipitation total; surface pressure; diurnal temperature range; and 2 m dew/frost-point temperature. The Daily API and its query parameters are documented at power.larc.nasa.gov/docs/services/api/temporal/daily/.

POWER provides point values by spatially sampling and interpolating authoritative gridded sources. For the parameters used here: (i) the meteorological variables (e.g., pressure, humidity, temperatures, winds, precipitation) originate from the NASA GMAO MERRA-2 reanalysis [49]; (ii) the surface and TOA radiative fields are sourced from the CERES (Clouds and the Earth’s Radiant Energy System) family of radiation products (CERES-SYN1deg and FLASHFlux), specialized for near-surface energy applications [51]. At delivery time, POWER harmonizes units—e.g., solar irradiance in MJ m⁻² day⁻¹—and applies standard formatting and metadata. A concise provenance summary for POWER, which lists MERRA-2 for meteorology and CERES-SYN1deg/FLASHFlux for radiation, along with product lineages and links to upstream documentation, is provided at power.larc.nasa.gov/data-sources.

Within POWER, PRECTOTCORR is a daily bias-corrected precipitation total designed for surface applications; radiation variables include ALLSKY_SFC_SW_DWN (all-sky surface shortwave on a horizontal plane), CLRSKY_SFC_SW_DWN (clear-sky surface shortwave), and TOA_SW_DWN (top-of-atmosphere shortwave), all provided as daily totals in MJ m⁻² day⁻¹ Wind fields (WS10M, WS10M_MAX) are in m s⁻¹, humidity in %, pressure in kPa, diurnal temperature range in °C, and dew/frost-point in °C. The POWER Parameter Dictionary documents the full variable set, naming, units, and community-specific conventions, and is the reference for parameter semantics used in this work: power.larc.nasa.gov/docs/services/api/parameters/.

From a data-generation standpoint, MERRA-2 is NASA’s modern global atmospheric reanalysis produced by the Global Modelling and Assimilation Office (GMAO). It assimilates a wide range of satellite and in-situ observations into the GEOS modelling system to yield a spatially and temporally consistent estimate of atmospheric state and fluxes from 1980 to the present [49]; POWER serves MERRA-2-based surface meteorology at a daily cadence suitable for agrometeorological and energy applications [48]. The POWER data-source documentation summarizes this linkage and provides upstream references to GMAO and the CERES radiation suite (SYN1deg for consistent climate records; FLASHFlux for low-latency fluxes). Users should consult that page for details on spatial resolution, updates, and cross-product consistency: power.larc.nasa.gov/data-sources.

To characterize short-scale regime behavior consistently across variables and sites, each daily series is processed with a 50-day sliding window and step size of one day. For every window, we estimate the tail-heaviness $\alpha$ and the mean-reversion $\theta$. Windows are assigned to discrete bins of these two quantities, and absolute counts per category are accumulated over two multi-decadal intervals, 1985–2004 and 2005–2025. These intervals are chosen to be as long and symmetric as possible while maintaining uniform data availability for all locations and variables, so that changes in regime composition can be interpreted as temporal shifts rather than artifacts of coverage. The resulting distributions summarize how short-horizon dynamics change between the two periods, both for individual physical variables and for the full Austrian composite. A spatial overview of the analyzed locations is shown in Figure 14.

Following the spatial overview, we first present absolute counts by $\alpha$ and $\theta$ category for all variables and both periods (Table 8 and Table 9). These tables provide the primary basis for comparing shifts in short-horizon regimes between the time periods 1985 to 2004 and 2005 to 2025. Complementary binary summaries that collapse categories into Gaussian vs. Lévy for $\alpha$ and mean-reverting vs. non-mean-reverting for $\theta$ are reported in the Appendix (Table A6 and Table A7); they corroborate the main count patterns and are referenced in the discussion below.

Figure 15 highlights a marked redistribution in CLRSKY_SFC_SW_DWN around the year 2000: the annual share of Gaussian windows increases while Lévy windows drop; in parallel, mean-reverting windows become more dominant relative to non-mean-reverting ones. This visual change is consistent with the pooled counts in Table 8 and Table 9 and with the binary summaries in the Appendix (Table A6, Table A7), where CLRSKY shows a higher $\alpha =2.0$ mass and a relative tilt toward stronger reversion in the timespan 2005 to 2025 compared to 1985 to 2004. The fact that these differences appear clearly in $(\alpha ,\theta )$ alone is useful for future applications: it implies that, when only short-window statistics are available, changes in the joint $(\alpha ,\theta )$ regime can be used as indicators of underlying shifts in the clear-sky radiative environment.

The aggregated counts across all Austrian locations reveal coherent and physically interpretable patterns. Radiative variables show strong contrasts between deterministic and weather-modulated components. The top-of-atmosphere shortwave flux (TOA_SW_DWN) concentrates almost entirely in the heaviest-tail side of the $\alpha$ grid and in the weakest-reversion bins of $\theta$ in both periods, reflecting the smooth but non-stationary seasonal cycle: windows drift slowly along the annual irradiance envelope with little tendency to revert to a local mean, producing very low $\theta$ values and an apparent departure from Gaussianity. In contrast, the surface-level shortwave components (ALLSKY_SFC_SW_DWN, CLRSKY_SFC_SW_DWN) allocate much larger mass to $\alpha =2.0$ and to high $\theta$, consistent with faster short-term equilibration driven by day-to-day weather variability. Between periods, small but systematic redistributions are visible: for instance, $\alpha =2.0$ counts rise for CLRSKY (395,689 vs. 310,660) and $\theta =32$ decreases slightly for ALLSKY (434,606 vs. 456,012), indicating a modest rebalancing between quickly mean-reverting, cloud-dominated variability and slowly varying seasonal components. These shifts are mirrored in the Appendix’s binary summaries, where CLRSKY shows more Gaussian and mean-reverting windows after 2005 (Table A6, Table A7). Because these shifts are captured directly in $(\alpha ,\theta )$, they demonstrate that the framework can distinguish changes in short-horizon radiative behavior between multi-decadal periods using only local window statistics.

Wind variables (WS10M, WS10M_MAX) remain dominated by $\alpha =2.0$ and strong mean reversion in both intervals, consistent with near-Gaussian daily increments and rapid synoptic-scale adjustment typical of mid-latitude surface winds. The mild shift from $\alpha =2.0$ toward $\alpha =1.5$ after 2005 points to slightly heavier tails—occasional stronger gust events—but the $\theta$ profile remains heavily concentrated in the fastest reversion class, implying that such excursions decay quickly at the 50-day scale.

Precipitation (PRECTOTCORR) is distinct in showing persistent heavy tails with $\alpha$ concentrated around 1.0–1.5 and appreciable counts in 0.5. This reflects the intermittency and skewed nature of rainfall series, where dry spells alternate with bursts of wet days. Its $\theta$-profile, while still anchored by high-reversion bins, includes non-trivial counts at weaker $\theta$, consistent with multi-day weather regimes or blocked patterns that maintain elevated precipitation probabilities over several days.

Thermodynamic fields—surface pressure, relative humidity, dew-point temperature, and diurnal temperature range—show the opposite signature. These variables concentrate almost entirely in the Gaussian, fast-reverting regime. Surface pressure and humidity display high $\theta$ and $\alpha =2.0$, reflecting small, symmetric day-to-day fluctuations about a slowly varying mean. Dew-point temperature behaves similarly but with slightly broader tails, while diurnal range exhibits the strongest mean-reversion of all variables, with more than $5.1\times {10}^5$ windows classified at $\theta =32$ in each period, indicating extremely rapid recovery of short-term anomalies.

Across all variables, three consistent themes emerge. Deterministically forced quantities such as top-of-atmosphere radiation remain structurally heavy-tailed and weakly reverting at this window scale. Weather-modulated fields—radiation near the surface, wind, and precipitation—occupy intermediate regimes where the balance between intermittency and reversion depends on local processes and changes only subtly between the two periods. Thermodynamic quantities remain stably Gaussian and fast-reverting, with minimal redistribution across decades. These patterns conform to established statistical characteristics of daily meteorological data and show that NASA POWER, with harmonized MERRA-2 inputs, reproduces these relations consistently across space and time. In this sense, the $(\alpha ,\theta )$ representation functions as a compact tool for analyzing short-horizon variability and for inferring regime shifts from local window statistics alone.

9. Discussion

Our goal is to recover two local, window–scale descriptors of dynamics, the stability $\alpha$ (tail heaviness) and the mean-reversion rate $\theta$, and to show that the same trained classifiers are useful across synthetic benchmarks and three real-world domains. Methodologically, the approach is best viewed as a tool: given a univariate time series and a choice of window length, it produces a categorical description of tail behavior and recovery speed without re-estimating a full stochastic model for each segment. The synthetic experiments establish that this tool can separate classes on controlled data; the financial, solar, and climate applications then illustrate how the same instrument behaves in applied settings, where “ground truth’’ is partial and domain knowledge must guide interpretation.

On synthetic, dimensionless Ornstein Uhlenbeck processes with $\alpha$-stable increments, we train separate CatBoost models for $\alpha$ and $\theta$ on windows of length $w\in \{50,100,250\}$. The multi-class reports in Table 1, Table 2 and Table 3 and the confusion matrices in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 show consistent generalization: accuracies increase with w, errors are predominantly between adjacent classes, and extremes are easiest to identify (e.g., $\alpha =2$ and $\alpha =0.05$; $\theta ={10}^{-6}$ and $\theta =32$). The auxiliary binary evaluation (Appendix A) confirms that the same models, without retraining, act as reliable screeners: “Gaussian vs. Lévy’’ for $\alpha$ reaches roughly 94–96% accuracy at $w=50$ and higher at longer windows, while “no_mean_rev vs. mean_rev’’ for $\theta$ rises from about 89% ( $w=50$ ) to about 95% ( $w=250$ ). These outcomes provide a controlled baseline for interpreting the real-data applications.

Applied to daily financial series (AAPL, MSFT, DJI, GSPC; Section 6), the window-level $\alpha$ and $\theta$ maps indicate a post-2010 redistribution toward heavier effective tails and fewer occurrences of the strongest mean-reversion band in the broad indices, with more modest shifts for the two stocks. Concretely, $\alpha =2.0$ counts fall for GSPC and DJI between 1995–2004 and 2010–2024, with mass moving largely to $\alpha =1.5$ (Table 4); in parallel, $\theta =32$ counts decline for the indices and for AAPL, while MSFT remains closer to flat (Table 5). The annual traces in Figure 10 and Figure 11 visualize the same pattern. This is compatible with established “stylized facts’’ of asset returns, such as fat tails, volatility clustering, and scale-dependent dependence [5,6,52], and it does not contradict informational efficiency [42,43]. Rather, it echoes classic variance-ratio evidence that equity returns are positively autocorrelated at short horizons but negatively autocorrelated at multi-year horizons, implying a sizable transitory component and mean reversion over longer windows [8,52,53]. In our taxonomy, more windows labeled as heavy-tailed and fewer labeled as very fast-reverting mean that outlier moves occur more often and take longer to dissipate, especially in the indices after 2010.

A distinct episode visible across all annual $\alpha$ and $\theta$ traces is the COVID-19 period (2020–2021). In each of the MSFT, GSPC, DJI, and AAPL panels, the onset of the pandemic coincides with a sharp increase in windows classified as heavy-tailed and a concurrent reduction in strongly mean-reverting windows, reflecting the abrupt and persistent volatility associated with the global shutdown and subsequent market reactions. The $\alpha$-based counts show a temporary shift toward the Lévy side, while the $\theta$-based counts display fewer fast-reversion labels, consistent with slower dissipation of shocks during crisis conditions. This pattern is observable in all four assets despite their different long-term trajectories, indicating that the classifiers capture a coherent, system-wide regime disturbance. The clear COVID-19 signal illustrates that the method can detect short-horizon structural breaks and extreme periods directly from window-level dynamics, reinforcing its suitability as a tool for identifying regime changes in high-frequency financial time series.

Here the role of the tool is exploratory: it condenses large price archives into local “regime maps’’ that are consistent with known properties but are not tied to a single structural event.

For daily sunspot numbers (SILSO; Section 7), the two-interval comparison (1985–2004 vs. 2005–2024) shows a clear shift away from Gaussian-like windows toward heavier tails and a concurrent weakening of mean reversion. Interval counts move from 4805 Gaussian vs. 2500 Lévy in the time span from 1985 to 2004 to 3730 vs. 3575 in the next period from 2005 to 2024 (Table 6), while no_mean_rev rises from 2216 to 3311 windows (Table 7). The annual binary plots in Figure 12 and Figure 13 show that Gaussian and mean-reverting windows cluster near cycle maxima, whereas Lévy and non-reverting windows increase near minima. Interpreting the pair $(\alpha ,\theta )$ together: higher-activity cycles yield many small, quickly corrected fluctuations (more Gaussian, more mean reversion), while the weaker and more irregular cycles after 2005 present longer durations interrupted by intermittent bursts (heavier tails, weaker reversion).

For the Austrian NASA POWER dataset (Section 8), the accumulated counts across locations and variables (Table 8 and Table 9) separate three clear classes of behavior. The top-of-atmosphere shortwave downward irradiance (TOA_SW_DWN), which represents the incoming solar radiation before any atmospheric interaction, is governed almost entirely by deterministic orbital geometry. As a consequence, its 50-day windows fall predominantly into heavy-tailed, weak-reversion bins: within such short windows, the signal follows the seasonal cycle and exhibits slow drift with little short-term corrective dynamics. Surface-level shortwave irradiance under all-sky conditions (ALLSKY_SFC_SW_DWN) and clear-sky conditions (CLRSKY_SFC_SW_DWN), together with near-surface wind speed (WS10M, WS10M_MAX) and daily precipitation totals (PRECTOTCORR), reflect weather-modulated processes. Surface irradiance and wind fields show largely Gaussian, fast-reverting windows, indicating frequent short-term adjustments driven by synoptic variability. In contrast, precipitation exhibits persistent heavy tails and non-negligible counts in weaker-reversion classes, consistent with the intermittent alternation between dry spells and wet episodes typical of mid-latitude rainfall. Thermodynamic variables—surface pressure (PS), relative humidity at 2 m (RH2M), dew-point temperature (T2MDEW)—and the diurnal temperature range (T2M_RANGE) exhibit the opposite pattern. These quantities remain almost exclusively in the Gaussian, strongly mean-reverting regime, reflecting tight physical constraints and rapid synoptic adjustment. Their distributions show little redistribution across decades, consistent with their stable dynamical control. A notable finding is the redistribution in clear-sky surface shortwave irradiance (CLRSKY_SFC_SW_DWN) around the year 2000, visible in the annual binary splits (Figure 15). After 2000, windows labeled as Gaussian and mean-reverting become more frequent. This shift aligns with the period counts and binary summaries (Table A6 and Table A7). While not attributable to a single documented forcing event, the pattern illustrates how the method distinguishes between (i) deterministic seasonal drift (TOA irradiance), (ii) weather-driven variability (surface irradiance, winds, precipitation), and (iii) strongly regulated thermodynamic fields, using one consistent window-level diagnostic framework.

Across these three use cases, the role of the examples is therefore complementary. Financial series test whether the method recovers and refines well-known stylized facts. Sunspots provide a case where independent physical knowledge of cycle changes allows a qualitative external check. The Austrian climate fields illustrate how the same classifiers behave on multi-variable, spatially distributed data and how they separate deterministic and stochastic components. In all three settings, the tool does not “prove’’ specific causal stories; rather, it produces compact, interpretable summaries that can be compared with external knowledge and used to formulate hypotheses about regime changes, resilience, and control.

Across all domains, classifying $\theta$ and $\alpha$ provides a operational summary: treat $\theta$ as a recovery meter and $\alpha$ as a tail-risk meter. Trends toward lower $\theta$ with stable or heavier tails warn of slower recovery and more persistent deviations; trends toward higher $\theta$ with lighter tails indicate stabilization. Mixed signals call for domain checks (e.g., coupling changes, controls, measurement). This framing is consistent with early-warning theory in complex systems: as systems approach critical transitions, recovery slows, autocorrelation rises, and variance inflates [54,55,56]. Our window-level $\theta$ labels operationalize recovery speed, while $\alpha$ labels track the potential for large shocks; together they translate information about a system’s resilience into measurable segments. Finally, our sensitivity analyses on synthetic data highlight an inherent limitation: finite windows impose unavoidable uncertainty, particularly for intermediate classes where local statistics overlap; increasing w reduces but does not eliminate near-neighbor confusions. In sum, the synthetic benchmarks show separability and accuracy of our approach, and the three real-world applications show that the same $(\alpha ,\theta )$ classifiers produce interpretable, domain-consistent classifications of tail behavior and recovery speed across markets, solar activity, and surface climate.

9.1. Larger-Scope Interpretation: Linking $\alpha$ and $\theta$ to Real-World Systems

Our two parameters, the stability index $\alpha$ and the mean-reversion rate $\theta$, summarize complementary dimensions of local dynamics: how large typical fluctuations are, and how quickly a system returns toward equilibrium after a disturbance. Reading them together allows us to link the window-level patterns found in financial, solar, and climatic time series to broader mechanisms of variability, resilience, and control across complex systems.

When $\theta$ is large, perturbations decay rapidly and the system quickly reverts to baseline; when $\theta$ is small, deviations persist, and recovery slows. In the literature on critical slowing down, such declines in recovery speed and increases in short-lag autocorrelation have been interpreted as early-warning indicators of approaching tipping points [54,55,56]. In our analysis, extended runs of persistent mean reversion windows with declining effective $\theta$ therefore suggest weakening resilience, whereas stable or increasing $\theta$ implies a regime of fast equilibration. The results across all domains fit this intuition: synthetic and real systems with strong restoring forces, e.g., atmospheric pressure or diurnal temperature range, remain predominantly high-$\theta$; those with slower or state-dependent feedbacks, e.g. precipitation, financial volatility, or solar activity—show more mass in weak-reversion bins and longer persistence of deviations.

The stability index $\alpha$ quantifies the thickness of the tails: values near 2 indicate Gaussian-like noise with rare extremes, while smaller $\alpha$ mark increasingly heavy tails, with higher probability of large jumps [5,6]. Heavy-tailed windows capture episodes where the variance of innovations is effectively infinite, implying bursts, jumps, or clustered extremes rather than steady diffusion. In markets, this corresponds to volatility clustering and drawdowns; in solar or climatic series, it maps to intermittent high-amplitude events (solar storms, intense precipitation). A shift toward higher $\alpha$ may signal stabilization or averaging (e.g., due to stronger controls or measurement smoothing), while a decline toward lower $\alpha$ points to more erratic, bursty regimes. We therefore interpret $\alpha$ only jointly with $\theta$: persistent heavy tails with low $\theta$ imply systems that both generate and retain large deviations, whereas heavy tails with high $\theta$ indicate frequent shocks but rapid recovery.

The financial results follow this logic. Equity returns show the well-documented pattern of positive short-horizon autocorrelation and negative long-horizon autocorrelation, implying transitory components and mean reversion [8,42,43,52,53]. In our maps, such behavior appears as alternating patches of mean_rev and no_mean_rev windows. The increase in heavy-tailed, weakly reverting intervals in indices after 2010 (Section 6) corresponds to periods of higher volatility and slower dissipation of shocks—consistent with post-crisis market structure changes and the persistence of macro volatility. Individual stocks, particularly AAPL and MSFT, show similar but less pronounced effects, suggesting partial insulation due to firm-specific growth factors.

The sunspot record provides a natural analog with a clear physical driver. Around each 11-year solar cycle maximum, the series exhibits many small, quickly corrected fluctuations, producing high $\alpha$ and high $\theta$; near minima, activity becomes intermittent and clustered, lowering both parameters. Our counts and annual traces (Section 7) show that after 2005, the weaker, more irregular cycles coincide with a rise in Lévy and non-reverting windows, indicating a slower and more uneven recovery of solar activity after shocks. This matches known asymmetries and intermittency in modern solar cycles and supports interpreting $\theta$ as a dynamical recovery metric beyond classical variance measures.

The NASA POWER results (Section 8) extend this view to terrestrial climate fields. Deterministically forced quantities such as top-of-atmosphere shortwave flux concentrate in low-$\theta$, heavy-tail bins, reflecting smooth seasonal drift with minimal stochastic correction. Near-surface variables—irradiance under clouds, winds, precipitation—occupy intermediate regimes where fast meteorological feedbacks coexist with intermittent forcing. The sharp redistribution in clear-sky surface shortwave irradiance (CLRSKY_SFC_SW_DWN) around 2000 (Figure 15) shows how large-scale radiative changes project onto local tail and reversion patterns: the post-2000 era displays more Gaussian and mean-reverting windows, potentially reflecting shifts in aerosol loading or cloud climatology. Thermodynamic fields such as pressure and humidity remain stably Gaussian and rapidly reverting, as expected under strong synoptic control. The binary appendices (Table A6 and Table A7) confirm that these categories are consistent across space and decades.

Taken together, the cross-domain synthesis suggests a simple operational reading. Treat $\theta$ as a recovery meter and $\alpha$ as a tail-risk meter. Moving toward lower $\theta$ with stable or heavier tails signals slower recovery and higher persistence of shocks—an early-warning pattern observed in many complex systems [54,55]. Moving toward higher $\theta$ with lighter tails denotes stabilization and stronger feedback. Mixed signals, such as heavier tails but faster recovery, call for domain-specific interpretation (e.g., improved buffering or new controls masking underlying volatility). This joint perspective avoids over-interpreting any single metric and links our window-level classifications to the macroscopic behavior of markets, solar cycles, and climate systems. It also clarifies the role of the method in practice: not as a stand-alone detector of specific events, but as a compact diagnostic tool that turns long time series into interpretable maps of variability, recovery, and systemic risk.

10. Conclusions

We introduced a method to estimate the local mean-reversion rate $\theta$ and the stability index $\alpha$ from time series using a feature-based, supervised learning framework trained on synthetic Ornstein–Uhlenbeck paths with $\alpha$-stable increments. The approach combines interpretable statistical descriptors with gradient-boosted classifiers to assign each short window of a real series to discrete $(\alpha ,\theta )$ categories. Synthetic benchmarks confirm that both parameters can be recovered with high accuracy and that confusion is largely confined to adjacent classes across window lengths.

Applying the same models to financial, solar, and climatic datasets shows that $\alpha$ and $\theta$ capture different aspects of system dynamics in a consistent way. In financial returns, we find a redistribution between 1995–2004 and 2010–2024 toward heavier effective tails and fewer occurrences of the strongest mean-reversion band in the broad indices, with smaller changes for AAPL and MSFT. In solar activity, the comparison between 1985–2004 and 2005–2024 reveals that the weaker and more irregular recent cycles exhibit more Lévy-type and non-mean-reverting windows, while Gaussian and mean-reverting windows cluster around earlier cycle maxima. In the Austrian NASA POWER climate data, atmosphere radiation remains structurally heavy-tailed and weakly reverting, near-surface radiative and wind fields are largely Gaussian and fast-reverting, and thermodynamic variables are stably Gaussian and strongly mean-reverting; a clear redistribution in clear-sky shortwave irradiance around 2000 illustrates how large-scale radiative changes project onto $(\alpha ,\theta )$ respectively. In all three use cases, the method does not replace detailed physical or econometric modeling but acts as a compact diagnostic that produces domain-consistent categorizations of tail behavior and recovery speed.

Conceptually, the two parameters provide a simple operational reading: $\theta$ acts as a recovery meter and $\alpha$ as a tail-risk meter. Trends toward lower $\theta$ with stable or heavier tails signal slower recovery and more persistent deviations, whereas higher $\theta$ with lighter tails indicates faster equilibration and reduced local tail risk. Because the classifiers operate directly on short windows and require only a generic feature representation, the framework is portable across domains and suitable for scanning long observational records for changes across short-horizon regimes.

This study has several limitations. First, the classifiers are trained on a specific synthetic family (dimensionless Ornstein–Uhlenbeck processes with symmetric $\alpha$-stable noise), so applications to systems with qualitatively different dynamics or strong asymmetries should be interpreted as phenomenological mappings rather than structural parameter recovery. Second, we work with finite window lengths and a discrete grid of $(\alpha ,\theta )$ classes; windows near class boundaries are intrinsically hard to distinguish, and some loss of information relative to continuous parameter estimates is unavoidable. Third, the feature set and CatBoost models were chosen for robustness and practicality, not optimality; other feature constructions or architectures might capture additional structure or yield different decision boundaries. Finally, all real-data case studies rely on external data products (financial indices, SILSO sunspot numbers, NASA POWER fields) with their own uncertainties and revisions, so our labels reflect both underlying dynamics and data-processing choices. These constraints do not invalidate the results but should be kept in mind when transferring the method to new domains or when using $(\alpha ,\theta )$ labels for operational decisions.

Future work should extend this approach to multivariate and spatiotemporal settings, integrate uncertainty estimates for the $(\alpha ,\theta )$ assignments, and investigate whether systematic shifts in the joint distribution of these labels can serve as early indicators of emerging transitions in economic, heliophysical, or climate-related systems.

The full code to reproduce our findings together with the trained models is available at

https://github.com/Raubkatz/HeavyTailsMeanReversion.

Appendix A.0.1. Binary Classification Reports 

We mirror the formatting of the multi-class section: for each window size, $\alpha$ (left) and $\theta$ (right) reports are shown side-by-side.

Table A1. Binary reports for 50-point windows.

Table A1. Binary reports for 50-point windows.

Table A2. Binary reports for 100-point windows.

Table A2. Binary reports for 100-point windows.

Table A3. Binary reports for 250-point windows.

Table A3. Binary reports for 250-point windows.

Interpretation (binary).

For $\alpha$, the classifier is a reliable Gaussian screener: already at $w=50$ it distinguishes Gaussian windows from heavy-tailed alternatives with $\approx 94\%$ accuracy and balanced macro-averages, and it improves slightly with w. This indicates that the window-level complexity features capture distributional shape robustly, with longer windows providing more opportunities to observe tail events and intermittency. For $\theta$, the binary split cleanly separates “presence” vs. “absence” of mean reversion by $w=250$ (accuracy $\approx 95\%$); at small w, weak reversion is harder to verify (lower recall for no_mean_rev), which is consistent with the temporal nature of reversion signals. Practically, these binary tests validate that the trained models can deliver crisp, interpretable flags (Gaussian vs. Lévy; strong vs. weak reversion) that align with the most consequential regime distinctions in downstream analyses.