Implied Volatility in Python: Why Newton Fails

You pull an option chain, loop over the quotes, call a Newton-Raphson solver for implied volatility, and most rows come back fine. Then a handful return nan, a negative number, or a volatility of 1,400%. You add a try/except, drop the bad rows, and fit your smile to what is left.

That is how most implied volatility code is written, and it is exactly the code that breaks on the strikes a desk cares about most: the short-dated wings. This guide shows why the textbook Newton iteration fails, how to write an implied volatility solver in Python that converges for every arbitrage-free price, and which market-data checks matter more than the root-finder itself.

1. The problem: inverting Black-Scholes

Implied volatility is the single number \(\sigma\) that makes the Black-Scholes (Black-76, when written on the forward) price equal the observed market price. With forward \(F\), strike \(K\), expiry \(T\) and discount factor \(D\):

There is no closed-form inverse, so we solve \(C(\sigma) = C_{\text{mkt}}\) numerically. The good news is that the problem is well posed: the price is strictly increasing in volatility (vega is positive), so a solution exists and is unique whenever the price lies inside the no-arbitrage bounds

If a solver fails on a price inside those bounds, the solver is the bug, not the market.

2. Why naive Newton-Raphson fails

The textbook iteration is \(\sigma_{n+1} = \sigma_n - \big(C(\sigma_n) - C_{\text{mkt}}\big)/\text{vega}(\sigma_n)\), started from something like 20%:

def naive_newton(price, F, K, T, df=1.0, call=True, sigma0=0.2, tol=1e-10):
          sigma = sigma0
          for _ in range(100):
              diff = bs_price(F, K, T, sigma, df, call) - price
              if abs(diff) < tol:
                  return sigma
              sigma -= diff / bs_vega(F, K, T, sigma, df)   # nothing stops this going negative
          raise RuntimeError("no convergence")

Try it on an ordinary wing option: forward 100, strike 130, 0.1 years to expiry, discount factor 0.99, true volatility 35%. The price is about 0.0371. Starting from 20%, this is what happens:

Iteration	\(\sigma_n\)	Vega	Price error	Next \(\sigma\)
0	0.2000	0.0026	−0.0371	14.41
1	14.41	1.06	+96.4	−76.35

At 20% the option is so far out of the money that vega is almost zero, so the Newton step (price error divided by vega) is enormous: the solver jumps to 1,441% volatility. There the price is close to its upper bound, so the next step overshoots to a meaningless negative volatility (−76), and within a few more steps vega hits zero and the run ends in nan. Nothing is wrong with the option. The root-finder was simply started on the wrong side of a curve it cannot handle.

The underlying reason is shape. As a function of \(\sigma\), the option price is convex, then concave. Differentiating vega gives volga (vomma):

So the curvature changes sign exactly where \(\sigma^2 T = 2|x|\), at the inflection point

Newton's method is only guaranteed to behave on a region where the function does not change curvature. Start too low on the convex part and a tiny vega throws you far to the right; start too high on the concave part and you overshoot to the left. For the example above, \(\sigma^{*} = \sqrt{2\ln(1.3)/0.1} \approx 2.29\): the true root 0.35 and the 0.20 starting guess are both on the convex side, far below it, which is why the first step explodes.

3. The fix: start at the inflection point and keep a bracket

Two ideas make the solver robust.

1. Start at \(\sigma^{*}\). Manaster and Koehler (1982) showed that Newton's iteration started from this point converges monotonically to the implied volatility. It is also a far better guess than a constant 20% for short-dated wings.
2. Keep a bracket and fall back to bisection. Because the price is monotone in \(\sigma\), every evaluation tells you which side of the root you are on. Keep an interval \([\sigma_{lo}, \sigma_{hi}]\) that always contains the root; accept a Newton step only if it lands inside it, otherwise bisect. This is the classic safeguarded Newton, and it cannot diverge.

Two more details matter in practice:

• Solve on the out-of-the-money option. An in-the-money call is intrinsic value plus a small amount of time value, and all the volatility information lives in that small part. Converting to the OTM put with put-call parity, \(C - P = D(F-K)\), makes the problem better conditioned.
• Use a tolerance relative to the OTM price. An absolute tolerance of \(10^{-10}\) is meaningless for a wing option worth \(10^{-9}\).

4. A robust implied volatility solver in Python

import math
      from scipy.stats import norm

      def bs_price(F, K, T, sigma, df=1.0, call=True):
          """Black (forward) price. df = discount factor to expiry."""
          if sigma <= 0 or T <= 0:
              intrinsic = max(F - K, 0.0) if call else max(K - F, 0.0)
              return df * intrinsic
          v = sigma * math.sqrt(T)
          d1 = math.log(F / K) / v + 0.5 * v
          d2 = d1 - v
          if call:
              return df * (F * norm.cdf(d1) - K * norm.cdf(d2))
          return df * (K * norm.cdf(-d2) - F * norm.cdf(-d1))

      def bs_vega(F, K, T, sigma, df=1.0):
          v = sigma * math.sqrt(T)
          d1 = math.log(F / K) / v + 0.5 * v
          return df * F * norm.pdf(d1) * math.sqrt(T)

      def implied_vol(price, F, K, T, df=1.0, call=True,
                      tol=1e-12, max_iter=100, lo=1e-6, hi=5.0):
          """Safeguarded Newton: Newton steps inside a bisection bracket."""
          # 1. No-arbitrage bounds (undiscounted)
          p = price / df
          intrinsic = max(F - K, 0.0) if call else max(K - F, 0.0)
          upper = F if call else K
          if not (intrinsic < p < upper):
              raise ValueError(f"price {price:.6g} outside no-arbitrage bounds "
                               f"({df*intrinsic:.6g}, {df*upper:.6g})")
          # 2. Use the OTM side via put-call parity (better conditioned)
          if call and K < F:
              p, call = p - (F - K), False
          elif not call and K > F:
              p, call = p + (F - K), True
          f = lambda s: bs_price(F, K, T, s, 1.0, call) - p
          # 3. Start at the inflection point of price(sigma)
          x = math.log(F / K)
          sigma = math.sqrt(2.0 * abs(x) / T) if x != 0 else 0.2
          sigma = min(max(sigma, lo), hi)
          while f(hi) < 0:                  # widen bracket for extreme prices
              hi *= 2.0
          for i in range(max_iter):
              diff = f(sigma)
              if abs(diff) <= tol * p:          # relative to the OTM price
                  return sigma
              if diff > 0: hi = sigma
              else:        lo = sigma
              vega = bs_vega(F, K, T, sigma)
              step_ok = vega > 1e-14
              if step_ok:
                  cand = sigma - diff / vega
                  step_ok = lo < cand < hi
              sigma = cand if step_ok else 0.5 * (lo + hi)
              if hi - lo < 1e-15:
                  return sigma
          return sigma

A round-trip test prices options across strikes from 40 to 250 (forward 100), expiries from one day to five years and volatilities from 5% to 150%, then solves for implied volatility and compares. For every out-of-the-money quote, the recovered volatility matches the true one to within about \(10^{-12}\), and there were no failures. When the same deep in-the-money options are passed as ITM prices, the error grows to around \(10^{-7}\). That is not the solver: a 150.00 option whose time value is \(1.5\times10^{-10}\) simply does not carry enough significant digits in double precision. Quote the OTM side.

If you only need something quick, scipy.optimize.brentq on \(C(\sigma) - C_{\text{mkt}}\) over a bracket such as \([10^{-6}, 5]\) is also robust, just slower than a safeguarded Newton. For production speed and machine-precision accuracy across all regimes, Peter Jäckel's Let's Be Rational (2015) algorithm reaches full double precision in at most two iterations; the open-source py_vollib library wraps an implementation of it.

5. A good closed-form first guess

For near-the-money options, Brenner and Subrahmanyam (1988) give an approximation from the fact that \(N(d_1) - N(d_2) \approx \sigma\sqrt{T}/\sqrt{2\pi}\) at the money:

For an at-the-money option with 25% volatility and six months to expiry it gives 24.97%. It is excellent as a starting point near the money and poor in the wings, which is why the inflection point is the safer default.

6. The market-data problems that matter more than the solver

Once the root-finder is reliable, almost every remaining "bad implied vol" comes from the inputs.

Symptom	Usual cause	Fix
Price below intrinsic, solver refuses	Stale or crossed quote, wrong forward, or an American option priced as European	Check bounds first and report the row; do not clip it silently
Calls and puts at the same strike give different vols	Wrong forward or dividend/borrow assumption	Imply the forward from put-call parity near the money
Wings jump around between snapshots	Using mid when the bid is zero or the spread is wide	Compute bid and ask vols; flag rows where the band is wider than your tolerance
Short-dated vols look too high	Calendar-day vs trading-day time, or wrong expiry time of day	Measure \(T\) consistently, down to the expiry cut-off
ITM single-stock options give odd vols	Early-exercise premium in American options	Use OTM options, or de-Americanise with a tree before inverting

A useful rule is to report an implied volatility band, not a point: the bid vol and the ask vol. If a smile fit sits inside every band, it is consistent with the market, even if it misses every mid.

import pandas as pd

      def chain_iv(chain: pd.DataFrame, F: float, T: float, df: float) -> pd.DataFrame:
          """chain columns: strike, bid, ask, is_call. Uses OTM quotes only."""
          otm = chain[(chain.is_call & (chain.strike >= F)) |
                      (~chain.is_call & (chain.strike <= F))].copy()
          for side in ("bid", "mid", "ask"):
              px = otm[side] if side != "mid" else 0.5 * (otm.bid + otm.ask)
              otm[f"iv_{side}"] = [
                  implied_vol(p, F, k, T, df, c) if p > 0 else float("nan")
                  for p, k, c in zip(px, otm.strike, otm.is_call)
              ]
          return otm

7. The interview version

"How would you compute implied volatility?" is a common quant interview question, and "Newton-Raphson" is only the start of a good answer. A strong answer covers:

1. The solution exists and is unique inside the no-arbitrage bounds, because vega is positive.
2. Plain Newton can diverge because price is convex then concave in \(\sigma\) and vega vanishes in the wings.
3. Start at \(\sigma^{*} = \sqrt{2|\ln(F/K)|/T}\) and safeguard the step with a bisection bracket.
4. Invert OTM options, use relative tolerances, and treat bid/ask as a band.
5. Know that Jäckel's method is the production standard when speed and precision both matter.

8. Checklist

• Compute the forward and discount factor first; solve in forward (Black) terms.
• Reject prices outside \(\big(D\max(F-K,0),\,DF\big)\) and log them.
• Convert to the OTM option with put-call parity.
• Start from the inflection point; keep a bracket; bisect when Newton leaves it.
• Use a tolerance relative to the OTM price.
• Round-trip test across strikes, expiries and volatilities before trusting the code.
• Produce bid and ask vols, not just mid.

Sources and further reading

• S. Manaster and G. Koehler (1982), "The Calculation of Implied Variances from the Black-Scholes Model: A Note", Journal of Finance 37(1) — the inflection-point starting value and monotone convergence of Newton's method.
• M. Brenner and M. Subrahmanyam (1988), "A Simple Formula to Compute the Implied Standard Deviation", Financial Analysts Journal 44(5) — the at-the-money approximation.
• P. Jäckel (2015), "Let's Be Rational", Wilmott — machine-precision implied volatility in two iterations.
• SciPy documentation: scipy.optimize.brentq — a bracketed root-finder you can use as a simple robust fallback.
• py_vollib on GitHub — Python option pricing and implied volatility built on Jäckel's algorithm.

Continue on Desk2Quant

Implied volatility is the input to everything downstream. The Vol Surface Construction Playbook takes these implied vols and builds an arbitrage-free SVI/SSVI surface, Numerical Methods for Quants covers root-finding, stability and calibration in depth, and Python for Quants is the interview-ready companion for the code. For what vega's own derivatives do to a book, read The Greeks You Were Never Taught.

Educational content only. Code is provided as-is for learning; validate it against your own data and systems before any production use.