5 Python libraries that do the math quant finance actually runs on

Search "top quant Python libraries" and you'll get the same list every time. Pandas. NumPy. yfinance. A few mostly-abandoned wrappers around someone's PhD thesis. None of them actually do the math.

Strip away the jargon and quant finance is a small number of hard mathematical problems. You're solving differential equations. You're solving optimization problems. You're doing Bayesian inference over latent state. You're pricing derivatives. That's most of the work. Pick a real research paper from a derivatives desk, an HFT shop, or a systematic fund and you'll see one of these four problems sitting underneath it.

This article is about the five Python libraries that handle those problems seriously. Not toys. Not 50-line wrappers. The libraries that production quants actually use when they need to ship.

• Diffrax for stochastic differential equations, with reverse-mode autodiff through the solver.
• py-pde for finite-difference PDEs with symbolic equation specification and numba JIT.
• CVXPY for convex optimization, all the way up to mixed-integer quadratic programs.
• NumPyro for Bayesian inference with NUTS, JAX-compiled, fast enough for real time-series posteriors.
• FinancePy for derivatives pricing across equity, rates, FX, and credit, with the models desks actually use.

For each one, I'll show what it is, the problem it solves, and the core code. The visuals in this article are all real output from the libraries themselves. The carousel on Instagram walked through them as a teaser; this is the deep version with code you can actually run.

Setup

Everything below runs in a single Python 3.13 environment. The JAX-based libraries (Diffrax, NumPyro) install CPU JAX by default. CVXPY needs the SCIP backend for the mixed-integer step. py-pde brings numba in as a dep.

1. Diffrax: differentiable SDEs in JAX

Diffrax is a JAX-native library for numerical differential equations, written by Patrick Kidger (the same person who introduced Neural SDEs in his thesis). It covers ODEs, SDEs, and controlled differential equations, and the defining feature is that every solver is reverse-mode autodifferentiable. You can compute exact gradients of the solution with respect to your model parameters by backpropagating through diffeqsolve itself. It composes with jit and vmap and runs identically on CPU, GPU, and TPU.

The problem

The Heston model couples a risk-neutral price $S_t$ to a mean-reverting instantaneous variance $v_t$:

with correlated drivers $\langle dW^{S},\, dW^{v}\rangle = \rho\, dt$. There are two things you want to do with this model. The first is the forward problem: simulate paths, take expectations, price an option as

That part every quant library can do. Here it is in Diffrax:

The backward problem, and why this library exists

The second thing you want to do is calibrate: given a quoted strike chain $\{C_k^{\text{mkt}}\}$, find the parameters $\Theta^\star$ that minimize the model-market squared error,

For gradient-based calibration you need $\nabla_\Theta C$. Finite differences re-price the chain twice per parameter and the gradient gets dominated by Monte-Carlo variance, so calibration becomes flaky and slow. Characteristic-function pricing works for Heston specifically but doesn't generalize to rough volatility, jumps, or any model where you can't write the MGF in closed form.

Diffrax computes the gradient exactly by backpropagating through the SDE solver itself. The same price function below is differentiated with jax.grad, then optimized with optax.adam.

The second image is what Diffrax was built to enable. Gradient-based calibration of an SDE model with no finite differences and no model-specific Fourier shortcut. Swap the Heston dynamics for rough Bergomi or a jump-diffusion and the calibration code is the same.

2. py-pde: finite-difference PDEs in pure Python

py-pde is a finite-difference PDE solver written in Python with a numba-compiled hot loop. The standout feature is the PDE class. You write the evolution equation as a string in mathematical notation, the library parses it, and numba compiles a fast solver. It supports coupled multi-field systems, custom and time-dependent boundary conditions, stochastic PDEs (SPDEs via SDEBase), and 2D/3D grids in Cartesian, polar, cylindrical, and spherical geometries.

The problem

Two species $u$ and $v$ diffuse and react on a periodic 2D plane:

$F$ feeds species $u$ from a reservoir; $k$ removes species $v$. For certain regions of the $(F, k)$ plane, a uniform mixture is Turing-unstable: tiny perturbations grow into self-replicating spots, stripes, or labyrinths. The parameters used here, $(F, k) = (0.054,\, 0.063)$, sit in the labyrinthine regime, and the filaments above are the steady state.

This is not finance, but the machinery is. The same finite-difference solver that runs Gray-Scott runs the Black-Scholes PDE backward in time, the Hamilton-Jacobi-Bellman equation for optimal execution, and the Dupire local-volatility equation. The PDE here just happens to be visually extraordinary.

Five lines specify the PDE, one line solves it, the MemoryStorage tracker stores intermediate frames so you can render the evolution as an animation. A from-scratch numpy port is fast to write but roughly 10x slower at runtime, and you'd reinvent the storage and animation pipeline yourself.

3. CVXPY: convex optimization as a DSL

CVXPY is a domain-specific language for convex optimization, developed by Stephen Boyd's group at Stanford. You write the problem in math notation, the library verifies convexity through the Disciplined Convex Programming (DCP) ruleset, transforms it into standard cone form, and dispatches it to a solver. It covers LP, QP, SOCP, SDP, MILP, MIQP, robust optimization, and chance constraints in the same syntax. The solver list is plug-and-play: ECOS, SCS, OSQP, CLARABEL, SCIP, CBC, plus the commercial ones (GUROBI, MOSEK).

The problem

Classical mean-variance portfolio selection picks weights $w \in \mathbb{R}^{n}$ to trade off expected return against variance:

That's a quadratic program. Three lines in CVXPY, solved in milliseconds.

Real portfolios face an additional operational, tax, or regulatory constraint: hold at most $k$ non-zero positions. Adding $\|w\|_{0} \le k$ turns the QP into a mixed-integer quadratic program (MIQP), enforced through binary indicators $z \in \{0,1\}^{n}$ with $w_i \le z_i$ and $\sum z_i \le k$.

What CVXPY lets you do is add the integer constraint by changing one keyword. Same problem object, same solve call, the library re-dispatches to SCIP under the hood.

The story the visual tells: at $k = 2$ the frontier sits well inside the unconstrained one, at $k = 3$ slightly less so, and by $k = 5$ the curve has essentially caught up. The cardinality constraint stops binding around five names for this universe. Translation: you can hold a 30-stock universe with five concentrated positions and lose almost no Sharpe. Going to three costs real return at every risk level. Going to two is a serious haircut.

The whole optimization is about thirteen lines. Re-implementing this against a raw MIQP solver (Gurobi's Python API, say) is five times the code and pinned to one backend.

4. NumPyro: gradient-based MCMC for real time-series posteriors

NumPyro is a JAX-backed probabilistic programming language. You write a generative model with numpyro.sample(...) statements and numpyro.contrib.control_flow.scan for time-series; NumPyro gives you the No-U-Turn variant of HMC (NUTS, the adaptive-trajectory sampler that powers most modern Bayesian inference) plus variational inference (SVI) and a large distributions library. NUTS through JAX autodiff is one of the fastest mainstream gradient-based MCMC implementations available.

The problem

A latent log-variance $h_t$ follows an AR(1) process, and conditional on $h_t$ the returns $y_t$ are Gaussian with time-varying scale:

The inference target is the joint posterior $p(\mu,\, \phi,\, \sigma_{\eta},\, \mathbf{h} \mid \mathbf{y})$ over both the structural parameters and the full latent volatility path $\mathbf{h} = (h_1, \ldots, h_T)$. For $T = 500$ that's a 503-dimensional posterior with very strong correlations between consecutive $h_t$. Random-walk samplers stall on this geometry. NUTS uses the gradient of the log-posterior with respect to all 503 latent variables to pick its step direction, and JAX provides that gradient through autodiff.

The same autodiff machinery that backpropagates through diffeqsolve in section 1 lets NUTS compute the gradient of the log-posterior with respect to all 503 latent variables in a single backward pass. That gradient is what tells the sampler which way to move. Stan does the same thing but in its own language. PyMC does it about 5x slower for high-dim time-series posteriors. For Bayesian SV models, NumPyro is the live option.

5. FinancePy: derivatives pricing and the models desks actually trade

FinancePy is a pure-Python derivatives pricing and risk library by Dominic O'Kane (ex-Lehman, now NYU). The coverage is broad. Equity options, fixed income, FX, credit, inflation, structured products. And it ships with the models that desks actually use: Black-Scholes, Heston, SABR, Bachelier, Hull-White, LMM. Pricing methods include closed-form (Heston via characteristic function), Monte Carlo, PDE, and tree-based. The library bootstraps yield curves and fits implied vol surfaces out of the box.

The problem

SABR (Hagan et al., 2002) parameterizes the implied-volatility surface through four parameters $(\alpha, \beta, \rho, \nu)$. Under the SABR dynamics

Hagan's asymptotic expansion gives a closed-form Black implied volatility $\sigma_{\text{BS}}(K, T;\, \alpha, \beta, \rho, \nu, F)$. Fast enough to calibrate a smile in milliseconds. With $\beta$ fixed at $1$ (lognormal SABR), each maturity slice has three free parameters: $\alpha$ for the ATM vol level, $\rho$ for the skew, $\nu$ for the vol of vol.

Implementing the Hagan formula from the paper is about 50 lines including edge cases (the ATM limit, $\beta = 0$ and $\beta = 1$ branches, near-money handling). FinancePy gives you the same thing, debugged and tested, in one line.

Two states of the same model, rendered on the same axes. ATM vol decreases with maturity in both (the surface flattens going back). The smile is steepest at short maturity (steep ramp on the left at $T \to 0$). Skew is negative throughout, sharply so in the stressed regime. These are the features every options desk monitors. SABR is also the model that actually trades on every interest-rate, FX, and equity-vol desk worldwide. FinancePy's SABR(alpha, beta, rho, nu).black_vol(F, K, T) is one line.

Note the connection back to section 1. The Diffrax visual fit Heston by backpropagating through a Monte Carlo SDE solver. This visual fits SABR using a closed-form expansion plus scipy.optimize.minimize. Same fundamental problem (calibrate a vol model to market quotes), two completely different solver paths, two different libraries. Both are how it really gets done.

Putting the toolkit together

If you took every quant paper published in the last five years and looked at the libraries actually used, this is most of what you'd see, with PyTorch occupying the rest of the JAX share for ML-heavy work. The five solve four different problems:

• Differential equations, including stochastic ones, where the model parameters need calibration. Diffrax.
• Coupled finite-difference PDEs for option pricing, optimal control, local-vol fitting. py-pde.
• Convex optimization with portfolio constraints, hedging, execution, risk budgeting. CVXPY.
• Bayesian inference over latent state in a non-trivial geometry. NumPyro.
• Derivatives pricing with the standard models, calibration utilities, and risk metrics. FinancePy.

Real workflows use them in combination. A research desk might calibrate Heston with Diffrax (because they want to extend to a rough-volatility model where Fourier methods fail), price the calibrated model with FinancePy's analytic Heston implementation for sanity checks, then run portfolio-level risk through CVXPY. None of these libraries are toys, and none of them are doing what NumPy alone can do.

That's the deep version of the carousel. If you want more of this kind of breakdown, with code and visuals from real libraries, follow along on Instagram. The next one is already in progress.