new example: Asian option pricing using 2D PDE (Pawel Magnuszewski's MSc project draft, to be submitted to arXiv soon)

docs/bibliography.json

@@ -186,5 +186,12 @@
     ],
     "label": "Smolarkiewicz and Pudykiewicz 1992 (J. Atmos. Sci. 49)",
     "title": "A Class of Semi-Lagrangian Approximations for Fluids"
+  },
+  "https://doi.org/10.1017/CBO9781139017404": {
+    "usages": [
+      "examples/PyMPDATA_examples/Magnuszewski_et_al_2025/monte_carlo.py"
+    ],
+    "label": "Capiński and Zastawniak 2012 (Cambridge University Press)",
+    "title": "Numerical Methods in Finance with C++"
   }
 }

examples/PyMPDATA_examples/Magnuszewski_et_al_2025/asian_option.py

@@ -0,0 +1,117 @@
+from functools import cached_property
+
+import numpy as np
+from PyMPDATA_examples.utils.discretisation import discretised_analytical_solution
+
+from PyMPDATA import ScalarField, Solver, Stepper, VectorField
+from PyMPDATA.boundary_conditions import Extrapolated
+from PyMPDATA.impl.enumerations import INNER, OUTER
+
+
+# pylint: disable=too-few-public-methods
+class Settings:
+    def __init__(self, T, sgma, r, K, S_min, S_max):
+        self.T = T
+        self.sgma = sgma
+        self.r = r
+        self.K = K
+        self.S_min = S_min
+        self.S_max = S_max
+        self.rh = None
+
+    def payoff(self, A: np.ndarray, da: np.float32 = 0, variant="call"):
+        def call(x):
+            return np.maximum(0, x - self.K)
+
+        def put(x):
+            return np.maximum(0, self.K - x)
+
+        if variant == "call":
+            payoff_func = call
+        else:
+            payoff_func = put
+
+        self.rh = np.linspace(A[0] - da / 2, A[-1] + da / 2, len(A) + 1)
+        output = discretised_analytical_solution(self.rh, payoff_func)
+        return output
+
+
+class Simulation:
+    def __init__(self, settings, *, nx, ny, nt, options, variant="call"):
+        self.nx = nx
+        self.nt = nt
+        self.settings = settings
+        self.ny = ny
+        self.dt = settings.T / self.nt
+        log_s_min = np.log(settings.S_min)
+        log_s_max = np.log(settings.S_max)
+        self.S = np.exp(np.linspace(log_s_min, log_s_max, self.nx))
+        self.A, self.dy = np.linspace(
+            0, settings.S_max, self.ny, retstep=True, endpoint=True
+        )
+        self.dx = (log_s_max - log_s_min) / self.nx
+        self.settings = settings
+        sigma_squared = pow(settings.sgma, 2)
+        courant_number_x = -(0.5 * sigma_squared - settings.r) * (-self.dt) / self.dx
+        self.l2 = self.dx * self.dx / sigma_squared / self.dt
+        self.mu_coeff = (0.5 / self.l2, 0)
+        assert (
+            self.l2 > 2
+        ), f"Lambda squared should be more than 2 for stability {self.l2}"
+        self.payoff = settings.payoff(A=self.A, da=self.dy, variant=variant)
+        stepper = Stepper(options=options, n_dims=2)
+        x_dim_advector = np.full(
+            (self.nx + 1, self.ny),
+            courant_number_x,
+            dtype=options.dtype,
+        )
+        cfl_condition = np.max(np.abs(self.a_dim_advector)) + np.max(
+            np.abs(x_dim_advector)
+        )
+        assert cfl_condition < 1, f"CFL condition not met {cfl_condition}"
+        self.solver = Solver(
+            stepper=stepper,
+            advectee=ScalarField(
+                self.payoff_2d.astype(dtype=options.dtype)
+                * np.exp(-self.settings.r * self.settings.T),
+                halo=options.n_halo,
+                boundary_conditions=self.boundary_conditions,
+            ),
+            advector=VectorField(
+                (x_dim_advector, self.a_dim_advector),
+                halo=options.n_halo,
+                boundary_conditions=self.boundary_conditions,
+            ),
+        )
+        self.rhs = np.zeros((self.nx, self.ny))
+
+    @property
+    def payoff_2d(self):
+        raise NotImplementedError()
+
+    @property
+    def a_dim_advector(self):
+        raise NotImplementedError()
+
+    @property
+    def boundary_conditions(self):
+        return (
+            Extrapolated(OUTER),
+            Extrapolated(INNER),
+        )
+
+    def step(self, nt=1):
+        self.solver.advance(n_steps=nt, mu_coeff=self.mu_coeff)
+
+
+class AsianArithmetic(Simulation):
+    @cached_property
+    def payoff_2d(self):
+        return np.repeat([self.payoff], self.nx, axis=0)
+
+    @property
+    def a_dim_advector(self):
+        a_dim_advector = np.zeros((self.nx, self.ny + 1))
+        for i in range(self.ny + 1):
+            a_dim_advector[:, i] = -self.dt / self.dy * self.S / self.settings.T
+        return a_dim_advector

examples/PyMPDATA_examples/Magnuszewski_et_al_2025/monte_carlo.py

@@ -0,0 +1,124 @@
+"""
+This code is a Python numba-fied implementation of the Monte Carlo method
+for pricing Asian options taken from
+[Numerical Methods in Finance with C++](https://doi.org/10.1017/CBO9781139017404)
+"""
+
+from functools import cached_property, lru_cache, partial
+from typing import Callable
+
+import numba
+import numpy as np
+
+jit = partial(numba.jit, fastmath=True, error_model="numpy", cache=True, nopython=True)
+
+# pylint: disable=too-few-public-methods
+
+
+class BSModel:
+    def __init__(self, S0, r, sigma, T, M, seed):
+        self.S0 = S0
+        self.r = r
+        self.sigma = sigma
+        self.sigma2 = sigma * sigma
+        self.b = r - 0.5 * self.sigma2
+        self.T = T
+        self.M = M
+        self.t = np.linspace(0, T, M)
+        self.bt = self.b * self.t
+        self.sqrt_tm = np.sqrt(T / M)
+        self.seed = seed
+
+    @cached_property
+    def generate_path(self):
+        M = self.M
+        S0 = self.S0
+        bt = self.bt
+        sigma = self.sigma
+        sqrt_tm = self.sqrt_tm
+        seed = self.seed
+
+        @jit
+        def numba_seed():
+            np.random.seed(seed)
+
+        if seed is not None:
+            numba_seed()
+
+        @jit
+        def body(path):
+            path[:] = S0 * np.exp(
+                bt + sigma * np.cumsum(np.random.standard_normal(M)) * sqrt_tm
+            )
+
+        return body
+
+
+class PathDependentOption:
+    def __init__(self, T, model, N):
+        self.T = T
+        self.model = model
+        self.N = N
+        self.payoff: Callable[[np.ndarray], float] = lambda path: 0.0
+
+    @cached_property
+    def price_by_mc(self):
+        T = self.T
+        model_generate_path = self.model.generate_path
+        model_r = self.model.r
+        payoff = self.payoff
+        M = self.model.M
+        N = self.N
+
+        @jit
+        def body():
+            sum_ct = 0.0
+            path = np.empty(M)
+            for _ in range(N):
+                model_generate_path(path)
+                sum_ct += payoff(path)
+            return np.exp(-model_r * T) * (sum_ct / N)
+
+        return body
+
+
+@lru_cache
+def make_payoff(K: float, option_type: str, average_type: str = "arithmetic"):
+    assert average_type in ["arithmetic", "geometric"]
+    if average_type != "arithmetic":
+        raise NotImplementedError("Only arithmetic average is supported")
+    if option_type == "call":
+
+        @jit
+        def payoff(path):
+            return max(np.mean(path) - K, 0)
+
+    elif option_type == "put":
+
+        @jit
+        def payoff(path):
+            return max(K - np.mean(path), 0)
+
+    else:
+        raise ValueError("Invalid option")
+    return payoff
+
+
+class FixedStrikeArithmeticAsianOption(PathDependentOption):
+    def __init__(self, T, K, variant, model, N):
+        super().__init__(T, model, N)
+        self.K = K
+        self.payoff = make_payoff(K, variant)
+
+
+class FixedStrikeGeometricAsianOption(PathDependentOption):
+    def __init__(self, T, K, variant, model, N):
+        super().__init__(T, model, N)
+        self.K = K
+
+        if variant == "call":
+            self.payoff = lambda path: max(np.exp(np.mean(np.log(path))) - K, 0)
+        elif variant == "put":
+            self.payoff = lambda path: max(K - np.exp(np.mean(np.log(path))), 0)
+        else:
+            raise ValueError("Invalid option type")