American Put Option example, Extension of #2056

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/american.put
 
 - `Beltrami flow <https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/Beltrami_flow.py>`_
 - `Wave propagation with spatio-temporal multi-scale Fourier feature architecture <https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/wave_1d.py>`_

docs/demos/pinn_forward/american.put.rst

@@ -0,0 +1,205 @@
+===================
+American Put Option
+===================
+
+Financial Context: American vs. European Options
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Before diving into the math, it is helpful to understand the core difference: exercise flexibility.
+
+* European Put: Can only be exercised at the expiration date (:math:`T`).
+* American Put: Can be exercised at any time :math:`t \le T`.
+
+Because an American option gives the holder more rights (the right to "get out early" if it is profitable), it must always be worth at least as much as a European option. This extra value is called the Early Exercise Premium.
+
+Problem Setup
+-------------
+
+We solve the Black-Scholes partial differential equation for pricing an American put option.
+Unlike European options, American options can be exercised at any time before maturity. This leads to a free-boundary problem, which is formulated as a Linear Complementarity Problem (LCP).
+The problem is formulated in terms of time-to-maturity :math:`\tau = T - t`.
+
+The LCP is given by:
+
+.. math::
+
+    L(V) &= -\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 \le 0 \\
+    V(S, \tau) &\ge \max(K - S, 0) \\
+    (V(S, \tau) - \max(K - S, 0)) \cdot L(V) &= 0
+
+Financial Context: The "No-Arbitrage" Constraints (The LCP)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The three equations above satisfy the Linear Complementarity Problem. Financially, they represent the rules of a market without free money:
+
+* The Payoff Floor (:math:`V(S, \tau) \ge \max(K - S, 0)`): An American option can never be worth less than its intrinsic value (what you'd get if you exercised right now). If it were, an arbitrageur could buy the option and exercise it instantly for an immediate risk-free profit.
+* The Time-Value Decay (:math:`L(V) \le 0`): The Black-Scholes operator :math:`L(V)` represents the option's value decay. In the region where you don't exercise, the option follows the standard Black-Scholes PDE (:math:`L(V) = 0`). In the region where you do exercise early, the option's value is decaying faster than the risk-free rate (:math:`L(V) < 0`).
+* Complementarity Condition: At any given time, exactly one of these must be true: either the option behaves exactly like a PDE (holding is optimal), or the option is worth exactly its exercise value (exercising is optimal).
+
+where:
+
+* :math:`V(S, \tau)` is the option price
+* :math:`S` is the underlying asset price
+* :math:`\tau` is the time to maturity
+* :math:`K=1.0` is the strike price
+* :math:`\sigma=0.2` is the volatility
+* :math:`r=0.05` is the risk-free interest rate
+
+Financial Context: The Role of the Parameters
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+* Volatility (:math:`\sigma = 0.2`): This represents the fear or uncertainty in the market. Higher volatility usually makes the American option more valuable because there is a higher chance the stock price will swing deeply into the profit zone.
+* Risk-free Rate (:math:`r = 0.05`): This is the opportunity cost of money. For a Put option, a higher interest rate actually makes early exercise more attractive. When you exercise a put, you receive cash (:math:`K`). The sooner you get that cash, the sooner you can put it in the bank to earn interest (:math:`r`).
+
+The boundary and initial conditions are:
+
+1. Initial Condition (at maturity :math:`\tau=0`):
+
+.. math:: 
+
+    V(S, 0) = \max(K - S, 0)
+
+2. Boundary Conditions:
+
+.. math::
+
+    V(0, \tau) &= K \\
+    V(S_{max}, \tau) &= 0
+
+where :math:`S_{max} = 4.0` limits the domain to a reasonable multiple of the strike price.
+
+Implementation
+--------------
+
+This example uses a Physics-Informed Neural Network (PINN) to solve the American put LCP using the Fischer-Burmeister function.
+
+First, we import DeepXDE, PyTorch, and define the problem parameters:
+
+.. code-block:: python
+
+    import deepxde as dde
+    import numpy as np
+    import torch
+
+    K = 1.0       # Strike price
+    r = 0.05      # Risk-free rate
+    sigma = 0.2   # Volatility
+    T = 1.0       # Time to maturity
+    S_max = 4.0   # Maximum stock price
+
+We define the computational geometry. We use ``GeometryXTime`` to combine the asset price interval :math:`[0, 4.0]` and the time domain :math:`[0, 1.0]`:
+
+.. code-block:: python
+
+    geom = dde.geometry.Interval(0, S_max)
+    timedomain = dde.geometry.TimeDomain(0, T)
+    geomtime = dde.geometry.GeometryXTime(geom, timedomain)
+
+For the American put, the PDE residual must account for the early exercise constraint. This is very important to make sure the model does not learn European Options.
+We formulate this using soft penalties and the Fischer-Burmeister equation :math:`\phi(a, b) = \sqrt{a^2 + b^2} - (a + b) = 0` to satisfy the complementarity condition:
+
+.. code-block:: python
+
+    def pde(x, y):
+        s = x[:, 0:1]
+        v = y
+
+        dv_dtau = dde.grad.jacobian(y, x, i=0, j=1) 
+        dv_ds = dde.grad.jacobian(y, x, i=0, j=0)
+        dv_ds2 = dde.grad.hessian(y, x, i=0, j=0)
+
+        # Black-Scholes Operator
+        bs_op = -dv_dtau + 0.5 * (sigma**2) * (s**2) * dv_ds2 + r * s * dv_ds - r * v
+
+        # Intrinsic Value (Payoff)
+        payoff = dde.backend.relu(K - s)
+
+        # LCP Constraints
+        res_pde = dde.backend.relu(bs_op)            # Soft penalty for L(V) > 0
+        penalty = dde.backend.relu(payoff - v)       # Soft penalty for V < Payoff
+
+        # Fischer-Burmeister Complementarity
+        a = v - payoff
+        b = -bs_op 
+        # Add 1e-8 for numerical stability of the derivative of sqrt at 0
+        res_fb = torch.sqrt(a**2 + b**2 + 1e-8) - (a + b)
+
+        return [res_pde, penalty, res_fb]
+
+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
+
+    def initial_condition(x):
+        s = x[:, 0:1]
+        return np.maximum(K - s, 0)
+
+    IC = dde.icbc.IC(geomtime, initial_condition, lambda _, on_initial: on_initial)
+
+    LB = dde.icbc.DirichletBC(
+        geomtime, 
+        lambda x: K, 
+        lambda x, on_boundary: on_boundary and np.isclose(x[0], 0), 
+        component=0
+    )
+    HB = dde.icbc.DirichletBC(
+        geomtime, 
+        lambda x: 0, 
+        lambda x, on_boundary: on_boundary and np.isclose(x[0], S_max), 
+        component=0
+    )
+
+Sampling Strategy:
+~~~~~~~~~~~~~~~~~~
+To perfectly capture the non-differentiable kink at the strike price, we densely sample anchor points around :math:`S=1.0`. These are combined with the standard domain and boundary points:
+
+.. code-block:: python
+
+    S_anchors = np.random.uniform(0.5 * K, 1.5 * K, (1000, 1))
+    tau_anchors = np.random.uniform(0, T, (1000, 1))
+    anchors = np.hstack((S_anchors, tau_anchors))
+
+    data = dde.data.TimePDE(
+        geomtime, 
+        pde, 
+        [IC, LB, HB],
+        num_domain=6000, 
+        num_boundary=1500, 
+        num_initial=1500,
+        anchors=anchors
+    )
+
+Network Architecture:
+~~~~~~~~~~~~~~~~~~~~~
+We use a fully connected neural network with 4 hidden layers and 128 neurons per layer, increasing the width compared to the European option to handle the sharper Fischer-Burmeister gradients:
+
+.. code-block:: python
+
+    net = dde.nn.FNN([2] + [128] * 4 + [1], "tanh", "Glorot normal")
+
+Training Strategy:
+~~~~~~~~~~~~~~~~~~
+We employ a two-stage training process. First, we use the Adam optimizer for 12,000 iterations, heavily weighting the boundary conditions and the intrinsic payoff penalty. Then, we use L-BFGS with a ``PDEPointResampler`` to adaptively refine the solution.
+
+.. code-block:: python
+
+    model = dde.Model(data, net)
+
+    # Heavily weight the initial/boundary conditions and the payoff penalty
+    model.compile("adam", lr=0.001, loss_weights=[1, 10, 5, 100, 100, 100])
+    model.train(iterations=12000)
+
+    model.compile("L-BFGS")
+    checker = dde.callbacks.PDEPointResampler(period=1000)
+    losshistory, train_state = model.train(callbacks=[checker])
+
+Validation
+----------
+
+The model is validated against a standard Binomial Tree pricing model. Evaluating the results alongside the analytical binomial ground truth allows us to compute the relative L2 error and visually confirm the model captures the early exercise premium.
+
+Complete code
+-------------
+
+.. literalinclude:: ../../../examples/pinn_forward/american_put.py
+   :language: python

examples/pinn_forward/american_put.py

@@ -0,0 +1,150 @@
+"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch, paddle
+
+Black-Scholes PDE for American Put Option:
+∂V/∂t + 0.5 * σ^2 * S^2 * ∂^2V/∂S^2 + r * S * ∂V/∂S - r * V = 0, for S > 0, t < T
+With the early exercise constraint:
+V(S, t) >= max(K - S, 0) for all S, t
+And the complementarity condition:
+(V(S, t) - max(K - S, 0)) * (∂V/∂t + 0.5 * σ^2 * S^2 * ∂^2V/∂S^2 + r * S * ∂V/∂S - r * V) = 0   
+"""
+import deepxde as dde
+import numpy as np
+import matplotlib.pyplot as plt
+import math 
+
+# --- 1. Financial Parameters ---
+r = 0.05
+sigma = 0.2
+K = 1.0 
+# Anything too high above strike price will have zero payoff, so we can limit the domain to a few multiples of K
+S_max = K * 4.0 
+T = 1.0
+
+# --- 2. PDE Definition (Fischer-Burmeister LCP) ---
+def pde(x, y):
+    s = x[:, 0:1]
+    v = y
+
+    dv_dtau = dde.grad.jacobian(y, x, i=0, j=1) 
+    dv_ds = dde.grad.jacobian(y, x, i=0, j=0)
+    dv_ds2 = dde.grad.hessian(y, x, i=0, j=0)
+
+    # Black-Scholes Operator (Time reversed: tau = T - t)
+    bs_op = -dv_dtau + 0.5 * (sigma**2) * (s**2) * dv_ds2 + r * s * dv_ds - r * v
+
+    # Intrinsic Value (Payoff)
+    payoff = dde.backend.relu(K - s)
+
+    # LCP Constraints
+    # 1. Soft penalty for L(V) > 0
+    res_pde = dde.backend.relu(bs_op) 
+
+    # 2. Soft penalty for V < Payoff
+    penalty = dde.backend.relu(payoff - v)
+
+    # 3. Fischer-Burmeister Complementarity
+    # phi(a, b) = sqrt(a^2 + b^2) - (a + b) = 0
+    a = v - payoff
+    b = -bs_op 
+    # Add 1e-8 for numerical stability of the derivative of sqrt at 0
+    res_fb = dde.backend.pow(a**2 + b**2 + 1e-8, 0.5) - (a + b)
+
+    return [res_pde, penalty, res_fb]
+
+# --- 3. Geometry and Time Domain ---
+geom = dde.geometry.Interval(0, S_max)
+timedomain = dde.geometry.TimeDomain(0, T)
+geomtime = dde.geometry.GeometryXTime(geom, timedomain)
+
+# Create anchor points dense around the strike price (S=1.0) 
+# This forces the network to learn the kink perfectly.
+S_anchors = np.random.uniform(0.5 * K, 1.5 * K, (1000, 1))
+tau_anchors = np.random.uniform(0, T, (1000, 1))
+anchors = np.hstack((S_anchors, tau_anchors))
+
+# --- 4. Boundary and Initial Conditions ---
+def initial_condition(x):
+    s = x[:, 0:1]
+    return np.maximum(K - s, 0)
+
+IC = dde.icbc.IC(geomtime, initial_condition, lambda _, on_initial: on_initial)
+
+def low_boundary(x, on_boundary):
+    return on_boundary and np.isclose(x[0], 0)
+
+def high_boundary(x, on_boundary):
+    return on_boundary and np.isclose(x[0], S_max)
+
+LB = dde.icbc.DirichletBC(geomtime, lambda x: K, low_boundary, component=0)
+HB = dde.icbc.DirichletBC(geomtime, lambda x: 0, high_boundary, component=0)
+
+# Compile Data
+data = dde.data.TimePDE(
+    geomtime, 
+    pde, 
+    [IC, LB, HB],
+    num_domain=6000, 
+    num_boundary=1500, 
+    num_initial=1500,
+    anchors=anchors # Pass the focus points
+)
+
+# --- 5. Network Architecture ---
+# Increased width slightly to handle the sharper FB gradients
+net = dde.nn.FNN([2] + [128] * 4 + [1], "tanh", "Glorot normal")
+model = dde.Model(data, net)
+
+# --- 6. Training ---
+print("Stage 1: Training with Adam...")
+# Weights: [res_pde, penalty, res_fb, IC, LB, HB]
+# Heavily weight the initial/boundary conditions and the payoff penalty
+model.compile("adam", lr=0.001, loss_weights=[1, 10, 5, 100, 100, 100])
+model.train(iterations=12000)
+
+print("Stage 2: Fine-tuning with L-BFGS and Adaptive Refinement...")
+model.compile("L-BFGS")
+checker = dde.callbacks.PDEPointResampler(period=1000)
+losshistory, train_state = model.train(callbacks=[checker])
+
+dde.saveplot(losshistory, train_state, issave=True, isplot=True)
+
+# --- 7. Ground Truth Evaluation ---
+def binomial_american_put(S, K, T, r, sigma, N=2000):
+    dt = T / N
+    u = np.exp(sigma * np.sqrt(dt))
+    d = 1 / u
+    p = (np.exp(r * dt) - d) / (u - d)
+
+    S_tree = S * d**(np.arange(N, -1, -1)) * u**(np.arange(0, N + 1, 1))
+    V = np.maximum(K - S_tree, 0)
+
+    for i in range(N - 1, -1, -1):
+        V = (p * V[1:] + (1 - p) * V[:-1]) * np.exp(-r * dt)
+        S_i = S * d**(np.arange(i, -1, -1)) * u**(np.arange(0, i + 1, 1))
+        V = np.maximum(V, K - S_i) 
+
+    return V[0]
+
+S_eval = np.linspace(0.1, 2.0, 100).reshape(-1, 1) 
+tau_eval = np.ones_like(S_eval) * T
+X_test = np.hstack((S_eval, tau_eval))
+
+y_pred = model.predict(X_test)
+
+print("Calculating Binomial Ground Truth (this may take a moment)...")
+y_true = [binomial_american_put(s[0], K, T, r, sigma) for s in S_eval]
+
+plt.figure(figsize=(10, 6))
+plt.plot(S_eval, y_true, 'r-', label='Ground Truth (Binomial Tree)', linewidth=2)
+plt.plot(S_eval, y_pred, 'b--', label='PINN Prediction', linewidth=2)
+plt.axvline(K, color='black', linestyle=':', label='Strike Price (K)')
+plt.xlabel('Stock Price (S)')
+plt.ylabel('Option Value (V)')
+plt.title('American Put Option: PINN vs. Ground Truth')
+plt.legend()
+plt.grid(True, alpha=0.3)
+plt.savefig("american_put_comparison.png")
+plt.show()
+
+error_l2 = np.linalg.norm(y_true - y_pred.flatten()) / np.linalg.norm(y_true)
+print(f"Relative L2 Error: {error_l2:.4e}")