Black–Scholes from First Principles: Deriving the Formula by Replication

We price a European call on a non-dividend stock. The model rests on a short list of assumptions, and every one of them is a place the real world will later disagree: the stock follows a geometric Brownian motion with constant volatility σ, the risk-free rate r is constant, trading is continuous and frictionless, and — the load-bearing one — no arbitrage is allowed.

Under those rules the stock’s dynamics are dS = μS·dt + σS·dW. Notice what is not there: the drift μ. The replication argument is about to delete it, and that deletion is the whole magic trick.

Before believing “constant σ”, ask the data. QQQ’s trailing one-year realised volatility over 2017-01-03 → 2024-12-31 ranged from 10.4% to 36.3% — the “constant” more than tripled inside one sample (the trailing one-month vol peaked far higher still around 2020-04-06). Keep that crack in mind; the smile in section 5 grows out of it.

Hold one option worth V(S, t) and short Δ shares. Apply Itô’s lemma to V, then choose Δ = ∂V/∂S so that the random dW term cancels exactly. The portfolio is now instantaneously riskless — and a riskless portfolio, by no-arbitrage, must earn precisely the risk-free rate. Equate the two and the drift falls out, leaving a deterministic equation in V.

The economic content is the cancellation: because we can hedge the risk away, the option’s fair value cannot depend on how bullish or bearish anyone feels about the stock. Only volatility and the rate survive.

Solving the PDE with the call payoff max(S − K, 0) gives the closed form. Everything funnels through two standardised distances, d₁ andd₂:

Read it as a probability-weighted payoff: N(d₂) is (roughly) the chance the option finishes in the money, and S·N(d₁) is the expected stock you receive, both under the risk-neutral measure. Plotted across spot for our real contract (K = 505, T = 1y, σ = 18.0%, r = 4.21%), the formula smooths the hockey-stick payoff — the vertical gap between the curves is time value, and the slope of the smooth one is delta:

The implementation is a direct transcription, and this time every input is observed: spot is QQQ’s last close of $507.45, the strike is the nearest listed 505, σ is the trailing one-year realised vol of 18.04%, and r is the 13-week T-bill yield of 4.21% (2024-12-31, ^IRX). The Greeks — the sensitivities that tell a desk how its book moves — are just the analytic derivatives of C, so they come almost for free.

import numpy as np
from scipy.stats import norm

def bs_price(S, K, T, r, sigma, kind="call"):
    """European Black-Scholes price. S spot, K strike, T years, r rate, sigma vol."""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    if kind == "call":
        return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)

def bs_greeks(S, K, T, r, sigma):
    """Delta, Gamma, Vega (per 1% vol), Theta (per day), Rho (per 1% rate)."""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    pdf = norm.pdf(d1)
    return {
        "delta": norm.cdf(d1),
        "gamma": pdf / (S * sigma * np.sqrt(T)),
        "vega":  S * pdf * np.sqrt(T) / 100,
        "theta": (-(S * pdf * sigma) / (2 * np.sqrt(T))
                  - r * K * np.exp(-r * T) * norm.cdf(d2)) / 365,
        "rho":   K * T * np.exp(-r * T) * norm.cdf(d2) / 100,
    }

S0, K, T = 507.45, 505, 1.0     # QQQ close + nearest strike
r, sigma = 0.0421, 0.1804       # ^IRX and 1y realised vol
print(bs_price(S0, K, T, r, sigma))   # 48.45

For this contract the desk sees (put–call parity holds to < 1e-8, the free correctness check):

One assumption is a known fiction: constant volatility. If it were true, every strike on a name would imply the same σ. To see why it fails, look at the returns themselves: QQQ’s daily log returns carry an excess kurtosis of 6.78 (a normal distribution has 0), and the worst day in the sample — -12.8% on 2020-03-16 — was a 8.8σ event under the model’s own calibration. A two-state Gaussian mixture fits those returns far better than one lognormal: a calm regime (77% of days, σ ≈ 14.4%) and a stress regime (23% of days, σ ≈ 39.0%). Price options under that mixture and invert each price back through Black–Scholes, and the implied σ is no longer flat:

The numbers off that curve: 20.1% at the money, 22% at 85% moneyness, 20.5% at 115% — a genuine smile generated by nothing more exotic than two volatility regimes. Real equity smiles are steeper still and asymmetric — OTM puts trade extra rich because crashes are fatter and faster than any symmetric mixture allows — but the mechanism is exactly this one: fat tails force the wings up.

So is the model useless? No — it is the language. Traders quote in implied vol rather than price precisely because Black–Scholes gives an invertible, one-number translation. You stop believing the assumptions and start using the formula as a coordinate system. That shift — from belief to convention — is the real lesson, and the gateway to local-vol, stochastic-vol, and everything that repairs the smile.