Matthew JacksonMJ

Project 01 · Derivatives pricing

Options pricing & volatility

Black-Scholes-Merton and Cox-Ross-Rubinstein pricers built from first principles, analytical Greeks, and a Newton-Raphson implied volatility solver that turns a live options chain into a volatility surface.

Code on GitHub
Fig. 1. European call under BSM, \(K = 100\), decaying toward its payoff as expiry nears. Move across the figure to change spot; up and down changes volatility.

(1)Black-Scholes-Merton

1.1

Paths to a distribution

The model assumes the underlying follows geometric Brownian motion, the same process drawn on the home page. Pricing happens under the risk-neutral measure, where the drift is the risk-free rate:

\[ dS_t = r S_t\, dt + \sigma S_t\, dW_t \quad\Longrightarrow\quad S_T = S_0\, e^{(r - \frac{1}{2}\sigma^2)T + \sigma\sqrt{T}\,Z} \]

Every path ends somewhere, and the endpoints pile up into a lognormal distribution.

1.2

Payoff, average, discount

A call pays \((S_T - K)^+\) at expiry: nothing for the paths that finish below the strike, the excess for the rest. Its price today is the risk-neutral expectation of that payoff, discounted at the risk-free rate:

\[ C = e^{-rT}\, \mathbb{E}^{\mathbb{Q}}\big[(S_T - K)^+\big] \]

Averaging the simulated payoffs and discounting gives a Monte Carlo price. The closed form does the same integral exactly.

1.3

Hedge away the randomness

Why is the price an expectation under the risk-free drift? Itô's lemma gives the change in an option's value \(V(S, t)\) when \(S\) follows GBM with any drift \(\mu\):

\[ dV = \Big(\frac{\partial V}{\partial t} + \mu S \frac{\partial V}{\partial S} + \tfrac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\Big)dt + \sigma S \frac{\partial V}{\partial S}\, dW \]

Hold the option and sell \(\Delta = \partial V/\partial S\) shares. The \(dW\) terms cancel, and so does \(\mu\):

\[ \Pi = V - \Delta S, \qquad d\Pi = \Big(\frac{\partial V}{\partial t} + \tfrac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\Big)dt \]

The figure runs it on one path that drifts at 8%, not 4%: the option alone swings by dollars while the hedged book barely moves. What's left is the cost of rebalancing only once a day, and it grows when the hedge is adjusted weekly.

1.4

The Black-Scholes equation

A riskless portfolio has to earn the risk-free rate, \(d\Pi = r\Pi\,dt\). Setting the two expressions equal gives the Black-Scholes partial differential equation:

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

Each term has a meaning: time decay, the value of convexity in a moving market, growth at the risk-free rate, and the cost of financing the position. They cancel exactly at every point, which the waterfall checks as \(S\) sweeps across. The real-world drift never appears, which is why pricing can use \(r\) in its place.

1.5

Solving it backward

The payoff fixes the value at expiry, \(V(S, T) = (S - K)^+\), and the equation carries it back to today. In time to expiry \(\tau = T - t\) it's a diffusion, close cousin of the heat equation, so the payoff's sharp kink smooths out as \(\tau\) grows.

\[ \frac{\partial V}{\partial \tau} = \tfrac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV, \qquad V(S, 0) = (S - K)^+ \]

The figure solves it with Crank-Nicolson: on a grid in \(S\) and \(\tau\), each new value comes from three neighbours one step back, averaging an explicit and an implicit step, which leaves a tridiagonal system to solve at every step:

\[ \frac{V^{n+1}_i - V^{n}_i}{\Delta\tau} = \tfrac{1}{2}\big(\mathcal{L}V^{n+1}\big)_i + \tfrac{1}{2}\big(\mathcal{L}V^{n}\big)_i \]
1.6

The closed form

The equation can also be solved exactly. Its solution is the discounted risk-neutral expectation, and with \(S_T = S e^{(r - \sigma^2/2)T + \sigma\sqrt{T}Z}\) the payoff is positive exactly when \(Z > -d_2\). The integral splits in two:

\[ C = e^{-rT}\,\mathbb{E}\big[S_T\,\mathbf{1}_{Z > -d_2}\big] - K e^{-rT}\,\mathbb{P}(Z > -d_2) \]

The second piece is the strike times \(N(d_2)\), the risk-neutral probability of exercise. In the first, \(e^{-rT}S_T\,\varphi(z) = S\,\varphi(z - \sigma\sqrt{T})\): weighting by the stock price slides the bell curve right by \(\sigma\sqrt{T}\), and the same region's area becomes \(N(d_1)\).

\[ \begin{gathered} C = S e^{-qT} N(d_1) - K e^{-rT} N(d_2) \\[4pt] d_1 = \frac{\ln(S/K) + \big(r - q + \tfrac{1}{2}\sigma^2\big)T}{\sigma\sqrt{T}}, \qquad d_2 = d_1 - \sigma\sqrt{T} \end{gathered} \]

Puts follow from put-call parity, which the implementation satisfies to better than \(10^{-10}\).

(2)Binomial tree

2.1

Build the lattice

The Cox-Ross-Rubinstein tree splits the time to expiry into \(N\) steps of \(\Delta t = T/N\). At each step the price moves up by \(u\) or down by \(d\), with risk-neutral probability \(p\) of an up move:

\[ u = e^{\sigma\sqrt{\Delta t}}, \qquad d = \frac{1}{u}, \qquad p = \frac{e^{(r-q)\Delta t} - d}{u - d} \]

Because \(ud = 1\) the tree recombines, leaving \(N + 1\) terminal prices instead of \(2^N\).

2.2

Payoffs at expiry

At the last step the option's value is known exactly. For a put, it's the intrinsic value at each terminal price:

\[ V_{N,j} = \max\big(K - S_0\, u^{j} d^{N-j},\; 0\big) \]
2.3

Backward induction

Each earlier node is worth the discounted expectation of its two children, so values flow back through the tree to today:

\[ V_{n,j} = e^{-r\Delta t}\big[\,p\, V_{n+1,\,j+1} + (1-p)\, V_{n+1,\,j}\,\big] \]
2.4

Early exercise

An American option can be exercised at any node, so each node takes whichever is worth more, holding on or exercising now:

\[ V_{n,j} = \max\Big(e^{-r\Delta t}\big[p\, V_{n+1,\,j+1} + (1-p)\, V_{n+1,\,j}\big],\; K - S_{n,j}\Big) \]

The highlighted nodes are where exercising wins: deep in the money, where the interest on the strike is worth more than the option's remaining time value. On the engine's at-the-money benchmark the early-exercise premium is about $0.52.

2.5

Convergence

As \(N\) grows the tree converges to the Black-Scholes price, zigzagging as the strike falls between terminal nodes. The error shrinks like \(1/N\); on the engine's own benchmark it's under \(10^{-2}\) at \(N = 200\).

\[ \big|\,V^{\text{tree}}_N - V^{\text{BSM}}\big| = O(1/N) \]

The page navigator is drawn as one of these trees. Each section is a node, and the path through it is one possible price path.

(3)Greeks

3.1

Delta is a slope

Each Greek is a derivative of the price. Delta is the slope of the price against the spot, the number of shares that hedges one option:

\[ \Delta = \frac{\partial C}{\partial S} = e^{-qT} N(d_1) \]
3.2

Gamma is curvature

Gamma is how fast delta changes. It peaks at the money and sharpens as expiry approaches, when the option is about to become either worthless or a share:

\[ \Gamma = \frac{\partial^2 C}{\partial S^2} = \frac{e^{-qT} N'(d_1)}{S\sigma\sqrt{T}} \]
GreekDefinitionWhat it measures
Delta \(\Delta\)\(\partial V/\partial S\)The hedge ratio
Gamma \(\Gamma\)\(\partial^2 V/\partial S^2\)How fast delta changes
Vega \(\nu\)\(\partial V/\partial \sigma\)Sensitivity to volatility
Theta \(\Theta\)\(\partial V/\partial t\)Time decay
Rho \(\rho\)\(\partial V/\partial r\)Interest rate sensitivity
3.3

Checking with finite differences

All five Greeks are computed analytically and checked against central finite differences. As the step \(h\) shrinks the secant closes on the tangent and the error falls like \(h^2\), until floating-point rounding takes over and it climbs again:

\[ \Delta \approx \frac{C(S+h) - C(S-h)}{2h}, \qquad \text{error} \approx \underbrace{\tfrac{h^2}{6}\,C'''}_{\text{truncation}} + \underbrace{\tfrac{\varepsilon\, C}{h}}_{\text{rounding}} \]

The engine's analytical and numerical Greeks agree to better than \(10^{-7}\).

(4)Implied volatility

4.1

Inverting the price

The BSM price rises monotonically with \(\sigma\), so for any market price there is exactly one volatility that reproduces it. That volatility is the implied volatility:

\[ \operatorname{BSM}(\sigma^{*}) = C_{\text{mkt}} \]
4.2

Newton-Raphson

Newton's method follows the tangent of the price curve down to the market price. Vega is that tangent's slope, so each step is one division:

\[ \sigma_{n+1} = \sigma_n - \frac{\operatorname{BSM}(\sigma_n) - C_{\text{mkt}}}{\nu(\sigma_n)} \]

Convergence is quadratic: the number of correct digits roughly doubles each step, typically reaching machine precision in 5 to 10 iterations. A bisection fallback handles contracts with near-zero vega, and the solver is vectorized with NumPy broadcasting across whole chains.

4.3

The smile and the surface

Live chains come from yfinance, filtered for liquidity and priced at the bid-ask midpoint, and every contract is inverted. Under BSM the result would be flat across strikes. It isn't: equity index options such as SPY show a pronounced skew, with lower strikes carrying higher implied volatility, which reflects demand for downside protection.

Across expiries the smiles form a surface, which the dashboard draws in 3D. The figure's market is synthetic, but each point on it is a real Newton inversion.

(5)Results

Live SPY chain
< 4 s3,800+ contracts across 32 expiries, fetched and IV-solved
Put-call parity
\(< 10^{-10}\)Holds to machine precision
Greeks vs finite differences
\(< 10^{-7}\)All five, analytical against numerical
IV round trip
\(< 10^{-8}\)Recovery error across all strikes and volatility levels

(6)Code

Python, NumPy, SciPy, pandas, yfinance, Streamlit, and Plotly. The dashboard has two tabs: a pricer with all five Greeks plotted against spot and time to expiry, and a market data tab that fetches a ticker's chain and draws its implied volatility surface.

Next: Heston stochastic volatility to fit the surface, Dupire local volatility, delta-hedging P&L attribution, and a strategy builder.

Code on GitHub