Binomial lattice pricing

Deep dive

How it works

Parts:

• Step parameters — Δt = T/N, u, d, p, discount e^{-rΔt}; match the lognormal mean and variance.
• Rolling value array — v[j] = step-i node with j up moves, overwritten in place going backwards.
• Exercise rule — European: continuation; American: max(continuation, intrinsic).

1. Compute Δt, u, d, p; require 0 ≤ p ≤ 1.
2. Fill leaves v[j] = payoff(S·u^{2j-N}), j = 0..N.
3. For i = N-1 down to 0 and j = 0..i: v[j] = e^{-rΔt}(p·v[j+1] + (1-p)·v[j]); American: max with payoff(S·u^{2j-i}).
4. Save v[0..2] at i=1 and v[0..3] at i=2 (Greek taps).
5. Price = v[0]. Delta = (V₁ᵤ - V₁_d)/(S·u - S·d); gamma = difference of the two step-2 deltas over (S·u² - S·d²)/2; theta = (V₂_mid - V₀)/(2Δt).

S0 ─ u ─ uu ─ uuu        step i has i+1 nodes: j up moves, price S·u^(2j-i)
   \  \  / \  / \
    d   ud   uud          node (i,j) -> children (i+1,j+1) up, (i+1,j) down
      \  / \  /
       dd   udd           ud = du: the tree recombines, O(N²) nodes in all

Code

#[derive(Clone, Copy)]
pub struct Opt { pub s: f64, pub k: f64, pub r: f64, pub q: f64, pub sigma: f64, pub t: f64, pub call: bool, pub american: bool }
pub struct Out { pub price: f64, pub delta: f64, pub gamma: f64 }

/// CRR tree, n >= 3 steps, one rolling array.
pub fn crr(o: &Opt, n: usize) -> Out {
    assert!(n >= 3);
    let dt = o.t / n as f64;
    let u = (o.sigma * dt.sqrt()).exp();
    let d = 1.0 / u;
    let p = (((o.r - o.q) * dt).exp() - d) / (u - d);     // risk-neutral up probability
    assert!((0.0..=1.0).contains(&p), "p outside [0,1]: use more steps");
    let disc = (-o.r * dt).exp();
    let (dp, dq) = (disc * p, disc * (1.0 - p));            // precomputed once
    let sign = if o.call { 1.0 } else { -1.0 };
    let payoff = |s: f64| (sign * (s - o.k)).max(0.0);
    let u2 = u * u;
    let mut v = Vec::with_capacity(n + 1);
    let mut s = o.s * d.powi(n as i32);                     // lowest leaf: all down moves
    for _ in 0..=n { v.push(payoff(s)); s *= u2; }          // each extra up move multiplies by u^2
    let (mut v1, mut v2) = ([0.0; 2], [0.0; 3]);
    for i in (0..n).rev() {
        let mut s = o.s * d.powi(i as i32);
        for j in 0..=i {                                    // node (i, j): j up moves
            let cont = dp * v[j + 1] + dq * v[j];           // v[j+1] still holds step i+1
            v[j] = if o.american { cont.max(payoff(s)) } else { cont };
            s *= u2;
        }
        if i == 2 { v2.copy_from_slice(&v[..3]); }
        if i == 1 { v1.copy_from_slice(&v[..2]); }
    }
    let (suu, sdd) = (o.s * u2, o.s * d * d);
    let delta = (v1[1] - v1[0]) / (o.s * (u - d));
    let (d_up, d_dn) = ((v2[2] - v2[1]) / (suu - o.s), (v2[1] - v2[0]) / (o.s - sdd));
    Out { price: v[0], delta, gamma: (d_up - d_dn) / (0.5 * (suu - sdd)) }
}

/// Abramowitz & Stegun 7.1.26 (|erf error| <= 1.5e-7) and closed-form Black-Scholes-Merton.
pub fn norm_cdf(x: f64) -> f64 {
    let z = x.abs() / std::f64::consts::SQRT_2;
    let t = 1.0 / (1.0 + 0.3275911 * z);
    let poly = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429))));
    let half = 1.0 - 0.5 * poly * (-z * z).exp();
    if x >= 0.0 { half } else { 1.0 - half }
}
pub fn bsm(o: &Opt) -> f64 {
    let v = o.sigma * o.t.sqrt();
    let d1 = ((o.s / o.k).ln() + (o.r - o.q + 0.5 * o.sigma * o.sigma) * o.t) / v;
    let (fs, fk) = (o.s * (-o.q * o.t).exp(), o.k * (-o.r * o.t).exp());
    if o.call { fs * norm_cdf(d1) - fk * norm_cdf(d1 - v) } else { fk * norm_cdf(v - d1) - fs * norm_cdf(-d1) }
}

fn main() {
    let c = Opt { s: 100.0, k: 100.0, r: 0.05, q: 0.02, sigma: 0.2, t: 1.0, call: true, american: false };
    let p = Opt { call: false, ..c };
    let (ec, ep) = (crr(&c, 2000), crr(&p, 2000));
    assert!((ec.price - bsm(&c)).abs() < 5e-3 && (ep.price - bsm(&p)).abs() < 5e-3); // converges to BS
    let parity = ec.price - ep.price - 100.0 * ((-0.02f64).exp() - (-0.05f64).exp());
    assert!(parity.abs() < 1e-9);                                        // tree mean is exact: parity holds
    let nq = Opt { q: 0.0, ..c };
    let am = |o: Opt| crr(&Opt { american: true, ..o }, 2000).price;
    assert!((crr(&nq, 2000).price - am(nq)).abs() < 1e-12);              // no dividend: never exercise a call early
    assert!(am(p) > ep.price + 0.1);                                     // early-exercise premium on the put
    println!("call {:.4} (BS {:.4}) delta {:.4} gamma {:.4}", ec.price, bsm(&c), ec.delta, ec.gamma);
}

Compiled and run, asserts pass. S=K=100, r=5%, q=2%, σ=20%, T=1, N=2000: tree call 9.2260 vs closed form 9.2270; American put 6.660 vs European 6.329.

Complexity and constants

N(N+1)/2 internal nodes, each two multiplies, an add and an optional max. Memory 8(N+1) bytes: 80 KB at N=10,000, in cache, so the loop is compute-bound. Error is O(1/N): measured 1.9e-2, 1.9e-3, 3.8e-4 at N=100, 1,000, 5,000 on the call above; the sign alternates with the parity of N because the strike sits differently relative to the terminal nodes.

Efficiency

• Cost: Apple M4 Pro, one core, release build, put, N=100 / 1,000 / 10,000 (measured): European 0.001 / 0.063 / 5.2 ms (0.1-0.2 ns/node); American 0.004 / 0.49 / 51 ms (~1 ns/node; serial s *= u² chain plus max). Compute-bound.
• Levers:
  • Average N and N+1 — error 5.5e-4 at N=100/101 vs 1.9e-3 at N=1,000, about 50x fewer nodes (measured) — any option whose payoff is kinked at one strike.
  • Rolling array + precomputed disc·p, disc·(1-p) — about 4x vs the naive row below (measured) — always.
  • Early-exercise boundary — put nodes far above the boundary are never exercised, so skip payoff and max there (gain depends on parameters, not measured) — American puts.
  • Smoothing — closed form at the last step (Broadie-Detemple) or a Leisen-Reimer tree: no oscillation — accuracy at small N.
  • Batching — independent options in the same loop; SIMD lanes across strikes — option chains, portfolios.
• Naive → standard → state of the art:
  • Full 2-D tree, powi per node, fresh Vec per level · ~3.7-4.5 ns/node American (measured)
  • Rolling array, precomputed factors, running price product · ~1 ns/node
  • Smoothed or Leisen-Reimer tree, N of order 100, batched with SIMD · accuracy of N of order 1,000 at a fraction of the nodes

Implementation notes

• In place, ascending j: v[j+1] is still the old layer when v[j] is overwritten.
• Node prices are a running product S·d^i·(u²)^j; no table, no per-node powf.
• Delta from the step-1 pair, gamma from the step-2 triple; QuantLib's BinomialVanillaEngine does the same and derives theta from delta and gamma via the Black-Scholes PDE. Vega, rho: bump and reprice.
• p ∈ [0,1] needs σ√Δt ≥ |r-q|Δt.
• Proportional dividends: multiply node prices by (1-δ) after the ex-date (still recombines). Cash dividends: lattice the price minus the present value of dividends still to come.
• Jarrow-Rudd uses equal probabilities with drift-adjusted u, d; trinomial trees add a middle node (Hull-White rate trees use them) and suit mean-reverting processes.

Pitfalls

• Trusting one N: oscillation hides error; compare N and N+1.
• With q > 0 American and European calls differ.
• Low σ with few steps: p leaves [0,1].
• Rate lattices must be fitted to today's curve (yield curve bootstrapping).

Variants and alternatives

In production

• Equity/FX: QuantLib's BinomialVanillaEngine with CoxRossRubinstein, JarrowRudd, Tian, LeisenReimer, Joshi4 trees.
• Corporate finance: real options on energy, mining and pharma projects; employee stock options.
• Fixed income, insurance: convertible and callable bonds, surrender-right guarantees.

Implementations

Sources

• Cox, Ross, Rubinstein, Option pricing: a simplified approach, Journal of Financial Economics 7(3), 1979, 229-263. paper
• Leisen, Reimer, Binomial models for option valuation - examining and improving convergence, Applied Mathematical Finance 3(4), 1996, 319-346. paper
• Broadie, Detemple, American option valuation: new bounds, approximations, and a comparison of existing methods, Review of Financial Studies 9(4), 1996, 1211-1250. paper
• Black, Derman, Toy, A one-factor model of interest rates and its application to Treasury bond options, Financial Analysts Journal 46(1), 1990, 33-39. paper
• Hull, White, Numerical procedures for implementing term structure models I: single-factor models, Journal of Derivatives 2(1), 1994, 7-16. paper
• QuantLib binomialengine.hpp, binomialtree.hpp