Making the model agree with the market
A pricing model has parameters, and the right values are not found in a textbook but in today's market. Calibration chooses the parameters so that the model's prices are as close as possible to observed ones, usually the option quotes of the day. The recipe has three parts. An objective that measures the misfit, almost always a sum of squared differences, possibly weighted so that liquid quotes count more. A search that moves the parameters to reduce it. And checks that the result is sensible: stable from day to day, inside plausible bounds and good on quotes that were not used in the fit.
# The objective: how far the model's prices are from the market's, as a (weighted) sum of squares
def sum_squared_error(model, market, weights=None):
if weights is None:
weights = [1.0] * len(model)
return sum(w * (m - y) ** 2 for w, m, y in zip(weights, model, market))
market = [4.20, 8.60, 14.10]
model = [4.00, 8.80, 14.20]
print(round(sum_squared_error(model, market), 4)) # 0.09 0.2^2 + 0.2^2 + 0.1^2
print(round(sum_squared_error(model, market, [1, 4, 1]), 4)) # 0.21 the middle quote counts four times as muchTwo search tools: finding a root and finding a minimum
Two tasks come up constantly. First, implied volatility: given a market price, find the volatility at which the model gives it, that is a root of model(sigma) - price = 0. Because the price rises steadily with volatility, bisection works: start with an interval where the difference changes sign, halve it, keep the half where the sign still changes, and repeat; each step doubles the accuracy and it cannot fail. Second, a minimum: for a function with a single valley, golden-section search shrinks the interval by the factor 0.618 each step while reusing one of the two function values, so it needs one evaluation per step, and it does not need derivatives. Libraries (scipy.optimize) are faster and handle many parameters, but these two show the ideas.
# Implied volatility by bisection: the root of (model price - market price) in sigma
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def bs_call(s, k, t, r, sigma):
d1 = (math.log(s / k) + (r + sigma ** 2 / 2) * t) / (sigma * math.sqrt(t))
return s * norm_cdf(d1) - k * math.exp(-r * t) * norm_cdf(d1 - sigma * math.sqrt(t))
def bisect_root(f, lo, hi, tol=1e-10):
f_lo = f(lo)
while hi - lo > tol:
mid = (lo + hi) / 2
f_mid = f(mid)
if f_lo * f_mid <= 0:
hi = mid
else:
lo, f_lo = mid, f_mid
return (lo + hi) / 2
market_price = bs_call(100.0, 105.0, 1.0, 0.02, 0.27) # a quote generated with a volatility of 27%
implied = bisect_root(lambda sigma: bs_call(100.0, 105.0, 1.0, 0.02, sigma) - market_price, 0.01, 2.0)
print(round(implied, 6)) # 0.27 recovered# Golden-section search: one function evaluation per step, the interval shrinks by 0.618 each time
import math
def golden_section_minimize(f, lo, hi, tol=1e-8):
ratio = (math.sqrt(5) - 1) / 2
a, b = lo, hi
c, d = b - ratio * (b - a), a + ratio * (b - a)
fc, fd = f(c), f(d)
while b - a > tol:
if fc < fd:
b, d, fd = d, c, fc
c = b - ratio * (b - a)
fc = f(c)
else:
a, c, fc = c, d, fd
d = a + ratio * (b - a)
fd = f(d)
return (a + b) / 2
evaluations = []
def valley(x):
evaluations.append(x)
return (x - 1.3) ** 2 + 0.5
print(round(golden_section_minimize(valley, 0.0, 5.0), 5), len(evaluations) < 60) # 1.3 True about 40 evaluations for 1e-8Calibrating a volatility, and doing it responsibly
Putting the pieces together, calibrating a single volatility to several option prices means minimising the sum of squared errors between price_fn(strike, sigma) and the market prices over sigma, with price_fn passed in so the same code works for any model. For a model with several parameters (Heston has five, SABR three after fixing beta) the search is multidimensional and can land in a local minimum: start from several initial guesses, bound the parameters, and prefer relative or implied-volatility errors to raw price errors so that cheap out-of-the-money options are not ignored. Then check that the fitted parameters are stable from one day to the next and that the model prices options it was not fitted to well. A fit that is perfect but unstable is overfitting.
# Calibrate one volatility to three option prices: minimise the sum of squared errors over sigma
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def bs_call(s, k, t, r, sigma):
d1 = (math.log(s / k) + (r + sigma ** 2 / 2) * t) / (sigma * math.sqrt(t))
return s * norm_cdf(d1) - k * math.exp(-r * t) * norm_cdf(d1 - sigma * math.sqrt(t))
def golden_section_minimize(f, lo, hi, tol=1e-8):
ratio = (math.sqrt(5) - 1) / 2
a, b = lo, hi
c, d = b - ratio * (b - a), a + ratio * (b - a)
fc, fd = f(c), f(d)
while b - a > tol:
if fc < fd:
b, d, fd = d, c, fc
c = b - ratio * (b - a)
fc = f(c)
else:
a, c, fc = c, d, fd
d = a + ratio * (b - a)
fd = f(d)
return (a + b) / 2
def calibrate_sigma(strikes, market_prices, price_fn, lo, hi):
def error(sigma):
return sum((price_fn(k, sigma) - p) ** 2 for k, p in zip(strikes, market_prices))
return golden_section_minimize(error, lo, hi)
strikes = [90.0, 100.0, 110.0]
quotes = [bs_call(100.0, k, 1.0, 0.02, 0.23) for k in strikes] # market quotes made with sigma = 23%
sigma = calibrate_sigma(strikes, quotes, lambda k, s: bs_call(100.0, k, 1.0, 0.02, s), 0.01, 1.5)
print(round(sigma, 4)) # 0.23- Calibration picks model parameters that minimise a (weighted) sum of squared price errors; check stability across days and accuracy on quotes not used in the fit.
- Bisection finds a root (an implied volatility) by halving an interval with a sign change; golden section finds the minimum of a single-valley function with one evaluation per step.
- Pass the pricing function in, bound the parameters, start from several guesses for multi-parameter models, and reject nonsense fits.
Next up: neural networks for forecasting, learning a prediction from data instead of writing the formula.
Lab · what's graded
- Measure the misfit with
sum_squared_error25 pts - Solve f(x) = 0 with
bisect_root25 pts - Find the lowest point with
golden_section_minimize25 pts - Fit the market with
calibrate_sigma25 pts