Numerical Methods for Derivative Pricing and Calibration

This notebook develops and compares several numerical methods used in quantitative finance for pricing European options.

We proceed progressively:

1. The Black–Scholes model in continuous time
2. Closed-form pricing and Greeks
3. Implied volatility inversion
4. Monte Carlo pricing and variance reduction
5. Finite difference methods for the Black–Scholes PDE
6. Generic Parabolic PDE Form
7. Calibration on a synthetic implied volatility surface

The objective is to combine mathematical rigor with clear numerical implementation and empirical validation.

1. The Black–Scholes Model in Continuous Time

We consider a filtered probability space

$$ (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P}) $$

The underlying asset price process $S_t$ is assumed to follow a geometric Brownian motion under the physical measure $\mathbb{P}$:

$$ dS_t = \mu S_t \, dt + \sigma S_t \, dW_t $$

where:

• $\mu \in \mathbb{R}$ is the drift,
• $\sigma > 0$ is the volatility,
• $W_t$ is a standard Brownian motion.

Under standard no-arbitrage assumptions, there exists an equivalent martingale measure $\mathbb{Q}$ such that discounted asset prices are martingales.

Under the risk-neutral measure $\mathbb{Q}$, the dynamics become:

$$ dS_t = r S_t \, dt + \sigma S_t \, dW_t^{\mathbb{Q}} $$

where $r$ denotes the constant risk-free rate.

The stochastic differential equation admits the explicit solution:

$$ S_T = S_0 \exp \left( \left( r - \frac{1}{2}\sigma^2 \right) T + \sigma \sqrt{T} Z \right) $$

with

$$ Z \sim \mathcal{N}(0,1) $$

Hence, under $\mathbb{Q}$, the log-price is normally distributed:

$$ \ln S_T \sim \mathcal{N} \left( \ln S_0 + \left( r - \frac{1}{2}\sigma^2 \right)T, \sigma^2 T \right) $$

1.1 Risk-Neutral Valuation and Option Pricing

Let $\Phi(S_T)$ denote the payoff of a European contingent claim at maturity $T$.

Under the risk-neutral measure $\mathbb{Q}$, the arbitrage-free price at time $t$ is given by the discounted conditional expectation:

$$ V(S_t,t) = e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}} \left[ \Phi(S_T) \mid S_t \right] $$

In particular, for a European call option with strike $K$:

$$ \Phi(S_T) = (S_T - K)^+ $$

and for a European put option:

$$ \Phi(S_T) = (K - S_T)^+ $$

Since $S_T$ is lognormally distributed under $\mathbb{Q}$, the expectation can be computed in closed form, leading to the Black–Scholes pricing formulas.

The pricing problem can therefore be approached in two equivalent ways:

• via probabilistic valuation (risk-neutral expectation),
• via partial differential equations (Black–Scholes PDE).

Both approaches will be studied and compared.

2. Closed-Form Black–Scholes Pricing

We now state the closed-form formulas for European options in the Black–Scholes model.

Define

$$ d_1 = \frac{\ln(S_0/K) + \left(r + \frac{1}{2}\sigma^2\right)T}{\sigma\sqrt{T}}, \qquad d_2 = d_1 - \sigma\sqrt{T} $$

Let $N(\cdot)$ denote the cumulative distribution function of a standard normal variable.

The Black–Scholes price of a European call option is:

$$ C_{\mathrm{BS}}(S_0,K,T,r,\sigma) = S_0\,N(d_1) - K e^{-rT} N(d_2) $$

and the price of a European put option is:

$$ P_{\mathrm{BS}}(S_0,K,T,r,\sigma) = K e^{-rT} N(-d_2) - S_0\,N(-d_1) $$

These formulas will serve as a benchmark to validate numerical methods (Monte Carlo and finite differences).

2.1 Greeks (Closed-Form)

In addition to option prices, the Black–Scholes model provides closed-form expressions for the main sensitivities with respect to model parameters and market inputs.

Let $N(\cdot)$ be the standard normal CDF and let $\varphi(\cdot)$ denote the standard normal PDF:

$$ \varphi(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2} $$

Using the notation introduced previously ($d_1$, $d_2$), the Greeks of a European call option are:

Delta $$ \Delta_C = \frac{\partial C}{\partial S_0} = N(d_1) $$

Gamma $$ \Gamma = \frac{\partial^2 C}{\partial S_0^2} = \frac{\varphi(d_1)}{S_0 \sigma \sqrt{T}} $$

Vega $$ \nu = \frac{\partial C}{\partial \sigma} = S_0 \varphi(d_1)\sqrt{T} $$

Theta (one common convention) $$ \Theta_C = \frac{\partial C}{\partial T} = -\frac{S_0 \varphi(d_1)\sigma}{2\sqrt{T}} - r K e^{-rT} N(d_2) $$

These sensitivities will be useful later to:

• understand the behavior of option prices,
• implement implied volatility inversion (through the Vega in Newton's method),
• validate numerical methods against analytical benchmarks.

import sys
from pathlib import Path

# Allow importing from ../src
SRC_DIR = Path("..") / "src"
if str(SRC_DIR.resolve()) not in sys.path:
    sys.path.append(str(SRC_DIR.resolve()))

from bs import bs_price, bs_greeks

from pathlib import Path
import matplotlib.pyplot as plt
FIG_DIR = Path("..") / "reports" / "figures"
FIG_DIR.mkdir(parents=True, exist_ok=True)

def savefig(name):
    path = FIG_DIR / f"{name}.png"
    plt.tight_layout()
    plt.savefig(path, dpi=200)

import numpy as np
import matplotlib.pyplot as plt

# Parameters
S0 = 100
K = 100
T = 1.0
r = 0.05
sigma = 0.2

call_price = bs_price(S0, K, T, r, sigma, option="call")
put_price = bs_price(S0, K, T, r, sigma, option="put")

print("Call price:", round(call_price, 4))
print("Put price :", round(put_price, 4))

# Price as a function of S
S_grid = np.linspace(20, 180, 200)

call_values = bs_price(S_grid, K, T, r, sigma, option="call")
put_values = bs_price(S_grid, K, T, r, sigma, option="put")

plt.figure()
plt.plot(S_grid, call_values, label="Call")
plt.plot(S_grid, put_values, label="Put")
plt.axvline(K, linestyle="--", color="black", alpha=0.6)
plt.title("Black–Scholes Price as a Function of S")
plt.xlabel("Underlying price S")
plt.ylabel("Option price")
plt.legend()

savefig("bs_price_vs_S")
plt.show()

Call price: 10.4506
Put price : 5.5735

greeks_call = bs_greeks(S_grid, K, T, r, sigma, option="call")

delta = greeks_call["delta"]
gamma = greeks_call["gamma"]
vega = greeks_call["vega"]
theta = greeks_call["theta"]

fig, axes = plt.subplots(1, 4, figsize=(15, 4))

axes[0].plot(S_grid, delta)
axes[0].set_title("Delta")
axes[0].set_xlabel("S")
axes[0].set_ylabel("Delta")

axes[1].plot(S_grid, gamma)
axes[1].set_title("Gamma")
axes[1].set_xlabel("S")
axes[1].set_ylabel("Gamma")

axes[2].plot(S_grid, vega)
axes[2].set_title("Vega")
axes[2].set_xlabel("S")
axes[2].set_ylabel("Vega")

axes[3].plot(S_grid, theta)
axes[3].set_title("Theta")
axes[3].set_xlabel("S")
axes[3].set_ylabel("Theta")

plt.tight_layout()
savefig("bs_greeks_vs_S")
plt.show()

Interpretation of the Greeks

The plots above illustrate how the sensitivities evolve with the underlying price.

• Delta:
  Around-the-money, Delta is close to 0.5, meaning that a small increase of 1 unit in $S$ increases the call price by approximately 0.5 units.
  As $S$ becomes large, Delta approaches 1, meaning the option behaves like the underlying asset.
  As $S$ becomes small, Delta approaches 0.

• Gamma:
  Gamma is highest near-the-money.
  This means Delta changes most rapidly when the option is close to the strike.
  Deep in-the-money or out-of-the-money options have very small Gamma.

• Vega:
  Vega is maximal near-the-money.
  This indicates that option prices are most sensitive to volatility when $S \approx K$.
  Far from the strike, volatility has much less impact.

• Theta:
  Theta is typically negative for calls.
  This reflects time decay: as maturity approaches, the option loses time value.
  The decay is strongest near-the-money.

3. Implied Volatility Inversion

In practice, option prices are observed in the market, while the volatility parameter is not directly tradable. Given a market price $C_{\text{mkt}}$, the implied volatility $\sigma_{\text{imp}}$ is defined as the solution of:

$$ C_{\mathrm{BS}}(S_0, K, T, r, \sigma_{\text{imp}}) = C_{\text{mkt}} $$

Equivalently, define the function

$$ f(\sigma) = C_{\mathrm{BS}}(S_0, K, T, r, \sigma) - C_{\text{mkt}} $$

and solve:

$$ f(\sigma) = 0 $$

This is a one-dimensional root-finding problem. We will implement and compare:

• the bisection method (robust but slower),
• the Newton method (faster but sensitive to initialization).

3.1 Bisection Method

Assume that $\sigma_{\min} < \sigma_{\max}$ are such that $f(\sigma_{\min})$ and $f(\sigma_{\max})$ have opposite signs. Bisection iteratively halves the interval until convergence.

3.2 Newton Method

Newton's update rule is:

$$ \sigma_{n+1} = \sigma_n - \frac{f(\sigma_n)}{f'(\sigma_n)} $$

In our case,

$$ f'(\sigma) = \frac{\partial C_{\mathrm{BS}}}{\partial \sigma}(S_0, K, T, r, \sigma) $$

which is precisely the Vega. Newton's method is therefore efficient when Vega is not too small.

from bs import implied_vol_bisection, implied_vol_newton

# Synthetic market price generated from a known volatility
true_sigma = 0.25
market_call = bs_price(S0, K, T, r, true_sigma, option="call")

sigma_bis = implied_vol_bisection(S0, K, T, r, market_call, option="call",
                                  sigma_min=1e-6, sigma_max=5.0, tol=1e-10, max_iter=200)

sigma_newton = implied_vol_newton(S0, K, T, r, market_call, option="call",
                                  sigma_init=1, tol=1e-10, max_iter=50)

print("True sigma   :", true_sigma)
print("Bisection    :", sigma_bis)
print("Newton       :", sigma_newton)
print("Abs error (bisection):", abs(sigma_bis - true_sigma))
print("Abs error (newton)   :", abs(sigma_newton - true_sigma))

True sigma   : 0.25
Bisection    : 0.250000000002051
Newton       : 0.24999999999999986
Abs error (bisection): 2.051026015692514e-12
Abs error (newton)   : 1.3877787807814457e-16

from bs import implied_vol_bisection_path, implied_vol_newton_path

sig_bis, err_bis = implied_vol_bisection_path(
    S0, K, T, r, market_call, option="call",
    sigma_min=1e-6, sigma_max=5.0, max_iter=60
)

sig_new, err_new = implied_vol_newton_path(
    S0, K, T, r, market_call, option="call",
    sigma_init=0.2, max_iter=30
)

plt.figure()
plt.plot(err_bis, label="Bisection")
plt.plot(err_new, label="Newton")
plt.yscale("log")
plt.title("Implied Volatility Convergence (Absolute Pricing Error)")
plt.xlabel("Iteration")
plt.ylabel("Absolute pricing error (log scale)")
plt.legend()

savefig("implied_vol_convergence")
plt.show()

Interpretation of the Convergence Plot

The figure compares the convergence of bisection and Newton's method using the absolute pricing error $|C_{\mathrm{BS}}(\sigma_n) - C_{\text{mkt}}|$ at each iteration (log scale).

• Bisection converges steadily but slowly.
  Each iteration halves the interval that brackets the solution, so the error decreases at a predictable rate. This robustness is obtained at the cost of speed.

• Newton's method converges much faster when the initial guess is reasonable and Vega is not too small.
  Near the solution, Newton exhibits quadratic convergence: the number of correct digits typically grows very quickly. This explains the sharp drop observed in the first few iterations.

In practice, a common robust approach is to start with a bracketing method (such as bisection) and then switch to Newton once the iterate is close enough to the root.

4. Monte Carlo Pricing

Under the risk-neutral measure, the price of a European call is

$$ C = e^{-rT} \mathbb{E}^{\mathbb{Q}}[(S_T - K)^+] $$

Using the explicit solution of the Black–Scholes model:

$$ S_T = S_0 \exp\left(\left(r - \frac{1}{2}\sigma^2\right)T + \sigma \sqrt{T} Z \right), \quad Z \sim \mathcal{N}(0,1) $$

A Monte Carlo estimator with $N$ simulations is:

$$ \hat{C}_N = e^{-rT} \frac{1}{N} \sum_{i=1}^{N} (S_T^{(i)} - K)^+ $$

By the Law of Large Numbers:

$$ \hat{C}_N \to C $$

and by the Central Limit Theorem:

$$ \sqrt{N}(\hat{C}_N - C) \xrightarrow{d} \mathcal{N}(0, \sigma^2_{\text{MC}}) $$

Thus, the standard error decreases at rate:

$$ \mathcal{O}\left(\frac{1}{\sqrt{N}}\right) $$

We will implement:

• a plain Monte Carlo estimator,
• variance reduction techniques (antithetic variates).

Plain Monte Carlo vs Antithetic Variates

We consider two Monte Carlo estimators.

Plain Monte Carlo

We simulate $N$ independent standard normal variables:

$$ Z_1, \dots, Z_N \sim \mathcal{N}(0,1) $$

and compute:

$$ \widehat{C}_N = e^{-rT}\frac{1}{N}\sum_{i=1}^N \Phi(S_T^{(i)}) $$

This estimator is unbiased and its standard error decreases at rate:

$$ \mathcal{O}(N^{-1/2}) $$

Antithetic Variates

Since $Z \sim \mathcal{N}(0,1)$ implies $-Z \sim \mathcal{N}(0,1)$, we can use pairs $(Z, -Z)$ to construct the estimator:

$$ \frac{1}{2}\left[\Phi(S_T(Z)) + \Phi(S_T(-Z))\right] $$

This often reduces variance because the two terms are negatively correlated.

We compare both estimators in terms of:

• estimated price,
• standard error,
• convergence behavior.

from mc import mc_price_plain, mc_price_antithetic, simulate_gbm_paths

# Benchmark (analytical Black–Scholes)
bs_call = bs_price(S0, K, T, r, sigma, option="call")

# Monte Carlo settings
n_paths = 200_000
rng = np.random.default_rng(42)

mc_call_plain, se_plain = mc_price_plain(S0, K, T, r, sigma, n_paths, option="call", rng=rng)

rng = np.random.default_rng(42)  # same seed for fair comparison
mc_call_anti, se_anti = mc_price_antithetic(S0, K, T, r, sigma, n_paths, option="call", rng=rng)

def ci95(est, se):
    return (est - 1.96 * se, est + 1.96 * se)

print("Black–Scholes call:", bs_call)

print("\nPlain MC call:", mc_call_plain)
print("Plain MC SE  :", se_plain)
print("Plain MC 95% CI:", ci95(mc_call_plain, se_plain))

print("\nAntithetic MC call:", mc_call_anti)
print("Antithetic MC SE  :", se_anti)
print("Antithetic MC 95% CI:", ci95(mc_call_anti, se_anti))

Black–Scholes call: 10.450583572185565

Plain MC call: 10.463413536194789
Plain MC SE  : 0.033117555441239864
Plain MC 95% CI: (np.float64(10.39850312752996), np.float64(10.528323944859618))

Antithetic MC call: 10.476284520066631
Antithetic MC SE  : 0.033089502857212114
Antithetic MC 95% CI: (np.float64(10.411429094466495), np.float64(10.541139945666767))

Simulating Sample Paths

So far we only simulated the terminal value $S_T$ to price the option. To build intuition, we now simulate full paths $(S_t)_{0 \le t \le T}$ under the risk-neutral dynamics:

$$ dS_t = r S_t\,dt + \sigma S_t\,dW_t $$

Using the exact discretization on a uniform time grid $t_k = k\Delta t$:

$$ S_{t_{k+1}} = S_{t_k}\exp\left(\left(r - \frac{1}{2}\sigma^2\right)\Delta t + \sigma\sqrt{\Delta t}\,Z_k\right), \quad Z_k \sim \mathcal{N}(0,1) $$

This is not needed for plain European pricing (since $S_T$ is enough), but it becomes essential for:

• path-dependent options (barriers, Asians),
• visual diagnostics and intuition about volatility and drift.

from mc import simulate_gbm_paths

n_steps = 252
n_plot_paths = 30
rng = np.random.default_rng(123)

t_grid, S_paths = simulate_gbm_paths(S0, T, r, sigma, n_steps=n_steps, n_paths=n_plot_paths, rng=rng)

plt.figure()
for i in range(n_plot_paths):
    plt.plot(t_grid, S_paths[i, :])

plt.title("Sample Risk-Neutral GBM Paths")
plt.xlabel("Time t")
plt.ylabel("Underlying price S(t)")

savefig("mc_gbm_paths")
plt.show()

# Large simulation to study terminal distribution
n_paths_large = 200_000
rng = np.random.default_rng(321)

_, S_paths_large = simulate_gbm_paths(
    S0, T, r, sigma,
    n_steps=1,          # only need S_T
    n_paths=n_paths_large,
    rng=rng
)

S_T = S_paths_large[:, -1]

plt.figure()
plt.hist(S_T, bins=100, density=True, alpha=0.6, label="MC histogram")

# Theoretical lognormal density
x = np.linspace(min(S_T), max(S_T), 500)

mu = np.log(S0) + (r - 0.5 * sigma**2) * T
var = sigma**2 * T
std = np.sqrt(var)

from scipy.stats import lognorm
theoretical_pdf = lognorm.pdf(x, s=std, scale=np.exp(mu))

plt.plot(x, theoretical_pdf, linewidth=2, label="Theoretical density")

plt.title("Distribution of $S_T$ (Monte Carlo vs Theory)")
plt.xlabel("$S_T$")
plt.ylabel("Density")
plt.legend()

savefig("mc_terminal_distribution")
plt.show()

# Convergence study using estimated standard error

N_values = np.logspace(2, 5, 50, dtype=int)

se_plain_list = []
se_anti_list = []

for N in N_values:

    rng = np.random.default_rng(123)
    _, se_plain = mc_price_plain(S0, K, T, r, sigma, N, option="call", rng=rng)

    rng = np.random.default_rng(123)
    _, se_anti = mc_price_antithetic(S0, K, T, r, sigma, N, option="call", rng=rng)

    se_plain_list.append(se_plain)
    se_anti_list.append(se_anti)


plt.figure()
plt.loglog(N_values, se_plain_list, marker="o", label="Plain MC SE")
plt.loglog(N_values, se_anti_list, marker="o", label="Antithetic MC SE")

# theoretical slope
ref = se_plain_list[0] * np.sqrt(N_values[0]) / np.sqrt(N_values)
plt.loglog(N_values, ref, linestyle="--", label="$O(N^{-1/2})$")

plt.xlabel("Number of paths N")
plt.ylabel("Standard error")
plt.title("Monte Carlo Standard Error Convergence")
plt.legend()

savefig("mc_convergence_rate_clean")
plt.show()

Monte Carlo Pricing — Summary and Interpretation

We have implemented and analyzed Monte Carlo pricing under the risk-neutral Black–Scholes framework.

Main results:

• The simulated terminal distribution of $S_T$ matches the theoretical lognormal density.
• The Monte Carlo estimator is unbiased.
• The empirical standard error decreases at the theoretical rate:

$$ \text{SE} \sim \mathcal{O}(N^{-1/2}) $$

• Antithetic variates reduce variance by introducing negative correlation between paired samples.

In practice:

• Plain Monte Carlo is simple but noisy.
• Antithetic Monte Carlo improves stability at almost no additional cost.
• The $N^{-1/2}$ convergence rate is slow compared to PDE-based methods, which motivates deterministic numerical schemes.

Monte Carlo is therefore:

• Flexible (dimension-independent),
• Easy to implement,
• Well-suited for path-dependent derivatives,
• But relatively slow in terms of convergence.

5. Finite Difference Methods for the Black–Scholes PDE

Monte Carlo methods are flexible but converge slowly at rate:

$$ \mathcal{O}(N^{-1/2}) $$

We now turn to deterministic numerical methods based on partial differential equations.

Under the risk-neutral measure, the Black–Scholes PDE is:

$$ \frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + r S \frac{\partial V}{\partial S} - rV = 0 $$

with terminal condition:

$$ V(S,T) = \Phi(S) $$

For a European call:

$$ \Phi(S) = (S-K)^+ $$

Our goal is to approximate the solution on a discrete grid:

• Time grid: $t_k = k\Delta t$
• Space grid: $S_i = i\Delta S$

We will implement:

• Explicit scheme
• Implicit scheme
• Crank–Nicolson scheme

and compare:

• Accuracy
• Stability
• Convergence
• Computational cost

5.1 Grid Construction

To solve the Black–Scholes PDE numerically, we discretize both time and space.

We define:

• Time grid:

$$ t_k = k \Delta t, \quad k = 0,1,\dots,M $$

where

$$ \Delta t = \frac{T}{M} $$

• Space grid:

$$ S_i = i \Delta S, \quad i = 0,1,\dots,N $$

where

$$ \Delta S = \frac{S_{\max}}{N} $$

We approximate the solution by:

$$ V_i^k \approx V(S_i, t_k) $$

Our objective is to compute the solution backward in time, starting from the terminal condition:

$$ V_i^M = \Phi(S_i) $$

# Grid parameters
S_max = 3 * K      # upper bound in space
N = 200            # number of space steps
M = 200            # number of time steps

dS = S_max / N
dt = T / M

# Space and time grids
S_grid = np.linspace(0, S_max, N + 1)
t_grid = np.linspace(0, T, M + 1)

# Initialize solution matrix
V = np.zeros((N + 1, M + 1))

# Terminal condition (European call)
V[:, -1] = np.maximum(S_grid - K, 0)

import matplotlib.pyplot as plt
import numpy as np

# Small illustrative grid (just for visualization)
N_demo = 8
M_demo = 6

S_demo = np.linspace(0, 8, N_demo)
t_demo = np.linspace(0, 6, M_demo)

S_mesh_demo, t_mesh_demo = np.meshgrid(S_demo, t_demo)

plt.figure(figsize=(6, 5))
plt.scatter(S_mesh_demo, t_mesh_demo, color="black")

# Highlight a central node
i0 = 4
k0 = 3

plt.scatter(S_demo[i0], t_demo[k0], color="red", s=100)

# Neighboring points (Crank-Nicolson style stencil)
neighbors = [
    (i0-1, k0),
    (i0+1, k0),
    (i0, k0-1),
    (i0, k0+1)
]

for (i, k) in neighbors:
    plt.scatter(S_demo[i], t_demo[k], color="blue", s=60)
    plt.plot([S_demo[i0], S_demo[i]],
             [t_demo[k0], t_demo[k]],
             linestyle="--")

plt.title("Finite Difference Grid and Local Stencil")
plt.xlabel("Space index i (Underlying price S)")
plt.ylabel("Time index k")

plt.gca().invert_yaxis()  # time goes backward visually
plt.grid(True)
savefig("pde_grid_stencil")
plt.text(S_demo[i0], t_demo[k0]-0.3, r"$V_i^k$", color="red", ha="center")

plt.show()

5.2 Boundary Conditions

Since the Black–Scholes PDE is posed on $S \in (0,\infty)$, we truncate the domain to $[0, S_{\max}]$. To obtain a well-defined numerical scheme, we must specify boundary conditions at $S=0$ and $S=S_{\max}$.

For a European call option:

• At $S=0$, the option is worthless:

$$ V(0,t) = 0 $$

• For large $S$, the call behaves like $S - K e^{-r(T-t)}$:

$$ V(S_{\max}, t) \approx S_{\max} - K e^{-r(T-t)} $$

These boundary conditions will be enforced at each time step.

# Boundary conditions across time
V[0, :] = 0.0
V[-1, :] = S_max - K * np.exp(-r * (T - t_grid))

5.2.1 Derivation of boundary conditions (call / put)

The Black–Scholes PDE is posed on an unbounded domain (S \in (0,\infty)). In practice we truncate the grid to ([0,S_{\max}]), so we must impose boundary conditions consistent with the option payoff.

Left boundary: (S=0)

• Call: if (S=0), the option is worthless for all (t \le T): [ V(0,t)=0. ]

• Put: if (S=0), the put is certainly exercised at maturity, hence: [ V(0,t)=K e^{-r(T-t)}. ]

Right boundary: (S=S_{\max}) (large underlying)

For large (S), a European call behaves like its intrinsic value discounted for the strike: [ V(S,t) \sim S - K e^{-r(T-t)}. ] Thus we impose: [ V(S_{\max},t) = S_{\max} - K e^{-r(T-t)}. ]

For a European put, when (S) is large it is almost surely out-of-the-money: [ V(S_{\max},t) \approx 0. ]

These conditions are enforced at each time step and ensure the finite-difference system is well-posed.

5.3 Finite Difference Approximations

For interior grid points $S_i$ (with $i=1,\dots,N-1$), we approximate derivatives by finite differences.

Second derivative (central difference):

$$ \frac{\partial^2 V}{\partial S^2}(S_i,t_k) \approx \frac{V_{i+1}^k - 2V_i^k + V_{i-1}^k}{(\Delta S)^2} $$

First derivative (central difference):

$$ \frac{\partial V}{\partial S}(S_i,t_k) \approx \frac{V_{i+1}^k - V_{i-1}^k}{2\Delta S} $$

Time derivative (backward stepping in time):

$$ \frac{\partial V}{\partial t}(S_i,t_k) \approx \frac{V_i^{k+1} - V_i^{k}}{\Delta t} $$

We will compute the solution backward in time from $t=T$ to $t=0$.

5.3.1 Grid notation and backward time stepping

We use a uniform grid in space and time: [ S_i = i\Delta S,\quad i=0,\dots,N, \qquad t_k = k\Delta t,\quad k=0,\dots,M, ] and denote [ V_i^k \approx V(S_i,t_k). ]

Because the terminal condition is given at maturity, [ V(S,T)=\Phi(S), ] we compute the solution backward in time, i.e. from (t_M=T) down to (t_0=0). In practice, each scheme will express the unknown values at time (t_k) using already known values at time (t_{k+1}).

5.4 Explicit Scheme

We rewrite the Black–Scholes PDE as:

$$ \frac{\partial V}{\partial t} = -\frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} - r S \frac{\partial V}{\partial S} + rV $$

Using the finite difference approximations from the previous section, we obtain an explicit update formula:

$$ V_i^{k} = a_i V_{i-1}^{k+1} + b_i V_i^{k+1} + c_i V_{i+1}^{k+1} $$

for $i=1,\dots,N-1$, where the coefficients are:

$$ a_i = \frac{1}{2}\Delta t\left(\sigma^2 i^2 - r i\right), \qquad b_i = 1 - \Delta t\left(\sigma^2 i^2 + r\right), \qquad c_i = \frac{1}{2}\Delta t\left(\sigma^2 i^2 + r i\right) $$

This scheme is simple but conditionally stable, which means the time step $\Delta t$ must be small enough relative to $\Delta S$.

5.4.1 Derivation of the explicit update and CFL-type stability condition

Starting from the Black–Scholes PDE written backward in time: $$ \frac{\partial V}{\partial t} = -\frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} - rS \frac{\partial V}{\partial S} + rV, $$ we substitute the finite-difference approximations at $(S_i,t_k)$: $$ \frac{V_i^{k+1}-V_i^{k}}{\Delta t} = -\frac{1}{2}\sigma^2 S_i^2 \frac{V_{i+1}^{k+1}-2V_i^{k+1}+V_{i-1}^{k+1}}{(\Delta S)^2} - rS_i \frac{V_{i+1}^{k+1}-V_{i-1}^{k+1}}{2\Delta S} + rV_i^{k+1}. $$

Rearranging yields an explicit recursion of the form: $$ V_i^{k} = a_i V_{i-1}^{k+1} + b_i V_i^{k+1} + c_i V_{i+1}^{k+1}, \qquad i=1,\dots,N-1, $$ with coefficients depending on $\Delta t$, $\Delta S$, and $S_i$.

Stability (CFL-type condition)

The explicit scheme is conditionally stable: $\Delta t$ must be sufficiently small relative to $\Delta S$. A common sufficient condition is obtained by comparing the diffusion term to the explicit heat equation stability bound: $$ \boxed{ \Delta t \;\lesssim\; \frac{(\Delta S)^2}{\sigma^2 S_{\max}^2} } $$ (up to a constant factor depending on the discretization convention).

This explains why refining the spatial grid (smaller $\Delta S$) forces a quadratically smaller time step to keep the explicit method stable.

from pde import solve_bs_explicit
from bs import bs_price

# Solve PDE (explicit scheme)
S_grid_fd, t_grid_fd, V_fd = solve_bs_explicit(
    K=K, T=T, r=r, sigma=sigma,
    S_max_mult=3.0, N=200, M=1500, option="call"
)


# Extract price at t=0 by interpolating on the S grid
V_t0 = V_fd[:, 0]
price_fd_S0 = float(np.interp(S0, S_grid_fd, V_t0))

# Analytical Black–Scholes benchmark
price_bs = float(bs_price(S0, K, T, r, sigma, option="call"))

print("Explicit FD price at S0:", price_fd_S0)
print("Black–Scholes price      :", price_bs)
print("Absolute error           :", abs(price_fd_S0 - price_bs))

Explicit FD price at S0: 10.45513019275099
Black–Scholes price      : 10.450583572185565
Absolute error           : 0.004546620565426096

Interpretation.
The explicit finite-difference scheme is conditionally stable. With a coarse time grid (large $\Delta t$), the solution may explode numerically and produce meaningless values.

After refining the time discretization (increasing $M$, hence reducing $\Delta t$), the scheme becomes stable and the price at $S_0$ matches the analytical Black–Scholes benchmark up to a small discretization error.

# Compare the full curve V(S,0) (FD) vs Black–Scholes closed-form

V_t0_fd = V_fd[:, 0]
V_t0_bs = bs_price(S_grid_fd[1:], K, T, r, sigma, option="call")  # avoid S=0 for log

plt.figure()
plt.plot(S_grid_fd, V_t0_fd, label="Explicit FD (t=0)")
plt.plot(S_grid_fd[1:], V_t0_bs, linestyle="--", label="Black–Scholes (t=0)")

plt.axvline(S0, linestyle="--", alpha=0.6)
plt.title("European Call Price at t=0: Explicit FD vs Black–Scholes")
plt.xlabel("Underlying price S")
plt.ylabel("Option price")
plt.legend()

savefig("pde_explicit_vs_bs_t0")
plt.show()

# Global error metrics at t=0

# Avoid S=0 for BS evaluation
S_eval = S_grid_fd[1:]
V_fd_eval = V_t0_fd[1:]
V_bs_eval = bs_price(S_eval, K, T, r, sigma, option="call")

abs_errors = np.abs(V_fd_eval - V_bs_eval)

max_error = np.max(abs_errors)
rmse = np.sqrt(np.mean((V_fd_eval - V_bs_eval)**2))

print("Max error (sup norm):", max_error)
print("RMSE:", rmse)

Max error (sup norm): 0.0013053334225117674
RMSE: 0.00038071660065808877

import numpy as np
import matplotlib.pyplot as plt

# 3D surface plot of V(S,t) from Explicit scheme
S_mesh_fd, t_mesh_fd = np.meshgrid(S_grid_fd, t_grid_fd, indexing="ij")

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

ax.plot_surface(S_mesh_fd, t_mesh_fd, V_fd,
                linewidth=0, antialiased=True)

ax.set_title("Explicit solution surface V(S,t)")
ax.set_xlabel("Underlying price S")
ax.set_ylabel("Time t")
ax.set_zlabel("Option value V")

ax.view_init(elev=25, azim=-135)

savefig("pde_explicit_surface_V_S_t")
plt.show()

5.5 Implicit Scheme (Backward Euler)

The implicit (Backward Euler) scheme evaluates the spatial derivatives at the unknown time level $t_k$.

Using the same finite-difference approximations in space, the implicit time discretization writes:

$$ \frac{V_i^{k+1} - V_i^{k}}{\Delta t} = -\frac{1}{2}\sigma^2 S_i^2 \frac{\partial^2 V}{\partial S^2}(S_i,t_k) - r S_i \frac{\partial V}{\partial S}(S_i,t_k) + rV_i^{k} $$

This leads to a linear system at each step:

$$ A_i V_{i-1}^{k} + B_i V_i^{k} + C_i V_{i+1}^{k} = V_i^{k+1} $$

for $i=1,\dots,N-1$.

Key point:

• the method is unconditionally stable, but
• we must solve a tridiagonal system at each time step.

5.5.1 Tridiagonal system and coefficients (Backward Euler)

With the implicit (Backward Euler) discretization, the unknown values at time $t_k$ appear inside the spatial derivatives. After substitution of finite differences, we obtain for each interior node $i=1,\dots,N-1$: $$ A_i V_{i-1}^{k} + B_i V_i^{k} + C_i V_{i+1}^{k} = V_i^{k+1}, $$ where the coefficients take the generic form: $$ A_i = -\Delta t\left(\frac{1}{2}\sigma^2 S_i^2\frac{1}{(\Delta S)^2} - \frac{rS_i}{2\Delta S}\right), $$ $$ B_i = 1 + \Delta t\left(\sigma^2 S_i^2\frac{1}{(\Delta S)^2} + r\right), $$ $$ C_i = -\Delta t\left(\frac{1}{2}\sigma^2 S_i^2\frac{1}{(\Delta S)^2} + \frac{rS_i}{2\Delta S}\right). $$

Thus at each time step we solve a tridiagonal linear system: $$ \mathbf{A} V^{k} = V^{k+1}, $$ which can be solved efficiently with the Thomas algorithm in $\mathcal{O}(N)$.

from pde import solve_bs_implicit
from bs import bs_price

# Implicit FD solve (Backward Euler)
S_grid_imp, t_grid_imp, V_imp = solve_bs_implicit(
    K=K, T=T, r=r, sigma=sigma,
    S_max_mult=3.0, N=200, M=200, option="call"
)

# Price at t=0 by interpolation
V_t0_imp = V_imp[:, 0]
price_imp_S0 = float(np.interp(S0, S_grid_imp, V_t0_imp))

# Black–Scholes benchmark
price_bs = float(bs_price(S0, K, T, r, sigma, option="call"))

print("Implicit FD price at S0:", price_imp_S0)
print("Black–Scholes price      :", price_bs)
print("Absolute error           :", abs(price_imp_S0 - price_bs))

Implicit FD price at S0: 10.449196020476712
Black–Scholes price      : 10.450583572185565
Absolute error           : 0.001387551708852186

# Compare the whole curve V(S,0): Implicit FD vs Black–Scholes

V_t0_imp = V_imp[:, 0]

# Avoid S=0 for BS evaluation (log)
S_eval = S_grid_imp[1:]
V_bs_eval = bs_price(S_eval, K, T, r, sigma, option="call")

plt.figure()
plt.plot(S_grid_imp, V_t0_imp, label="Implicit FD (t=0)")
plt.plot(S_eval, V_bs_eval, linestyle="--", label="Black–Scholes (t=0)")

plt.axvline(S0, linestyle="--", alpha=0.6)
plt.title("European Call Price at t=0: Implicit FD vs Black–Scholes")
plt.xlabel("Underlying price S")
plt.ylabel("Option price")
plt.legend()

savefig("pde_implicit_vs_bs_t0")
plt.show()

# Global error metrics at t=0 (Implicit scheme)

S_eval = S_grid_imp[1:]   # avoid S=0
V_fd_eval = V_t0_imp[1:]
V_bs_eval = bs_price(S_eval, K, T, r, sigma, option="call")

abs_errors = np.abs(V_fd_eval - V_bs_eval)

max_error = np.max(abs_errors)
rmse = np.sqrt(np.mean((V_fd_eval - V_bs_eval)**2))

print("Implicit - Max error (sup norm):", max_error)
print("Implicit - RMSE:", rmse)

Implicit - Max error (sup norm): 0.006168423536504264
Implicit - RMSE: 0.002048469421116114

# 3D surface plot of V(S,t) from Implicit scheme
S_mesh_imp, t_mesh_imp = np.meshgrid(S_grid_imp, t_grid_imp, indexing="ij")

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

ax.plot_surface(S_mesh_imp, t_mesh_imp, V_imp,
                linewidth=0, antialiased=True)

ax.set_title("Implicit solution surface V(S,t)")
ax.set_xlabel("Underlying price S")
ax.set_ylabel("Time t")
ax.set_zlabel("Option value V")

ax.view_init(elev=25, azim=-135)

savefig("pde_implicit_surface_V_S_t")
plt.show()

5.6 Crank–Nicolson Scheme

The Crank–Nicolson method is obtained by averaging the explicit and implicit schemes. It corresponds to a trapezoidal discretization in time.

For the Black–Scholes PDE:

$$ \frac{\partial V}{\partial t} = \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + r S \frac{\partial V}{\partial S} - rV $$

The Crank–Nicolson discretization reads:

$$ \frac{V_i^{k+1} - V_i^k}{\Delta t} = \frac{1}{2} \left[ \mathcal{L}V_i^{k+1} + \mathcal{L}V_i^k \right] $$

This leads to a tridiagonal linear system at each time step:

$$ M V^{k} = N V^{k+1} $$

The method is:

• unconditionally stable
• second-order accurate in time and space

5.6.1 Crank–Nicolson as a matrix system $M V^{k} = N V^{k+1}$

Crank–Nicolson averages the explicit and implicit discretizations: $$ \frac{V_i^{k+1}-V_i^k}{\Delta t}=\frac{1}{2}\left(\mathcal{L}V^{k+1}_i+\mathcal{L}V^k_i\right), $$ which leads to a tridiagonal system coupling the vectors (V^k) and (V^{k+1}).

For interior points $i=1,\dots,N-1$, the scheme can be written: $$ \alpha_i V_{i-1}^{k} + \beta_i V_i^{k} + \gamma_i V_{i+1}^{k} = \tilde{\alpha}_i V_{i-1}^{k+1} + \tilde{\beta}_i V_i^{k+1} + \tilde{\gamma}_i V_{i+1}^{k+1}, $$ where the coefficients come from splitting the operator terms between time levels (k) and (k+1).

In matrix form: $$ \mathbf{M}\,V^{k}=\mathbf{N}\,V^{k+1}, $$ with $\mathbf{M}$ and $\mathbf{N}$ tridiagonal matrices, so each time step requires solving a tridiagonal linear system. Crank–Nicolson is typically second-order accurate in time and more accurate than Backward Euler.

from pde import solve_bs_crank_nicolson
from bs import bs_price

S_grid_cn, t_grid_cn, V_cn = solve_bs_crank_nicolson(
    K=K, T=T, r=r, sigma=sigma,
    S_max_mult=3.0, N=200, M=200, option="call"
)

V_t0_cn = V_cn[:, 0]
price_cn_S0 = float(np.interp(S0, S_grid_cn, V_t0_cn))

price_bs = float(bs_price(S0, K, T, r, sigma, option="call"))

print("Crank–Nicolson price at S0:", price_cn_S0)
print("Black–Scholes price       :", price_bs)
print("Absolute error            :", abs(price_cn_S0 - price_bs))

Crank–Nicolson price at S0: 10.454438036590357
Black–Scholes price       : 10.450583572185565
Absolute error            : 0.0038544644047924237

import numpy as np
import matplotlib.pyplot as plt

# 3D surface plot of V(S,t) from Crank–Nicolson
# V_cn has shape (N+1, M+1) with:
#   rows  -> S grid
#   cols  -> t grid

S_mesh, t_mesh = np.meshgrid(S_grid_cn, t_grid_cn, indexing="ij")  # both shape (N+1, M+1)

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

# plot_surface expects 2D arrays
ax.plot_surface(S_mesh, t_mesh, V_cn, linewidth=0, antialiased=True)

ax.set_title("Crank–Nicolson solution surface V(S,t)")
ax.set_xlabel("Underlying price S")
ax.set_ylabel("Time t")
ax.set_zlabel("Option value V")

# a nicer viewing angle
ax.view_init(elev=25, azim=-135)

savefig("pde_cn_surface_V_S_t")
plt.show()

# Plot several time slices V(S, t_k)

plt.figure()

time_indices = np.linspace(0, len(t_grid_cn)-1, 10, dtype=int)

for k in time_indices:
    plt.plot(S_grid_cn, V_cn[:, k], alpha=0.8)

plt.title("Crank-Nicolson: Time evolution V(S,t)")
plt.xlabel("Underlying price S")
plt.ylabel("Option value V")
plt.show()

5.7 Convergence Study

We study the convergence of the finite-difference schemes by increasing the spatial resolution N (and M accordingly).

We compute the pricing error at S₀ for increasing grid sizes.

from pde import solve_bs_crank_nicolson
from bs import bs_price

N_values = [50, 100, 150, 200, 300, 400]
errors = []

for N_test in N_values:
    M_test = N_test  

    Sg, tg, Vcn = solve_bs_crank_nicolson(
        K=K, T=T, r=r, sigma=sigma,
        S_max_mult=3.0, N=N_test, M=M_test,
        option="call"
    )

    V0 = Vcn[:, 0]
    price_fd = float(np.interp(S0, Sg, V0))
    price_exact = float(bs_price(S0, K, T, r, sigma, option="call"))

    errors.append(abs(price_fd - price_exact))

plt.figure()
plt.loglog(N_values, errors, marker="o", label="Crank–Nicolson error")

# theoretical slope ~ O(N^-2)
ref = errors[0] * (np.array(N_values[0]) / np.array(N_values))**2
plt.loglog(N_values, ref, linestyle="--", label=r"$O(N^{-2})$")

plt.xlabel("Number of space steps N")
plt.ylabel("Absolute error at S0")
plt.title("Convergence of Crank–Nicolson Scheme")
plt.legend()

savefig("pde_cn_convergence")
plt.show()

import pandas as pd

results = pd.DataFrame({
    "Method": ["Explicit", "Implicit", "Crank-Nicolson"],
    "Price at S0": [price_fd_S0, price_imp_S0, price_cn_S0],
    "BS Price": [price_bs, price_bs, price_bs],
    "Abs Error": [
        abs(price_fd_S0 - price_bs),
        abs(price_imp_S0 - price_bs),
        abs(price_cn_S0 - price_bs)
    ]
})

results

	Method	Price at S0	BS Price	Abs Error
0	Explicit	10.455130	10.450584	0.004547
1	Implicit	10.449196	10.450584	0.001388
2	Crank-Nicolson	10.454438	10.450584	0.003854

from pde import solve_bs_explicit, solve_bs_implicit, solve_bs_crank_nicolson
from bs import bs_price

N_values = [50, 100, 150, 200, 300, 400]

errors_exp = []
errors_imp = []
errors_cn  = []

for N_test in N_values:

    M_test = N_test  

    # Explicit
    Sg, tg, Vexp = solve_bs_explicit(
        K=K, T=T, r=r, sigma=sigma,
        S_max_mult=3.0, N=N_test, M=M_test, option="call"
    )
    price_exp = float(np.interp(S0, Sg, Vexp[:, 0]))

    # Implicit
    Sg, tg, Vimp = solve_bs_implicit(
        K=K, T=T, r=r, sigma=sigma,
        S_max_mult=3.0, N=N_test, M=M_test, option="call"
    )
    price_imp = float(np.interp(S0, Sg, Vimp[:, 0]))

    # Crank-Nicolson
    Sg, tg, Vcn = solve_bs_crank_nicolson(
        K=K, T=T, r=r, sigma=sigma,
        S_max_mult=3.0, N=N_test, M=M_test, option="call"
    )
    price_cn = float(np.interp(S0, Sg, Vcn[:, 0]))

    price_exact = float(bs_price(S0, K, T, r, sigma, option="call"))

    errors_exp.append(abs(price_exp - price_exact))
    errors_imp.append(abs(price_imp - price_exact))
    errors_cn.append(abs(price_cn - price_exact))

plt.figure()

plt.loglog(N_values, errors_exp, marker="o", label="Explicit")
plt.loglog(N_values, errors_imp, marker="o", label="Implicit")
plt.loglog(N_values, errors_cn, marker="o", label="Crank-Nicolson")

plt.xlabel("Number of space steps N")
plt.ylabel("Absolute error at S0")
plt.title("Convergence comparison of FD schemes")
plt.legend()

savefig("pde_convergence_comparison")
plt.show()

C:\Users\mathi\Desktop\GIT\numerical-methods-finance\src\pde.py:65: RuntimeWarning: overflow encountered in multiply
  V[i, k] = a * V[i - 1, k + 1] + b * V[i, k + 1] + c * V[i + 1, k + 1]
C:\Users\mathi\Desktop\GIT\numerical-methods-finance\src\pde.py:65: RuntimeWarning: overflow encountered in add
  V[i, k] = a * V[i - 1, k + 1] + b * V[i, k + 1] + c * V[i + 1, k + 1]
C:\Users\mathi\Desktop\GIT\numerical-methods-finance\src\pde.py:65: RuntimeWarning: invalid value encountered in add
  V[i, k] = a * V[i - 1, k + 1] + b * V[i, k + 1] + c * V[i + 1, k + 1]

Interpretation: stability vs convergence

In this experiment we increase the spatial resolution $N$ (and set $M=N$ for the time grid). This choice is not appropriate for the explicit scheme: the explicit finite-difference method is conditionally stable and requires a sufficiently small time step $\Delta t$ relative to $\Delta S$ (CFL-type condition).

As $N$ increases, $\Delta S$ decreases, but with $M=N$ we have $\Delta t = T/M = O(1/N)$, which is not small enough compared to $(\Delta S)^2 = O(1/N^2)$. As a result, the explicit solution becomes unstable and produces numerical overflow (huge errors).

By contrast, the implicit (Backward Euler) and Crank–Nicolson schemes are unconditionally stable and remain well-behaved as we refine the grid. Therefore, the convergence comparison should be interpreted as:

• Explicit: can diverge if the stability condition is violated.
• Implicit / Crank–Nicolson: stable and convergent under grid refinement.

Section 5 wrap-up

We implemented three finite-difference solvers for the Black–Scholes PDE:

• Explicit scheme: simple and fast per step, but conditionally stable. If the time step is not small enough, the numerical solution can diverge.
• Implicit (Backward Euler): unconditionally stable, but only first-order accurate in time.
• Crank–Nicolson: unconditionally stable and typically more accurate (second-order in time and space), at the cost of solving a tridiagonal system at each time step.

Overall, Crank–Nicolson provides the best accuracy–stability trade-off for European option pricing in this setting.

5.8 Computational Cost Analysis

We now compare the computational cost of:

• Monte Carlo
• Explicit finite differences
• Implicit finite differences
• Crank–Nicolson scheme

For each method we measure:

• Price at S0
• Absolute error vs Black–Scholes
• Execution time

import time
import pandas as pd
from mc import mc_price_plain
from pde import solve_bs_explicit, solve_bs_implicit, solve_bs_crank_nicolson
from bs import bs_price

results = []

# Reference price
price_exact = bs_price(S0, K, T, r, sigma, option="call")

# --- Monte Carlo ---
start = time.perf_counter()
price_mc, _ = mc_price_plain(S0, K, T, r, sigma, n_paths=200_000)
elapsed = time.perf_counter() - start

results.append([
    "Monte Carlo",
    price_mc,
    abs(price_mc - price_exact),
    elapsed
])

# --- Explicit FD ---
start = time.perf_counter()
Sg, tg, Vexp = solve_bs_explicit(K, T, r, sigma, N=200, M=200)
price_exp = float(np.interp(S0, Sg, Vexp[:, 0]))
elapsed = time.perf_counter() - start

results.append([
    "Explicit FD",
    price_exp,
    abs(price_exp - price_exact),
    elapsed
])

# --- Implicit FD ---
start = time.perf_counter()
Sg, tg, Vimp = solve_bs_implicit(K, T, r, sigma, N=200, M=200)
price_imp = float(np.interp(S0, Sg, Vimp[:, 0]))
elapsed = time.perf_counter() - start

results.append([
    "Implicit FD",
    price_imp,
    abs(price_imp - price_exact),
    elapsed
])

# --- Crank-Nicolson ---
start = time.perf_counter()
Sg, tg, Vcn = solve_bs_crank_nicolson(K, T, r, sigma, N=200, M=200)
price_cn = float(np.interp(S0, Sg, Vcn[:, 0]))
elapsed = time.perf_counter() - start

results.append([
    "Crank-Nicolson",
    price_cn,
    abs(price_cn - price_exact),
    elapsed
])

df_cost = pd.DataFrame(results, columns=[
    "Method",
    "Price at S0",
    "Abs Error",
    "Time (sec)"
])

df_cost

	Method	Price at S0	Abs Error	Time (sec)
0	Monte Carlo	1.041414e+01	3.644827e-02	0.032200
1	Explicit FD	-1.076981e+115	1.076981e+115	0.007625
2	Implicit FD	1.044920e+01	1.387552e-03	0.053650
3	Crank-Nicolson	1.045444e+01	3.854464e-03	0.050157

6. Generic Parabolic PDE Form

The Black–Scholes PDE belongs to the class of linear parabolic PDEs of the form: $$ ∂u/∂t = a(x,t) ∂²u/∂x² + b(x,t) ∂u/∂x + c(x,t) u $$

with appropriate boundary and terminal conditions.

The numerical schemes implemented above (Explicit, Implicit, Crank–Nicolson) can be applied to a wide class of parabolic PDEs beyond Black–Scholes.

This highlights that option pricing is a special case of a broader class of diffusion-type equations.

7. Calibration on a synthetic implied volatility surface

import numpy as np
import matplotlib.pyplot as plt


try:
    from plots import savefig
except Exception:
    def savefig(name):
        pass


S0 = 100.0
r = 0.02

# Grid for strikes and maturities
K_grid = np.linspace(60, 140, 41)          # strikes
T_grid = np.linspace(0.05, 2.0, 25)        # maturities in years

K_mesh, T_mesh = np.meshgrid(K_grid, T_grid, indexing="ij")


def sigma_synth(S0, K, T):
    x = np.log(K / S0)  # log-moneyness

    base = 0.18
    term = 0.06 * np.sqrt(T)                # term structure (higher vol for longer T)
    skew = -0.25 * x * np.exp(-0.7 * T)     # skew that fades with maturity
    smile = 0.45 * (x**2)                   # convex smile

    sig = base + term + skew + smile
    return np.clip(sig, 0.05, 1.50)

IV_true = sigma_synth(S0, K_mesh, T_mesh)

# 3D plot 
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(K_mesh, T_mesh, IV_true, linewidth=0, antialiased=True)

ax.set_title("Synthetic Implied Volatility Surface σ(K,T)")
ax.set_xlabel("Strike K")
ax.set_ylabel("Maturity T")
ax.set_zlabel("Implied vol σ")

savefig("iv_surface_3d")
plt.show()

# Heatmap 
plt.figure()
plt.imshow(
    IV_true.T,
    origin="lower",
    aspect="auto",
    extent=[K_grid.min(), K_grid.max(), T_grid.min(), T_grid.max()],
)
plt.colorbar()
plt.title("Synthetic IV Heatmap σ(K,T)")
plt.xlabel("Strike K")
plt.ylabel("Maturity T")
savefig("iv_surface_heatmap")
plt.show()

from bs import bs_price

# Build synthetic market prices using BS with the synthetic IV surface
C_mkt = np.zeros_like(IV_true)

for iK, K in enumerate(K_grid):
    for iT, T in enumerate(T_grid):
        sigma = IV_true[iK, iT]
        C_mkt[iK, iT] = bs_price(S0, K, T, r, sigma, option="call")

# 3D price surface
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(K_mesh, T_mesh, C_mkt, linewidth=0, antialiased=True)

ax.set_title("BS Price Surface C(K,T) from Synthetic IV")
ax.set_xlabel("Strike K")
ax.set_ylabel("Maturity T")
ax.set_zlabel("Call price C")

savefig("price_surface_3d")
plt.show()

# Heatmap
plt.figure()
plt.imshow(
    C_mkt.T,
    origin="lower",
    aspect="auto",
    extent=[K_grid.min(), K_grid.max(), T_grid.min(), T_grid.max()],
)
plt.colorbar()
plt.title("Call Price Heatmap C(K,T)")
plt.xlabel("Strike K")
plt.ylabel("Maturity T")
savefig("price_surface_heatmap")
plt.show()

from bs import implied_vol_newton, implied_vol_bisection, bs_price

IV_newton = np.full_like(IV_true, np.nan, dtype=float)
IV_bis    = np.full_like(IV_true, np.nan, dtype=float)


def implied_vol_safe(S0, K, T, r, price, option="call"):
    
    sigma_init = 0.2

    # Try Newton first
    try:
        sig_n = implied_vol_newton(
            S0, K, T, r, price, option=option,
            sigma_init=sigma_init, tol=1e-10, max_iter=100
        )
        # sanity clamp 
        if not np.isfinite(sig_n) or sig_n <= 0:
            raise ValueError("Newton produced invalid sigma.")
        return sig_n, "newton"
    except Exception:
        # Fallback to bisection 
        sig_b = implied_vol_bisection(
            S0, K, T, r, price, option=option,
            sigma_min=1e-6, sigma_max=5.0, tol=1e-10, max_iter=300
        )
        return sig_b, "bisection"

methods_used = {"newton": 0, "bisection": 0}

for iK, K in enumerate(K_grid):
    for iT, T in enumerate(T_grid):
        price = C_mkt[iK, iT]

        sig, method = implied_vol_safe(S0, K, T, r, price, option="call")
        methods_used[method] += 1

        # store in both arrays: Newton only when it truly worked
        if method == "newton":
            IV_newton[iK, iT] = sig
        IV_bis[iK, iT] = sig  # final "robust" IV 

# Errors
err_newton = np.abs(IV_newton - IV_true)   # will be NaN where Newton failed
err_robust = np.abs(IV_bis - IV_true)      # always defined

# Heatmap Newton error (mask NaNs)
plt.figure()
plt.imshow(
    np.nan_to_num(err_newton.T, nan=np.nan),
    origin="lower",
    aspect="auto",
    extent=[K_grid.min(), K_grid.max(), T_grid.min(), T_grid.max()],
)
plt.colorbar()
plt.title("IV Abs Error Heatmap |σ_newton - σ_true| (NaN where Newton fails)")
plt.xlabel("Strike K")
plt.ylabel("Maturity T")
savefig("iv_error_heatmap_newton_masked")
plt.show()

# Heatmap robust error
plt.figure()
plt.imshow(
    err_robust.T,
    origin="lower",
    aspect="auto",
    extent=[K_grid.min(), K_grid.max(), T_grid.min(), T_grid.max()],
)
plt.colorbar()
plt.title("IV Abs Error Heatmap |σ_recovered - σ_true| (Newton + fallback)")
plt.xlabel("Strike K")
plt.ylabel("Maturity T")
savefig("iv_error_heatmap_robust")
plt.show()

print("Points recovered with Newton   :", methods_used["newton"])
print("Points recovered with Bisection:", methods_used["bisection"])

# Summary stats
print("Robust max abs error :", float(np.max(err_robust)))
print("Robust mean abs error:", float(np.mean(err_robust)))

# Newton-only stats (ignore NaNs)
print("Newton max abs error (where it works):", float(np.nanmax(err_newton)))
print("Newton mean abs error (where it works):", float(np.nanmean(err_newton)))

Points recovered with Newton   : 978
Points recovered with Bisection: 47
Robust max abs error : 0.03868945691741307
Robust mean abs error: 0.00021972901359723775
Newton max abs error (where it works): 0.03868945691741307
Newton mean abs error (where it works): 0.00023027882544505255

The Newton method is extremely fast but locally unstable when Vega is small (deep ITM/OTM, short maturities). A hybrid Newton–Bisection approach ensures global robustness while preserving quadratic convergence where possible.

import numpy as np
import pandas as pd
import time

mask_ok = np.isfinite(IV_newton)
success_rate = 100.0 * mask_ok.mean()

# Timing (assumes you still have the loops available; otherwise skip timing)
# If you already computed IV_newton / IV_bis, we can still summarize without timing.

summary = pd.DataFrame({
    "Metric": [
        "Newton success rate (%)",
        "Newton mean abs err (ok only)",
        "Newton max abs err (ok only)",
        "Robust mean abs err",
        "Robust max abs err",
    ],
    "Value": [
        success_rate,
        float(np.nanmean(np.abs(IV_newton - IV_true))),
        float(np.nanmax(np.abs(IV_newton - IV_true))),
        float(np.mean(np.abs(IV_bis - IV_true))),
        float(np.max(np.abs(IV_bis - IV_true))),
    ]
})

summary

	Metric	Value
0	Newton success rate (%)	95.414634
1	Newton mean abs err (ok only)	0.000230
2	Newton max abs err (ok only)	0.038689
3	Robust mean abs err	0.000220
4	Robust max abs err	0.038689

General Conclusion

This project provided a unified numerical study of European option pricing under the Black–Scholes framework.

We began with the continuous-time stochastic model and derived both the risk-neutral valuation formula and the associated Black–Scholes partial differential equation. The closed-form solution served as a theoretical benchmark to validate all numerical implementations.

Three main computational approaches were developed and compared:

• Monte Carlo simulation, highlighting its flexibility and dimension-independence, but also its slow $ \mathcal{O}(N^{-1/2}) $ convergence rate.
• Explicit finite-difference schemes, illustrating the importance of stability conditions and CFL-type constraints.
• Implicit and Crank–Nicolson schemes, demonstrating unconditional stability and superior robustness under grid refinement.

The experiments confirmed several key theoretical insights:

• Monte Carlo estimators are unbiased and converge at the expected statistical rate.
• Explicit PDE schemes may diverge when stability conditions are violated.
• Implicit and Crank–Nicolson methods remain stable and converge deterministically as the grid is refined.

From a computational perspective, the study highlights the fundamental trade-offs between flexibility, stability, convergence speed, and numerical cost. Monte Carlo methods are well-suited for high-dimensional or path-dependent derivatives, while PDE-based approaches provide faster deterministic convergence for low-dimensional problems.

Overall, this project illustrates the deep connection between probabilistic valuation and partial differential equation methods in quantitative finance. Mastering both perspectives is essential for understanding modern derivative pricing, model calibration, and the numerical challenges encountered in real-world quantitative trading and risk management.