Running the formula backwards
The Black–Scholes formula turns a volatility into a price. In the market you have the opposite: the price is quoted and the volatility is unknown. The implied volatility is the sigma for which the formula reproduces the market price. There is no algebraic inverse, so you search: compute the price for a trial sigma, look at the price error (model minus market), and move sigma until the error is zero. The simplest search tries a grid of values and keeps the best.
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def bs_call(spot, strike, rate, sigma, t):
d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))
market = 8.0 # the quoted price of the 100-strike, 6-month call on a stock at 100
for sigma in (0.10, 0.15, 0.20, 0.25, 0.30):
model = bs_call(100, 100, 0.05, sigma, 0.5)
print(sigma, round(model, 4), "error", round(model - market, 4))Why there is exactly one answer
The call price rises steadily with sigma: vega is positive, so the error crosses zero only once. That makes the answer unique, and it tells you which way to move: if the model is too expensive, lower sigma. A solution only exists if the market price lies strictly between the *intrinsic floor* max(S - K * exp(-rT), 0) (volatility zero) and the spot S (volatility infinite). Outside those bounds the quote contradicts no-arbitrage, and the right answer is an error, not a number.
# The bounds: no volatility can explain a price outside them
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def bs_call(spot, strike, rate, sigma, t):
d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))
spot, strike, rate, t = 100, 90, 0.05, 1
floor = max(spot - strike * math.exp(-rate * t), 0)
print(round(floor, 4), spot) # 14.3894 100
print(round(bs_call(spot, strike, rate, 0.0001, t), 4)) # 14.3894 tiny volatility: the floor
print(round(bs_call(spot, strike, rate, 5.0, t), 4)) # 98.8513 huge volatility: close to the spotHalving the interval
A grid is crude: the answer is only as fine as the grid. Better use the sign of the error. Start with sigma between a very low and a very high value, take the midpoint, and look at the error there: positive means the midpoint is too high, so keep the lower half; negative, keep the upper. Every step halves the interval, so about thirty steps give eight digits. The next topic names and generalises this; here it is the tool that reads the market.
# Implied volatility by halving the interval, then reprice to check it
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def bs_call(spot, strike, rate, sigma, t):
d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))
def implied_vol(price, spot, strike, rate, t, tol=1e-8):
low, high = 0.0001, 5.0
while high - low > tol:
mid = (low + high) / 2
if bs_call(spot, strike, rate, mid, t) > price:
high = mid
else:
low = mid
return (low + high) / 2
iv = implied_vol(8.0, 100, 100, 0.05, 0.5)
print(round(iv, 6)) # 0.240534
print(round(bs_call(100, 100, 0.05, iv, 0.5), 6)) # 8.0The smile
If Black–Scholes were exactly right, every option on the same stock and date would imply the same volatility. They do not: low strikes usually imply more than high ones (the *skew*), or both wings more than the centre (the *smile*). The pattern shows what the model misses (jumps, fat tails, the demand for protection), which is why options are quoted in implied volatility: a single, comparable number across strikes and expiries, from which the price follows.
# A smile: market prices at five strikes, converted to implied volatilities
import math
def norm_cdf(x):
return 0.5 * (1 + math.erf(x / math.sqrt(2)))
def bs_call(spot, strike, rate, sigma, t):
d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))
def implied_vol(price, spot, strike, rate, t, tol=1e-8):
low, high = 0.0001, 5.0
while high - low > tol:
mid = (low + high) / 2
low, high = (low, mid) if bs_call(spot, strike, rate, mid, t) > price else (mid, high)
return (low + high) / 2
strikes = [90, 95, 100, 105, 110]
quotes = [bs_call(100, k, 0.05, v, 0.5) for k, v in zip(strikes, [0.28, 0.24, 0.21, 0.22, 0.25])] # prices from a made-up smile
for strike, quote in zip(strikes, quotes):
print(strike, round(quote, 4), round(implied_vol(quote, 100, strike, 0.05, 0.5), 4))- Implied volatility is the
sigmathat makesbs_callequal the market price: found by search, because there is no inverse formula. - The price rises with
sigma, so there is one answer and the sign of the error says which way to move; prices outside the bounds have none. - Halving an interval converges fast (about 30 steps for 8 digits); the smile is the implied volatility plotted against strike.
Next up: root finding with bisection and Newton, the general method behind this search and a much faster one.
Lab · what's graded
- Compare model and market with
call_price_error25 pts - Search a grid of volatilities with
implied_vol_scan25 pts - Invert the formula with
implied_vol25 pts - Read a whole smile with
implied_vol_smile25 pts