Merge pull request #44 from chemardes/feature/adding-greek-plots

FEATURE: adding greek plots - delta / gamma / theta

Author: chemardes

.github/workflows/ci.yml

@@ -14,6 +14,6 @@ jobs:
         with:
           python-version: '3.13'
           architecture: 'x64'
-      - name: Run Script
+      - name: Run python tests
         run: |
-          bash ./ci/script.sh
+          bash ./ci/run_python_tests.sh

.gitignore

@@ -16,3 +16,4 @@
 
 # Other files
 *.log
+/.coverage

ci/script.sh ci/run_python_tests.sh

@@ -11,4 +11,4 @@ source venv/bin/activate
 
 pip install -r requirements.txt
 
-pytest
+pytest --cov=pdesolvers pdesolvers/tests/

pdesolvers/main.py

@@ -19,18 +19,12 @@ def main():
     equation2 = pde.BlackScholesEquation('call', 300, 1, 0.2, 0.05, 100, 100, 20000)
 
     solver1 = pde.BlackScholesCNSolver(equation2)
-    solver2 = pde.BlackScholesExplicitSolver(equation2)
-    res1 = solver1.solve().get_result()
-    res2 = solver2.solve().get_result()
+    # solver2 = pde.BlackScholesExplicitSolver(equation2)
+    sol1 = solver1.solve()
 
-    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())
+    sol1.plot_greek('gamma')
+    sol1.plot()
 
 
-# print(res.shape)
-#     res1.plot()
-
 if __name__ == "__main__":
     main()

pdesolvers/pdes/black_scholes.py

@@ -26,11 +26,9 @@ def __init__(self, option_type, S_max, expiry, sigma, r, K, s_nodes=1, t_nodes=N
         self.__t_nodes = t_nodes
         self.__V = None
 
-    def generate_asset_grid(self):
-        return np.linspace(0, self.__S_max, self.__s_nodes+1)
-
-    def generate_time_grid(self):
-        return np.linspace(0, self.__expiry, self.__t_nodes+1)
+    @staticmethod
+    def generate_grid(value, nodes):
+        return np.linspace(0, value, nodes + 1)
 
     @property
     def s_nodes(self):

pdesolvers/pdes/heat_1d.py

@@ -38,22 +38,20 @@ def __check_conditions(self):
 
         if self.__left_boundary_temp is not None:
             err = np.abs(self.__left_boundary_temp(0) - self.__initial_temp(0))
-            assert err < 1e-12
+            assert err < 1e-12, f"Left boundary condition at t=0 does not match the initial condition."
 
         if self.__right_boundary_temp is not None:
-            err = np.abs(self.__right_boundary_temp(0) - self.__initial_temp(0))
-            assert err < 1e-12
+            err = np.abs(self.__right_boundary_temp(0) - self.__initial_temp(self.__length))
+            assert err < 1e-12, f"Right boundary condition at t=0 does not match the initial condition."
 
     @staticmethod
     def __validate_callable(func):
         if not callable(func):
             raise ValueError("Temperature conditions must be a callable function")
 
-    def generate_x_grid(self):
-        return np.linspace(0, self.__length, self.__x_nodes)
-
-    def generate_t_grid(self):
-        return np.linspace(0, self.__time, self.__t_nodes)
+    @staticmethod
+    def generate_grid(value, nodes):
+        return np.linspace(0, value, nodes)
 
     @property
     def length(self):

pdesolvers/solution/solution.py

@@ -44,10 +44,13 @@ def get_result(self):
         return self.result
 
 class SolutionBlackScholes:
-    def __init__(self, result, s_grid, t_grid):
+    def __init__(self, result, s_grid, t_grid, delta, gamma, theta):
         self.result = result
         self.s_grid = s_grid
         self.t_grid = t_grid
+        self.delta = delta
+        self.gamma = gamma
+        self.theta = theta
 
     def plot(self):
         """
@@ -77,4 +80,29 @@ def get_result(self):
 
         :return: grid result
         """
-        return self.result
+        return self.result
+
+    def plot_greek(self, greek_type='delta', time_step=0):
+
+        greek_types = {
+            'delta': {'data': self.delta, 'title': 'Delta'},
+            'gamma': {'data': self.gamma, 'title': 'Gamma'},
+            'theta': {'data': self.theta, 'title': 'Theta'}
+        }
+
+        if greek_type.lower() not in greek_types:
+            raise ValueError("Invalid greek type - please choose between delta/gamma/theta.")
+
+        chosen_greek = greek_types[greek_type.lower()]
+        greek_data = chosen_greek['data'][:, time_step]
+        plt.figure(figsize=(8, 6))
+        plt.plot(self.s_grid, greek_data, label=f"Delta at t={self.t_grid[time_step]:.4f}", color="blue")
+
+        plt.title(f"{chosen_greek['title']} vs. Stock Price at t={self.t_grid[time_step]:.4f}")
+        plt.xlabel("Stock Price (S)")
+        plt.ylabel(chosen_greek['title'])
+        plt.grid()
+        plt.legend()
+
+        plt.show()
+

pdesolvers/solvers/black_scholes_solvers.py

@@ -16,23 +16,22 @@ def solve(self):
         :return: the solver instance with the computed option values
         """
 
-        S = self.equation.generate_asset_grid()
-        T = self.equation.generate_time_grid()
+        S = self.equation.generate_grid(self.equation.S_max, self.equation.s_nodes)
+        T = self.equation.generate_grid(self.equation.expiry, self.equation.t_nodes)
 
         dt_max = 1/((self.equation.s_nodes**2) * (self.equation.sigma**2)) # cfl condition to ensure stability
 
         if self.equation.t_nodes is None:
             dt = 0.9 * dt_max
             self.equation.t_nodes = int(self.equation.expiry/dt)
-            dt = self.equation.expiry / self.equation.t_nodes # to ensure that the expiration time is integer time steps away
+            dt = self.equation.expiry / self.equation.t_nodes
         else:
-            # possible fix - set a check to see that user-defined value is within cfl condition
             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!")
 
-        dS = S[1] - S[0]
+        ds = S[1] - S[0]
 
         V = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
 
@@ -44,19 +43,25 @@ def solve(self):
         else:
             raise ValueError("Invalid option type - please choose between call/put")
 
+        delta = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
+        gamma = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
+        theta = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
+
         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)
-                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)
+                delta[i, tau] = (V[i+1, tau+1] - V[i-1, tau+1]) / (2 * ds)
+                gamma[i, tau] = (V[i+1, tau+1] - 2 * V[i,tau+1] + V[i-1, tau+1]) / (ds ** 2)
+                theta[i, tau] = -0.5 * (self.equation.sigma ** 2) * (S[i] ** 2) * gamma[i, tau] - self.equation.rate * S[i] * delta[i, tau] + self.equation.rate * V[i, tau+1]
+                V[i, tau] = V[i, tau + 1] - (theta[i, tau] * dt)
 
             # setting boundary conditions
             lower, upper = self.__set_boundary_conditions(T, tau)
             V[0, tau] = lower
             V[self.equation.s_nodes, tau] = upper
 
-        return sol.SolutionBlackScholes(V,S,T)
+            delta, gamma, theta = self.__calculate_greeks_at_boundary(delta, gamma, theta, tau, V, S, ds)
+
+        return sol.SolutionBlackScholes(V,S,T, delta, gamma, theta)
 
     def __set_boundary_conditions(self, T, tau):
         """
@@ -78,6 +83,20 @@ def __set_boundary_conditions(self, T, tau):
 
         return lower_boundary, upper_boundary
 
+    def __calculate_greeks_at_boundary(self, delta, gamma, theta, tau, V, S, ds):
+        delta[0, tau] = (V[1, tau+1] - V[0, tau+1]) / ds  # Forward difference for lower boundary
+        delta[self.equation.s_nodes, tau] = (V[self.equation.s_nodes, tau+1] - V[self.equation.s_nodes-1, tau+1]) / ds  # Backward difference for upper boundary
+
+        gamma[0, tau] = (V[2, tau+1] - 2*V[1, tau+1] + V[0, tau+1]) / (ds**2)  # Forward approximation
+        gamma[self.equation.s_nodes, tau] = (V[self.equation.s_nodes, tau+1] - 2*V[self.equation.s_nodes-1, tau+1] + V[self.equation.s_nodes-2, tau+1]) / (ds**2)  # Backward approximation
+
+        # Calculate theta for boundary points using the same formula
+        theta[0, tau] = -0.5 * (self.equation.sigma**2) * (S[0]**2) * gamma[0, tau] - self.equation.rate * S[0] * delta[0, tau] + self.equation.rate * V[0, tau+1]
+        theta[self.equation.s_nodes, tau] = -0.5 * (self.equation.sigma**2) * (S[-1]**2) * gamma[self.equation.s_nodes, tau] - self.equation.rate * S[-1] * delta[self.equation.s_nodes, tau] + self.equation.rate * V[self.equation.s_nodes, tau+1]
+
+        return delta, gamma, theta
+
+
 class BlackScholesCNSolver:
 
     def __init__(self, equation: bse.BlackScholesEquation):
@@ -90,22 +109,27 @@ def solve(self):
         :return: the solver instance with the computed option values
         """
 
-        S = self.equation.generate_asset_grid()
-        T = self.equation.generate_time_grid()
+        S = self.equation.generate_grid(self.equation.S_max, self.equation.s_nodes)
+        T = self.equation.generate_grid(self.equation.expiry, self.equation.t_nodes)
 
-        dS = S[1] - S[0]
-        dT = T[1] - T[0]
+        ds = S[1] - S[0]
+        dt = T[1] - T[0]
 
-        alpha = 0.25 * dT * ((self.equation.sigma**2) * (S**2) / (dS**2) - self.equation.rate * S / dS)
-        beta = -dT * 0.5 * (self.equation.sigma**2 * (S**2) / (dS**2) + self.equation.rate)
-        gamma = 0.25 * dT * (self.equation.sigma**2 * (S**2) / (dS**2) + self.equation.rate * S / dS)
+        a = 0.25 * dt * ((self.equation.sigma**2) * (S**2) / (ds**2) - self.equation.rate * S / ds)
+        b = -dt * 0.5 * (self.equation.sigma**2 * (S**2) / (ds**2) + self.equation.rate)
+        c = 0.25 * dt * (self.equation.sigma**2 * (S**2) / (ds**2) + self.equation.rate * S / ds)
 
-        lhs = sparse.diags([-alpha[2:], 1-beta[1:], -gamma[1:-1]], [-1, 0, 1], shape = (self.equation.s_nodes - 1, self.equation.s_nodes - 1), format='csr')
-        rhs = sparse.diags([alpha[2:], 1+beta[1:], gamma[1:-1]], [-1, 0, 1], shape = (self.equation.s_nodes - 1, self.equation.s_nodes - 1) , format='csr')
+        lhs = sparse.diags([-a[2:], 1-b[1:], -c[1:-1]], [-1, 0, 1], shape = (self.equation.s_nodes - 1, self.equation.s_nodes - 1), format='csr')
+        rhs = sparse.diags([a[2:], 1+b[1:], c[1:-1]], [-1, 0, 1], shape = (self.equation.s_nodes - 1, self.equation.s_nodes - 1) , format='csr')
 
         V = np.zeros((self.equation.s_nodes+1, self.equation.t_nodes+1))
 
-        # setting terminal condition (for all values of S at time T)
+        delta = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
+        gamma = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
+        theta = np.zeros((self.equation.s_nodes + 1, self.equation.t_nodes + 1))
+
+
+# setting terminal condition (for all values of S at time T)
         if self.equation.option_type == 'call':
             V[:,-1] = np.maximum((S - self.equation.strike_price), 0)
 
@@ -123,10 +147,29 @@ def solve(self):
             rhs_vector = rhs @ V[1:-1, tau + 1]
 
             # Apply boundary conditions to the RHS vector
-            rhs_vector[0] += alpha[1] * (V[0, tau + 1] + V[0, tau])
-            rhs_vector[-1] += gamma[self.equation.s_nodes-1] *(V[-1, tau+1] + V[-1, tau])
+            rhs_vector[0] += a[1] * (V[0, tau + 1] + V[0, tau])
+            rhs_vector[-1] += c[self.equation.s_nodes-1] *(V[-1, tau+1] + V[-1, tau])
 
             # Solve the linear system for interior points
             V[1:-1, tau] = spsolve(lhs, rhs_vector)
 
-        return sol.SolutionBlackScholes(V,S,T)
+            # Calculate Greeks for interior points
+            delta[1:-1, tau] = (V[2:, tau] - V[:-2, tau]) / (2 * ds)
+            gamma[1:-1, tau] = (V[2:, tau] - 2 * V[1:-1, tau] + V[:-2, tau]) / (ds**2)
+            theta[1:-1, tau] = -0.5 * (self.equation.sigma**2) * (S[1:-1]**2) * gamma[1:-1, tau] - self.equation.rate * S[1:-1] * delta[1:-1, tau] + self.equation.rate * V[1:-1, tau]
+
+            delta, gamma, theta = self.__calculate_greeks_at_boundary(delta, gamma, theta, tau, V, S, ds)
+
+        return sol.SolutionBlackScholes(V,S,T, delta, gamma, theta)
+
+    def __calculate_greeks_at_boundary(self, delta, gamma, theta, tau, V, S, ds):
+        delta[0, tau] = (V[1, tau+1] - V[0, tau+1]) / ds
+        delta[self.equation.s_nodes, tau] = (V[self.equation.s_nodes, tau+1] - V[self.equation.s_nodes-1, tau+1]) / ds
+
+        gamma[0, tau] = (V[2, tau+1] - 2*V[1, tau+1] + V[0, tau+1]) / (ds**2)
+        gamma[self.equation.s_nodes, tau] = (V[self.equation.s_nodes, tau+1] - 2*V[self.equation.s_nodes-1, tau+1] + V[self.equation.s_nodes-2, tau+1]) / (ds**2)
+
+        theta[0, tau] = -0.5 * (self.equation.sigma**2) * (S[0]**2) * gamma[0, tau] - self.equation.rate * S[0] * delta[0, tau] + self.equation.rate * V[0, tau+1]
+        theta[self.equation.s_nodes, tau] = -0.5 * (self.equation.sigma**2) * (S[-1]**2) * gamma[self.equation.s_nodes, tau] - self.equation.rate * S[-1] * delta[self.equation.s_nodes, tau] + self.equation.rate * V[self.equation.s_nodes, tau+1]
+
+        return delta, gamma, theta

pdesolvers/solvers/heat_solvers.py

@@ -16,23 +16,29 @@ def solve(self):
         :return: the solver instance with the computed temperature values
         """
 
-        x = self.equation.generate_x_grid()
-        dx = x[1] - x[0]
+        x = self.equation.generate_grid(self.equation.length, self.equation.x_nodes)
+        t = self.equation.generate_grid(self.equation.time, self.equation.t_nodes)
 
+        dx = x[1] - x[0]
         dt_max = 0.5 * (dx**2) / self.equation.k
-        dt = 0.8 * dt_max
-        time_step = int(self.equation.time/dt)
-        self.equation.t_nodes = time_step
 
-        t = np.linspace(0, self.equation.time, self.equation.t_nodes)
+        if self.equation.t_nodes is None:
+            dt = 0.8 * dt_max
+            self.equation.t_nodes = int(self.equation.time/dt)
+            dt = self.equation.time / self.equation.t_nodes
+        else:
+            dt = t[1] - t[0]
 
-        u = np.zeros((time_step, self.equation.x_nodes))
+        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!")
+
+        u = np.zeros((self.equation.t_nodes, self.equation.x_nodes))
 
         u[0, :] = self.equation.get_initial_temp(x)
         u[:, 0] = self.equation.get_left_boundary(t)
         u[:, -1] = self.equation.get_right_boundary(t)
 
-        for tau in range(0, time_step-1):
+        for tau in range(0, self.equation.t_nodes-1):
             for i in range(1, self.equation.x_nodes - 1):
                 u[tau+1,i] = u[tau, i] + (dt * self.equation.k * (u[tau, i-1] - 2 * u[tau, i] + u[tau, i+1]) / dx**2)
 
@@ -49,8 +55,8 @@ def solve(self):
         :return: the solver instance with the computed temperature values
         """
 
-        x = self.equation.generate_x_grid()
-        t = self.equation.generate_t_grid()
+        x = self.equation.generate_grid(self.equation.length, self.equation.x_nodes)
+        t = self.equation.generate_grid(self.equation.time, self.equation.t_nodes)
 
         dx = x[1] - x[0]
         dt = t[1] - t[0]

pdesolvers/tests/__init__.py

@@ -19,7 +19,7 @@ def test_check_lower_boundary_for_call_explicit(self):
     def test_check_terminal_condition_for_call_explicit(self):
         result = solver.BlackScholesExplicitSolver(self.equation).solve().get_result()
 
-        test_asset_grid = self.equation.generate_asset_grid()
+        test_asset_grid = self.equation.generate_grid(self.equation.S_max, self.equation.s_nodes)
         test_strike_price = self.equation.strike_price
         expected_payoff = np.maximum(test_asset_grid - test_strike_price, 0)
 
@@ -29,7 +29,7 @@ def test_check_terminal_condition_for_put_explicit(self):
         self.equation.option_type = 'put'
         result = solver.BlackScholesExplicitSolver(self.equation).solve().get_result()
 
-        test_asset_grid = self.equation.generate_asset_grid()
+        test_asset_grid = self.equation.generate_grid(self.equation.S_max, self.equation.s_nodes)
         test_strike_price = self.equation.strike_price
         expected_payoff = np.maximum(test_strike_price - test_asset_grid, 0)
 
@@ -50,7 +50,7 @@ def test_check_lower_boundary_for_call_cn(self):
     def test_check_terminal_condition_for_call_cn(self):
         result = solver.BlackScholesCNSolver(self.equation).solve().get_result()
 
-        test_asset_grid = self.equation.generate_asset_grid()
+        test_asset_grid = self.equation.generate_grid(self.equation.S_max, self.equation.s_nodes)
         test_strike_price = self.equation.strike_price
         expected_payoff = np.maximum(test_asset_grid - test_strike_price, 0)
 
@@ -60,7 +60,7 @@ def test_check_terminal_condition_for_put_cn(self):
         self.equation.option_type = 'put'
         result = solver.BlackScholesCNSolver(self.equation).solve().get_result()
 
-        test_asset_grid = self.equation.generate_asset_grid()
+        test_asset_grid = self.equation.generate_grid(self.equation.S_max, self.equation.s_nodes)
         test_strike_price = self.equation.strike_price
         expected_payoff = np.maximum(test_strike_price - test_asset_grid, 0)
 
@@ -86,7 +86,7 @@ def test_convergence_between_interpolated_data(self):
 
         diff = np.abs(data1 - data2)
 
-        assert diff < 1e-4
+        assert np.max(diff) < 1e-4