Alphanume
Guest
DatasetsToolsDocsPricingContact
Sign inSign up

Insights

Black-Scholes in Python (From Scratch)

Alphanume Team · June 8, 2026

A clean, vectorized implementation you can trust.

The Black-Scholes model is the bedrock of modern options pricing, yet most implementations floating around the internet cut corners in ways that bite you the moment you feed in edge cases or try to price an entire strike chain at once. This post builds black scholes python code the right way: a scalar version you can read and verify, followed by a vectorized version that handles arrays without a single Python loop. We will also confirm correctness with put-call parity, so you know the math is right before you do anything else with it. If you want to jump straight to pricing without the code, our options pricing calculator does it interactively.

The Black-Scholes formula

Before writing any code, it helps to pin down exactly what we are computing. The price of a European call option under Black-Scholes is:

C = S · N(d1) − K · e−rT · N(d2)

and the corresponding put is:

P = K · e−rT · N(−d2) − S · N(−d1)

where the two auxiliary quantities are:

d1 = [ln(S/K) + (r + σ2/2) · T] / (σ · √T)

d2 = d1 − σ · √T

Here S is the spot price, K is the strike, r is the continuously compounded risk-free rate, σ is implied volatility, and T is time to expiry in years. N(·) is the standard normal CDF. A fuller derivation of each term lives in the companion post on the Black-Scholes formula; here we focus on turning that math into reliable code.

Scalar implementation

The cleanest starting point is a pure-scalar function that takes five numbers and returns two — one call price, one put price. scipy.stats.norm.cdf gives us N(·) without any hand-rolled approximation, and math.log and math.sqrt keep the intent obvious.

import math
from scipy.stats import norm


def bs_price(S, K, T, r, sigma):
    """Black-Scholes call and put prices for a single option.

    Parameters
    ----------
    S     : float  spot price
    K     : float  strike price
    T     : float  time to expiry in years
    r     : float  continuously compounded risk-free rate
    sigma : float  implied volatility (annualised)

    Returns
    -------
    call : float
    put  : float
    """
    if T <= 0:
        call = max(S - K, 0.0)
        put = max(K - S, 0.0)
        return call, put

    if sigma <= 0:
        pv_k = K * math.exp(-r * T)
        call = max(S - pv_k, 0.0)
        put = max(pv_k - S, 0.0)
        return call, put

    sqrt_T = math.sqrt(T)
    d1 = (math.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sqrt_T)
    d2 = d1 - sigma * sqrt_T

    call = S * norm.cdf(d1) - K * math.exp(-r * T) * norm.cdf(d2)
    put = K * math.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    return call, put

Two edge cases are handled explicitly at the top. When T <= 0 the option has expired and its value is just intrinsic — the payoff you would collect by exercising right now. When sigma <= 0 there is no uncertainty, so the option price collapses to the present-value comparison of spot and discounted strike. Both checks prevent a division-by-zero inside the d1 formula and make the function safe to call in any loop over expiry dates.

Vectorized implementation with NumPy

Pricing one option at a time is fine for exploration, but pricing a whole chain — say, every strike at a given expiry — calls for vectorization. NumPy's np.where handles the edge-case branching element-wise, and scipy.stats.norm.cdf already accepts arrays, so the structure mirrors the scalar version exactly.

import numpy as np
from scipy.stats import norm as sp_norm


def bs_price_vec(S, K, T, r, sigma):
    """Vectorized Black-Scholes prices.

    All arguments may be scalars or broadcastable NumPy arrays.

    Returns
    -------
    call : ndarray
    put  : ndarray
    """
    S = np.asarray(S, dtype=float)
    K = np.asarray(K, dtype=float)
    T = np.asarray(T, dtype=float)
    r = np.asarray(r, dtype=float)
    sigma = np.asarray(sigma, dtype=float)

    sqrt_T = np.sqrt(np.maximum(T, 0.0))

    # Avoid log(0) and division by zero by working only where valid
    safe_T = np.where(T > 0, T, 1.0)
    safe_sig = np.where(sigma > 0, sigma, 1.0)
    safe_sqrt_T = np.sqrt(safe_T)

    d1 = (
        np.log(S / K) + (r + 0.5 * safe_sig ** 2) * safe_T
    ) / (safe_sig * safe_sqrt_T)
    d2 = d1 - safe_sig * safe_sqrt_T

    pv_k = K * np.exp(-r * T)

    call_bs = S * sp_norm.cdf(d1) - pv_k * sp_norm.cdf(d2)
    put_bs = pv_k * sp_norm.cdf(-d2) - S * sp_norm.cdf(-d1)

    # Intrinsic values for edge cases
    call_intrinsic = np.maximum(S - K, 0.0)
    put_intrinsic = np.maximum(K - S, 0.0)
    call_zero_sig = np.maximum(S - pv_k, 0.0)
    put_zero_sig = np.maximum(pv_k - S, 0.0)

    call = np.where(T <= 0, call_intrinsic,
           np.where(sigma <= 0, call_zero_sig, call_bs))
    put = np.where(T <= 0, put_intrinsic,
          np.where(sigma <= 0, put_zero_sig, put_bs))

    return call, put

The safe_T and safe_sig substitutions let NumPy compute d1 and d2 everywhere without branching — the results are only used where T > 0 and sigma > 0, guarded by the final np.where calls. Broadcasting rules mean you can pass a scalar spot price and an array of strikes and get back a price array for the full chain in one call.

Worked example and put-call parity check

Concretely: AAPL is at $180, we want the 185-strike option expiring in 45 days, with a risk-free rate of 5.25% and implied volatility of 22%.

S = 180.0
K = 185.0
T = 45 / 365          # ~0.1233 years
r = 0.0525
sigma = 0.22

call, put = bs_price(S, K, T, r, sigma)
print(f"Call: ${call:.4f}")
print(f"Put:  ${put:.4f}")

# Put-call parity: C - P == S - K * exp(-r * T)
parity_lhs = call - put
parity_rhs = S - K * math.exp(-r * T)
print(f"C - P:           {parity_lhs:.6f}")
print(f"S - K*exp(-rT):  {parity_rhs:.6f}")
print(f"Parity holds: {abs(parity_lhs - parity_rhs) < 1e-10}")

Running that block should print something close to Call: $3.0268, Put: $7.8476, and Parity holds: True. Put-call parity — C − P = S − K · e−rT — must hold algebraically for any Black-Scholes implementation because it follows directly from the formula structure, not from any numerical approximation. If it fails, something is wrong in your d1/d2 signs or your norm.cdf arguments. Treat it as the first unit test for any option pricing code you write.

Assumptions and extensions

The formulas above assume a European-exercise option on a non-dividend-paying stock. Both assumptions are worth keeping in mind before applying this to real data.

Dividends. For a stock paying a continuous dividend yield q, replace every bare S in the formula with S * exp(-q * T) — the spot discounted by the foregone dividends over the option's life. Concretely, d1 becomes:

d1 = [ln(S/K) + (r − q + σ2/2) · T] / (σ · √T)

and the call price becomes S * exp(-q*T) * N(d1) - K * exp(-r*T) * N(d2). Adding a single q=0.0 default parameter to either function above is all the code change required.

American options. Black-Scholes has no closed form for American-exercise options. For those you need a binomial tree, finite-difference grid, or the Barone-Adesi–Whaley approximation — each a separate post.

Greeks. Once you have the d1 and d2 values in hand, delta, gamma, vega, theta, and rho all follow analytically from the same inputs. That is covered in detail in the post on computing option Greeks in Python.

The functions here are intentionally small — one computation each, no global state, no hidden dependencies. That makes them easy to unit-test, easy to vectorize further, and easy to slot into a larger pricing engine without surprises.