Realised-volatility estimators

"""M09 exercise 1: realised-volatility estimators (notes §9.2).

All estimators return a VARIANCE PER DAY (not annualised, not a vol) unless
the name says otherwise. Inputs are numpy arrays of daily open, high, low,
close prices (same length n). Where a previous close is needed, it is passed
explicitly as `prev_close`.

Notation (notes eq. 9.4):
    o_t = ln(O_t / C_{t-1})   overnight (opening) log return
    u_t = ln(H_t / O_t)       normalised high
    d_t = ln(L_t / O_t)       normalised low
    c_t = ln(C_t / O_t)       open-to-close log return

Run:  cd exercises && .venv/bin/python -m pytest modules/M09-volatility-trading/tests/test_m09_realised_vol.py
"""

from __future__ import annotations

import math

import numpy as np

LN2 = math.log(2.0)


# ------------------------------------------------------------------ simulation

def simulate_ohlc(
    n_days: int,
    sigma_daily: float,
    n_intraday: int = 390,
    overnight_frac: float = 0.0,
    drift_daily: float = 0.0,
    s0: float = 100.0,
    seed: int = 0,
) -> dict[str, np.ndarray]:
    """Simulate daily OHLC bars from a Brownian log price (notes §9.2.4).

    Total daily log variance is sigma_daily**2. A fraction `overnight_frac` of it
    is realised as a N(0, overnight_frac * sigma_daily**2) jump between yesterday's
    close and today's open; the remaining (1 - overnight_frac) * sigma_daily**2
    diffuses during the session as a random walk on `n_intraday` equal Gaussian
    steps (each with mean drift_daily / n_intraday). High and low are the max and
    min of the sampled intraday log path INCLUDING the open and the close.

    Use np.random.default_rng(seed). Draw the overnight jumps first (only if
    overnight_frac > 0), then the (n_days, n_intraday) step matrix.

    Returns dict with keys 'open', 'high', 'low', 'close' (arrays of length
    n_days) and 'prev_close' (length n_days; prev_close[0] = s0).
    """
    raise NotImplementedError


# ------------------------------------------------------------------ per-day pieces

def daily_close_to_close_sq(close: np.ndarray, prev_close: np.ndarray) -> np.ndarray:
    """Per-day close-to-close estimator r_t**2 with r_t = ln(C_t / C_{t-1}) (zero mean assumed)."""
    raise NotImplementedError


def daily_parkinson(high: np.ndarray, low: np.ndarray) -> np.ndarray:
    """Per-day Parkinson estimator (ln(H/L))**2 / (4 ln 2). Notes Prop. 9.2."""
    raise NotImplementedError


def daily_garman_klass(open_: np.ndarray, high: np.ndarray, low: np.ndarray, close: np.ndarray) -> np.ndarray:
    """Per-day Garman-Klass estimator 0.5 (ln H/L)**2 - (2 ln 2 - 1)(ln C/O)**2. Notes eq. (9.7)."""
    raise NotImplementedError


def daily_rogers_satchell(open_: np.ndarray, high: np.ndarray, low: np.ndarray, close: np.ndarray) -> np.ndarray:
    """Per-day Rogers-Satchell estimator u(u - c) + d(d - c). Notes eq. (9.8)."""
    raise NotImplementedError


# ------------------------------------------------------------------ window estimators

def close_to_close_var(close: np.ndarray) -> float:
    """Sample variance (ddof=1) of the daily log returns ln(C_t/C_{t-1}) within `close`
    (so n closes give n-1 returns)."""
    raise NotImplementedError


def parkinson_var(high: np.ndarray, low: np.ndarray) -> float:
    """Mean of daily Parkinson estimates over the window."""
    raise NotImplementedError


def garman_klass_var(open_, high, low, close) -> float:
    """Mean of daily Garman-Klass estimates over the window."""
    raise NotImplementedError


def rogers_satchell_var(open_, high, low, close) -> float:
    """Mean of daily Rogers-Satchell estimates over the window."""
    raise NotImplementedError


def yang_zhang_var(open_, high, low, close, prev_close) -> float:
    """Yang-Zhang (2000) variance per day over a window of n days (notes eq. 9.9).

        V_O  = sample var (ddof=1) of o_t = ln(O_t / C_{t-1})
        V_C  = sample var (ddof=1) of c_t = ln(C_t / O_t)
        V_RS = mean Rogers-Satchell
        k    = 0.34 / (1.34 + (n + 1) / (n - 1))
        V_YZ = V_O + k V_C + (1 - k) V_RS
    """
    raise NotImplementedError


def annualise_vol(var_per_day: float, periods: int = 252) -> float:
    """sqrt(var_per_day * periods)."""
    raise NotImplementedError


# ------------------------------------------------------------------ range moments

def range_moment_exact(p: float, sigma: float = 1.0, t: float = 1.0) -> float:
    """E[R^p] for the range R = max - min of a driftless BM with vol sigma on [0, t].

    Notes Prop. 9.1: E[R^p] = 4 * E|Z|^p * eta(p - 1) * (sigma * sqrt(t))**p, with
    E|Z|^p = 2^{p/2} Gamma((p+1)/2) / sqrt(pi) and eta(s) = sum_{k>=1} (-1)^{k-1} k^{-s}
    the Dirichlet eta function; eta(0) = 1/2 (Abel sum), eta(1) = ln 2, and for other
    s > 0 sum enough terms (average two consecutive partial sums to accelerate).
    Valid for p >= 1. Check: p = 2 gives 4 ln 2, p = 4 gives 9 zeta(3).
    """
    raise NotImplementedError


def relative_efficiency(per_day_estimates: np.ndarray, benchmark: np.ndarray) -> float:
    """Var(benchmark) / Var(estimator) (ddof=1), both arrays of per-day variance estimates."""
    raise NotImplementedError