Black-Scholes, Greeks and Implied Volatility by Bisection in C++

An option trades at four point eight one three six. Nobody tells you the volatility. The market gives you only the price, and you must work backwards to the one number that produced it. That backwards step is called implied volatility, and it looks like a perfect job for calculus.

• Market gives the price
• We need sigma
• Invert Black-Scholes

Yet Kestrel does not use the famous fast method for this. It uses bisection, the slowest, plainest root finder there is. In this episode we read the production Black-Scholes code, handle its edge cases, and see exactly why the boring algorithm is the right one. Last time we priced trading costs. Today we price options.

• Newton is fast
• Bisection is safe
• Kestrel picks safe

We will go in order. First the formulas and Greeks. Then the normal distribution function, then the edge cases. Then the solver, with a full worked example using real numbers, and finally where this code is used and how it is tested. This is educational material, not financial advice, and every number here is illustrative.

• Formulas and Greeks
• Edge cases
• Bisection solver
• Tests

Black-Scholes starts with two numbers. d one equals the natural log of S over K, plus r plus sigma squared over two times T, all divided by sigma times the square root of T. Then d two is simply d one minus sigma root T. Here S is spot, K is strike, T is time in years, r is the rate.

• d1 = [ln(S/K) + (r + σ²/2)T] / (σ√T)
• d2 = d1 − σ√T

Picture the standard normal bell curve. The cumulative area to the left of d one, called N of d one, and to the left of d two, are the two probability-like weights in the price. N of d two is roughly how likely the option finishes in the money under the pricing measure. N of d one is the hedge ratio, the delta.

• Area left of d1 is N(d1)
• Area left of d2 is N(d2)
• d1 − d2 = σ√T

The call price is S times N of d one, minus K times e to the minus r T, times N of d two. The put is K e to the minus r T times N of minus d two, minus S times N of minus d one. Delta is N of d one for the call, and N of d one minus one for the put.

• Call: S·N(d1) − K·e^(−rT)·N(d2)
• Put: K·e^(−rT)·N(−d2) − S·N(−d1)
• Delta: N(d1), N(d1) − 1

Gamma is phi of d one divided by S sigma root T, where phi is the normal density. Vega is S times phi of d one times root T. Theta is quoted per year, and contains the decay term minus S phi of d one sigma over two root T. Gamma and vega are shared by calls and puts.

• gamma = φ(d1)/(Sσ√T)
• vega = S·φ(d1)·√T
• theta per year

A detail that matters later. Kestrel's vega is per one point zero zero of volatility, not per one percent. If volatility moves from zero point twenty to zero point twenty-one, the price changes by about vega times zero point zero one. Mixing these units is a classic bug, and we come back to it in the pitfalls.

• Vega is per 1.00 vol
• Not per 1%
• Scale by 0.01 yourself

Now the picture that powers everything. Fix S, K, T and r, and plot price against sigma. The curve starts at the discounted intrinsic value and rises smoothly, never dipping. It is monotone increasing. Its slope at each point is exactly vega. Remember this shape, because the entire solver rests on it.

• Price rises with sigma
• Never decreases
• Slope equals vega

Everything needs N of x. Kestrel computes normCdf of x as zero point five times erfc of minus x over root two. The complementary error function is in the C++ standard library, so there is no polynomial approximation to maintain and no table to trust.

• 0.5 · erfc(−x/√2)
• Standard library only
• Simplified listing

Why not write one minus something? Because of cancellation. For x equal to minus eight, N of x is about six times ten to the minus sixteen. Computing one minus a number very close to one throws away almost every significant digit. Erfc evaluates the small tail directly, so relative accuracy stays good.

• 1 − big number loses digits
• erfc keeps the tail
• Deep options need tails

Deep out-of-the-money options live in those tails. Their prices are tiny, and the implied volatility solver compares them to a target again and again. If the tail were noisy, the comparison would be noisy, and the root finder would chase rounding error instead of the true answer.

• Deep OTM prices are tiny
• Solver compares them repeatedly
• Noise misleads the search

Real inputs are messy. The first guard is time. Year fraction is seconds divided by three hundred sixty-five point two five times eighty-six thousand four hundred, clamped to be at least zero. That is calendar time, including weekends and holidays, not trading days.

• Seconds ÷ (365.25 × 86400)
• Calendar time
• Clamped to ≥ 0

If T is zero or negative, the formula would divide by zero. So expiredGreeks takes over and returns intrinsic value. A call is worth max of S minus K, a put max of K minus S. Delta is exactly zero or one, and the other Greeks are zero. At expiry there is nothing left to model.

• T ≤ 0: intrinsic value
• Delta exactly 0 or 1
• Other Greeks zero

Next, zero volatility, or a non-positive spot or strike. The logarithm and the division blow up. So deterministicGreeks returns the discounted forward intrinsic, max of S minus K e to the minus r T, and zero. With no randomness, the future is known, so the price is just the forward value discounted.

• Zero vol: no randomness
• S or K ≤ 0 handled
• max(S − K·e^(−rT), 0)

The design rule is simple: the pricer always returns a finite, sensible number. Derived bars from episode six call it on every bar, so one NaN would poison a whole backtest. Guards live at the top, the clean formula lives below them, and every path returns something defined.

• Guards first
• Formula second
• Always finite

Now to implied volatility. Newton's method is the obvious tool. Start with a guess, measure the error in price, and divide by the slope, which is vega, to get the next guess. Near a good starting point it converges astonishingly fast, doubling correct digits each step.

• Step = error ÷ vega
• Fast near the answer
• Quadratic convergence

But look at the curve for a deep out-of-the-money option. It is nearly flat. Vega is almost zero. Dividing the price error by a tiny vega produces a gigantic step, and the next guess lands at an absurd or even negative volatility. Newton shoots off the curve instead of climbing it.

• Deep OTM: vega near zero
• Division by tiny number
• Step flies off the curve

Bisection needs only one property: the price is monotone increasing in volatility. If the price at the low end is below the target, and the price at the high end is above it, a root is trapped between. Check the middle, keep the half that still brackets the target, repeat. It cannot diverge.

• Monotone price in vol
• Keep the bracketing half
• Cannot diverge

Bisection gains about one binary digit per step, so it needs more iterations. But each iteration is just one price evaluation, and the whole search is bounded. In a pipeline that runs unattended on thousands of option bars, a solver that always finishes beats one that is sometimes brilliant and sometimes catastrophic.

• One digit per step
• Bounded work
• Predictable behaviour

Here is the function. It first rejects the hopeless: T less than or equal to zero, S or K non-positive, or a price that is not finite. Each returns std nullopt, an honest answer meaning there is no implied volatility. It never returns NaN, because a NaN spreads silently through later arithmetic.

• Reject T ≤ 0, S ≤ 0, K ≤ 0
• Reject non-finite price
• Return nullopt, never NaN

Now the bounds check. A call cannot be worth less than max of S minus K e to the minus r T and zero, nor more than S itself. A put lies between max of K e to the minus r T minus S and zero, and K e to the minus r T. If the quoted price is outside, with a small tolerance, no volatility can match it, so return nullopt.

• Call in [max(S − Ke^(−rT), 0), S]
• Put in [max(Ke^(−rT) − S, 0), Ke^(−rT)]
• Outside: nullopt

This check prevents a whole class of bad searches. Without it, bisection would grind through two hundred iterations chasing a price no model can produce, then return garbage. With it, a stale or crossed quote is rejected immediately and clearly, which is also useful data about the quality of an input.

• Skip impossible searches
• Fail fast and clearly
• Bad quotes are data too

Next, set lo to one times ten to the minus four and hi to five. Evaluate the model at lo. If its price is already at or above the target, the answer is lo itself. Then evaluate at hi. If even that is below the target, return hi only when the shortfall is within tolerance, otherwise nullopt. Now the bracket is guaranteed.

• Range [1e-4, 5.0]
• f(lo) ≥ 0: answer is lo
• f(hi) ≤ 0: hi or nullopt

The loop itself is tiny. Take mid as half of lo plus hi, price at mid, and subtract the target to get f. If the absolute value of f is below tolerance, or the interval is narrower than tolerance, return mid. If f is negative, the price is too low, so lo becomes mid. Otherwise hi becomes mid. At most two hundred rounds.

• mid = (lo + hi) / 2
• f < 0: lo = mid
• f ≥ 0: hi = mid

Each pass halves the interval, so after about forty halvings five point zero shrinks below one times ten to the minus eleven. The cap of two hundred iterations is far beyond what is needed, and exists purely as a safety net. Two exit tests, price error or interval width, guard against both flat and steep regions.

• Width = 5 / 2ⁿ
• About 40 halvings → under 1e-11
• Cap: 200 iterations

Let us run it. Spot and strike are both one hundred, T is zero point two five years, r is zero point zero six five. At volatility zero point two, the model gives a call price of four point eight one three six. Pretend that is the quote, and forget the volatility. The bracket starts at nearly zero to five.

• S = K = 100, T = 0.25
• r = 0.065
• Target price 4.8136

First mid is two point five zero zero zero five. The model price there is too high by forty-two point four two, so hi moves down. Next mid is one point two five zero zero eight, error plus twenty point three four. Then zero point six two five zero nine, plus eight point three three. Then zero point three one two five nine, plus two point two zero.

• 2.50005 → f = +42.42
• 1.25008 → f = +20.34
• 0.62509 → f = +8.33
• 0.31259 → f = +2.20

Now zero point one five six three five. The error is minus zero point eight five. The model price is too low, so lo moves up. The root is trapped between zero point one five six and zero point three one three. Next, zero point two three four four seven gives plus zero point six seven, so hi moves down again.

• 0.15635 → f = −0.85, lo
• 0.23447 → f = +0.67, hi
• Root is bracketed

Next mid is zero point one nine five four one, error minus zero point zero nine, so lo moves. Then zero point two one four nine four, error plus zero point two nine, so hi moves. The guesses alternate around the target and squeeze in, converging to zero point two zero zero zero. We recovered the volatility we started with.

• 0.19541 → f = −0.09
• 0.21494 → f = +0.29
• Converges to 0.2000

Where does this code live in the platform? In four places. Derived option bars, from episode six, price options from the underlying. Alphas use volatility as an input. The portfolio greeks function aggregates risk. And the IV surface in episode thirteen solves for volatility across a whole chain.

• Derived bars (ep 6)
• Alpha volatility
• portfolio_greeks
• IV surface (ep 13)

The portfolio greeks function sums delta, gamma, vega and theta per underlying across all open option positions. Those totals can feed optional risk dot greeks limits, so a strategy can be capped on aggregate exposure. Because the pricer never returns NaN, the totals stay finite even when one contract is expired or degenerate.

• Per-underlying sums
• Delta, gamma, vega, theta
• Optional risk.greeks limits

Numerical code needs strong tests. Kestrel checks reference values against known prices, put-call parity, and the Greeks against finite differences of the price. Implied volatility is tested by round trips, price to volatility and back, and also by out-of-bounds prices, which must come back as nullopt.

• Reference values
• Put-call parity
• Greeks vs finite differences
• IV round trips and bounds

Three pitfalls. First, vega units: ours is per one point zero zero, so a one percent move is vega times zero point zero one. Second, time: using trading days instead of calendar time changes T, and the model's rate and decay no longer match the code. Kestrel uses three hundred sixty-five point two five days, consistently.

• Vega scale: 1.00 vs 1%
• Calendar time, not trading days
• Be consistent

Third, failure signalling. Returning NaN for a failed solve looks harmless, but every comparison with NaN is false, so it slips past checks and contaminates averages. A std optional forces the caller to decide what a missing volatility means. That is the design: make failure visible in the type.

• NaN fails every comparison
• NaN spreads silently
• optional makes failure explicit

To recap. Black-Scholes gives price and Greeks from d one and d two. Kestrel computes the normal cdf with erfc to keep the tails accurate, and guards expiry, zero volatility and bad inputs with deterministic fallbacks. Implied volatility checks no-arbitrage bounds, checks the endpoints, and bisects, because price is monotone in volatility.

• Formulas and Greeks
• erfc tails and edge cases
• Bounds, endpoints, halving

Newton asks the curve to be well behaved. Bisection only asks it to go up, and that is the one promise the model always keeps. Next time, episode thirteen, we build an entire implied volatility surface from an option chain.

• Next: IV surface (ep 13)
• Subscribe for more