Stops and targets as barrier problems

"""M04A exercise 1: stops and targets as barrier problems (notes §4A.6).

Model: the trade's P&L per unit, measured from entry, is an arithmetic Brownian
motion X_t = mu*t + sigma*W_t, X_0 = 0 (t in days, X in rupees or index points).
A bracket order exits at the target +a or the stop -b (a, b > 0), whichever is
hit first; tau is that exit time.

Closed forms you will implement (Propositions 4A.6-4A.9):
  P(target first) = b/(a+b)                                          (mu = 0)
                  = (1 - exp(g*b)) / (exp(-g*a) - exp(g*b)),  g = 2mu/sigma^2
  E[tau]          = a*b/sigma^2                                       (mu = 0)
                  = (a*P - b*(1-P)) / mu
  E[X_tau]        = mu * E[tau]                                (optional stopping)
  P(min_{t<=T} X_t <= -b) = Phi((-b - mu T)/(s sqrt T)) + exp(-2 mu b/s^2) Phi((-b + mu T)/(s sqrt T))

Discrete monitoring every dt: a barrier monitored only on the grid behaves like a
continuously monitored barrier shifted away from the start by BETA*sigma*sqrt(dt),
BETA = -zeta(1/2)/sqrt(2*pi) ~ 0.5826 (Broadie, Glasserman and Kou, 1997).

Run (tests import this file; set COURSE_IMPL=solutions to run the reference):
    cd exercises && .venv/bin/python -m pytest modules/M04A-reading-markets-trade-management/tests/test_m04a_barriers.py
"""

from __future__ import annotations

import math

import numpy as np
from scipy.stats import norm

BETA = 0.5825971579390106  # -zeta(1/2)/sqrt(2 pi)


def two_barrier_hit_prob(a: float, b: float, mu: float, sigma: float) -> float:
    """P(X hits +a before -b), X_t = mu t + sigma W_t from 0, continuous monitoring.

    Raise ValueError unless a, b, sigma > 0. For |2 mu (a+b)/sigma^2| tiny, return
    b/(a+b) (use math.expm1 elsewhere to avoid cancellation).
    """
    raise NotImplementedError


def expected_exit_time(a: float, b: float, mu: float, sigma: float) -> float:
    """E[tau] for the two-barrier exit: a b / sigma^2 if mu = 0, else (a P - b (1-P)) / mu."""
    raise NotImplementedError


def bracket_expectancy(a: float, b: float, mu: float, sigma: float) -> float:
    """E[X_tau] = a P - b (1-P); it must equal mu * E[tau]."""
    raise NotImplementedError


def stop_hit_prob(b: float, T: float, mu: float, sigma: float) -> float:
    """P(min_{t<=T} X_t <= -b) for a single stop at -b over horizon T (reflection principle).

    Return 0.0 for T <= 0.
    """
    raise NotImplementedError


def discrete_barrier_shift(sigma: float, dt: float) -> float:
    """BETA * sigma * sqrt(dt)."""
    raise NotImplementedError


def simulate_two_barriers(
    a: float, b: float, mu: float, sigma: float, dt: float, n_paths: int, seed: int, max_steps: int = 100_000
) -> dict:
    """Monte Carlo of the bracket monitored on a grid of step dt.

    Use np.random.default_rng(seed). Exit at the first grid time k*dt with X >= a
    (target) or X <= -b (stop); the exit P&L is the observed X then (it overshoots
    the barrier). Paths still open after max_steps are closed at their value.

    Returns dict(p_target, mean_tau, mean_pnl, se_pnl) with se_pnl the standard
    error of mean_pnl (sample std with ddof=1 over sqrt(n_paths)).
    """
    raise NotImplementedError


def simulate_stop_hit(b: float, T: float, mu: float, sigma: float, n_steps: int, n_paths: int, seed: int) -> float:
    """Fraction of grid-monitored paths (n_steps equal steps over [0, T], X_0 = 0 not
    counted) whose minimum over the grid is <= -b. Use np.random.default_rng(seed)."""
    raise NotImplementedError