Merge pull request #49 from chemardes/feature/implement-monte-carlo-simulations

FEATURE: adding monte carlo simulations

Author: chemardes

.github/workflows/ci.yml

@@ -12,7 +12,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: actions/setup-python@v5
         with:
-          python-version: '3.13'
+          python-version: '3.10'
           architecture: 'x64'
       - name: Run python tests
         run: |

ci/run_python_tests.sh

@@ -6,7 +6,7 @@ export PYTHONPATH=$(pwd)
 
 # Set up your Python environment
 pip install virtualenv
-virtualenv -p python3.13 venv
+virtualenv -p python3.10 venv
 source venv/bin/activate
 
 pip install -r requirements.txt

pdesolvers/__init__.py

@@ -2,4 +2,5 @@
 from .solution import *
 from .solvers import *
 from .utils import *
-from .enums import *
+from .enums import *
+from .optionspricing import *

pdesolvers/main.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 import pdesolvers as pde
 
 def main():
@@ -14,14 +15,36 @@ def main():
     # solver1 = pde.Heat1DCNSolver(equation1)
     # solver2 = pde.Heat1DExplicitSolver(equation1)
 
-    # testing for bse
-    equation2 = pde.BlackScholesEquation(pde.OptionType.EUROPEAN_CALL, 300, 1, 0.2, 0.05, 100, 100, 20000)
+    # testing for monte carlo pricing
+
+    ticker = 'AAPL'
+
+    # STOCK
+    historical_data = pde.HistoricalStockData(ticker)
+    historical_data.fetch_stock_data( "2024-03-21","2025-03-21")
+    sigma, r = historical_data.estimate_metrics()
+    current_price = historical_data.get_latest_stock_price()
+
+    equation2 = pde.BlackScholesEquation(pde.OptionType.EUROPEAN_CALL, current_price, 1, sigma, r, 100, 100, 20000)
 
     solver1 = pde.BlackScholesCNSolver(equation2)
     solver2 = pde.BlackScholesExplicitSolver(equation2)
     sol1 = solver1.solve()
-    sol1.plot_greek(pde.Greeks.GAMMA)
-    sol1.plot()
+    # sol1.plot()
+
+    # COMPARISON
+    #  look to see the corresponding option price for the expiration date and strike price
+    pricing_1 = pde.BlackScholesFormula(pde.OptionType.EUROPEAN_CALL, current_price, 100, r, sigma, 1)
+    pricing_2 = pde.MonteCarloPricing(pde.OptionType.EUROPEAN_CALL, current_price, 100, r, sigma, 1, 365, 10000)
+
+    bs_price = pricing_1.get_black_scholes_merton_price()
+    monte_carlo_price = pricing_2.get_monte_carlo_option_price()
+
+    pde_price = sol1.get_result()[-1, 0]
+    print(f"PDE Price: {pde_price}")
+    print(f"Black-Scholes Price: {bs_price}")
+    print(f"Monte-Carlo Price: {monte_carlo_price}")
+    pricing_2.plot()
 
 
 if __name__ == "__main__":

pdesolvers/optionspricing/__init__.py

@@ -0,0 +1,3 @@
+from .black_scholes_formula import BlackScholesFormula
+from .monte_carlo import MonteCarloPricing
+from .market_data import HistoricalStockData

pdesolvers/optionspricing/black_scholes_formula.py

@@ -0,0 +1,26 @@
+import numpy as np
+from scipy.stats import norm
+from pdesolvers.enums.enums import OptionType
+
+class BlackScholesFormula:
+    def __init__(self, option_type: OptionType, S0, strike_price, r, sigma, expiry):
+        self.__option_type = option_type
+        self.__S0 = S0
+        self.__strike_price = strike_price
+        self.__r = r
+        self.__sigma = sigma
+        self.__expiry = expiry
+
+    def get_black_scholes_merton_price(self):
+
+        d1 = (np.log(self.__S0/self.__strike_price) + (self.__r + 0.5 * self.__sigma**2) * self.__expiry) / (self.__sigma * np.sqrt(self.__expiry))
+        d2 = d1 - (self.__sigma * np.sqrt(self.__expiry))
+
+        if self.__option_type == OptionType.EUROPEAN_CALL:
+            payoff = self.__S0 * norm.cdf(d1) - self.__strike_price * np.exp(-self.__r * self.__expiry) * norm.cdf(d2)
+        elif self.__option_type == OptionType.EUROPEAN_PUT:
+            payoff =  self.__strike_price * np.exp(-self.__r * self.__expiry) * norm.cdf(-d2) - self.__S0 * norm.cdf(-d1)
+        else:
+            raise ValueError(f'Unsupported option type: {self.__option_type}')
+
+        return payoff

pdesolvers/optionspricing/market_data.py

@@ -0,0 +1,36 @@
+import numpy as np
+import yfinance as yf
+
+class HistoricalStockData:
+    def __init__(self, ticker):
+        self.__stock_data = None
+        self.__options = None
+        self.__ticker = ticker
+
+    def fetch_stock_data(self, start_date, end_date):
+        self.__stock_data = yf.download(self.__ticker, start=start_date, end=end_date, interval='1d')
+        return self.__stock_data
+
+    def estimate_metrics(self):
+
+        if self.__stock_data is None:
+            raise ValueError(f'No data available - data must be fetched first.')
+
+        self.__stock_data["Log Returns"] = np.log(self.__stock_data["Close"] / self.__stock_data["Close"].shift(1))
+
+        # takes into account the gaps between trading days
+        self.__stock_data["Time Diff"] = self.__stock_data.index.to_series().diff().dt.days
+        self.__stock_data["Z"] = self.__stock_data["Log Returns"] / self.__stock_data["Time Diff"]
+
+        print(self.__stock_data)
+
+        annualized_sigma = self.__stock_data["Z"].std() * np.sqrt(252)
+        annualized_mu = self.__stock_data["Z"].mean() * 252
+
+        return annualized_sigma, annualized_mu
+
+    def get_latest_stock_price(self):
+        if self.__stock_data is None:
+            raise ValueError("No data available. Call fetch_data first.")
+
+        return self.__stock_data["Close"].iloc[-1].item()

pdesolvers/optionspricing/monte_carlo.py

@@ -0,0 +1,98 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pdesolvers.enums.enums import OptionType
+
+class MonteCarloPricing:
+
+    def __init__(self, option_type: OptionType, S0, strike_price, r, sigma, T, time_steps, sim):
+        """
+        Initialize the Geometric Brownian Motion model with the given parameters.
+
+        Parameters:
+        - S0: Initial stock price
+        - mu: Drift coefficient (expected return)
+        - sigma: Volatility coefficient (standard deviation of returns)
+        - T: Time period for the simulation (in years)
+        - time_steps: Number of time steps in the simulation
+        - sim: Number of simulations to run
+        """
+
+        self.__option_type = option_type
+        self.__S0 = S0
+        self.__strike_price = strike_price
+        self.__r = r
+        self.__sigma = sigma
+        self.__T = T
+        self.__time_steps = time_steps
+        self.__sim = sim
+        self.__S = None
+
+    def get_monte_carlo_option_price(self):
+
+        S = self.simulate_gbm()
+
+        if self.__option_type == OptionType.EUROPEAN_CALL:
+            payoff = np.maximum(S[:, -1] - self.__strike_price, 0)
+        elif self.__option_type == OptionType.EUROPEAN_PUT:
+            payoff = np.maximum(self.__strike_price - S[:, -1], 0)
+        else:
+            raise ValueError(f'Unsupported Option Type: {self.__option_type}')
+
+        option_price = np.exp(-self.__r * self.__T) * np.mean(payoff)
+
+        return option_price
+
+    def simulate_gbm(self):
+        """
+        Simulate the Geometric Brownian Motion for the given parameters.
+
+        This method calculates the stock prices at each time step for each simulation.
+        """
+
+        t = self.__generate_grid()
+        dt = t[1] - t[0]
+
+        B = np.zeros((self.__sim, self.__time_steps))
+        S = np.zeros((self.__sim, self.__time_steps))
+
+        # for all simulations at t = 0
+        S[:,0] = self.__S0
+        Z = np.random.normal(0, 1, (self.__sim, self.__time_steps))
+
+        for i in range(self.__sim):
+            for j in range (1, self.__time_steps):
+                # updates brownian motion
+                B[i,j] = B[i,j-1] + np.sqrt(dt) * Z[i,j-1]
+                # calculates stock price based on the incremental difference
+                S[i,j] = S[i, j-1] * np.exp((self.__r - 0.5*self.__sigma**2)*dt + self.__sigma * (B[i, j] - B[i, j - 1]))
+
+        self.__S = S
+
+        return self.__S
+
+    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 plot(self):
+        """
+        Plot the simulated stock prices for all simulations.
+        """
+
+        t = self.__generate_grid()
+
+        fig = plt.figure(figsize=(10,6))
+        for i in range(np.min([100, self.__sim])):
+            plt.plot(t, self.__S[i], color='grey', alpha=0.3)
+
+        plt.title("Simulated Geometric Brownian Motion")
+        plt.xlabel("Time (Years)")
+        plt.ylabel("Stock Price")
+        plt.show()

pdesolvers/utils/__init__.py

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

pdesolvers/utils/utility.py

@@ -1,4 +1,6 @@
 import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
 from scipy.sparse import csc_matrix
 
 import pdesolvers.enums.enums as enum
@@ -152,4 +154,42 @@ def interpolate(self, x, y):
         interpolated += rbf_weights[3] * self.__z[i_minus + 1, j_minus + 1]
         interpolated /= sum_rbf
 
-        return interpolated
+        return interpolated
+
+class GPUResults:
+
+    def __init__(self, file_path, s_max, expiry):
+        self.__file_path = file_path
+        self.__s_max = s_max
+        self.__expiry = expiry
+        self.__grid_data = None
+
+    def get_results(self):
+
+        # Load data
+        df = pd.read_csv(self.__file_path, header=None)
+        print(f"Data shape: {df.shape}")
+
+        self.__grid_data = df.values.T
+
+        return self.__grid_data
+
+    def plot_option_surface(self):
+        if self.__grid_data is None:
+            self.get_results()
+
+        price_grid = np.linspace(0, self.__s_max, self.__grid_data.shape[0])
+        time_grid = np.linspace(0, self.__expiry, self.__grid_data.shape[1])
+        X, Y = np.meshgrid(time_grid, price_grid)
+
+        fig = plt.figure(figsize=(10, 6))
+        ax = fig.add_subplot(111, projection='3d')
+        surf = ax.plot_surface(X, Y, self.__grid_data, cmap='viridis')
+
+        ax.set_xlabel('Time')
+        ax.set_ylabel('Asset Price')
+        ax.set_zlabel('Option Value')
+        ax.set_title('Option Value Surface Plot')
+        fig.colorbar(surf, shrink=0.5, aspect=5)
+
+        plt.show()