Newton-Raphson and Brent for Implied Volatility

Setup

Market Context

Implied volatility $\sigma^*(K, T)$ is the unique value of $\sigma$ such that the Black-Scholes price matches the observed market price of an option with strike $K$ and expiry $T$. It is the standard language in which option markets quote prices: traders say "I buy the 1-month 25-delta put at 18 vol", not "at price €3.47".

Computing implied vol is therefore an inversion problem, performed thousands of times per second in a live trading system. Correctness and speed both matter. A slow inverter causes pricing latency; a numerically fragile one produces garbage implied vols for deep OTM strikes, polluting the vol surface.

Conventions

Throughout this module:

• $C_{\mathrm{BS}}(\sigma; S, K, r, q, T)$: Black-Scholes call price, where $S$ is spot, $K$ is strike, $r$ is the continuously compounded risk-free rate, $q$ is the continuous dividend yield, $T$ is time to expiry in years (Act/365 unless stated).
• Volatility $\sigma$ is annualised, in decimal ($0.20 = 20\%$).
• $F = S e^{(r-q)T}$: the forward price.
• $d_1 = \frac{\ln(F/K) + \frac{1}{2}\sigma^2 T}{\sigma\sqrt{T}}$, $\quad d_2 = d_1 - \sigma\sqrt{T}$.

The Black-Scholes call price:

$C_{\mathrm{BS}}(\sigma) = e^{-rT}\bigl[F \Phi(d_1) - K \Phi(d_2)\bigr].$

Assumptions before use:

1. Constant volatility $\sigma$ over $[0, T]$.
2. Lognormal spot dynamics under the risk-neutral measure.
3. No early exercise (European option).
4. Liquid, arbitrage-free market quote $C_{\mathrm{mkt}}$ with $C_{\mathrm{mkt}} \in \bigl((F - K)^+\, e^{-rT},\; S\, e^{-qT}\bigr)$ — the no-arbitrage bounds.

Theory: Implied Vol as a Root-Finding Problem

Define the objective function:

$f(\sigma) = C_{\mathrm{BS}}(\sigma) - C_{\mathrm{mkt}}.$

We seek $\sigma^* > 0$ such that $f(\sigma^*) = 0$. This function has the following properties that govern algorithm design:

1. Monotonicity. For $\sigma > 0$, $C_{\mathrm{BS}}$ is strictly increasing in $\sigma$: $f'(\sigma) = \mathcal{V} > 0$, where $\mathcal{V} = S e^{-qT} \sqrt{T}\, \phi(d_1)$ is the Black-Scholes vega and $\phi = \Phi'$ is the standard normal density. Monotonicity guarantees at most one root.
2. Boundary behaviour. $C_{\mathrm{BS}}(0) = e^{-rT}(F - K)^+$ (intrinsic value) and $C_{\mathrm{BS}}(\infty) = S e^{-qT}$ (full asset price). For $C_{\mathrm{mkt}}$ in the arbitrage-free interior, a root exists.
3. Vega degeneracy. At extreme strikes ($|d_1| \to \infty$), $\phi(d_1) \to 0$, so $\mathcal{V} \to 0$. Newton-Raphson $\sigma_{n+1} = \sigma_n - f(\sigma_n)/f'(\sigma_n)$ divides by a near-zero vega and can produce catastrophic oscillations or divergence.

Newton-Raphson

Algorithm

Starting from an initial guess $\sigma_0$, Newton-Raphson iterates:

$\sigma_{n+1} = \sigma_n - \frac{f(\sigma_n)}{f'(\sigma_n)} = \sigma_n - \frac{C_{\mathrm{BS}}(\sigma_n) - C_{\mathrm{mkt}}}{\mathcal{V}(\sigma_n)}.$

Convergence rate. Near a simple root, Newton-Raphson converges quadratically: if $|\sigma_n - \sigma^*| = \varepsilon$, then $|\sigma_{n+1} - \sigma^*| = O(\varepsilon^2)$. The number of correct decimal digits doubles per iteration. For a good initial guess, 3–5 iterations suffice.

Formal proof. By Taylor expansion around $\sigma^*$:

$f(\sigma_n) = f(\sigma^*) + f'(\sigma^*)(\sigma_n - \sigma^*) + \tfrac{1}{2}f''(\xi)(\sigma_n - \sigma^*)^2 = f'(\sigma^*)(\sigma_n - \sigma^*) + O((\sigma_n - \sigma^*)^2),$

since $f(\sigma^*) = 0$. Thus:

$\sigma_{n+1} - \sigma^* = \sigma_n - \sigma^* - \frac{f(\sigma_n)}{f'(\sigma_n)} = -\frac{f''(\xi)}{2f'(\sigma_n)}(\sigma_n - \sigma^*)^2 + O((\sigma_n - \sigma^*)^3).$

The error is bounded by $|e_{n+1}| \leq C |e_n|^2$ where $C = |f''(\sigma^*)| / (2|f'(\sigma^*)|)$ is the asymptotic constant.

Initial Guess

A poor initial guess causes slow convergence or divergence. The Brenner-Subrahmanyam (1988) approximation for near-ATM options provides a first-order starting point:

$\sigma_0 \approx \sqrt{\frac{2\pi}{T}} \cdot \frac{C_{\mathrm{mkt}}}{S e^{-qT}}.$

This is a linearisation of the ATM Black-Scholes formula $C_{\mathrm{ATM}} \approx S e^{-qT} \sigma \sqrt{T/(2\pi)}$. For skewed strikes, a better guess uses the Corrado-Miller approximation.

Failure Modes

Deep OTM puts and calls. When $|d_1| > 4$, the vega is below $10^{-8}$ per unit notional. Newton-Raphson steps become enormous and overshoot wildly. The algorithm must detect this case and switch strategies.

Near-zero time to expiry. As $T \to 0$, the whole vol surface compresses. Vega scales as $\sqrt{T} \to 0$. Newton-Raphson becomes unreliable for $T < 0.001$ years (less than one trading day).

Near-intrinsic options. If $C_{\mathrm{mkt}} \approx (F - K)^+ e^{-rT}$ (i.e., the option is priced at near-intrinsic), the root is near $\sigma = 0$ and both $f(\sigma)$ and $\mathcal{V}(\sigma)$ are near zero simultaneously. l'Hôpital's rule applies in the limit but the numerical ratio is unstable.

Brent's Method

Brent's method (1973) combines three strategies: bisection, secant, and inverse quadratic interpolation (IQI). It is guaranteed to converge to machine precision within $O(\log_2(1/\varepsilon))$ function evaluations, making it the default fallback.

Bracketing

Brent requires a bracket $[\sigma_{\mathrm{lo}}, \sigma_{\mathrm{hi}}]$ with $f(\sigma_{\mathrm{lo}}) \cdot f(\sigma_{\mathrm{hi}}) < 0$. Since $f$ is monotone and bounded:

$\sigma_{\mathrm{lo}} = 10^{-8} \quad (\text{near zero but finite}), \qquad \sigma_{\mathrm{hi}} = 10.0 \quad (1000\%\ \text{vol}).$

This bracket is valid for any arbitrage-free market quote. Initial $f(\sigma_{\mathrm{lo}}) < 0$ (intrinsic value is below the market) and $f(\sigma_{\mathrm{hi}}) > 0$ (deep-ITM call at 1000 vol recovers full asset value), so the sign condition is satisfied.

Algorithm Sketch

Maintain three points $a$, $b$, $c$ with $f(a) \cdot f(b) \leq 0$, $|f(b)| \leq |f(a)|$. At each step:

1. IQI if the last two steps used secant and the three points are distinct: fit a quadratic through $(a, f(a))$, $(b, f(b))$, $(c, f(c))$ and take the root.
2. Secant if the last step was not secant: use the linear interpolant through $(a, f(a))$, $(b, f(b))$.
3. Bisection if the candidate from steps 1/2 lies outside $(b, (a+b)/2)$ or convergence is slow: take the midpoint.

Always bisect when in doubt. The convergence guarantee (superlinear, between order 1 and order 2 in the best case) comes from this fallback. In practice, Brent converges in 5–10 iterations for vol inversion.

Convergence Guarantee

Theorem (Brent 1973). Given a continuous function $f$ and bracket $[a, b]$ with $f(a) f(b) < 0$, Brent's method converges to a root in at most $O(\log_2((b-a)/\varepsilon))$ function evaluations, where $\varepsilon$ is the tolerance.

For $b - a = 10$ and $\varepsilon = 10^{-10}$, this is at most $\approx 33$ evaluations. In practice, far fewer are needed.

Hybrid Newton/Brent

Production implied vol solvers use a hybrid:

1. Initialise with a bracket $[\sigma_{\mathrm{lo}}, \sigma_{\mathrm{hi}}]$ and a starting guess $\sigma_0$ (Brenner-Subrahmanyam or similar).
2. Attempt Newton-Raphson if the candidate $\sigma_{n+1}$ lies within the current bracket and vega $\mathcal{V}(\sigma_n) > \mathcal{V}_{\min}$.
3. Fall back to bisection (or Brent) if Newton's candidate exits the bracket or vega is degenerate.
4. Update bracket: narrow $[\sigma_{\mathrm{lo}}, \sigma_{\mathrm{hi}}]$ each iteration using the sign of $f$.

This gives quadratic convergence near ATM (where Newton is well-conditioned) and guaranteed convergence everywhere else. A well-implemented hybrid solves the inverse problem to $10^{-10}$ tolerance in 3–7 function evaluations for the entire vol surface.

Implementation

#include <cmath>
#include <stdexcept>
#include <limits>
#include <algorithm>

// Black-Scholes cumulative normal and density
static double bs_phi(double x) noexcept {
    return std::exp(-0.5 * x * x) / std::sqrt(2.0 * M_PI);
}
static double bs_Phi(double x) noexcept {
    return 0.5 * std::erfc(-x / std::sqrt(2.0));
}

// Black-Scholes call price and vega
// F: forward price = S * exp((r-q)*T)
// K: strike, T: time to expiry (years), df: discount factor exp(-r*T)
struct BsResult { double price; double vega; };

static BsResult bs_call(double F, double K, double T, double df) noexcept {
    if (T <= 0.0) {
        double intrinsic = df * std::max(F - K, 0.0);
        return {intrinsic, 0.0};
    }
    // avoid division by zero in degenerate cases
    if (K <= 0.0) return {df * F, 0.0};
    return {df * F, 0.0};  // placeholder; full impl below
}

// Full implementation
struct ImpliedVolResult {
    double vol;             // implied vol (annualised)
    int    iterations;      // function evaluations used
    bool   converged;
};

// Hybrid Newton-Raphson / Brent implied vol inverter.
//
// Assumptions:
//   - European call option (put via put-call parity)
//   - Continuously compounded rates
//   - market_price must satisfy no-arbitrage bounds
//
// Parameters:
//   market_price: observed call price
//   S           : spot price
//   K           : strike
//   r           : risk-free rate (continuously compounded, annualised)
//   q           : continuous dividend yield (annualised)
//   T           : time to expiry (years)
//   tol         : target absolute tolerance on vol (default 1e-10)
//   max_iter    : maximum iterations
ImpliedVolResult implied_vol(
    double market_price,
    double S, double K,
    double r, double q, double T,
    double tol     = 1.0e-10,
    int    max_iter = 100
) {
    if (T <= 0.0)
        throw std::domain_error("T must be positive");
    if (S <= 0.0 || K <= 0.0)
        throw std::domain_error("S and K must be positive");

    const double F  = S * std::exp((r - q) * T);
    const double df = std::exp(-r * T);

    // No-arbitrage bounds check
    const double lo_bound = df * std::max(F - K, 0.0);
    const double hi_bound = S * std::exp(-q * T);
    if (market_price <= lo_bound || market_price >= hi_bound)
        throw std::domain_error("market_price violates no-arbitrage bounds");

    // BS price and vega as a function of sigma
    auto price_and_vega = [&](double sigma) -> BsResult {
        const double sqrtT = std::sqrt(T);
        const double d1 = (std::log(F / K) + 0.5 * sigma * sigma * T) / (sigma * sqrtT);
        const double d2 = d1 - sigma * sqrtT;
        const double price = df * (F * bs_Phi(d1) - K * bs_Phi(d2));
        const double vega  = S * std::exp(-q * T) * sqrtT * bs_phi(d1);
        return {price, vega};
    };

    // Bracket: [sigma_lo, sigma_hi]
    double s_lo = 1.0e-8;
    double s_hi = 10.0;   // 1000% vol — safe upper bound

    double f_lo = price_and_vega(s_lo).price - market_price;
    double f_hi = price_and_vega(s_hi).price - market_price;

    // Should always be true for valid inputs
    if (f_lo * f_hi >= 0.0)
        throw std::domain_error("Failed to bracket root");

    // Brenner-Subrahmanyam initial guess (near-ATM approximation)
    double sigma = std::sqrt(2.0 * M_PI / T) * (market_price / (S * std::exp(-q * T)));
    sigma = std::clamp(sigma, s_lo, s_hi);

    int iters = 0;
    for (; iters < max_iter; ++iters) {
        auto [price, vega] = price_and_vega(sigma);
        double f = price - market_price;

        // Update bracket using sign of f
        if (f * f_lo < 0.0) {
            s_hi = sigma; f_hi = f;
        } else {
            s_lo = sigma; f_lo = f;
        }

        // Convergence check
        if (std::abs(f) < tol && std::abs(s_hi - s_lo) < tol)
            return {sigma, iters + 1, true};

        // Newton step — only if vega is non-degenerate
        constexpr double vega_floor = 1.0e-12;
        if (vega > vega_floor) {
            double s_newton = sigma - f / vega;
            if (s_newton > s_lo && s_newton < s_hi) {
                sigma = s_newton;
                continue;  // accept Newton step
            }
        }

        // Fallback: bisection
        sigma = 0.5 * (s_lo + s_hi);
    }

    return {sigma, iters, false};  // did not converge
}

Put Options via Put-Call Parity

The inverter above handles calls. For a put:

$C_{\mathrm{BS}} = P_{\mathrm{BS}} + e^{-rT}(F - K),$

so convert $P_{\mathrm{mkt}}$ to a call price first:

// Convert put to synthetic call via put-call parity
double put_to_call(double put_price, double S, double K, double r, double q, double T) {
    double F  = S * std::exp((r - q) * T);
    double df = std::exp(-r * T);
    return put_price + df * (F - K);  // C = P + df*(F - K)
}

Validation

Round-Trip Test

The fundamental correctness test: generate a Black-Scholes price from a known vol, then invert it and verify recovery.

// Catch2 test — compile with -DCATCH_CONFIG_MAIN
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>

TEST_CASE("Implied vol: round-trip ATM", "[implied_vol]") {
    const double S = 100.0, K = 100.0, r = 0.05, q = 0.02, T = 1.0;
    const double true_vol = 0.20;
    // price at true_vol
    const double sqrtT = std::sqrt(T), F = S * std::exp((r-q)*T);
    const double d1 = (std::log(F/K) + 0.5*true_vol*true_vol*T) / (true_vol*sqrtT);
    const double d2 = d1 - true_vol * sqrtT;
    const double df = std::exp(-r*T);
    const double cprice = df * (F * bs_Phi(d1) - K * bs_Phi(d2));

    auto result = implied_vol(cprice, S, K, r, q, T);
    REQUIRE(result.converged);
    REQUIRE_THAT(result.vol, Catch::Matchers::WithinAbs(true_vol, 1.0e-8));
}

TEST_CASE("Implied vol: deep OTM call", "[implied_vol]") {
    const double S = 100.0, K = 150.0, r = 0.05, q = 0.0, T = 0.5;
    const double true_vol = 0.30;
    const double sqrtT = std::sqrt(T), F = S * std::exp(r*T);
    const double d1 = (std::log(F/K) + 0.5*true_vol*true_vol*T) / (true_vol*sqrtT);
    const double d2 = d1 - true_vol * sqrtT;
    const double df = std::exp(-r*T);
    const double cprice = df * (F * bs_Phi(d1) - K * bs_Phi(d2));

    auto result = implied_vol(cprice, S, K, r, q, T);
    REQUIRE(result.converged);
    REQUIRE_THAT(result.vol, Catch::Matchers::WithinAbs(true_vol, 1.0e-7));
}

Convergence Across the Surface

A robust inverter must work for the full range of strikes and maturities encountered on a real vol surface. Typical production test: a 10 × 10 grid of strikes $K/F \in [0.5, 2.0]$ and maturities $T \in [0.01, 5.0]$, at four vol levels (5%, 20%, 50%, 100%). The hybrid Newton/Brent solver described above converges to $10^{-10}$ tolerance in all cases within 10 iterations.

Limitations

Model dependency. Implied vol is defined relative to a pricing model. Black-Scholes implied vol is the conventional market standard, but it is a model-dependent quantity. Under a different model (e.g., normal Black model, shifted lognormal), implied vol means something different. Always state the model convention.

Non-uniqueness for exotic payoffs. For exotic options without closed-form pricing, implied vol is typically not well-defined as a scalar. A vol surface (rather than a single number) is used for calibration.

Negative vega for calendar spreads. Calendar spreads can have non-monotone value in vol, breaking the single-root guarantee. Implied vol is only well-defined for vanilla options (call, put, digital).

Arbitrage-free bounds violation. Market data occasionally contains crossed quotes (bid above no-arbitrage bound, or ask below intrinsic). The inverter must detect this and either error out or apply a bid-offer spread adjustment before inversion.

Numerical precision near intrinsic. When $C_{\mathrm{mkt}}$ is within a few cents of $e^{-rT}(F-K)^+$, the implied vol is near zero and the inversion is numerically ill-conditioned. Tolerances must be relaxed or an alternative formula (e.g., Bachelier normal vol) used.

Interview Angle

L1. Why is implied vol not the same as realised vol? What does it represent? Explain why implied vol is monotone in Black-Scholes — i.e., why does the call price always increase with $\sigma$?

Implied vol is the market-implied expectation of future vol embedded in option prices, under the Black-Scholes model assumption. Realised vol is an ex-post measurement of actual price moves. They differ because of the variance risk premium (option sellers demand compensation for variance risk), demand/supply imbalances, and jump risk not captured in GBM.

Monotonicity: $\partial C_{\mathrm{BS}}/\partial\sigma = \mathcal{V} = S e^{-qT} \sqrt{T} \phi(d_1) > 0$ for all $\sigma > 0$, $T > 0$, because $\phi > 0$. Higher vol increases the probability of large moves in both directions, always increasing option value.

L2. Why does Newton-Raphson fail for deep OTM options? How does the hybrid solver avoid this? Describe the Brent bracket and the bracketing invariant maintained throughout iteration.

NR failure: Vega $\mathcal{V} = S e^{-qT}\sqrt{T}\phi(d_1)$ is exponentially small for large $|d_1|$. The Newton step $-f/\mathcal{V}$ divides a finite residual by near-zero, producing a step that overshoots the bracket by orders of magnitude. The hybrid detects when the Newton candidate exits the bracket and replaces it with a bisection step, maintaining convergence.

Bracket invariant: At every iteration, $f(\sigma_{\mathrm{lo}}) < 0$ and $f(\sigma_{\mathrm{hi}}) > 0$ (or vice versa). This is preserved by updating the bracket endpoint whose $f$-value has the same sign as the current iterate.

L3. Analyse the convergence rate of the hybrid solver. Under what conditions does it achieve quadratic convergence, and when does it degrade to linear? Discuss the Corrado-Miller approximation as an alternative to Brenner-Subrahmanyam for off-ATM strikes, and its accuracy relative to the Newton iterate.

Convergence rate: Quadratic convergence is achieved when the Newton step remains within the bracket for every iteration — i.e., when the problem is well-conditioned (vega bounded away from zero). Brent's IQI achieves superlinear convergence (order ~1.84 asymptotically) in the general case. Pure bisection is linear (order 1), halving the bracket width each step.

The Corrado-Miller (1996) approximation improves on Brenner-Subrahmanyam by including a first-order correction for the moneyness $m = F e^{-rT} - K e^{-rT}$:

$\sigma_0 \approx \sqrt{\frac{2\pi}{T}} \cdot \frac{C_{\mathrm{mkt}} - m/2 + \sqrt{(C_{\mathrm{mkt}} - m/2)^2 - m^2/\pi}}{S e^{-qT} + K e^{-rT}},$

providing a starting point accurate to $O(\sigma^3)$ in moneyness. For typical equity options, this initial guess is within 1–2 vol points of $\sigma^*$, reducing Newton iterations from 5–7 to 2–3.