Skip to content

Commit 2279864

Browse files
authored
Merge pull request #48 from chemardes/feature/implementing-interpolation-between-grids
FEATURE: adding interpolation between two grids + refactoring Solution
2 parents 53eb1d2 + acd56df commit 2279864

11 files changed

Lines changed: 197 additions & 168 deletions

File tree

pdesolvers/enums/__init__.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
from .option_type import OptionType
1+
from .enums import OptionType, Greeks

pdesolvers/enums/enums.py

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
from enum import Enum
2+
3+
class OptionType(Enum):
4+
EUROPEAN_CALL = 'European Call'
5+
EUROPEAN_PUT = 'European Put'
6+
7+
class Greeks(Enum):
8+
DELTA = 'Delta'
9+
GAMMA = 'Gamma'
10+
THETA = 'Theta'

pdesolvers/enums/option_type.py

Lines changed: 0 additions & 5 deletions
This file was deleted.

pdesolvers/main.py

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -18,10 +18,9 @@ def main():
1818
equation2 = pde.BlackScholesEquation(pde.OptionType.EUROPEAN_CALL, 300, 1, 0.2, 0.05, 100, 100, 20000)
1919

2020
solver1 = pde.BlackScholesCNSolver(equation2)
21-
# solver2 = pde.BlackScholesExplicitSolver(equation2)
21+
solver2 = pde.BlackScholesExplicitSolver(equation2)
2222
sol1 = solver1.solve()
23-
24-
sol1.plot_greek('gamma')
23+
sol1.plot_greek(pde.Greeks.GAMMA)
2524
sol1.plot()
2625

2726

pdesolvers/pdes/black_scholes.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
import numpy as np
2-
from pdesolvers.enums.option_type import OptionType
2+
from pdesolvers.enums.enums import OptionType
33

44
class BlackScholesEquation:
55

pdesolvers/solution/solution.py

Lines changed: 77 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -1,35 +1,33 @@
11
import numpy as np
2+
import pdesolvers.utils.utility as utility
3+
import pdesolvers.enums.enums as enum
4+
25
from matplotlib import pyplot as plt
36

4-
class Solution1D:
7+
class Solution:
58

6-
def __init__(self, result, x_grid, t_grid):
9+
def __init__(self, result, x_grid, y_grid, dx, dy):
710
self.result = result
811
self.x_grid = x_grid
9-
self.t_grid = t_grid
12+
self.y_grid = y_grid
13+
self.dx = dx
14+
self.dy = dy
1015

1116
def plot(self):
1217
"""
13-
Generates a 3D surface plot of the temperature distribution across a grid of space and time
18+
Generates a 3D surface plot of the option values across a grid of asset prices and time
1419
1520
:return: 3D surface plot
1621
"""
1722

18-
if self.result is None:
19-
raise RuntimeError("Solution has not been computed - please run the solver.")
20-
21-
x_plot, t_plot = np.meshgrid(self.x_grid,self.t_grid)
23+
X, Y = np.meshgrid(self.x_grid, self.y_grid)
2224

2325
# plotting the 3d surface
2426
fig = plt.figure(figsize=(10,6))
2527
ax = fig.add_subplot(111, projection='3d')
26-
surf = ax.plot_surface(x_plot, t_plot, self.result, cmap='viridis')
28+
surf = ax.plot_surface(X, Y, self.result, cmap='viridis')
2729

28-
# set labels and title
29-
ax.set_xlabel('Space')
30-
ax.set_ylabel('Time')
31-
ax.set_zlabel('Temperature')
32-
ax.set_title('3D Surface Plot of 1D Heat Equation')
30+
self._set_plot_labels(ax)
3331

3432
fig.colorbar(surf, shrink=0.5, aspect=5)
3533
plt.show()
@@ -42,66 +40,93 @@ def get_result(self):
4240
"""
4341
return self.result
4442

45-
class SolutionBlackScholes:
46-
def __init__(self, result, s_grid, t_grid, delta, gamma, theta):
47-
self.result = result
48-
self.s_grid = s_grid
49-
self.t_grid = t_grid
50-
self.delta = delta
51-
self.gamma = gamma
52-
self.theta = theta
43+
def _set_plot_labels(self, ax):
44+
ax.set_xlabel('X-axis')
45+
ax.set_ylabel('Y-axis')
46+
ax.set_zlabel('Value')
47+
ax.set_title('3D Surface Plot')
5348

54-
def plot(self):
49+
def __sub__(self, other):
5550
"""
56-
Generates a 3D surface plot of the option values across a grid of asset prices and time
57-
58-
:return: 3D surface plot
51+
Compares two solutions by interpolating the sparse grid to the dense grid and computing the difference
52+
:param other: the grid to be compared against
53+
:return: the error (difference) between the two solutions
5954
"""
6055

61-
X, Y = np.meshgrid(self.t_grid, self.s_grid)
56+
sparser_grid = other
57+
denser_grid = self
6258

63-
# plotting the 3d surface
64-
fig = plt.figure(figsize=(10,6))
65-
ax = fig.add_subplot(111, projection='3d')
66-
surf = ax.plot_surface(X, Y, self.result, cmap='viridis')
59+
if self.x_grid.shape[0] < other.x_grid.shape[0] and self.y_grid.shape[0] < other.y_grid.shape[0]:
60+
sparser_grid = self
61+
denser_grid = other
6762

68-
ax.set_xlabel('Time')
69-
ax.set_ylabel('Asset Price')
70-
ax.set_zlabel('Option Value')
71-
ax.set_title('Option Value Surface Plot')
63+
interpolator_sparse = utility.RBFInterpolator(sparser_grid.result.T, sparser_grid.dx, sparser_grid.dy)
7264

73-
fig.colorbar(surf, shrink=0.5, aspect=5)
74-
plt.show()
65+
diff = 0
7566

76-
def get_result(self):
77-
"""
78-
Gets the grid of computed option prices
67+
for idx_x in range(denser_grid.x_grid.shape[0]):
68+
for idx_y in range(denser_grid.y_grid.shape[0]):
7969

80-
:return: grid result
81-
"""
82-
return self.result
70+
# Points (x, y) of the dense grid
71+
x = denser_grid.x_grid[idx_x]
72+
y = denser_grid.y_grid[idx_y]
73+
74+
# Value at (x, y) for the dense grid
75+
val_dense_x_y = denser_grid.result[idx_y, idx_x]
8376

84-
def plot_greek(self, greek_type='delta', time_step=0):
77+
# Interpolate the sparse grid at (x, y)
78+
val_sparse_x_y = interpolator_sparse.interpolate(x, y)
79+
80+
diff = np.max([diff, np.abs(val_dense_x_y - val_sparse_x_y)])
81+
82+
return diff
83+
84+
class Solution1D(Solution):
85+
86+
def __init__(self, result, x_grid, y_grid, dx, dy):
87+
super().__init__(result, x_grid, y_grid, dx, dy)
88+
89+
def _set_plot_labels(self, ax):
90+
ax.set_xlabel('Space')
91+
ax.set_ylabel('Time')
92+
ax.set_zlabel('Temperature')
93+
ax.set_title('3D Surface Plot of 1D Heat Equation')
94+
95+
96+
class SolutionBlackScholes(Solution):
97+
def __init__(self, result, x_grid, y_grid, dx, dy, delta, gamma, theta, option_type):
98+
super().__init__(result, x_grid, y_grid, dx, dy)
99+
self.option_type = option_type
100+
self.delta = delta
101+
self.gamma = gamma
102+
self.theta = theta
103+
104+
def plot_greek(self, greek_type=enum.Greeks.DELTA, time_step=0):
85105

86106
greek_types = {
87-
'delta': {'data': self.delta, 'title': 'Delta'},
88-
'gamma': {'data': self.gamma, 'title': 'Gamma'},
89-
'theta': {'data': self.theta, 'title': 'Theta'}
107+
enum.Greeks.DELTA : {'data': self.delta, 'title': enum.Greeks.DELTA.value},
108+
enum.Greeks.GAMMA : {'data': self.gamma, 'title': enum.Greeks.DELTA.value},
109+
enum.Greeks.THETA : {'data': self.theta, 'title': enum.Greeks.DELTA.value}
90110
}
91111

92-
if greek_type.lower() not in greek_types:
93-
raise ValueError("Invalid greek type - please choose between delta/gamma/theta.")
112+
# if greek_type.lower() not in greek_types:
113+
# raise ValueError("Invalid greek type - please choose between delta/gamma/theta.")
94114

95-
chosen_greek = greek_types[greek_type.lower()]
115+
chosen_greek = greek_types[greek_type]
96116
greek_data = chosen_greek['data'][:, time_step]
97117
plt.figure(figsize=(8, 6))
98-
plt.plot(self.s_grid, greek_data, label=f"Delta at t={self.t_grid[time_step]:.4f}", color="blue")
118+
plt.plot(self.y_grid, greek_data, label=f"Delta at t={self.x_grid[time_step]:.4f}", color="blue")
99119

100-
plt.title(f"{chosen_greek['title']} vs. Stock Price at t={self.t_grid[time_step]:.4f}")
120+
plt.title(f"{chosen_greek['title']} vs. Stock Price at t={self.x_grid[time_step]:.4f}")
101121
plt.xlabel("Stock Price (S)")
102122
plt.ylabel(chosen_greek['title'])
103123
plt.grid()
104124
plt.legend()
105125

106126
plt.show()
107127

128+
def _set_plot_labels(self, ax):
129+
ax.set_xlabel('Time')
130+
ax.set_ylabel('Asset Price')
131+
ax.set_zlabel(f'{self.option_type.value} Option Value')
132+
ax.set_title(f'{self.option_type.value} Option Value Surface Plot')

pdesolvers/solvers/black_scholes_solvers.py

Lines changed: 10 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,11 @@
1-
from scipy import sparse
2-
from scipy.sparse.linalg import spsolve
31
import numpy as np
42
import pdesolvers.solution as sol
53
import pdesolvers.pdes.black_scholes as bse
6-
import pdesolvers.enums.option_type as enum
4+
import pdesolvers.enums.enums as enum
5+
import pdesolvers.utils.utility as utility
6+
7+
from scipy import sparse
8+
from scipy.sparse.linalg import spsolve
79

810
class BlackScholesExplicitSolver:
911

@@ -56,46 +58,13 @@ def solve(self):
5658
V[i, tau] = V[i, tau + 1] - (theta[i, tau] * dt)
5759

5860
# setting boundary conditions
59-
lower, upper = self.__set_boundary_conditions(T, tau)
61+
lower, upper = utility.BlackScholesHelper.set_boundary_conditions(self.equation, T, tau)
6062
V[0, tau] = lower
6163
V[self.equation.s_nodes, tau] = upper
6264

63-
delta, gamma, theta = self.__calculate_greeks_at_boundary(delta, gamma, theta, tau, V, S, ds)
64-
65-
return sol.SolutionBlackScholes(V,S,T, delta, gamma, theta)
66-
67-
def __set_boundary_conditions(self, T, tau):
68-
"""
69-
Sets the boundary conditions for the Black-Scholes Equation based on option type
70-
71-
:param T: grid of time steps
72-
:param tau: index of current time step
73-
:return: a tuple representing the boundary values for the given time step
74-
"""
75-
76-
lower_boundary = None
77-
upper_boundary = None
78-
if self.equation.option_type == enum.OptionType.EUROPEAN_CALL:
79-
lower_boundary = 0
80-
upper_boundary = self.equation.S_max - self.equation.strike_price * np.exp(-self.equation.rate * (self.equation.expiry - T[tau]))
81-
elif self.equation.option_type == enum.OptionType.EUROPEAN_PUT:
82-
lower_boundary = self.equation.strike_price * np.exp(-self.equation.rate * (self.equation.expiry - T[tau]))
83-
upper_boundary = 0
84-
85-
return lower_boundary, upper_boundary
86-
87-
def __calculate_greeks_at_boundary(self, delta, gamma, theta, tau, V, S, ds):
88-
delta[0, tau] = (V[1, tau+1] - V[0, tau+1]) / ds # Forward difference for lower boundary
89-
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
65+
delta, gamma, theta = utility.BlackScholesHelper.calculate_greeks_at_boundary(self.equation, delta, gamma, theta, tau, V, S, ds)
9066

91-
gamma[0, tau] = (V[2, tau+1] - 2*V[1, tau+1] + V[0, tau+1]) / (ds**2) # Forward approximation
92-
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
93-
94-
# Calculate theta for boundary points using the same formula
95-
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]
96-
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]
97-
98-
return delta, gamma, theta
67+
return sol.SolutionBlackScholes(V, T, S, dt, ds, delta, gamma, theta, self.equation.option_type)
9968

10069

10170
class BlackScholesCNSolver:
@@ -159,18 +128,7 @@ def solve(self):
159128
gamma[1:-1, tau] = (V[2:, tau] - 2 * V[1:-1, tau] + V[:-2, tau]) / (ds**2)
160129
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]
161130

162-
delta, gamma, theta = self.__calculate_greeks_at_boundary(delta, gamma, theta, tau, V, S, ds)
163-
164-
return sol.SolutionBlackScholes(V,S,T, delta, gamma, theta)
165-
166-
def __calculate_greeks_at_boundary(self, delta, gamma, theta, tau, V, S, ds):
167-
delta[0, tau] = (V[1, tau+1] - V[0, tau+1]) / ds
168-
delta[self.equation.s_nodes, tau] = (V[self.equation.s_nodes, tau+1] - V[self.equation.s_nodes-1, tau+1]) / ds
169-
170-
gamma[0, tau] = (V[2, tau+1] - 2*V[1, tau+1] + V[0, tau+1]) / (ds**2)
171-
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)
131+
delta, gamma, theta = utility.BlackScholesHelper.calculate_greeks_at_boundary(self.equation, delta, gamma, theta, tau, V, S, ds)
172132

173-
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]
174-
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]
133+
return sol.SolutionBlackScholes(V, T, S, dt, ds, delta, gamma, theta, self.equation.option_type)
175134

176-
return delta, gamma, theta

pdesolvers/solvers/heat_solvers.py

Lines changed: 6 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
import numpy as np
2-
from scipy.sparse.linalg import spsolve
3-
from scipy.sparse import csc_matrix
4-
52
import pdesolvers.solution as sol
63
import pdesolvers.pdes.heat_1d as heat
4+
import pdesolvers.utils.utility as utility
5+
6+
from scipy.sparse.linalg import spsolve
77

88
class Heat1DExplicitSolver:
99
def __init__(self, equation: heat.HeatEquation):
@@ -42,7 +42,7 @@ def solve(self):
4242
for i in range(1, self.equation.x_nodes - 1):
4343
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)
4444

45-
return sol.Solution1D(u, x, t)
45+
return sol.Solution1D(u, x, t, dx, dt)
4646

4747
class Heat1DCNSolver:
4848
def __init__(self, equation: heat.HeatEquation):
@@ -72,7 +72,7 @@ def solve(self):
7272
u[:, 0] = self.equation.get_left_boundary(t)
7373
u[:, -1] = self.equation.get_right_boundary(t)
7474

75-
lhs = self.__build_tridiagonal_matrix(a, b, c, self.equation.x_nodes - 2)
75+
lhs = utility.Heat1DHelper.build_tridiagonal_matrix(a, b, c, self.equation.x_nodes - 2)
7676
rhs = np.zeros(self.equation.x_nodes - 2)
7777

7878
for tau in range(0, self.equation.t_nodes - 1):
@@ -85,25 +85,4 @@ def solve(self):
8585

8686
u[tau+1, 1:-1] = spsolve(lhs, rhs)
8787

88-
return sol.Solution1D(u, x, t)
89-
90-
@staticmethod
91-
def __build_tridiagonal_matrix(a, b, c, nodes):
92-
"""
93-
Initialises the tridiagonal matrix on the LHS of the equation
94-
95-
:param a: the coefficient of U @ (t = tau + 1 & x = i-1)
96-
:param b: the coefficient of U @ (t = tau + 1 & x = i)
97-
:param c: the coefficient of U @ (t = tau + 1 & x = i+1)
98-
:param nodes: number of spatial nodes ( used to initialise the size of the tridiagonal matrix)
99-
:return: the tridiagonal matrix consisting of coefficients
100-
"""
101-
102-
matrix = np.zeros((nodes, nodes))
103-
np.fill_diagonal(matrix, b)
104-
np.fill_diagonal(matrix[1:], a)
105-
np.fill_diagonal(matrix[:, 1:], c)
106-
107-
matrix = csc_matrix(matrix)
108-
109-
return matrix
88+
return sol.Solution1D(u, x, t, dx, dt)

0 commit comments

Comments
 (0)