Skip to content

Vectorized BSM IV Solver#

High-performance, vectorized Newton-Raphson solver for BSM implied volatility.

This solver is designed to calculate the implied volatility for a large number of options simultaneously, leveraging NumPy for vectorized operations.

Source code in src/quantfin/calibration/vectorized_bsm_iv.py 
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
class BSMIVSolver:
    """
    High-performance, vectorized Newton-Raphson solver for BSM implied volatility.

    This solver is designed to calculate the implied volatility for a large
    number of options simultaneously, leveraging NumPy for vectorized operations.
    """

    def __init__(
        self,
        max_iter: int = 20,
        tolerance: float = 1e-6,
    ):
        """
        Initializes the BSM implied volatility solver.

        Parameters
        ----------
        max_iter : int, optional
            The maximum number of iterations for the Newton-Raphson method,
            by default 20.
        tolerance : float, optional
            The error tolerance to determine convergence, by default 1e-6.
        """
        self.max_iter = max_iter
        self.tolerance = tolerance

    def solve(
        self,
        target_prices: np.ndarray,
        options: pd.DataFrame,
        stock: Stock,
        rate: Rate,
    ) -> np.ndarray:
        """
        Calculates implied volatility for an array of options.

        Parameters
        ----------
        target_prices : np.ndarray
            An array of market prices for which to find the implied volatility.
        options : pd.DataFrame
            A DataFrame of option contracts, must include 'strike', 'maturity',
            and 'optionType' columns.
        stock : Stock
            The underlying asset's properties.
        rate : Rate
            The risk-free rate structure.

        Returns
        -------
        np.ndarray
            An array of calculated implied volatilities corresponding to the
            target prices.
        """
        S, q = stock.spot, stock.dividend
        K, T = options["strike"].values, options["maturity"].values
        r = rate.get_rate(T)  # Use get_rate for term structure
        is_call = options["optionType"].values == "call"
        iv = np.full_like(target_prices, 0.20)

        for _ in range(self.max_iter):
            with np.errstate(all="ignore"):
                sqrt_T = np.sqrt(T)
                d1 = (np.log(S / K) + (r - q + 0.5 * iv**2) * T) / (iv * sqrt_T)
                d2 = d1 - iv * sqrt_T
                vega = S * np.exp(-q * T) * sqrt_T * norm.pdf(d1)
                call_prices = S * np.exp(-q * T) * norm.cdf(d1) - K * np.exp(
                    -r * T
                ) * norm.cdf(d2)
                put_prices = K * np.exp(-r * T) * norm.cdf(-d2) - S * np.exp(
                    -q * T
                ) * norm.cdf(-d1)
            model_prices = np.where(is_call, call_prices, put_prices)
            error = model_prices - target_prices
            if np.all(np.abs(error) < self.tolerance):
                break
            iv = iv - error / np.maximum(vega, 1e-8)
        return iv

__init__(max_iter: int = 20, tolerance: float = 1e-06) #

Initializes the BSM implied volatility solver.

Parameters:

Name Type Description Default
max_iter int

The maximum number of iterations for the Newton-Raphson method, by default 20.

 20
tolerance float

The error tolerance to determine convergence, by default 1e-6.

 1e-06
Source code in src/quantfin/calibration/vectorized_bsm_iv.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def __init__(
    self,
    max_iter: int = 20,
    tolerance: float = 1e-6,
):
    """
    Initializes the BSM implied volatility solver.

    Parameters
    ----------
    max_iter : int, optional
        The maximum number of iterations for the Newton-Raphson method,
        by default 20.
    tolerance : float, optional
        The error tolerance to determine convergence, by default 1e-6.
    """
    self.max_iter = max_iter
    self.tolerance = tolerance

solve(target_prices: np.ndarray, options: pd.DataFrame, stock: Stock, rate: Rate) -> np.ndarray #

Calculates implied volatility for an array of options.

Parameters:

Name Type Description Default
target_prices ndarray

An array of market prices for which to find the implied volatility.

 required
options DataFrame

A DataFrame of option contracts, must include 'strike', 'maturity', and 'optionType' columns.

 required
stock Stock

The underlying asset's properties.

 required
rate Rate

The risk-free rate structure.

 required

Returns:

Type Description
ndarray

An array of calculated implied volatilities corresponding to the target prices.
Source code in src/quantfin/calibration/vectorized_bsm_iv.py 
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
def solve(
    self,
    target_prices: np.ndarray,
    options: pd.DataFrame,
    stock: Stock,
    rate: Rate,
) -> np.ndarray:
    """
    Calculates implied volatility for an array of options.

    Parameters
    ----------
    target_prices : np.ndarray
        An array of market prices for which to find the implied volatility.
    options : pd.DataFrame
        A DataFrame of option contracts, must include 'strike', 'maturity',
        and 'optionType' columns.
    stock : Stock
        The underlying asset's properties.
    rate : Rate
        The risk-free rate structure.

    Returns
    -------
    np.ndarray
        An array of calculated implied volatilities corresponding to the
        target prices.
    """
    S, q = stock.spot, stock.dividend
    K, T = options["strike"].values, options["maturity"].values
    r = rate.get_rate(T)  # Use get_rate for term structure
    is_call = options["optionType"].values == "call"
    iv = np.full_like(target_prices, 0.20)

    for _ in range(self.max_iter):
        with np.errstate(all="ignore"):
            sqrt_T = np.sqrt(T)
            d1 = (np.log(S / K) + (r - q + 0.5 * iv**2) * T) / (iv * sqrt_T)
            d2 = d1 - iv * sqrt_T
            vega = S * np.exp(-q * T) * sqrt_T * norm.pdf(d1)
            call_prices = S * np.exp(-q * T) * norm.cdf(d1) - K * np.exp(
                -r * T
            ) * norm.cdf(d2)
            put_prices = K * np.exp(-r * T) * norm.cdf(-d2) - S * np.exp(
                -q * T
            ) * norm.cdf(-d1)
        model_prices = np.where(is_call, call_prices, put_prices)
        error = model_prices - target_prices
        if np.all(np.abs(error) < self.tolerance):
            break
        iv = iv - error / np.maximum(vega, 1e-8)
    return iv