SVI → Transformed PDE → Crank–Nicolson

1) Fit the implied-vol surface with SVI

Goal. Replace noisy per-quote implied vols with a smooth, arbitrage-aware surface spanning the pricing domain.

Raw SVI (per maturity)

w(k, T) = a + b · (ρ(k - m) + √{(k - m)² + σ²}),    IV(k,T)=√{w/T}.

• x-axis: log-moneyness k=ℓn(K/S), built from delta inversion.
• Symmetric delta grid (e.g., 0.10…0.90) for consistent wings/ATM coverage.

Delta → log-moneyness

• Δ_c=Φ(d₁), Δ_p=Φ(d₁)-1 ⇒ d₁=Φ⁻¹(Δ_c) or Φ⁻¹(Δ_p+1)
• k=ℓn(K/S)=(r-q-½σ²)T+σ√T·d₁

2) Feed the SVI surface into the pricing PDE

Mapping. Use Δ ↔ d_1 ↔ k=ln K/S to connect Δ to x. Interpolate SVI to get σ(x,τ) on PDE nodes.

3) Transformed PDE in log-price & backward time

Original Black-Scholes PDE

Standard form in forward time t and asset price S:

Subject to terminal payoff V(S,T) at expiration and appropriate boundary conditions.

Variable transformation

Apply coordinate transforms to simplify the PDE:

Transformed PDE

Define time-and-space dependent coefficients:

The transformed PDE becomes:

Initial & boundary conditions

Payoff at τ = 0 (expiration, normalized by K)

• Call: V(x,0) = max(e^x - 1, 0)
• Put: V(x,0) = max(1 - e^x, 0)

Boundary conditions at each time step

• Deep ITM call: V ≈ e^-qτe^x - e^-rτ
• Deep OTM: V ≈ 0
• For puts: swap ITM/OTM conditions

All values normalized by strike K for numerical stability. Multiply by K to recover dollar prices.

4) Crank–Nicolson finite differences

Discretization

• Space: x_i=x_min+iΔx, i=0..M
• Time: τ_n=nΔτ, march backward from n=0 (payoff) to valuation time.

Interior coefficients at node i

λ_i = a_i Δτ/Δx², ν_i = b_i Δτ/(2Δx), a_i=½σ²(x_i,τ), b_i=(r-q)-½σ²(x_i,τ)

• α_i = -½λ_i + ½ν_i
• β_i = 1 + ½ rΔτ + λ_i
• γ_i = -½λ_i - ½ν_i

RHS_i = (1 - ½ rΔτ - λ_i) V_iⁿ + (½λ_i + ½ν_i) V_i-1ⁿ + (½λ_i - ½ν_i) V_i+1ⁿ

Stability & accuracy

• Rannacher start: 1–2 backward-Euler steps to damp payoff kink.
• Vol floor: clamp σ ≥ 0.001.
• Ensure x_min,x_max far enough to minimize boundary influence.

5) How the PDE becomes a matrix–vector linear system

Crank–Nicolson averages the spatial operator between time levels n and n+1. Discretizing interior nodes i=1..M-1 yields a tridiagonal relation:

α_i V_i-1ⁿ⁺¹ + β_i V_iⁿ⁺¹ + γ_i V_i+1ⁿ⁺¹ = RHS_i(Vⁿ).

Vector notation

Collect interior values into vectors Vⁿ = [V₁ⁿ,...,V_M-1ⁿ]^T and Vⁿ⁺¹ of the same size. Define three diagonal vectors α,β,γ from the coefficients above.

Left-hand matrix A and right-hand matrix D

Form a tridiagonal matrix A ∈ ℝ^(M-1)×(M-1) with bands α,β,γ:

A = tridiag(lower = α2..αM-1,
            diag   = β1..βM-1,
            upper  = γ1..γM-2)

Similarly, the RHS can be written as DVⁿ plus boundary contributions bⁿ, where D is another tridiagonal matrix with bands

• α̃_i = ½λ_i + ½ν_i
• β̃_i = 1 - ½ rΔτ - λ_i
• γ̃_i = ½λ_i - ½ν_i

D = tridiag(lower = α̃2..α̃M-1,
            diag   = β̃1..β̃M-1,
            upper  = γ̃1..γ̃M-2)

Boundary vector bⁿ

The first and last interior equations reference the boundaries V₀ⁿ⁺¹, V_Mⁿ⁺¹ (on LHS) and V₀ⁿ, V_Mⁿ (on RHS). Move known boundary values to a vector bⁿ:

• Row 1 (i=1) gets -α₁ V₀ⁿ⁺¹ on LHS and +(α̃₁ V₀ⁿ) on RHS.
• Row M−1 (i=M−1) gets -γ_M-1 V_Mⁿ⁺¹ on LHS and +(γ̃_M-1 V_Mⁿ) on RHS.

Since boundaries are imposed from asymptotics each step, the terms above are known numbers.

Final linear system

A Vⁿ⁺¹ = D Vⁿ + bⁿ.

Solve for Vⁿ⁺¹ using the Thomas algorithm (tri-diagonal solver). Repeat for each time step, marching backward.

Thomas algorithm (outline)

# Given tridiagonal A with (lower l, diag d, upper u) and RHS y:
# 1) Forward sweep: modify coefficients in-place
for i in 2..N:
    w = l[i] / d[i-1]
    d[i] = d[i] - w * u[i-1]
    y[i] = y[i] - w * y[i-1]

# 2) Back substitution
x[N] = y[N] / d[N]
for i in N-1..1:
    x[i] = (y[i] - u[i] * x[i+1]) / d[i]

Matrix assembly summary

1. Compute a_i, b_i at each interior node (often at half time-step).
2. Form λ_i, ν_i then bands α,β,γ (for A) and α̃,β̃,γ̃ (for D).
3. Build bⁿ from boundary values at steps n and n+1.
4. Solve AVⁿ⁺¹ = DVⁿ + bⁿ.

For puts, the same construction applies; only payoff and boundary formulas change.

6) Diagnostics & plotting

• spread = ask − bid
• diff = model − mid, abs_diff = |model − mid|
• within_bid_ask flag (model in [bid, ask])

plot_bid_ask_check(df,
                   price_col="dollar_option_price",
                   mid_col="mid", bid_col="bid", ask_col="ask",
                   y_col="spread")

7) Practical notes & pitfalls

• Normalize columns early (dte, dollar_option_price).
• Ensure SVI coverage across the PDE domain; control extrapolation.
• Use a small vol floor; consider Rannacher start if oscillations appear.
• Place x_min,x_max far enough for negligible boundary impact around quoted strikes.

Appendix — key formulas

Black–Scholes

• k = ℓn(K/S)
• d₁ = Φ⁻¹(Δ_c) (calls), or Φ⁻¹(Δ_p + 1) (puts)
• k = (r - q - ½σ²) T + σ √T · d₁

CN coefficients (interior)

• α_i = −½ λ_i + ½ ν_i
• β_i = 1 + ½ r Δτ + λ_i
• γ_i = −½ λ_i − ½ ν_i
• α̃_i = ½ λ_i + ½ ν_i, β̃_i = 1 − ½ r Δτ − λ_i, γ̃_i = ½ λ_i − ½ ν_i

Matrix form: A Vⁿ⁺¹ = D Vⁿ + bⁿ with tridiagonal A,D and boundary vector bⁿ.