Chapter 13 — Volatility
"Volatility is the only thing in finance that mean-reverts reliably — everything else is wishful thinking."
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.
$$\hat{\sigma}{\text{close}} = \sqrt{\frac{252}{n-1} \sum{i=1}^n (r_i - \bar{r})^2}, \quad r_i = \ln(S_i / S_{i-1})$$
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:
$$\sigma_t^2 = \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma_{t-1}^2, \quad \epsilon_t = r_t - \mu$$
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:
$$\sigma_{\text{local}}^2(T, K) = \frac{\frac{\partial C}{\partial T} + (r-q) K \frac{\partial C}{\partial K} + q C }{\frac{1}{2} K^2 \frac{\partial^2 C}{\partial K^2}}$$
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:
$$dS_t = (r - q) S_t\cdot dt + \sqrt{V_t} S_t\cdot dW_t^S$$ $$dV_t = \kappa(\theta - V_t)\cdot dt + \xi \sqrt{V_t}\cdot dW_t^V, \quad \langle dW^S, dW^V \rangle = \rho\cdot dt$$
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
let _g = g in
(** Full formula left as exercise — real implementations use
Albrecher et al. (2007) formulation for numerical stability *)
let _ = (e_dt, g) in
(cos (u *. rate *. tau), sin (u *. rate *. tau)) (* placeholder *)
(** 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:
$$\text{Payoff} = N \cdot (\sigma_{\text{realised}}^2 - K_{\text{var}})$$
The fair strike can be replicated by a log contract portfolio:
$$K_{\text{var}} = \frac{2}{T}\int_0^\infty \frac{C(K) \mathbf{1}{K>F} + P(K) \mathbf{1}{K<F}}{K^2}\cdot dK$$
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.
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.