Gaussiaanse processen voor niet-parametrische prijsmodellering

Onderdeel van de serie "Classic ML baselines".

Er zijn twee dingen die een Gaussiaans proces een apart artikel op deze blog waard maken, en geen van beide is 'het geeft je onzekerheid'.

De eerste is kernelontwerp. De hele inductieve vooringenomenheid van een huisarts leeft in één functie $k(x, x')$, en die functie schrijf je doelbewust op: hoe ruw het pad is, of het zich herhaalt, of de herhaling vervalt. Met niets anders in de standaardtoolkit kunt u een structurele hypothese over de marktdynamiek formuleren die er expliciet bij past. De tweede is de marginale waarschijnlijkheid: een trainingsdoel met een complexiteitsboete die is afgeleid van het model zelf, en niet van een vastgehouden set. Elk ander artikel in de overfitting-boog op deze blog (plateau-analyse, PBO, deflated Sharpe) bestaat omdat de regularisatie van de validatieset kwetsbaar is bij zoeken. Een huisarts beweert er geen nodig te hebben. Die bewering is toetsbaar, en het testen ervan is interessanter dan weer een tutorial over de omvang van de onzekerheid.

Over de onzekerheid zelf: de posterieure variantie van de GP is structureel – zij valt buiten dezelfde gevolgtrekking die het gemiddelde oplevert, in plaats van achteraf om een ​​passend model te worden gewikkeld. Dat is het echte contrast met conformal forecast, waarin op deze blog al wordt besproken waarom onzekerheid de juiste input is voor positiebepaling en wat je met een interval moet doen als je er eenmaal een hebt. In dit artikel wordt deze zaak niet opnieuw beargumenteerd; het gaat achter het model aan.

Wat volgt zijn de kernel- en gevolgtrekkingsmachines, een GPyTorch-implementatie en – duidelijk vermeld aan het einde – de metingen die dit artikel nog niet heeft.

Wat is een Gaussiaans proces?

Huisartsen verschijnen al op deze blog als een surrogaat voor Bayesiaanse optimalisatie in Optuna vs. coördinatenafdaling, met dezelfde notatie en hetzelfde voorbehoud met betrekking tot de lage dimensionaliteit. Hier is de huisarts het model zelf, aangepast aan marktgegevens in plaats van aan een zoekoppervlak met hyperparameters, dus de behandeling gaat dieper.

Een Gaussiaans proces is een verzameling willekeurige variabelen, waarvan elk eindig aantal een gezamenlijke Gaussiaanse verdeling heeft. Het is een verdeling over functies, niet een verdeling over parameters.

Formeel een functie $f: \mathcal{X} \to \mathbb{R}$ wordt ontleend aan een GP indien voor een eindige reeks invoergegevens $\{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)$

waar $m(x) = \mathbb{E}[f(x)]$ is de gemiddelde functie en $k(x, x') = \text{Cov}(f(x), f(x'))$ is de covariantie (kernel) functie. We schrijven dit compact als:

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

De gemiddelde functie codeert voor eerdere overtuigingen over het gemiddelde gedrag van $f$. Bij het handelen stellen we doorgaans in $m(x) = 0$, wat de veronderstelling codeert dat we vooraf geen directionele bias op het gebied van rendement hebben. Alle structuur gaat in de kernel.

Waarom niet-parametrisch?

Een lineair regressiemodel met 5 kenmerken heeft 6 parameters. Een neuraal netwerk met twee verborgen lagen van 64 eenheden heeft duizenden. Een huisarts heeft geen vast aantal parameters; de complexiteit van het model groeit met de gegevens. Met 10 waarnemingen definieert de huisarts een 10-dimensionale Gaussiaans. Met 10.000 waarnemingen definieert het een 10.000-dimensionale Gaussiaans.

Dit betekent niet dat huisartsen geen hyperparameters hebben. De kernelfunctie heeft hyperparameters (lengteschalen, amplitudes, periodiciteiten) die de eigenschappen besturen van functies die uit de vorige zijn ontleend. Maar de functionele vorm zelf staat nooit vast. De GP kan elke continue functie vertegenwoordigen, gegeven voldoende gegevens en de juiste kernel. Dit is wat 'niet-parametrisch' betekent: het model is niet beperkt tot een parametrische familie zoals lineaire functies of polynomen.

Voor financiële modellering is dit waardevol. Markten veranderen. De relatie tussen kenmerken en rendementen is niet-lineair, niet-stationair en afhankelijk van het regime. Parametrische modellen leggen structuur op die mogelijk niet overeenkomt met de werkelijkheid. Huisartsen laten de data spreken.

Kernelfuncties: coderen van marktstructuur

De kernelfunctie $k(x, x')$ is de ziel van een Gaussiaans proces. Het definieert welke functies a priori waarschijnlijk zijn door de covariantie tussen functiewaarden op twee willekeurige invoerpunten te specificeren. Verschillende kernels coderen verschillende aannames over gladheid, periodiciteit en gedrag op lange termijn.

Radiale basisfunctie (RBF) / kwadratisch exponentieel

De RBF-kernel is het meest voorkomende startpunt:

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

waar $\sigma^2$ is de signaalvariantie (uitvoerschaal) en $\ell$ is de lengteschaal. Functies afkomstig van een GP met een RBF-kernel zijn oneindig differentieerbaar - zeer soepel.

Handelsinterpretatie: De lengteschaal $\ell$ bepaalt hoe ver twee gegevenspunten uit elkaar kunnen liggen en nog steeds gecorreleerd kunnen zijn. Een korte lengteschaal betekent dat het model reageert op lokale patronen; een lange lengteschaal betekent dat het brede trends vastlegt. De signaalvariantie $\sigma^2$ regelt de amplitude van de functie: hoe groot het voorspelde rendement kan zijn.

Financieel probleem: Oneindige gladheid is onrealistisch. Financiële rendementen kennen sprongen, regimeveranderingen en discontinuïteiten. De RBF-kernel kan deze kenmerken verzachten en voorspellingen produceren die te conservatief zijn in de buurt van structurele breuken.

Moederkernel

De Matern-klasse generaliseert de RBF door een gladheidsparameter te introduceren $\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)$

waar $K_{\nu}$ is de gemodificeerde Bessel-functie van de tweede soort. Als $\nu \to \infty$, convergeert de Matern-kernel naar de RBF. Veel voorkomende keuzes:

• $\nu = 1/2$: Equivalent aan het Ornstein-Uhlenbeck-proces. Functies zijn continu maar niet differentieerbaar – ruw, zoals de Brownse beweging.
• $\nu = 3/2$: Functies zijn ooit differentieerbaar. Een goede balans tussen soepelheid en flexibiliteit.
• $\nu = 5/2$: Functies zijn tweemaal differentieerbaar. Gladder dan $3/2$ maar minder rigide dan RBF.

Handelsinterpretatie: De Moeder-$3/2$ kernel is misschien wel de beste standaard voor financiële tijdreeksen. Het maakt het soort ruwheid mogelijk dat echte prijspaden vertonen, zonder zo grillig te zijn als $\nu = 1/2$. Dit komt overeen met de ‘volatiliteit is ruw’-literatuur (Gatheral, Jaisson, & Rosenbaum, 2018), die empirisch aantoont dat volatiliteitspaden Hurst-exponenten om zich heen hebben. $H \approx 0.1$, veel ruwer dan de Brownse beweging.

Het belangrijkste gepubliceerde bewijs dat Matern-kernels de klassieke volatiliteitsmodellen verslaan, is Rizvi et al. (2017), die grofweg 20% ​​betere MSE rapporteren dan een willekeurige wandeling en 50% beter dan GARCH - op basis van dagelijkse valutapaargegevens uit 2017, niet crypto, en hier niet weergegeven. Behandel het als motivatie om de kernel te proberen, niet als een benchmark. De eigen GARCH(1,1) van deze blog past bij echte dagelijkse BTC-gegevens, met Ljung-Box- en ARCH-LM-diagnostiek, in GARCH-volatiliteitsvoorspelling voor crypto; een onderlinge strijd tegen een moederarts uit hetzelfde monster zou een eerlijke vergelijking zijn, en die is nog niet uitgevoerd.

Periodieke kernel

Financiële markten kennen cyclische patronen: intraday-volumecurves, dag-van-week-effecten, maandelijkse herbalanceringsstromen, kwartaalwinstseizoenen. De periodieke kernel legt deze vast:

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

waar $p$ is de periode. Functies uit deze kernel worden herhaald met een punt $p$, gemoduleerd door de lengteschaal $\ell$ die bepaalt hoe snel de correlatie binnen een periode vervalt.

Handelsinterpretatie: Set $p = 24$ (uren) om intraday-patronen vast te leggen, of $p = 5$ (handelsdagen) voor wekelijkse seizoensinvloeden. In tegenstelling tot Fourier-kenmerken gaat de periodieke kernel niet uit van een vast aantal harmonischen; de huisarts leert de vorm van de cyclus uit gegevens.

Kernels combineren: additieve en multiplicatieve samenstelling

De echte kracht van GP-kernels ligt in de samenstelling. Als $k_1$ En $k_2$ zijn geldige kernels, en dat geldt ook voor:

• Som: $k_1 + k_2$ — de functie is de som van onafhankelijke componenten (additieve ontleding)
• Product: $k_1 \times k_2$ — interacties tussen componenten (bijvoorbeeld lokaal periodiek gedrag)

Een nuttige samengestelde kernel voor financieel rendement:

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

Dit ontleedt het signaal in:

1. Een niet-vloeiende, aperiodieke trendcomponent (Matern-3/2)
2. Een periodieke component waarvan de amplitude in de loop van de tijd afneemt (Periodiek $\times$ RBF)

Het product $k_{\text{Periodic}} \cdot k_{\text{RBF}}$ creëert een lokaal periodieke kern: het patroon herhaalt zich, maar herhalingen op afstand hebben minder invloed dan herhalingen in de buurt. Dit is precies goed voor de financiële seizoensgebondenheid, die in de loop van de tijd verandert naarmate de microstructuur van de markt evolueert.

Spectrale mengselkernels

Voor maximale flexibiliteit parameteriseert de spectrale mengselkernel (SM) (Wilson & Adams, 2013) de spectrale dichtheid van de kernel als een mengsel van Gaussianen:

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

waar $w_q$ zijn mengselgewichten, $v_q$ zijn spectrale varianties, en $\mu_q$ zijn spectrale gemiddelden (frequenties). Volgens de stelling van Bochner kan elke stationaire kern op deze manier worden weergegeven. De SM-kernel kan tegelijkertijd periodieke componenten, langetermijntrends en kortetermijncorrelaties ontdekken – allemaal op basis van gegevens.

Handelsinterpretatie: De SM-kernel is handig als u niet weet welke patronen er in de gegevens voorkomen. Het kan verborgen periodiciteiten in retourreeksen identificeren (bijvoorbeeld subtiele cycli van vier uur in cryptomarkten, aangedreven door geautomatiseerde herbalancering). Het nadeel is meer hyperparameters en het risico van overfitting met kleine datasets.

Posterieure gevolgtrekking: van vóór naar voorspelling

Gezien trainingsgegevens $\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n$ waar $y_i = f(x_i) + \epsilon_i$ En $\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$, de huisarts posterieur op de testpunten $X_*$ heeft een gesloten oplossing. Dit is het belangrijkste rekenvoordeel van huisartsen ten opzichte van de meeste Bayesiaanse modellen.

De posterieure vergelijkingen

Laat $K = k(X, X)$ wees de $n \times n$ training covariantiematrix, $K_* = k(X_*, X)$ wees de $m \times n$ cross-covariantiematrix, en $K_{**} = k(X_*, X_*)$ wees de $m \times m$ covariantiematrix testen. Het achterste is:

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

waar:

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

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

Het achterste gemiddelde $\bar{f}_*$ is een lineaire combinatie van de trainingsdoelen, gewogen op basis van de kernovereenkomst tussen test- en trainingspunten. De posterieure covariantie $\text{Cov}(f_*)$ begint bij de eerdere covariantie $K_{**}$ en trekt de verkregen informatie af van de trainingsgegevens. Waar de trainingsgegevens compact zijn, is de posterieure variantie klein. Waar trainingsgegevens schaars zijn, keert de posterieure variantie terug naar de eerdere.

De marginale waarschijnlijkheid en de claim die het testen waard is

De kernel-hyperparameters $\theta$ (lengteschalen, varianties, ruisniveau) worden geleerd door de log marginale waarschijnlijkheid te maximaliseren:

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

De eerste term is een data-fit-term (bestraft voorspellingen die verre van observaties zijn). De tweede term is een complexiteitsstraf (bestraft modellen die te flexibel zijn – d.w.z. waarbij de kernelmatrix een grote determinant heeft). De derde term is een normalisatieconstante.

Dit is het automatische scheermes van Occam, en het is het meest interessante wat een huisarts in een handelspijplijn brengt. De sterke versie van de bewering is dat er geen aparte validatieset nodig is voor regularisatie – de complexiteitsboete ligt binnen de doelstelling, dus het model kan niet gratis met flexibiliteit worden aangepast.

Die bewering verdient vooral op deze blog scepticisme. Plateau-analyse laat zien dat een validatiescore op één punt een slecht selectiecriterium is en dat robuustheid in de vorm van de buurt leeft; PBO kwantificeert hoe vaak de winnaar in de steekproef buiten de steekproef verliest; leeggelopen Sharpe prijzen in het aantal proeven. Een marginale waarschijnlijkheid is nog steeds een doelstelling binnen de steekproef waarover gemaximaliseerd wordt $\theta$ — de Occam-factor bestraft de capaciteit van het model, en niet de selectie over veel aangepaste kernels. Als je twaalf kandidaat-kernels past en degene kiest met de beste marginale waarschijnlijkheid, ben je terug in het gebied van meerdere tests en is de deflatielogica onveranderd van toepassing. De falsifieerbare versie: produceert selectie op marginale waarschijnlijkheid een kleinere kloof in de steekproef/buiten de steekproef dan selectie op basis van validatiesets op dezelfde gegevens en dezelfde kernelfamilie? Dat is meetbaar en wordt hieronder niet gemeten.

Het marginale waarschijnlijkheidsoppervlak heeft ook lokale optima. Meerdere willekeurige herstarts of zorgvuldige initialisatie zijn van belang: het initialiseren van de lengteschaal op de mediaan paarsgewijze afstand van de trainingsinputs en de ruisvariantie op de steekproefvariantie van de doelen is een redelijk uitgangspunt.

Computationele kosten en schaalbaarheid

De bottleneck is aan het omkeren $(K + \sigma_n^2 I)$, wat kost $O(n^3)$ op tijd en $O(n^2)$ in het geheugen. De kubusvormige muur wordt op deze blog al vanuit de andere richting beargumenteerd: zoekmethode vs. evaluatiekosten diskwalificeert GP-gebaseerde Bayesiaanse optimalisatie regelrecht wanneer het doel goedkoop is, omdat het surrogaat meer kost dan de evaluaties die het bespaart. De rekenkunde is hier hetzelfde; de toepassing is anders. Als een model dat is afgestemd op marktgegevens, is de kubieke term geen diskwalificatie maar een begroting: het stelt een hard plafond aan het trainingsvenster.

Schaalbaarheidsstrategieën voor handelsapplicaties:

1. Schaarse huisartsen (inducerende punten). Vervang de $n \times n$ matrix met een $m \times m$ matrix waar $m \ll n$. De inducerende punten $Z = \{z_1, \ldots, z_m\}$ zijn pseudo-invoer die de trainingsgegevens samenvat. De SVGP-formulering (Stochastic Variational GP) van Hensman et al. (2013) maakt mini-batchtraining tegen kosten mogelijk $O(nm^2)$ per iteratie. GPyTorch ondersteunt dit native.

2. Gestructureerde kernelinterpolatie (SKI/KISS-GP). Maakt gebruik van de Kronecker- en Toeplitz-structuur in de kernelmatrix wanneer invoer op een raster ligt. Verlaagt de kosten tot $O(n + g \log g)$ waar $g$ is de rastergrootte. Ideaal voor regelmatig gesamplede tijdreeksen (bijvoorbeeld maten van 1 minuut).

3. Lokale huisartsen op schuiframen. Afzonderlijke huisartsen trainen alleen op recente gegevens. Dit is waar de GP-specifieke beperking bijt: standaardkernels zijn stationair ($k(x,x')$ hangt er alleen van af $x - x'$) en markten niet, dus het gebruikelijke antwoord is een rollend venster – maar de vensterlengte wordt hierboven begrensd door $O(n^3)$, niet alleen door statistieken. Walk-forward-optimalisatie omvat verankerde versus rollende vensters, trein-/testlengtes en heroptimalisatiefrequentie in het algemeen; voor een huisarts is de venstergrootte zowel een rekenbeslissing als een statistische beslissing, en een verankerd (steeds groter wordend) venster is eenvoudigweg niet beschikbaar voorbij een paar duizend punten zonder benadering. Die herkadering is de praktische consequentie: bij huisartsen kun je niet kiezen voor 'gebruik de hele geschiedenis'.

GP voor rendementsvoorspelling: een praktisch raamwerk

Functieontwerp: waarom 5-20 functies, niet 500

De algemene kenmerktaxonomie voor ML van marktgegevens – onevenwicht in het orderboek, boekdruk, VPIN, Kyle’s lambda, kenmerken van gerealiseerde volatiliteit, cyclische tijdcodering, cross-asset- en financieringspercentagesignalen – is al uiteengezet in spread modeling with machine learning; gebruik die lijst.

Wat huisartsspecifiek is, is de grootte van de lijst. Kernelmethoden worden slechter in hoge dimensies: afstanden concentreren zich, en een stationaire kernel verliest discriminatie. Het praktische bereik voor een financiële huisarts is 5-20 inputs, niet de honderden die een boom met een gradiëntboost graag eet. Het mechanisme dat dit mogelijk maakt is ARD (Automatic Relevance Determination): geef elke invoerdimensie zijn eigen lengteschaal $\ell_d$en trainingen met marginale waarschijnlijkheid $\ell_d \to \infty$ voor dimensies die geen signaal dragen, aangezien een oneindige lengteschaal betekent dat de kernel die coördinaat negeert. Functieselectie wordt een bijproduct van het aanpassen en het geleerde $\ell_d$ zijn direct leesbaar als een relevantierangschikking - wat het ook falsifieerbaar maakt: pas de samengestelde kern op echte staven, druk de lengteschalen af ​​en kijk of de kenmerken waarin u gelooft, dezelfde zijn die het model behoudt.

Positiebepaling en onthouding

De GP posterior geeft $\sigma_*(x)$ direct, dus het valt rechtstreeks in het edge-ratio sizing en no-trade filter ontwikkeld in Conformal Prediction for Risk-Aware Position Sizing, met de intervalbreedte $w_t = 2\kappa\sigma_*(x)$. Het enige verschil is de herkomst: de huisarts haalt de breedte uit het model zelf in plaats van uit een kalibratieset, dus deze varieert afhankelijk van waar het testpunt zich bevindt ten opzichte van de trainingsgegevens in plaats van met een globaal kwantiel van residuen uit het verleden.

Implementatie met GPyTorch

GPyTorch is een op PyTorch gebaseerde bibliotheek voor schaalbare GP-inferentie. Het maakt gebruik van GPU-versnelling, automatische differentiatie en moderne lineaire algebratechnieken (conjugaatgradiënten, Lanczos-decompositie) om huisartsen verder te schalen dan het naïeve $O(n^3)$ beperken.

Basis Exact GP voor retourvoorspelling

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 wat de hierboven besproken lengteschalen per dimensie mogelijk maakt. Na de training, model.covar_matern.base_kernel.lengthscale is de vector die moet worden uitgelezen.

Trainingslus

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

Voorspelling met onzekerheid

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

De fast_pred_var() contextmanager gebruikt het LOVE-algoritme (Lanczos Variance Estimates) om voorspellende varianties te berekenen $O(n)$ tijd in plaats van $O(n^2)$.

End-to-end handelspijplijn

Opmerking over de functiebouwer hieronder: elke voortschrijdende statistiek moet strikt terugkijkend zijn, en de invoerstandaardisatie mag alleen op het trainingssegment worden aangebracht. Dat tweede punt is geen algemene hygiëne – het is precies het normalisatielekkage-kanaal dat is ontleed en gemeten in de look-ahead bias-taxonomie, dat de Sharpe-inflatie rapporteert die elk lektype produceert. Een huisarts standaardiseert zijn input per constructie, dus dit is het lek waaraan hij het meest wordt blootgesteld.

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

Schalen naar grotere datasets met schaarse huisartsen

Wanneer het trainingsvenster een paar duizend punten overschrijdt, wordt de exacte GP-afleiding traag. Schakel over naar een afwisselende, schaarse huisarts:

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

Met 128 inducerende punten zijn de kosten per batch gelijk $O(128^2 \times \text{batch\_size})$ — ongeveer 4 miljoen bewerkingen per batch. Dit verwerkt gemakkelijk datasets van meer dan 100.000 observaties op één GPU.

Diep kernelleren: huisarts ontmoet neurale netwerken

Wanneer de invoerruimte hoogdimensionaal is of de relatie tussen kenmerken en rendementen zeer niet-lineair is, kan een gewone kernel het moeilijk hebben. Deep kernel learning (DKL) geeft de invoer door via een neuraal netwerk voordat de GP-kernel wordt toegepast:

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

waar $g_\phi$ is een neuraal netwerk met parameters $\phi$ En $k_{\text{base}}$ is een standaardkernel (bijvoorbeeld Matern-3/2). Het netwerk leert een functierepresentatie kennen waarbij de GP-kernel het meest effectief is. Het hele model – netwerkparameters en kernel-hyperparameters – wordt end-to-end getraind door de marginale waarschijnlijkheid te maximaliseren.

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 combineert het representatieleren van neurale netwerken met de onzekerheidskwantificering van huisartsen. De GP-laag zorgt ervoor dat voorspellingen die ver verwijderd zijn van de trainingsgegevens een hoge onzekerheid hebben – iets waar standaard neurale netwerken om bekend niet in slagen te bieden. Merk ook op dat DKL precies de flexibiliteit herintroduceert die de marginale waarschijnlijkheid moest controleren: de Occam-factor bestraft de kernel, maar het netwerk ervoor heeft duizenden vrije parameters en een dergelijke straf bestaat niet.

Ludkovski en Risk (2025), Gaussian Process Models for Quantitative Finance, onderzoeken DKL in een bredere kwantitatieve financiële context, inclusief de prijsstelling van opties en portefeuilleoptimalisatie; die resultaten zijn van hen, op basis van hun problemen, en vormen geen bewijs voor de voorspelling van crypto-rendementen.

Diagnostiek en valkuilen

Kalibratie

Een goed gekalibreerde huisarts heeft voorspellende intervallen die overeenkomen met de empirische dekking. De controle is dezelfde als die wordt gebruikt in conformal forecast – die ook de theorie omvat waarom marginale dekking geen voorwaardelijke dekking is, een beperking die net zo goed van toepassing is op GP-intervallen:

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

De verwachte overschrijding bedraagt ​​0,317 / 0,046 / 0,003 bij 1/2/3 sigma. Het getal dat er toe doet, is wat dit afdrukt bij een echte walk-forward run over cryptobars, en die tabel staat nog niet in dit artikel. De eerdere verwachting is dat de huisarts overmoedig is – een Gaussiaanse waarschijnlijkheid op dikstaartige, niet-stationaire rendementen zou bij 3 sigma slecht undercover moeten zijn – maar ‘zou moeten’ is geen maatstaf.

Overige valkuilen

1. Invoerschaling. Lengteschalen zijn relatief ten opzichte van de invoerschaal, dus een kenmerk heeft een bereik $[0, 10000]$ en één varieerde $[0, 1]$ kan geen betekenisvolle delen $\ell$; ARD is wat heterogene bereiken redt, en standaardisatie is wat de initialisatie van ARD gezond maakt.

2. Overfitting van de marginale waarschijnlijkheid. Bij veel kernel-hyperparameters (vooral samengestelde of spectrale mengselkernels) kan de marginale waarschijnlijkheid nog steeds overfit zijn. Gebruik priors: een lognormale prior op lengteschalen gecentreerd op de mediaan paarsgewijze afstand, een halfnormale prior op varianties.

3. Covariantiematrixconditionering. $(K + \sigma_n^2 I)$ kan numeriek singulier worden als het geluid $\sigma_n^2$ te klein is of wanneer de trainingspunten bijna gedupliceerd zijn. GPyTorch voegt toe $10^{-6}$ jitter naar de diagonaal; slecht geconditioneerde financiële gegevens hebben vaak meer nodig.

4. Look-ahead bias. Behandeld hierboven en, in gemeten details, in de look-ahead bias-taxonomie.

Wanneer huisartsen versus andere modellen gebruiken

* Voor de XGBoost- en neurale-netwerkkolommen kun je het beste verwijzen naar de gepubliceerde vergelijkingen in plaats van naar een nieuwe bewering hier: spread modeling with machine learning heeft de gradiënt-boosting-vs-deep-learning tabel voor financiële tabelgegevens (drempels voor gegevensgrootte, interpreteerbaarheid, aanpassing van het regime, latentie) en temporele fusie transformers heeft TFT versus LSTM versus vanille Transformer.

Gebruik huisartsen wanneer:

• U beschikt over kleine tot middelgrote datasets (minder dan ~10.000 observaties per trainingsvenster)
• U wilt een expliciete, controleerbare structurele hypothese over het signaal (trend + seizoensinvloeden + ruis)
• De onzekerheid moet variëren met de afstand tot de trainingsgegevens, en niet alleen met een globaal restkwantiel

Gebruik geen huisarts wanneer:

• Je hebt inferentie van minder dan een milliseconde nodig over miljoenen waarnemingen
• De ingangsdimensionaliteit overschrijdt ~50
• Het signaal leeft in complexe functie-interacties van hoge orde die stationaire kernels niet kunnen vertegenwoordigen (DKL helpt, ten koste van de Occam-eigenschap)

Wat dit artikel nog niet meet

Alles hierboven zijn modelmachines. Niets hiervan is een bewijs dat een huisarts geld verdient met crypto, en de standaard van deze blog is dat het artikel zijn eigen cijfers bevat. De open items, in de volgorde waarin ze moeten worden uitgevoerd:

1. ARD-lengteschalen op echte BTCUSDT-staven van 1 m. Monteer de composietkernel, print $\ell_d$ per functie. Dit test direct de claim "ARD is ingebouwde functieselectie" en levert een falsifieerbare rangschikking op van de relevantie van functies.
2. De kalibratietabel van een walk-forward run. Verwachte vs. empirische overschrijding bij 1/2/3 sigma, duidelijk vermeld, inclusief het geval waarin de huisarts overmoedig blijkt te zijn.
3. De op vertrouwen gebaseerde strategie, walk-forward, op dezelfde vijf majors als het eerlijke negatieve, gedefleerd voor het aantal proeven. Als het mislukt, wordt het in die reeks opgenomen als een ander negatief resultaat.
4. De wandklok $O(n^3)$ curve voor $n = 250 / 500 / 1000 / 2000 / 5000$, exact vs. SVGP, dus de schaalbaarheidssectie rust op een diagram in plaats van op een bewering.

Conclusie

Het argument voor Gaussiaanse handelsprocessen is niet "ze geven je foutbalken" - conforme voorspelling geeft je foutbalken met minder verdelingsaannames en een dekkingsgarantie die de huisarts niet heeft. Het geval is dat een huisarts je je structurele hypothese laat opschrijven als een kernel, deze vergelijkt met een doelstelling die rekening houdt met zijn eigen complexiteit, en je vervolgens vertelt welke van je kenmerken hij daadwerkelijk heeft gebruikt via de ARD-lengteschalen. Dat is een ongebruikelijk leesbaar model.

De kosten zijn even concreet: kubieke schaalvergroting die je trainingsvenster afsluit, aannames over stationariteit die een rollend venster slechts gedeeltelijk repareert, een Gaussiaanse waarschijnlijkheid die met dikke staartrendementen zal worden geschonden, en – als je eenmaal reikt tot diepgaand kernelleren – het stille verlies van de Occam-eigenschap die in de eerste plaats de marginale waarschijnlijkheid motiveerde. Of wat overblijft verhandelbaar is, is een empirische vraag, en de vier hierboven genoemde metingen zouden deze vraag kunnen beantwoorden.

Referenties

• Rasmussen, C.E., & Williams, C.K.I. (2006). Gaussiaanse processen voor machinaal leren. MIT-pers.
• Wilson, A., en Adams, R. (2013). Gaussiaanse proceskernels voor patroondetectie en extrapolatie. ICML.
• Hensman, J., Fusi, N., & Lawrence, N.D. (2013). Gaussiaanse processen voor big data. UAI.
• Gatheral, J., Jaisson, T., en Rosenbaum, M. (2018). De volatiliteit is ruw. Kwantitatieve financiën, 18(6).
• Rizvi, S.A.A., Roberts, S.J., Osborne, M.A., & Nyikosa, F. (2017). Een nieuwe benadering voor het voorspellen van financiële volatiliteit met Gaussiaanse procesenveloppen. arXiv:1705.00891. – Ludkovski, M., & Risk, J. (2025). Gaussiaanse procesmodellen voor kwantitatieve financiën. Springer.
• Gardner, JR, Pleiss, G., Bindel, D., Weinberger, KQ, & Wilson, A.G. (2018). GPyTorch: Blackbox Matrix-Matrix Gaussiaanse procesinferentie met GPU-versnelling. NeurIPS.