The return distribution of systematic option selling

"""M09 exercise 4: the return distribution of systematic option selling (notes §9.9).

We sell a weekly strangle on a synthetic index whose daily log returns come from a
GARCH(1,1) with Student-t shocks and rare negative jumps (toolkit.data.synthetic_returns).
The strangle is struck k implied standard deviations either side of spot, priced with
Black-Scholes at an implied vol equal to `iv_markup` times a trailing realised-vol
forecast, and held to expiry (no hedging, no early exit).

`iv_markup` is a crude knob for the variance risk premium: markup 1.0 prices options
at the forecast of realised vol; markup 1.2 sets implied vol 20% above it.

Provided: bs_price(S, K, tau, r, sigma, kind="call") from toolkit.options (you wrote it
in M08; it returns intrinsic value at tau <= 0), and scipy.stats.skew (population-moment
skewness, the statistic you derived in M02) for the shape of the P&L distribution.

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

from __future__ import annotations

import math

import numpy as np
from scipy.stats import skew  # noqa: F401  (provided, M02)

from toolkit.options import bs_price  # noqa: F401  (provided, M08)


# ------------------------------------------------------------------ exercises

def short_strangle_pnl(log_returns, hold=5, k=1.0, iv_markup=1.1, vol_window=20, periods=252, r=0.0):
    """Per-trade P&L of a sequence of non-overlapping short strangles, as a fraction of spot.

    Prices: S_t = exp(sum of the first t log returns), so S_0 = 1.
    Trades open at t = vol_window, vol_window + hold, ... while t + hold <= n.
    At open time t:
        sigma_f = sqrt(periods) * std(log_returns[t - vol_window : t], ddof=1)   (forecast)
        sigma_i = iv_markup * sigma_f                                           (implied)
        tau     = hold / periods
        K_put   = S_t exp(-k sigma_i sqrt(tau)),  K_call = S_t exp(+k sigma_i sqrt(tau))
        premium = BS put(K_put) + BS call(K_call), both at sigma_i, rate r
    At t + hold: payoff = max(K_put - S, 0) + max(S - K_call, 0).
    P&L / S_t = (premium * e^{r tau} - payoff) / S_t.

    Returns a numpy array of P&L fractions, one per trade.
    """
    raise NotImplementedError


def gaussian_kelly(x) -> float:
    """Mean-variance Kelly fraction mean(x) / var(x) (ddof=0): the f* = mu / sigma^2 of M03."""
    raise NotImplementedError


def growth_rate(x, f) -> float:
    """Empirical log growth per trade g(f) = mean(log(1 + f x)); return -inf if any 1 + f x <= 0."""
    raise NotImplementedError


def kelly_fraction(x, tol=1e-10) -> float:
    """Exact empirical Kelly fraction argmax_f mean(log(1 + f x)) over 0 <= f < 1/max(-x).

    g is concave in f, so solve g'(f) = mean(x / (1 + f x)) = 0 by bisection on
    [0, (1 - 1e-12) / max(-x)]. Return 0 if mean(x) <= 0. Assumes min(x) < 0.
    Notes Prop. 9.11.
    """
    raise NotImplementedError