spb/wp7_uqo Public

1# =============================================================================2# Author: Simon-Pierre Boucher3# Contact: contact@spboucher.ai4# =============================================================================5"""Econometric primitives shared by every analysis script.67All estimators are deliberately implemented with :func:`numpy.linalg.lstsq`8and explicit sandwich formulas — exactly as in the original study — so that9re-running the pipeline reproduces the published numbers bit-for-bit. Do not10"modernize" the numerics (e.g. switch to statsmodels) without re-validating11every result table against ``results/``.12"""1314import numpy as np15import pandas as pd16from numpy.linalg import lstsq171819# --------------------------------------------------------------------------20# Transformations21# --------------------------------------------------------------------------22def winsorize(s: pd.Series, q: float = 0.01) -> pd.Series:23    """Clip a series at its ``q`` and ``1-q`` empirical quantiles."""24    lo, hi = s.quantile(q), s.quantile(1 - q)25    return s.clip(lo, hi)262728def standardize(X: np.ndarray) -> np.ndarray:29    """Column-standardize (population std, zero-variance columns left at 1)."""30    means = X.mean(axis=0)31    stds = X.std(axis=0)32    stds[stds == 0] = 133    return (X - means) / stds343536def add_constant(X: np.ndarray) -> np.ndarray:37    """Prepend an intercept column of ones."""38    return np.column_stack([np.ones(len(X)), X])394041def design_matrix(data: pd.DataFrame, features: list) -> np.ndarray:42    """Standardized feature matrix with intercept, ready for ``lstsq``."""43    return add_constant(standardize(data[features].values))444546# --------------------------------------------------------------------------47# OLS with robust inference48# --------------------------------------------------------------------------49def ols(X: np.ndarray, y: np.ndarray) -> np.ndarray:50    """Least-squares coefficients."""51    coefs, _, _, _ = lstsq(X, y, rcond=None)52    return coefs535455def r_squared(y: np.ndarray, y_pred: np.ndarray) -> float:56    ss_res = np.sum((y - y_pred) ** 2)57    ss_tot = np.sum((y - y.mean()) ** 2)58    return 1 - ss_res / ss_tot if ss_tot > 0 else 0596061def adjusted_r_squared(r2: float, n: int, k: int) -> float:62    return 1 - (1 - r2) * (n - 1) / (n - k - 1)636465def _xtx_inverse(X: np.ndarray) -> np.ndarray:66    try:67        return np.linalg.inv(X.T @ X)68    except np.linalg.LinAlgError:69        return np.linalg.pinv(X.T @ X)707172def hc1_tstats(X: np.ndarray, y: np.ndarray, coefs: np.ndarray,73               n_features: int) -> tuple:74    """Heteroskedasticity-robust (HC1) standard errors and t-statistics.7576    ``n_features`` excludes the intercept; the small-sample factor is77    ``n / (n - k - 1)`` as in the original study.78    """79    n = len(y)80    e = y - X @ coefs81    XtX_inv = _xtx_inverse(X)82    # X' diag(e²) X computed by broadcasting (identical to the dense-diagonal83    # formula, but O(nk) memory instead of O(n²)).84    S = (X * (e ** 2)[:, None]).T @ X * n / (n - n_features - 1)85    se = np.sqrt(np.abs(np.diag(XtX_inv @ S @ XtX_inv)))86    t_stats = coefs / np.where(se > 0, se, 1)87    return se, t_stats888990def newey_west_tstats(X: np.ndarray, y: np.ndarray, coefs: np.ndarray,91                      n_lags: int = 5) -> tuple:92    """Newey-West HAC standard errors with Bartlett kernel weights."""93    n, k = X.shape94    e = y - X @ coefs95    XtX_inv = _xtx_inverse(X)9697    S = np.zeros((k, k))98    for j in range(n):99        S += e[j] ** 2 * np.outer(X[j], X[j])100    for lag in range(1, n_lags + 1):101        w = 1 - lag / (n_lags + 1)102        for j in range(lag, n):103            cross = e[j] * e[j - lag] * np.outer(X[j], X[j - lag])104            S += w * (cross + cross.T)105106    V = XtX_inv @ S @ XtX_inv107    se = np.sqrt(np.abs(np.diag(V)))108    t_stats = coefs / np.where(se > 0, se, 1)109    return se, t_stats110111112def double_clustered_tstats(X: np.ndarray, y: np.ndarray, coefs: np.ndarray,113                            tickers: np.ndarray, dates: np.ndarray) -> tuple:114    """Two-way (ticker + date) clustered SEs, Cameron-Gelbach-Miller (2011)."""115    k = X.shape[1]116    e = y - X @ coefs117    XtX_inv = _xtx_inverse(X)118119    def cluster_meat(labels: np.ndarray) -> np.ndarray:120        S = np.zeros((k, k))121        for g in np.unique(labels):122            mask = labels == g123            u = (X[mask].T * e[mask]).sum(axis=1, keepdims=True)124            S += u @ u.T125        return S126127    S_hc = (X * (e ** 2)[:, None]).T @ X128    V = XtX_inv @ (cluster_meat(tickers) + cluster_meat(dates) - S_hc) @ XtX_inv129    se = np.sqrt(np.abs(np.diag(V)))130    t_stats = coefs / np.where(se > 0, se, 1)131    return se, t_stats132133134def quantile_regression(X: np.ndarray, y: np.ndarray, tau: float,135                        max_iter: int = 50) -> np.ndarray:136    """Quantile regression via iteratively reweighted least squares."""137    coefs = ols(X, y)138    for _ in range(max_iter):139        residuals = y - X @ coefs140        weights = np.where(residuals >= 0, tau, 1 - tau)141        weights = np.maximum(weights / (np.abs(residuals) + 1e-6), 1e-6)142        Xw = X * weights[:, None]143        try:144            coefs = np.linalg.solve(Xw.T @ X, Xw.T @ y)145        except np.linalg.LinAlgError:146            break147    return coefs148149150def pooled_regression_summary(data: pd.DataFrame, features: list, target: str,151                              label: str = "", min_obs: int = 100):152    """Winsorize → standardize → pooled OLS with HC1; compact summary dict.153154    Returns ``None`` when fewer than ``min_obs`` complete observations exist.155    """156    sub = data[features + [target]].dropna()157    if len(sub) < min_obs:158        return None159    for f in features:160        sub[f] = winsorize(sub[f])161    sub[target] = winsorize(sub[target])162163    X = design_matrix(sub, features)164    y = sub[target].values165    coefs = ols(X, y)166    r2 = r_squared(y, X @ coefs)167    _, t = hc1_tstats(X, y, coefs, len(features))168    return {169        'label': label,170        'n_obs': len(sub),171        'r2': r2,172        'adj_r2': adjusted_r_squared(r2, len(y), len(features)),173        'n_significant': int(np.sum(np.abs(t[1:]) > 1.96)),174    }175176177# --------------------------------------------------------------------------178# Fama-MacBeth179# --------------------------------------------------------------------------180def fama_macbeth(data: pd.DataFrame, features: list, target: str,181                 min_obs_per_day: int = 10, min_periods: int = 30):182    """Daily cross-sectional regressions; time-series means and t-statistics.183184    Returns ``(results_df, n_periods, n_tickers)`` or ``None`` when the panel185    is too small.186    """187    sub = data[features + [target, 'ticker', 'trade_date']].dropna()188    if len(sub) < 100:189        return None190    for f in features:191        sub[f] = winsorize(sub[f])192    sub[target] = winsorize(sub[target])193194    coef_series = []195    for dt in sub['trade_date'].unique():196        day = sub[sub['trade_date'] == dt]197        if len(day) < min_obs_per_day:198            continue199        X = add_constant(standardize(day[features].values))200        try:201            coef_series.append(ols(X, day[target].values))202        except np.linalg.LinAlgError:203            continue204205    if len(coef_series) < min_periods:206        return None207208    coef_array = np.array(coef_series)209    avg = coef_array.mean(axis=0)210    se = coef_array.std(axis=0) / np.sqrt(len(coef_series))211    t = avg / se212    results = pd.DataFrame({213        'variable': ['const'] + features,214        'fm_coefficient': avg,215        'fm_std_error': se,216        'fm_t_stat': t,217        'fm_significant_5pct': np.abs(t) > 1.96,218        'fm_significant_1pct': np.abs(t) > 2.576,219    })220    return results, len(coef_series), sub['ticker'].nunique()221222223# --------------------------------------------------------------------------224# Time-series: Granger causality and bivariate VAR225# --------------------------------------------------------------------------226def granger_test(y: np.ndarray, x: np.ndarray, max_lag: int = 5):227    """F-test of the null that ``x`` does not Granger-cause ``y``."""228    n = len(y)229    if n < max_lag + 50:230        return None231232    Y = y[max_lag:]233    X_r = np.column_stack([y[max_lag - i - 1:n - i - 1] for i in range(max_lag)])234    X_r = add_constant(X_r)235    ssr_r = np.sum((Y - X_r @ ols(X_r, Y)) ** 2)236237    X_u = np.column_stack([X_r] + [x[max_lag - i - 1:n - i - 1] for i in range(max_lag)])238    ssr_u = np.sum((Y - X_u @ ols(X_u, Y)) ** 2)239240    n_eff = len(Y)241    f_stat = ((ssr_r - ssr_u) / max_lag) / (ssr_u / (n_eff - 2 * max_lag - 1))242243    from scipy import stats244    try:245        p_value = 1 - stats.f.cdf(f_stat, max_lag, n_eff - 2 * max_lag - 1)246    except Exception:247        p_value = np.nan248    return {'f_stat': f_stat, 'p_value': p_value, 'n_obs': n_eff}249250251def estimate_var(y1: np.ndarray, y2: np.ndarray, lags: int = 5,252                 irf_periods: int = 20):253    """Bivariate VAR(p) by equation-wise OLS, with companion-form IRFs."""254    n = len(y1)255    if n < lags + 50:256        return None257258    Y1, Y2 = y1[lags:], y2[lags:]259    X_lags = []260    for lag in range(1, lags + 1):261        X_lags.append(y1[lags - lag:n - lag])262        X_lags.append(y2[lags - lag:n - lag])263    X = np.column_stack([np.ones(len(Y1))] + X_lags)264265    coefs1 = ols(X, Y1)266    r2_1 = 1 - np.sum((Y1 - X @ coefs1) ** 2) / np.sum((Y1 - Y1.mean()) ** 2)267    coefs2 = ols(X, Y2)268    r2_2 = 1 - np.sum((Y2 - X @ coefs2) ** 2) / np.sum((Y2 - Y2.mean()) ** 2)269270    # Companion matrix271    A = np.zeros((2 * lags, 2 * lags))272    for lag in range(lags):273        A[0, 2 * lag] = coefs1[1 + 2 * lag]274        A[0, 2 * lag + 1] = coefs1[1 + 2 * lag + 1]275        A[1, 2 * lag] = coefs2[1 + 2 * lag]276        A[1, 2 * lag + 1] = coefs2[1 + 2 * lag + 1]277    for i in range(2, 2 * lags):278        A[i, i - 2] = 1.0279280    irf = np.zeros((irf_periods, 2, 2))281    for h in range(irf_periods):282        power = np.eye(2 * lags) if h == 0 else power @ A283        irf[h] = power[:2, :2]284285    return {'coefs_eq1': coefs1, 'coefs_eq2': coefs2,286            'r2_eq1': r2_1, 'r2_eq2': r2_2,287            'irf': irf, 'n_obs': len(Y1)}288