Building long-short factor returns from a cross-section, and checking exposures

"""M06A · Building long-short factor returns from a cross-section, and checking exposures.

Exercise stubs. Replace each ``raise NotImplementedError`` with your code, then run the tests in
../tests. See notes.md §6A.12 (Propositions 6A.19–6A.20).

Timing convention (it is the whole game in factor research):
* Row t of every panel is the END of month t. A characteristic in row t is known at the end of t.
* A portfolio formed on characteristics in row t earns the returns in row t + 1.
* The factor return series is therefore indexed by the holding month t + 1.
"""

from __future__ import annotations

import numpy as np
import pandas as pd
import statsmodels.api as sm  # noqa: F401  (for factor_regression)


# ---------------------------------------------------------------------------------------------
# Provided helper (not a stub): a synthetic cross-section with known premia.
# ---------------------------------------------------------------------------------------------
def make_cross_section(n_stocks: int = 400, n_months: int = 240, seed: int = 0,
                       size_premium: float = -0.002, value_premium: float = 0.004,
                       alpha_sd: float = 0.006, alpha_rho: float = 0.9):
    """Synthetic monthly panel.

    Characteristics (cross-sectionally standardised each month, persistent over time):
        size_z (think log market cap), bm_z (think log book-to-market).
    Latent persistent expected-return component a_{i,t} (AR(1), coefficient alpha_rho,
    stationary sd alpha_sd) that is NOT observed: past returns are a noisy proxy for it, which is
    what creates momentum.
    Returns:
        r_{i,t+1} = beta_i * m_{t+1} + lam_s,t+1 * size_z_{i,t} + lam_v,t+1 * bm_z_{i,t}
                    + a_{i,t} + e_{i,t+1}
    with m ~ N(0.008, 0.05^2), lam_s ~ N(size_premium, 0.01^2), lam_v ~ N(value_premium, 0.01^2),
    e ~ N(0, 0.08^2).
    Returns a dict of DataFrames (index = month-end dates, columns = stock ids):
        'returns', 'size', 'bm', 'alpha'  and a Series 'market'.
    """
    rng = np.random.default_rng(seed)
    T, N = n_months, n_stocks
    idx = pd.period_range("2000-01", periods=T, freq="M").to_timestamp("M")
    cols = [f"S{i:03d}" for i in range(N)]

    def persistent(rho):
        x = np.empty((T, N))
        x[0] = rng.standard_normal(N)
        for t in range(1, T):
            x[t] = rho * x[t - 1] + np.sqrt(1 - rho ** 2) * rng.standard_normal(N)
        x = (x - x.mean(axis=1, keepdims=True)) / x.std(axis=1, keepdims=True)
        return x

    size = persistent(0.98)
    bm = persistent(0.97)
    a = np.empty((T, N))
    a[0] = alpha_sd * rng.standard_normal(N)
    for t in range(1, T):
        a[t] = alpha_rho * a[t - 1] + alpha_sd * np.sqrt(1 - alpha_rho ** 2) * rng.standard_normal(N)
    beta = 1.0 + 0.2 * rng.standard_normal(N)
    m = rng.normal(0.008, 0.05, T)
    lam_s = rng.normal(size_premium, 0.01, T)
    lam_v = rng.normal(value_premium, 0.01, T)
    r = np.zeros((T, N))
    r[0] = beta * m[0] + rng.normal(0, 0.08, N)
    for t in range(1, T):
        r[t] = (beta * m[t] + lam_s[t] * size[t - 1] + lam_v[t] * bm[t - 1] + a[t - 1]
                + rng.normal(0, 0.08, N))
    return {
        "returns": pd.DataFrame(r, idx, cols),
        "size": pd.DataFrame(size, idx, cols),
        "bm": pd.DataFrame(bm, idx, cols),
        "alpha": pd.DataFrame(a, idx, cols),
        "market": pd.Series(m, idx, name="market"),
    }


# ---------------------------------------------------------------------------------------------
# Exercises
# ---------------------------------------------------------------------------------------------
def momentum_signal(returns: pd.DataFrame, lookback: int = 11, skip: int = 1) -> pd.DataFrame:
    """Momentum characteristic known at the end of month t:

        mom_{i,t} = prod_{s = t-skip-lookback+1}^{t-skip} (1 + r_{i,s}) - 1.

    With lookback=11, skip=1 this is the return over months t-11..t-1, skipping month t
    (the IIMA library's definition: 11 months from the end of t-12 to the end of t-1).
    Rows without a full window are NaN. Must not use any return from row t or later.
    Hint: np.log1p, DataFrame.shift, DataFrame.rolling(...).sum(), np.expm1.
    """
    raise NotImplementedError


def long_short(char: pd.DataFrame, returns: pd.DataFrame, lo: float = 0.3, hi: float = 0.7,
               direction: int = 1) -> pd.Series:
    """Equal-weighted top-minus-bottom portfolio return.

    For each month t with characteristic row t: long stocks with char >= quantile(hi), short
    stocks with char <= quantile(lo) (quantiles over stocks with a non-NaN characteristic and a
    non-NaN next-month return), and record the return difference in month t+1 (indexed by t+1).
    Skip months with fewer than 10 valid stocks. direction=-1 returns bottom-minus-top (use it for
    size: small minus big).
    """
    raise NotImplementedError


def ff_2x3(size: pd.DataFrame, char: pd.DataFrame, returns: pd.DataFrame) -> pd.DataFrame:
    """Fama-French style 2x3 double sort, equal-weighted within each of the six portfolios.

    Each month t: Small = size below the cross-sectional median, Big = the rest; Low/Mid/High =
    char <= 30th percentile / between / >= 70th (breakpoints over all valid stocks). Skip months
    with fewer than 30 valid stocks. Holding month t+1 returns:
        SMB = (SL + SM + SH)/3 - (BL + BM + BH)/3
        HML = (SH + BH)/2 - (SL + BL)/2
    Returns a DataFrame with columns ['SMB', 'HML'] indexed by holding month.
    """
    raise NotImplementedError


def factor_regression(y: pd.Series, X: pd.DataFrame, nw_lags: int = 0) -> dict:
    """Time-series regression y_t = alpha + sum_j b_j X_{j,t} + e_t on the common non-NaN dates.

    Standard errors: OLS if nw_lags == 0, else Newey-West with Bartlett weights
    w_l = 1 - l/(nw_lags+1). Returns a dict with keys
        'alpha', 'betas' (Series indexed by X.columns), 't_alpha', 't_betas' (Series), 'r2', 'n'.

    The point is the factor model, not the linear algebra: use statsmodels,
        sm.OLS(y, sm.add_constant(X)).fit()                                  (classical), or
        .fit(cov_type="HAC", cov_kwds={"maxlags": nw_lags, "use_correction": True})  (Newey-West).
    """
    raise NotImplementedError


def annualised_stats(f: pd.Series) -> dict:
    """Annualised mean (12 x monthly mean), volatility (sqrt(12) x monthly sd, ddof=1), Sharpe
    (mean/vol), t-statistic of the monthly mean (mean / (sd/sqrt(n))), and maximum drawdown of the
    compounded series prod(1 + f). Returns a dict with keys 'mean', 'vol', 'sharpe', 't', 'max_dd'
    (max_dd as a negative fraction)."""
    raise NotImplementedError