Merge pull request #38 from chemardes/feature/adding-rbf-interpolation

FEATURE: adding rbf interpolation

Author: chemardes

pdesolvers/__init__.py

@@ -1,3 +1,4 @@
 from .pdes import *
 from .solution import *
-from .solvers import *
+from .solvers import *
+from .utils import *

pdesolvers/main.py

@@ -19,11 +19,18 @@ def main():
     equation2 = pde.BlackScholesEquation('call', 300, 1, 0.2, 0.05, 100, 100, 20000)
 
     solver1 = pde.BlackScholesCNSolver(equation2)
-    res1 = solver1.solve()
-    # res2 = solver2.solve()
+    solver2 = pde.BlackScholesExplicitSolver(equation2)
+    res1 = solver1.solve().get_result()
+    res2 = solver2.solve().get_result()
 
-    print(res1.get_result())
-    res1.plot()
+    interpolator1 = pde.RBFInterpolator(res1, 0.8, 200,0.1, 0.03)
+    interpolator2 = pde.RBFInterpolator(res2, 0.8, 200,0.1, 0.03)
+    print(interpolator1.rbf_interpolate())
+    print(interpolator2.rbf_interpolate())
+
+
+# print(res.shape)
+#     res1.plot()
 
 if __name__ == "__main__":
     main()

pdesolvers/solution/solution.py

@@ -10,6 +10,11 @@ def __init__(self, result, x_grid, t_grid):
         self.t_grid = t_grid
 
     def plot(self):
+        """
+        Generates a 3D surface plot of the temperature distribution across a grid of space and time
+
+        :return: 3D surface plot
+        """
 
         if self.result is None:
             raise RuntimeError("Solution has not been computed - please run the solver.")
@@ -31,6 +36,11 @@ def plot(self):
         plt.show()
 
     def get_result(self):
+        """
+        Gets the grid of computed temperature values
+
+        :return: grid result
+        """
         return self.result
 
 class SolutionBlackScholes:
@@ -40,6 +50,11 @@ def __init__(self, result, s_grid, t_grid):
         self.t_grid = t_grid
 
     def plot(self):
+        """
+        Generates a 3D surface plot of the option values across a grid of asset prices and time
+
+        :return: 3D surface plot
+        """
 
         X, Y = np.meshgrid(self.t_grid, self.s_grid)
 
@@ -57,4 +72,9 @@ def plot(self):
         plt.show()
 
     def get_result(self):
+        """
+        Gets the grid of computed option prices
+
+        :return: grid result
+        """
         return self.result

pdesolvers/solvers/black_scholes_solvers.py

@@ -12,9 +12,13 @@ def __init__(self, equation: bse.BlackScholesEquation):
     def solve(self):
         """
         This method solves the Black-Scholes equation using the explicit finite difference method
+
         :return: the solver instance with the computed option values
         """
 
+        S = self.equation.generate_asset_grid()
+        T = self.equation.generate_time_grid()
+
         dt_max = 1/((self.equation.s_nodes**2) * (self.equation.sigma**2)) # cfl condition to ensure stability
 
         if self.equation.t_nodes is None:
@@ -23,15 +27,11 @@ def solve(self):
             dt = self.equation.expiry / self.equation.t_nodes # to ensure that the expiration time is integer time steps away
         else:
             # possible fix - set a check to see that user-defined value is within cfl condition
-            dt = self.equation.expiry / self.equation.t_nodes
+            dt = T[1] - T[0]
 
             if dt > dt_max:
                 raise ValueError("User-defined t nodes is too small and exceeds the CFL condition. Possible action: Increase number of t nodes for stability!")
 
-
-        S = self.equation.generate_asset_grid()
-        T = self.equation.generate_time_grid()
-
         dS = S[1] - S[0]
 
         V = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
@@ -47,7 +47,7 @@ def solve(self):
         for tau in reversed(range(self.equation.t_nodes)):
             for i in range(1, self.equation.s_nodes):
                 delta = (V[i+1, tau+1] - V[i-1, tau+1]) /  (2 * dS)
-                gamma = (V[i+1, tau+1] - 2 * V[i,tau+1] + V[i-1, tau+1]) / dS ** 2
+                gamma = (V[i+1, tau+1] - 2 * V[i,tau+1] + V[i-1, tau+1]) / (dS ** 2)
                 theta = -0.5 * ( self.equation.sigma ** 2) * (S[i] ** 2) * gamma - self.equation.rate * S[i] * delta + self.equation.rate * V[i, tau+1]
                 V[i, tau] = V[i, tau + 1] - (theta * dt)
 
@@ -59,6 +59,14 @@ def solve(self):
         return sol.SolutionBlackScholes(V,S,T)
 
     def __set_boundary_conditions(self, T, tau):
+        """
+        Sets the boundary conditions for the Black-Scholes Equation based on option type
+
+        :param T: grid of time steps
+        :param tau: index of current time step
+        :return: a tuple representing the boundary values for the given time step
+        """
+
         lower_boundary = None
         upper_boundary = None
         if self.equation.option_type == 'call':
@@ -76,6 +84,11 @@ def __init__(self, equation: bse.BlackScholesEquation):
         self.equation = equation
 
     def solve(self):
+        """
+        This method solves the Black-Scholes equation using the Crank-Nicolson method
+
+        :return: the solver instance with the computed option values
+        """
 
         S = self.equation.generate_asset_grid()
         T = self.equation.generate_time_grid()

pdesolvers/solvers/heat_solvers.py

@@ -10,6 +10,11 @@ def __init__(self, equation: heat.HeatEquation):
         self.equation = equation
 
     def solve(self):
+        """
+        This method solves the heat (diffusion) equation using the explicit finite difference method
+
+        :return: the solver instance with the computed temperature values
+        """
 
         x = self.equation.generate_x_grid()
         dx = x[1] - x[0]
@@ -38,6 +43,11 @@ def __init__(self, equation: heat.HeatEquation):
         self.equation = equation
 
     def solve(self):
+        """
+        This method solves the heat (diffusion) equation using the Crank Nicolson method
+
+        :return: the solver instance with the computed temperature values
+        """
 
         x = self.equation.generate_x_grid()
         t = self.equation.generate_t_grid()

pdesolvers/tests/test_black_scholes.py

@@ -3,11 +3,12 @@
 
 import pdesolvers.pdes.black_scholes as bse
 import pdesolvers.solvers.black_scholes_solvers as solver
+import pdesolvers.utils.utility as utility
 
 class TestBlackScholesSolvers:
 
     def setup_method(self):
-        self.equation = bse.BlackScholesEquation('call', 300, 1, 0.2, 0.05, 100, 500, 20000)
+        self.equation = bse.BlackScholesEquation('call', 300, 1, 0.2, 0.05, 100, 100, 2000)
 
     # explicit method tests
 
@@ -72,10 +73,21 @@ def test_check_absolute_difference_between_two_results(self):
         u2 = result2.get_result()
         diff = u1 - u2
 
-        # X, Y = np.meshgrid(result1.t_grid, result1.s_grid)
+        assert np.max(np.abs(diff)) < 1e-2
+
+    def test_convergence_between_interpolated_data(self):
+        result1 = solver.BlackScholesExplicitSolver(self.equation).solve()
+        result2 = solver.BlackScholesCNSolver(self.equation).solve()
+        u1 = result1.get_result()
+        u2 = result2.get_result()
+
+        data1 = utility.RBFInterpolator(u1, 0.1, 0.03).interpolate(4,20)
+        data2 = utility.RBFInterpolator(u2,0.1, 0.03).interpolate(4,20)
+
+        diff = np.abs(data1 - data2)
+
+        assert diff < 1e-4
+
+
+
 
-        # fig = plt.figure(figsize=(10,6))
-        # ax = fig.add_subplot(111, projection='3d')
-        # surf = ax.plot_surface(X, Y, diff, cmap='viridis')
-        print(np.max(np.abs(diff)))
-        # plt.show()

pdesolvers/utils/__init__.py

@@ -0,0 +1 @@
+from .utility import RBFInterpolator

pdesolvers/utils/utility.py

@@ -0,0 +1,105 @@
+import numpy as np
+
+
+class RBFInterpolator:
+
+    def __init__(self, z, hx, hy):
+        """
+        Initializes the RBF Interpolator.
+
+        :param z: 2D array of values at the grid points.
+        :param x: x-coordinate of the point to interpolate.
+        :param y: y-coordinate of the point to interpolate.
+        :param hx: Grid spacing in the x-direction.
+        :param hy: Grid spacing in the y-direction.
+        """
+
+        self.__z = z
+        self.__hx = hx
+        self.__hy = hy
+        self.__nx, self._ny = z.shape
+
+    def __get_coordinates(self, x, y):
+        """
+        Determines the x and y coordinates of the bottom-left corner of the grid cell
+
+        :return: A tuple containing the coordinates and its corresponding indices
+        """
+
+        # gets the grid steps to x
+        i_minus_star = int(np.floor(x / self.__hx))
+        if i_minus_star > self.__nx - 1:
+            raise Exception("x is out of bounds")
+
+        # final i index for interpolation
+        i_minus = i_minus_star if i_minus_star < self.__nx - 1 else self.__nx - 1
+
+        # gets the grid steps to y
+        j_minus_star = int(np.floor(y / self.__hy))
+        if j_minus_star > self._ny - 1:
+            raise Exception("y is out of bounds")
+
+        # final j index for interpolation
+        j_minus = j_minus_star if j_minus_star < self._ny - 1 else self._ny - 1
+
+        # computes the coordinates at the computed indices
+        x_minus = i_minus * self.__hx
+        y_minus = j_minus * self.__hy
+
+        return x_minus, y_minus, i_minus, j_minus
+
+    def __euclidean_distances(self, x_minus, y_minus, x, y):
+        """
+        Calculates Euclidean distances between (x,y) and the surrounding grid points in the unit cell
+
+        :param x_minus: x-coordinate of the bottom-left corner of the grid
+        :param y_minus: y-coordinate of the bottom-left corner of the grid
+        :return: returns tuple with the Euclidean distances to the surrounding grid points:
+                [bottom left, top left, bottom right, top right]
+        """
+
+        bottom_left = np.sqrt((x_minus - x) ** 2 + (y_minus - y) ** 2)
+        top_left = np.sqrt((x_minus - x) ** 2 + (y_minus + self.__hy - y) ** 2)
+        bottom_right = np.sqrt((x_minus + self.__hx - x) ** 2 + (y_minus - y) ** 2)
+        top_right = np.sqrt((x_minus + self.__hx - x) ** 2 + (y_minus + self.__hy - y) ** 2)
+
+        return bottom_left, top_left, bottom_right, top_right
+
+    @staticmethod
+    def __rbf(d, gamma):
+        """
+        Computes the Radial Basis Function (RBF) for a given distance and gamma
+
+        :param d: the Euclidean distance to a grid point
+        :param gamma: gamma parameter
+        :return: the RBF value for the distance d
+        """
+        return np.exp(-gamma * d ** 2)
+
+    def interpolate(self, x, y):
+        """
+        Performs the Radial Basis function (RBF) interpolation for the point (x,y)
+
+        :return: the interpolated value at (x,y)
+        """
+
+        x_minus, y_minus, i_minus, j_minus = self.__get_coordinates(x, y)
+
+        distances = self.__euclidean_distances(x_minus, y_minus, x, y)
+
+        h_diag_squared = self.__hx ** 2 + self.__hy ** 2
+        gamma = -np.log(0.005) / h_diag_squared
+
+        rbf_weights = [self.__rbf(d, gamma) for d in distances]
+
+        sum_rbf = np.sum(rbf_weights)
+        interpolated = rbf_weights[0] * self.__z[i_minus, j_minus]
+        interpolated += rbf_weights[1] * self.__z[i_minus, j_minus + 1]
+        interpolated += rbf_weights[2] * self.__z[i_minus + 1, j_minus]
+        interpolated += rbf_weights[3] * self.__z[i_minus + 1, j_minus + 1]
+        interpolated /= sum_rbf
+
+        return interpolated
+
+    def interpolate_all(self):
+        pass