Deriving the Black-Scholes PDE from First Principles

Introduction

The Black-Scholes model is arguably the most important result in quantitative finance. Published by Fischer Black and Myron Scholes in 1973 (with Robert Merton's parallel independent derivation), it provided the first rigorous, no-arbitrage framework for pricing European options.

In this article we derive the Black-Scholes PDE from scratch — using only geometric Brownian motion, Itô's lemma, and a delta-hedging argument. No shortcuts. We then solve the PDE analytically and implement the full pricing model in Python, including the Greeks.

The Model: Geometric Brownian Motion

We model the stock price $S_t$ as a geometric Brownian motion (GBM):

dS_t = \mu S_t \, dt + \sigma S_t \, dW_t

where:

$\mu \in \mathbb{R}$ is the drift (expected return per unit time)
$\sigma > 0$ is the volatility (standard deviation of returns)
$W_t$ is a standard Brownian motion on a filtered probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}, \mathbb{P})$

The solution to this SDE is:

S_t = S_0 \exp\!\left[\left(\mu - \frac{\sigma^2}{2}\right)t + \sigma W_t\right]

The term $-\frac{\sigma^2}{2}$ is the Itô correction — it arises because log-returns are normally distributed but the exponential introduces a Jensen's inequality gap.

Itô's Lemma

Let $V(t, S_t)$ be a twice continuously differentiable function of time and the stock price (i.e., an option price). Itô's lemma gives us the stochastic differential of $V$ :

dV = \frac{\partial V}{\partial t} dt + \frac{\partial V}{\partial S} dS_t + \frac{1}{2} \frac{\partial^2 V}{\partial S^2} (dS_t)^2

Using the multiplication rules $dt \cdot dt = 0$ , $dt \cdot dW_t = 0$ , and $(dW_t)^2 = dt$ :

(dS_t)^2 = \sigma^2 S_t^2 \, dt

Substituting the GBM dynamics:

dV = \left(\frac{\partial V}{\partial t} + \mu S \frac{\partial V}{\partial S} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\right) dt + \sigma S \frac{\partial V}{\partial S} \, dW_t

The Delta-Hedging Argument

Construct a self-financing portfolio $\Pi$ consisting of one option and $-\Delta$ shares of stock:

\Pi = V - \Delta \cdot S

The portfolio change over an infinitesimal interval $dt$ is:

d\Pi = dV - \Delta \, dS_t

Substituting the expressions derived above:

d\Pi = \left(\frac{\partial V}{\partial t} + \mu S \frac{\partial V}{\partial S} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\right) dt + \sigma S \frac{\partial V}{\partial S} \, dW_t - \Delta\left(\mu S \, dt + \sigma S \, dW_t\right)

Key insight: choose $\Delta = \frac{\partial V}{\partial S}$ (the option's delta). This cancels all stochastic terms:

d\Pi = \left(\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\right) dt

The portfolio is now instantaneously riskless — it has no exposure to $dW_t$ .

No-Arbitrage and the PDE

A riskless portfolio must earn the risk-free rate $r$ to prevent arbitrage:

d\Pi = r \Pi \, dt

Substituting $\Pi = V - \frac{\partial V}{\partial S} \cdot S$ :

\left(\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\right) dt = r\left(V - S \frac{\partial V}{\partial S}\right) dt

Dividing through by $dt$ and rearranging gives the Black-Scholes PDE:

\boxed{\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0}

This PDE holds for any European-style derivative $V(t, S)$ — regardless of $\mu$ ! The drift has disappeared entirely. This is the risk-neutral pricing miracle.

Boundary Conditions

For a European call option with strike $K$ and expiry $T$ :

V(T, S) = \max(S - K, 0)

For a European put:

V(T, S) = \max(K - S, 0)

Additional boundary conditions (for the call):

V(t, 0) = 0, \qquad V(t, S) \approx S - Ke^{-r(T-t)} \text{ as } S \to \infty

The Closed-Form Solution

Via a change of variables that reduces the BS PDE to the heat equation, the closed-form price for a European call is:

C(S, t) = S \, \mathcal{N}(d_1) - K e^{-r(T-t)} \mathcal{N}(d_2)

where $\mathcal{N}(\cdot)$ is the standard normal CDF and:

d_1 = \frac{\ln(S/K) + \left(r + \frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}, \qquad d_2 = d_1 - \sigma\sqrt{T-t}

By put-call parity $C - P = Se^{-q(T-t)} - Ke^{-r(T-t)}$ (where $q$ is the continuous dividend yield), the put price is:

P(S, t) = Ke^{-r(T-t)}\mathcal{N}(-d_2) - S \, \mathcal{N}(-d_1)

Python Implementation

import numpy as np
from scipy.stats import norm
from dataclasses import dataclass
from typing import Literal
 
@dataclass
class BSMOption:
    S: float          # Spot price
    K: float          # Strike price
    T: float          # Time to expiry (years)
    r: float          # Risk-free rate (annualised)
    sigma: float      # Volatility (annualised)
    q: float = 0.0    # Continuous dividend yield
    kind: Literal["call", "put"] = "call"
 
    def _d1_d2(self) -> tuple[float, float]:
        d1 = (np.log(self.S / self.K) + (self.r - self.q + 0.5 * self.sigma**2) * self.T) / (
            self.sigma * np.sqrt(self.T)
        )
        d2 = d1 - self.sigma * np.sqrt(self.T)
        return d1, d2
 
    def price(self) -> float:
        d1, d2 = self._d1_d2()
        df = np.exp(-self.r * self.T)
        fwd = self.S * np.exp(-self.q * self.T)
        if self.kind == "call":
            return fwd * norm.cdf(d1) - self.K * df * norm.cdf(d2)
        return self.K * df * norm.cdf(-d2) - fwd * norm.cdf(-d1)
 
    # ── Greeks ──────────────────────────────────────────────
    def delta(self) -> float:
        d1, _ = self._d1_d2()
        disc = np.exp(-self.q * self.T)
        if self.kind == "call":
            return disc * norm.cdf(d1)
        return -disc * norm.cdf(-d1)
 
    def gamma(self) -> float:
        d1, _ = self._d1_d2()
        return (
            np.exp(-self.q * self.T)
            * norm.pdf(d1)
            / (self.S * self.sigma * np.sqrt(self.T))
        )
 
    def vega(self) -> float:
        d1, _ = self._d1_d2()
        return (
            self.S
            * np.exp(-self.q * self.T)
            * norm.pdf(d1)
            * np.sqrt(self.T)
            / 100  # per 1% move in vol
        )
 
    def theta(self) -> float:
        d1, d2 = self._d1_d2()
        fwd = self.S * np.exp(-self.q * self.T)
        term1 = -fwd * norm.pdf(d1) * self.sigma / (2 * np.sqrt(self.T))
        if self.kind == "call":
            return (term1 - self.r * self.K * np.exp(-self.r * self.T) * norm.cdf(d2)
                    + self.q * fwd * norm.cdf(d1)) / 365
        return (term1 + self.r * self.K * np.exp(-self.r * self.T) * norm.cdf(-d2)
                - self.q * fwd * norm.cdf(-d1)) / 365
 
    def rho(self) -> float:
        _, d2 = self._d1_d2()
        if self.kind == "call":
            return self.K * self.T * np.exp(-self.r * self.T) * norm.cdf(d2) / 100
        return -self.K * self.T * np.exp(-self.r * self.T) * norm.cdf(-d2) / 100

Let's price an at-the-money call and examine the Greeks:

opt = BSMOption(S=100, K=100, T=0.5, r=0.05, sigma=0.20, kind="call")
 
print(f"Price  : {opt.price():.4f}")
print(f"Delta  : {opt.delta():.4f}")
print(f"Gamma  : {opt.gamma():.6f}")
print(f"Vega   : {opt.vega():.4f}  (per 1% σ)")
print(f"Theta  : {opt.theta():.4f}  (per calendar day)")
print(f"Rho    : {opt.rho():.4f}   (per 1% r)")

Output:

Price  : 6.8887
Delta  : 0.5377
Gamma  : 0.019919
Vega   : 0.2801  (per 1% σ)
Theta  : -0.0188  (per calendar day)
Rho    : 0.2352   (per 1% r)

Key Intuitions

Why does $\mu$ vanish?

The drift $\mu$ cancels because we're pricing under the risk-neutral measure $\mathbb{Q}$ (also called the equivalent martingale measure). Under $\mathbb{Q}$ , all assets grow at the risk-free rate $r$ . The derivation above is the discrete-time hedging intuition for this profound result.

The $\frac{1}{2}\sigma^2$ term

This is the Itô correction — the quadratic variation of GBM contributes directly to the drift of any smooth function of $S$ . A financial consequence: the expected return on an option exceeds what a naive $\Delta$ -based argument would suggest.

PDE vs. Risk-Neutral Pricing

The PDE approach and the risk-neutral expectations approach are equivalent:

V(t, S) = e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}\!\left[\,V(T, S_T) \;\middle|\; \mathcal{F}_t\right]

The Feynman-Kac theorem formalises this connection — it states that the solution to the BS PDE is this conditional expectation under $\mathbb{Q}$ .

Limitations and Extensions

The Black-Scholes model makes assumptions that fail in real markets:

Assumption	Reality	Extension
Constant $\sigma$	Volatility smiles/skews	Local Vol, Heston, SABR
Continuous trading	Jump risk exists	Jump-diffusion (Merton, Kou)
No transaction costs	Friction in hedging	Leland's model
European exercise	American options exist	Free-boundary PDE, LSM
Log-normal returns	Fat tails observed	Variance-Gamma, NIG

Despite its limitations, Black-Scholes remains the market's lingua franca. Implied volatility — the $\sigma$ that makes the BS formula match the market price — is how options are quoted and compared globally.

In a future article we'll derive the Heston stochastic volatility model and its semi-analytic characteristic function pricing formula.