Processi gaussiani per la modellazione dei prezzi non parametrica

Parte della serie "Linee di base ML classiche".

Due cose rendono un processo gaussiano degno di un articolo separato su questo blog, e nessuno dei due è "ti dà incertezza".

Il primo è il design del kernel. L'intero pregiudizio induttivo di un medico di famiglia risiede in un'unica funzione $k(x, x')$, e quella funzione è qualcosa che scrivi deliberatamente: quanto è accidentato il percorso, se si ripete, se la ripetizione decade. Nient'altro nel toolkit standard consente di formulare un'ipotesi strutturale sulle dinamiche di mercato che espliciti e quindi adattarla. Il secondo è la probabilità marginale: un obiettivo formativo con una penalità di complessità derivata dal modello stesso, non da un insieme di resistenze. Ogni altro articolo nell'arco di overfitting su questo blog (analisi del plateau, PBO, deflated Sharpe) esiste perché la regolarizzazione del set di convalida è fragile durante la ricerca. Un medico di famiglia afferma di non averne bisogno. Questa affermazione è verificabile e testarla è più interessante di un altro tutorial sul dimensionamento dell’incertezza.

Sull’incertezza stessa: la varianza a posteriori GP è strutturale – cade dalla stessa inferenza che produce la media, piuttosto che essere avvolta successivamente attorno a un modello adattato. Questo è il vero contrasto con la previsione conforme, che spiega già in questo blog perché l'incertezza è l'input giusto per dimensionare la posizione e cosa fare con un intervallo una volta che ne hai uno. Questo articolo non riargomenta questo caso; va dietro al modello.

Quello che segue è il kernel e il meccanismo di inferenza, un'implementazione GPyTorch e, come affermato chiaramente alla fine, le misurazioni che questo articolo non ha ancora.

Cos'è un processo gaussiano?

I GP appaiono già su questo blog come surrogato dell'ottimizzazione bayesiana in Optuna vs. coordinate descends, con la stessa notazione e lo stesso avvertimento di bassa dimensionalità. Qui il GP è il modello stesso, adattato ai dati di mercato piuttosto che a una superficie di ricerca iperparametrica, quindi la trattazione va più in profondità.

Un processo gaussiano è un insieme di variabili casuali, un numero finito delle quali ha una distribuzione gaussiana congiunta. È una distribuzione sulle funzioni, non una distribuzione sui parametri.

Formalmente, una funzione $f: \mathcal{X} \to \mathbb{R}$ viene estratto da un GP se per qualsiasi insieme finito di input $\{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)$

Dove $m(x) = \mathbb{E}[f(x)]$ è la funzione media e $k(x, x') = \text{Cov}(f(x), f(x'))$ è la funzione di covarianza (kernel). Lo scriviamo in modo compatto come:

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

La funzione media codifica la convinzione precedente sul comportamento medio di $f$. Nel trading, in genere fissiamo $m(x) = 0$, codificando il presupposto che non vi siano precedenti distorsioni direzionali sui rendimenti. Tutta la struttura va nel kernel.

Perché non parametrico?

Un modello di regressione lineare con 5 caratteristiche ha 6 parametri. Una rete neurale con due strati nascosti di 64 unità ne ha migliaia. Un GP non ha un numero fisso di parametri: la complessità del modello cresce con i dati. Con 10 osservazioni, il GP definisce una gaussiana a 10 dimensioni. Con 10.000 osservazioni, definisce una gaussiana a 10.000 dimensioni.

Ciò non significa che i medici di base non abbiano iperparametri. La funzione del kernel ha iperparametri (scale di lunghezza, ampiezze, periodicità) che controllano le proprietà delle funzioni tratte dalla precedente. Ma la forma funzionale in sé non è mai fissa. Il GP può rappresentare qualsiasi funzione continua, dati abbastanza dati e il giusto kernel. Questo è ciò che significa "non parametrico": il modello non è vincolato a una famiglia parametrica come le funzioni lineari o i polinomi.

Per la modellazione finanziaria, questo è prezioso. I mercati cambiano. La relazione tra caratteristiche e rendimenti è non lineare, non stazionaria e dipendente dal regime. I modelli parametrici impongono una struttura che potrebbe non corrispondere alla realtà. I medici di base lasciano parlare i dati.

Funzioni del kernel: codifica della struttura del mercato

La funzione del kernel $k(x, x')$ è l'anima di un processo gaussiano. Definisce quali funzioni sono probabili a priori specificando la covarianza tra i valori della funzione in due punti di input qualsiasi. Kernel diversi codificano presupposti diversi su fluidità, periodicità e comportamento a lungo raggio.

Funzione a base radiale (RBF) / Esponenziale quadrato

Il kernel RBF è il punto di partenza più comune:

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

Dove $\sigma^2$ è la varianza del segnale (scala di uscita) e $\ell$ è la scala della lunghezza. Le funzioni tratte da un GP con un kernel RBF sono infinitamente differenziabili: molto fluide.

Interpretazione del trading: La scala della lunghezza $\ell$ controlla quanto distanti due punti dati possono essere ed essere comunque correlati. Una scala di breve durata significa che il modello reagisce a modelli locali; una scala di lunga durata significa che cattura tendenze generali. La varianza del segnale $\sigma^2$ controlla l'ampiezza della funzione: quanto possono essere grandi i rendimenti previsti.

Problema finanziario: La fluidità infinita non è realistica. I rendimenti finanziari presentano salti, cambiamenti di regime e discontinuità. Il kernel RBF può attenuare eccessivamente queste caratteristiche, producendo previsioni troppo prudenti in prossimità di rotture strutturali.

Kernel materno

La classe Matern generalizza l'RBF introducendo un parametro di fluidità $\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)$

Dove $K_{\nu}$ è la funzione di Bessel modificata del secondo tipo. COME $\nu \to \infty$, il nucleo Matern converge al RBF. Scelte comuni:

• $\nu = 1/2$: Equivalente al processo Ornstein-Uhlenbeck. Le funzioni sono continue ma non differenziabili — grezze, come il moto browniano.
• $\nu = 3/2$: Le funzioni sono differenziabili una volta. Un buon equilibrio tra morbidezza e flessibilità.
• $\nu = 5/2$: Le funzioni sono due volte differenziabili. Più liscio di $3/2$ ma meno rigido di RBF.

Interpretazione commerciale: Il Matern-$3/2$ kernel è probabilmente il miglior valore predefinito per le serie temporali finanziarie. Permette il tipo di ruvidità che i percorsi dei prezzi reali mostrano senza essere così frastagliati $\nu = 1/2$. Ciò è in linea con la letteratura "la volatilità è approssimativa" (Gatheral, Jaisson e Rosenbaum, 2018), che mostra empiricamente che i percorsi di volatilità hanno esponenti di Hurst intorno a $H \approx 0.1$, molto più approssimativo del moto browniano.

La principale prova pubblicata che i kernel Matern battono i modelli classici di volatilità è Rizvi et al. (2017), che riportano un MSE migliore di circa il 20% rispetto a una passeggiata casuale e del 50% migliore di GARCH - sui dati giornalieri sulle coppie di valute del 2017, non crittografici, e non riprodotti qui. Consideralo come motivazione per provare il kernel, non come un punto di riferimento. L'adattamento GARCH(1,1) di questo blog sui dati giornalieri BTC reali, con la diagnostica Ljung-Box e ARCH-LM, è in Previsione della volatilità GARCH per criptovalute; un testa a testa contro un GP di Matern sullo stesso campione sarebbe il confronto onesto, e non è stato effettuato.

Kernel periodico

I mercati finanziari hanno schemi ciclici: curve di volume intraday, effetti giornalieri della settimana, flussi di ribilanciamento mensili, stagioni degli utili trimestrali. Il kernel periodico cattura questi:

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

Dove $p$ è il periodo. Le funzioni estratte da questo kernel si ripetono con il punto $p$, modulato dalla scala di lunghezza $\ell$ che controlla la velocità con cui la correlazione decade all'interno di un periodo.

Interpretazione commerciale: Imposta $p = 24$ (ore) per acquisire modelli intraday, o $p = 5$ (giorni di Borsa) per la stagionalità settimanale. A differenza delle caratteristiche di Fourier, il nucleo periodico non presuppone un numero fisso di armoniche: il GP apprende la forma del ciclo dai dati.

Combinazione di noccioli: composizione additiva e moltiplicativa

Il vero potere dei kernel GP risiede nella composizione. Se $k_1$ E $k_2$ sono kernel validi, così come:

• Somma: $k_1 + k_2$ — la funzione è la somma di componenti indipendenti (scomposizione additiva)
• Prodotto: $k_1 \times k_2$ — interazioni tra componenti (ad esempio, comportamento localmente periodico)

Un utile nucleo composito per i rendimenti finanziari:

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

Questo decompone il segnale in:

1. Una componente di tendenza non uniforme e aperiodica (Matern-3/2)
2. Una componente periodica la cui ampiezza decade nel tempo (Periodic $\times$ RF)

Il prodotto $k_{\text{Periodic}} \cdot k_{\text{RBF}}$ crea un nucleo localmente periodico: lo schema si ripete, ma le ripetizioni distanti hanno meno influenza di quelle vicine. Ciò è esattamente vero per la stagionalità finanziaria, che varia nel tempo con l’evoluzione della microstruttura del mercato.

Noccioli di miscela spettrale

Per la massima flessibilità, il kernel della miscela spettrale (SM) (Wilson & Adams, 2013) parametrizza la densità spettrale del kernel come una miscela di gaussiane:

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

Dove $w_q$ sono pesi misti, $v_q$ sono varianze spettrali e $\mu_q$ sono mezzi spettrali (frequenze). Per il teorema di Bochner, qualsiasi nucleo stazionario può essere rappresentato in questo modo. Il kernel SM è in grado di scoprire simultaneamente componenti periodici, trend a lungo termine e correlazioni a breve termine, il tutto a partire dai dati.

Interpretazione del trading: Il kernel SM è utile quando non si sanno quali modelli esistono nei dati. Può identificare periodicità nascoste nelle serie di rendimenti (ad esempio, cicli sottili di 4 ore nei mercati delle criptovalute guidati dal ribilanciamento automatizzato). Lo svantaggio è rappresentato dal maggior numero di iperparametri e dal rischio di sovradimensionamento con set di dati di piccole dimensioni.

Inferenza posteriore: da prima della previsione

Dati i dati di allenamento $\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n$ Dove $y_i = f(x_i) + \epsilon_i$ E $\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$, il GP posteriore nei punti di test $X_*$ ha una soluzione in forma chiusa. Questo è il vantaggio computazionale chiave dei GP rispetto alla maggior parte dei modelli bayesiani.

Le equazioni posteriori

Lasciamo $K = k(X, X)$ essere il $n \times n$ matrice di covarianza dell'allenamento, $K_* = k(X_*, X)$ essere il $m \times n$ matrice di covarianza incrociata e $K_{**} = k(X_*, X_*)$ essere il $m \times m$ matrice di covarianza del test. Il posteriore è:

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

Dove:

$\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 a posteriori $\bar{f}_*$ è una combinazione lineare degli obiettivi di allenamento, ponderata in base alla somiglianza del kernel tra i punti di test e quelli di allenamento. La covarianza posteriore $\text{Cov}(f_*)$ parte dalla covarianza precedente $K_{**}$ e sottrae le informazioni ottenute dai dati di addestramento. Laddove i dati di addestramento sono densi, la varianza a posteriori è piccola. Laddove i dati di addestramento sono scarsi, la varianza a posteriori ritorna a quella a priori.

La verosimiglianza marginale e l'affermazione che vale la pena verificare

Gli iperparametri del kernel $\theta$ (scale di lunghezza, varianze, livello di rumore) vengono appresi massimizzando la verosimiglianza marginale logaritmica:

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

Il primo termine è un termine di data fit (penalizza le previsioni lontane dalle osservazioni). Il secondo termine è una penalità di complessità (penalizza i modelli troppo flessibili, ovvero in cui la matrice del kernel ha un determinante grande). Il terzo termine è una costante di normalizzazione.

Questo è il rasoio di Occam automatico, ed è la cosa più interessante che un GP porta in una pipeline di trading. La versione forte dell'affermazione è che non è necessario alcun set di convalida separato per la regolarizzazione: la penalità di complessità è all'interno dell'obiettivo, quindi il modello non può adattarsi gratuitamente alla flessibilità.

Tale affermazione merita scetticismo in particolare su questo blog. L'analisi del plateau mostra che un punteggio di convalida a punto singolo è un cattivo criterio di selezione e che la robustezza risiede nella forma del quartiere; PBO quantifica la frequenza con cui il vincitore nel campione perde fuori dal campione; deflated Sharpe prezzi nel numero di prove. Una probabilità marginale è ancora un obiettivo interno al campione che viene massimizzato $\theta$ — il fattore Occam penalizza la capacità del modello, non la selezione su molti chicchi adattati. Se si adattano dodici nuclei candidati e si sceglie quello con la migliore probabilità marginale, si torna nel territorio dei test multipli e la logica di deflazione si applica invariata. La versione falsificabile: la selezione della verosimiglianza marginale produce un divario all’interno del campione/fuori campione più piccolo rispetto alla selezione del set di validazione sugli stessi dati e sulla stessa famiglia di kernel? Questo è misurabile e non viene misurato di seguito.

La superficie di verosimiglianza marginale ha anche ottimi locali. Riavvii casuali multipli o un'inizializzazione attenta sono importanti: inizializzare la scala di lunghezza sulla distanza mediana a coppie degli input di addestramento e la varianza del rumore sulla varianza del campione dei target è un punto di partenza ragionevole.

Costo computazionale e scalabilità

Il collo di bottiglia si sta invertendo $(K + \sigma_n^2 I)$, che costa $O(n^3)$ nel tempo e $O(n^2)$ nella memoria. Il muro cubico è già discusso su questo blog dalla direzione opposta: metodo di ricerca vs. costo di valutazione squalifica completamente l'ottimizzazione bayesiana basata su GP quando l'obiettivo è economico, perché il surrogato costa più delle valutazioni che salva. Qui l'aritmetica è la stessa; l'applicazione è diversa. In quanto modello adattato ai dati di mercato, il termine cubico non è una squalifica ma un budget: fissa un tetto rigido alla finestra di formazione.

Strategie di scalabilità per le applicazioni di trading:

1. GP sparsi (induzione di punti). Sostituisci il $n \times n$ matrice con an $m \times m$ matrice dove $m \ll n$. I punti induttori $Z = \{z_1, \ldots, z_m\}$ sono pseudo-input che riassumono i dati di addestramento. La formulazione SVGP (Stochastic Variational GP) di Hensman et al. (2013) consente l'addestramento in mini-batch a pagamento $O(nm^2)$ per iterazione. GPyTorch lo supporta in modo nativo.

2. Interpolazione del kernel strutturato (SKI/KISS-GP). Sfrutta la struttura di Kronecker e Toeplitz nella matrice del kernel quando gli input si trovano su una griglia. Riduce i costi a $O(n + g \log g)$ Dove $g$ è la dimensione della griglia. Ideale per serie temporali campionate regolarmente (ad esempio, barre da 1 minuto).

3. GP locali su finestre scorrevoli. Addestra GP separati solo sui dati recenti. È qui che entra in gioco il vincolo specifico del GP: i kernel standard sono stazionari ($k(x,x')$ dipende solo da $x - x'$) e i mercati no, quindi la risposta abituale è una finestra mobile, ma la lunghezza della finestra è limitata sopra $O(n^3)$, non solo dalle statistiche. Ottimizzazione walk-forward copre le finestre ancorate e quelle rotanti, la durata dei treni/test e la frequenza di riottimizzazione in generale; per un medico di famiglia, il dimensionamento della finestra è una decisione di calcolo tanto quanto statistica, e una finestra ancorata (in continua crescita) semplicemente non è disponibile oltre poche migliaia di punti senza approssimazione. Questa ristrutturazione è la conseguenza pratica: con i medici di base non puoi scegliere di "utilizzare tutta la cronologia".

GP per la previsione del rendimento: un quadro pratico

Design delle funzionalità: perché 5-20 funzionalità e non 500

La tassonomia generale delle caratteristiche generali per il machine learning dei dati di mercato (sbilanciamento del book degli ordini, pressione del book, VPIN, lambda di Kyle, caratteristiche della volatilità realizzata, codifica temporale ciclica, segnali tra asset e tasso di finanziamento) è già illustrata in spread modeling with machine learning; usa quella lista.

Ciò che è specifico del GP è la dimensione dell'elenco. I metodi del kernel si degradano in dimensioni elevate: le distanze si concentrano e un kernel stazionario perde la discriminazione. Il range pratico per un GP finanziario è di 5-20 input, non le centinaia che un albero potenziato dal gradiente mangia felicemente. Il meccanismo che rende tutto ciò sostenibile è ARD (Automatic Relevance Determination): assegna a ciascuna dimensione di input la propria scala di lunghezza $\ell_d$e spinte alla formazione a probabilità marginale $\ell_d \to \infty$ per dimensioni che non portano alcun segnale, poiché una scala di lunghezza infinita significa che il kernel ignora quella coordinata. La selezione delle caratteristiche diventa un sottoprodotto dell'adattamento e dell'apprendimento $\ell_d$ sono direttamente leggibili come classifica di pertinenza, che è anche ciò che lo rende falsificabile: adatta il kernel composito su barre reali, stampa le scale di lunghezza e vedi se le caratteristiche in cui credi sono quelle conservate dal modello.

Ridimensionamento della posizione e astensione

Il GP posteriore dà $\sigma_*(x)$ direttamente, quindi cade direttamente nel dimensionamento del rapporto edge e nel filtro no-trade sviluppato in Predizione conforme per dimensionamento della posizione consapevole del rischio, con la larghezza dell'intervallo $w_t = 2\kappa\sigma_*(x)$. L'unica differenza è la provenienza: il GP produce la larghezza dal modello stesso piuttosto che da un set di calibrazione, quindi varia in base alla posizione del punto di test rispetto ai dati di addestramento piuttosto che con un quantile globale di residui passati.

Implementazione con GPyTorch

GPyTorch è una libreria basata su PyTorch per l'inferenza GP scalabile. Sfrutta l'accelerazione GPU, la differenziazione automatica e le moderne tecniche di algebra lineare (gradienti coniugati, decomposizione di Lanczos) per scalare i GP oltre l'ingenuo $O(n^3)$ limite.

GP esatto di base per la previsione del rendimento

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 è ciò che consente le scale di lunghezza per dimensione discusse sopra. Dopo l'allenamento, model.covar_matern.base_kernel.lengthscale è il vettore da leggere.

Ciclo di allenamento

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

Previsione con incertezza

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

IL fast_pred_var() il gestore del contesto utilizza l'algoritmo LOVE (Lanczos Variance Estimates) per calcolare le varianze predittive $O(n)$ tempo invece di $O(n^2)$.

Pipeline di trading end-to-end

Nota sul generatore di funzionalità di seguito: ogni statistica mobile deve essere rigorosamente guardando al passato e la standardizzazione dell'input deve essere adattata solo alla sezione di addestramento. Questo secondo punto non è un problema di igiene generica: è proprio il canale di normalizzazione delle perdite sezionato e misurato nella tassonomia dei pregiudizi look-ahead, che riporta l'inflazione di Sharpe prodotta da ciascun tipo di perdita. Un GP standardizza i suoi input per costruzione, quindi questa è la perdita a cui è maggiormente esposto.

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

Adattamento a set di dati più grandi con GP sparsi

Quando la finestra di allenamento supera alcune migliaia di punti, l'esatta inferenza GP diventa lenta. Passa a un GP sparse variazionale:

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 punti di induzione, il costo per lotto è $O(128^2 \times \text{batch\_size})$ — circa 4 milioni di operazioni per lotto. Questo gestisce facilmente set di dati di oltre 100.000 osservazioni su una singola GPU.

Apprendimento profondo del kernel: GP incontra le reti neurali

Quando lo spazio di input è altamente dimensionale o la relazione tra caratteristiche e rendimenti è altamente non lineare, un kernel semplice potrebbe avere difficoltà. Il deep kernel learning (DKL) passa gli input attraverso una rete neurale prima di applicare il kernel GP:

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

Dove $g_\phi$ è una rete neurale con parametri $\phi$ E $k_{\text{base}}$ è un kernel standard (ad esempio, Matern-3/2). La rete apprende una rappresentazione delle funzionalità in cui il kernel GP è più efficace. L'intero modello (parametri di rete e iperparametri del kernel) viene addestrato end-to-end massimizzando la probabilità marginale.

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 l'apprendimento della rappresentazione delle reti neurali con la quantificazione dell'incertezza dei GP. Il livello GP garantisce che le previsioni lontane dai dati di addestramento abbiano un’elevata incertezza, qualcosa che notoriamente le reti neurali standard non riescono a fornire. Si noti inoltre che DKL reintroduce esattamente la flessibilità che la probabilità marginale avrebbe dovuto controllare: il fattore Occam penalizza il kernel, ma la rete di fronte ad esso ha migliaia di parametri liberi e nessuna penalità del genere.

Ludkovski e Risk (2025), Gaussian Process Models for Quantitative Finance, esaminano DKL in un contesto di finanza quantitativa più ampio, comprendente la determinazione del prezzo delle opzioni e l'ottimizzazione del portafoglio; quei risultati appartengono a loro, ai loro problemi, e non costituiscono una prova della previsione del rendimento delle criptovalute.

Diagnostica e insidie

Calibrazione

Un medico di famiglia ben calibrato ha intervalli predittivi che corrispondono alla copertura empirica. Il controllo è lo stesso utilizzato in previsione conforme — che copre anche la teoria del perché la copertura marginale non è una copertura condizionale, una limitazione che si applica altrettanto agli intervalli 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}"
        )

Il superamento previsto è 0,317/0,046/0,003 a 1/2/3 sigma. Il numero che conta è ciò che viene stampato su una vera corsa in avanti sulle barre crittografiche, e quella tabella non è ancora in questo articolo. L’aspettativa preliminare è che il medico di famiglia sia troppo sicuro di sé – una probabilità gaussiana su rendimenti non stazionari e a coda grassa dovrebbe essere mal coperta a 3 sigma – ma “dovrebbe” non è una misura.

Insidie rimanenti

1. Ridimensionamento dell'input. Le scale di lunghezza sono relative alla scala di input, quindi una funzionalità variava $[0, 10000]$ e uno a distanza $[0, 1]$ non può condividere un significato $\ell$; ARD è ciò che salva intervalli eterogenei e la standardizzazione è ciò che rende sensata l'inizializzazione di ARD.

2. Adattamento eccessivo della probabilità marginale. Con molti iperparametri del kernel (in particolare kernel compositi o di miscele spettrali), la probabilità marginale può ancora adattarsi eccessivamente. Utilizzare i valori a priori: un valore a priori log-normale sulle scale di lunghezza centrata sulla distanza mediana a coppie, un valore a priori semi-normale sulle varianze.

3. Condizionamento della matrice di covarianza. $(K + \sigma_n^2 I)$ può diventare numericamente singolare quando il rumore $\sigma_n^2$ è troppo piccolo o quando i punti di allenamento sono quasi duplicati. GPyTorch aggiunge $10^{-6}$ jitter sulla diagonale; i dati finanziari mal condizionati spesso hanno bisogno di qualcosa di più.

4. Pregiudizio di previsione. Trattato sopra e, in modo dettagliato, nella tassonomia del pregiudizio di previsione.

Quando utilizzare i medici di famiglia rispetto ad altri modelli

* Per le colonne XGBoost e rete neurale, rimandare ai confronti pubblicati piuttosto che a una nuova affermazione qui: spread modeling with machine learning ha la tabella gradient-boosting-vs-deep-learning per i dati tabulari finanziari (soglie delle dimensioni dei dati, interpretabilità, adattamento del regime, latenza) e fusione temporale trasformatori ha TFT, LSTM e Vanilla Transformer.

Utilizzare i medici di famiglia quando:

• Hai set di dati di piccole e medie dimensioni (meno di circa 10.000 osservazioni per finestra di addestramento)
• Vuoi un'ipotesi strutturale esplicita e ispezionabile sul segnale (tendenza + stagionalità + rumore)
• L'incertezza deve variare con la distanza dai dati di addestramento, non solo con un quantile residuo globale

Non utilizzare i medici di famiglia quando:

• Hai bisogno di un'inferenza inferiore al millisecondo su milioni di osservazioni
• La dimensionalità dell'input supera ~50
• Il segnale vive in complesse interazioni di caratteristiche di ordine elevato che i kernel stazionari non possono rappresentare (DKL aiuta, a scapito della proprietà Occam)

Ciò che questo articolo non misura ancora

Tutto quanto sopra è un modello di macchinario. Niente di tutto ciò è una prova che un medico di famiglia guadagni con le criptovalute, e lo standard di questo blog è che l'articolo riporti i propri numeri. Gli elementi aperti, nell'ordine in cui devono essere eseguiti:

1. La lunghezza dell'ARD è scalabile su barre reali da 1 milione di BTCUSDT. Adatta il kernel composito, stampa $\ell_d$ per funzione. Ciò verifica direttamente l'affermazione "ARD è una selezione di funzionalità integrata" e produce una classificazione falsificabile della pertinenza delle funzionalità.
2. La tabella di calibrazione da una corsa walk-forward. Superamento previsto rispetto a quello empirico a 1/2/3 sigma, dichiarato chiaramente, incluso il caso in cui il medico di famiglia risulta essere troppo sicuro di sé.
3. La strategia basata sulla fiducia, walk-forward, sugli stessi cinque major di the onesto negativo, sgonfiata per il conteggio delle prove. Se fallisce, viene inserita in quella serie come un altro risultato negativo.
4. L'orologio da parete $O(n^3)$ curva per $n = 250 / 500 / 1000 / 2000 / 5000$, esatto rispetto a SVGP, quindi la sezione sulla scalabilità si basa su un grafico anziché su un'asserzione.

Conclusione

Il caso dei processi gaussiani nel trading non è "ti danno barre di errore" - la previsione conforme ti dà barre di errore con meno presupposti distribuzionali e una garanzia di copertura che il GP non ha. Il caso è che un GP ti fa scrivere la tua ipotesi strutturale come un nucleo, la adatta a un obiettivo che valuta nella sua stessa complessità, e poi ti dice quale delle tue caratteristiche ha effettivamente utilizzato tramite le scale di lunghezza ARD. Questo è un modello insolitamente leggibile.

I costi sono ugualmente concreti: ridimensionamento cubico che limita la finestra di addestramento, ipotesi di stazionarietà che una finestra mobile ripara solo in parte, una probabilità gaussiana che i rendimenti a coda grassa violeranno e – una volta raggiunto l’apprendimento profondo del kernel – la silenziosa perdita della stessa proprietà di Occam che ha motivato in primo luogo la probabilità marginale. Se ciò che resta sia commerciabile è una questione empirica, e le quattro misurazioni sopra elencate sono ciò che potrebbe rispondere.

Riferimenti

• Rasmussen, C.E., & Williams, C.K.I. (2006). Processi gaussiani per l'apprendimento automatico. Stampa del MIT.
• Wilson, A. e Adams, R. (2013). Kernel di processo gaussiano per la scoperta e l'estrapolazione di pattern. ICML.
• Hensman, J., Fusi, N. e Lawrence, N. D. (2013). Processi gaussiani per i Big Data. UAI.
• Gatheral, J., Jaisson, T. e Rosenbaum, M. (2018). La volatilità è approssimativa. Finanza quantitativa, 18(6).
• Rizvi, S. A. A., Roberts, S. J., Osborne, M. A. e Nyikosa, F. (2017). Un nuovo approccio alla previsione della volatilità finanziaria con inviluppi di processi gaussiani. arXiv:1705.00891.
• Ludkovski, M., & Risk, J. (2025). Modelli di processo gaussiani per la finanza quantitativa. Springer.
• Gardner, J. R., Pleiss, G., Bindel, D., Weinberger, K. Q., & Wilson, A. G. (2018). GPyTorch: inferenza del processo gaussiano Matrix-Matrix Blackbox con accelerazione GPU. NeurIPS.