diff --git a/docs/demos/pinn_forward.rst b/docs/demos/pinn_forward.rst index 5a7df612b..a2d750677 100644 --- a/docs/demos/pinn_forward.rst +++ b/docs/demos/pinn_forward.rst @@ -52,6 +52,7 @@ Time-dependent PDEs pinn_forward/burgers.rar pinn_forward/allen.cahn pinn_forward/klein.gordon + pinn_forward/black_scholes pinn_forward/american.put - `Beltrami flow `_ diff --git a/docs/demos/pinn_forward/black_scholes.rst b/docs/demos/pinn_forward/black_scholes.rst new file mode 100644 index 000000000..62acf3dc5 --- /dev/null +++ b/docs/demos/pinn_forward/black_scholes.rst @@ -0,0 +1,168 @@ +Black-Scholes Equation +====================== + +Problem Setup +------------- + +We solve the Black-Scholes partial differential equation for pricing a European call option. +The problem is formulated in terms of **time-to-maturity** :math:`\tau = T - t`. + +The PDE is given by: + +.. math:: + \frac{\partial V}{\partial \tau} = \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + r S \frac{\partial V}{\partial S} - r V, \qquad S \in [0, S_{max}], \quad \tau \in [0, T] + +where: + - :math:`V(S, \tau)` is the option price + - :math:`S` is the underlying asset price + - :math:`\tau` is the time to maturity + - :math:`\sigma=0.2` is the volatility + - :math:`r=0.05` is the risk-free interest rate + +The boundary and initial conditions are: + +1. **Initial Condition** (at maturity :math:`\tau=0`): + +.. math:: V(S, 0) = \max(S - K, 0) + +where :math:`K=50` is the strike price. + +2. **Boundary Conditions**: + +.. math:: + V(0, \tau) &= 0 \\ + V(S_{max}, \tau) &= S_{max} - K e^{-r\tau} + +Implementation +-------------- + +This example uses a Physics-Informed Neural Network (PINN) to solve the equation. + +First, we import DeepXDE and define the problem parameters: + +.. code-block:: python + + import deepxde as dde + import numpy as np + + K = 50.0 # Strike price + r = 0.05 # Risk-free rate + sigma = 0.2 # Volatility + T = 1.0 # Time to maturity + S_max = 150.0 # Maximum stock price + +We define the computational geometry. We use ``GeometryXTime`` to combine the asset price interval :math:`[0, 150]` and the time domain :math:`[0, 1]`: + +.. code-block:: python + + geom = dde.geometry.Interval(0, S_max) + timedomain = dde.geometry.TimeDomain(0, T) + geomtime = dde.geometry.GeometryXTime(geom, timedomain) + +The PDE residual is defined below. Here ``x[:, 0]`` is the asset price +:math:`S`, ``x[:, 1]`` is the time to maturity :math:`\tau`, and ``y`` +is the option price :math:`V(S, \tau)` output by the network. + +.. code-block:: python + + def pde(x, y): + # x[:, 0]: asset price S + # x[:, 1]: time to maturity tau + # y: option price V(S, tau) + S = x[:, 0:1] + dV_tau = dde.grad.jacobian(y, x, i=0, j=1) + dV_S = dde.grad.jacobian(y, x, i=0, j=0) + dV_SS = dde.grad.hessian(y, x, i=0, j=0) + return dV_tau - (0.5 * sigma**2 * S**2 * dV_SS + r * S * dV_S - r * y) + +For the boundary conditions, we define the payoff function at maturity (:math:`\tau=0`) and the Dirichlet conditions at :math:`S=0` and :math:`S=S_{max}`: + +.. code-block:: python + + ic = dde.icbc.IC( + geomtime, + lambda x: np.maximum(x[:, 0:1] - K, 0.0), + lambda _, on_initial: on_initial, + ) + bc_lower = dde.icbc.DirichletBC( + geomtime, + lambda x: np.zeros((len(x), 1)), + lambda x, on_boundary: on_boundary and np.isclose(x[0], 0.0), + ) + bc_upper = dde.icbc.DirichletBC( + geomtime, + lambda x: x[:, 0:1] - K * np.exp(-r * x[:, 1:2]), + lambda x, on_boundary: on_boundary and np.isclose(x[0], S_max), + ) + +**Sampling Strategy**: +To improve domain coverage and training stability, we use the **Sobol sequence** +for sampling points. Sobol sequences are low-discrepancy quasi-random sequences +that provide more uniform coverage of the domain than pseudorandom sampling [2]_. +According to the SciPy documentation, Sobol sequences lose their balance +properties when the sample size is not a power of 2 [3]_. Therefore, we choose +sample counts close to powers of two in this implementation: + +.. code-block:: python + + data = dde.data.TimePDE( + geomtime, + pde, + [ic, bc_lower, bc_upper], + num_domain=2047, + num_boundary=63, + num_initial=127, + num_test=2047, + train_distribution="Sobol", + ) + +**Network Architecture**: +We use a fully connected neural network with **4 hidden layers** and **28 neurons** per layer, using ``tanh`` activation and ``Glorot normal`` initialization: + +.. code-block:: python + + net = dde.nn.FNN([2] + [28] * 4 + [1], "tanh", "Glorot normal") + +**Training Strategy**: +We employ a two-stage training process: first using the **Adam** optimizer for 20,000 iterations, followed by **L-BFGS** to fine-tune the solution and achieve lower residual error. + +.. code-block:: python + + model = dde.Model(data, net) + model.compile("adam", lr=1e-3) + model.train(iterations=20000) + model.compile("L-BFGS") + losshistory, train_state = model.train() + +Validation +---------- + +The model is validated against the analytical Black-Scholes formula. In our +runs, we observed an L2 relative error around ``1e-3`` and a mean PDE residual +around ``0.05``. These hyperparameters (network width, depth, learning rate, +and iteration count) were selected empirically based on observed convergence +and were not subject to a systematic hyperparameter search. Further tuning +may improve accuracy. + +References +---------- + +.. [1] Tanios, R. (2021). "Physics Informed Neural Networks in Computational Finance: + High Dimensional Forward & Inverse Option Pricing". ETH Zurich Master's Thesis. + https://doi.org/10.3929/ethz-b-000491555 + +.. [2] Sobol', I. M. (1967). "On the distribution of points in a cube and the + approximate evaluation of integrals." USSR Computational Mathematics and + Mathematical Physics, 7(4), 86–112. + https://www.sciencedirect.com/science/article/abs/pii/0041555367901449 + +.. [3] SciPy documentation: ``scipy.stats.qmc.Sobol``. + https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.qmc.Sobol.html + + +Complete code +------------- + +.. literalinclude:: ../../../examples/pinn_forward/black_scholes_call.py + :language: python + diff --git a/examples/pinn_forward/black_scholes_call.py b/examples/pinn_forward/black_scholes_call.py new file mode 100644 index 000000000..bcd146d58 --- /dev/null +++ b/examples/pinn_forward/black_scholes_call.py @@ -0,0 +1,128 @@ +"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch, paddle + +Black-Scholes PDE for European call option pricing using PINN. + +This example solves the Black-Scholes PDE in time-to-maturity form: + dV/dtau = 0.5*sigma^2*S^2*d^2V/dS^2 + r*S*dV/dS - r*V + +where: + S : underlying asset price (x[:, 0]) + tau : time to maturity (x[:, 1]) + V : option price (network output y) + +with boundary conditions: + V(S, 0) = max(S - K, 0) [Initial condition at maturity] + V(0, tau) = 0 [Lower boundary] + V(S_max, tau) = S_max - K*exp(-r*tau) [Upper boundary] + +Reference: + Tanios, R. (2021). "Physics Informed Neural Networks in Computational Finance: + High Dimensional Forward & Inverse Option Pricing". ETH Zurich Master's Thesis. + https://doi.org/10.3929/ethz-b-000491555 +""" + +import deepxde as dde +import numpy as np +from scipy.stats import norm + + +# Problem parameters +K = 50.0 # Strike price +r = 0.05 # Risk-free rate +sigma = 0.2 # Volatility +T = 1.0 # Time to maturity +S_max = 150.0 # Maximum stock price + + +def pde(x, y): + """Black-Scholes PDE residual. + + Args: + x: Input tensor. + x[:, 0]: asset price S + x[:, 1]: time to maturity tau + y: Network output representing the option price V(S, tau). + + Returns: + PDE residual: dV/dtau - (0.5*sigma^2*S^2*d^2V/dS^2 + r*S*dV/dS - r*V) + """ + S = x[:, 0:1] + dV_tau = dde.grad.jacobian(y, x, i=0, j=1) + dV_S = dde.grad.jacobian(y, x, i=0, j=0) + dV_SS = dde.grad.hessian(y, x, i=0, j=0) + return dV_tau - (0.5 * sigma**2 * S**2 * dV_SS + r * S * dV_S - r * y) + + +def black_scholes_call(S, tau): + """Analytical Black-Scholes formula.""" + tau = np.maximum(tau, 1e-8) + S = np.maximum(S, 1e-10) + d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * tau) / (sigma * np.sqrt(tau)) + d2 = d1 - sigma * np.sqrt(tau) + return S * norm.cdf(d1) - K * np.exp(-r * tau) * norm.cdf(d2) + + +def gen_testdata(): + """Generate test data with analytical solution.""" + S = np.linspace(0, S_max, 256).reshape(256, 1) + tau = np.linspace(0, T, 201).reshape(201, 1) + SS, TT = np.meshgrid(S, tau) + X = np.vstack((np.ravel(SS), np.ravel(TT))).T + y = np.array([black_scholes_call(X[i, 0], X[i, 1]) for i in range(len(X))]) + return X, y[:, None] + + +# Computational geometry +geom = dde.geometry.Interval(0, S_max) +timedomain = dde.geometry.TimeDomain(0, T) +geomtime = dde.geometry.GeometryXTime(geom, timedomain) + +# Boundary and initial conditions +ic = dde.icbc.IC( + geomtime, + lambda x: np.maximum(x[:, 0:1] - K, 0.0), + lambda _, on_initial: on_initial, +) +bc_lower = dde.icbc.DirichletBC( + geomtime, + lambda x: np.zeros((len(x), 1)), + lambda x, on_boundary: on_boundary and np.isclose(x[0], 0.0), +) +bc_upper = dde.icbc.DirichletBC( + geomtime, + lambda x: x[:, 0:1] - K * np.exp(-r * x[:, 1:2]), + lambda x, on_boundary: on_boundary and np.isclose(x[0], S_max), +) + +# Define PDE problem +data = dde.data.TimePDE( + geomtime, + pde, + [ic, bc_lower, bc_upper], + num_domain=2047, + num_boundary=63, + num_initial=127, + num_test=2047, + train_distribution="Sobol", +) + +# Neural network +net = dde.nn.FNN([2] + [28] * 4 + [1], "tanh", "Glorot normal") +model = dde.Model(data, net) + +# Training +model.compile("adam", lr=1e-3) +model.train(iterations=20000) +model.compile("L-BFGS") +losshistory, train_state = model.train() + +# Plot and save results +dde.saveplot(losshistory, train_state, issave=True, isplot=True) + +# Validation +X, y_true = gen_testdata() +y_pred = model.predict(X) +f = model.predict(X, operator=pde) +print("Mean residual:", np.mean(np.absolute(f))) +print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred)) +np.savetxt("test.dat", np.hstack((X, y_true, y_pred)))