pde/black_scholes.viv

K = 100.0           #: K      strike
RATE = 0.05         #: RATE   risk-free interest rate, per year
SIGMA = 0.2         #: SIGMA  volatility, per root-year
T = 1.0             #: T      time to expiry, years
SMAX = 300.0        #: SMAX   where the grid stops: three strikes out
SMOOTHED = 2.0      #: SMOOTHED  how many opening steps `smooth` replaces by two half-steps each

#| The standard normal cumulative distribution.
def normal_cdf(x) {
    return 0.5 * (1.0 + erf(x / sqrt(2.0)))
}

#| The Black-Scholes closed form for a European call at stock price s.
def call_exact(s) {
    d1 = (log(s / K) + (RATE + 0.5 * SIGMA * SIGMA) * T) / (SIGMA * sqrt(T))
    d2 = d1 - SIGMA * sqrt(T)
    return s * normal_cdf(d1) - K * exp(-RATE * T) * normal_cdf(d2)
}

#| Its delta, dV/dS = N(d1).
def delta_exact(s) {
    d1 = (log(s / K) + (RATE + 0.5 * SIGMA * SIGMA) * T) / (SIGMA * sqrt(T))
    return normal_cdf(d1)
}

#| One theta-step of length dt from time-to-expiry tau: theta = 1 is implicit Euler,
#| theta = 1/2 Crank-Nicolson. `s` holds the n interior stock prices, spacing h.
def theta_step(v, s, h, tau, dt, theta) {
    n = len(v)
    diff = 0.5 * SIGMA * SIGMA * s * s / (h * h)
    drift = 0.5 * RATE * s / h
    low = diff - drift
    mid = -2.0 * diff - RATE
    high = diff + drift
    e = (1.0 - theta) * dt
    rhs = v + e * (low * roll(v, 1.0) + mid * v + high * roll(v, -1.0))
    # roll wraps around; the two edge rows see the boundary values instead
    edge_old = SMAX - K * exp(-RATE * tau)
    edge_new = SMAX - K * exp(-RATE * (tau + dt))
    last = n - 1.0
    rhs[0] = v[0] + e * (mid[0] * v[0] + high[0] * v[1])
    rhs[last] = v[last] + e * (low[last] * v[last - 1.0] + mid[last] * v[last] + high[last] * edge_old) + theta * dt * high[last] * edge_new
    return tridiag_solve(-theta * dt * low, 1.0 - theta * dt * mid, -theta * dt * high, rhs)
}

#| Price the call on n interior points in `steps` time steps. With `smooth`, the first two
#| steps are each replaced by two implicit-Euler half-steps.
def price_call(n, steps, smooth: Bool) {
    h = SMAX / (n + 1.0)
    s = linspace(h, SMAX - h, n)
    v = max(s - K, 0.0)
    dt = T / steps
    tau = 0.0
    first = 0.0
    if smooth {
        for q to 2.0 * SMOOTHED {
            v = theta_step(v, s, h, tau, 0.5 * dt, 1.0)
            tau = tau + 0.5 * dt
        }
        first = SMOOTHED
    }
    for m to steps - first {
        v = theta_step(v, s, h, tau, dt, 0.5)
        tau = tau + dt
    }
    return v
}

#| The worst error in the delta (a centered difference of the prices) within 20 of the strike.
def delta_error(n, steps, smooth: Bool) {
    h = SMAX / (n + 1.0)
    v = price_call(n, steps, smooth)
    worst = 0.0
    for i to n - 2.0 {
        s = (i + 2.0) * h
        if abs(s - K) < 20.0 {
            d = (v[i + 2.0] - v[i]) / (2.0 * h)
            worst = max(worst, abs(d - delta_exact(s)))
        }
    }
    return worst
}

#| The error in the at-the-money price. The strike sits on a grid point: (n + 1) / 3 is whole.
def price_error(n, steps, smooth: Bool) {
    v = price_call(n, steps, smooth)
    return abs(v[(n + 1.0) / 3.0 - 1.0] - call_exact(K))
}

N = 599
STEPS = 200
S = linspace(SMAX / (N + 1.0), SMAX - SMAX / (N + 1.0), N)
V = price_call(N, STEPS, true)
exact = call_exact(K)
computed = V[(N + 1.0) / 3.0 - 1.0]
show exact
show computed

# few, long time steps make dt / h^2 large: that is where the kink shows
rough = delta_error(599.0, 25.0, false)
smoothed = delta_error(599.0, 25.0, true)
show rough
show smoothed

order = log(price_error(149.0, 50.0, true) / price_error(299.0, 100.0, true)) / log(2.0)
show order

figure("A European call, one year from expiry", "stock price S", "value")
plot(S, max(S - K, 0.0), "payoff at expiry")
plot(S, V, "value today")

test "the at-the-money price matches the closed form to a tenth of a cent" {
    assert_near(computed, exact, 1.0e-3)
}
test "halving h and the time step together cuts the error by four: second order" {
    assert_near(order, 2.0, 0.1)
}
test "with long time steps, plain Crank-Nicolson leaves the kink ringing in the delta" {
    assert_gt(rough, 10.0 * smoothed)
}
test "four implicit half-steps at the start remove it" {
    assert_lt(smoothed, 1.0e-3)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors