FEATURE: fixing greek plots + adding execution times

Author: chemardes

pdesolvers/optionspricing/brownian_motion.py

@@ -0,0 +1,61 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+class BrownianMotion:
+
+    def __init__(self, T, time_steps, sim, mu=0, sigma=1):
+        self.__mu = mu
+        self.__sigma = sigma
+        self.__T = T
+        self.__time_steps = time_steps
+        self.__sim = sim
+        self.__B = None
+
+    def simulate_bm(self):
+        t = self.__generate_grid()
+        dt = t[1] - t[0]
+
+        Z = np.random.normal(0,1, (self.__sim, self.__time_steps))
+        B = np.zeros((self.__sim, self.__time_steps))
+
+        for i in range(self.__sim):
+            for j in range(1, self.__time_steps):
+                B[i,j] = B[i,j-1] + self.__mu * dt + self.__sigma * np.sqrt(dt) * Z[i,j]
+
+        self.__B = B
+        return self
+
+
+    def __generate_grid(self):
+        """
+        Generate a time grid from 0 to T with `time_steps` intervals.
+
+        Returns:
+        - A numpy array representing the time grid.
+        """
+
+        return np.linspace(0, self.__T, self.__time_steps)
+
+    def __get_paths(self):
+        return self.__B
+
+    def plot(self):
+        t = self.__generate_grid()
+        B = self.__get_paths()
+
+        plt.figure(figsize=(10,6))
+        for i in range(self.__sim):
+            plt.plot(t, B[i], color='grey', alpha=0.3)
+
+        plt.title("Simulated Brownian Motion Paths")
+        plt.xlabel("Time (Years)")
+        plt.ylabel("B(t)")
+        # plt.show()
+        plt.savefig("brownian_plot.pdf", bbox_inches='tight')
+
+
+def main():
+    BrownianMotion(T=1, time_steps=365, sim=200).simulate_bm().plot()
+
+if __name__ == "__main__":
+    main()

pdesolvers/solution/solution.py

@@ -6,12 +6,13 @@
 
 class Solution:
 
-    def __init__(self, result, x_grid, y_grid, dx, dy):
+    def __init__(self, result, x_grid, y_grid, dx, dy, duration):
         self.result = result
         self.x_grid = x_grid
         self.y_grid = y_grid
         self.dx = dx
         self.dy = dy
+        self.duration = duration
 
     def plot(self):
         """
@@ -46,6 +47,9 @@ def _set_plot_labels(self, ax):
         ax.set_zlabel('Value')
         ax.set_title('3D Surface Plot')
 
+    def get_execution_time(self):
+        return self.duration
+
     def __sub__(self, other):
         """
         Compares two solutions by interpolating the sparse grid to the dense grid and computing the difference
@@ -83,8 +87,8 @@ def __sub__(self, other):
 
 class Solution1D(Solution):
 
-    def __init__(self, result, x_grid, y_grid, dx, dy):
-        super().__init__(result, x_grid, y_grid, dx, dy)
+    def __init__(self, result, x_grid, y_grid, dx, dy, duration):
+        super().__init__(result, x_grid, y_grid, dx, dy, duration)
 
     def _set_plot_labels(self, ax):
         ax.set_xlabel('Space')
@@ -94,8 +98,8 @@ def _set_plot_labels(self, ax):
 
 
 class SolutionBlackScholes(Solution):
-    def __init__(self, result, x_grid, y_grid, dx, dy, delta, gamma, theta, option_type):
-        super().__init__(result, x_grid, y_grid, dx, dy)
+    def __init__(self, result, x_grid, y_grid, dx, dy, duration, delta, gamma, theta, option_type):
+        super().__init__(result, x_grid, y_grid, dx, dy, duration)
         self.option_type = option_type
         self.delta = delta
         self.gamma = gamma
@@ -105,8 +109,8 @@ def plot_greek(self, greek_type=enum.Greeks.DELTA, time_step=0):
 
         greek_types = {
             enum.Greeks.DELTA : {'data': self.delta, 'title': enum.Greeks.DELTA.value},
-            enum.Greeks.GAMMA : {'data': self.gamma, 'title': enum.Greeks.DELTA.value},
-            enum.Greeks.THETA : {'data': self.theta, 'title': enum.Greeks.DELTA.value}
+            enum.Greeks.GAMMA : {'data': self.gamma, 'title': enum.Greeks.GAMMA.value},
+            enum.Greeks.THETA : {'data': self.theta, 'title': enum.Greeks.THETA.value}
         }
 
         # if greek_type.lower() not in greek_types:
@@ -115,14 +119,15 @@ def plot_greek(self, greek_type=enum.Greeks.DELTA, time_step=0):
         chosen_greek = greek_types[greek_type]
         greek_data = chosen_greek['data'][:, time_step]
         plt.figure(figsize=(8, 6))
-        plt.plot(self.y_grid, greek_data, label=f"Delta at t={self.x_grid[time_step]:.4f}", color="blue")
+        plt.plot(self.y_grid, greek_data, label=f"{chosen_greek['title']} at t={self.x_grid[time_step]:.4f}", color="grey", alpha=0.8)
 
         plt.title(f"{chosen_greek['title']} vs. Stock Price at t={self.x_grid[time_step]:.4f}")
         plt.xlabel("Stock Price (S)")
         plt.ylabel(chosen_greek['title'])
         plt.grid()
         plt.legend()
 
+        np.savetxt(f"option_{chosen_greek['title']}.csv", np.column_stack((self.y_grid, greek_data)), delimiter=',', header=f"StockPrice,{chosen_greek['title']}", comments='')
         plt.show()
 
     def _set_plot_labels(self, ax):

pdesolvers/solvers/black_scholes_solvers.py

@@ -3,6 +3,7 @@
 import pdesolvers.pdes.black_scholes as bse
 import pdesolvers.enums.enums as enum
 import pdesolvers.utils.utility as utility
+import time
 
 from scipy import sparse
 from scipy.sparse.linalg import spsolve
@@ -19,6 +20,8 @@ def solve(self):
         :return: the solver instance with the computed option values
         """
 
+        start = time.perf_counter()
+
         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)
 
@@ -64,7 +67,10 @@ def solve(self):
 
             delta, gamma, theta = utility.BlackScholesHelper.calculate_greeks_at_boundary(self.equation, delta, gamma, theta, tau, V, S, ds)
 
-        return sol.SolutionBlackScholes(V, T, S, dt, ds, delta, gamma, theta, self.equation.option_type)
+        end = time.perf_counter()
+        duration = end - start
+
+        return sol.SolutionBlackScholes(V, T, S, dt, ds, duration, delta, gamma, theta, self.equation.option_type)
 
 
 class BlackScholesCNSolver:
@@ -79,6 +85,7 @@ def solve(self):
         :return: the solver instance with the computed option values
         """
 
+        start = time.perf_counter()
         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)
 
@@ -130,5 +137,8 @@ def solve(self):
 
             delta, gamma, theta = utility.BlackScholesHelper.calculate_greeks_at_boundary(self.equation, delta, gamma, theta, tau, V, S, ds)
 
-        return sol.SolutionBlackScholes(V, T, S, dt, ds, delta, gamma, theta, self.equation.option_type)
+        end = time.perf_counter()
+        duration = end - start
+
+        return sol.SolutionBlackScholes(V, T, S, dt, ds, duration, delta, gamma, theta, self.equation.option_type)
 

pdesolvers/solvers/heat_solvers.py

@@ -1,3 +1,5 @@
+import time
+
 import numpy as np
 import pdesolvers.solution as sol
 import pdesolvers.pdes.heat_1d as heat
@@ -16,6 +18,8 @@ def solve(self):
         :return: the solver instance with the computed temperature values
         """
 
+        start = time.perf_counter()
+
         x = self.equation.generate_grid(self.equation.length, self.equation.x_nodes)
         t = self.equation.generate_grid(self.equation.time, self.equation.t_nodes)
 
@@ -42,7 +46,10 @@ def solve(self):
             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)
 
-        return sol.Solution1D(u, x, t, dx, dt)
+        end = time.perf_counter()
+        duration = end - start
+
+        return sol.Solution1D(u, x, t, dx, dt, duration)
 
 class Heat1DCNSolver:
     def __init__(self, equation: heat.HeatEquation):
@@ -55,6 +62,8 @@ def solve(self):
         :return: the solver instance with the computed temperature values
         """
 
+        start = time.perf_counter()
+
         x = self.equation.generate_grid(self.equation.length, self.equation.x_nodes)
         t = self.equation.generate_grid(self.equation.time, self.equation.t_nodes)
 
@@ -85,4 +94,7 @@ def solve(self):
 
             u[tau+1, 1:-1] = spsolve(lhs, rhs)
 
-        return sol.Solution1D(u, x, t, dx, dt)
+        end = time.perf_counter()
+        duration = end - start
+
+        return sol.Solution1D(u, x, t, dx, dt, duration)

playground/geometric_brownian_motion.py

@@ -82,7 +82,7 @@ def plot(self):
 
         fig = plt.figure(figsize=(10,6))
         for i in range(self.__sim):
-            plt.plot(t, S[i])
+            plt.plot(t, S[i], color='grey', alpha=0.3)
 
         plt.title("Simulated Geometric Brownian Motion")
         plt.xlabel("Time (Years)")
@@ -91,7 +91,7 @@ def plot(self):
 
 
 def main():
-    GeometricBrownianMotion(100, 0.05, 0.03, 1, 365, 100).simulate_gbm().plot()
+    GeometricBrownianMotion(100, 2, 0.3, 1, 365, 100).simulate_gbm().plot()
 
 if __name__ == "__main__":
     main()