Regression Models for Sports Betting: From Linear to Logistic to Ridge

"""
Regression models for sports betting agents.
Covers: OLS, Logistic, Ridge, Poisson regression with
full training, evaluation, and prediction pipelines.
"""

import numpy as np
import pandas as pd
from sklearn.linear_model import (
    LinearRegression,
    LogisticRegressionCV,
    Ridge,
    Lasso,
    ElasticNet,
    PoissonRegressor,
)
from sklearn.model_selection import cross_val_score, TimeSeriesSplit
from sklearn.metrics import (
    brier_score_loss,
    log_loss,
    mean_squared_error,
    r2_score,
)
from sklearn.preprocessing import StandardScaler
from dataclasses import dataclass
from typing import Optional


@dataclass
class RegressionPrediction:
    """Output from a sports regression model."""
    predicted_value: float       # point spread, win prob, or expected goals
    model_type: str              # "ols", "logistic", "ridge", "poisson"
    confidence_interval: tuple[float, float]
    features_used: list[str]
    r2_or_accuracy: float


def compute_vif(X: pd.DataFrame) -> pd.DataFrame:
    """
    Compute Variance Inflation Factor for each feature.
    VIF > 5 indicates concerning multicollinearity.
    VIF > 10 indicates severe multicollinearity.
    """
    from sklearn.linear_model import LinearRegression

    vif_data = []
    for i, col in enumerate(X.columns):
        y_vif = X[col].values
        X_vif = X.drop(columns=[col]).values

        reg = LinearRegression().fit(X_vif, y_vif)
        r2 = reg.score(X_vif, y_vif)

        vif = 1.0 / (1.0 - r2) if r2 < 1.0 else np.inf
        vif_data.append({"feature": col, "VIF": round(vif, 2)})

    return pd.DataFrame(vif_data).sort_values("VIF", ascending=False)


def build_nfl_logistic_model(
    df: pd.DataFrame,
    features: list[str],
    target: str = "home_win",
    regularization: str = "elasticnet",
    n_splits: int = 5,
) -> dict:
    """
    Train a logistic regression model for NFL win probability.

    Args:
        df: DataFrame with game-level features and outcomes.
            Must have columns matching `features` plus `target`.
        features: List of feature column names.
        target: Binary target column (1 = home win, 0 = away win).
        regularization: One of 'l1', 'l2', 'elasticnet'.
        n_splits: Number of cross-validation folds.

    Returns:
        Dict with model, scaler, metrics, and coefficients.
    """
    X = df[features].copy()
    y = df[target].values

    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # TimeSeriesSplit respects temporal ordering — critical for sports
    tscv = TimeSeriesSplit(n_splits=n_splits)

    if regularization == "elasticnet":
        solver = "saga"
        penalty = "elasticnet"
        l1_ratios = [0.1, 0.3, 0.5, 0.7, 0.9]
    elif regularization == "l1":
        solver = "saga"
        penalty = "l1"
        l1_ratios = None
    else:
        solver = "lbfgs"
        penalty = "l2"
        l1_ratios = None

    model = LogisticRegressionCV(
        penalty=penalty,
        solver=solver,
        cv=tscv,
        Cs=20,
        l1_ratios=l1_ratios,
        max_iter=5000,
        scoring="neg_log_loss",
        random_state=42,
    )
    model.fit(X_scaled, y)

    # Cross-validated metrics
    cv_logloss = -cross_val_score(
        model, X_scaled, y, cv=tscv, scoring="neg_log_loss"
    )
    cv_brier = -cross_val_score(
        model, X_scaled, y, cv=tscv, scoring="neg_brier_score"
    )
    cv_accuracy = cross_val_score(
        model, X_scaled, y, cv=tscv, scoring="accuracy"
    )

    # Coefficient table
    coef_df = pd.DataFrame({
        "feature": features,
        "coefficient": model.coef_[0],
        "odds_ratio": np.exp(model.coef_[0]),
    }).sort_values("coefficient", ascending=False)

    return {
        "model": model,
        "scaler": scaler,
        "best_C": model.C_[0],
        "coefficients": coef_df,
        "cv_logloss_mean": cv_logloss.mean(),
        "cv_logloss_std": cv_logloss.std(),
        "cv_brier_mean": cv_brier.mean(),
        "cv_accuracy_mean": cv_accuracy.mean(),
    }


def predict_win_probability(
    model_dict: dict,
    game_features: dict[str, float],
) -> RegressionPrediction:
    """
    Predict win probability for a single game.

    Args:
        model_dict: Output from build_nfl_logistic_model.
        game_features: Dict of feature_name -> value for one game.

    Returns:
        RegressionPrediction with win probability and metadata.
    """
    model = model_dict["model"]
    scaler = model_dict["scaler"]
    features = model_dict["coefficients"]["feature"].tolist()

    X = np.array([[game_features[f] for f in features]])
    X_scaled = scaler.transform(X)

    prob = model.predict_proba(X_scaled)[0, 1]

    # Bootstrap-style confidence interval approximation
    # using coefficient standard errors
    z = model.decision_function(X_scaled)[0]
    se_approx = 0.15  # typical SE for NFL logistic models
    ci_low = 1 / (1 + np.exp(-(z - 1.96 * se_approx)))
    ci_high = 1 / (1 + np.exp(-(z + 1.96 * se_approx)))

    return RegressionPrediction(
        predicted_value=prob,
        model_type="logistic",
        confidence_interval=(round(ci_low, 4), round(ci_high, 4)),
        features_used=features,
        r2_or_accuracy=model_dict["cv_accuracy_mean"],
    )


def build_ridge_spread_model(
    df: pd.DataFrame,
    features: list[str],
    target: str = "point_diff",
    alphas: Optional[list[float]] = None,
) -> dict:
    """
    Train a Ridge regression model for point spread prediction.

    Args:
        df: DataFrame with game features and point differential.
        features: Feature column names.
        target: Continuous target (home points - away points).
        alphas: List of regularization strengths to try.

    Returns:
        Dict with model, scaler, and evaluation metrics.
    """
    from sklearn.linear_model import RidgeCV

    if alphas is None:
        alphas = np.logspace(-2, 4, 50)

    X = df[features].copy()
    y = df[target].values

    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    tscv = TimeSeriesSplit(n_splits=5)

    model = RidgeCV(alphas=alphas, cv=tscv, scoring="neg_mean_squared_error")
    model.fit(X_scaled, y)

    y_pred = model.predict(X_scaled)
    r2 = r2_score(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))

    cv_rmse = np.sqrt(-cross_val_score(
        model, X_scaled, y, cv=tscv, scoring="neg_mean_squared_error"
    ))

    coef_df = pd.DataFrame({
        "feature": features,
        "coefficient": model.coef_,
    }).sort_values("coefficient", ascending=False)

    return {
        "model": model,
        "scaler": scaler,
        "best_alpha": model.alpha_,
        "r2": r2,
        "rmse": rmse,
        "cv_rmse_mean": cv_rmse.mean(),
        "cv_rmse_std": cv_rmse.std(),
        "coefficients": coef_df,
    }


def build_poisson_model(
    df: pd.DataFrame,
    features: list[str],
    target: str = "goals_scored",
) -> dict:
    """
    Train a Poisson regression model for count outcomes
    (goals, runs, touchdowns).

    Args:
        df: DataFrame with features and count target.
        features: Feature column names.
        target: Count target column.

    Returns:
        Dict with model, scaler, and metrics.
    """
    X = df[features].copy()
    y = df[target].values

    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    model = PoissonRegressor(alpha=1.0, max_iter=1000)
    model.fit(X_scaled, y)

    y_pred = model.predict(X_scaled)

    # Overdispersion test: compare variance to mean
    residuals = y - y_pred
    pearson_chi2 = np.sum(residuals**2 / y_pred)
    dispersion = pearson_chi2 / (len(y) - len(features) - 1)

    coef_df = pd.DataFrame({
        "feature": features,
        "coefficient": model.coef_,
        "rate_ratio": np.exp(model.coef_),
    }).sort_values("coefficient", ascending=False)

    return {
        "model": model,
        "scaler": scaler,
        "coefficients": coef_df,
        "dispersion": dispersion,
        "overdispersed": dispersion > 1.5,
        "mean_predicted": y_pred.mean(),
        "mean_actual": y.mean(),
    }


def calibration_table(
    y_true: np.ndarray,
    y_prob: np.ndarray,
    n_bins: int = 10,
) -> pd.DataFrame:
    """
    Build a reliability diagram table for model calibration assessment.

    A well-calibrated model has observed frequency ≈ predicted probability
    in each bin. Use this to diagnose systematic over/under-confidence.

    Args:
        y_true: Binary outcomes (0 or 1).
        y_prob: Predicted probabilities.
        n_bins: Number of bins.

    Returns:
        DataFrame with bin_center, mean_predicted, observed_frequency, count.
    """
    bins = np.linspace(0, 1, n_bins + 1)
    table = []

    for i in range(n_bins):
        mask = (y_prob >= bins[i]) & (y_prob < bins[i + 1])
        if mask.sum() == 0:
            continue

        table.append({
            "bin_center": (bins[i] + bins[i + 1]) / 2,
            "mean_predicted": y_prob[mask].mean(),
            "observed_frequency": y_true[mask].mean(),
            "count": mask.sum(),
        })

    return pd.DataFrame(table)


def calculate_clv(
    model_prob: float,
    opening_implied: float,
    closing_implied: float,
) -> dict:
    """
    Calculate Closing Line Value.

    CLV measures whether your model's edge was real by comparing
    against the closing line (the sharpest estimate of true probability).

    Args:
        model_prob: Your model's predicted probability.
        opening_implied: Sportsbook opening implied probability.
        closing_implied: Sportsbook closing implied probability.

    Returns:
        Dict with CLV metrics.
    """
    # Model edge vs opening
    edge_vs_open = model_prob - opening_implied

    # Model edge vs closing
    edge_vs_close = model_prob - closing_implied

    # Did model beat the closing line?
    beat_close = model_prob > closing_implied

    # CLV in percentage points
    clv = closing_implied - opening_implied

    return {
        "model_prob": model_prob,
        "opening_implied": opening_implied,
        "closing_implied": closing_implied,
        "edge_vs_open": round(edge_vs_open, 4),
        "edge_vs_close": round(edge_vs_close, 4),
        "clv_points": round(clv, 4),
        "beat_closing_line": beat_close,
    }


# ------------------------------------------------------------------
# Example: End-to-end NFL model pipeline
# ------------------------------------------------------------------
if __name__ == "__main__":
    np.random.seed(42)

    # Simulate 500 NFL games with realistic feature distributions
    n_games = 500
    df = pd.DataFrame({
        "home_off_dvoa": np.random.normal(0, 15, n_games),
        "away_off_dvoa": np.random.normal(0, 15, n_games),
        "home_def_dvoa": np.random.normal(0, 12, n_games),
        "away_def_dvoa": np.random.normal(0, 12, n_games),
        "home_adv": np.ones(n_games),
        "rest_diff": np.random.choice([-3, 0, 3, 7], n_games, p=[0.1, 0.6, 0.2, 0.1]),
    })

    # Generate outcomes from a known logistic model
    true_z = (
        0.35 * df["home_adv"]
        + 0.04 * df["home_off_dvoa"]
        - 0.035 * df["away_off_dvoa"]
        - 0.03 * df["home_def_dvoa"]
        + 0.03 * df["away_def_dvoa"]
        + 0.08 * df["rest_diff"]
    )
    true_prob = 1 / (1 + np.exp(-true_z))
    df["home_win"] = (np.random.rand(n_games) < true_prob).astype(int)

    # Point differential (correlated with win outcome)
    df["point_diff"] = (
        true_z * 7
        + np.random.normal(0, 10, n_games)
    ).round(0)

    features = [
        "home_off_dvoa", "away_off_dvoa",
        "home_def_dvoa", "away_def_dvoa",
        "home_adv", "rest_diff",
    ]

    # 1. Check multicollinearity
    print("=== VIF Analysis ===")
    vif = compute_vif(df[features])
    print(vif.to_string(index=False))

    # 2. Train logistic model
    print("\n=== Logistic Regression (Elastic Net) ===")
    result = build_nfl_logistic_model(df, features, regularization="elasticnet")
    print(f"Best C (inverse regularization): {result['best_C']:.4f}")
    print(f"CV Log-Loss: {result['cv_logloss_mean']:.4f} +/- {result['cv_logloss_std']:.4f}")
    print(f"CV Brier: {result['cv_brier_mean']:.4f}")
    print(f"CV Accuracy: {result['cv_accuracy_mean']:.1%}")
    print("\nCoefficients:")
    print(result["coefficients"].to_string(index=False))

    # 3. Predict a specific game
    game = {
        "home_off_dvoa": 18.2,
        "away_off_dvoa": 12.5,
        "home_def_dvoa": -8.3,
        "away_def_dvoa": -12.4,
        "home_adv": 1.0,
        "rest_diff": 0.0,
    }
    pred = predict_win_probability(result, game)
    print(f"\nPredicted win prob: {pred.predicted_value:.1%}")
    print(f"95% CI: ({pred.confidence_interval[0]:.1%}, {pred.confidence_interval[1]:.1%})")

    # 4. Calibration check
    print("\n=== Calibration Table ===")
    model = result["model"]
    scaler = result["scaler"]
    X_scaled = scaler.transform(df[features])
    probs = model.predict_proba(X_scaled)[:, 1]
    cal = calibration_table(df["home_win"].values, probs)
    print(cal.to_string(index=False))

    # 5. CLV calculation
    print("\n=== CLV Example ===")
    clv = calculate_clv(
        model_prob=0.815,
        opening_implied=0.75,
        closing_implied=0.78,
    )
    for k, v in clv.items():
        print(f"  {k}: {v}")

    # 6. Ridge spread model
    print("\n=== Ridge Regression (Point Spread) ===")
    spread_result = build_ridge_spread_model(df, features, target="point_diff")
    print(f"Best alpha: {spread_result['best_alpha']:.2f}")
    print(f"R²: {spread_result['r2']:.3f}")
    print(f"RMSE: {spread_result['rmse']:.1f}")
    print(f"CV RMSE: {spread_result['cv_rmse_mean']:.1f} +/- {spread_result['cv_rmse_std']:.1f}")