Option Pricing II: Black-Scholes, Greeks and Implied Volatility
Implement the Black-Scholes formula with std::erfc, verify every Greek by bump-and-reprice, and invert prices to implied volatility with bisection and safeguarded Newton.
From a thousand-step tree to one line of arithmetic
The binomial tree in the previous lesson priced a European call at 10.4486 with 1000 steps and 500,500 node updates, and was still 0.002 away from its limit. For European options on a stock without dividends, that limit has a closed form, published by Fischer Black and Myron Scholes in 1973 and extended by Robert Merton the same year. It evaluates in a few dozen floating-point operations.
Speed is only part of why it matters. A closed form can be differentiated, which gives the Greeks: how the price responds to the stock price, volatility, time and interest rates. Those sensitivities are what a risk system aggregates across thousands of positions. And the formula can be run backwards: given a price, find the volatility that produces it. That number, the implied volatility, is how option prices are commonly quoted and compared.
This lesson implements all three: the price, the five standard Greeks, and implied volatility. Each comes with a numerical check, because each has a way of being silently wrong: a Greek in the wrong units, a solver that diverges on one input in a thousand, a formula that returns NaN at expiry. The parameters are illustrative; nothing here describes any real instrument or is advice of any kind.
The formula and what each piece means
For a European call with spot price S, strike K, continuously compounded risk-free rate r, volatility sigma and time to expiry T in years, on a stock that pays no dividends:
d1 = (log(S / K) + (r + sigma^2 / 2) * T) / (sigma * sqrt(T)), d2 = d1 - sigma * sqrt(T)
C = S * N(d1) - K * exp(-r * T) * N(d2), P = K * exp(-r * T) * N(-d2) - S * N(-d1)
where N is the cumulative distribution function of the standard normal distribution. The formula is the risk-neutral expectation from the binomial lesson, carried to its continuous limit: under the risk-neutral probability the stock follows geometric Brownian motion with drift r, log(S_T) is normal with standard deviation sigma * sqrt(T), and C = exp(-r * T) * E[max(S_T - K, 0)] integrates to the expression above. Two pieces have a direct reading:
N(d2)is the risk-neutral probability that the call finishes in the money (S_T > K), soK * exp(-r * T) * N(d2)is the discounted expected payment of the strike.N(d1)is the call's delta, the number of shares in the replicating portfolio, as in the one-step example of the previous lesson.
N is not in the C++ standard library under that name, but std::erfc, the complementary error function, is (in cmath since C++11), and N(x) = 0.5 * erfc(-x / sqrt(2)). Use erfc, not 0.5 * (1 + erf(x / sqrt(2))): for very negative x, erf is close to -1 and the sum cancels to rounding noise, while erfc keeps full relative precision in the tail. Deep out-of-the-money options live in that tail.
Dividends. For a stock with a continuous dividend yield q, Merton's extension replaces S by S * exp(-q * T) in the price and r by r - q in d1. This lesson uses q = 0 throughout.
The first program implements the formula and reproduces a published example: Example 15.6 in John Hull's "Options, Futures, and Other Derivatives" (9th edition, 2014), a six-month option with S = 42, K = 40, r = 10%, sigma = 20%, for which Hull reports a call worth 4.76 and a put worth 0.81 (with d1 and d2 of about 0.7693 and 0.6278). The second case is the at-the-money option the binomial tree converged towards.
#include <cmath>
#include <numbers>
#include <print>
// Standard normal cumulative distribution function.
// erfc is accurate in the far left tail, where 1 + erf(x) would cancel.
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
struct Inputs {
double spot, strike, rate, vol, expiry; // rate: cont. compounded; vol, expiry: years
};
struct Prices { double d1, d2, call, put; };
Prices black_scholes(const Inputs& in) {
const double sd = in.vol * std::sqrt(in.expiry); // std. dev. of log(S_T)
const double drift = (in.rate + 0.5 * in.vol * in.vol) * in.expiry;
const double d1 = (std::log(in.spot / in.strike) + drift) / sd;
const double d2 = d1 - sd;
const double disc_k = in.strike * std::exp(-in.rate * in.expiry); // K e^{-rT}
return {d1, d2, in.spot * norm_cdf(d1) - disc_k * norm_cdf(d2),
disc_k * norm_cdf(-d2) - in.spot * norm_cdf(-d1)};
}
void show(const char* label, const Inputs& in) {
Prices p = black_scholes(in);
std::println("{}: S={} K={} r={} vol={} T={}", label, in.spot, in.strike, in.rate,
in.vol, in.expiry);
std::println(" d1 = {:.4f} d2 = {:.4f}", p.d1, p.d2);
std::println(" call = {:.4f} put = {:.4f}", p.call, p.put);
std::println(" parity: C - P = {:.6f}, S - K e^-rT = {:.6f}", p.call - p.put,
in.spot - in.strike * std::exp(-in.rate * in.expiry));
}
int main() {
show("Hull, Example 15.6", {42, 40, 0.10, 0.20, 0.5});
show("tree test case ", {100, 100, 0.05, 0.20, 1.0});
}Press Run to compile and run this example.
Hull, Example 15.6: S=42 K=40 r=0.1 vol=0.2 T=0.5
d1 = 0.7693 d2 = 0.6278
call = 4.7594 put = 0.8086
parity: C - P = 3.950823, S - K e^-rT = 3.950823
tree test case : S=100 K=100 r=0.05 vol=0.2 T=1
d1 = 0.3500 d2 = 0.1500
call = 10.4506 put = 5.5735
parity: C - P = 4.877058, S - K e^-rT = 4.877058Every number matches the published example, and the second case gives 10.4506, the value the binomial tree was approaching. Independently of both, the program checks put-call parity, C - P = S - K * exp(-r * T), which holds to every printed digit. Three different kinds of evidence (a textbook value, a second pricing method, a model-free identity) agree; that is what "the formula is implemented correctly" should mean. (The reference values were also recomputed for this lesson by numerical integration of the lognormal payoff, independently of the closed form.)
Seeing the price as a curve
A single price hides the shape that the Greeks describe. This program draws the call price against the spot price for two expiries, with K = 50, r = 5%, sigma = 30%, along with the payoff at expiry. SVG's y-axis points down, so the program maps higher values to smaller y coordinates.
#include <algorithm>
#include <cmath>
#include <format>
#include <numbers>
#include <print>
#include <string>
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
double bs_call(double s, double k, double r, double vol, double t) {
const double sd = vol * std::sqrt(t);
const double d1 = (std::log(s / k) + (r + 0.5 * vol * vol) * t) / sd;
return s * norm_cdf(d1) - k * std::exp(-r * t) * norm_cdf(d1 - sd);
}
int main() {
const double k = 50, r = 0.05, vol = 0.30;
const double s_min = 25, s_max = 75, v_max = 28;
// Map spot to x and value to y; SVG's y axis points down, so subtract.
auto x = [&](double s) { return 30 + 360 * (s - s_min) / (s_max - s_min); };
auto y = [&](double v) { return 190 - 180 * v / v_max; };
auto curve = [&](auto value) {
std::string pts;
for (double s = s_min; s <= s_max + 1e-9; s += 0.5)
pts += std::format("{:.1f},{:.1f} ", x(s), y(value(s)));
return pts;
};
std::println(R"(<svg viewBox="0 0 400 220" font-size="9" fill="#1f2937">)");
std::println(R"(<line x1="30" y1="190" x2="390" y2="190" stroke="#1f2937"/>)");
std::println(R"(<line x1="30" y1="10" x2="30" y2="190" stroke="#1f2937"/>)");
std::println(R"(<polyline points="{}" fill="none" stroke="#1f2937" )"
R"(stroke-width="1.5"/>)",
curve([&](double s) { return std::max(s - k, 0.0); }));
std::println(R"(<polyline points="{}" fill="none" stroke="#2563eb" )"
R"(stroke-dasharray="4 3"/>)",
curve([&](double s) { return bs_call(s, k, r, vol, 0.1); }));
std::println(R"(<polyline points="{}" fill="none" stroke="#2563eb" )"
R"(stroke-width="1.5"/>)",
curve([&](double s) { return bs_call(s, k, r, vol, 1.0); }));
std::println(R"(<text x="{:.1f}" y="205">K = 50</text>)", x(k) - 14);
std::println(R"(<text x="{:.1f}" y="205">spot 25</text>)", x(25) - 4);
std::println(R"(<text x="{:.1f}" y="205">75</text>)", x(75) - 10);
std::println(R"(<text x="40" y="20">payoff max(S - K, 0)</text>)");
std::println(R"(<text x="40" y="32" fill="#2563eb">T = 1 year (solid), )"
R"(T = 0.1 year (dashed)</text>)");
std::println("</svg>");
}Press Run to compile and run this example.
The black line is the payoff max(S - K, 0), with its kink at the strike. Both blue curves lie above it: an option before expiry is worth at least its exercise value, plus the value of the chance that things improve. The one-year curve (solid) is higher and smoother than the 0.1-year curve (dashed), which hugs the payoff except near the strike. The slope of each curve is delta, rising from 0 on the left to 1 on the right; the curvature is gamma, concentrated around the strike and sharper for the shorter expiry; the gap between the two curves is what the passage of time takes away, which is theta.
The Greeks in closed form
Each Greek is a partial derivative of the price with respect to one input, with the others held fixed. For a European call without dividends, with n(x) = exp(-x^2 / 2) / sqrt(2 * pi) the standard normal density:
| Greek | Derivative | Call | Put |
|---|---|---|---|
| delta | dV/dS | N(d1) | N(d1) - 1 |
| gamma | d2V/dS2 | n(d1) / (S * sigma * sqrt(T)) | same as call |
| vega | dV/dsigma | S * n(d1) * sqrt(T) | same as call |
| theta | dV/dt | -S * n(d1) * sigma / (2 * sqrt(T)) - r * K * exp(-r * T) * N(d2) | -S * n(d1) * sigma / (2 * sqrt(T)) + r * K * exp(-r * T) * N(-d2) |
| rho | dV/dr | K * T * exp(-r * T) * N(d2) | -K * T * exp(-r * T) * N(-d2) |
The put's gamma and vega equal the call's because, by put-call parity, C - P is linear in S and does not depend on sigma.
Units cause more bugs than formulas. Theta here is the derivative with respect to calendar time t, which shortens T, in units of value per year; to report it per calendar day divide by 365, per trading day by 252, and say which. Vega is per 1.00 of volatility (from 20% to 120%); the usual quote "per vol point" (20% to 21%) is vega divided by 100. Rho is per 1.00 of rate, so per 1% it is rho divided by 100.
Bump and reprice. The most direct check of a Greek is to move the input a little and reprice. The central difference (V(x + h) - V(x - h)) / (2h) approximates the first derivative with an error proportional to h^2, much better than the one-sided (V(x + h) - V(x)) / h, whose error is proportional to h. For gamma, the second difference (V(S + h) - 2V(S) + V(S - h)) / h^2 also has an h^2 error. The program uses the parameters of the delta-hedging example in Hull's chapter "The Greek Letters": S = 49, K = 50, r = 5%, sigma = 20%, 20 weeks to expiry. Hull values the call at 2.40 (as $240,000 for 100,000 options) with a delta of 0.522.
#include <cmath>
#include <numbers>
#include <print>
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
double norm_pdf(double x) {
return std::exp(-0.5 * x * x) / std::sqrt(2.0 * std::numbers::pi);
}
struct Inputs { double spot, strike, rate, vol, expiry; };
double bs_call(const Inputs& in) {
const double sd = in.vol * std::sqrt(in.expiry);
const double drift = (in.rate + 0.5 * in.vol * in.vol) * in.expiry;
const double d1 = (std::log(in.spot / in.strike) + drift) / sd;
return in.spot * norm_cdf(d1) -
in.strike * std::exp(-in.rate * in.expiry) * norm_cdf(d1 - sd);
}
struct Greeks { double delta, gamma, vega, theta, rho; };
// Closed-form Greeks of a European call (no dividends). Units: per 1.00 change in
// spot, vol (i.e. 100 vol points) and rate; theta per year of calendar time.
Greeks call_greeks(const Inputs& in) {
const double sq = std::sqrt(in.expiry), sd = in.vol * sq;
const double drift = (in.rate + 0.5 * in.vol * in.vol) * in.expiry;
const double d1 = (std::log(in.spot / in.strike) + drift) / sd, d2 = d1 - sd;
const double disc_k = in.strike * std::exp(-in.rate * in.expiry);
return {norm_cdf(d1),
norm_pdf(d1) / (in.spot * sd),
in.spot * norm_pdf(d1) * sq,
-in.spot * norm_pdf(d1) * in.vol / (2 * sq) - in.rate * disc_k * norm_cdf(d2),
in.expiry * disc_k * norm_cdf(d2)};
}
// Central differences: bump one input up and down by h and reprice.
Greeks bumped_greeks(Inputs in) {
auto central = [&](double Inputs::*field, double h) {
Inputs up = in, down = in;
up.*field += h;
down.*field -= h;
return (bs_call(up) - bs_call(down)) / (2 * h);
};
const double hs = 0.001 * in.spot; // 0.1% of spot
Inputs up = in, down = in;
up.spot += hs;
down.spot -= hs;
return {central(&Inputs::spot, hs),
(bs_call(up) - 2 * bs_call(in) + bs_call(down)) / (hs * hs),
central(&Inputs::vol, 1e-4),
-central(&Inputs::expiry, 1e-4), // theta: time passing shortens T
central(&Inputs::rate, 1e-4)};
}
int main() {
const Inputs in{49, 50, 0.05, 0.20, 20.0 / 52}; // 20 weeks to expiry
const Greeks a = call_greeks(in), b = bumped_greeks(in);
std::println("call price: {:.4f}", bs_call(in));
std::println("{:<6} {:>11} {:>11}", "greek", "closed form", "bumped");
std::println("{:<6} {:>11.4f} {:>11.4f}", "delta", a.delta, b.delta);
std::println("{:<6} {:>11.5f} {:>11.5f}", "gamma", a.gamma, b.gamma);
std::println("{:<6} {:>11.4f} {:>11.4f}", "vega", a.vega, b.vega);
std::println("{:<6} {:>11.4f} {:>11.4f}", "theta", a.theta, b.theta);
std::println("{:<6} {:>11.4f} {:>11.4f}", "rho", a.rho, b.rho);
}Press Run to compile and run this example.
call price: 2.4005
greek closed form bumped
delta 0.5216 0.5216
gamma 0.06554 0.06554
vega 12.1055 12.1055
theta -4.3053 -4.3053
rho 8.9070 8.9070The price, 2.4005, and delta, 0.5216, match Hull's 2.40 and 0.522. For every Greek, the closed form and the bumped reprice agree to the printed precision. Theta of -4.3053 per year is about -0.0118 per calendar day. Vega of 12.1055 means a one-point rise in volatility (20% to 21%) adds about 0.121 to a price of 2.40.
The bump-and-reprice function is worth keeping even after the closed forms are trusted. It needs nothing but the pricer, so it works for any model, including ones where no closed-form Greeks exist, such as the American options of the previous lesson or the Monte Carlo prices of the next.
Choosing the bump size
The bump h cannot simply be made as small as possible. Two errors pull in opposite directions:
- Truncation error comes from the Taylor expansion the difference formula is built on. For central differences it is about
c * h^2for some constantcthat depends on the higher derivatives; it shrinks ashshrinks. - Rounding error comes from the prices themselves. Each price carries a relative error around machine epsilon (
2.2e-16fordouble), and the second difference divides that byh^2. It grows ashshrinks.
The program estimates the constant c for gamma from one bump where rounding is negligible, then compares the predicted truncation error with the rounding scale eps * C / h^2 for bumps from 10% of spot down to 1e-8 of spot.
#include <cmath>
#include <limits>
#include <numbers>
#include <print>
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
double norm_pdf(double x) {
return std::exp(-0.5 * x * x) / std::sqrt(2.0 * std::numbers::pi);
}
const double kStrike = 50, kRate = 0.05, kVol = 0.20, kExpiry = 20.0 / 52;
double d1_of(double spot) {
const double drift = (kRate + 0.5 * kVol * kVol) * kExpiry;
return (std::log(spot / kStrike) + drift) / (kVol * std::sqrt(kExpiry));
}
double call(double spot) {
const double d1 = d1_of(spot), d2 = d1 - kVol * std::sqrt(kExpiry);
return spot * norm_cdf(d1) - kStrike * std::exp(-kRate * kExpiry) * norm_cdf(d2);
}
double gamma_exact(double spot) {
return norm_pdf(d1_of(spot)) / (spot * kVol * std::sqrt(kExpiry));
}
int main() {
const double spot = 49.0, eps = std::numeric_limits<double>::epsilon();
auto second_difference = [&](double h) {
return (call(spot + h) - 2 * call(spot) + call(spot - h)) / (h * h);
};
// Truncation error of the central second difference is about c * h^2.
// Measure c once, at a bump where rounding is negligible.
const double h_ref = spot * 1e-2;
const double c =
std::abs(second_difference(h_ref) - gamma_exact(spot)) / (h_ref * h_ref);
std::println("{:>6} {:>11} {:>14} {}", "h / S", "truncation", "rounding scale",
"measured error");
for (int k = 1; k <= 8; ++k) {
const double h = spot * std::pow(10.0, -k);
const double truncation = c * h * h;
// Each price carries a relative error near eps; the second difference
// divides that by h^2. The exact noise varies from machine to machine.
const double rounding = eps * call(spot) / (h * h);
const double error = std::abs(second_difference(h) - gamma_exact(spot));
if (truncation > 1000 * rounding)
std::println("{:>6.0e} {:>11.0e} {:>14.0e} {:.0e}", std::pow(10.0, -k),
truncation, rounding, error);
else
std::println("{:>6.0e} {:>11.0e} {:>14.0e} noise (varies by machine)",
std::pow(10.0, -k), truncation, rounding);
}
}Press Run to compile and run this example.
h / S truncation rounding scale measured error
1e-01 3e-03 2e-17 3e-03
1e-02 3e-05 2e-15 3e-05
1e-03 3e-07 2e-13 3e-07
1e-04 3e-09 2e-11 noise (varies by machine)
1e-05 3e-11 2e-09 noise (varies by machine)
1e-06 3e-13 2e-07 noise (varies by machine)
1e-07 3e-15 2e-05 noise (varies by machine)
1e-08 3e-17 2e-03 noise (varies by machine)Down to h / S = 1e-3, truncation dominates by a wide margin, and the measured error falls by a factor of 100 for every factor of 10 in h, the h^2 law, landing at 3e-7. At 1e-4 truncation (3e-9) is only about 150 times the rounding scale (2e-11), so the last digits of the measurement already depend on the machine. From 1e-5 on, the rounding scale exceeds the truncation error and grows by 100 per row: at h / S = 1e-8 it is 2e-3, a 3% error in a gamma of 0.066. The measured error in those rows is noise, and its exact value depends on the last bits of exp, log and erfc, which differ between math libraries; that is why the program refuses to print it. (The "Try it" section shows how to see the noise yourself.) The best h sits where the two curves cross, around h / S = 3e-5 here.
The general rule of thumb: for a central first difference the best relative bump is around eps^(1/3), about 6e-6; for a second difference around eps^(1/4), about 1e-4. Choose bumps relative to the size of the input, and prefer bumps that are exact in binary when you can (powers of two), so that x + h and x - h are exactly 2h apart.
Implied volatility: running the formula backwards
Every input to the formula is observable except sigma. The implied volatility of an option is the sigma at which the formula reproduces the option's market price. Because traders compare options with different strikes and expiries, they commonly quote volatility rather than price: "20 vol" carries the same information as 4.7594 for the Hull example, but is comparable across contracts.
Inverting is well posed. Vega is positive, so the call price is strictly increasing in sigma. As sigma goes to zero the price falls to max(S - K * exp(-r * T), 0), and as sigma grows without bound it rises to S. So a solution exists and is unique exactly when the price lies strictly between those two bounds, and a price outside them is not a volatility problem: it is an arbitrage or bad data, and the solver must reject it rather than return something.
Two solvers, with opposite strengths:
- Bisection keeps an interval
[lo, hi]known to contain the answer and halves it each step. With the bracket[1e-6, 5](from essentially zero to 500% volatility) and a tolerance of1e-10, it always takes 36 steps, becauselog2(5 / 1e-10)is about 35.5. It cannot fail on a valid price. - Newton-Raphson steps
sigma = sigma - (C(sigma) - price) / vega(sigma). Near the answer the number of correct digits roughly doubles per step. Far from it, where vega is tiny, one step can throwsigmaanywhere.
A safeguarded Newton combines them: keep the bracket, take the Newton step when it lands inside the bracket, bisect when it does not. It also starts from a good point: Manaster and Koehler (1982) proposed sigma0 = sqrt(2 * |log(S / K) + r * T| / T), the volatility at which the price, as a function of volatility, changes from convex to concave.
#include <algorithm>
#include <cmath>
#include <expected>
#include <numbers>
#include <print>
#include <string>
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
double norm_pdf(double x) {
return std::exp(-0.5 * x * x) / std::sqrt(2.0 * std::numbers::pi);
}
struct Inputs { double spot, strike, rate, expiry; };
struct PriceVega { double price, vega; };
PriceVega call(const Inputs& in, double vol) {
const double sq = std::sqrt(in.expiry), sd = vol * sq;
const double drift = (in.rate + 0.5 * vol * vol) * in.expiry;
const double d1 = (std::log(in.spot / in.strike) + drift) / sd;
const double disc_k = in.strike * std::exp(-in.rate * in.expiry);
return {in.spot * norm_cdf(d1) - disc_k * norm_cdf(d1 - sd),
in.spot * norm_pdf(d1) * sq};
}
struct Solved { double vol; int iterations; };
using Result = std::expected<Solved, std::string>;
// A call price must lie strictly between its no-arbitrage bounds to have a volatility.
bool in_bounds(const Inputs& in, double price) {
double lower = std::max(in.spot - in.strike * std::exp(-in.rate * in.expiry), 0.0);
return price > lower && price < in.spot;
}
Result bisection(const Inputs& in, double price) {
if (!in_bounds(in, price)) return std::unexpected("price outside arbitrage bounds");
double lo = 1e-6, hi = 5.0; // stated bracket: 0.0001% to 500% volatility
int it = 0;
while (hi - lo > 1e-10) {
double mid = 0.5 * (lo + hi);
(call(in, mid).price > price ? hi : lo) = mid; // the price rises with vol
++it;
}
return Solved{0.5 * (lo + hi), it};
}
Result newton(const Inputs& in, double price, double vol) {
for (int it = 1; it <= 50; ++it) {
auto [p, vega] = call(in, vol);
vol -= (p - price) / vega;
if (!std::isfinite(vol) || vol <= 0 || vol > 100)
return std::unexpected(std::format("diverged at iteration {}, vol = {:.1e}",
it, vol));
if (std::abs(p - price) < 1e-12) return Solved{vol, it};
}
return std::unexpected("no convergence in 50 iterations");
}
// Newton steps inside a shrinking bracket; bisect whenever Newton would leave it.
Result safeguarded(const Inputs& in, double price) {
if (!in_bounds(in, price)) return std::unexpected("price outside arbitrage bounds");
double lo = 1e-6, hi = 5.0;
// Manaster and Koehler's starting point: sqrt(2 |log(S/K) + rT| / T).
double m = std::log(in.spot / in.strike) + in.rate * in.expiry;
double vol = std::clamp(std::sqrt(2 * std::abs(m) / in.expiry), lo, hi);
for (int it = 1; it <= 100; ++it) {
auto [p, vega] = call(in, vol);
if (std::abs(p - price) < 1e-12) return Solved{vol, it};
(p > price ? hi : lo) = vol;
double step = vol - (p - price) / vega;
vol = (step > lo && step < hi) ? step : 0.5 * (lo + hi);
}
return std::unexpected("no convergence in 100 iterations");
}
void report(const char* name, const Result& r) {
if (r) std::println(" {:<12} vol = {:.6f} after {:>2} iterations", name, r->vol,
r->iterations);
else std::println(" {:<12} failed: {}", name, r.error());
}
int main() {
const Inputs hull{42, 40, 0.10, 0.5}, far_otm{100, 200, 0.05, 0.25};
struct Case { const char* label; Inputs in; double price; };
const Case cases[] = {
{"Hull 15.6 inputs, true vol 30%", hull, call(hull, 0.30).price},
{"deep out of the money, true vol 60%", far_otm, call(far_otm, 0.60).price},
{"call quoted above the stock price", hull, 43.0},
};
for (const auto& c : cases) {
std::println("{} (price {:.6f})", c.label, c.price);
report("bisection", bisection(c.in, c.price));
report("Newton, 0.2", newton(c.in, c.price, 0.2));
report("safeguarded", safeguarded(c.in, c.price));
}
}Press Run to compile and run this example.
Hull 15.6 inputs, true vol 30% (price 5.714711)
bisection vol = 0.300000 after 36 iterations
Newton, 0.2 vol = 0.300000 after 5 iterations
safeguarded vol = 0.300000 after 5 iterations
deep out of the money, true vol 60% (price 0.167532)
bisection vol = 0.600000 after 36 iterations
Newton, 0.2 failed: diverged at iteration 1, vol = 6.9e+07
safeguarded vol = 0.600000 after 9 iterations
call quoted above the stock price (price 43.000000)
bisection failed: price outside arbitrage bounds
Newton, 0.2 failed: diverged at iteration 5, vol = 1.8e+06
safeguarded failed: price outside arbitrage boundsOn the Hull inputs, all three methods recover 30%; Newton needs 5 steps where bisection needs 36. On the deep out-of-the-money call (spot 100, strike 200, three months), plain Newton started at 20% fails on its first step. At 20% volatility the option is worth about 1e-11 (essentially nothing), vega is about 2.4e-9, and the step 0.17 / 2.4e-9 sends sigma to 69 million. Bisection is unaffected, and the safeguarded version starts at the Manaster-Koehler point (about 233%), stays inside its bracket and converges in 9 steps. The last case is a call quoted at 43 on a stock trading at 42. No volatility produces that price; plain Newton wanders off, while the two bracketed solvers reject it before iterating.
One more property of implied volatility deserves a warning. Deep in the money, the price barely depends on sigma (vega is tiny relative to the price), so a small error in the price becomes a large error in the implied volatility. Practitioners therefore usually back out implied volatility from the out-of-the-money option at each strike, the put below the forward and the call above it, using put-call parity to translate where needed.
The volatility smile. If Black-Scholes described markets exactly, every option on the same stock would have the same implied volatility. Plotting implied volatility against strike for one expiry typically gives a curve instead, called a smile or skew depending on its shape, and the curve changes with expiry. That is direct evidence against the assumption of a constant volatility and lognormal prices. Black-Scholes survives in practice as a quoting convention and a local model: prices are quoted in implied volatility, and models that fit the whole curve (local and stochastic volatility models) are calibrated to those quotes.