The Heston model

When volatility moves too

Black-Scholes assumes constant volatility, but real volatility rises in crashes and fades in calm, and the options market shows it: implied volatility differs by strike (the smile or skew). The Heston model keeps the geometric price process but makes the *variance* v random too: it is pulled back toward a long-run level theta at speed kappa, it is shaken by noise of size xi that scales with sqrt(v), and the shocks to price and variance have correlation rho. A negative rho (price falls as volatility rises, the leverage effect) produces the downward skew seen in equity options.

# The Heston equations, in words:  dS = mu S dt + sqrt(v) S dW      dv = kappa (theta - v) dt + xi sqrt(v) dZ      corr(dW, dZ) = rho
kappa, theta, xi = 2.0, 0.04, 0.3             # reversion speed, long-run variance (a 20% volatility), volatility of variance
print(round(theta ** 0.5, 2))                 # 0.2    the long-run volatility
print(2 * kappa * theta, xi ** 2)             # 0.16 0.09   the Feller condition 2 kappa theta >= xi^2 holds: the variance stays positive

Simulating it: truncation and correlation

To simulate, step through small intervals dt. The variance update is v + kappa * (theta - v) * dt + xi * sqrt(v * dt) * dz, but a plain Euler step can push v below zero, where the square root fails. The full-truncation scheme uses max(v, 0) wherever v appears on the right-hand side (so a negative value acts as zero) and keeps the raw result for the next step. The Feller condition 2 * kappa * theta >= xi^2 is the parameter regime where the true process never hits zero; if it fails, truncation matters even more. The two shocks must be correlated: draw independent normals z1, z2, set dz = z1 and dw = rho * z1 + sqrt(1 - rho^2) * z2.

# One variance step with full truncation, and the correlated shocks for the price
import math

def variance_step(v, kappa, theta, xi, dz, dt):
    vp = max(v, 0.0)
    return v + kappa * (theta - vp) * dt + xi * math.sqrt(vp * dt) * dz

def correlated_normals(z1, z2, rho):
    return z1, rho * z1 + math.sqrt(1 - rho ** 2) * z2

print(round(variance_step(0.04, 2.0, 0.04, 0.3, 0.5, 1 / 252), 6))      # 0.04189   at the long-run level only the noise moves it
dz, dw = correlated_normals(1.0, 1.0, -0.7)
print(round(dz, 4), round(dw, 4))             # 1.0 0.0141   a negative rho: the price shock pulls against the variance shock
print(variance_step(-0.01, 2.0, 0.04, 0.3, 1.0, 1 / 252))               # -0.00968...   still negative, but it moved up by kappa * theta * dt: the negative value acted as 0

A full path and a sanity check

One step updates both state variables together: the price uses the variance at the *start* of the step, S * exp((mu - v/2) * dt + sqrt(v * dt) * dw), then the variance moves. Repeat for the number of steps in the horizon and you have a path; many paths give a distribution of terminal prices, and averaging a discounted payoff over them prices an option. A good sanity check: with the drift mu, the average terminal price should be close to S0 * exp(mu * T), whatever the volatility parameters. If it is not, there is a bug in the price step. Calibrating kappa, theta, xi, rho to option prices comes in the topics ahead.

# A short Monte Carlo: simulate many paths, and check the average terminal price against S0 * exp(mu * T)
import math
import random

def simulate(s0, v0, mu, kappa, theta, xi, rho, steps, years, rng):
    dt = years / steps
    s, v = s0, v0
    for _ in range(steps):
        z1, z2 = rng.gauss(0, 1), rng.gauss(0, 1)
        dz, dw = z1, rho * z1 + math.sqrt(1 - rho ** 2) * z2
        vp = max(v, 0.0)
        s *= math.exp((mu - vp / 2) * dt + math.sqrt(vp * dt) * dw)
        v = v + kappa * (theta - vp) * dt + xi * math.sqrt(vp * dt) * dz
    return s

rng = random.Random(11)
finals = [simulate(100.0, 0.04, 0.05, 2.0, 0.04, 0.3, -0.7, 50, 1.0, rng) for _ in range(4000)]
mean = sum(finals) / len(finals)
print(abs(mean / (100.0 * math.exp(0.05)) - 1) < 0.02)       # True   the average terminal price follows S0 * exp(mu * T)
print(min(finals) > 0)                                       # True   prices stay positive

# A check that costs nothing: with xi = 0 and v0 = theta the variance never moves (Black-Scholes with volatility sqrt(theta))
import math

kappa, theta, xi = 2.0, 0.04, 0.0
v = theta
for z in (1.5, -2.0, 0.3, 0.9):
    vp = max(v, 0.0)
    v = v + kappa * (theta - vp) * (1 / 252) + xi * math.sqrt(vp / 252) * z
print(v == theta)                                    # True   no noise and already at the long-run level
print(round(math.sqrt(v), 2))                        # 0.2    a constant 20% volatility

• Heston makes the variance random: it mean-reverts to theta at speed kappa, has noise xi * sqrt(v), and is correlated with the price by rho (negative for equities).
• Simulate with full truncation (max(v, 0) on the right-hand side), correlate shocks with dw = rho z1 + sqrt(1 - rho^2) z2, and check the Feller condition 2 kappa theta >= xi^2.
• The price step uses the variance at the start of the step; check that the average terminal price follows S0 * exp(mu * T) and that xi = 0 reproduces Black-Scholes.

Next up: local volatility, a volatility surface that depends on price and time.

Lab · what's graded

• Keep the variance positive with feller_condition25 pts
• Step the variance with variance_step25 pts
• Correlate the shocks with correlated_normals25 pts
• Move price and variance with heston_step25 pts