Koopman Operators and DMD: Do Market Modes Survive Out-of-Sample?

Dynamic Mode Decomposition hands you a spectrum: a handful of complex eigenvalues, each with a growth rate and a frequency, each attached to a spatial mode over your asset cross-section. It looks like structure. The entire question is whether it is structure or whether it is a linear operator dutifully memorising one window of noise.

That question has two testable halves, and this article is built around them:

1. Mode persistence. Fit DMD on window $t$ and on window $t+1$. Do the dominant modes span the same subspace, or do they reshuffle every refit? If they reshuffle, DMD is an in-sample decomposition and nothing more — and saying that plainly is more useful than another tutorial.
2. Spectral radius as a leading indicator. The largest eigenvalue modulus $\rho = \max_j |\lambda_j|$ is a single scalar summarising how explosive the fitted dynamics are. Does it lead realised volatility, lag it, or merely restate it? Either answer is publishable; only the first is tradeable.

That markets are non-stationary nonlinear systems is assumed here, not argued — the blog has already argued it, with phase-space geometry in attractors in algotrading and with measured per-regime BTC statistics in regime detection with HMMs. What follows is the narrower question of what non-stationarity does to a fitted Koopman operator, and how to measure it.

1. The Core Idea: Linearizing Nonlinear Dynamics

Consider a discrete-time dynamical system on a state space $\mathcal{M} \subseteq \mathbb{R}^n$:

$x_{k+1} = F(x_k)$

where $F: \mathcal{M} \to \mathcal{M}$ is a possibly nonlinear map. For markets, $x_k$ is a vector of asset returns, volatilities, or order-book imbalances at time step $k$.

The Koopman operator $\mathcal{K}$ acts not on the state $x$ directly, but on scalar-valued observable functions $g: \mathcal{M} \to \mathbb{C}$:

$[\mathcal{K} g](x) = g(F(x))$

The key property: $\mathcal{K}$ is linear, even when $F$ is not. The price is dimensionality — $\mathcal{K}$ acts on an infinite-dimensional function space. A good finite-dimensional approximation buys you the expressiveness of nonlinear dynamics with the tractability of linear algebra: prediction becomes matrix exponentiation, and every mode is inspectable rather than buried in network weights.

2. Spectral Decomposition of the Koopman Operator

If $\mathcal{K}$ has eigenvalues $\lambda_j$ and eigenfunctions $\varphi_j$, then $\mathcal{K} \varphi_j = \lambda_j \varphi_j$, and any observable $g$ in the span of these eigenfunctions decomposes as:

$g(x_k) = \sum_{j=1}^{\infty} \varphi_j(x_0) \, \lambda_j^{k} \, v_j$

where $v_j$ are the Koopman modes — vector-valued coefficients describing how each eigenfunction contributes to the full observable vector.

Each eigenvalue $\lambda_j = |\lambda_j| e^{i\omega_j}$ encodes a growth or decay rate ($|\lambda_j|$) and an oscillation frequency ($\omega_j$):

This table is the promise. Section 4 is where the promise gets checked against data.

3. Dynamic Mode Decomposition (DMD)

DMD is the workhorse algorithm for approximating $\mathcal{K}$ from data. Given snapshots arranged into matrices:

$X = \begin{bmatrix} x_0 & x_1 & \cdots & x_{m-1} \end{bmatrix}, \quad X' = \begin{bmatrix} x_1 & x_2 & \cdots & x_m \end{bmatrix}$

DMD seeks a best-fit linear operator $A$ with $X' \approx AX$:

1. Compute the SVD: $X = U \Sigma V^*$
2. Project: $\tilde{A} = U^* X' V \Sigma^{-1}$
3. Eigendecompose: $\tilde{A} W = W \Lambda$
4. Recover full-space modes: $\Phi = X' V \Sigma^{-1} W$

The columns of $\Phi$ are the DMD modes; the diagonal of $\Lambda$ holds the DMD eigenvalues.

import numpy as np
from numpy.linalg import svd, eig, lstsq

def dmd(X: np.ndarray, rank: int | None = None) -> tuple:
    """
    Dynamic Mode Decomposition.

    Parameters
    ----------
    X : np.ndarray, shape (n_features, n_snapshots)
        Data matrix where each column is a state snapshot.
    rank : int or None
        Truncation rank for the SVD. None = no truncation.

    Returns
    -------
    eigenvalues : np.ndarray, shape (r,)
        DMD eigenvalues (approximating Koopman eigenvalues).
    modes : np.ndarray, shape (n_features, r)
        DMD modes (columns), L2-normalised.
    amplitudes : np.ndarray, shape (r,)
        Mode amplitudes fitted to the FINAL snapshot, so that a one-step
        forecast is simply modes @ (eigenvalues * amplitudes).
    """
    X0 = X[:, :-1]
    X1 = X[:, 1:]

    U, S, Vh = svd(X0, full_matrices=False)

    if rank is not None:
        U = U[:, :rank]
        S = S[:rank]
        Vh = Vh[:rank, :]

    S_inv = np.diag(1.0 / S)

    A_tilde = U.conj().T @ X1 @ Vh.conj().T @ S_inv
    eigenvalues, W = eig(A_tilde)

    modes = X1 @ Vh.conj().T @ S_inv @ W

    norms = np.linalg.norm(modes, axis=0)
    norms[norms == 0] = 1.0
    modes = modes / norms

    amplitudes = lstsq(modes, X[:, -1].astype(complex), rcond=None)[0]

    return eigenvalues, modes, amplitudes

Two deliberate choices here. Modes are L2-normalised, because section 4.2 compares mode subspaces across windows and unnormalised amplitudes would swamp the comparison. Amplitudes are fitted to the last snapshot rather than the first, which means a forecast never raises an eigenvalue to a large power — the source of the numerical blow-up that makes naive DMD signals look like dynamics when they are floating-point overflow.

4. The Measurement

This is the part of the article that is not textbook. Everything above is a fitting procedure; below is the protocol for finding out whether the fit means anything.

Data. Use the project's own exchange data — a cross-section of BTC, ETH and liquid alts on a consistent minute or tick grid, not a daily download of equity ETFs. The blog is crypto-first and the microstructure argument does not transfer. Build $X$ with assets on rows and time on columns, on log returns, demeaned per asset within each window.

4.1 The eigenvalue spectrum

Report, for a representative window: how many of the $r$ eigenvalues land within a tolerance of the unit circle, what oscillation period each corresponds to in hours (period $= 2\pi / |\omega_j|$ bars, converted), and what fraction of return variance the top $r$ modes reconstruct.

def spectrum_report(eigenvalues: np.ndarray, bar_minutes: float,
                    tol: float = 0.05) -> list[dict]:
    """
    Turn a DMD spectrum into human-readable rows: modulus, period in hours,
    and whether the eigenvalue sits on the unit circle within `tol`.
    """
    rows = []
    for lam in eigenvalues:
        modulus = float(np.abs(lam))
        omega = float(np.angle(lam))
        period_hours = (2 * np.pi / abs(omega)) * bar_minutes / 60 if omega else np.inf
        rows.append({
            "modulus": modulus,
            "period_hours": period_hours,
            "on_unit_circle": abs(modulus - 1.0) < tol,
            "regime": "unstable" if modulus > 1 + tol
                      else "persistent" if abs(modulus - 1.0) <= tol
                      else "decaying",
        })
    return sorted(rows, key=lambda r: -r["modulus"])


def reconstruction_r2(X: np.ndarray, modes: np.ndarray,
                      eigenvalues: np.ndarray, amplitudes: np.ndarray) -> float:
    """Fraction of in-window return variance captured by the truncated modes."""
    n_steps = X.shape[1]
    powers = eigenvalues[:, None] ** np.arange(-(n_steps - 1), 1)
    X_hat = (modes @ (amplitudes[:, None] * powers)).real
    resid = np.var(X - X_hat)
    return 1.0 - resid / np.var(X)

An honest spectrum report is already worth the article. If nothing sits near the unit circle, there are no persistent cycles to trade and the "annual seasonality" story from the equities literature simply does not carry over to 24/7 crypto.

4.2 Mode stability across adjacent windows

The decisive test. Fit DMD on window $t$, then on window $t+1$, and measure how much of the dominant mode subspace survives. The right statistic is not a naive correlation of mode vectors — mode ordering and complex phase are arbitrary — but the principal angles between the two subspaces.

def subspace_stability(modes_a: np.ndarray, modes_b: np.ndarray,
                       k: int = 3) -> float:
    """
    Overlap between the leading-k DMD mode subspaces of two adjacent windows.

    Returns the mean cosine of the principal angles: 1.0 = identical subspace,
    0.0 = orthogonal. Immune to mode reordering and complex phase, both of
    which are arbitrary in a DMD fit.
    """
    Qa, _ = np.linalg.qr(modes_a[:, :k])
    Qb, _ = np.linalg.qr(modes_b[:, :k])
    sing = np.linalg.svd(Qa.conj().T @ Qb, compute_uv=False)
    return float(np.mean(np.clip(sing, 0.0, 1.0)))


def stability_curve(returns: np.ndarray, window: int, step: int,
                    rank: int, k: int = 3) -> np.ndarray:
    """Subspace overlap between every pair of adjacent windows."""
    fits = []
    for t_end in range(window, returns.shape[1], step):
        evals, modes, _ = dmd(returns[:, t_end - window:t_end], rank=rank)
        order = np.argsort(-np.abs(evals))
        fits.append(modes[:, order])
    return np.array([subspace_stability(fits[i], fits[i + 1], k=k)
                     for i in range(len(fits) - 1)])

Report the distribution of this overlap, and report it against a null: the same statistic computed on phase-randomised surrogates of the same returns. Overlap that is high but no higher than the surrogate null means the modes are tracking the covariance structure, not the dynamics.

4.3 Rolling spectral radius against realised volatility

The claim to test: $\rho_t = \max_j |\lambda_j|$, recomputed on rolling windows, moves before realised volatility rather than with it. The scalar is mechanically new — the blog has monitored geometric scalars in real time before, notably Kobayashi curvature in complex manifolds for algorithmic trading — but a Koopman eigenvalue modulus is a different quantity with a different failure mode, and it deserves its own lead-lag test rather than an inherited pitch.

def rolling_spectral_radius(returns: np.ndarray, window: int = 1440,
                            step: int = 60, rank: int = 5) -> dict:
    """
    Rolling DMD spectrum for regime monitoring.

    returns : np.ndarray, shape (n_assets, n_timesteps)
    window  : rolling window length in bars
    step    : bars between refits
    """
    idx, radii, dom_freq = [], [], []

    for t_end in range(window, returns.shape[1], step):
        X_win = returns[:, t_end - window:t_end]
        try:
            evals, _, _ = dmd(X_win, rank=rank)
        except np.linalg.LinAlgError:
            continue

        idx.append(t_end)
        radii.append(float(np.max(np.abs(evals))))

        on_circle = np.abs(np.abs(evals) - 1.0) < 0.1
        if on_circle.any():
            sel = evals[on_circle]
            dom_freq.append(float(np.abs(np.angle(sel[np.argmax(np.abs(sel))])) / (2 * np.pi)))
        else:
            dom_freq.append(0.0)

    return {"index": np.array(idx),
            "spectral_radius": np.array(radii),
            "dominant_frequency": np.array(dom_freq)}


def lead_lag(signal: np.ndarray, target: np.ndarray, max_lag: int = 24) -> dict:
    """
    Cross-correlation of `signal` against `target` over +/- max_lag steps.
    A peak at negative lag means the signal LEADS the target.
    """
    s = (signal - signal.mean()) / (signal.std() + 1e-12)
    y = (target - target.mean()) / (target.std() + 1e-12)
    lags = np.arange(-max_lag, max_lag + 1)
    corrs = []
    for L in lags:
        if L < 0:
            corrs.append(float(np.corrcoef(s[:L], y[-L:])[0, 1]))
        elif L > 0:
            corrs.append(float(np.corrcoef(s[L:], y[:-L])[0, 1]))
        else:
            corrs.append(float(np.corrcoef(s, y)[0, 1]))
    corrs = np.array(corrs)
    return {"lags": lags, "corr": corrs, "peak_lag": int(lags[np.argmax(np.abs(corrs))])}

Align spectral_radius to realised volatility computed on the same grid and read off peak_lag. A peak at lag 0 means $\rho_t$ is a volatility restatement with extra steps. A peak at negative lag, stable across the sample and across ranks, is the only version of this article with a tradeable claim in it.

A note on what comes after a positive result

If $\rho_t$ does lead, the obvious next move is the cross-sectional construction: rank assets by DMD-predicted next-step return, go long the predicted winners and short the predicted losers. That construction is not new here — it is the factor-residual trade already covered in statistical arbitrage and pairs trading in crypto and section 4 of complex arbitrage with vectors and matrices; the only genuinely Koopman-specific twist is that the eigenportfolios carry time-varying eigenvalues rather than static PCA loadings.

The strategy section is deliberately absent from this article because it has not been backtested here with fees and slippage. When it is, it has to clear the blog's own bar: the significance test in the Deflated Sharpe Ratio and multiple testing, against the standing counter-example of the honest negative. A DMD spectrum plot is not a result.

For the one-step forecast itself, the correct implementation predicts from the last snapshot rather than propagating an eigenvalue to the power of the window length:

def dmd_one_step(returns: np.ndarray, rank: int = 4) -> np.ndarray:
    """
    One-step-ahead prediction from the final snapshot of the window.

    Never raise eigenvalues to the window length: any |lambda| != 1 then
    overflows or underflows and the "signal" becomes numerical garbage.
    """
    evals, modes, amplitudes = dmd(returns, rank=rank)
    return (modes @ (evals * amplitudes)).real

5. Extended DMD (EDMD): Nonlinear Observables

Standard DMD operates on the raw state vector. EDMD lifts the data through a dictionary of nonlinear basis functions first.

Given a dictionary $\mathbf{D} = \{d_1, \ldots, d_p\}$ of scalar functions $d_i: \mathbb{R}^n \to \mathbb{R}$, define the lifted state:

$z_k = \begin{bmatrix} d_1(x_k) \\ d_2(x_k) \\ \vdots \\ d_p(x_k) \end{bmatrix} \in \mathbb{R}^p$

EDMD seeks $K \in \mathbb{R}^{p \times p}$ with $z_{k+1} \approx K z_k$. In the convention used by the code below — states as columns, $K$ acting on the left:

$K = A G^{-1}, \quad G = \sum_{k=0}^{m-1} z_k z_k^{T}, \quad A = \sum_{k=0}^{m-1} z_{k+1} z_k^{T}$

The dictionary is where domain knowledge enters, and the time-delay row is the Koopman-specific reason to care about delay coordinates at all: they are not a separate technique bolted on, they are one more block of the lifting map. The embedding itself — choice of delay $\tau$, embedding dimension $d$, and the reconstruction theorem behind it — is already introduced and coded in complex manifolds for algorithmic trading; take the delay vectors from there and feed them straight into build_financial_dictionary as extra rows.

import numpy as np
from itertools import combinations_with_replacement

def build_financial_dictionary(X: np.ndarray, max_poly_degree: int = 2,
                               include_volatility: bool = True,
                               delay_steps: int = 0) -> np.ndarray:
    """
    Build a dictionary of nonlinear observables for EDMD.

    X : np.ndarray, shape (n_features, n_snapshots)
    Returns Z of shape (n_dict, n_snapshots - delay_steps).
    """
    n_features, n_snapshots = X.shape
    offset = max(delay_steps, 0)
    X_eff = X[:, offset:]
    n_eff = X_eff.shape[1]

    lifted = [X_eff]  # degree-1 terms (identity)

    if max_poly_degree >= 2:
        for deg in range(2, max_poly_degree + 1):
            for combo in combinations_with_replacement(range(n_features), deg):
                term = np.ones(n_eff)
                for idx in combo:
                    term *= X_eff[idx]
                lifted.append(term.reshape(1, -1))

    if include_volatility:
        lifted.append(np.abs(X_eff))   # absolute returns
        lifted.append(X_eff ** 2)      # squared returns

    for d in range(1, delay_steps + 1):
        lifted.append(X[:, offset - d : n_snapshots - d])

    return np.vstack(lifted)


def edmd(X: np.ndarray, dictionary_fn=None, reg: float = 1e-8,
         **dict_kwargs) -> tuple:
    """
    Extended Dynamic Mode Decomposition.

    Solves Z1 ~= K @ Z0 in the least-squares sense. Uses a least-squares
    solve rather than an explicit Gram inverse: `inv` on a near-singular
    dictionary Gram matrix is how EDMD spectra get silently corrupted.
    """
    if dictionary_fn is None:
        dictionary_fn = lambda x: build_financial_dictionary(x, **dict_kwargs)

    Z = dictionary_fn(X)
    Z0, Z1 = Z[:, :-1], Z[:, 1:]

    p = Z0.shape[0]
    G = Z0 @ Z0.T + reg * np.eye(p)   # regularised Gram matrix
    A = Z1 @ Z0.T

    K = np.linalg.solve(G, A.T).T

    eigenvalues, eigenvectors = np.linalg.eig(K)
    return K, eigenvalues, eigenvectors

6. Deep Koopman Networks

The EDMD dictionary is hand-crafted, which is a real limitation when the Koopman-invariant subspace is unknown. Deep Koopman networks learn the lifting and the operator jointly.

The architecture is an autoencoder — encoder, latent code, decoder, reconstruction loss, all as introduced in anomaly detection in algotrading — with one addition that is the entire point of the section: a linearity loss forcing the latent dynamics through a single learned matrix $K$.

x_k  -->  [Encoder φ]  -->  z_k  -->  [Linear K]  -->  z_{k+1}  -->  [Decoder ψ]  -->  x̂_{k+1}
                              |                            |
                              +--- Linearity loss: ‖z_{k+1} - K z_k‖ ---+

$\mathcal{L} = \underbrace{\lVert x_{k+1} - \psi(K \varphi(x_k)) \rVert^2}_{\text{prediction}} + \alpha \underbrace{\lVert \varphi(x_{k+1}) - K \varphi(x_k) \rVert^2}_{\text{linearity}} + \beta \underbrace{\lVert x_k - \psi(\varphi(x_k)) \rVert^2}_{\text{reconstruction}}$

Without the $\alpha$ term you have an ordinary autoencoder whose latent space happens to be followed by a matrix multiply. With it, the network is penalised for any latent representation whose evolution is not linear — which is what makes the learned $K$ a Koopman approximation and its eigenvalues comparable to the DMD spectrum of section 4.

import torch
import torch.nn as nn

class DeepKoopman(nn.Module):
    """Deep Koopman autoencoder: learned lifting + linear latent dynamics."""

    def __init__(self, input_dim: int, latent_dim: int, hidden_dim: int = 128):
        super().__init__()

        self.encoder = nn.Sequential(
            nn.Linear(input_dim, hidden_dim), nn.ReLU(),
            nn.Linear(hidden_dim, hidden_dim), nn.ReLU(),
            nn.Linear(hidden_dim, latent_dim),
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, hidden_dim), nn.ReLU(),
            nn.Linear(hidden_dim, hidden_dim), nn.ReLU(),
            nn.Linear(hidden_dim, input_dim),
        )
        self.K = nn.Linear(latent_dim, latent_dim, bias=False)

    def encode(self, x: torch.Tensor) -> torch.Tensor:
        return self.encoder(x)

    def decode(self, z: torch.Tensor) -> torch.Tensor:
        return self.decoder(z)

    def forward(self, x_k: torch.Tensor) -> dict:
        z_k = self.encode(x_k)
        z_k1_pred = self.K(z_k)
        return {
            "z_k": z_k,
            "z_k1_pred": z_k1_pred,
            "x_k1_pred": self.decode(z_k1_pred),
            "x_k_recon": self.decode(z_k),
        }

    def multi_step_predict(self, x_0: torch.Tensor, n_steps: int) -> torch.Tensor:
        """Roll out by repeated application of the linear operator."""
        z = self.encode(x_0)
        preds = []
        for _ in range(n_steps):
            z = self.K(z)
            preds.append(self.decode(z))
        return torch.stack(preds, dim=1)

    def latent_spectrum(self) -> np.ndarray:
        """Eigenvalues of the learned K — directly comparable to DMD's."""
        return np.linalg.eigvals(self.K.weight.detach().cpu().numpy())


def koopman_loss(model: DeepKoopman, x_k: torch.Tensor, x_k1: torch.Tensor,
                 alpha: float = 1.0, beta: float = 0.5) -> torch.Tensor:
    out = model(x_k)
    z_k1_true = model.encode(x_k1)

    prediction = nn.functional.mse_loss(out["x_k1_pred"], x_k1)
    linearity = nn.functional.mse_loss(out["z_k1_pred"], z_k1_true)
    reconstruction = nn.functional.mse_loss(out["x_k_recon"], x_k)

    return prediction + alpha * linearity + beta * reconstruction

latent_spectrum is what makes this worth training rather than reaching for a sequence model: the learned operator is still a matrix, so the section 4.2 stability test and the section 4.3 lead-lag test apply unchanged to a deep model.

7. Practical Considerations and Pitfalls

Rank selection. The truncation rank $r$ is a bias-variance dial — too low misses dynamics, too high fits noise. Do not eyeball a singular-value elbow; the blog already answers "how many components before you are fitting noise" properly, with the Marchenko-Pastur bound from Random Matrix Theory in complex arbitrage with vectors and matrices. Keep the components whose singular values exceed the Marchenko-Pastur edge for your window shape, and state the resulting $r$ explicitly alongside every spectrum you publish.

Window length. Koopman theory assumes a fixed $F$; markets do not supply one. Rolling refits are mandatory, and section 4.2's stability curve is precisely the diagnostic for whether the chosen window is long enough to estimate and short enough to stay inside one regime.

Noise sensitivity. Financial data has a low signal-to-noise ratio and standard DMD is biased by noise in $X$. Remedies worth trying before concluding that modes are unstable:

• Total DMD (TDMD) — treats both $X$ and $X'$ as noisy via total least squares.
• Optimized DMD — directly optimises the eigenvalue-mode decomposition against the residual Frobenius norm.
• Kernel EDMD — works implicitly in a high-dimensional feature space without building the dictionary.

If mode stability rises materially under TDMD, the instability was measurement noise. If it does not, it was the market.

8. Where DMD Sits

For volatility-specific forecasting the comparison point is the GARCH family — see GARCH volatility forecasting for crypto. For fully nonlinear sequence models and the interpretability and error-accumulation trade-offs they carry, see the Temporal Fusion Transformer in trading.

The niche DMD occupies is narrow but real: multi-step forecasts generated by a single matrix power rather than autoregressive rollout, with every mode inspectable. Whether that niche contains alpha is section 4's question, not this table's.

Conclusion

Koopman theory is a genuinely elegant way to look at market dynamics, and elegance is exactly why it needs the harshest test available. The takeaways:

1. Fit DMD, then immediately test the fit. Subspace overlap between adjacent windows, measured against a phase-randomised null, tells you within an afternoon whether you have found structure or memorised a window.
2. The rolling spectral radius is the one scalar worth monitoring, and its value depends entirely on the lead-lag result. At lag 0 it is a volatility proxy; at negative lag it is a regime warning.
3. Anchor amplitudes at the last snapshot and never raise eigenvalues to the window length. A large fraction of DMD "signals" in the wild are floating-point artifacts.
4. Choose rank by Marchenko-Pastur, not by eyeballing an elbow, and publish the rank with every spectrum.
5. A spectrum plot is not a result. Any strategy built on top of this has to survive fees, slippage and a deflated Sharpe test before it counts.

For implementations beyond these, the PyDMD library covers the DMD variants comprehensively, and Mallen et al.'s reference code covers the deep Koopman side.

Markets will stay messy, non-stationary and partially observed. Koopman theory offers a principled lens for extracting structure from that mess — provided you check that the structure is still there next week.

References and further reading:

• B. O. Koopman, "Hamiltonian Systems and Transformation in Hilbert Space," Proceedings of the National Academy of Sciences, 1931.
• J. H. Tu et al., "On Dynamic Mode Decomposition: Theory and Applications," Journal of Computational Dynamics, 2014.
• M. O. Williams, I. G. Kevrekidis, C. W. Rowley, "A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition," Journal of Nonlinear Science, 2015.
• B. Lusch, J. N. Kutz, S. L. Brunton, "Deep Learning for Universal Linear Embeddings of Nonlinear Dynamics," Nature Communications, 2018.
• J. Mann and J. N. Kutz, "Dynamic Mode Decomposition for Financial Trading Strategies," Quantitative Finance, 2016.
• A. Mallen et al., "Koopman Neural Forecaster for Time Series with Temporal Distribution Shifts," ICML, 2023.
• E. Gonzalez and M. Generelo, "Analysis of chaotic economic models through Koopman operators, EDMD, Takens' theorem and Machine Learning," Data Science in Finance and Economics, 2022.