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.