Procesos gaussianos para modelado de precios no paramétrico

Parte de la serie "Líneas básicas de ML clásico".

Hay dos cosas que hacen que un proceso gaussiano merezca un artículo aparte en este blog, y ninguna de ellas es "da incertidumbre".

El primero es diseño del kernel. Todo el sesgo inductivo de un médico de cabecera reside en una función $k(x, x')$, y esa función es algo que se escribe deliberadamente: qué tan accidentado es el camino, si se repite, si la repetición decae. Nada más en el conjunto de herramientas estándar permite plantear una hipótesis estructural sobre la dinámica del mercado que se ajuste explícitamente y luego se ajuste a ella. El segundo es la probabilidad marginal: un objetivo de entrenamiento con una penalización de complejidad derivada del modelo en sí, no de un conjunto reservado. Todos los demás artículos en el arco de sobreajuste en este blog (análisis de plateau, PBO, deflated Sharpe) existen porque la regularización del conjunto de validación es frágil bajo la búsqueda. Un médico de cabecera afirma que no lo necesita. Esa afirmación es comprobable, y probarla es más interesante que otro tutorial sobre dimensionamiento de la incertidumbre.

Sobre la incertidumbre misma: la varianza posterior del GP es estructural: surge de la misma inferencia que produce la media, en lugar de envolverse posteriormente alrededor de un modelo ajustado. Ese es el verdadero contraste con la predicción conforme, que ya explica en este blog por qué la incertidumbre es el factor correcto para determinar el tamaño de una posición y qué hacer con un intervalo una vez que se tiene uno. Este artículo no vuelve a argumentar ese caso; va tras el modelo.

Lo que sigue es el núcleo y la maquinaria de inferencia, una implementación de GPyTorch y, como se indica claramente al final, las medidas que este artículo aún no tiene.

¿Qué es un proceso gaussiano?

Los médicos de cabecera ya aparecen en este blog como un sustituto de optimización bayesiana en Optuna vs. descenso de coordenadas, con la misma notación y la misma advertencia de baja dimensionalidad. Aquí el GP es el modelo en sí, ajustado a datos de mercado en lugar de a una superficie de búsqueda de hiperparámetros, por lo que el tratamiento es más profundo.

Un proceso gaussiano es una colección de variables aleatorias, cualquier número finito de las cuales tiene una distribución gaussiana conjunta. Es una distribución sobre funciones, no una distribución sobre parámetros.

Formalmente, una función $f: \mathcal{X} \to \mathbb{R}$ se extrae de un médico de cabecera si para cualquier conjunto finito de entradas $\{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)$

dónde $m(x) = \mathbb{E}[f(x)]$ es la función media y $k(x, x') = \text{Cov}(f(x), f(x'))$ es la función de covarianza (kernel). Escribimos esto de forma compacta como:

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

La función media codifica la creencia previa sobre el comportamiento promedio de $f$. En el comercio, normalmente establecemos $m(x) = 0$, codificando el supuesto de que no tenemos un sesgo direccional previo en los rendimientos. Toda la estructura va al núcleo.

¿Por qué no paramétrico?

Un modelo de regresión lineal con 5 características tiene 6 parámetros. Una red neuronal con dos capas ocultas de 64 unidades tiene miles. Un médico de cabecera no tiene un número fijo de parámetros: la complejidad del modelo crece con los datos. Con 10 observaciones, el médico de cabecera define una gaussiana de 10 dimensiones. Con 10.000 observaciones, define una gaussiana de 10.000 dimensiones.

Esto no significa que los médicos de cabecera no tengan hiperparámetros. La función kernel tiene hiperparámetros (escalas de longitud, amplitudes, periodicidades) que controlan las propiedades de las funciones extraídas de la anterior. Pero la forma funcional en sí nunca es fija. El GP puede representar cualquier función continua, con suficientes datos y el núcleo correcto. Esto es lo que significa "no paramétrico": el modelo no está restringido a una familia paramétrica como funciones lineales o polinomios.

Para los modelos financieros, esto es valioso. Los mercados cambian. La relación entre características y rendimientos no es lineal, no estacionaria y depende del régimen. Los modelos paramétricos imponen una estructura que puede no coincidir con la realidad. Los médicos de cabecera dejan que los datos hablen.

Funciones del kernel: codificación de la estructura del mercado

La función del núcleo $k(x, x')$ es el alma de un proceso gaussiano. Define qué funciones son probables a priori especificando la covarianza entre los valores de la función en dos puntos de entrada cualesquiera. Diferentes núcleos codifican diferentes suposiciones sobre suavidad, periodicidad y comportamiento a largo plazo.

Función de base radial (RBF) / Exponencial al cuadrado

El kernel RBF es el punto de partida más común:

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

dónde $\sigma^2$ es la varianza de la señal (escala de salida) y $\ell$ es la escala de longitud. Las funciones extraídas de un GP con un kernel RBF son infinitamente diferenciables, muy fluidas.

Interpretación comercial: La escala de longitud $\ell$ controla qué tan separados pueden estar dos puntos de datos y aún estar correlacionados. Una escala de corta duración significa que el modelo reacciona a patrones locales; una escala larga significa que captura tendencias amplias. La variación de la señal. $\sigma^2$ controla la amplitud de la función: qué tan grandes pueden ser los rendimientos previstos.

Problema financiero: La suavidad infinita no es realista. Los rendimientos financieros tienen saltos, cambios de régimen y discontinuidades. El núcleo RBF puede suavizar demasiado estas características, produciendo predicciones demasiado conservadoras cerca de las rupturas estructurales.

Núcleo materno

La clase Matern generaliza el RBF introduciendo un parámetro de suavidad $\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)$

dónde $K_{\nu}$ es la función de Bessel modificada de segundo tipo. Como $\nu \to \infty$, el núcleo de Matern converge al RBF. Opciones comunes:

• $\nu = 1/2$: Equivalente al proceso de Ornstein-Uhlenbeck. Las funciones son continuas pero no diferenciables, aproximadas, como el movimiento browniano.
• $\nu = 3/2$: Las funciones alguna vez son diferenciables. Un buen equilibrio entre suavidad y flexibilidad.
• $\nu = 5/2$: Las funciones son dos veces diferenciables. Más suave que $3/2$ pero menos rígido que RBF.

Interpretación comercial: The Matern-$3/2$ Se puede decir que el kernel es el mejor valor predeterminado para series temporales financieras. Permite el tipo de aspereza que exhiben las trayectorias de precios reales sin ser tan irregular como $\nu = 1/2$. Esto se alinea con la literatura sobre "la volatilidad es aproximada" (Gatheral, Jaisson y Rosenbaum, 2018), que muestra empíricamente que las trayectorias de volatilidad tienen exponentes de Hurst alrededor $H \approx 0.1$, mucho más brusco que el movimiento browniano.

La principal evidencia publicada de que los kernels de Matern superan los modelos de volatilidad clásicos es Rizvi et al. (2017), quienes reportan aproximadamente un 20% mejor de MSE que un paseo aleatorio y un 50% mejor que GARCH - sobre datos diarios de pares de divisas de 2017, no criptográficos, y no se reproducen aquí. Trátelo como una motivación para probar el kernel, no como un punto de referencia. El ajuste GARCH(1,1) de este blog con datos diarios reales de BTC, con diagnósticos Ljung-Box y ARCH-LM, se encuentra en Pronóstico de volatilidad GARCH para criptomonedas; un enfrentamiento directo contra un médico de cabecera de Maternidad con la misma muestra sería una comparación honesta, y no se ha realizado.

Núcleo periódico

Los mercados financieros tienen patrones cíclicos: curvas de volumen intradía, efectos de día de la semana, flujos de reequilibrio mensuales, temporadas de ganancias trimestrales. El núcleo periódico captura estos:

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

dónde $p$ es el periodo. Las funciones extraídas de este núcleo se repiten con punto $p$, modulado por la escala de longitud $\ell$ que controla la rapidez con la que la correlación decae dentro de un período.

Interpretación comercial: Conjunto $p = 24$ (horas) para capturar patrones intradiarios, o $p = 5$ (días hábiles) para la estacionalidad semanal. A diferencia de las características de Fourier, el núcleo periódico no asume un número fijo de armónicos: el GP aprende la forma del ciclo a partir de los datos.

Combinando núcleos: composición aditiva y multiplicativa

El verdadero poder de los núcleos GP reside en la composición. Si $k_1$ y $k_2$ son núcleos válidos, también lo son:

• Suma: $k_1 + k_2$ — la función es la suma de componentes independientes (descomposición aditiva)
• Producto: $k_1 \times k_2$ — interacciones entre componentes (por ejemplo, comportamiento periódico local)

Un núcleo compuesto útil para los rendimientos financieros:

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

Esto descompone la señal en:

1. Un componente de tendencia aperiódico y no uniforme (Matern-3/2)
2. Un componente periódico cuya amplitud decae con el tiempo (periódico $\times$ FBR)

el producto $k_{\text{Periodic}} \cdot k_{\text{RBF}}$ crea un núcleo localmente periódico: el patrón se repite, pero las repeticiones distantes tienen menos influencia que las cercanas. Esto es exactamente cierto en el caso de la estacionalidad financiera, que varía con el tiempo a medida que evoluciona la microestructura del mercado.

Núcleos de mezcla espectral

Para una máxima flexibilidad, el núcleo de mezcla espectral (SM) (Wilson & Adams, 2013) parametriza la densidad espectral del núcleo como una mezcla de gaussianos:

$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)$

dónde $w_q$ son pesos de mezcla, $v_q$ son variaciones espectrales, y $\mu_q$ son medias espectrales (frecuencias). Según el teorema de Bochner, cualquier núcleo estacionario se puede representar de esta manera. El núcleo SM puede descubrir componentes periódicos, tendencias de largo alcance y correlaciones de corto alcance simultáneamente, todo a partir de datos.

Interpretación comercial: El kernel SM es útil cuando no se sabe qué patrones existen en los datos. Puede identificar periodicidades ocultas en series de retornos (por ejemplo, ciclos sutiles de 4 horas en criptomercados impulsados ​​por reequilibrio automatizado). La desventaja es que hay más hiperparámetros y el riesgo de sobreajustar con conjuntos de datos pequeños.

Inferencia posterior: desde antes de la predicción

Datos de entrenamiento dados $\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n$ dónde $y_i = f(x_i) + \epsilon_i$ y $\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$, el GP posterior en los puntos de prueba $X_*$ tiene una solución de forma cerrada. Ésta es la ventaja computacional clave de los médicos de cabecera sobre la mayoría de los modelos bayesianos.

Las ecuaciones posteriores

dejar $K = k(X, X)$ ser el $n \times n$ matriz de covarianza de entrenamiento, $K_* = k(X_*, X)$ ser el $m \times n$ matriz de covarianza cruzada, y $K_{**} = k(X_*, X_*)$ ser el $m \times m$ prueba de matriz de covarianza. La parte posterior es:

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

dónde:

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

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

La media posterior $\bar{f}_*$ es una combinación lineal de los objetivos de entrenamiento, ponderada por la similitud del núcleo entre los puntos de prueba y de entrenamiento. La covarianza posterior $\text{Cov}(f_*)$ comienza desde la covarianza anterior $K_{**}$ y resta la información obtenida de los datos de entrenamiento. Cuando los datos de entrenamiento son densos, la varianza posterior es pequeña. Cuando los datos de entrenamiento son escasos, la varianza posterior vuelve a la anterior.

La probabilidad marginal y la afirmación que vale la pena probar

Los hiperparámetros del kernel. $\theta$ (escalas de longitud, variaciones, nivel de ruido) se aprenden maximizando la probabilidad marginal logarítmica:

$\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$

El primer término es un término de ajuste de datos (penaliza las predicciones alejadas de las observaciones). El segundo término es una penalización por complejidad (penaliza los modelos que son demasiado flexibles, es decir, donde la matriz del núcleo tiene un determinante grande). El tercer término es una constante de normalización.

Esta es la navaja automática de Occam y es lo más interesante que aporta un médico de cabecera a un canal de negociación. La versión sólida de la afirmación es que no se necesita un conjunto de validación separado para la regularización: la penalización por complejidad está dentro del objetivo, por lo que el modelo no puede comprar ajuste con flexibilidad de forma gratuita.

Esa afirmación merece escepticismo específicamente en este blog. Análisis de meseta muestra que una puntuación de validación de un solo punto es un mal criterio de selección y que la solidez depende de la forma del vecindario; PBO cuantifica con qué frecuencia el ganador dentro de la muestra pierde fuera de la muestra; Precios de Sharpe deflactado en el número de ensayos. Una probabilidad marginal sigue siendo un objetivo dentro de la muestra que se maximiza durante $\theta$ — el factor Occam penaliza la capacidad del modelo, no la selección entre muchos núcleos ajustados. Si se ajustan doce núcleos candidatos y se elige el que tiene la mejor probabilidad marginal, se vuelve a estar en territorio de pruebas múltiples y la lógica de la deflación se aplica sin cambios. La versión falsificable: ¿la selección de verosimilitud marginal produce una brecha dentro/fuera de la muestra más pequeña que la selección del conjunto de validación con los mismos datos y la misma familia de núcleos? Eso es mensurable y no se mide a continuación.

La superficie de probabilidad marginal también tiene óptimos locales. Múltiples reinicios aleatorios o una inicialización cuidadosa son importantes: inicializar la escala de longitud a la distancia mediana por pares de las entradas de entrenamiento y la varianza del ruido a la varianza de la muestra de los objetivos es un punto de partida razonable.

Costo computacional y escalabilidad

El cuello de botella se está invirtiendo $(K + \sigma_n^2 I)$, que cuesta $O(n^3)$ en el tiempo y $O(n^2)$ en la memoria. La pared cúbica ya se argumenta en este blog desde la otra dirección: método de búsqueda versus costo de evaluación descalifica por completo la optimización bayesiana basada en GP cuando el objetivo es barato, porque el sustituto cuesta más que las evaluaciones que ahorra. La aritmética es la misma aquí; La aplicación es diferente. Como modelo ajustado a los datos del mercado, el término cúbico no es una descalificación sino un presupuesto: establece un límite estricto en la ventana de capacitación.

Estrategias de escalabilidad para aplicaciones comerciales:

1. GP escasos (puntos inductores). Reemplace el $n \times n$ matriz con una $m \times m$ matriz donde $m \ll n$. Los puntos inductores $Z = \{z_1, \ldots, z_m\}$ Son pseudoentradas que resumen los datos de entrenamiento. La formulación SVGP (GP estocástica variacional) de Hensman et al. (2013) permite el entrenamiento en mini lotes con costo $O(nm^2)$ por iteración. GPyTorch admite esto de forma nativa.

2. Interpolación de kernel estructurada (SKI/KISS-GP). Explota la estructura de Kronecker y Toeplitz en la matriz del kernel cuando las entradas se encuentran en una cuadrícula. Reduce el costo de $O(n + g \log g)$ dónde $g$ es el tamaño de la cuadrícula. Ideal para series temporales muestreadas regularmente (por ejemplo, barras de 1 minuto).

3. ** Médicos de cabecera locales en ventanas deslizantes. ** Capacite a médicos de cabecera separados únicamente con datos recientes. Aquí es donde aparece la restricción específica de GP: los núcleos estándar son estacionarios ($k(x,x')$ depende sólo de $x - x'$) y los mercados no, por lo que la respuesta habitual es una ventana móvil, pero la duración de la ventana está limitada arriba por $O(n^3)$, no sólo por las estadísticas. Optimización de avance cubre ventanas ancladas versus ventanas móviles, longitudes de tren/prueba y frecuencia de reoptimización en general; para un médico de cabecera, el tamaño de la ventana es una decisión informática tanto como estadística, y una ventana anclada (en constante crecimiento) simplemente no está disponible más allá de unos pocos miles de puntos sin una aproximación. Ese replanteamiento es la consecuencia práctica: con los médicos de cabecera no puedes elegir "usar todo el historial".

GP para la predicción del retorno: un marco práctico

Diseño de funciones: por qué entre 5 y 20 funciones, no 500

La taxonomía de características generales para el ML de datos de mercado (desequilibrio de la cartera de pedidos, presión contable, VPIN, lambda de Kyle, características de volatilidad realizada, codificación de tiempo cíclico, señales de tasas de financiación y activos cruzados) ya está expuesta en modelado de spread con aprendizaje automático; usa esa lista.

Lo que es específico del médico de cabecera es el tamaño de la lista. Los métodos de kernel se degradan en dimensiones altas: las distancias se concentran y un kernel estacionario pierde discriminación. El rango práctico para un médico de cabecera financiero es de 5 a 20 entradas, no los cientos que un árbol impulsado por gradiente come felizmente. El mecanismo que hace que esto pueda sobrevivir es ARD (Determinación automática de relevancia): le da a cada dimensión de entrada su propia escala de longitud. $\ell_d$, y el entrenamiento de probabilidad marginal empuja $\ell_d \to \infty$ para dimensiones que no transportan señal, ya que una escala de longitud infinita significa que el núcleo ignora esa coordenada. La selección de características se convierte en un subproducto de la adaptación, y lo aprendido $\ell_d$ son directamente legibles como una clasificación de relevancia, lo que también lo hace falsificable: ajuste el núcleo compuesto en barras reales, imprima las escalas de longitud y vea si las características en las que cree son las que conserva el modelo.

Dimensionamiento de posiciones y abstención

El GP posterior da $\sigma_*(x)$ directamente, por lo que cae directamente en el filtro de tamaño de relación de borde y no comercio desarrollado en Predicción conforme para el tamaño de posiciones consciente del riesgo, con el ancho del intervalo $w_t = 2\kappa\sigma_*(x)$. La única diferencia es la procedencia: el GP produce el ancho a partir del propio modelo en lugar de un conjunto de calibración, por lo que varía según la ubicación del punto de prueba en relación con los datos de entrenamiento en lugar de con un cuantil global de residuos pasados.

Implementación con GPyTorch

GPyTorch es una biblioteca basada en PyTorch para inferencia de GP escalable. Aprovecha la aceleración de GPU, la diferenciación automática y las técnicas modernas de álgebra lineal (gradientes conjugados, descomposición de Lanczos) para escalar los GP más allá de lo ingenuo. $O(n^3)$ límite.

GP exacto básico para la predicción de retorno

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 es lo que permite las escalas de longitud por dimensión discutidas anteriormente. Después del entrenamiento, model.covar_matern.base_kernel.lengthscale es el vector a leer.

Bucle de entrenamiento

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

Predicción con incertidumbre

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(),
    }

El fast_pred_var() El administrador de contexto utiliza el algoritmo LOVE (Lanczos Variance Estimates) para calcular las variaciones predictivas en $O(n)$ tiempo en lugar de $O(n^2)$.

Canal de negociación de extremo a extremo

Nota sobre el generador de funciones a continuación: cada estadística móvil debe ser estrictamente retrospectiva y la estandarización de entrada debe ajustarse únicamente al segmento de entrenamiento. Ese segundo punto no es higiene genérica: es precisamente el canal de fuga de normalización analizado y medido en la taxonomía de sesgo anticipado, que informa la inflación de Sharpe que produce cada tipo de fuga. Un médico de cabecera estandariza sus entradas por construcción, por lo que esta es la fuga a la que está más expuesto.

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

Escalar a conjuntos de datos más grandes con médicos de cabecera dispersos

Cuando la ventana de entrenamiento excede unos pocos miles de puntos, la inferencia exacta del GP se vuelve lenta. Cambie a un GP disperso variante:

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

Con 128 puntos inductores, el costo por lote es $O(128^2 \times \text{batch\_size})$ – aproximadamente 4 millones de operaciones por lote. Esto maneja fácilmente conjuntos de datos de más de 100.000 observaciones en una sola GPU.

Aprendizaje profundo del kernel: GP se encuentra con las redes neuronales

Cuando el espacio de entrada es de alta dimensión o la relación entre características y retornos es altamente no lineal, un núcleo simple puede tener dificultades. El aprendizaje profundo del kernel (DKL) pasa las entradas a través de una red neuronal antes de aplicar el kernel GP:

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

dónde $g_\phi$ es una red neuronal con parámetros $\phi$ y $k_{\text{base}}$ es un kernel estándar (por ejemplo, Matern-3/2). La red aprende una representación de características donde el kernel GP es más efectivo. Todo el modelo (parámetros de red e hiperparámetros del kernel) se entrena de un extremo a otro maximizando la probabilidad marginal.

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 combina el aprendizaje de representación de redes neuronales con la cuantificación de incertidumbre de los médicos de cabecera. La capa GP garantiza que las predicciones alejadas de los datos de entrenamiento tengan una alta incertidumbre, algo que las redes neuronales estándar no logran proporcionar. Tenga en cuenta también que DKL reintroduce exactamente la flexibilidad que se suponía debía controlar la probabilidad marginal: el factor Occam penaliza al núcleo, pero la red que tiene delante tiene miles de parámetros libres y no tiene tal penalización.

Ludkovski y Risk (2025), Gaussian Process Models for Quantitative Finance, analizan DKL en un contexto financiero cuantitativo más amplio que incluye precios de opciones y optimización de carteras; esos resultados son suyos, de sus problemas, y no son evidencia sobre la predicción del retorno de las criptomonedas.

Diagnósticos y trampas

Calibración

Un médico de cabecera bien calibrado tiene intervalos predictivos que coinciden con la cobertura empírica. La verificación es la misma que se usa en predicción conforme, que también cubre la teoría de por qué la cobertura marginal no es cobertura condicional, una limitación que se aplica igualmente a los intervalos de GP:

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}"
        )

La superación esperada es 0,317/0,046/0,003 en 1/2/3 sigma. El número que importa es lo que esto imprime en un recorrido real sobre barras criptográficas, y esa tabla aún no está en este artículo. La expectativa anterior es que el médico de cabecera está demasiado confiado (una probabilidad gaussiana sobre rendimientos no estacionarios de cola gruesa debería ocultar mal en 3 sigma), pero "debería" no es una medida.

Errores restantes

1. Escala de entrada. Las escalas de longitud son relativas a la escala de entrada, por lo que una característica varía $[0, 10000]$ y uno a distancia $[0, 1]$ no puedo compartir un significado $\ell$; ARD es lo que rescata rangos heterogéneos y la estandarización es lo que hace que la inicialización de ARD sea sensata.

2. Sobreajuste de la probabilidad marginal. Con muchos hiperparámetros del núcleo (especialmente núcleos compuestos o de mezcla espectral), la probabilidad marginal aún puede sobreajustarse. Utilice valores previos: un valor previo logarítmico normal en escalas de longitud centrado en la distancia mediana por pares, un valor previo medio normal en variaciones.

3. Condicionamiento de matriz de covarianza. $(K + \sigma_n^2 I)$ puede volverse numéricamente singular cuando el ruido $\sigma_n^2$ es demasiado pequeño o cuando los puntos de entrenamiento están casi duplicados. GPyTorch agrega $10^{-6}$ inquietud hacia la diagonal; los datos financieros mal condicionados a menudo necesitan más.

4. Sesgo de anticipación. Cubierto anteriormente y, con detalle, en la taxonomía de sesgo de anticipación.

Cuándo utilizar médicos de cabecera frente a otros modelos

* Para las columnas XGBoost y redes neuronales, consulte las comparaciones publicadas en lugar de una nueva afirmación aquí: modelado extendido con aprendizaje automático tiene la tabla de aumento de gradiente versus aprendizaje profundo para datos tabulares financieros (umbrales de tamaño de datos, interpretabilidad, adaptación de régimen, latencia) y fusión temporal transformadores tiene TFT frente a LSTM frente a Vanilla Transformer.

Utilice médicos de cabecera cuando:

• Tiene conjuntos de datos pequeños y medianos (menos de ~10 000 observaciones por ventana de entrenamiento)
• Quiere una hipótesis estructural explícita e inspeccionable sobre la señal (tendencia + estacionalidad + ruido)
• La incertidumbre debe variar con la distancia desde los datos de entrenamiento, no solo con un cuantil residual global.

No utilice médicos de cabecera cuando:

• Necesita inferencias en submilisegundos sobre millones de observaciones
• La dimensionalidad de entrada supera ~50
• La señal vive en interacciones complejas de características de alto orden que los núcleos estacionarios no pueden representar (DKL ayuda, a costa de la propiedad de Occam)

Lo que este artículo aún no mide

Todo lo anterior es maquinaria modelo. Nada de esto es evidencia de que un médico de cabecera gane dinero con las criptomonedas, y el estándar de este blog es que el artículo incluya sus propios números. Los elementos abiertos, en el orden en que deben ejecutarse:

1. Escalas de longitud ARD en barras BTCUSDT reales de 1 m. Coloque el núcleo compuesto, imprima $\ell_d$ por característica. Esto prueba directamente la afirmación "ARD es una selección de funciones integrada" y produce una clasificación falsificable de relevancia de funciones.
2. La tabla de calibración de un análisis previo. Exceso esperado versus empírico en 1/2/3 sigma, expresado claramente, incluido el caso en el que el médico de cabecera resulta tener exceso de confianza.
3. La estrategia basada en la confianza, caminar hacia adelante, en las mismas cinco principales que la negativa honesta, desinflada para el recuento de pruebas. Si falla, se ubica en esa serie como otro resultado negativo.
4. El reloj de pared $O(n^3)$ curva para $n = 250 / 500 / 1000 / 2000 / 5000$, exacto frente a SVGP, por lo que la sección de escalabilidad se basa en un gráfico en lugar de una afirmación.

Conclusión

El argumento a favor de los procesos gaussianos en el trading no es que "te dan barras de error": la predicción conforme te da barras de error con menos supuestos distributivos y una garantía de cobertura que el GP no tiene. El caso es que un médico de cabecera te hace escribir tu hipótesis estructural como un núcleo, la ajusta a un objetivo que valora su propia complejidad y luego te dice cuál de tus características utilizó realmente a través de las escalas de longitud de ARD. Se trata de un modelo inusualmente legible.

Los costos son igualmente concretos: escalamiento cúbico que limita su ventana de entrenamiento, supuestos de estacionariedad que una ventana móvil sólo repara parcialmente, una probabilidad gaussiana de que los retornos de cola gruesa se violen y, una vez que se alcanza el aprendizaje profundo del núcleo, la pérdida silenciosa de la misma propiedad de Occam que motivó la probabilidad marginal en primer lugar. Si lo que queda es comercializable es una cuestión empírica, y las cuatro medidas enumeradas anteriormente son las que la responderían.

Referencias

• Rasmussen, CE y Williams, CKI (2006). Procesos Gaussianos para Aprendizaje Automático. Prensa del MIT.
• Wilson, A. y Adams, R. (2013). Núcleos del proceso gaussiano para el descubrimiento y la extrapolación de patrones. ICML.
• Hensman, J., Fusi, N. y Lawrence, ND (2013). Procesos Gaussianos para Big Data. UAI.
• Gatheral, J., Jaisson, T. y Rosenbaum, M. (2018). La volatilidad es dura. Finanzas cuantitativas, 18(6).
• Rizvi, S. A. A., Roberts, S. J., Osborne, M. A. y Nyikosa, F. (2017). Un enfoque novedoso para pronosticar la volatilidad financiera con envolventes de procesos gaussianos. arXiv:1705.00891.
• Ludkovski, M. y Risk, J. (2025). Modelos de Procesos Gaussianos para Finanzas Cuantitativas. Saltador.
• Gardner, J. R., Pleiss, G., Bindel, D., Weinberger, K. Q. y Wilson, A. G. (2018). GPyTorch: Inferencia de procesos gaussianos Blackbox Matrix-Matrix con aceleración de GPU. NeurIPS.