new formulae-only math-finance example notebook: analytic solution for valuation of fixed-strike geometric Asian options

examples/PyMPDATA_examples/asian_option/Mkhize_2007.py

@@ -0,0 +1,11 @@
+import numpy as np
+from scipy.stats import norm
+
+def geometric_mkhize(s_t=100, K=100, r=0.008, sigma=0.2, T=30, T_0=0):
+    d_1 = (np.log(s_t / K) + 0.5 * (r + (sigma ** 2) / 6) * (T - T_0)) / (sigma * np.sqrt((T - T_0) / 3))
+    d_2 = d_1 - sigma * np.sqrt((T - T_0) / 3)
+    C_0 = s_t * np.exp(-0.5 * (r + (sigma ** 2) / 6) * (T - T_0)) * norm.cdf(d_1) - K * np.exp(
+        -r * (T - T_0)) * norm.cdf(d_2)
+    P_0 = K * np.exp(-r * (T - T_0)) * norm.cdf(-d_2) - s_t * np.exp(
+        -0.5 * (r + (sigma ** 2) / 6) * (T - T_0)) * norm.cdf(-d_1)
+    return C_0

examples/PyMPDATA_examples/asian_option/__init__.py

@@ -0,0 +1 @@
+from .simulation import Simulation

examples/PyMPDATA_examples/asian_option/colors.py

@@ -0,0 +1 @@
+colors = ("purple", "teal", "turquoise")

examples/PyMPDATA_examples/asian_option/options.py

@@ -0,0 +1,8 @@
+OPTIONS = {
+    "n_iters": 2,
+    "infinite_gauge": True,
+    "nonoscillatory": True,
+    "divergent_flow": True,
+    "third_order_terms": True,
+    "non_zero_mu_coeff": True,
+}

examples/PyMPDATA_examples/asian_option/setup3_asian_option.py

@@ -0,0 +1,39 @@
+import numpy as np
+import PyMPDATA_examples.asian_option.Mkhize_2007 as MKH
+from pystrict import strict
+from scipy.stats import norm
+
+
+
+
+@strict
+class Settings:
+    # S0 = 55
+    T = 0.5
+    # amer = False
+    S_min = 10
+    S_max = 2000
+    sigma = 0.6
+    r = 0.008
+    K1 = 75
+    # K2 = 175
+    S_match = 175
+
+    def __init__(self, *, n_iters: int = 2, l2_opt: int = 2, C_opt: float = 0.034):
+        self.n_iters = n_iters
+        self.l2_opt = l2_opt
+        self.C_opt = C_opt
+
+    def payoff(self, A: np.ndarray):
+        return np.exp(-self.r * self.T) * np.maximum(0, A - self.K1)
+
+        # return np.exp(-self.r * self.T) * (
+        #     np.maximum(0, self.K2 - S) - np.maximum(0, self.K1 - S)
+        # )
+
+    def analytical_solution(self, S: np.ndarray):
+        return MKH.geometric_mkhize(s_t=S, K=self.K1, r=self.r, sigma=self.sigma, T=self.T, T_0=0)
+
+        # return BS73.p_euro(
+        #     S, K=self.K2, T=self.T, r=self.r, b=self.r, sgma=self.sigma
+        # ) - BS73.p_euro(S, K=self.K1, T=self.T, r=self.r, b=self.r, sgma=self.sigma)

examples/PyMPDATA_examples/asian_option/simulation.py

@@ -0,0 +1,119 @@
+import numba
+import numpy as np
+from PyMPDATA_examples.asian_option.options import OPTIONS
+
+from PyMPDATA import Options, ScalarField, Solver, Stepper, VectorField
+from PyMPDATA.boundary_conditions import Extrapolated, Constant, Periodic
+
+
+class Simulation:
+    @staticmethod
+    def _factory(
+        *, options: Options, advectee: np.ndarray, advector_s: float, advector_a: np.ndarray, boundary_conditions
+    ):
+        stepper = Stepper(
+            options=options, n_dims=len(advectee.shape), non_unit_g_factor=False
+        )
+        a_dim_advector = np.multiply.outer(np.ones(advectee.shape[0]+1, dtype=options.dtype), advector_a)
+        x_dim_advector = np.full((advectee.shape[0], advectee.shape[1]+1), advector_s, dtype=options.dtype)
+        print(f"{x_dim_advector.shape=}", f"{a_dim_advector.shape=}")
+        advector_values = (a_dim_advector, x_dim_advector)
+        # print(advector_values)
+        return Solver(
+            stepper=stepper,
+            advectee=ScalarField(
+                advectee.astype(dtype=options.dtype),
+                halo=options.n_halo,
+                boundary_conditions=boundary_conditions,
+            ),
+            advector=VectorField(
+                advector_values,
+                halo=options.n_halo,
+                boundary_conditions=boundary_conditions,
+            ),
+        )
+
+    def __init__(self, settings):
+        self.settings = settings
+
+        sigma2 = pow(settings.sigma, 2)
+        dx_opt = abs(
+            settings.C_opt / (0.5 * sigma2 - settings.r) * settings.l2_opt * sigma2
+        )
+        dt_opt = pow(dx_opt, 2) / sigma2 / settings.l2_opt
+
+        # adjusting dt so that nt is integer
+        self.dt = settings.T
+        self.nt = 0
+        while self.dt > dt_opt:
+            self.nt += 1
+            self.dt = settings.T / self.nt
+
+        # adjusting dx to match requested l^2
+        dx = np.sqrt(settings.l2_opt * self.dt) * settings.sigma
+
+        # calculating actual u number and lambda
+        self.C = -(0.5 * sigma2 - settings.r) * (-self.dt) / dx
+
+        self.l2 = dx * dx / sigma2 / self.dt
+
+        # adjusting nx and setting S_beg, S_end
+        S_beg = settings.S_match
+        self.nx = 1
+        while S_beg > settings.S_min:
+            self.nx += 1
+            S_beg = np.exp(np.log(settings.S_match) - self.nx * dx)
+
+        self.na = 15 # TODO: why?
+
+
+        self.ix_match = self.nx
+
+        S_end = settings.S_match
+        while S_end < settings.S_max:
+            self.nx += 1
+            S_end = np.exp(np.log(S_beg) + (self.nx - 1) * dx)
+
+        # asset price
+        self.S = np.exp(np.log(S_beg) + np.arange(self.nx) * dx)
+        self.A, self.da = np.linspace(0, S_end, self.na, retstep=True)
+        print(f"{self.S.shape=}, {self.A.shape=}")
+
+
+        # a advector
+        self.C_a = (self.S / settings.T) * (-self.dt)/self.da
+        try:
+            assert np.max(np.abs(self.C_a)) < 1
+        except AssertionError:
+            print(f"{np.max(np.abs(self.C_a))=}")
+            raise
+        # meshgrid
+        self.S_mesh, self.A_mesh = np.meshgrid(self.S, self.A)
+        print(f"{self.S_mesh.shape=}, {self.A_mesh.shape=}")
+
+        self.mu_coeff = (0.5 / self.l2, 0)
+        self.solvers = {}
+        # self.solvers[1] = self._factory(
+        #     advectee=settings.payoff(self.A_mesh),
+        #     advector=self.C,
+        #     options=Options(n_iters=1, non_zero_mu_coeff=True),
+        #     boundary_conditions=(Extrapolated(), Extrapolated()),
+        #     time_to_maturity=settings.T,
+        #     advectee_x_values=self.S,
+        # )
+        self.solvers[2] = self._factory(
+            advectee=settings.payoff(self.A_mesh),
+            advector_s=self.C,
+            advector_a=self.C_a,
+            options=Options(**OPTIONS),
+            boundary_conditions=(Periodic(), Extrapolated()),
+            # time_to_maturity=settings.T,
+            # advectee_x_values=self.S,
+        )
+
+    def run(self, n_iters: int):
+        self.solvers[n_iters].advance(self.nt, self.mu_coeff)
+        return self.solvers[n_iters].advectee.get()
+
+    # def terminal_value(self):
+    #     return self.solvers[1].advectee.get()