"""Spurious regression check for marketing time series (numpy only).

From the Vexpower MMM library: vexpower.com/mmm/marketing-econometrics

Part 1 reproduces the Granger & Newbold (1974) experiment: regress one random
walk on another, independent one, and count how often the slope looks
"significant" at 5%, in levels and in first differences.

Part 2 checks your own weekly data before you model it:
  * a Dickey-Fuller style unit-root test on every column (is it trending or
    wandering, rather than fluctuating around a stable level?),
  * the same regression in levels and in weekly changes, with R-squared and
    the Durbin-Watson statistic, so you can see whether a relationship
    survives differencing.

Usage:
    python spurious-regression-check.py                      # simulation + practice CSV
    python spurious-regression-check.py your.csv revenue     # your file, outcome column
The CSV needs a header row; non-numeric columns (like a date) are skipped.

This is a screening tool, not a replacement for a proper ADF/KPSS test from
statsmodels or a full MMM. Critical values are approximate.
"""
import csv
import sys
from pathlib import Path

import numpy as np

# Approximate 5% critical value for the Dickey-Fuller t-statistic with a
# constant (no trend), for samples of roughly 50-250 observations.
DF_CRIT_5 = -2.88


def ols(X, y):
    """OLS with intercept already in X. Returns coef, se, r2, residuals."""
    XtX = X.T @ X
    coef = np.linalg.solve(XtX, X.T @ y)
    resid = y - X @ coef
    n, p = X.shape
    s2 = resid @ resid / (n - p)
    se = np.sqrt(np.diag(np.linalg.inv(XtX)) * s2)
    r2 = 1 - (resid @ resid) / ((y - y.mean()) @ (y - y.mean()))
    return coef, se, r2, resid


def durbin_watson(e):
    return float(np.sum(np.diff(e) ** 2) / np.sum(e ** 2))


def simple_regression(x, y):
    X = np.column_stack([np.ones_like(x), x])
    coef, se, r2, e = ols(X, y)
    return {"slope": coef[1], "t": coef[1] / se[1], "r2": r2, "dw": durbin_watson(e)}


def dickey_fuller(y, lags=1):
    """Augmented Dickey-Fuller t-stat with a constant: dy_t = a + g*y_{t-1} + sum b_i dy_{t-i}."""
    y = np.asarray(y, float)
    dy = np.diff(y)
    rows = []
    target = []
    for t in range(lags, len(dy)):
        rows.append([1.0, y[t]] + [dy[t - i] for i in range(1, lags + 1)])
        target.append(dy[t])
    X, z = np.array(rows), np.array(target)
    coef, se, _, _ = ols(X, z)
    return coef[1] / se[1]


def simulate(n=50, reps=1000, seed=1974):
    rng = np.random.default_rng(seed)
    sig_levels = sig_diffs = 0
    r2s = []
    for _ in range(reps):
        a = np.cumsum(rng.normal(size=n))
        b = np.cumsum(rng.normal(size=n))
        lv = simple_regression(a, b)
        df = simple_regression(np.diff(a), np.diff(b))
        sig_levels += abs(lv["t"]) > 1.96
        sig_diffs += abs(df["t"]) > 1.96
        r2s.append(lv["r2"])
    print(f"Part 1: {reps} pairs of independent random walks, {n} periods each")
    print(f"  'Significant' at 5% in levels:      {sig_levels / reps:6.1%}  (average R2 {np.mean(r2s):.2f})")
    print(f"  'Significant' at 5% in differences: {sig_diffs / reps:6.1%}")
    print("  Granger & Newbold (1974) reported 76% and 8% with 100 replications.\n")


def load_numeric(path):
    with open(path, newline="") as f:
        rows = list(csv.DictReader(f))
    cols = {}
    for name in rows[0].keys():
        try:
            cols[name] = np.array([float(r[name]) for r in rows])
        except ValueError:
            continue
    return cols


def check_file(path, outcome):
    cols = load_numeric(path)
    if outcome not in cols:
        sys.exit(f"Outcome column '{outcome}' not found or not numeric. Columns: {', '.join(cols)}")
    print(f"Part 2: {Path(path).name}, {len(cols[outcome])} rows, outcome = {outcome}")
    print(f"  Unit-root screen (Dickey-Fuller t, 1 lag; below {DF_CRIT_5} suggests stationary):")
    for name, v in cols.items():
        if np.std(v) == 0 or set(np.unique(v)) <= {0.0, 1.0}:
            print(f"    {name:<16} skipped (constant or 0/1 flag)")
            continue
        t = dickey_fuller(v)
        verdict = "looks stationary" if t < DF_CRIT_5 else "may be trending or wandering: check levels vs changes"
        print(f"    {name:<16} t = {t:6.2f}   {verdict}")
    print("\n  Regressing the outcome on each input, in levels and in weekly changes:")
    print(f"    {'input':<16} {'levels R2':>9} {'levels t':>9} {'DW':>6}   {'changes R2':>10} {'changes t':>9} {'DW':>6}")
    y = cols[outcome]
    for name, x in cols.items():
        if name == outcome or np.std(x) == 0:
            continue
        lv = simple_regression(x, y)
        ch = simple_regression(np.diff(x), np.diff(y))
        flag = "  <- levels only" if abs(lv["t"]) > 1.96 and abs(ch["t"]) <= 1.96 else ""
        print(f"    {name:<16} {lv['r2']:9.2f} {lv['t']:9.1f} {lv['dw']:6.2f}   {ch['r2']:10.2f} {ch['t']:9.1f} {ch['dw']:6.2f}{flag}")
    print("\n  A relationship that is strong in levels, with a Durbin-Watson well below 2, and vanishes")
    print("  in changes deserves suspicion. It may still be real (effects can be slow, and adstock")
    print("  smooths spend), but it needs a reason beyond a shared trend before it goes in a model.")


if __name__ == "__main__":
    simulate()
    here = Path(__file__).resolve().parent
    path = sys.argv[1] if len(sys.argv) > 1 else str(here / "mmm-practice-synthetic.csv")
    outcome = sys.argv[2] if len(sys.argv) > 2 else "revenue"
    if Path(path).exists():
        check_file(path, outcome)
    else:
        print(f"No CSV at {path}; download mmm-practice-synthetic.csv from vexpower.com/mmm to run Part 2.")