core.py

"""
Core functions for Nested Multilevel Monte Carlo Simulation.

Based on the Master Thesis by Nicolai Ehrhardt, TU Kaiserslautern, 2018.
Supervised by Prof. Dr. Klaus Ritter, AG Computational Stochastics.
"""

import numpy as np
from scipy.stats import norm


def sim_brownian_motion(T, l, M):
    """Simulate M independent Brownian motion paths on [0, T] with 2^l steps.

    Parameters
    ----------
    T : float
        Terminal time.
    l : int
        Refinement level (number of steps = 2^l).
    M : int
        Number of independent paths.

    Returns
    -------
    bm : ndarray, shape (2^l + 1, M)
        Brownian motion paths, starting at 0.
    """
    n_steps = 2**l
    dt = T / n_steps
    increments = np.random.normal(0, np.sqrt(dt), size=(n_steps, M))
    bm = np.zeros((n_steps + 1, M))
    bm[1:, :] = np.cumsum(increments, axis=0)
    return bm


def euler_gbm(s0, mu, sigma, T, l, bm):
    """Euler-Maruyama approximation of geometric Brownian motion.

    Parameters
    ----------
    s0 : float
        Initial value.
    mu : float
        Drift coefficient.
    sigma : float
        Volatility.
    T : float
        Terminal time.
    l : int
        Refinement level (2^l time steps).
    bm : ndarray, shape (2^l + 1, M)
        Brownian motion paths.

    Returns
    -------
    gbm : ndarray, shape (2^l + 1, M)
        Euler approximation paths.
    """
    n_steps = 2**l
    M = bm.shape[1] if bm.ndim > 1 else 1
    if bm.ndim == 1:
        bm = bm[:, np.newaxis]
    dt = T / n_steps
    gbm = np.zeros_like(bm)
    gbm[0, :] = s0
    for k in range(1, n_steps + 1):
        dW = bm[k, :] - bm[k - 1, :]
        gbm[k, :] = gbm[k - 1, :] + mu * gbm[k - 1, :] * dt + sigma * gbm[k - 1, :] * dW
    return gbm


def milstein_gbm(s0, mu, sigma, T, l, bm):
    """Milstein approximation of geometric Brownian motion.

    Parameters
    ----------
    s0 : float
        Initial value.
    mu : float
        Drift coefficient.
    sigma : float
        Volatility.
    T : float
        Terminal time.
    l : int
        Refinement level (2^l time steps).
    bm : ndarray, shape (2^l + 1, M)
        Brownian motion paths.

    Returns
    -------
    gbm : ndarray, shape (2^l + 1, M)
        Milstein approximation paths.
    """
    n_steps = 2**l
    if bm.ndim == 1:
        bm = bm[:, np.newaxis]
    dt = T / n_steps
    gbm = np.zeros_like(bm)
    gbm[0, :] = s0
    for k in range(1, n_steps + 1):
        dW = bm[k, :] - bm[k - 1, :]
        gbm[k, :] = (gbm[k - 1, :]
                      + mu * gbm[k - 1, :] * dt
                      + sigma * gbm[k - 1, :] * dW
                      + 0.5 * sigma**2 * gbm[k - 1, :] * (dW**2 - dt))
    return gbm


def exact_gbm(s0, mu, sigma, t, bm):
    """Exact geometric Brownian motion using the analytic formula.

    S(t) = s0 * exp((mu - sigma^2/2) * t + sigma * W(t))

    Parameters
    ----------
    s0 : float
        Initial value.
    mu : float
        Drift coefficient.
    sigma : float
        Volatility.
    t : ndarray, shape (N,)
        Time points.
    bm : ndarray, shape (N,) or (N, M)
        Brownian motion values at time points t.

    Returns
    -------
    gbm : ndarray
        Exact GBM paths.
    """
    if bm.ndim == 1:
        return s0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * bm)
    else:
        return s0 * np.exp((mu - 0.5 * sigma**2) * t[:, np.newaxis] + sigma * bm)


def construct_smoothing_polynomial(r):
    """Construct the smoothing polynomial p of degree r+1.

    The polynomial satisfies the moment conditions (S4) from Chapter 4.1:
    - p(1) = 0, p(-1) = 1
    - integral conditions for moments j = 0, ..., r-1

    Parameters
    ----------
    r : int
        Smoothness parameter (density is r-times continuously differentiable).

    Returns
    -------
    coeffs : ndarray
        Polynomial coefficients in descending order (for use with np.polyval).
    """
    n = r + 2
    A = np.zeros((n, n))
    b = np.zeros(n)

    # p(1) = sum(a_k) = 0
    A[0, :] = 1.0
    b[0] = 0.0

    # p(-1) = sum(a_k * (-1)^k) = 1
    for k in range(n):
        A[1, k] = (-1)**k
    b[1] = 1.0

    # Moment conditions: integral from -1 to 1 of s^j * (1_{s<=0} - p(s)) ds = 0
    # for j = 0, ..., r-1
    for j in range(r):
        b[2 + j] = (-1)**j / (2 * (j + 1))
        for k in range(n):
            if (k + j + 1) % 2 == 1:
                A[2 + j, k] = 1.0 / (k + j + 1)
            else:
                A[2 + j, k] = 0.0

    coeffs = np.linalg.solve(A, b)
    # Flip to get descending order for np.polyval
    return coeffs[::-1]


def smoothing_function(t, s, delta, coeffs):
    """Evaluate the smoothing function Psi((t-s)/delta).

    Psi approximates the indicator function 1_{(-inf, s]}.
    - Psi(x) = 1 for x < -1
    - Psi(x) = p(x) for |x| <= 1
    - Psi(x) = 0 for x > 1

    Parameters
    ----------
    t : float or ndarray
        Points to evaluate.
    s : float or ndarray
        Threshold points.
    delta : float
        Smoothing parameter.
    coeffs : ndarray
        Polynomial coefficients (descending order).

    Returns
    -------
    result : ndarray
        Values of the smoothing function.
    """
    t = np.atleast_1d(np.asarray(t, dtype=float))
    s = np.atleast_1d(np.asarray(s, dtype=float))

    # Handle broadcasting for (t, s) pairs
    if t.shape != s.shape:
        t, s = np.broadcast_arrays(t, s)

    x = (t - s) / delta
    result = np.where(x < -1, 1.0, np.where(x > 1, 0.0, np.polyval(coeffs, x)))
    return result


def european_call_payoff(S_T, K):
    """European call option payoff: max(S_T - K, 0)."""
    return np.maximum(S_T - K, 0)


def lookback_call_payoff(path, K=0):
    """Lookback call option payoff: max(max(path) - K, 0)."""
    return np.maximum(np.max(path, axis=0) - K, 0)


def black_scholes_call(s0, K, T, mu, sigma):
    """Black-Scholes formula for European call option price (undiscounted).

    Parameters
    ----------
    s0 : float or ndarray
        Initial stock price.
    K : float
        Strike price.
    T : float
        Time to maturity.
    mu : float
        Drift (risk-neutral rate).
    sigma : float
        Volatility.

    Returns
    -------
    price : float or ndarray
        Option price.
    """
    d1 = (np.log(s0 / K) + (mu + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return s0 * norm.cdf(d1) - np.exp(-mu * T) * K * norm.cdf(d2)


def mlmc_european_call(x, mu, sigma, T0, T1, L, K):
    """Inner MLMC estimator for European call option price.

    Computes Z^L(x) = sum_{l=0}^{L} Z_l, where Z_l uses M_l = ceil(2^{L-l}/L)
    replications with coupled Euler fine/coarse simulations.

    Parameters
    ----------
    x : float
        Initial value for the inner simulation (= risk scenario at T0).
    mu : float
        Drift.
    sigma : float
        Volatility.
    T0 : float
        Start time of inner simulation.
    T1 : float
        Maturity.
    L : int
        Maximum level.
    K : float
        Strike price.

    Returns
    -------
    Z : float
        MLMC estimate of E[Phi(S(T1)) | S(T0) = x].
    """
    Z_levels = np.zeros(L + 1)
    dt_inner = T1 - T0

    for l in range(L + 1):
        M_l = max(1, int(np.ceil(2**(L - l) / L)))
        n_fine = 2**l
        dt_fine = dt_inner / n_fine

        bm = sim_brownian_motion(dt_inner, l, M_l)

        # Fine Euler simulation
        gbm_fine = np.zeros((n_fine + 1, M_l))
        gbm_fine[0, :] = x
        for k in range(1, n_fine + 1):
            dW = bm[k, :] - bm[k - 1, :]
            gbm_fine[k, :] = (gbm_fine[k - 1, :]
                              + mu * gbm_fine[k - 1, :] * dt_fine
                              + sigma * gbm_fine[k - 1, :] * dW)

        Phi_fine = european_call_payoff(gbm_fine[-1, :], K)

        if l > 0:
            # Coarse Euler simulation using same Brownian motion
            n_coarse = 2**(l - 1)
            dt_coarse = dt_inner / n_coarse
            bm_coarse = bm[::2, :]  # subsample

            gbm_coarse = np.zeros((n_coarse + 1, M_l))
            gbm_coarse[0, :] = x
            for k in range(1, n_coarse + 1):
                dW_c = bm_coarse[k, :] - bm_coarse[k - 1, :]
                gbm_coarse[k, :] = (gbm_coarse[k - 1, :]
                                    + mu * gbm_coarse[k - 1, :] * dt_coarse
                                    + sigma * gbm_coarse[k - 1, :] * dW_c)

            Phi_coarse = european_call_payoff(gbm_coarse[-1, :], K)
        else:
            Phi_coarse = np.zeros(M_l)

        Z_levels[l] = np.mean(Phi_fine - Phi_coarse)

    return np.sum(Z_levels)


def mlmc_milstein_european_call(x, mu, sigma, T0, T1, L, K):
    """Inner MLMC estimator using Milstein method for European call.

    Uses M_l = ceil(2^{L-3/2*l}) replications (Theorem 5.14).

    Parameters
    ----------
    x : float
        Initial value.
    mu, sigma : float
        SDE coefficients.
    T0, T1 : float
        Time interval.
    L : int
        Maximum level.
    K : float
        Strike price.

    Returns
    -------
    Z : float
        MLMC Milstein estimate.
    """
    Z_levels = np.zeros(L + 1)
    dt_inner = T1 - T0

    for l in range(L + 1):
        M_l = max(1, int(np.ceil(2**(L - 1.5 * l))))
        n_fine = 2**l
        dt_fine = dt_inner / n_fine

        bm = sim_brownian_motion(dt_inner, l, M_l)

        # Fine Milstein simulation
        gbm_fine = np.zeros((n_fine + 1, M_l))
        gbm_fine[0, :] = x
        for k in range(1, n_fine + 1):
            dW = bm[k, :] - bm[k - 1, :]
            gbm_fine[k, :] = (gbm_fine[k - 1, :]
                              + mu * gbm_fine[k - 1, :] * dt_fine
                              + sigma * gbm_fine[k - 1, :] * dW
                              + 0.5 * sigma**2 * gbm_fine[k - 1, :] * (dW**2 - dt_fine))

        Phi_fine = european_call_payoff(gbm_fine[-1, :], K)

        if l > 0:
            n_coarse = 2**(l - 1)
            dt_coarse = dt_inner / n_coarse
            bm_coarse = bm[::2, :]

            gbm_coarse = np.zeros((n_coarse + 1, M_l))
            gbm_coarse[0, :] = x
            for k in range(1, n_coarse + 1):
                dW_c = bm_coarse[k, :] - bm_coarse[k - 1, :]
                gbm_coarse[k, :] = (gbm_coarse[k - 1, :]
                                    + mu * gbm_coarse[k - 1, :] * dt_coarse
                                    + sigma * gbm_coarse[k - 1, :] * dW_c
                                    + 0.5 * sigma**2 * gbm_coarse[k - 1, :] * (dW_c**2 - dt_coarse))

            Phi_coarse = european_call_payoff(gbm_coarse[-1, :], K)
        else:
            Phi_coarse = np.zeros(M_l)

        Z_levels[l] = np.mean(Phi_fine - Phi_coarse)

    return np.sum(Z_levels)