Physics-Informed Neural Networks for Options Pricing

The interesting part of applying PINNs to derivatives is not Black-Scholes. Black-Scholes has a closed form; a network that reproduces it is a sanity check, not a result. The interesting part is everything that does not have a closed form and where finite differences run out of road: the Heston two-factor PDE with its mixed partial $\partial^2 V / \partial S \partial v$, the American free-boundary problem, and Merton's partial integro-differential equation whose integral term couples every point in the domain to every other point.

This post is about how each of those goes into a loss function. Concretely: the log-price residual that makes the Black-Scholes PINN trainable, the three-input Heston network whose cross-derivative falls out of autograd for free, the penalty-versus-two-network treatment of early exercise, and a Gauss-Hermite quadrature that puts the jump integral inside a PDE residual.

What it is not is a benchmark. Every number below that looks like a result is either a target or a citation, and it is marked as such. The measurement pass is scoped at the end.

PINNs were introduced by Raissi et al. (2019), and the mechanics are already covered on this blog in The Navier-Stokes Problem, where the same autograd-residual trick hunts for singularities in fluid equations: the PDE becomes a loss term, autodiff supplies exact derivatives of the network with respect to its inputs, and no grid is ever built. Here the only thing that changes is which PDE goes into the loss.

The Black-Scholes PDE as a Physics Constraint

The Black-Scholes PDE, its terminal payoff $\max(S-K,0)$, and the no-arbitrage argument behind it are covered in The Black-Scholes Formula — this section assumes them and goes straight to the form a network can actually be trained on.

Log-price transformation

Working directly with $S$ is problematic: the domain is $[0, \infty)$ and the PDE has variable coefficients. The standard transformation $x = \ln(S)$ converts it into a constant-coefficient form:

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

This is a convection-diffusion-reaction equation on $x \in (-\infty, \infty)$. The variable coefficients $S^2$ and $S$ are gone, which matters for a network: a residual whose magnitude scales with $S^2$ is dominated by the far end of the domain and the optimizer spends its budget there.

PINN loss for Black-Scholes

Let $u_\theta(x, t)$ be the network. Sample $N_r$ collocation points $(x_i, t_i)$ in the interior, $N_b$ on the boundaries, and $N_0$ at terminal time:

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

Every partial derivative is computed via torch.autograd.grad with create_graph=True, so gradients flow through the derivative computation during backpropagation. That single flag is the whole implementation trick.

PyTorch Implementation

Parameters below are crypto-flavored rather than textbook: $\sigma = 0.7$ and $T = 0.08$ (roughly a 30-day BTC option) instead of the equity-desk $\sigma = 0.2$, $T = 1.0$. Short maturity plus high vol is the regime where the payoff kink is sharpest and where the PINN is hardest to train — which is the point.

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

Two structural choices worth flagging:

• Weight on the terminal condition. The payoff gets weight 10x because it defines the problem. Without strong enforcement the network can satisfy the PDE trivially by outputting zero everywhere — the homogeneous PDE has infinitely many solutions and the terminal condition is what selects one.
• Resampling collocation points every epoch. Points are drawn fresh each step, which acts as stochastic regularization on the residual. The alternative — a fixed point set — lets the network overfit the residual at those exact coordinates.

Extending to the Heston Stochastic Volatility Model

Black-Scholes prices with a single constant $\sigma$, which is exactly the assumption GARCH(1,1) volatility forecasting is built to reject. Heston makes volatility a second stochastic factor, the instantaneous variance $v$:

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

where $\kappa$ is the mean-reversion speed, $\theta$ the long-run variance, $\xi$ the vol-of-vol, and $dW_1 \cdot dW_2 = \rho\,dt$.

The pricing PDE is two-dimensional in state:

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

This is where a meshfree method starts to look attractive. A finite difference scheme needs a grid in $(S, v, t)$ — a typical $200 \times 100 \times 500$ is $10^7$ nodes. Add a third stochastic factor, such as stochastic rates, and the grid becomes impractical. Monte Carlo scales better in dimension but converges slowly, especially for Greeks.

PINN architecture for Heston

The network takes three inputs $(x, v, t)$ with $x = \ln S$. The mixed partial $\partial^2 u / \partial x \partial v$ — the term that makes ADI schemes awkward — is one more autograd.grad call:

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

Note that u_xv is obtained by differentiating u_x with respect to v, reusing the graph built by the first call. Symmetry of the mixed partial means differentiating u_v with respect to x should give the same tensor; in float32 it will not, exactly, and the gap is a cheap diagnostic for whether the graph is behaving.

American Options and Free Boundary Problems

American options add an early exercise constraint: value must never fall below intrinsic. This turns the PDE into a free boundary problem, equivalently a linear complementarity problem (LCP):

$\frac{\partial V}{\partial t} + \mathcal{L}V \leq 0, \quad V \geq \Phi(S), \quad \left(\frac{\partial V}{\partial t} + \mathcal{L}V\right)(V - \Phi(S)) = 0$

where $\mathcal{L}$ is the Black-Scholes operator and $\Phi(S)$ the payoff. Two ways to put this in a loss.

Penalty method

Replace the complementarity constraint with a smooth penalty:

$\frac{\partial V}{\partial t} + \mathcal{L}V + \rho_p \cdot \max(\Phi(S) - V, 0) = 0$

with $\rho_p$ large (typically $10^4$ to $10^6$). When $V$ drops below intrinsic, the penalty forces it back up, and the PINN loss becomes:

$\mathcal{L}_{\text{PDE}}^{\text{American}} = \frac{1}{N_r}\sum_i \left[\frac{\partial u_\theta}{\partial t} + \mathcal{L}u_\theta + \rho_p \cdot \max(\Phi - u_\theta, 0)\right]^2$

No architecture change, just a modified residual. The cost is a new hyperparameter with a bad conditioning tradeoff: too small and the constraint is violated, too large and the loss landscape is dominated by the penalty term.

Direct free boundary approach

Train two networks jointly: one for the price $u_\theta(S, t)$, one for the optimal exercise boundary $S^*_\phi(t)$. The loss carries the PDE in the continuation region, the smooth-pasting condition at the boundary, and the payoff in the exercise region. The boundary comes out as a first-class output, which is what you actually want for hedging an American book.

Jump-Diffusion Models with PINNs

Merton's model adds Poisson jumps to geometric Brownian motion:

$dS = (r - \lambda \bar{k})S\,dt + \sigma S\,dW + S\,dJ$

where $J$ is a compound Poisson process with intensity $\lambda$ and log-normal jump sizes. Pricing becomes a partial integro-differential equation (PIDE):

$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + (r - \lambda\bar{k})S\frac{\partial V}{\partial S} - (r + \lambda)V + \lambda \int_0^\infty V(Sy, t) f(y)\,dy = 0$

with $f(y)$ the density of the jump multiplier.

The integral is what breaks finite differences — it couples every point in the domain to every other point, destroying the banded structure the solver relies on. A PINN does not have that structure to lose. The integral is just another term in the residual, evaluated by quadrature at each collocation point, and because the network is defined everywhere, $V(Sy, t)$ is a free forward pass rather than an interpolation:

def jump_integral(model, x, t, lam, mu_j, sigma_j, n_quad=32):
    """Approximate jump integral using Gauss-Hermite quadrature."""
    nodes, weights = np.polynomial.hermite.hermgauss(n_quad)
    nodes = torch.tensor(nodes, dtype=torch.float32, device=x.device)
    weights = torch.tensor(weights, dtype=torch.float32, device=x.device)

    y_nodes = mu_j + sigma_j * np.sqrt(2) * nodes
    integral = torch.zeros_like(x)

    for i in range(n_quad):
        x_shifted = x + y_nodes[i]
        v_shifted = model(x_shifted, t)
        integral += weights[i] * v_shifted

    integral *= 1.0 / np.sqrt(np.pi)
    return integral

The PIDE residual is then:

$\text{Residual} = \frac{\partial u}{\partial t} + \frac{1}{2}\sigma^2 u_{xx} + (r - \lambda\bar{k} - \frac{\sigma^2}{2})u_x - (r + \lambda)u + \lambda \cdot I[u]$

with $I[u]$ the quadrature approximation. The loop costs n_quad extra forward passes per training step, all of which land in the autograd graph — the real price of jumps is memory, not math.

Comparison with Traditional Methods

The structural differences are real and can be stated without a benchmark:

The performance rows that would normally live here — training time, inference latency, accuracy — are deliberately absent, because measuring them is the open work, not a table entry.

Where PINNs plausibly win

Real-time repricing. Once trained, evaluation is a forward pass over a batch, so an entire book reprices in one kernel launch rather than one solve per parameter set. Whether that beats a well-tuned Crank-Nicolson solver at realistic book sizes is an empirical question and unmeasured here.

High-dimensional models. Stochastic volatility plus stochastic rates plus jumps is 4+ dimensions, where finite differences are effectively dead. PINNs degrade gracefully in dimension, subject to training difficulty.

Continuous Greeks. The network is smooth and differentiable, so Delta, Gamma, Theta, and Vega come from the same autograd machinery as the PDE residual — no bump-and-revalue, no finite-difference noise. This is the strongest claim in the post and also the one most in need of a measured Gamma surface behind it.

Parametric solutions. Feeding model parameters ($\sigma$, $\kappa$, $\theta$) in as network inputs lets one PINN cover a family of models rather than a single calibration.

Where PINNs struggle

Training cost. For one option at one parameter set, a finite difference solve is over before a PINN finishes its first thousand epochs.

No error bounds. Finite differences with Richardson extrapolation reach machine precision with a known convergence order. A PINN reports a loss value, which is not an error bound. Uncertainty-aware variants (Bai et al., 2025) attach confidence intervals, but the field is young.

Optimization difficulty. The loss landscape is badly non-convex and the weighting between PDE, boundary, and terminal terms is a tuning problem. The characteristic failure is a network that drives the PDE residual to near zero while ignoring the boundary conditions entirely — a perfectly valid solution to the wrong problem.

Reproducibility. Different seeds, collocation distributions, and optimizer settings can land on meaningfully different solutions. Ensembling helps and multiplies cost.

Those last two are the parts of this post most worth turning into plots, because they are failure modes anyone reimplementing this will hit.

What Still Needs Measuring

The honest state of this article: the formulations are correct and the code runs, but nothing here has been benchmarked on this desk's hardware or on this desk's data. The pass that would make it a result rather than a derivation:

1. Sup-norm and RMSE against closed-form Black-Scholes across the $(S, t)$ grid, per seed, five seeds. The usual folklore target is agreement to 4+ decimal places; the point is to test that at $\sigma = 0.7$, $T = 0.08$ and publish the seeds where it fails.
2. Delta and Gamma error, not just price error. A network matching prices to four digits with a noisy $\Gamma$ is useless for hedging, so Gamma is the acceptance criterion.
3. Wall-clock training time and per-option inference latency for a 10,000-option book, against a Crank-Nicolson solve of the identical problem on the same GPU.
4. The failure modes, reproduced deliberately. Down-weight the boundary loss until the network satisfies the PDE and misses the boundary; run five seeds and plot the spread. Both are cheap to produce and more useful than another correct-looking price surface.
5. Fit to real Deribit BTC quotes, or at minimum keep the crypto-realistic parameters above rather than reverting to $\sigma = 0.2$, $T = 1.0$.
6. Architecture and optimizer choices as findings, not advice. Wider-versus-deeper, tanh versus ReLU (ReLU's discontinuous second derivative should visibly wreck the residual), $\lambda_{\text{IC}}$ sweeps, and Adam-then-L-BFGS versus Adam alone are all one-line experiments. Until they are run, they are folklore repeated from the literature, and this post declines to repeat them as recommendations.

The two structural claims that do survive without measurement, because they are properties of the formulation rather than of a run: use log-price coordinates (the PDE becomes constant-coefficient, which is a fact about the algebra), and weight the terminal condition above the PDE term (the homogeneous PDE admits the zero solution, which is a fact about the problem).

What Else Is Out There

Directions worth knowing about, none of them tested here:

• Residual-adaptive collocation sampling — place points where the residual is large, near the strike for Black-Scholes or near the Feller boundary for Heston, rather than uniformly. The implementation is short:

def adaptive_resample(model, x_pool, t_pool, n_select):
    """Select collocation points with highest PDE residual."""
    with torch.no_grad():
        residuals = compute_pde_residual(model, x_pool, t_pool).abs()
    probs = residuals.squeeze() / residuals.sum()
    indices = torch.multinomial(probs, n_select, replacement=False)
    return x_pool[indices], t_pool[indices]

• Transfer learning across strikes and maturities — fine-tune from a neighboring $(K, T)$ instead of training from scratch, to maintain a library covering the option chain.
• Physics-Informed Extreme Learning Machines — freeze hidden weights, train only the output layer, reducing training to one linear solve. Less expressive, but Black-Scholes is low-dimensional enough that it may not matter.
• Operator learning (DeepONet, FNO) — learn the solution operator rather than the solution, mapping $(\sigma(\cdot), r(\cdot), \text{payoff}(\cdot)) \to V(\cdot,\cdot)$ so new payoff structures need no retraining. See DeepSVM.

Conclusion

PINNs do not replace finite differences or Monte Carlo. For a single low-dimensional option at a single parameter set, finite differences win and it is not close. The case for PINNs is narrow and specific: high-dimensional models where grids die, parametric solutions reused across a book, and Greeks that come from the same autodiff pass as the price.

What this post delivers is the translation layer — how the Heston mixed partial, the American free boundary, and the Merton jump integral each become a term in a loss function. What it does not yet deliver is evidence that the resulting networks are accurate enough to hedge with. That is the next pass, and until it runs, treat everything above as a derivation rather than a recommendation.

References and Further Reading

• Raissi, M., Perdikaris, P., Karniadakis, G.E. (2019). Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378, 686-707.
• Dhiman, A., Kaur, D. (2023). Physics Informed Neural Network for Option Pricing. arXiv preprint.
• Salvador, B. et al. (2023). Meshless methods for American option pricing through Physics-Informed Neural Networks. Engineering Analysis with Boundary Elements.
• Bai, G. et al. (2025). An uncertainty-aware physics-informed neural network solution for the Black-Scholes equation. arXiv preprint.
• Chen, Y. et al. (2025). Stochastic jump diffusion process informed neural networks for American option pricing under data scarcity. Applied Soft Computing.
• Benth, F.E. et al. (2025). DeepSVM: Learning Stochastic Volatility Models with Physics-Informed Deep Operator Networks. arXiv preprint.
• He, Y. et al. (2024). Finance-Informed Neural Network: Learning the Geometry of Option Pricing. arXiv preprint.
• MathWorks: Physics-Informed Neural Networks (PINNs) for Option Pricing. MATLAB Finance Blog, 2025.
• GitHub: BlackScholesPINN (PyTorch implementation).
• Raissi PINNs GitHub repository.