Gaussian Processes for Non-Parametric Price Modeling

Part of the "Classic ML baselines" series.

Two things make a Gaussian process worth a separate article on this blog, and neither of them is "it gives you uncertainty."

The first is kernel design. A GP's entire inductive bias lives in one function $k(x, x')$, and that function is something you write down deliberately: how rough the path is, whether it repeats, whether the repetition decays. Nothing else in the standard toolkit lets you state a structural hypothesis about market dynamics that explicitly and then fit it. The second is the marginal likelihood — a training objective with a complexity penalty derived from the model itself, not from a held-out set. Every other article in the overfitting arc on this blog (plateau analysis, PBO, deflated Sharpe) exists because validation-set regularization is fragile under search. A GP claims not to need one. That claim is testable, and testing it is more interesting than another uncertainty-sizing tutorial.

On uncertainty itself: the GP posterior variance is structural — it falls out of the same inference that produces the mean, rather than being wrapped around a fitted model afterwards. That is the real contrast with conformal prediction, which already covers on this blog why uncertainty is the right input to position sizing and what to do with an interval once you have one. This article does not re-argue that case; it goes after the model.

What follows is the kernel and inference machinery, a GPyTorch implementation, and — stated plainly at the end — the measurements this article does not yet have.

What Is a Gaussian Process?

GPs already appear on this blog as a Bayesian-optimization surrogate in Optuna vs. coordinate descent, with the same notation and the same low-dimensionality caveat. Here the GP is the model itself, fitted to market data rather than to a hyperparameter search surface, so the treatment goes deeper.

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution. It is a distribution over functions, not a distribution over parameters.

Formally, a function $f: \mathcal{X} \to \mathbb{R}$ is drawn from a GP if for any finite set of inputs $\{x_1, x_2, \ldots, x_n\} \subset \mathcal{X}$:

$\begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{pmatrix} \sim \mathcal{N}\left(\begin{pmatrix} m(x_1) \\ m(x_2) \\ \vdots \\ m(x_n) \end{pmatrix}, \begin{pmatrix} k(x_1, x_1) & \cdots & k(x_1, x_n) \\ \vdots & \ddots & \vdots \\ k(x_n, x_1) & \cdots & k(x_n, x_n) \end{pmatrix}\right)$

where $m(x) = \mathbb{E}[f(x)]$ is the mean function and $k(x, x') = \text{Cov}(f(x), f(x'))$ is the covariance (kernel) function. We write this compactly as:

$f \sim \mathcal{GP}(m(\cdot), k(\cdot, \cdot))$

The mean function encodes prior belief about the average behavior of $f$. In trading, we typically set $m(x) = 0$, encoding the assumption that we have no prior directional bias on returns. All the structure goes into the kernel.

Why Non-Parametric?

A linear regression model with 5 features has 6 parameters. A neural network with two hidden layers of 64 units has thousands. A GP has no fixed number of parameters — the complexity of the model grows with the data. With 10 observations, the GP defines a 10-dimensional Gaussian. With 10,000 observations, it defines a 10,000-dimensional Gaussian.

This does not mean GPs have no hyperparameters. The kernel function has hyperparameters (length scales, amplitudes, periodicities) that control the properties of functions drawn from the prior. But the functional form itself is never fixed. The GP can represent any continuous function, given enough data and the right kernel. This is what "non-parametric" means — the model is not constrained to a parametric family like linear functions or polynomials.

For financial modeling, this is valuable. Markets change. The relationship between features and returns is nonlinear, non-stationary, and regime-dependent. Parametric models impose structure that may not match reality. GPs let the data speak.

Kernel Functions: Encoding Market Structure

The kernel function $k(x, x')$ is the soul of a Gaussian process. It defines which functions are probable a priori by specifying the covariance between function values at any two input points. Different kernels encode different assumptions about smoothness, periodicity, and long-range behavior.

Radial Basis Function (RBF) / Squared Exponential

The RBF kernel is the most common starting point:

$k_{\text{RBF}}(x, x') = \sigma^2 \exp\left(-\frac{\|x - x'\|^2}{2\ell^2}\right)$

where $\sigma^2$ is the signal variance (output scale) and $\ell$ is the length scale. Functions drawn from a GP with an RBF kernel are infinitely differentiable — very smooth.

Trading interpretation: The length scale $\ell$ controls how far apart two data points can be and still be correlated. A short length scale means the model reacts to local patterns; a long length scale means it captures broad trends. The signal variance $\sigma^2$ controls the amplitude of the function — how large the predicted returns can be.

Problem for finance: Infinite smoothness is unrealistic. Financial returns have jumps, regime changes, and discontinuities. The RBF kernel can oversmooth these features, producing predictions that are too conservative near structural breaks.

Matern Kernel

The Matern class generalizes the RBF by introducing a smoothness parameter $\nu$:

$k_{\text{Matern}}(x, x') = \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{\sqrt{2\nu}\|x - x'\|}{\ell}\right)^{\nu} K_{\nu}\left(\frac{\sqrt{2\nu}\|x - x'\|}{\ell}\right)$

where $K_{\nu}$ is the modified Bessel function of the second kind. As $\nu \to \infty$, the Matern kernel converges to the RBF. Common choices:

• $\nu = 1/2$: Equivalent to the Ornstein-Uhlenbeck process. Functions are continuous but not differentiable — rough, like Brownian motion.
• $\nu = 3/2$: Functions are once differentiable. A good balance between smoothness and flexibility.
• $\nu = 5/2$: Functions are twice differentiable. Smoother than $3/2$ but less rigid than RBF.

Trading interpretation: The Matern-$3/2$ kernel is arguably the best default for financial time series. It allows the kind of roughness that real price paths exhibit without being as jagged as $\nu = 1/2$. This aligns with the "volatility is rough" literature (Gatheral, Jaisson, & Rosenbaum, 2018), which empirically shows that volatility paths have Hurst exponents around $H \approx 0.1$, far rougher than Brownian motion.

The main published evidence for Matern kernels beating classical volatility models is Rizvi et al. (2017), who report roughly 20% better MSE than a random walk and 50% better than GARCH — on 2017 daily currency-pair data, not crypto, and not reproduced here. Treat it as motivation for trying the kernel, not as a benchmark. This blog's own GARCH(1,1) fit on real BTC daily data, with Ljung-Box and ARCH-LM diagnostics, is in GARCH volatility forecasting for crypto; a head-to-head against a Matern GP on the same sample would be the honest comparison, and it has not been run.

Periodic Kernel

Financial markets have cyclical patterns: intraday volume curves, day-of-week effects, monthly rebalancing flows, quarterly earnings seasons. The periodic kernel captures these:

$k_{\text{Periodic}}(x, x') = \sigma^2 \exp\left(-\frac{2\sin^2\left(\pi|x - x'|/p\right)}{\ell^2}\right)$

where $p$ is the period. Functions drawn from this kernel repeat with period $p$, modulated by the length scale $\ell$ which controls how quickly the correlation decays within a period.

Trading interpretation: Set $p = 24$ (hours) to capture intraday patterns, or $p = 5$ (trading days) for weekly seasonality. Unlike Fourier features, the periodic kernel does not assume a fixed number of harmonics — the GP learns the shape of the cycle from data.

Combining Kernels: Additive and Multiplicative Composition

The real power of GP kernels lies in composition. If $k_1$ and $k_2$ are valid kernels, so are:

• Sum: $k_1 + k_2$ — the function is the sum of independent components (additive decomposition)
• Product: $k_1 \times k_2$ — interactions between components (e.g., locally periodic behavior)

A useful composite kernel for financial returns:

$k(x, x') = k_{\text{Matern-3/2}}(x, x') + k_{\text{Periodic}}(x, x') \cdot k_{\text{RBF}}(x, x')$

This decomposes the signal into:

1. A non-smooth, aperiodic trend component (Matern-3/2)
2. A periodic component whose amplitude decays over time (Periodic $\times$ RBF)

The product $k_{\text{Periodic}} \cdot k_{\text{RBF}}$ creates a locally periodic kernel: the pattern repeats, but distant repetitions have less influence than nearby ones. This is exactly right for financial seasonality, which drifts over time as market microstructure evolves.

Spectral Mixture Kernels

For maximum flexibility, the spectral mixture (SM) kernel (Wilson & Adams, 2013) parameterizes the spectral density of the kernel as a mixture of Gaussians:

$k_{\text{SM}}(x, x') = \sum_{q=1}^{Q} w_q \exp\left(-2\pi^2 \|x - x'\|^2 v_q\right) \cos\left(2\pi \|x - x'\| \mu_q\right)$

where $w_q$ are mixture weights, $v_q$ are spectral variances, and $\mu_q$ are spectral means (frequencies). By Bochner's theorem, any stationary kernel can be represented this way. The SM kernel can discover periodic components, long-range trends, and short-range correlations simultaneously — all from data.

Trading interpretation: The SM kernel is useful when you do not know what patterns exist in the data. It can identify hidden periodicities in return series (e.g., subtle 4-hour cycles in crypto markets driven by automated rebalancing). The downside is more hyperparameters and the risk of overfitting with small datasets.

Posterior Inference: From Prior to Prediction

Given training data $\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n$ where $y_i = f(x_i) + \epsilon_i$ and $\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$, the GP posterior at test points $X_*$ has a closed-form solution. This is the key computational advantage of GPs over most Bayesian models.

The Posterior Equations

Let $K = k(X, X)$ be the $n \times n$ training covariance matrix, $K_* = k(X_*, X)$ be the $m \times n$ cross-covariance matrix, and $K_{**} = k(X_*, X_*)$ be the $m \times m$ test covariance matrix. The posterior is:

$f_* | X, y, X_* \sim \mathcal{N}(\bar{f}_*, \text{Cov}(f_*))$

where:

$\bar{f}_* = K_* (K + \sigma_n^2 I)^{-1} y$

$\text{Cov}(f_*) = K_{**} - K_* (K + \sigma_n^2 I)^{-1} K_*^T$

The posterior mean $\bar{f}_*$ is a linear combination of the training targets, weighted by the kernel similarity between test and training points. The posterior covariance $\text{Cov}(f_*)$ starts from the prior covariance $K_{**}$ and subtracts the information gained from the training data. Where training data is dense, the posterior variance is small. Where training data is sparse, the posterior variance reverts to the prior.

The Marginal Likelihood, and the Claim Worth Testing

The kernel hyperparameters $\theta$ (length scales, variances, noise level) are learned by maximizing the log marginal likelihood:

$\log p(y | X, \theta) = -\frac{1}{2} y^T (K + \sigma_n^2 I)^{-1} y - \frac{1}{2} \log |K + \sigma_n^2 I| - \frac{n}{2} \log 2\pi$

The first term is a data fit term (penalizes predictions far from observations). The second term is a complexity penalty (penalizes models that are too flexible — i.e., where the kernel matrix has large determinant). The third term is a normalization constant.

This is the automatic Occam's razor, and it is the most interesting thing a GP brings to a trading pipeline. The strong version of the claim is that no separate validation set is needed for regularization — the complexity penalty is inside the objective, so the model cannot buy fit with flexibility for free.

That claim deserves scepticism on this blog specifically. Plateau analysis shows that a single-point validation score is a bad selection criterion and that robustness lives in the shape of the neighborhood; PBO quantifies how often the in-sample winner loses out-of-sample; deflated Sharpe prices in the number of trials. A marginal likelihood is still an in-sample objective being maximized over $\theta$ — the Occam factor penalizes model capacity, not selection over many fitted kernels. If you fit twelve candidate kernels and pick the one with the best marginal likelihood, you are back in multiple-testing territory and the deflation logic applies unchanged. The falsifiable version: does marginal-likelihood selection produce a smaller in-sample/out-of-sample gap than validation-set selection on the same data and the same kernel family? That is measurable and is not measured below.

The marginal likelihood surface also has local optima. Multiple random restarts or careful initialization matter — initializing the length scale to the median pairwise distance of the training inputs and the noise variance to the sample variance of the targets is a reasonable starting point.

Computational Cost and Scalability

The bottleneck is inverting $(K + \sigma_n^2 I)$, which costs $O(n^3)$ in time and $O(n^2)$ in memory. The cubic wall is already argued on this blog from the other direction: search method vs. evaluation cost disqualifies GP-based Bayesian optimization outright when the objective is cheap, because the surrogate costs more than the evaluations it saves. The arithmetic is the same here; the application is different. As a model fitted to market data, the cubic term is not a disqualification but a budget: it sets a hard ceiling on the training window.

Scalability strategies for trading applications:

1. Sparse GPs (inducing points). Replace the $n \times n$ matrix with an $m \times m$ matrix where $m \ll n$. The inducing points $Z = \{z_1, \ldots, z_m\}$ are pseudo-inputs that summarize the training data. The SVGP (Stochastic Variational GP) formulation by Hensman et al. (2013) allows mini-batch training with cost $O(nm^2)$ per iteration. GPyTorch supports this natively.

2. Structured kernel interpolation (SKI/KISS-GP). Exploits Kronecker and Toeplitz structure in the kernel matrix when inputs lie on a grid. Reduces cost to $O(n + g \log g)$ where $g$ is the grid size. Ideal for regularly-sampled time series (e.g., 1-minute bars).

3. Local GPs on sliding windows. Train separate GPs on recent data only. This is where the GP-specific constraint bites: standard kernels are stationary ($k(x,x')$ depends only on $x - x'$) and markets are not, so the usual answer is a rolling window — but the window length is bounded above by $O(n^3)$, not just by statistics. Walk-forward optimization covers anchored vs. rolling windows, train/test lengths, and re-optimization frequency in general; for a GP, window sizing is a compute decision as much as a statistical one, and an anchored (ever-growing) window is simply not available past a few thousand points without approximation. That reframing is the practical consequence: with GPs you do not get to choose "use all the history."

GP for Return Prediction: A Practical Framework

Feature Design: Why 5-20 Features, Not 500

The general feature taxonomy for market-data ML — order book imbalance, book pressure, VPIN, Kyle's lambda, realized-volatility features, cyclical time encoding, cross-asset and funding-rate signals — is already laid out in spread modeling with machine learning; use that list.

What is GP-specific is the size of the list. Kernel methods degrade in high dimensions: distances concentrate, and a stationary kernel loses discrimination. Practical range for a financial GP is 5-20 inputs, not the hundreds a gradient-boosted tree happily eats. The mechanism that makes this survivable is ARD (Automatic Relevance Determination): give each input dimension its own length scale $\ell_d$, and marginal-likelihood training pushes $\ell_d \to \infty$ for dimensions that carry no signal, since an infinite length scale means the kernel ignores that coordinate. Feature selection becomes a byproduct of fitting, and the learned $\ell_d$ are directly readable as a relevance ranking — which is also what makes it falsifiable: fit the composite kernel on real bars, print the length scales, and see whether the features you believe in are the ones the model keeps.

Position Sizing and Abstention

The GP posterior gives $\sigma_*(x)$ directly, so it drops straight into the edge-ratio sizing and no-trade filter developed in Conformal Prediction for Risk-Aware Position Sizing, with the interval width $w_t = 2\kappa\sigma_*(x)$. The only difference is provenance: the GP produces the width from the model itself rather than from a calibration set, so it varies with where the test point sits relative to the training data rather than with a global quantile of past residuals.

Implementation with GPyTorch

GPyTorch is a PyTorch-based library for scalable GP inference. It leverages GPU acceleration, automatic differentiation, and modern linear algebra techniques (conjugate gradients, Lanczos decomposition) to scale GPs beyond the naive $O(n^3)$ limit.

Basic Exact GP for Return Prediction

import torch
import gpytorch
import numpy as np
from torch.utils.data import TensorDataset, DataLoader


class ExactGPModel(gpytorch.models.ExactGP):
    """Exact GP with a composite kernel for financial returns."""

    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()

        self.covar_matern = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.MaternKernel(nu=1.5, ard_num_dims=train_x.shape[1])
        )
        self.covar_periodic = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.PeriodicKernel()
        )
        self.covar_rbf_decay = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )

    def forward(self, x):
        mean = self.mean_module(x)
        covar = self.covar_matern(x) + self.covar_periodic(x) * self.covar_rbf_decay(x)
        return gpytorch.distributions.MultivariateNormal(mean, covar)

ard_num_dims is what enables the per-dimension length scales discussed above. After training, model.covar_matern.base_kernel.lengthscale is the vector to read out.

Training Loop

def train_gp(train_x, train_y, n_epochs=200, lr=0.05, device=None):
    """Train the GP by maximizing the marginal log-likelihood.

    Returns the device alongside the model so that callers move test
    tensors to the same place -- otherwise prediction crashes on GPU.
    """
    if device is None:
        device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

    train_x = train_x.to(device)
    train_y = train_y.to(device)

    likelihood = gpytorch.likelihoods.GaussianLikelihood().to(device)
    model = ExactGPModel(train_x, train_y, likelihood).to(device)

    model.train()
    likelihood.train()

    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

    losses = []
    for epoch in range(n_epochs):
        optimizer.zero_grad()
        output = model(train_x)
        loss = -mll(output, train_y)
        loss.backward()
        optimizer.step()
        losses.append(loss.item())

        if (epoch + 1) % 50 == 0:
            noise = likelihood.noise.item()
            print(
                f"Epoch {epoch+1}/{n_epochs} | "
                f"Loss: {loss.item():.4f} | "
                f"Noise: {noise:.6f}"
            )

    return model, likelihood, device, losses

Prediction with Uncertainty

def predict_with_uncertainty(model, likelihood, test_x, device):
    """Generate predictions with uncertainty estimates."""
    model.eval()
    likelihood.eval()

    test_x = test_x.to(device)

    with torch.no_grad(), gpytorch.settings.fast_pred_var():
        posterior = likelihood(model(test_x))

        mean = posterior.mean
        variance = posterior.variance
        lower, upper = posterior.confidence_region()  # 2-sigma bounds

    return {
        "mean": mean.cpu().numpy(),
        "std": variance.sqrt().cpu().numpy(),
        "lower_2sigma": lower.cpu().numpy(),
        "upper_2sigma": upper.cpu().numpy(),
    }

The fast_pred_var() context manager uses the LOVE (Lanczos Variance Estimates) algorithm to compute predictive variances in $O(n)$ time instead of $O(n^2)$.

End-to-End Trading Pipeline

Note on the feature builder below: every rolling statistic must be strictly backward-looking, and the input standardization must be fitted on the training slice only. That second point is not generic hygiene — it is precisely the normalization leakage channel dissected and measured in the look-ahead bias taxonomy, which reports the Sharpe inflation each leak type produces. A GP standardizes its inputs by construction, so this is the leak it is most exposed to.

import pandas as pd


def build_features(df: pd.DataFrame, lookback: int = 10) -> pd.DataFrame:
    """Build features for GP-based return prediction."""
    features = pd.DataFrame(index=df.index)

    for lag in range(1, lookback + 1):
        features[f"ret_lag_{lag}"] = df["close"].pct_change().shift(lag)

    ret = df["close"].pct_change()
    features["vol_ratio"] = (
        ret.rolling(10).std() / ret.rolling(50).std()
    )

    features["vol_zscore"] = (
        (df["volume"] - df["volume"].rolling(50).mean())
        / df["volume"].rolling(50).std()
    )

    if "bid_vol" in df.columns and "ask_vol" in df.columns:
        features["obi"] = (
            (df["bid_vol"] - df["ask_vol"])
            / (df["bid_vol"] + df["ask_vol"])
        )

    if hasattr(df.index, "hour"):
        hours = df.index.hour + df.index.minute / 60.0
        features["time_sin"] = np.sin(2 * np.pi * hours / 24)
        features["time_cos"] = np.cos(2 * np.pi * hours / 24)

    features.dropna(inplace=True)
    return features


def run_gp_strategy(
    df: pd.DataFrame,
    train_window: int = 500,
    retrain_every: int = 50,
    confidence_threshold: float = 1.0,
    risk_fraction: float = 0.02,
    max_leverage: float = 1.0,
):
    """Walk-forward GP trading strategy with uncertainty-based sizing."""
    features = build_features(df)
    returns = df["close"].pct_change().reindex(features.index)
    target = returns.shift(-1)  # predict next-bar return

    mask = features.notna().all(axis=1) & target.notna()
    features = features[mask]
    target = target[mask]

    positions = pd.Series(0.0, index=features.index)
    predictions = pd.DataFrame(
        index=features.index, columns=["mean", "std"], dtype=float
    )

    model = likelihood = device = None
    x_mean = x_std = y_mean = y_std = None

    for i in range(train_window, len(features)):
        if model is None or (i - train_window) % retrain_every == 0:
            train_x = torch.tensor(
                features.iloc[i - train_window : i].values,
                dtype=torch.float32,
            )
            train_y = torch.tensor(
                target.iloc[i - train_window : i].values,
                dtype=torch.float32,
            )

            x_mean, x_std = train_x.mean(0), train_x.std(0) + 1e-8
            y_mean, y_std = train_y.mean(), train_y.std() + 1e-8
            train_x_norm = (train_x - x_mean) / x_std
            train_y_norm = (train_y - y_mean) / y_std

            model, likelihood, device, _ = train_gp(
                train_x_norm, train_y_norm, n_epochs=100
            )

        test_x = torch.tensor(
            features.iloc[i : i + 1].values, dtype=torch.float32
        )
        test_x_norm = (test_x - x_mean) / x_std

        pred = predict_with_uncertainty(model, likelihood, test_x_norm, device)

        pred_mean = pred["mean"][0] * y_std.item() + y_mean.item()
        pred_std = pred["std"][0] * y_std.item()

        predictions.iloc[i] = [pred_mean, pred_std]

        z_score = abs(pred_mean) / (pred_std + 1e-8)
        if z_score > confidence_threshold:
            size = min(z_score * risk_fraction, max_leverage)
            positions.iloc[i] = np.sign(pred_mean) * size

    strategy_returns = positions.shift(1) * returns
    return strategy_returns, predictions, positions

Scaling to Larger Datasets with Sparse GPs

When the training window exceeds a few thousand points, exact GP inference becomes slow. Switch to a variational sparse GP:

class SparseGPModel(gpytorch.models.ApproximateGP):
    """Sparse variational GP for large-scale return prediction."""

    def __init__(self, inducing_points):
        variational_distribution = (
            gpytorch.variational.CholeskyVariationalDistribution(
                inducing_points.size(0)
            )
        )
        variational_strategy = (
            gpytorch.variational.VariationalStrategy(
                self,
                inducing_points,
                variational_distribution,
                learn_inducing_locations=True,
            )
        )
        super().__init__(variational_strategy)
        self.mean_module = gpytorch.means.ZeroMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.MaternKernel(nu=1.5)
        )

    def forward(self, x):
        mean = self.mean_module(x)
        covar = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean, covar)


def train_sparse_gp(train_x, train_y, n_inducing=128, n_epochs=50, batch_size=256):
    """Train sparse GP with mini-batch stochastic variational inference."""
    indices = torch.randperm(train_x.size(0))[:n_inducing]
    inducing_points = train_x[indices]

    model = SparseGPModel(inducing_points)
    likelihood = gpytorch.likelihoods.GaussianLikelihood()

    model.train()
    likelihood.train()

    optimizer = torch.optim.Adam(
        [{"params": model.parameters()}, {"params": likelihood.parameters()}],
        lr=0.01,
    )
    mll = gpytorch.mlls.VariationalELBO(
        likelihood, model, num_data=train_y.size(0)
    )

    dataset = TensorDataset(train_x, train_y)
    loader = DataLoader(dataset, batch_size=batch_size, shuffle=True)

    for epoch in range(n_epochs):
        for x_batch, y_batch in loader:
            optimizer.zero_grad()
            output = model(x_batch)
            loss = -mll(output, y_batch)
            loss.backward()
            optimizer.step()

    return model, likelihood

With 128 inducing points, the per-batch cost is $O(128^2 \times \text{batch\_size})$ — roughly 4 million operations per batch. This easily handles datasets of 100,000+ observations on a single GPU.

Deep Kernel Learning: GP Meets Neural Networks

When the input space is high-dimensional or the relationship between features and returns is highly nonlinear, a plain kernel may struggle. Deep kernel learning (DKL) passes the inputs through a neural network before applying the GP kernel:

$k_{\text{DKL}}(x, x') = k_{\text{base}}(g_\phi(x), g_\phi(x'))$

where $g_\phi$ is a neural network with parameters $\phi$ and $k_{\text{base}}$ is a standard kernel (e.g., Matern-3/2). The network learns a feature representation where the GP kernel is most effective. The entire model — network parameters and kernel hyperparameters — is trained end-to-end by maximizing the marginal likelihood.

class DeepKernelGP(gpytorch.models.ExactGP):
    """GP with a neural network feature extractor."""

    def __init__(self, train_x, train_y, likelihood, input_dim):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()

        self.feature_extractor = torch.nn.Sequential(
            torch.nn.Linear(input_dim, 8),
            torch.nn.ReLU(),
            torch.nn.Linear(8, 4),
            torch.nn.ReLU(),
            torch.nn.Linear(4, 2),
        )

        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.MaternKernel(nu=1.5, ard_num_dims=2)
        )

    def forward(self, x):
        features = self.feature_extractor(x)
        mean = self.mean_module(features)
        covar = self.covar_module(features)
        return gpytorch.distributions.MultivariateNormal(mean, covar)

DKL combines the representation learning of neural networks with the uncertainty quantification of GPs. The GP layer ensures that predictions far from the training data have high uncertainty — something that standard neural networks famously fail to provide. Note also that DKL reintroduces exactly the flexibility the marginal likelihood was supposed to police: the Occam factor penalizes the kernel, but the network in front of it has thousands of free parameters and no such penalty.

Ludkovski and Risk (2025), Gaussian Process Models for Quantitative Finance, survey DKL in a broader quantitative-finance context including options pricing and portfolio optimization; those results are theirs, on their problems, and are not evidence about crypto return prediction.

Diagnostics and Pitfalls

Calibration

A well-calibrated GP has predictive intervals that match empirical coverage. The check is the same one used in conformal prediction — which also covers the theory of why marginal coverage is not conditional coverage, a limitation that applies to GP intervals just as much:

from scipy.stats import norm

def calibration_report(predictions, actuals):
    """Check if GP uncertainty is well-calibrated."""
    z_scores = (actuals - predictions["mean"]) / (predictions["std"] + 1e-8)

    for sigma in [1, 2, 3]:
        expected_outside = 2 * (1 - norm.cdf(sigma))
        actual_outside = (np.abs(z_scores) > sigma).mean()
        print(
            f"{sigma}-sigma | Expected outside: {expected_outside:.3f} | "
            f"Actual outside: {actual_outside:.3f}"
        )

Expected exceedance is 0.317 / 0.046 / 0.003 at 1/2/3 sigma. The number that matters is what this prints on a real walk-forward run over crypto bars, and that table is not in this article yet. The prior expectation is that the GP is overconfident — a Gaussian likelihood on fat-tailed, non-stationary returns should undercover badly at 3 sigma — but "should" is not a measurement.

Remaining Pitfalls

1. Input scaling. Length scales are relative to input scale, so a feature ranged $[0, 10000]$ and one ranged $[0, 1]$ cannot share a meaningful $\ell$; ARD is what rescues heterogeneous ranges, and standardization is what makes ARD's initialization sane.

2. Overfitting the marginal likelihood. With many kernel hyperparameters (especially composite or spectral mixture kernels), the marginal likelihood can still overfit. Use priors: a log-normal prior on length scales centered at the median pairwise distance, a half-normal prior on variances.

3. Covariance matrix conditioning. $(K + \sigma_n^2 I)$ can become numerically singular when the noise $\sigma_n^2$ is too small or when training points are nearly duplicated. GPyTorch adds $10^{-6}$ jitter to the diagonal; poorly-conditioned financial data often needs more.

4. Look-ahead bias. Covered above and, in measured detail, in the look-ahead bias taxonomy.

When to Use GPs vs. Other Models

* For the XGBoost and neural-network columns, defer to the published comparisons rather than to a fresh assertion here: spread modeling with machine learning has the gradient-boosting-vs-deep-learning table for financial tabular data (data-size thresholds, interpretability, regime adaptation, latency), and temporal fusion transformers has TFT vs. LSTM vs. vanilla Transformer.

Use GPs when:

• You have small-to-medium datasets (under ~10,000 observations per training window)
• You want an explicit, inspectable structural hypothesis about the signal (trend + seasonality + noise)
• Uncertainty must vary with distance from the training data, not just with a global residual quantile

Do not use GPs when:

• You need sub-millisecond inference over millions of observations
• Input dimensionality exceeds ~50
• The signal lives in complex high-order feature interactions that stationary kernels cannot represent (DKL helps, at the cost of the Occam property)

What This Article Does Not Yet Measure

Everything above is model machinery. None of it is evidence that a GP makes money on crypto, and this blog's standard is that the article carries its own numbers. The open items, in the order they should be run:

1. ARD length scales on real BTCUSDT 1m bars. Fit the composite kernel, print $\ell_d$ per feature. This directly tests the "ARD is built-in feature selection" claim and produces a falsifiable feature-relevance ranking.
2. The calibration table from a walk-forward run. Expected vs. empirical exceedance at 1/2/3 sigma, stated plainly, including the case where the GP turns out to be overconfident.
3. The confidence-gated strategy, walk-forward, on the same five majors as the honest negative, deflated for trial count. If it fails, it slots into that series as another negative result.
4. The wall-clock $O(n^3)$ curve for $n = 250 / 500 / 1000 / 2000 / 5000$, exact vs. SVGP, so the scalability section rests on a chart instead of an assertion.

Conclusion

The case for Gaussian processes in trading is not "they give you error bars" — conformal prediction gives you error bars with fewer distributional assumptions and a coverage guarantee the GP does not have. The case is that a GP makes you write your structural hypothesis down as a kernel, fits it against an objective that prices in its own complexity, and then tells you which of your features it actually used via the ARD length scales. That is an unusually legible model.

The costs are equally concrete: cubic scaling that caps your training window, stationarity assumptions that a rolling window only partly repairs, a Gaussian likelihood that fat-tailed returns will violate, and — once you reach for deep kernel learning — the quiet loss of the very Occam property that motivated the marginal likelihood in the first place. Whether what remains is tradeable is an empirical question, and the four measurements listed above are what would answer it.

References

• Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
• Wilson, A., & Adams, R. (2013). Gaussian Process Kernels for Pattern Discovery and Extrapolation. ICML.
• Hensman, J., Fusi, N., & Lawrence, N. D. (2013). Gaussian Processes for Big Data. UAI.
• Gatheral, J., Jaisson, T., & Rosenbaum, M. (2018). Volatility is rough. Quantitative Finance, 18(6).
• Rizvi, S. A. A., Roberts, S. J., Osborne, M. A., & Nyikosa, F. (2017). A Novel Approach to Forecasting Financial Volatility with Gaussian Process Envelopes. arXiv:1705.00891.
• Ludkovski, M., & Risk, J. (2025). Gaussian Process Models for Quantitative Finance. Springer.
• Gardner, J. R., Pleiss, G., Bindel, D., Weinberger, K. Q., & Wilson, A. G. (2018). GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. NeurIPS.