Touched up notebooks + webapp
This commit is contained in:
@@ -0,0 +1,37 @@
|
||||
"""Public API for the Adaptive Barrier Monitor engine."""
|
||||
|
||||
__version__ = "0.4.7"
|
||||
|
||||
from adaptive_barrier.engine import (
|
||||
DROP_FRACTION,
|
||||
LOG_BARRIER,
|
||||
LOG_UPPER_BARRIER,
|
||||
RISE_FRACTION,
|
||||
TRADING_DAYS_PER_YEAR,
|
||||
TRADING_MINUTES_PER_YEAR,
|
||||
TYPICAL_ANNUAL_VOL,
|
||||
WINDOW_MINUTES,
|
||||
WINDOW_YEARS,
|
||||
adaptive_schedule,
|
||||
barrier_miss_prob,
|
||||
brownian_bridge,
|
||||
brownian_motion,
|
||||
brownian_motion_cholesky,
|
||||
detect_barrier_event,
|
||||
estimate_bridge_breach_prob,
|
||||
estimate_first_passage_prob,
|
||||
estimate_first_passage_prob_bb,
|
||||
first_passage_cdf,
|
||||
first_passage_cdf_zero_drift,
|
||||
fixed_cadence_indices,
|
||||
fixed_uniform_indices,
|
||||
geometric_brownian_motion,
|
||||
horizon_vol,
|
||||
max_safe_dt,
|
||||
merton_jump_diffusion,
|
||||
minutes_to_years,
|
||||
run_monte_carlo_simulation,
|
||||
time_grid,
|
||||
)
|
||||
|
||||
__all__ = [name for name in globals() if not name.startswith("_")]
|
||||
@@ -0,0 +1,781 @@
|
||||
"""Core stochastic-process and barrier-monitoring utilities.
|
||||
|
||||
The closed-form Brownian-motion results are exact under their stated diffusion
|
||||
assumptions. The adaptive schedule is deliberately described as a *local
|
||||
heuristic*: before the next observation is known, it substitutes the current
|
||||
barrier distance for both Brownian-bridge endpoint distances. Consequently,
|
||||
``eps`` is a design parameter, not an unconditional real-time miss guarantee,
|
||||
and it does not control jump risk.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from typing import Optional
|
||||
|
||||
import numpy as np
|
||||
from scipy.stats import norm
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Trading-time and volatility conventions
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
TRADING_DAYS_PER_YEAR = 252
|
||||
TRADING_HOURS_PER_DAY = 6.5
|
||||
TRADING_MINUTES_PER_HOUR = 60
|
||||
TRADING_MINUTES_PER_DAY = TRADING_HOURS_PER_DAY * TRADING_MINUTES_PER_HOUR
|
||||
TRADING_MINUTES_PER_YEAR = TRADING_DAYS_PER_YEAR * TRADING_MINUTES_PER_DAY
|
||||
|
||||
DROP_FRACTION = 0.10
|
||||
RISE_FRACTION = 0.10
|
||||
WINDOW_MINUTES = 5
|
||||
WINDOW_YEARS = WINDOW_MINUTES / TRADING_MINUTES_PER_YEAR
|
||||
LOG_BARRIER = float(np.log1p(-DROP_FRACTION))
|
||||
LOG_UPPER_BARRIER = float(np.log1p(RISE_FRACTION))
|
||||
TYPICAL_ANNUAL_VOL = 0.30
|
||||
|
||||
|
||||
def _require_positive(name: str, value: float) -> None:
|
||||
if not np.isfinite(value) or value <= 0:
|
||||
raise ValueError(f"{name} must be finite and positive")
|
||||
|
||||
|
||||
def _require_positive_int(name: str, value: int) -> None:
|
||||
if isinstance(value, bool) or int(value) != value or value <= 0:
|
||||
raise ValueError(f"{name} must be a positive integer")
|
||||
|
||||
|
||||
def minutes_to_years(minutes: float | np.ndarray) -> float | np.ndarray:
|
||||
"""Convert trading minutes to years."""
|
||||
values = np.asarray(minutes, dtype=float)
|
||||
if np.any(~np.isfinite(values)) or np.any(values < 0):
|
||||
raise ValueError("minutes must be finite and non-negative")
|
||||
result = values / TRADING_MINUTES_PER_YEAR
|
||||
return float(result) if result.ndim == 0 else result
|
||||
|
||||
|
||||
def horizon_vol(sigma_annual: float, minutes: float | np.ndarray) -> float | np.ndarray:
|
||||
"""Return the diffusion standard deviation over ``minutes`` of trading time."""
|
||||
_require_positive("sigma_annual", sigma_annual)
|
||||
result = sigma_annual * np.sqrt(minutes_to_years(minutes))
|
||||
return float(result) if np.ndim(result) == 0 else result
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Process samplers
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
|
||||
def _as_rng(rng: Optional[np.random.Generator] = None) -> np.random.Generator:
|
||||
return np.random.default_rng() if rng is None else rng
|
||||
|
||||
|
||||
def time_grid(T: float, n_steps: int) -> np.ndarray:
|
||||
"""Uniform grid from 0 to ``T`` with ``n_steps + 1`` points."""
|
||||
_require_positive("T", T)
|
||||
_require_positive_int("n_steps", n_steps)
|
||||
return np.linspace(0.0, T, int(n_steps) + 1)
|
||||
|
||||
|
||||
def brownian_motion(
|
||||
T: float,
|
||||
n_steps: int,
|
||||
n_paths: int = 1,
|
||||
drift: float = 0.0,
|
||||
sigma: float = 1.0,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Sample ``dX = drift dt + sigma dW`` using independent Gaussian increments."""
|
||||
_require_positive("sigma", sigma)
|
||||
_require_positive_int("n_paths", n_paths)
|
||||
if not np.isfinite(drift):
|
||||
raise ValueError("drift must be finite")
|
||||
rng = _as_rng(rng)
|
||||
t = time_grid(T, n_steps)
|
||||
dt = T / n_steps
|
||||
increments = rng.normal(
|
||||
loc=drift * dt,
|
||||
scale=sigma * np.sqrt(dt),
|
||||
size=(int(n_paths), int(n_steps)),
|
||||
)
|
||||
paths = np.zeros((int(n_paths), int(n_steps) + 1))
|
||||
paths[:, 1:] = np.cumsum(increments, axis=1)
|
||||
return t, paths
|
||||
|
||||
|
||||
def brownian_motion_cholesky(
|
||||
T: float,
|
||||
n_steps: int,
|
||||
n_paths: int = 1,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Sample standard Brownian motion from the ``min(s,t)`` covariance matrix."""
|
||||
_require_positive_int("n_paths", n_paths)
|
||||
rng = _as_rng(rng)
|
||||
t = time_grid(T, n_steps)
|
||||
inner = t[1:]
|
||||
covariance = np.minimum(inner[:, None], inner[None, :])
|
||||
factor = np.linalg.cholesky(covariance)
|
||||
normals = rng.standard_normal(size=(int(n_paths), int(n_steps)))
|
||||
paths = np.zeros((int(n_paths), int(n_steps) + 1))
|
||||
paths[:, 1:] = normals @ factor.T
|
||||
return t, paths
|
||||
|
||||
|
||||
def geometric_brownian_motion(
|
||||
S0: float,
|
||||
T: float,
|
||||
n_steps: int,
|
||||
n_paths: int = 1,
|
||||
mu: float = 0.0,
|
||||
sigma: float = 1.0,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Sample GBM exactly in log space, without Euler discretisation error."""
|
||||
_require_positive("S0", S0)
|
||||
if not np.isfinite(mu):
|
||||
raise ValueError("mu must be finite")
|
||||
_require_positive("sigma", sigma)
|
||||
t, log_returns = brownian_motion(
|
||||
T,
|
||||
n_steps,
|
||||
n_paths,
|
||||
drift=mu - 0.5 * sigma**2,
|
||||
sigma=sigma,
|
||||
rng=rng,
|
||||
)
|
||||
return t, S0 * np.exp(log_returns)
|
||||
|
||||
|
||||
def brownian_bridge(
|
||||
T: float,
|
||||
n_steps: int,
|
||||
n_paths: int = 1,
|
||||
start: float = 0.0,
|
||||
end: float = 0.0,
|
||||
sigma: float = 1.0,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Sample Brownian motion conditioned on its two endpoints."""
|
||||
if not np.isfinite(start) or not np.isfinite(end):
|
||||
raise ValueError("bridge endpoints must be finite")
|
||||
t, motion = brownian_motion(T, n_steps, n_paths, sigma=sigma, rng=rng)
|
||||
linear = start + (end - start) * (t / T)
|
||||
bridge = motion - np.outer(motion[:, -1], t / T) + linear[None, :]
|
||||
return t, bridge
|
||||
|
||||
|
||||
def merton_jump_diffusion(
|
||||
S0: float,
|
||||
T: float,
|
||||
n_steps: int,
|
||||
n_paths: int = 1,
|
||||
mu: float = 0.0,
|
||||
sigma: float = 0.2,
|
||||
jump_intensity: float = 0.0,
|
||||
jump_mean: float = -0.10,
|
||||
jump_sigma: float = 0.15,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Sample a Merton jump diffusion with annualised Poisson intensity."""
|
||||
_require_positive("S0", S0)
|
||||
_require_positive("sigma", sigma)
|
||||
_require_positive_int("n_paths", n_paths)
|
||||
if not np.isfinite(mu):
|
||||
raise ValueError("mu must be finite")
|
||||
if not np.isfinite(jump_intensity) or jump_intensity < 0:
|
||||
raise ValueError("jump_intensity must be finite and non-negative")
|
||||
if not np.isfinite(jump_mean):
|
||||
raise ValueError("jump_mean must be finite")
|
||||
if not np.isfinite(jump_sigma) or jump_sigma < 0:
|
||||
raise ValueError("jump_sigma must be finite and non-negative")
|
||||
|
||||
rng = _as_rng(rng)
|
||||
dt = T / n_steps
|
||||
compensator = np.exp(jump_mean + 0.5 * jump_sigma**2) - 1.0
|
||||
log_drift = mu - 0.5 * sigma**2 - jump_intensity * compensator
|
||||
t, diffusion = brownian_motion(
|
||||
T,
|
||||
n_steps,
|
||||
n_paths,
|
||||
drift=log_drift,
|
||||
sigma=sigma,
|
||||
rng=rng,
|
||||
)
|
||||
|
||||
counts = rng.poisson(jump_intensity * dt, size=(int(n_paths), int(n_steps)))
|
||||
jump_log = np.zeros((int(n_paths), int(n_steps)))
|
||||
max_count = int(counts.max()) if counts.size else 0
|
||||
if max_count:
|
||||
log_jumps = rng.normal(
|
||||
jump_mean,
|
||||
jump_sigma,
|
||||
size=(int(n_paths), int(n_steps), max_count),
|
||||
)
|
||||
mask = np.arange(max_count) < counts[:, :, None]
|
||||
jump_log = (log_jumps * mask).sum(axis=2)
|
||||
|
||||
cumulative_jumps = np.zeros((int(n_paths), int(n_steps) + 1))
|
||||
cumulative_jumps[:, 1:] = np.cumsum(jump_log, axis=1)
|
||||
return t, S0 * np.exp(diffusion + cumulative_jumps)
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Exact diffusion barrier formulae
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
|
||||
def first_passage_cdf_zero_drift(B: float, T: float, sigma: float) -> float:
|
||||
"""Return ``P(inf_{s<=T} X_s <= B)`` for zero-drift BM started at zero.
|
||||
|
||||
``B`` must be negative. For extremely remote barriers the floating-point
|
||||
result can underflow to zero; this correctly means "below machine precision".
|
||||
"""
|
||||
if not np.isfinite(B) or B >= 0:
|
||||
raise ValueError("B must be a finite negative lower barrier")
|
||||
_require_positive("T", T)
|
||||
_require_positive("sigma", sigma)
|
||||
return float(2.0 * norm.cdf(B / (sigma * np.sqrt(T))))
|
||||
|
||||
|
||||
def first_passage_cdf(B: float, T: float, nu: float, sigma: float) -> float:
|
||||
"""Bachelier--Lévy lower-barrier CDF for drifted BM started at zero."""
|
||||
if not np.isfinite(B) or B >= 0:
|
||||
raise ValueError("B must be a finite negative lower barrier")
|
||||
_require_positive("T", T)
|
||||
_require_positive("sigma", sigma)
|
||||
if not np.isfinite(nu):
|
||||
raise ValueError("nu must be finite")
|
||||
scale = sigma * np.sqrt(T)
|
||||
probability = norm.cdf((B - nu * T) / scale) + np.exp(
|
||||
2.0 * nu * B / sigma**2
|
||||
) * norm.cdf((B + nu * T) / scale)
|
||||
return float(np.clip(probability, 0.0, 1.0))
|
||||
|
||||
|
||||
def barrier_miss_prob(
|
||||
x0: float,
|
||||
xT: float,
|
||||
B: float,
|
||||
sigma: float,
|
||||
dt: float,
|
||||
) -> float:
|
||||
"""Conditional Brownian-bridge probability of crossing a lower barrier."""
|
||||
for name, value in (("x0", x0), ("xT", xT), ("B", B)):
|
||||
if not np.isfinite(value):
|
||||
raise ValueError(f"{name} must be finite")
|
||||
_require_positive("sigma", sigma)
|
||||
_require_positive("dt", dt)
|
||||
if x0 <= B or xT <= B:
|
||||
return 1.0
|
||||
exponent = -2.0 * (x0 - B) * (xT - B) / (sigma**2 * dt)
|
||||
return float(np.exp(exponent))
|
||||
|
||||
|
||||
def max_safe_dt(
|
||||
D: float | np.ndarray,
|
||||
sigma: float,
|
||||
epsilon: float,
|
||||
) -> float | np.ndarray:
|
||||
"""Invert the symmetric-endpoint bridge formula for ``dt``.
|
||||
|
||||
This is exact *conditional on both endpoint distances being ``D``*. In a
|
||||
live scheduler the future endpoint is unknown, so using the current distance
|
||||
for both endpoints is a local design approximation rather than a guarantee.
|
||||
"""
|
||||
_require_positive("sigma", sigma)
|
||||
if not np.isfinite(epsilon) or not 0.0 < epsilon < 1.0:
|
||||
raise ValueError("epsilon must lie strictly between 0 and 1")
|
||||
distances = np.asarray(D, dtype=float)
|
||||
if np.any(~np.isfinite(distances)) or np.any(distances < 0):
|
||||
raise ValueError("D must be finite and non-negative")
|
||||
result = 2.0 * distances**2 / (sigma**2 * np.log(1.0 / epsilon))
|
||||
return float(result) if result.ndim == 0 else result
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Monte Carlo estimators used by the notebooks
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
|
||||
def estimate_first_passage_prob(
|
||||
B: float,
|
||||
T: float,
|
||||
nu: float,
|
||||
sigma: float,
|
||||
n_paths: int,
|
||||
n_steps: int,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> float:
|
||||
"""Naive grid estimator; it has downward discretisation bias."""
|
||||
_, paths = brownian_motion(T, n_steps, n_paths, drift=nu, sigma=sigma, rng=rng)
|
||||
return float(np.mean(np.any(paths <= B, axis=1)))
|
||||
|
||||
|
||||
def estimate_first_passage_prob_bb(
|
||||
B: float,
|
||||
T: float,
|
||||
nu: float,
|
||||
sigma: float,
|
||||
n_paths: int,
|
||||
n_steps: int,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> float:
|
||||
"""Brownian-bridge-corrected first-passage Monte Carlo estimator."""
|
||||
rng = _as_rng(rng)
|
||||
_, paths = brownian_motion(T, n_steps, n_paths, drift=nu, sigma=sigma, rng=rng)
|
||||
dt = T / n_steps
|
||||
x0, x1 = paths[:, :-1], paths[:, 1:]
|
||||
both_above = (x0 > B) & (x1 > B)
|
||||
probabilities = np.where(
|
||||
both_above,
|
||||
np.exp(-2.0 * (x0 - B) * (x1 - B) / (sigma**2 * dt)),
|
||||
1.0,
|
||||
)
|
||||
uniforms = rng.uniform(size=probabilities.shape)
|
||||
return float(np.mean(np.any(uniforms < probabilities, axis=1)))
|
||||
|
||||
|
||||
def estimate_bridge_breach_prob(
|
||||
x0: float,
|
||||
xT: float,
|
||||
B: float,
|
||||
sigma: float,
|
||||
dt: float,
|
||||
n_paths: int,
|
||||
n_inner: int,
|
||||
rng: Optional[np.random.Generator] = None,
|
||||
) -> float:
|
||||
"""Monte Carlo check of :func:`barrier_miss_prob`."""
|
||||
_, paths = brownian_bridge(
|
||||
dt,
|
||||
n_inner,
|
||||
n_paths,
|
||||
start=x0,
|
||||
end=xT,
|
||||
sigma=sigma,
|
||||
rng=rng,
|
||||
)
|
||||
return float(np.mean(np.any(paths <= B, axis=1)))
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Adaptive schedule and detector evaluation
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
|
||||
def adaptive_schedule(
|
||||
t: np.ndarray,
|
||||
X: np.ndarray,
|
||||
B: float,
|
||||
sigma: float,
|
||||
eps: float,
|
||||
dt_cap: Optional[float] = None,
|
||||
B_upper: Optional[float] = None,
|
||||
) -> np.ndarray:
|
||||
"""Choose sample indices from a pre-generated path using a local proxy.
|
||||
|
||||
``t`` and ``sigma`` must use matching units: when ``t`` is measured in
|
||||
minutes, ``sigma`` must be the diffusion standard deviation per
|
||||
``sqrt(minute)``. At each observation, the current distance is substituted
|
||||
for both bridge endpoint distances in :func:`max_safe_dt`.
|
||||
"""
|
||||
times = np.asarray(t, dtype=float).ravel()
|
||||
values = np.asarray(X, dtype=float)
|
||||
if values.ndim > 1:
|
||||
if values.shape[0] != 1:
|
||||
raise ValueError("X must be one-dimensional or contain one path")
|
||||
values = values[0]
|
||||
values = values.ravel()
|
||||
if len(times) != len(values) or len(times) < 2:
|
||||
raise ValueError("t and X must have the same length of at least two")
|
||||
if np.any(~np.isfinite(times)) or np.any(np.diff(times) <= 0):
|
||||
raise ValueError("t must be finite and strictly increasing")
|
||||
if np.any(~np.isfinite(values)) or not np.isfinite(B):
|
||||
raise ValueError("X and B must be finite")
|
||||
_require_positive("sigma", sigma)
|
||||
if not 0.0 < eps < 1.0:
|
||||
raise ValueError("eps must lie strictly between 0 and 1")
|
||||
if dt_cap is not None:
|
||||
_require_positive("dt_cap", dt_cap)
|
||||
if B_upper is not None and (not np.isfinite(B_upper) or B_upper <= B):
|
||||
raise ValueError("B_upper must be finite and greater than B")
|
||||
|
||||
indices = [0]
|
||||
i = 0
|
||||
while i < len(times) - 1:
|
||||
lower_distance = values[i] - B
|
||||
distance = lower_distance
|
||||
if B_upper is not None:
|
||||
distance = min(lower_distance, B_upper - values[i])
|
||||
interval = max_safe_dt(max(float(distance), 1e-12), sigma, eps)
|
||||
if dt_cap is not None:
|
||||
interval = min(interval, dt_cap)
|
||||
target = times[i] + interval
|
||||
j = int(np.searchsorted(times, target, side="left"))
|
||||
j = min(max(j, i + 1), len(times) - 1)
|
||||
indices.append(j)
|
||||
i = j
|
||||
return np.asarray(indices, dtype=int)
|
||||
|
||||
|
||||
def fixed_uniform_indices(n: int, k: int) -> np.ndarray:
|
||||
"""Return exactly ``k`` approximately uniform indices from ``0`` to ``n-1``."""
|
||||
_require_positive_int("n", n)
|
||||
_require_positive_int("k", k)
|
||||
if k > n:
|
||||
raise ValueError("k cannot exceed n")
|
||||
if k == 1:
|
||||
return np.array([0], dtype=int)
|
||||
return np.rint(np.linspace(0, n - 1, k)).astype(int)
|
||||
|
||||
|
||||
def fixed_cadence_indices(t: np.ndarray, cadence: float) -> np.ndarray:
|
||||
"""Return grid indices for a fixed monitoring cadence.
|
||||
|
||||
The first and final grid points are always included. Intermediate target
|
||||
times are mapped to the first available grid point at or after each cadence
|
||||
tick. If the requested cadence is finer than the simulation grid, the
|
||||
resulting schedule is limited to one observation per grid point.
|
||||
"""
|
||||
times = np.asarray(t, dtype=float).ravel()
|
||||
if times.size < 2:
|
||||
raise ValueError("t must contain at least two points")
|
||||
if np.any(~np.isfinite(times)) or np.any(np.diff(times) <= 0):
|
||||
raise ValueError("t must be finite and strictly increasing")
|
||||
_require_positive("cadence", cadence)
|
||||
|
||||
targets = np.arange(times[0], times[-1] + cadence, cadence, dtype=float)
|
||||
targets = targets[targets <= times[-1] + 1e-12]
|
||||
indices = np.searchsorted(times, targets, side="left")
|
||||
indices = np.clip(indices, 0, len(times) - 1)
|
||||
indices = np.unique(indices.astype(int))
|
||||
if indices[0] != 0:
|
||||
indices = np.insert(indices, 0, 0)
|
||||
if indices[-1] != len(times) - 1:
|
||||
indices = np.append(indices, len(times) - 1)
|
||||
return indices
|
||||
|
||||
|
||||
def detect_barrier_event(
|
||||
sample_indices: np.ndarray,
|
||||
values: np.ndarray,
|
||||
breach_idx: Optional[int],
|
||||
barrier: float,
|
||||
direction: str,
|
||||
max_lag_steps: Optional[int] = None,
|
||||
) -> tuple[bool, Optional[int], Optional[int]]:
|
||||
"""Check whether a sampled point confirms a barrier event in time.
|
||||
|
||||
Returns ``(detected, lag_steps, detection_index)``. A detection must still
|
||||
lie beyond the barrier. When ``max_lag_steps`` is set, later observations
|
||||
do not count as catching the original event.
|
||||
"""
|
||||
if breach_idx is None:
|
||||
return False, None, None
|
||||
indices = np.asarray(sample_indices, dtype=int).ravel()
|
||||
path = np.asarray(values, dtype=float).ravel()
|
||||
if direction not in {"down", "up"}:
|
||||
raise ValueError("direction must be 'down' or 'up'")
|
||||
if max_lag_steps is not None:
|
||||
if isinstance(max_lag_steps, bool) or int(max_lag_steps) != max_lag_steps or max_lag_steps < 0:
|
||||
raise ValueError("max_lag_steps must be a non-negative integer or None")
|
||||
deadline = breach_idx + int(max_lag_steps)
|
||||
else:
|
||||
deadline = len(path) - 1
|
||||
|
||||
for sample_idx in indices:
|
||||
if sample_idx < breach_idx:
|
||||
continue
|
||||
if sample_idx > deadline:
|
||||
break
|
||||
beyond = path[sample_idx] <= barrier if direction == "down" else path[sample_idx] >= barrier
|
||||
if beyond:
|
||||
return True, int(sample_idx - breach_idx), int(sample_idx)
|
||||
return False, None, None
|
||||
|
||||
|
||||
def run_monte_carlo_simulation(
|
||||
S0: float = 100.0,
|
||||
sigma_annual: float = 0.30,
|
||||
mu_annual: float = 0.07,
|
||||
window_minutes: float = 1950.0,
|
||||
n_paths: int = 20,
|
||||
n_steps: int = 500,
|
||||
use_jumps: bool = False,
|
||||
jump_intensity: float = 25.0,
|
||||
jump_mean: float = -0.02,
|
||||
jump_sigma: float = 0.05,
|
||||
eps: float = 1e-3,
|
||||
dt_cap_minutes: Optional[float] = None,
|
||||
rng_seed: Optional[int] = None,
|
||||
drop_fraction: float = DROP_FRACTION,
|
||||
rise_fraction: float = RISE_FRACTION,
|
||||
max_detection_lag_steps: Optional[int] = 3,
|
||||
comparison_mode: str = "equal_budget",
|
||||
fixed_cadence_minutes: float = 60.0,
|
||||
) -> dict:
|
||||
"""Compare adaptive and fixed samplers on identical simulated paths.
|
||||
|
||||
``comparison_mode='equal_budget'`` gives the fixed baseline exactly the
|
||||
adaptive schedule's sample count on each path, isolating placement quality.
|
||||
``comparison_mode='fixed_cadence'`` samples independently at the requested
|
||||
cadence, exposing the detection-versus-observation-cost trade-off. Lower and
|
||||
upper breaches are evaluated independently, and a detection must occur
|
||||
within ``max_detection_lag_steps`` grid steps when that limit is not ``None``.
|
||||
"""
|
||||
_require_positive("S0", S0)
|
||||
_require_positive("sigma_annual", sigma_annual)
|
||||
_require_positive("window_minutes", window_minutes)
|
||||
_require_positive_int("n_paths", n_paths)
|
||||
_require_positive_int("n_steps", n_steps)
|
||||
if not np.isfinite(mu_annual):
|
||||
raise ValueError("mu_annual must be finite")
|
||||
if not 0.0 < drop_fraction < 1.0:
|
||||
raise ValueError("drop_fraction must lie between 0 and 1")
|
||||
if rise_fraction <= 0 or not np.isfinite(rise_fraction):
|
||||
raise ValueError("rise_fraction must be finite and positive")
|
||||
if max_detection_lag_steps is not None:
|
||||
if (
|
||||
isinstance(max_detection_lag_steps, bool)
|
||||
or int(max_detection_lag_steps) != max_detection_lag_steps
|
||||
or max_detection_lag_steps < 0
|
||||
):
|
||||
raise ValueError("max_detection_lag_steps must be a non-negative integer or None")
|
||||
if comparison_mode not in {"equal_budget", "fixed_cadence"}:
|
||||
raise ValueError("comparison_mode must be 'equal_budget' or 'fixed_cadence'")
|
||||
_require_positive("fixed_cadence_minutes", fixed_cadence_minutes)
|
||||
|
||||
rng = np.random.default_rng(rng_seed)
|
||||
horizon_years = minutes_to_years(window_minutes)
|
||||
if use_jumps:
|
||||
t_years, prices = merton_jump_diffusion(
|
||||
S0,
|
||||
horizon_years,
|
||||
n_steps,
|
||||
n_paths,
|
||||
mu=mu_annual,
|
||||
sigma=sigma_annual,
|
||||
jump_intensity=jump_intensity,
|
||||
jump_mean=jump_mean,
|
||||
jump_sigma=jump_sigma,
|
||||
rng=rng,
|
||||
)
|
||||
else:
|
||||
t_years, prices = geometric_brownian_motion(
|
||||
S0,
|
||||
horizon_years,
|
||||
n_steps,
|
||||
n_paths,
|
||||
mu=mu_annual,
|
||||
sigma=sigma_annual,
|
||||
rng=rng,
|
||||
)
|
||||
|
||||
times_minutes = np.asarray(t_years) * TRADING_MINUTES_PER_YEAR
|
||||
if dt_cap_minutes is None:
|
||||
dt_cap_minutes = max(window_minutes / 15.0, 5.0)
|
||||
_require_positive("dt_cap_minutes", dt_cap_minutes)
|
||||
|
||||
lower_price = float(S0 * (1.0 - drop_fraction))
|
||||
upper_price = float(S0 * (1.0 + rise_fraction))
|
||||
sigma_per_sqrt_minute = horizon_vol(sigma_annual, 1.0)
|
||||
|
||||
path_records: list[dict] = []
|
||||
adaptive_lags_steps: list[int] = []
|
||||
fixed_lags_steps: list[int] = []
|
||||
adaptive_lags_minutes: list[float] = []
|
||||
fixed_lags_minutes: list[float] = []
|
||||
adaptive_lags_by_direction: dict[str, list[int]] = {"lower": [], "upper": []}
|
||||
fixed_lags_by_direction: dict[str, list[int]] = {"lower": [], "upper": []}
|
||||
counts = {
|
||||
"lower_events": 0,
|
||||
"upper_events": 0,
|
||||
"adaptive_lower": 0,
|
||||
"adaptive_upper": 0,
|
||||
"fixed_lower": 0,
|
||||
"fixed_upper": 0,
|
||||
}
|
||||
|
||||
for path_prices in prices:
|
||||
log_relative = np.log(path_prices / path_prices[0])
|
||||
lower_log = float(np.log1p(-drop_fraction))
|
||||
upper_log = float(np.log1p(rise_fraction))
|
||||
adaptive_idx = adaptive_schedule(
|
||||
times_minutes,
|
||||
log_relative,
|
||||
lower_log,
|
||||
sigma_per_sqrt_minute,
|
||||
eps,
|
||||
dt_cap=dt_cap_minutes,
|
||||
B_upper=upper_log,
|
||||
)
|
||||
if comparison_mode == "equal_budget":
|
||||
fixed_idx = fixed_uniform_indices(len(times_minutes), len(adaptive_idx))
|
||||
else:
|
||||
fixed_idx = fixed_cadence_indices(times_minutes, fixed_cadence_minutes)
|
||||
|
||||
below = np.flatnonzero(path_prices <= lower_price)
|
||||
above = np.flatnonzero(path_prices >= upper_price)
|
||||
lower_breach_idx = int(below[0]) if below.size else None
|
||||
upper_breach_idx = int(above[0]) if above.size else None
|
||||
|
||||
lower_adapt = detect_barrier_event(
|
||||
adaptive_idx,
|
||||
path_prices,
|
||||
lower_breach_idx,
|
||||
lower_price,
|
||||
"down",
|
||||
max_detection_lag_steps,
|
||||
)
|
||||
upper_adapt = detect_barrier_event(
|
||||
adaptive_idx,
|
||||
path_prices,
|
||||
upper_breach_idx,
|
||||
upper_price,
|
||||
"up",
|
||||
max_detection_lag_steps,
|
||||
)
|
||||
lower_fixed = detect_barrier_event(
|
||||
fixed_idx,
|
||||
path_prices,
|
||||
lower_breach_idx,
|
||||
lower_price,
|
||||
"down",
|
||||
max_detection_lag_steps,
|
||||
)
|
||||
upper_fixed = detect_barrier_event(
|
||||
fixed_idx,
|
||||
path_prices,
|
||||
upper_breach_idx,
|
||||
upper_price,
|
||||
"up",
|
||||
max_detection_lag_steps,
|
||||
)
|
||||
|
||||
event_specs = [
|
||||
("lower", lower_breach_idx, lower_adapt, lower_fixed),
|
||||
("upper", upper_breach_idx, upper_adapt, upper_fixed),
|
||||
]
|
||||
for label, breach_idx, adaptive_result, fixed_result in event_specs:
|
||||
if breach_idx is None:
|
||||
continue
|
||||
counts[f"{label}_events"] += 1
|
||||
if adaptive_result[0]:
|
||||
counts[f"adaptive_{label}"] += 1
|
||||
adaptive_lags_steps.append(adaptive_result[1])
|
||||
adaptive_lags_by_direction[label].append(adaptive_result[1])
|
||||
adaptive_lags_minutes.append(
|
||||
float(times_minutes[adaptive_result[2]] - times_minutes[breach_idx])
|
||||
)
|
||||
if fixed_result[0]:
|
||||
counts[f"fixed_{label}"] += 1
|
||||
fixed_lags_steps.append(fixed_result[1])
|
||||
fixed_lags_by_direction[label].append(fixed_result[1])
|
||||
fixed_lags_minutes.append(
|
||||
float(times_minutes[fixed_result[2]] - times_minutes[breach_idx])
|
||||
)
|
||||
|
||||
first_candidates = [
|
||||
(idx, direction)
|
||||
for idx, direction in ((lower_breach_idx, "down"), (upper_breach_idx, "up"))
|
||||
if idx is not None
|
||||
]
|
||||
if first_candidates:
|
||||
first_breach_idx, first_breach_dir = min(first_candidates, key=lambda item: item[0])
|
||||
else:
|
||||
first_breach_idx, first_breach_dir = None, None
|
||||
|
||||
if first_breach_dir == "down":
|
||||
first_adapt, first_fixed = lower_adapt, lower_fixed
|
||||
elif first_breach_dir == "up":
|
||||
first_adapt, first_fixed = upper_adapt, upper_fixed
|
||||
else:
|
||||
first_adapt = first_fixed = (False, None, None)
|
||||
|
||||
path_records.append(
|
||||
{
|
||||
"prices": path_prices.tolist(),
|
||||
"log_prices": log_relative.tolist(),
|
||||
"sample_indices": adaptive_idx.tolist(),
|
||||
"sample_times": times_minutes[adaptive_idx].tolist(),
|
||||
"sample_prices": path_prices[adaptive_idx].tolist(),
|
||||
"fixed_sample_indices": fixed_idx.tolist(),
|
||||
"fixed_sample_times": times_minutes[fixed_idx].tolist(),
|
||||
"fixed_sample_prices": path_prices[fixed_idx].tolist(),
|
||||
"breach_idx": first_breach_idx,
|
||||
"breach_dir": first_breach_dir,
|
||||
"lower_breach_idx": lower_breach_idx,
|
||||
"upper_breach_idx": upper_breach_idx,
|
||||
"adaptive_detected": first_adapt[0],
|
||||
"fixed_detected": first_fixed[0],
|
||||
"adaptive_detection_lag": first_adapt[1],
|
||||
"fixed_detection_lag": first_fixed[1],
|
||||
"adaptive_lower_detected": lower_adapt[0],
|
||||
"adaptive_upper_detected": upper_adapt[0],
|
||||
"fixed_lower_detected": lower_fixed[0],
|
||||
"fixed_upper_detected": upper_fixed[0],
|
||||
"adaptive_lower_detection_lag": lower_adapt[1],
|
||||
"adaptive_upper_detection_lag": upper_adapt[1],
|
||||
"fixed_lower_detection_lag": lower_fixed[1],
|
||||
"fixed_upper_detection_lag": upper_fixed[1],
|
||||
}
|
||||
)
|
||||
|
||||
adaptive_total_samples = sum(len(record["sample_indices"]) for record in path_records)
|
||||
fixed_total_samples = sum(len(record["fixed_sample_indices"]) for record in path_records)
|
||||
n_events = counts["lower_events"] + counts["upper_events"]
|
||||
n_paths_with_any_breach = sum(record["breach_idx"] is not None for record in path_records)
|
||||
adaptive_detections = counts["adaptive_lower"] + counts["adaptive_upper"]
|
||||
fixed_detections = counts["fixed_lower"] + counts["fixed_upper"]
|
||||
|
||||
def _mean(values: list[float | int]) -> Optional[float]:
|
||||
return float(np.mean(values)) if values else None
|
||||
|
||||
return {
|
||||
"paths": path_records,
|
||||
"times_minutes": times_minutes.tolist(),
|
||||
"lower_barrier_price": lower_price,
|
||||
"upper_barrier_price": upper_price,
|
||||
"S0": S0,
|
||||
"sigma_annual": sigma_annual,
|
||||
"mu_annual": mu_annual,
|
||||
"window_minutes": window_minutes,
|
||||
"use_jumps": use_jumps,
|
||||
"eps": eps,
|
||||
"drop_fraction": float(drop_fraction),
|
||||
"rise_fraction": float(rise_fraction),
|
||||
"dt_cap_minutes": float(dt_cap_minutes),
|
||||
"max_detection_lag_steps": max_detection_lag_steps,
|
||||
"comparison_mode": comparison_mode,
|
||||
"fixed_cadence_minutes": float(fixed_cadence_minutes),
|
||||
"grid_step_minutes": float(times_minutes[1] - times_minutes[0]),
|
||||
"n_paths": n_paths,
|
||||
"n_paths_with_any_breach": n_paths_with_any_breach,
|
||||
"n_barrier_events": n_events,
|
||||
"n_lower_events": counts["lower_events"],
|
||||
"n_upper_events": counts["upper_events"],
|
||||
# Backward-compatible aliases retained for existing clients.
|
||||
"n_breaches": n_events,
|
||||
"n_lower_breaches": counts["lower_events"],
|
||||
"n_upper_breaches": counts["upper_events"],
|
||||
"adaptive_detections": adaptive_detections,
|
||||
"adaptive_lower_detections": counts["adaptive_lower"],
|
||||
"adaptive_upper_detections": counts["adaptive_upper"],
|
||||
"fixed_detections": fixed_detections,
|
||||
"fixed_lower_detections": counts["fixed_lower"],
|
||||
"fixed_upper_detections": counts["fixed_upper"],
|
||||
"adaptive_total_samples": adaptive_total_samples,
|
||||
"fixed_total_samples": fixed_total_samples,
|
||||
"mean_detection_lag": _mean(adaptive_lags_steps),
|
||||
"mean_lower_detection_lag": _mean(adaptive_lags_by_direction["lower"]),
|
||||
"mean_upper_detection_lag": _mean(adaptive_lags_by_direction["upper"]),
|
||||
"mean_fixed_detection_lag": _mean(fixed_lags_steps),
|
||||
"mean_fixed_lower_detection_lag": _mean(fixed_lags_by_direction["lower"]),
|
||||
"mean_fixed_upper_detection_lag": _mean(fixed_lags_by_direction["upper"]),
|
||||
"mean_detection_lag_minutes": _mean(adaptive_lags_minutes),
|
||||
"mean_fixed_detection_lag_minutes": _mean(fixed_lags_minutes),
|
||||
"model_scope": (
|
||||
"eps controls a local Brownian-diffusion scheduling proxy; it is not an "
|
||||
"unconditional miss guarantee and does not control jump risk."
|
||||
),
|
||||
}
|
||||
Reference in New Issue
Block a user