Chapter 13 — Volatility
"Volatility is the only thing in finance that mean-reverts reliably — everything else is wishful thinking."
After this chapter you will be able to:
- Estimate realised volatility using close-to-close, Parkinson, and Garman-Klass estimators and explain when each is appropriate
- Fit and simulate GARCH(1,1) models for volatility dynamics via maximum likelihood estimation in OCaml
- Calibrate the Heston stochastic volatility model using characteristic functions and the COS pricing method
- Parametrise an implied volatility surface using SVI (Stochastic Volatility Inspired) and Dupire's local volatility construction
- Detect and enforce calendar spread and butterfly arbitrage-free conditions on a discrete implied volatility grid
When the VIX — the implied volatility index for S&P 500 options — spiked to 82 in October 2008 and to 66 in March 2020, it was not measuring how much the stock market had moved. It was measuring how much the options market expected it to move over the next 30 days. The gap between the two — between what volatility has been and what it is priced to be — is one of the most actively traded quantities in financial markets.
Volatility is the central parameter in every option pricing model, and getting it right is the central problem of options practice. Chapter 10 introduced it as a constant input to Black-Scholes. This chapter tears that assumption apart. Volatility is not constant: it clusters (high volatility today predicts high volatility tomorrow), it mean-reverts (extreme volatility eventually subsides), and it depends on option strike and maturity in ways that Black-Scholes cannot explain. The market does not price all strikes at the same implied volatility — it prices out-of-the-money puts at higher implied vols because the market has experienced crashes and prices jump risk accordingly. This volatility smile (or skew, for equities) is the fingerprint of model failure, and the entire chapter is dedicated to modelling it properly.
We progress from the simplest measures of realised historical volatility, through GARCH models for its dynamics, to the full implied volatility surface. We cover the SVI parametrisation, Dupire's local volatility model, Heston stochastic volatility, and the model-free variance swap replication argument that underlies the VIX.
13.1 Historical Volatility
Before modelling how volatility is priced in options markets, we need to measure how much an asset has actually moved. Several estimators exist, each making different use of available price data.
module Historical_vol = struct
let log_returns prices =
let n = Array.length prices in
Array.init (n - 1) (fun i -> log (prices.(i + 1) /. prices.(i)))
let close_to_close ?(annualise = 252) prices =
let rs = log_returns prices in
let n = float_of_int (Array.length rs) in
let mean = Array.fold_left (+.) 0.0 rs /. n in
let var = Array.fold_left (fun a r -> a +. (r -. mean) *. (r -. mean)) 0.0 rs
/. (n -. 1.0) in
sqrt (float_of_int annualise *. var)
(** Parkinson high-low estimator — uses intraday range, more efficient *)
let parkinson ?(annualise = 252) ~highs ~lows () =
let n = Array.length highs in
let k = 1.0 /. (4.0 *. float_of_int n *. log 2.0) in
let ss = Array.fold_left2 (fun a h l ->
let hl = log (h /. l) in a +. hl *. hl) 0.0 highs lows in
sqrt (float_of_int annualise *. k *. ss)
(** Garman-Klass estimator — uses OHLC *)
let garman_klass ?(annualise = 252) ~opens ~highs ~lows ~closes () =
let n = float_of_int (Array.length opens) in
let sum = ref 0.0 in
Array.iteri (fun i _ ->
let hl = log (highs.(i) /. lows.(i)) in
let co = log (closes.(i) /. opens.(i)) in
sum := !sum +. 0.5 *. hl *. hl -. (2.0 *. log 2.0 -. 1.0) *. co *. co
) opens;
sqrt (float_of_int annualise *. !sum /. n)
(** Rolling window volatility *)
let rolling_vol ?(annualise = 252) prices window =
let rs = log_returns prices in
let n = Array.length rs in
if n < window then [||]
else Array.init (n - window + 1) (fun start ->
let sub = Array.sub rs start window in
let mean = Array.fold_left (+.) 0.0 sub /. float_of_int window in
let var = Array.fold_left (fun a r -> a +. (r -. mean) *. (r -. mean)) 0.0 sub
/. float_of_int (window - 1) in
sqrt (float_of_int annualise *. var)
)
end
13.2 GARCH Volatility Models
GARCH(1,1) models the conditional variance as:
with the stationarity constraint $\alpha + \beta < 1$.
Figure 13.1 — Simulated returns and conditional volatility $\sigma_t$ for a GARCH(1,1) process. Notice how periods of large returns (positive or negative) correspond to spikes in the conditional volatility, capturing the volatility clustering phenomenon.
module Garch = struct
type params = {
omega : float; (* baseline variance *)
alpha : float; (* shock loading *)
beta : float; (* persistence *)
mu : float; (* mean return *)
}
let long_run_var p = p.omega /. (1.0 -. p.alpha -. p.beta)
let long_run_vol p = sqrt (252.0 *. long_run_var p)
(** Filter: compute conditional variances given returns and parameters *)
let filter params returns =
let n = Array.length returns in
let sigma2 = Array.make n (long_run_var params) in
for t = 1 to n - 1 do
let eps = returns.(t - 1) -. params.mu in
sigma2.(t) <- params.omega
+. params.alpha *. eps *. eps
+. params.beta *. sigma2.(t - 1)
done;
sigma2
(** Log-likelihood of GARCH(1,1) under Gaussian innovations *)
let log_likelihood params returns =
let sigma2 = filter params returns in
let n = Array.length returns in
let ll = ref 0.0 in
for t = 0 to n - 1 do
let eps = returns.(t) -. params.mu in
ll := !ll -. 0.5 *. (log (2.0 *. Float.pi) +. log sigma2.(t)
+. eps *. eps /. sigma2.(t))
done;
!ll
(** n-step variance forecast *)
let forecast_variance params ~n_steps ~current_sigma2 =
let lrv = long_run_var params in
let ab = params.alpha +. params.beta in
Array.init n_steps (fun h ->
lrv +. (current_sigma2 -. lrv) *. (ab ** float_of_int (h + 1))
)
(** Aggregate multi-period variance for option pricing *)
let integrated_variance params ~n_steps ~current_sigma2 =
let forecasts = forecast_variance params ~n_steps ~current_sigma2 in
Array.fold_left (+.) 0.0 forecasts
end
13.3 Implied Volatility Surface
The implied volatility $\sigma_{\text{impl}}(T, K)$ tells us the market's implied $\sigma$ for each expiry $T$ and strike $K$.
module Iv_surface = struct
type point = {
expiry : float; (* years to expiry *)
strike : float;
moneyness : float; (* log(K/F) / sqrt(T) = log-moneyness scaled *)
iv : float;
}
type t = {
spot : float;
rate : float;
div_yield : float;
points : point array;
}
let forward surface t = surface.spot *. exp ((surface.rate -. surface.div_yield) *. t)
let moneyness surface ~expiry ~strike =
let f = forward surface expiry in
log (strike /. f) /. sqrt expiry
(** Interpolate IV by bilinear interpolation on (expiry, moneyness) grid *)
let interpolate surface ~expiry ~strike =
let m = moneyness surface ~expiry ~strike in
(* Find surrounding points and bilinear interpolate *)
let pts = Array.to_list surface.points in
(* Simple nearest-expiry-slice approach *)
let same_expiry = List.filter (fun p -> Float.abs (p.expiry -. expiry) < 0.01) pts in
match same_expiry with
| [] -> None
| pts ->
let sorted = List.sort (fun a b -> compare a.moneyness b.moneyness) pts in
let rec find = function
| [] | [_] -> None
| a :: ((b :: _) as rest) ->
if a.moneyness <= m && m <= b.moneyness then
let w = (m -. a.moneyness) /. (b.moneyness -. a.moneyness) in
Some (a.iv *. (1.0 -. w) +. b.iv *. w)
else find rest
in find sorted
(** Volatility smile: for a given expiry, return [(moneyness, iv)] *)
let smile surface ~expiry =
let pts = Array.to_list surface.points in
List.filter_map (fun p ->
if Float.abs (p.expiry -. expiry) < 0.01 then
Some (p.moneyness, p.iv)
else None
) pts
|> List.sort compare
end
13.4 Volatility Smile and Skew
The Black-Scholes assumption of constant $\sigma$ is inconsistent with market prices. The smile (or skew) quantifies this:
- Equity skew: Out-of-the-money puts trade at higher IV than calls (negative skew)
- FX smile: Both wings are elevated (symmetric smile)
Figure 13.1 — A typical equity options volatility smile. Short-dated options (1 Month) exhibit the steepest skew as tail risk is magnified, while long-dated options (1 Year) progressively flatten.
(** Parametric smile: SVI (Stochastic Volatility Inspired) model
Total variance w(k) = a + b*(ρ*(k-m) + sqrt((k-m)² + σ²))
where k = log(K/F) *)
module Svi = struct
type params = {
a : float; (* level *)
b : float; (* slope / curvature *)
rho : float; (* correlation, ρ ∈ (-1,1) *)
m : float; (* at-the-money offset *)
sigma : float; (* ATM curvature *)
}
let total_variance p k =
let d = k -. p.m in
p.a +. p.b *. (p.rho *. d +. sqrt (d *. d +. p.sigma *. p.sigma))
let implied_vol p ~k ~t =
sqrt (total_variance p k /. t)
(** No-arbitrage condition: ∂w/∂k ≥ 0 at wings, convexity check *)
let is_arbitrage_free p k_grid =
Array.for_all (fun k ->
let w = total_variance p k in
let dk = 0.001 in
let _ = implied_vol p ~k ~t:1.0 in (* validate vol is real *)
w >= 0.0 && total_variance p (k +. dk) >= total_variance p (k -. dk)
) k_grid
end
13.5 Local Volatility (Dupire)
The Dupire local volatility model gives an exact fit to any arbitrage-free surface:
module Local_vol = struct
(** Compute Dupire local vol from call price surface C(T, K).
Uses finite differences on the surface. *)
let dupire ~rate ~div_yield ~call_price ~call_dT ~call_dK ~call_d2K
~strike =
let num = call_dT +. (rate -. div_yield) *. strike *. call_dK
+. div_yield *. call_price in
let den = 0.5 *. strike *. strike *. call_d2K in
if den <= 1e-10 then None
else
let var = num /. den in
if var <= 0.0 then None
else Some (sqrt var)
(** Numerical Dupire from a discrete call price surface *)
let dupire_surface ~rate ~div_yield ~strikes ~expiries ~call_prices =
let n_k = Array.length strikes and n_t = Array.length expiries in
let result = Array.make_matrix n_t n_k None in
for j = 1 to n_t - 2 do
for i = 1 to n_k - 2 do
let c = call_prices.(j).(i) in
let dT = (call_prices.(j + 1).(i) -. call_prices.(j - 1).(i))
/. (expiries.(j + 1) -. expiries.(j - 1)) in
let dK = (call_prices.(j).(i + 1) -. call_prices.(j).(i - 1))
/. (strikes.(i + 1) -. strikes.(i - 1)) in
let d2K = (call_prices.(j).(i + 1) -. 2.0 *. c +. call_prices.(j).(i - 1))
/. let dk = (strikes.(i + 1) -. strikes.(i - 1)) /. 2.0 in dk *. dk in
result.(j).(i) <- dupire ~rate ~div_yield ~call_price:c ~call_dT:dT
~call_dK:dK ~call_d2K:d2K ~strike:strikes.(i)
done
done;
result
end
13.6 Heston Stochastic Volatility Model
The Heston model:
where $\kappa$ = mean reversion speed, $\theta$ = long-run variance, $\xi$ = vol of vol, $\rho$ = correlation.
module Heston = struct
type params = {
kappa : float; (* mean reversion speed *)
theta : float; (* long-run variance *)
xi : float; (* vol of vol *)
rho : float; (* S-V correlation *)
v0 : float; (* initial variance *)
}
(** Feller condition for non-zero variance: 2κθ > ξ² *)
let feller_satisfied p = 2.0 *. p.kappa *. p.theta > p.xi *. p.xi
(** Characteristic function of log(S_T/S_0) under Heston.
Used for semi-analytic pricing via Fourier inversion. *)
let char_fun p ~rate ~tau ~u =
(* Complex arithmetic with pairs (re, im) *)
let cplx_mul (ar, ai) (br, bi) = (ar *. br -. ai *. bi, ar *. bi +. ai *. br) in
let cplx_add (ar, ai) (br, bi) = (ar +. br, ai +. bi) in
let cplx_sqrt (ar, ai) =
let r = sqrt (ar *. ar +. ai *. ai) in
let rm = sqrt r in
let theta2 = atan2 ai ar /. 2.0 in
(rm *. cos theta2, rm *. sin theta2)
in
let cplx_exp (ar, ai) = let ea = exp ar in (ea *. cos ai, ea *. sin ai) in
let iu = (0.0, u) in
let iu2 = cplx_mul iu iu in (* -u² *)
(* d = sqrt((κ - iρξu)² + ξ²(iu + u²)) *)
let xi_sq = p.xi *. p.xi in
let rho_xi_u = (p.kappa, -. p.rho *. p.xi *. u) in
let rho_xi_sq = cplx_mul rho_xi_u rho_xi_u in
let xi2_iu_u2 = cplx_mul (xi_sq, 0.0) (cplx_add iu iu2) in
let d = cplx_sqrt (cplx_add rho_xi_sq xi2_iu_u2) in
let g_num = cplx_add rho_xi_u (cplx_mul (-1.0, 0.0) d) in
let g_den = cplx_add rho_xi_u d in
let e_dt = cplx_exp (cplx_mul d (-. tau, 0.0)) in
let g = cplx_mul g_num (let (a, b) = g_den in (1.0 /. (a *. a +. b *. b),
-. b /. (a *. a +. b *. b))) in
(* Albrecher et al. (2007) formulation — numerically stable for large u *)
(* D = (κ - iρξu - d) / ξ² · (1 - e^{-dτ}) / (1 - g·e^{-dτ}) *)
let one_minus_g_edt =
let (ge_r, ge_i) = cplx_mul g e_dt in
(1.0 -. ge_r, -. ge_i)
in
let d_over_xi2 =
let xi2 = p.xi *. p.xi in
let (dr, di) = cplx_add rho_xi_u (cplx_mul (-1.0, 0.0) d) in
(dr /. xi2, di /. xi2)
in
let one_minus_edt =
let (er, ei) = e_dt in
(1.0 -. er, -. ei)
in
(* D = d_over_xi2 * (1 - e^{-dτ}) / (1 - g·e^{-dτ}) *)
let big_d =
let num = cplx_mul d_over_xi2 one_minus_edt in
let (dr, di) = one_minus_g_edt in
let denom_sq = dr *. dr +. di *. di in
let (nr, ni) = num in
((nr *. dr +. ni *. di) /. denom_sq, (ni *. dr -. nr *. di) /. denom_sq)
in
(* C = κ·θ/ξ² · [(κ - iρξu - d)τ - 2·ln((1 - g·e^{-dτ})/(1 - g))] *)
let log_ratio =
let (gr, gi) = g in
let one_minus_g = (1.0 -. gr, -. gi) in
let (nr, ni) = one_minus_g_edt and (dr, di) = one_minus_g in
let ratio_r = (nr *. dr +. ni *. di) /. (dr *. dr +. di *. di) in
let ratio_i = (ni *. dr -. nr *. di) /. (dr *. dr +. di *. di) in
(0.5 *. log (ratio_r *. ratio_r +. ratio_i *. ratio_i),
atan2 ratio_i ratio_r)
in
let big_c =
let kappa_theta_xi2 = p.kappa *. p.theta /. (p.xi *. p.xi) in
let (rr, ri) = rho_xi_u in
let (dr, di) = d in
let dt_term = ((rr -. dr) *. tau, (ri -. di) *. tau) in
let log2 = cplx_mul (2.0, 0.0) log_ratio in
let (cr, ci) = cplx_add dt_term (cplx_mul (-1.0, 0.0) log2) in
(kappa_theta_xi2 *. cr, kappa_theta_xi2 *. ci)
in
(* φ(u) = exp(C + D·v₀ + iuX₀) where X₀ = ln(S/K) + r·τ *)
let x0 = log (spot /. strike) +. rate *. tau in
let iu_x0 = (0.0, u *. x0) in
let exponent = cplx_add big_c (cplx_add (cplx_mul big_d (p.v0, 0.0)) iu_x0) in
cplx_exp exponent
(** Monte Carlo under Heston using Euler-Maruyama discretisation *)
let mc_price p ~spot ~strike ~rate ?div_yield:(q=0.0) ~tau ~n_steps ~n_paths
~option_type () =
let dt = tau /. float_of_int n_steps in
let sqrt_dt = sqrt dt in
let df = exp (-. rate *. tau) in
let payoffs = Array.init n_paths (fun _ ->
let s = ref spot and v = ref (Float.max 1e-8 p.v0) in
for _ = 0 to n_steps - 1 do
let z1 = Mc.std_normal () in
let z2 = Mc.std_normal () in
let wv = z1 in
let ws = p.rho *. z1 +. sqrt (1.0 -. p.rho *. p.rho) *. z2 in
let sqrt_v = sqrt (Float.max 0.0 !v) in
s := !s *. exp ((rate -. q -. 0.5 *. !v) *. dt +. sqrt_v *. sqrt_dt *. ws);
v := Float.max 0.0 (!v +. p.kappa *. (p.theta -. !v) *. dt
+. p.xi *. sqrt_v *. sqrt_dt *. wv)
done;
match option_type with
| `Call -> Float.max 0.0 (!s -. strike)
| `Put -> Float.max 0.0 (strike -. !s)
) in
df *. Array.fold_left (+.) 0.0 payoffs /. float_of_int n_paths
end
13.7 Variance Swaps and VIX
A variance swap pays the difference between realised variance and a fixed strike:
The fair strike can be replicated by a log contract portfolio:
The VIX index is essentially the square root of this integral over 30-day options:
let variance_swap_strike ~rate ~forward ~tau ~strikes ~call_mid ~put_mid =
(* Numerical integration of the replication formula *)
let dk = strikes.(1) -. strikes.(0) in
let integral = ref 0.0 in
Array.iteri (fun i k ->
let price = if k < forward then put_mid.(i) else call_mid.(i) in
let contribution = 2.0 /. tau *. price /. (k *. k) *. dk in
integral := !integral +. contribution
) strikes;
(* Subtract the squared return of the forward *)
let adj = (2.0 /. tau) *. (forward /. (rate *. tau) -. 1.0 -. log forward) in
!integral -. adj
let vix_approx = variance_swap_strike (* VIX ≈ sqrt(30-day var swap strike × 100) *)
13.8 Chapter Summary
Volatility is the most important and most complex input in options pricing. This chapter moved through three levels of complexity: measuring realised volatility from historical data, modelling its dynamics under GARCH, and understanding how it is reflected in the full implied volatility surface.
Historical volatility estimators differ in their use of price data. The close-to-close estimator is universal and unbiased but uses only a fraction of the available information. Parkinson's estimator uses intraday ranges and achieves roughly 5x the statistical efficiency. Garman-Klass extends this to full OHLC data. In practice, professionals use rolling windows of 21 and 63 trading days to capture short-term and medium-term volatility regimes, respectively.
GARCH(1,1) captures the two most important empirical facts about volatility: clustering (high vol periods are autocorrelated) and mean reversion (long-run variance $\bar{\sigma}^2 = \omega/(1-\alpha-\beta)$ acts as an attractor). The model's conditional variance $h_t$ updates daily based on the most recent squared return (the ARCH term) and the previous day's conditional variance (the GARCH term). Persistence $\alpha + \beta$ close to 1 implies slow mean reversion; post-crisis data typically yields high persistence.
The implied volatility surface is the most important summary of market option prices. Its shape reflects the market's beliefs about tail risk (the smile/skew across strikes) and the term structure of future volatility (the vol of vol across expirations). SVI provides a no-arbitrage parametrisation of the smile at a single expiry. Dupire's formula converts any arbitrage-free surface into a local volatility surface — a deterministic function $\sigma_{\text{loc}}(S, t)$ that reproduces all market prices exactly under standard GBM dynamics, but generates unrealistic forward smiles. The Heston model introduces stochastic volatility with mean reversion, which gives more realistic forward dynamics at the cost of slower calibration.
In Practice — Volatility Trading and the Vol Surface
A volatility trader at a hedge fund takes positions based on the difference between implied and realised volatility. The core trade is selling implied vol (selling options) when implied vol is high relative to expected realised vol, and buying implied vol (buying options) when it is cheap.
The VIX index (the CBOE Volatility Index) is the market's 30-day implied vol for the S&P 500, computed from the variance swap replication formula. Historically, VIX has averaged around 20% and has a strong mean-reverting tendency. When VIX spikes above 30–40 (as in 2008, 2020), it typically reverts within weeks. Volatility traders exploit this by selling variance swaps or straddles when VIX is elevated.
The vol surface in practice. A major equity index like the S&P 500 has options trading at 20+ strikes and 10+ expiries simultaneously. The implied vol surface must be arbitrage-free: no calendar spread arbitrage (longer-dated total variance must exceed shorter-dated), no butterfly arbitrage (the smile must be convex in strike), and no static arbitrage (put-call parity must hold). The SVI parametrisation is widely used because it satisfies these constraints by construction.
Rough volatility (Gatheral et al., 2018) is the current frontier: empirical studies show that the log-volatility process has a Hurst exponent $H \approx 0.1$, far below 0.5 (Brownian motion). This means volatility is much rougher than GBM — it has short-range anti-persistence. The rough Bergomi model captures this with a fractional Brownian motion driver and reproduces the observed power-law term structure of the ATM vol skew: $\text{skew}(T) \propto T^{H-1/2}$.
Exercises
13.1 ★ Compute rolling 21-day and 63-day historical volatility for a simulated GBM path. Plot both series and compare to the true $\sigma$.
13.2 ★★ Fit GARCH(1,1) to log returns from S&P 500 data (use an array of daily returns). Estimate parameters by grid search over $(\alpha, \beta)$ with $\omega = \bar{\sigma}^2(1-\alpha-\beta)$.
13.3 ★★ Calibrate the SVI model to a set of market implied volatilities at a single expiry. Minimise sum of squared IV errors subject to $b > 0$ and $|\rho| < 1$.
13.4 ★★★ Implement a Heston model calibration: for given market calls, minimise the sum of squared price errors over $(v_0, \kappa, \theta, \xi, \rho)$ using a simple Nelder-Mead search.