GARCH(1,1) by maximum likelihood, from scratch

"""M13 · GARCH(1,1) by maximum likelihood, from scratch. Read notes §13.4 first.

Model:
    r_t = eps_t,  eps_t = sigma_t z_t,  z_t iid N(0, 1)
    sigma_t^2 = omega + alpha eps_{t-1}^2 + beta sigma_{t-1}^2
with omega > 0, alpha >= 0, beta >= 0, alpha + beta < 1.

Initialisation ("backcast"): set the unobserved pre-sample eps_0^2 and sigma_0^2
both equal to a number b, so sigma_1^2 = omega + (alpha + beta) b. This is the
convention of the `arch` package; the tests pass arch's own backcast and expect
your log-likelihood to agree with arch's to ~1e-8 relative.

Run:  cd exercises && .venv/bin/python -m pytest modules/M13-forecasting/tests/test_m13_garch.py
"""

from __future__ import annotations

import numpy as np


def simulate_garch11(omega: float, alpha: float, beta: float, n: int, seed: int = 0,
                     burn: int = 500) -> np.ndarray:
    """Simulate n returns from a zero-mean Gaussian GARCH(1,1), discarding `burn` start-up values. (Given.)"""
    rng = np.random.default_rng(seed)
    z = rng.standard_normal(n + burn)
    r = np.empty(n + burn)
    s2 = omega / (1.0 - alpha - beta)
    for t in range(n + burn):
        if t > 0:
            s2 = omega + alpha * r[t - 1] ** 2 + beta * s2
        r[t] = np.sqrt(s2) * z[t]
    return r[burn:]


def default_backcast(r: np.ndarray) -> float:
    """Exponentially weighted mean of the first min(75, n) squared returns, weights 0.94^j, normalised. (Given.)"""
    r = np.asarray(r, dtype=float)
    tau = min(75, r.size)
    w = 0.94 ** np.arange(tau)
    w = w / w.sum()
    return float(np.sum(w * r[:tau] ** 2))


def garch11_variance(r: np.ndarray, omega: float, alpha: float, beta: float,
                     backcast: float | None = None) -> np.ndarray:
    """Conditional variances sigma_t^2 for t = 1..n (notes eq. 13.16).

    Inputs: returns r (length n), parameters, backcast b (default: default_backcast(r)).
    Output: array of length n with sigma_1^2 = omega + (alpha + beta) b and
            sigma_t^2 = omega + alpha r_{t-1}^2 + beta sigma_{t-1}^2 for t >= 2.
    Hint: a Python loop is fine; scipy.signal.lfilter is faster.
    """
    raise NotImplementedError


def garch11_loglik(params, r: np.ndarray, backcast: float | None = None) -> float:
    """Full Gaussian log-likelihood (notes eq. 13.20), including the constant:

        L(theta) = -1/2 sum_{t=1}^n [ ln(2 pi) + ln sigma_t^2 + r_t^2 / sigma_t^2 ]

    params = (omega, alpha, beta).
    """
    raise NotImplementedError


def fit_garch11(r: np.ndarray, backcast: float | None = None) -> dict:
    """Maximum-likelihood estimates. Return dict(omega, alpha, beta, loglik).

    Enforce omega > 0, alpha >= 0, beta >= 0, alpha + beta < 1, either with
    bounds/constraints or (recommended) an unconstrained reparametrisation such as
    omega = s exp(t0), alpha + beta = logistic(t1), alpha / (alpha + beta) = logistic(t2).
    Use scipy.optimize.minimize from more than one starting point.
    """
    raise NotImplementedError


def unconditional_variance(omega: float, alpha: float, beta: float) -> float:
    """sigma_bar^2 = omega / (1 - alpha - beta) (Proposition 13.10). Raise ValueError if alpha + beta >= 1."""
    raise NotImplementedError


def garch_half_life(alpha: float, beta: float) -> float:
    """Half-life of a variance shock: h solving (alpha + beta)^h = 1/2 (eq. 13.18).
    (Not the OU half-life ln 2 / kappa of M10/M13A, toolkit.ou.ou_half_life.)"""
    raise NotImplementedError


def garch11_forecast(omega: float, alpha: float, beta: float, sigma2_next: float, h: int) -> np.ndarray:
    """Multi-step forecasts E_T[sigma_{T+k}^2], k = 1..h (eq. 13.17):

        sigma_bar^2 + (alpha + beta)^{k-1} (sigma_{T+1}^2 - sigma_bar^2).
    """
    raise NotImplementedError


def mc_forecast(omega: float, alpha: float, beta: float, sigma2_next: float, h: int,
                n_paths: int = 200_000, seed: int = 0) -> np.ndarray:
    """Monte Carlo estimate of E_T[eps_{T+k}^2], k = 1..h.

    Start every path at sigma_{T+1}^2 = sigma2_next; at each step draw
    eps^2 = sigma^2 z^2, record its cross-path mean, then update
    sigma^2 <- omega + alpha eps^2 + beta sigma^2. Vectorise over paths.
    """
    raise NotImplementedError