implied-vol-surface-calibration/dupire.py at main · Alex-Rib/implied-vol-surface-calibration · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
"""Dupire local volatility derived from the calibrated SVI surface

Local variance is computed from Gatheral's formula in total-variance and
log-moneyness form. Because the surface is SVI (analytic in k), the k-derivatives
are closed-form and interpolated linearly in total-variance space. The tau-derivative
uses central finite differences. The result is a pre-computed 2D grid
sigma_loc[i, j] over (k_i, tau_j), aligned with the PDE mesh
"""

import numpy as np

from vol_surface import SVIParams, VolSurface, svi_total_variance


def svi_dk(k: np.ndarray, params: SVIParams) -> np.ndarray:
    """First derivative dw/dk of the SVI total variance"""
    y = k - params.m
    return params.b * (params.rho + y / np.sqrt(y**2 + params.sigma**2))


def svi_dkk(k: np.ndarray, params: SVIParams) -> np.ndarray:
    """Second derivative d2w/dk2 of the SVI total variance"""
    y = k - params.m
    return params.b * params.sigma**2 / (y**2 + params.sigma**2) ** 1.5


def _get_interpolated_components(
    surface: VolSurface, k_grid: np.ndarray, tau: float
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Interpolated w, dw/dk and d2w/dk2 at maturity tau, vectorised over k_grid"""
    taus = surface.taus

    # left extrapolation: implied vol -> w proportional to tau
    if tau <= taus[0]:
        ratio = tau / taus[0]
        p = surface.slices[0]
        w = svi_total_variance(k_grid, p) * ratio
        wk = svi_dk(k_grid, p) * ratio
        wkk = svi_dkk(k_grid, p) * ratio
        return w, wk, wkk

    # right extrapolation: implied vol
    if tau >= taus[-1]:
        ratio = tau / taus[-1]
        p = surface.slices[-1]
        w = svi_total_variance(k_grid, p) * ratio
        wk = svi_dk(k_grid, p) * ratio
        wkk = svi_dkk(k_grid, p) * ratio
        return w, wk, wkk

    # exact linear interpolation between calibrated slices
    idx = np.searchsorted(taus, tau)
    tau_1, tau_2 = taus[idx - 1], taus[idx]
    p_1, p_2 = surface.slices[idx - 1], surface.slices[idx]

    alpha = (tau_2 - tau) / (tau_2 - tau_1)
    beta = 1.0 - alpha

    w = alpha * svi_total_variance(k_grid, p_1) + beta * svi_total_variance(k_grid, p_2)
    wk = alpha * svi_dk(k_grid, p_1) + beta * svi_dk(k_grid, p_2)
    wkk = alpha * svi_dkk(k_grid, p_1) + beta * svi_dkk(k_grid, p_2)

    return w, wk, wkk


def local_vol_grid(
    surface: VolSurface,
    k_grid: np.ndarray,
    tau_grid: np.ndarray,
    floor: float = 1e-4,
    dtau: float = 1e-4,
) -> np.ndarray:
    """Pre-compute the Dupire local-volatility grid sigma_loc[i, j]"""
    nk, nt = len(k_grid), len(tau_grid)
    sigma_loc = np.empty((nk, nt))

    for j, tau in enumerate(tau_grid):
        # numerator: tau-derivative by central finite difference
        tau_up = tau + dtau
        tau_dn = max(tau - dtau, 1e-8)

        w_up, _, _ = _get_interpolated_components(surface, k_grid, tau_up)
        w_dn, _, _ = _get_interpolated_components(surface, k_grid, tau_dn)
        wt = (w_up - w_dn) / (tau_up - tau_dn)

        # denominator: guard against the w=0 singularity with a tau floor
        safe_tau = max(tau, 1e-4)
        w, wk, wkk = _get_interpolated_components(surface, k_grid, safe_tau)

        # absolute safeguard against division by zero
        w = np.maximum(w, 1e-9)

        denom = (
            1.0
            - (k_grid / w) * wk
            + 0.25 * (-0.25 - 1.0 / w + (k_grid**2) / w**2) * wk**2
            + 0.5 * wkk
        )

        local_var = wt / denom
        sigma_loc[:, j] = np.sqrt(np.maximum(local_var, floor))

    return sigma_loc