The Lo-MacKinlay variance ratio

"""Exercise: the Lo-MacKinlay variance ratio (notes §2.3, Propositions 2.6-2.8).

Implement every function, then run from the repo root:

    cd exercises && .venv/bin/python -m pytest modules/M02-statistics-of-returns/tests/test_m02_variance_ratio.py

numpy only. Returns are 1-D arrays (or anything np.asarray accepts).
"""

from __future__ import annotations

import numpy as np


def kperiod_variance_ar1(phi: float, sigma2: float, k: int) -> float:
    """Variance of the k-period return r_1 + ... + r_k when one-period returns
    follow a stationary AR(1), r_t = phi r_{t-1} + e_t, with Var(r_t) = sigma2.

    Closed form (notes eq. 2.6 with rho_j = phi^j; eq. 2.8 sums the series):
        Var = sigma2 * [ k + 2 * sum_{j=1}^{k-1} (k - j) * phi^j ].

    Inputs: |phi| < 1, sigma2 > 0, integer k >= 1. Output: float.
    """
    raise NotImplementedError


def vr_ar1_theory(phi: float, q: int) -> float:
    """Population variance ratio VR(q) = Var(q-period return) / (q Var(r_t))
    for AR(1) returns:
        VR(q) = 1 + 2 * sum_{j=1}^{q-1} (1 - j/q) * phi^j     (notes eqs. 2.7, 2.9).
    """
    raise NotImplementedError


def variance_ratio(r, q: int) -> tuple[float, float, float]:
    """Lo-MacKinlay (1988) variance-ratio test with overlapping q-period returns.

    With T returns r_1..r_T and mu = mean(r):
        sigma_a^2 = sum (r_t - mu)^2 / (T - 1)
        sigma_c^2 = (1/m) sum_{t=q}^{T} (r_t + ... + r_{t-q+1} - q mu)^2,
                    m = q (T - q + 1) (1 - q/T)
        VR = sigma_c^2 / sigma_a^2
    Homoskedastic statistic:
        z = (VR - 1) / sqrt( 2 (2q - 1)(q - 1) / (3 q T) )
    Heteroskedasticity-robust statistic:
        delta_j = T * sum_{t=j+1}^{T} d_t^2 d_{t-j}^2 / (sum_t d_t^2)^2,  d_t = r_t - mu
        theta   = sum_{j=1}^{q-1} (2 (q - j) / q)^2 delta_j
        z*      = sqrt(T) (VR - 1) / sqrt(theta)
    Both are asymptotically N(0,1) under their null (notes Prop. 2.8).

    Raise ValueError unless 2 <= q < T. Return (VR, z, z*) as floats.
    Hint: overlapping sums are differences of a cumulative sum.
    """
    raise NotImplementedError


def rejection_rate(sim_returns, q: int, n_sims: int, seed: int, robust: bool = True, level: float = 1.96) -> float:
    """Monte Carlo size/power of the variance-ratio test.

    Create rng = np.random.default_rng(seed). Repeat n_sims times: draw
    x = sim_returns(rng), compute variance_ratio(x, q), and count |z*| > level
    (or |z| > level if robust is False). Return the fraction rejected.
    Under a true null this is the test's size (about 0.05); under an
    alternative it is the power.
    """
    raise NotImplementedError


def simulate_ar1(n: int, phi: float, rng, sigma: float = 0.01) -> np.ndarray:
    """Simulate n draws of a stationary Gaussian AR(1) with unconditional sd sigma:
        r_0 ~ N(0, sigma^2),  r_t = phi r_{t-1} + e_t,  e_t ~ N(0, sigma^2 (1 - phi^2)).
    Draw the n innovations e first with rng.normal, then r_0, then iterate.
    """
    raise NotImplementedError