PyMPDATA_examples.Magnuszewski_et_al_2025.monte_carlo

  1"""
  2This code is a Python numba-fied implementation of the Monte Carlo method
  3for pricing Asian options taken from
  4[Numerical Methods in Finance with C++](https://doi.org/10.1017/CBO9781139017404)
  5"""
  6
  7from functools import cached_property, lru_cache, partial
  8from typing import Callable
  9
 10import numba
 11import numpy as np
 12
 13jit = partial(numba.jit, fastmath=True, error_model="numpy", cache=True, nopython=True)
 14
 15# pylint: disable=too-few-public-methods
 16
 17
 18class BSModel:
 19    def __init__(self, S0, r, sigma, T, M, seed):
 20        self.S0 = S0
 21        self.r = r
 22        self.sigma = sigma
 23        self.sigma2 = sigma * sigma
 24        self.b = r - 0.5 * self.sigma2
 25        self.T = T
 26        self.M = M
 27        self.t = np.linspace(0, T, M)
 28        self.bt = self.b * self.t
 29        self.sqrt_tm = np.sqrt(T / M)
 30        self.seed = seed
 31
 32    @cached_property
 33    def generate_path(self):
 34        M = self.M
 35        S0 = self.S0
 36        bt = self.bt
 37        sigma = self.sigma
 38        sqrt_tm = self.sqrt_tm
 39        seed = self.seed
 40
 41        @jit
 42        def numba_seed():
 43            np.random.seed(seed)
 44
 45        if seed is not None:
 46            numba_seed()
 47
 48        @jit
 49        def body(path):
 50            path[:] = S0 * np.exp(
 51                bt + sigma * np.cumsum(np.random.standard_normal(M)) * sqrt_tm
 52            )
 53
 54        return body
 55
 56
 57class PathDependentOption:
 58    def __init__(self, T, model, N):
 59        self.T = T
 60        self.model = model
 61        self.N = N
 62        self.payoff: Callable[[np.ndarray], float] = lambda path: 0.0
 63
 64    @cached_property
 65    def price_by_mc(self):
 66        T = self.T
 67        model_generate_path = self.model.generate_path
 68        model_r = self.model.r
 69        payoff = self.payoff
 70        M = self.model.M
 71        N = self.N
 72
 73        @jit
 74        def body():
 75            sum_ct = 0.0
 76            path = np.empty(M)
 77            for _ in range(N):
 78                model_generate_path(path)
 79                sum_ct += payoff(path)
 80            return np.exp(-model_r * T) * (sum_ct / N)
 81
 82        return body
 83
 84
 85@lru_cache
 86def make_payoff(K: float, option_type: str, average_type: str = "arithmetic"):
 87    assert average_type in ["arithmetic", "geometric"]
 88    if average_type != "arithmetic":
 89        raise NotImplementedError("Only arithmetic average is supported")
 90    if option_type == "call":
 91
 92        @jit
 93        def payoff(path):
 94            return max(np.mean(path) - K, 0)
 95
 96    elif option_type == "put":
 97
 98        @jit
 99        def payoff(path):
100            return max(K - np.mean(path), 0)
101
102    else:
103        raise ValueError("Invalid option")
104    return payoff
105
106
107class FixedStrikeArithmeticAsianOption(PathDependentOption):
108    def __init__(self, T, K, variant, model, N):
109        super().__init__(T, model, N)
110        self.K = K
111        self.payoff = make_payoff(K, variant)
112
113
114class FixedStrikeGeometricAsianOption(PathDependentOption):
115    def __init__(self, T, K, variant, model, N):
116        super().__init__(T, model, N)
117        self.K = K
118
119        if variant == "call":
120            self.payoff = lambda path: max(np.exp(np.mean(np.log(path))) - K, 0)
121        elif variant == "put":
122            self.payoff = lambda path: max(K - np.exp(np.mean(np.log(path))), 0)
123        else:
124            raise ValueError("Invalid option type")