5D Library › Credit Cycles, Self-Organised Criticality and the Statistical Physics of Market Phase Transitions · 中文
Chapter 5 Implied Volatility, Skew and the Fear Index
The Black–Scholes formula and implied volatility
In 1973 Black and Scholes gave a pricing formula for European options (Black and Scholes 1973). It assumes that the price of the underlying asset follows geometric Brownian motion with constant volatility and can be traded continuously without frictions. The price of a call is
where is the price of the underlying, the strike, the time to expiry, the risk-free rate and the standard normal distribution function. Only in the formula cannot be observed directly. Conversely, given a market price for an option, the that makes the formula’s price equal the market price is called the implied volatility. It has no closed form but can be found by Newton’s method: each step corrects using the derivative of the option price with respect to (called vega).
If the Black–Scholes assumptions held, options with the same expiry and different strikes would give the same implied volatility. In practice they do not.
Evidence
The skew after 1987. Rubinstein reports that after the stock-market crash of October 1987 the implied volatilities of S&P 500 index options took on a marked and persistent shape: the lower the strike (the further out of the money the put), the higher the implied volatility (Rubinstein 1994). The shape is usually called the volatility “skew” or “smile”.
How jumps produce a skew
One way to produce a skew is to let the price process make occasional downward jumps in addition to its continuous fluctuations. In 1976 Merton gave option prices under such a “jump-diffusion” process: a sum of Black–Scholes prices weighted by the Poisson probabilities of the number of jumps (Merton 1976).
Calculation
Underlying at 100, three months to expiry, volatility of the continuous part 15 %; in addition, on average one jump every two years, with mean size about −10 % and standard deviation 10 %. First check Newton’s method on a known price, then, at each strike, take the jump-diffusion price of the out-of-the-money option (puts below 100, calls above) and back out its Black–Scholes implied volatility.
# Black-Scholes (1973) prices, implied volatility by Newton's method, and the skew that appears when
# the true process has downward jumps (Merton 1976 jump-diffusion, priced exactly as a Poisson-
# weighted sum of Black-Scholes prices). Parameters are illustrative, not fitted to any market.
import math
Ncdf = lambda x: 0.5 * math.erfc(-x / math.sqrt(2))
npdf = lambda x: math.exp(-x * x / 2) / math.sqrt(2 * math.pi)
def bs(S, K, T, r, s, call=True):
d1 = (math.log(S / K) + (r + s * s / 2) * T) / (s * math.sqrt(T)); d2 = d1 - s * math.sqrt(T)
if call: return S * Ncdf(d1) - K * math.exp(-r * T) * Ncdf(d2)
return K * math.exp(-r * T) * Ncdf(-d2) - S * Ncdf(-d1)
def vega(S, K, T, r, s):
d1 = (math.log(S / K) + (r + s * s / 2) * T) / (s * math.sqrt(T))
return S * npdf(d1) * math.sqrt(T)
def implied_vol(price, S, K, T, r, call=True, s=0.3):
for it in range(50):
diff = bs(S, K, T, r, s, call) - price
if abs(diff) < 1e-10: return s, it
s = max(1e-4, s - diff / vega(S, K, T, r, s)) # Newton step
raise ValueError("no convergence")
def merton(S, K, T, r, s, lam, muj, dj, call=True, nmax=40):
k = math.exp(muj + dj * dj / 2) - 1 # mean relative jump size
lp = lam * (1 + k)
tot = 0.0
for n in range(nmax):
sn = math.sqrt(s * s + n * dj * dj / T)
rn = r - lam * k + n * math.log(1 + k) / T
tot += math.exp(-lp * T) * (lp * T) ** n / math.factorial(n) * bs(S, K, T, rn, sn, call)
return tot
S, T, r = 100.0, 0.25, 0.03
p = bs(S, 95, T, r, 0.22, call=False)
iv, it = implied_vol(p, S, 95, T, r, call=False)
print(f"round trip: put K = 95 priced at sigma 22 % -> {p:.4f}; Newton recovers {iv:.6f} in {it} steps")
print("strike option jump-diffusion price Black-Scholes implied vol")
ivs = {}
for K in (80, 85, 90, 95, 100, 105, 110, 115, 120):
call = K >= S # out-of-the-money option at each strike
pr = merton(S, K, T, r, 0.15, lam=0.5, muj=-0.10, dj=0.10, call=call)
ivs[K], _ = implied_vol(pr, S, K, T, r, call)
print(f" {K:4d} {'call' if call else 'put ':4s} {pr:8.4f} {ivs[K]:6.2%}")
print(f"skew: implied vol at 90 minus at 110 = {ivs[90] - ivs[110]:+.2%}"
f" (diffusion sd alone is 15.00 %; with no jumps every strike would show 15.00 %)")
# round trip: put K = 95 priced at sigma 22 % -> 1.9997; Newton recovers 0.220000 in 3 steps
# strike option jump-diffusion price Black-Scholes implied vol
# 80 put 0.0958 23.41%
# 85 put 0.2347 21.41%
# 90 put 0.5592 19.46%
# 95 put 1.3503 18.00%
# 100 call 3.7864 17.11%
# 105 call 1.6830 16.61%
# 110 call 0.6137 16.33%
# 115 call 0.1856 16.19%
# 120 call 0.0483 16.19%
# skew: implied vol at 90 minus at 110 = +3.13% (diffusion sd alone is 15.00 %; with no jumps every strike would show 15.00 %)
Newton’s method recovers the 22 % implied volatility in three steps. In the jump-diffusion model the put with strike 80 has an implied volatility of about 23.4 % and calls with strikes of 110 and above about 16.2–16.3 %; the difference between 90 and 110 is about 3.1 percentage points. Implied volatility exceeds 15 % at every strike because the jumps add to total volatility, and it is much higher on the left because the jumps are mostly downward, which makes deep out-of-the-money puts worth more than in a purely continuous model.
Note
The skew is not “the probability of a crash”. Option prices are risk-neutral prices: they reflect both market participants’ judgements of how likely a fall is and the insurance premium they are willing to pay to avoid it, that is, a risk premium. A high implied volatility for out-of-the-money puts may arise because a fall is thought more likely, or because demand for protection is high and those willing to sell such insurance demand ample compensation; usually both. So the skew cannot be converted directly into “the market’s probability of a crash next month”. Nor are jumps the only mechanism that produces a skew: stochastic-volatility models, in which volatility is negatively correlated with the price, produce a similar shape.
The fear index and mean reversion
The VIX index of the Chicago Board Options Exchange is often called the “fear index”. Whaley explains its history: when launched in 1993 it used the implied volatilities of at-the-money options on the S&P 100 index; in 2003 it was changed to compute the expected volatility over the next 30 days from the prices of S&P 500 index options across a range of strikes, no longer relying on the Black–Scholes model (Whaley 2009). It too measures a risk-neutral expectation, and it tends to jump when stocks fall.
A notable feature of volatility indices is mean reversion: after a spike they tend to fall back, and after a long lull they tend to rise. A common simplified description lets the logarithm of the index follow an Ornstein–Uhlenbeck process (Uhlenbeck and Ornstein 1930):
The larger , the faster the reversion; the “half-life” of a deviation is .
Calculation
With a simulated series (not the real VIX) whose half-life is 30 trading days and long-run level 18, first compute the path back after a spike to 40 with no new shocks (which is also the median path in this model); then see how accurately the half-life can be estimated from 1 year and from 10 years of daily data; finally compute two common indicators, “IV rank” and the percentile.
# Mean reversion of a volatility index, modelled as an Ornstein-Uhlenbeck process for its logarithm:
# x = ln(V), dx = theta (mu - x) dt + eta dW (daily steps of 1/252 year)
# mu = ln 18, half-life 30 trading days, eta = 1.1 per sqrt(year). A simulated series, not the VIX.
import math, random, statistics as st
MU, HALF, ETA, DT = math.log(18), 30, 1.1, 1 / 252
TH = math.log(2) / (HALF * DT)
def path(n, rnd, x0=MU):
phi, sd = math.exp(-TH * DT), ETA * math.sqrt((1 - math.exp(-2 * TH * DT)) / (2 * TH))
x, out = x0, []
for _ in range(n):
x = MU + phi * (x - MU) + sd * rnd.gauss(0, 1); out.append(x)
return out
def half_life(x):
"""Regress x[t+1] on x[t]; half-life = -ln 2 / ln(slope), in days."""
a, b = x[:-1], x[1:]
ma, mb = st.fmean(a), st.fmean(b)
s = sum((u - ma) * (v - mb) for u, v in zip(a, b)) / sum((u - ma) ** 2 for u in a)
return -math.log(2) / math.log(s) if 0 < s < 1 else float("inf")
rnd = random.Random(4)
print("level after a spike to 40 with no new shocks (the median path), days:", ", ".join(
f"{h}: {math.exp(MU + (math.log(40) - MU) * math.exp(-TH * DT * h)):.1f}" for h in (0, 10, 30, 60, 120)))
for years in (1, 10):
est = sorted(half_life(path(252 * years, rnd)) for _ in range(200))
print(f"half-life estimated from {years:2d} year(s) of daily data, 200 samples: median {est[100]:.0f} days,"
f" middle 80 % {est[20]:.0f}-{est[180]:.0f} (true {HALF})")
v = [math.exp(u) for u in path(252, rnd)]
cur, lo, hi = v[-1], min(v), max(v)
print(f"one simulated year: last {cur:.1f}, range {lo:.1f}-{hi:.1f};"
f" IV-rank-style position {(cur - lo) / (hi - lo):.0%}, percentile {sum(u < cur for u in v) / len(v):.0%}")
# level after a spike to 40 with no new shocks (the median path), days: 0: 40.0, 10: 33.9, 30: 26.8, 60: 22.0, 120: 18.9
# half-life estimated from 1 year(s) of daily data, 200 samples: median 18 days, middle 80 % 9-37 (true 30)
# half-life estimated from 10 year(s) of daily data, 200 samples: median 29 days, middle 80 % 23-36 (true 30)
# one simulated year: last 26.8, range 13.4-34.1; IV-rank-style position 65%, percentile 85%In this model, after a spike to 40 and with no new shocks the index is back to about 27 after 30 trading days and about 19 after 120; this is also the median of the distribution, and because the log-normal distribution is skewed to the right the mean is slightly higher. But this path is not guaranteed: mean reversion describes an average tendency, and one spike can perfectly well be followed by a higher one.
The estimates deserve more attention. With 1 year of data the median estimated half-life is only about 18 days (the true value is 30), scattered between 9 and 37 days; only with 10 years of data does it come close to the true value. Short samples systematically overstate the speed of mean reversion, a consequence of the downward bias of autoregressive coefficients in small samples.
The last line shows that two common indicators can differ a great deal: the same reading sits at 65 % of the past year’s range (the IV-rank calculation: subtract the low and divide by the range) but is above 85 % of the past year’s days (the percentile). IV rank is determined by just two days, the highest and the lowest, and one extreme reading changes it; the percentile uses all the data.
Summary
Implied volatility is the volatility that makes the Black–Scholes formula match a market price, and it can be backed out with Newton’s method. After 1987 index options showed a persistent skew, and a price process with downward jumps can produce such a shape. The skew is a risk-neutral price that includes a risk premium and cannot be read directly as a crash probability. Volatility indices tend to revert to the mean, but short samples overstate the speed of reversion, and the commonly used IV rank and percentile can give very different readings. The next chapter discusses two more ambitious approaches: recognising the signals of a crisis before it comes.
References
Fan, Y. (2026). Credit Cycles, Self-Organised Criticality and the Statistical Physics of Market Phase Transitions: From Minsky's Financial Instability Hypothesis and Sandpile Power-Law Avalanches to Implied-Volatility Skew, Read Beside Classical Chinese Thought on Rise and Decline [信贷周期、自组织临界性与市场相变的统计物理] (Version 1.0.0). 5D Library. https://doi.org/10.5281/zenodo.23754711 · Licensed under Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0).