Physikinformierte neuronale Netze für die Preisgestaltung von Optionen

Der interessante Teil der Anwendung von PINNs auf Derivate ist nicht Black-Scholes. Black-Scholes hat eine geschlossene Form; Ein Netzwerk, das es reproduziert, ist eine Plausibilitätsprüfung, kein Ergebnis. Der interessante Teil ist alles, was keine geschlossene Form hat und bei dem endliche Differenzen außer Kontrolle geraten: die Heston-Zweifaktor-PDE mit ihrem gemischten Teil $\partial^2 V / \partial S \partial v$, das amerikanische Problem der freien Grenzen und Mertons partielle Integro-Differentialgleichung, deren Integralterm jeden Punkt in der Domäne mit jedem anderen Punkt koppelt.

In diesem Beitrag geht es darum, wie jede davon in eine Verlustfunktion übergeht. Konkret: das Logarithmus-Preis-Residuum, das die Black-Scholes-PINN trainierbar macht, das Heston-Netzwerk mit drei Eingängen, dessen Kreuzableitung kostenlos aus Autograd herausfällt, die Strafe-gegen-zwei-Netzwerk-Behandlung der frühen Übung und eine Gauß-Hermite-Quadratur, die das Sprungintegral in ein PDE-Residuum einfügt.

Was es nicht ist, ist ein Maßstab. Jede Zahl darunter, die wie ein Ergebnis aussieht, ist entweder ein Ziel oder eine Zitierung und wird als solche gekennzeichnet. Am Ende wird der Messdurchgang skaliert.

PINNs wurden von Raissi et al. (2019), und die Mechanik wird bereits in diesem Blog in Das Navier-Stokes-Problem behandelt, wo der gleiche Autograd-Residuen-Trick nach Singularitäten in Fluidgleichungen sucht: Die PDE wird zu einem Verlustterm, Autodiff liefert exakte Ableitungen des Netzwerks in Bezug auf seine Eingaben und es wird nie ein Gitter erstellt. Hier ändert sich lediglich, welche PDE in den Verlust geht.

Die Black-Scholes-PDE als physikalische Einschränkung

Die Black-Scholes-PDE, ihre endgültige Auszahlung $\max(S-K,0)$, und das Argument der No-Arbitrage dahinter werden in Die Black-Scholes-Formel behandelt – dieser Abschnitt geht davon aus und geht direkt auf die Form ein, auf der ein Netzwerk tatsächlich trainiert werden kann.

Logarithmische Preistransformation

Direkte Zusammenarbeit mit $S$ ist problematisch: Die Domain ist $[0, \infty)$ und die PDE hat variable Koeffizienten. Die Standardtransformation $x = \ln(S)$ wandelt es in eine Form mit konstantem Koeffizienten um:

$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 \frac{\partial^2 V}{\partial x^2} + \left(r - \frac{1}{2}\sigma^2\right)\frac{\partial V}{\partial x} - rV = 0$

Dabei handelt es sich um eine Konvektions-Diffusions-Reaktionsgleichung $x \in (-\infty, \infty)$. Die variablen Koeffizienten $S^2$ Und $S$ sind weg, was für ein Netzwerk wichtig ist: ein Residuum, dessen Größe mit skaliert $S^2$ wird vom anderen Ende der Domäne dominiert und der Optimierer gibt sein Budget dort aus.

PINN-Verlust für Black-Scholes

Lass $u_\theta(x, t)$ sei das Netzwerk. Probe $N_r$ Kollokationspunkte $(x_i, t_i)$ im Innenraum, $N_b$ an den Grenzen, und $N_0$ zur Endzeit:

$\mathcal{L}_{\text{PDE}} = \frac{1}{N_r}\sum_{i=1}^{N_r}\left[\frac{\partial u_\theta}{\partial t} + \frac{1}{2}\sigma^2 \frac{\partial^2 u_\theta}{\partial x^2} + \left(r - \frac{\sigma^2}{2}\right)\frac{\partial u_\theta}{\partial x} - r u_\theta\right]^2$

$\mathcal{L}_{\text{IC}} = \frac{1}{N_0}\sum_{j=1}^{N_0}\left[u_\theta(x_j, T) - \max(e^{x_j} - K, 0)\right]^2$

$\mathcal{L}_{\text{BC}} = \frac{1}{N_b}\sum_{k=1}^{N_b}\left[u_\theta(x_{\min}, t_k)\right]^2 + \frac{1}{N_b}\sum_{k=1}^{N_b}\left[u_\theta(x_{\max}, t_k) - (e^{x_{\max}} - Ke^{-r(T-t_k)})\right]^2$

Jede partielle Ableitung wird über berechnet torch.autograd.grad mit create_graph=True, sodass während der Rückausbreitung Gradienten durch die Ableitungsberechnung fließen. Dieses einzelne Flag ist der ganze Implementierungstrick.

PyTorch-Implementierung

Die folgenden Parameter sind eher kryptotypisch als lehrbuchmäßig: $\sigma = 0.7$ Und $T = 0.08$ (ungefähr eine 30-Tage-BTC-Option) anstelle des Equity-Desk $\sigma = 0.2$, $T = 1.0$. Kurze Laufzeit plus hohe Volatilität ist das Regime, in dem der Auszahlungsknick am stärksten ist und die PINN am schwierigsten zu trainieren ist – und genau darum geht es.

import torch
import torch.nn as nn
import numpy as np

r = 0.05        # risk-free rate
sigma = 0.7     # volatility
K = 1.0         # strike (moneyness-normalized)
T = 0.08        # ~30 days
x_min, x_max = -1.0, 1.0  # log-price domain (S/K in [0.37, 2.72])

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")


class BSPINN(nn.Module):
    """Physics-Informed Neural Network for Black-Scholes PDE."""

    def __init__(self, hidden_dim=128, num_layers=4):
        super().__init__()
        layers = [nn.Linear(2, hidden_dim), nn.Tanh()]
        for _ in range(num_layers - 1):
            layers += [nn.Linear(hidden_dim, hidden_dim), nn.Tanh()]
        layers.append(nn.Linear(hidden_dim, 1))
        self.net = nn.Sequential(*layers)

    def forward(self, x, t):
        inputs = torch.cat([x, t], dim=1)
        return self.net(inputs)


def compute_pde_residual(model, x, t):
    """Compute Black-Scholes PDE residual using autodiff."""
    x.requires_grad_(True)
    t.requires_grad_(True)
    u = model(x, t)

    grads = torch.autograd.grad(u, [x, t], grad_outputs=torch.ones_like(u),
                                create_graph=True)
    u_x, u_t = grads[0], grads[1]

    u_xx = torch.autograd.grad(u_x, x, grad_outputs=torch.ones_like(u_x),
                               create_graph=True)[0]

    residual = u_t + 0.5 * sigma**2 * u_xx + (r - 0.5 * sigma**2) * u_x - r * u
    return residual


def terminal_condition(x):
    """European call payoff: max(S - K, 0) = max(exp(x) - K, 0)."""
    return torch.relu(torch.exp(x) - K)


def train_pinn(epochs=10000, lr=1e-3, n_interior=5000, n_boundary=500, n_terminal=1000):
    model = BSPINN().to(device)
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=epochs)

    for epoch in range(epochs):
        optimizer.zero_grad()

        x_int = (torch.rand(n_interior, 1, device=device)
                 * (x_max - x_min) + x_min)
        t_int = torch.rand(n_interior, 1, device=device) * T

        residual = compute_pde_residual(model, x_int, t_int)
        loss_pde = (residual ** 2).mean()

        x_tc = (torch.rand(n_terminal, 1, device=device)
                * (x_max - x_min) + x_min)
        t_tc = torch.ones(n_terminal, 1, device=device) * T
        u_tc = model(x_tc, t_tc)
        loss_ic = ((u_tc - terminal_condition(x_tc)) ** 2).mean()

        t_bc = torch.rand(n_boundary, 1, device=device) * T

        x_lo = torch.full((n_boundary, 1), x_min, device=device)
        loss_bc_lo = (model(x_lo, t_bc) ** 2).mean()

        x_hi = torch.full((n_boundary, 1), x_max, device=device)
        target_hi = torch.exp(x_hi) - K * torch.exp(-r * (T - t_bc))
        loss_bc_hi = ((model(x_hi, t_bc) - target_hi) ** 2).mean()

        loss = loss_pde + 10.0 * loss_ic + loss_bc_lo + loss_bc_hi

        loss.backward()
        optimizer.step()
        scheduler.step()

        if epoch % 1000 == 0:
            print(f"Epoch {epoch:5d} | PDE: {loss_pde:.2e} | "
                  f"IC: {loss_ic:.2e} | BC: {loss_bc_lo + loss_bc_hi:.2e}")

    return model

Zwei strukturelle Optionen, die es wert sind, hervorgehoben zu werden:

• Gewichtung der Endbedingung. Die Auszahlung erhält das 10-fache Gewicht, da sie das Problem definiert. Ohne starke Durchsetzung kann das Netzwerk die PDE trivial erfüllen, indem es überall Nullen ausgibt – die homogene PDE hat unendlich viele Lösungen und die Endbedingung ist es, die eine auswählt.
• Resampling von Kollokationspunkten in jeder Epoche. Punkte werden bei jedem Schritt neu gezeichnet, was als stochastische Regularisierung für das Residuum fungiert. Die Alternative – ein Festpunktsatz – ermöglicht dem Netzwerk eine Überanpassung des Residuums an genau diesen Koordinaten.

Erweiterung auf das stochastische Volatilitätsmodell von Heston

Black-Scholes-Preise mit einer einzigen Konstante $\sigma$, was genau die Annahme ist, die GARCH(1,1) Volatility Forecasting ablehnen soll. Heston macht die Volatilität zu einem zweiten stochastischen Faktor, der momentanen Varianz $v$:

$dS = rS\,dt + \sqrt{v}S\,dW_1$ $dv = \kappa(\theta - v)\,dt + \xi\sqrt{v}\,dW_2$

Wo $\kappa$ ist die mittlere Umkehrgeschwindigkeit, $\theta$ die langfristige Varianz, $\xi$ das Vol-of-Vol, und $dW_1 \cdot dW_2 = \rho\,dt$.

Die Preis-PDE hat einen zweidimensionalen Zustand:

$\frac{\partial V}{\partial t} + \frac{1}{2}vS^2\frac{\partial^2 V}{\partial S^2} + \rho\xi v S\frac{\partial^2 V}{\partial S \partial v} + \frac{1}{2}\xi^2 v\frac{\partial^2 V}{\partial v^2} + rS\frac{\partial V}{\partial S} + \kappa(\theta - v)\frac{\partial V}{\partial v} - rV = 0$

Hier beginnt eine netzfreie Methode attraktiv zu erscheinen. Ein Finite-Differenzen-Schema benötigt ein Gitter $(S, v, t)$ – ein typisches $200 \times 100 \times 500$ Ist $10^7$ Knoten. Fügen Sie einen dritten stochastischen Faktor hinzu, beispielsweise stochastische Raten, und das Raster wird unpraktisch. Monte Carlo lässt sich in der Dimension besser skalieren, konvergiert jedoch langsam, insbesondere für Griechen.

PINN-Architektur für Heston

Das Netzwerk benötigt drei Eingaben $(x, v, t)$ mit $x = \ln S$. Der gemischte Teil $\partial^2 u / \partial x \partial v$ – der Begriff, der ADI-Systeme umständlich macht – ist ein weiterer autograd.grad Anruf:

class HestonPINN(nn.Module):
    def __init__(self, hidden_dim=256, num_layers=5):
        super().__init__()
        layers = [nn.Linear(3, hidden_dim), nn.Tanh()]
        for _ in range(num_layers - 1):
            layers += [nn.Linear(hidden_dim, hidden_dim), nn.Tanh()]
        layers.append(nn.Linear(hidden_dim, 1))
        self.net = nn.Sequential(*layers)

    def forward(self, x, v, t):
        return self.net(torch.cat([x, v, t], dim=1))


def heston_pde_residual(model, x, v, t, kappa, theta, xi, rho, r):
    x.requires_grad_(True)
    v.requires_grad_(True)
    t.requires_grad_(True)
    u = model(x, v, t)

    u_x, u_v, u_t = torch.autograd.grad(
        u, [x, v, t], torch.ones_like(u), create_graph=True
    )
    u_xx = torch.autograd.grad(u_x, x, torch.ones_like(u_x), create_graph=True)[0]
    u_vv = torch.autograd.grad(u_v, v, torch.ones_like(u_v), create_graph=True)[0]
    u_xv = torch.autograd.grad(u_x, v, torch.ones_like(u_x), create_graph=True)[0]

    residual = (
        u_t
        + 0.5 * v * u_xx
        + rho * xi * v * u_xv
        + 0.5 * xi**2 * v * u_vv
        + (r - 0.5 * v) * u_x
        + kappa * (theta - v) * u_v
        - r * u
    )
    return residual

Beachten Sie, dass u_xv erhält man durch Differenzieren u_x in Bezug auf v, wobei der beim ersten Aufruf erstellte Graph wiederverwendet wird. Symmetrie des gemischten Teils bedeutet Differenzierung u_v in Bezug auf x sollte den gleichen Tensor ergeben; In float32 ist dies nicht genau der Fall, und die Lücke ist eine kostengünstige Diagnose dafür, ob sich das Diagramm verhält.