Correlation Analysis

#!/usr/bin/env python3
"""Multi-asset correlation matrix analysis with hierarchical clustering.

Fetches historical price data for multiple crypto assets, computes correlation
matrices (Pearson and Spearman), performs hierarchical clustering to identify
asset groups, and reports diversification metrics.

Usage:
    python scripts/correlation_matrix.py --demo
    python scripts/correlation_matrix.py --coins bitcoin,ethereum,solana,cardano

Dependencies:
    uv pip install pandas numpy scipy httpx

Environment Variables:
    None required for demo mode.
"""

import argparse
import sys
import time
from typing import Optional

import httpx
import numpy as np
import pandas as pd
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import squareform


# ── Configuration ───────────────────────────────────────────────────
COINGECKO_BASE = "https://api.coingecko.com/api/v3"
DEFAULT_COINS = [
    "bitcoin", "ethereum", "solana", "cardano",
    "avalanche-2", "polkadot", "chainlink", "dogecoin",
]
DAYS = 180
REQUEST_DELAY = 1.5  # CoinGecko free tier rate limit


# ── Data Fetching ───────────────────────────────────────────────────
def fetch_prices(
    coin_id: str,
    days: int = DAYS,
    client: Optional[httpx.Client] = None,
) -> Optional[pd.Series]:
    """Fetch daily closing prices from CoinGecko.

    Args:
        coin_id: CoinGecko coin identifier (e.g., 'bitcoin').
        days: Number of days of history to fetch.
        client: Optional reusable httpx client.

    Returns:
        Series of daily prices indexed by date, or None on failure.
    """
    url = f"{COINGECKO_BASE}/coins/{coin_id}/market_chart"
    params = {"vs_currency": "usd", "days": days, "interval": "daily"}
    try:
        if client:
            resp = client.get(url, params=params, timeout=15)
        else:
            resp = httpx.get(url, params=params, timeout=15)
        resp.raise_for_status()
        data = resp.json()
        prices = data.get("prices", [])
        if not prices:
            print(f"  Warning: no price data for {coin_id}")
            return None
        series = pd.Series(
            {pd.Timestamp(ts, unit="ms").normalize(): price for ts, price in prices},
            name=coin_id,
        )
        return series
    except httpx.HTTPStatusError as e:
        print(f"  HTTP error fetching {coin_id}: {e.response.status_code}")
        return None
    except Exception as e:
        print(f"  Error fetching {coin_id}: {e}")
        return None


def generate_demo_data(n_assets: int = 8, n_days: int = 180) -> pd.DataFrame:
    """Generate synthetic correlated return data for demo mode.

    Creates assets with realistic crypto correlation structure:
    - A market factor driving all assets
    - Sector factors for subgroups
    - Idiosyncratic noise per asset

    Args:
        n_assets: Number of synthetic assets.
        n_days: Number of days of data.

    Returns:
        DataFrame of daily returns with asset columns.
    """
    rng = np.random.default_rng(42)
    names = ["BTC", "ETH", "SOL", "AVAX", "DOT", "LINK", "DOGE", "SHIB"][:n_assets]

    # Market factor (drives all assets)
    market = rng.normal(0.0005, 0.025, n_days)

    # Sector factors
    l1_factor = rng.normal(0, 0.015, n_days)  # L1 factor
    meme_factor = rng.normal(0, 0.02, n_days)  # Meme factor

    # Asset-specific loadings on factors
    loadings = {
        "BTC":  {"market": 0.85, "l1": 0.10, "meme": 0.00, "idio": 0.015},
        "ETH":  {"market": 0.90, "l1": 0.30, "meme": 0.00, "idio": 0.018},
        "SOL":  {"market": 0.80, "l1": 0.50, "meme": 0.05, "idio": 0.025},
        "AVAX": {"market": 0.75, "l1": 0.45, "meme": 0.00, "idio": 0.022},
        "DOT":  {"market": 0.70, "l1": 0.40, "meme": 0.00, "idio": 0.020},
        "LINK": {"market": 0.65, "l1": 0.20, "meme": 0.00, "idio": 0.020},
        "DOGE": {"market": 0.50, "l1": 0.05, "meme": 0.70, "idio": 0.035},
        "SHIB": {"market": 0.45, "l1": 0.00, "meme": 0.75, "idio": 0.040},
    }

    returns_data: dict[str, np.ndarray] = {}
    for name in names:
        l = loadings[name]
        asset_return = (
            l["market"] * market
            + l["l1"] * l1_factor
            + l["meme"] * meme_factor
            + rng.normal(0, l["idio"], n_days)
        )
        returns_data[name] = asset_return

    dates = pd.bdate_range(end=pd.Timestamp.today().normalize(), periods=n_days)
    return pd.DataFrame(returns_data, index=dates)


# ── Correlation Analysis ────────────────────────────────────────────
def compute_correlation_matrices(
    returns: pd.DataFrame,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Compute Pearson and Spearman correlation matrices.

    Args:
        returns: DataFrame of asset returns.

    Returns:
        Tuple of (pearson_corr, spearman_corr) DataFrames.
    """
    pearson = returns.corr(method="pearson")
    spearman = returns.corr(method="spearman")
    return pearson, spearman


def find_extreme_pairs(
    corr_matrix: pd.DataFrame, n: int = 5
) -> tuple[list[tuple[str, str, float]], list[tuple[str, str, float]]]:
    """Find the most and least correlated pairs.

    Args:
        corr_matrix: Correlation matrix DataFrame.
        n: Number of pairs to return.

    Returns:
        Tuple of (strongest_pairs, weakest_pairs), each a list of
        (asset_a, asset_b, correlation) tuples.
    """
    pairs: list[tuple[str, str, float]] = []
    cols = corr_matrix.columns.tolist()
    for i in range(len(cols)):
        for j in range(i + 1, len(cols)):
            pairs.append((cols[i], cols[j], corr_matrix.iloc[i, j]))

    pairs_sorted = sorted(pairs, key=lambda x: x[2], reverse=True)
    strongest = pairs_sorted[:n]
    weakest = pairs_sorted[-n:]
    return strongest, weakest


# ── Hierarchical Clustering ────────────────────────────────────────
def cluster_assets(
    corr_matrix: pd.DataFrame, max_clusters: int = 5
) -> tuple[np.ndarray, np.ndarray, list[int]]:
    """Perform hierarchical clustering on assets based on correlation.

    Args:
        corr_matrix: Correlation matrix DataFrame.
        max_clusters: Maximum number of clusters.

    Returns:
        Tuple of (linkage_matrix, distance_matrix, cluster_labels).
    """
    # Convert correlation to distance
    dist_matrix = np.sqrt(2.0 * (1.0 - corr_matrix.values))
    np.fill_diagonal(dist_matrix, 0.0)

    # Ensure symmetry (numerical precision)
    dist_matrix = (dist_matrix + dist_matrix.T) / 2.0

    condensed = squareform(dist_matrix)
    linkage_matrix = linkage(condensed, method="ward")

    # Determine optimal number of clusters (up to max_clusters)
    labels = fcluster(linkage_matrix, t=max_clusters, criterion="maxclust")
    return linkage_matrix, dist_matrix, labels.tolist()


def print_dendrogram_text(
    linkage_matrix: np.ndarray, labels: list[str]
) -> None:
    """Print a text-based representation of the clustering hierarchy.

    Args:
        linkage_matrix: Linkage matrix from scipy.
        labels: Asset names.
    """
    n = len(labels)
    clusters: dict[int, str] = {i: labels[i] for i in range(n)}

    print("\n  Clustering Hierarchy (Ward linkage):")
    print("  " + "-" * 55)
    for i, row in enumerate(linkage_matrix):
        left, right, dist, count = int(row[0]), int(row[1]), row[2], int(row[3])
        left_label = clusters.get(left, f"C{left}")
        right_label = clusters.get(right, f"C{right}")
        merged_label = f"({left_label} + {right_label})"
        clusters[n + i] = merged_label
        print(f"  d={dist:.3f}: {left_label}  <->  {right_label}")
    print()


# ── Diversification Metrics ────────────────────────────────────────
def diversification_metrics(corr_matrix: pd.DataFrame) -> dict[str, float]:
    """Compute portfolio diversification metrics from correlation matrix.

    Args:
        corr_matrix: Correlation matrix DataFrame.

    Returns:
        Dictionary with diversification metrics.
    """
    n = corr_matrix.shape[0]
    corr_values = corr_matrix.values

    # Average pairwise correlation (off-diagonal)
    mask = ~np.eye(n, dtype=bool)
    avg_corr = corr_values[mask].mean()

    # Effective N from average correlation
    if avg_corr < 1.0:
        effective_n_simple = 1.0 / (1.0 / n + (1.0 - 1.0 / n) * avg_corr)
    else:
        effective_n_simple = 1.0

    # Effective N from eigenvalue entropy
    eigenvalues = np.linalg.eigvalsh(corr_values)
    eigenvalues = eigenvalues[eigenvalues > 1e-10]
    proportions = eigenvalues / eigenvalues.sum()
    entropy = -np.sum(proportions * np.log(proportions))
    effective_n_eigen = np.exp(entropy)

    # First eigenvalue dominance (market factor strength)
    eigenvalues_sorted = np.sort(eigenvalues)[::-1]
    first_eigen_pct = eigenvalues_sorted[0] / eigenvalues_sorted.sum() * 100

    return {
        "n_assets": n,
        "avg_pairwise_corr": avg_corr,
        "effective_n_simple": effective_n_simple,
        "effective_n_eigenvalue": effective_n_eigen,
        "first_eigenvalue_pct": first_eigen_pct,
        "max_corr": corr_values[mask].max(),
        "min_corr": corr_values[mask].min(),
    }


# ── Eigenvalue Analysis ────────────────────────────────────────────
def eigenvalue_analysis(corr_matrix: pd.DataFrame) -> None:
    """Print eigenvalue decomposition of correlation matrix.

    Args:
        corr_matrix: Correlation matrix DataFrame.
    """
    eigenvalues, eigenvectors = np.linalg.eigh(corr_matrix.values)
    idx = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    total = eigenvalues.sum()
    cumulative = 0.0

    print("\n  Eigenvalue Decomposition:")
    print("  " + "-" * 55)
    print(f"  {'Factor':<10} {'Eigenvalue':>10} {'% Var':>8} {'Cum %':>8}")
    print("  " + "-" * 55)
    for i, ev in enumerate(eigenvalues):
        pct = ev / total * 100
        cumulative += pct
        label = "Market" if i == 0 else f"Factor {i}"
        print(f"  {label:<10} {ev:>10.3f} {pct:>7.1f}% {cumulative:>7.1f}%")
    print()


# ── Display Functions ───────────────────────────────────────────────
def print_correlation_heatmap(corr_matrix: pd.DataFrame, title: str) -> None:
    """Print a text-based correlation heatmap.

    Args:
        corr_matrix: Correlation matrix DataFrame.
        title: Title for the heatmap.
    """
    labels = corr_matrix.columns.tolist()
    max_label = max(len(l) for l in labels)

    print(f"\n  {title}")
    print("  " + "-" * (max_label + 2 + len(labels) * 7))

    # Header
    header = " " * (max_label + 2)
    for label in labels:
        header += f"{label:>7}"
    print(f"  {header}")

    # Rows
    for i, row_label in enumerate(labels):
        row_str = f"  {row_label:>{max_label}}  "
        for j in range(len(labels)):
            val = corr_matrix.iloc[i, j]
            if i == j:
                row_str += "   1.00"
            elif val >= 0.8:
                row_str += f"  {val:5.2f}"  # Strong positive
            elif val >= 0.5:
                row_str += f"  {val:5.2f}"
            elif val <= -0.3:
                row_str += f" {val:5.2f}"
            else:
                row_str += f"  {val:5.2f}"
        print(row_str)
    print()

    # Legend
    print("  Interpretation: >0.8 strong | 0.5-0.8 moderate | <0.5 weak")


def print_cluster_report(
    labels: list[str], clusters: list[int]
) -> None:
    """Print cluster membership report.

    Args:
        labels: Asset names.
        clusters: Cluster assignment for each asset.
    """
    print("\n  Asset Clusters:")
    print("  " + "-" * 40)
    cluster_map: dict[int, list[str]] = {}
    for asset, cluster in zip(labels, clusters):
        cluster_map.setdefault(cluster, []).append(asset)

    for cluster_id in sorted(cluster_map.keys()):
        members = cluster_map[cluster_id]
        print(f"  Cluster {cluster_id}: {', '.join(members)}")
    print()


# ── Main ────────────────────────────────────────────────────────────
def main() -> None:
    """Run multi-asset correlation analysis."""
    parser = argparse.ArgumentParser(
        description="Multi-asset correlation matrix analysis"
    )
    parser.add_argument(
        "--demo", action="store_true",
        help="Use synthetic demo data instead of fetching from API"
    )
    parser.add_argument(
        "--coins", type=str, default=None,
        help="Comma-separated CoinGecko coin IDs (e.g., bitcoin,ethereum,solana)"
    )
    parser.add_argument(
        "--days", type=int, default=DAYS,
        help=f"Days of history (default: {DAYS})"
    )
    args = parser.parse_args()

    print("=" * 65)
    print("  MULTI-ASSET CORRELATION ANALYSIS")
    print("=" * 65)

    if args.demo:
        print("\n  Mode: DEMO (synthetic data)")
        returns = generate_demo_data()
        asset_names = returns.columns.tolist()
    else:
        coins = args.coins.split(",") if args.coins else DEFAULT_COINS
        print(f"\n  Fetching {args.days}d prices for: {', '.join(coins)}")
        print("  (CoinGecko free tier — may be slow due to rate limits)")

        price_data: dict[str, pd.Series] = {}
        with httpx.Client() as client:
            for coin in coins:
                print(f"  Fetching {coin}...", end=" ", flush=True)
                series = fetch_prices(coin, days=args.days, client=client)
                if series is not None:
                    price_data[coin] = series
                    print("OK")
                else:
                    print("FAILED")
                time.sleep(REQUEST_DELAY)

        if len(price_data) < 3:
            print("\n  Error: need at least 3 assets. Try --demo mode.")
            sys.exit(1)

        prices = pd.DataFrame(price_data).dropna()
        returns = prices.pct_change().dropna()
        # Shorten names for display
        name_map = {c: c.upper()[:6] for c in returns.columns}
        returns = returns.rename(columns=name_map)
        asset_names = returns.columns.tolist()

    print(f"\n  Assets: {len(asset_names)}")
    print(f"  Observations: {len(returns)}")
    print(f"  Date range: {returns.index[0].strftime('%Y-%m-%d')} to "
          f"{returns.index[-1].strftime('%Y-%m-%d')}")

    # Compute correlation matrices
    pearson_corr, spearman_corr = compute_correlation_matrices(returns)

    # Display heatmaps
    print_correlation_heatmap(pearson_corr, "Pearson Correlation Matrix")
    print_correlation_heatmap(spearman_corr, "Spearman Rank Correlation Matrix")

    # Strongest and weakest pairs
    strongest, weakest = find_extreme_pairs(pearson_corr)
    print("\n  Strongest Correlated Pairs (Pearson):")
    print("  " + "-" * 40)
    for a, b, corr in strongest:
        print(f"  {a:>6} / {b:<6}  r = {corr:+.3f}")

    print("\n  Weakest Correlated Pairs (Pearson):")
    print("  " + "-" * 40)
    for a, b, corr in weakest:
        print(f"  {a:>6} / {b:<6}  r = {corr:+.3f}")

    # Eigenvalue analysis
    eigenvalue_analysis(pearson_corr)

    # Hierarchical clustering
    linkage_mat, dist_mat, clusters = cluster_assets(pearson_corr)
    print_dendrogram_text(linkage_mat, asset_names)
    print_cluster_report(asset_names, clusters)

    # Diversification metrics
    metrics = diversification_metrics(pearson_corr)
    print("\n  Diversification Metrics:")
    print("  " + "-" * 50)
    print(f"  Number of assets:             {metrics['n_assets']}")
    print(f"  Avg pairwise correlation:     {metrics['avg_pairwise_corr']:.3f}")
    print(f"  Effective N (simple):         {metrics['effective_n_simple']:.1f}")
    print(f"  Effective N (eigenvalue):     {metrics['effective_n_eigenvalue']:.1f}")
    print(f"  Market factor dominance:      {metrics['first_eigenvalue_pct']:.1f}%")
    print(f"  Max pairwise correlation:     {metrics['max_corr']:.3f}")
    print(f"  Min pairwise correlation:     {metrics['min_corr']:.3f}")

    # Assessment
    print("\n  Assessment:")
    print("  " + "-" * 50)
    avg = metrics["avg_pairwise_corr"]
    if avg > 0.7:
        print("  WARNING: High average correlation — portfolio is poorly diversified.")
        print("  Consider adding uncorrelated assets (stablecoins, different sectors).")
    elif avg > 0.4:
        print("  Moderate diversification. Some correlated clusters present.")
        print("  Consider reducing within-cluster allocations.")
    else:
        print("  Good diversification. Assets provide independent return streams.")

    mkt = metrics["first_eigenvalue_pct"]
    if mkt > 70:
        print(f"  Market factor explains {mkt:.0f}% of variance — all assets move together.")
    elif mkt > 50:
        print(f"  Market factor explains {mkt:.0f}% — significant common driver.")

    print("\n  Note: This is analysis output, not financial advice.")
    print("=" * 65)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Rolling correlation analysis with regime detection and tail dependence.

Computes rolling correlation between two assets at multiple windows, detects
correlation regime changes using z-score methodology, and estimates tail
dependence to assess crash co-movement.

Usage:
    python scripts/rolling_correlation.py --demo
    python scripts/rolling_correlation.py --coins bitcoin,ethereum --days 365

Dependencies:
    uv pip install pandas numpy httpx

Environment Variables:
    None required for demo mode.
"""

import argparse
import sys
import time
from typing import Optional

import httpx
import numpy as np
import pandas as pd


# ── Configuration ───────────────────────────────────────────────────
COINGECKO_BASE = "https://api.coingecko.com/api/v3"
WINDOWS = {"short": 20, "medium": 60, "long": 120}
TAIL_QUANTILE = 0.05
REQUEST_DELAY = 1.5


# ── Data Fetching ───────────────────────────────────────────────────
def fetch_prices(
    coin_id: str,
    days: int = 365,
    client: Optional[httpx.Client] = None,
) -> Optional[pd.Series]:
    """Fetch daily closing prices from CoinGecko.

    Args:
        coin_id: CoinGecko coin identifier.
        days: Number of days of history.
        client: Optional reusable httpx client.

    Returns:
        Series of daily prices indexed by date, or None on failure.
    """
    url = f"{COINGECKO_BASE}/coins/{coin_id}/market_chart"
    params = {"vs_currency": "usd", "days": days, "interval": "daily"}
    try:
        if client:
            resp = client.get(url, params=params, timeout=15)
        else:
            resp = httpx.get(url, params=params, timeout=15)
        resp.raise_for_status()
        data = resp.json()
        prices = data.get("prices", [])
        if not prices:
            return None
        return pd.Series(
            {pd.Timestamp(ts, unit="ms").normalize(): price for ts, price in prices},
            name=coin_id,
        )
    except Exception as e:
        print(f"  Error fetching {coin_id}: {e}")
        return None


def generate_demo_data(n_days: int = 365) -> tuple[pd.Series, pd.Series, str, str]:
    """Generate synthetic data showing correlation breakdown.

    Creates two assets that are normally correlated (~0.7) but experience
    a correlation spike mid-sample (simulating a market stress event)
    followed by a correlation breakdown (simulating structural change).

    Args:
        n_days: Number of days to generate.

    Returns:
        Tuple of (returns_a, returns_b, name_a, name_b).
    """
    rng = np.random.default_rng(123)
    dates = pd.bdate_range(end=pd.Timestamp.today().normalize(), periods=n_days)

    # Shared market factor
    market = rng.normal(0.0005, 0.02, n_days)

    # Asset A: consistent market exposure
    idio_a = rng.normal(0, 0.015, n_days)
    returns_a = 0.8 * market + idio_a

    # Asset B: changing correlation regime
    idio_b = rng.normal(0, 0.018, n_days)
    returns_b = np.zeros(n_days)

    # Phase 1 (days 0-120): normal correlation ~0.7
    phase1 = slice(0, 120)
    returns_b[phase1] = 0.7 * market[phase1] + idio_b[phase1]

    # Phase 2 (days 120-180): stress event — correlation spikes to ~0.95
    phase2 = slice(120, 180)
    stress = rng.normal(-0.02, 0.03, 60)  # negative drift during stress
    returns_b[phase2] = 0.95 * market[phase2] + 0.3 * stress + 0.5 * idio_b[phase2]

    # Phase 3 (days 180-280): recovery — correlation normalizes
    phase3 = slice(180, 280)
    returns_b[phase3] = 0.65 * market[phase3] + idio_b[phase3]

    # Phase 4 (days 280+): structural change — correlation breaks down
    phase4 = slice(280, n_days)
    sector_factor = rng.normal(0, 0.025, n_days)
    returns_b[phase4] = 0.3 * market[phase4] + 0.6 * sector_factor[phase4] + idio_b[phase4]

    series_a = pd.Series(returns_a, index=dates, name="Asset_A")
    series_b = pd.Series(returns_b, index=dates, name="Asset_B")
    return series_a, series_b, "Asset_A", "Asset_B"


# ── Rolling Correlation ────────────────────────────────────────────
def compute_rolling_correlations(
    returns_a: pd.Series,
    returns_b: pd.Series,
    windows: Optional[dict[str, int]] = None,
) -> pd.DataFrame:
    """Compute rolling correlation at multiple window sizes.

    Args:
        returns_a: First asset return series.
        returns_b: Second asset return series.
        windows: Dict of {label: window_size}. Defaults to WINDOWS.

    Returns:
        DataFrame with rolling correlation columns.
    """
    if windows is None:
        windows = WINDOWS

    result = pd.DataFrame(index=returns_a.index)
    for label, w in windows.items():
        result[f"corr_{label}"] = returns_a.rolling(w).corr(returns_b)
    return result


def ewma_correlation(
    x: pd.Series, y: pd.Series, span: int = 60
) -> pd.Series:
    """Compute EWMA correlation between two return series.

    Args:
        x: First return series.
        y: Second return series.
        span: EWMA span parameter.

    Returns:
        Series of EWMA correlation values.
    """
    mean_x = x.ewm(span=span).mean()
    mean_y = y.ewm(span=span).mean()
    cov = (x * y).ewm(span=span).mean() - mean_x * mean_y
    var_x = x.pow(2).ewm(span=span).mean() - mean_x.pow(2)
    var_y = y.pow(2).ewm(span=span).mean() - mean_y.pow(2)
    denom = var_x.pow(0.5) * var_y.pow(0.5)
    return cov / denom.replace(0, np.nan)


# ── Regime Detection ───────────────────────────────────────────────
def detect_correlation_regimes(
    rolling_corr: pd.Series,
    lookback: int = 120,
    z_threshold: float = 2.0,
) -> pd.DataFrame:
    """Detect correlation regime changes using z-score method.

    Args:
        rolling_corr: Series of rolling correlation values.
        lookback: Window for computing mean/std of correlation.
        z_threshold: Z-score threshold for regime flag.

    Returns:
        DataFrame with correlation, z-score, and regime columns.
    """
    clean = rolling_corr.dropna()
    mean = clean.rolling(lookback, min_periods=30).mean()
    std = clean.rolling(lookback, min_periods=30).std()
    zscore = (clean - mean) / std.replace(0, np.nan)

    regime = pd.Series("normal", index=clean.index)
    regime[zscore > z_threshold] = "HIGH_CORR"
    regime[zscore < -z_threshold] = "LOW_CORR"

    return pd.DataFrame({
        "correlation": clean,
        "rolling_mean": mean,
        "zscore": zscore,
        "regime": regime,
    })


def find_regime_changes(regime_df: pd.DataFrame) -> list[tuple[str, str, str]]:
    """Find dates where correlation regime changes.

    Args:
        regime_df: DataFrame from detect_correlation_regimes.

    Returns:
        List of (date_str, from_regime, to_regime) tuples.
    """
    changes: list[tuple[str, str, str]] = []
    prev_regime = "normal"
    for date, row in regime_df.iterrows():
        current = row["regime"]
        if current != prev_regime:
            changes.append((
                date.strftime("%Y-%m-%d"),  # type: ignore[union-attr]
                prev_regime,
                current,
            ))
        prev_regime = current
    return changes


# ── Tail Dependence ─────────────────────────────────────────────────
def compute_tail_dependence(
    returns_a: pd.Series,
    returns_b: pd.Series,
    quantile: float = TAIL_QUANTILE,
) -> dict[str, float]:
    """Estimate tail dependence coefficients.

    Args:
        returns_a: First asset returns.
        returns_b: Second asset returns.
        quantile: Quantile threshold for extreme events (default 5%).

    Returns:
        Dict with lower_tail, upper_tail, and joint_extreme_count.
    """
    n = len(returns_a)

    # Lower tail: crash together
    thresh_a_low = returns_a.quantile(quantile)
    thresh_b_low = returns_b.quantile(quantile)
    a_extreme_low = returns_a < thresh_a_low
    b_extreme_low = returns_b < thresh_b_low
    joint_low = (a_extreme_low & b_extreme_low).sum()
    marginal_low = a_extreme_low.sum()
    lower_tail = joint_low / marginal_low if marginal_low > 0 else 0.0

    # Upper tail: rally together
    thresh_a_high = returns_a.quantile(1 - quantile)
    thresh_b_high = returns_b.quantile(1 - quantile)
    a_extreme_high = returns_a > thresh_a_high
    b_extreme_high = returns_b > thresh_b_high
    joint_high = (a_extreme_high & b_extreme_high).sum()
    marginal_high = a_extreme_high.sum()
    upper_tail = joint_high / marginal_high if marginal_high > 0 else 0.0

    return {
        "lower_tail_dependence": lower_tail,
        "upper_tail_dependence": upper_tail,
        "joint_crash_events": int(joint_low),
        "joint_rally_events": int(joint_high),
        "total_observations": n,
        "extreme_threshold_pct": quantile * 100,
    }


# ── Summary Statistics ──────────────────────────────────────────────
def correlation_summary(
    returns_a: pd.Series,
    returns_b: pd.Series,
    name_a: str,
    name_b: str,
) -> dict[str, float]:
    """Compute full-sample correlation statistics.

    Args:
        returns_a: First asset returns.
        returns_b: Second asset returns.
        name_a: Name of first asset.
        name_b: Name of second asset.

    Returns:
        Dict of correlation statistics.
    """
    n = len(returns_a)
    pearson = returns_a.corr(returns_b, method="pearson")
    spearman = returns_a.corr(returns_b, method="spearman")
    kendall = returns_a.corr(returns_b, method="kendall")

    # Statistical significance (Pearson)
    if abs(pearson) < 1.0 and n > 2:
        t_stat = pearson * np.sqrt(n - 2) / np.sqrt(1 - pearson ** 2)
    else:
        t_stat = np.inf

    # Fisher z-transform confidence interval
    if abs(pearson) < 1.0 and n > 3:
        z = np.arctanh(pearson)
        se = 1.0 / np.sqrt(n - 3)
        ci_low = np.tanh(z - 1.96 * se)
        ci_high = np.tanh(z + 1.96 * se)
    else:
        ci_low, ci_high = pearson, pearson

    return {
        "pair": f"{name_a} / {name_b}",
        "n_obs": n,
        "pearson": pearson,
        "spearman": spearman,
        "kendall": kendall,
        "t_statistic": t_stat,
        "ci_95_low": ci_low,
        "ci_95_high": ci_high,
    }


# ── Display ─────────────────────────────────────────────────────────
def print_rolling_summary(
    rolling_df: pd.DataFrame, ewma_corr: pd.Series
) -> None:
    """Print summary of rolling correlation values.

    Args:
        rolling_df: DataFrame from compute_rolling_correlations.
        ewma_corr: Series of EWMA correlation values.
    """
    print("\n  Rolling Correlation Summary:")
    print("  " + "-" * 55)
    print(f"  {'Window':<12} {'Current':>8} {'Mean':>8} {'Min':>8} {'Max':>8} {'Std':>8}")
    print("  " + "-" * 55)

    for col in rolling_df.columns:
        s = rolling_df[col].dropna()
        if len(s) == 0:
            continue
        label = col.replace("corr_", "").capitalize()
        print(f"  {label:<12} {s.iloc[-1]:>8.3f} {s.mean():>8.3f} "
              f"{s.min():>8.3f} {s.max():>8.3f} {s.std():>8.3f}")

    ewma_clean = ewma_corr.dropna()
    if len(ewma_clean) > 0:
        print(f"  {'EWMA(60)':<12} {ewma_clean.iloc[-1]:>8.3f} "
              f"{ewma_clean.mean():>8.3f} {ewma_clean.min():>8.3f} "
              f"{ewma_clean.max():>8.3f} {ewma_clean.std():>8.3f}")
    print()


def print_regime_timeline(
    regime_df: pd.DataFrame, last_n: int = 20
) -> None:
    """Print recent correlation regime timeline.

    Args:
        regime_df: DataFrame from detect_correlation_regimes.
        last_n: Number of recent observations to show.
    """
    recent = regime_df.tail(last_n)
    print(f"\n  Recent Correlation Regime (last {last_n} observations):")
    print("  " + "-" * 55)
    print(f"  {'Date':<12} {'Corr':>7} {'Z-Score':>8} {'Regime':<12}")
    print("  " + "-" * 55)
    for date, row in recent.iterrows():
        date_str = date.strftime("%Y-%m-%d")  # type: ignore[union-attr]
        corr_val = row["correlation"]
        z_val = row["zscore"]
        regime = row["regime"]
        z_str = f"{z_val:>8.2f}" if not np.isnan(z_val) else "     N/A"
        marker = " !" if regime != "normal" else "  "
        print(f"  {date_str:<12} {corr_val:>7.3f} {z_str} {regime:<12}{marker}")
    print()


def print_correlation_assessment(
    summary: dict[str, float],
    tail_dep: dict[str, float],
    regime_df: pd.DataFrame,
) -> None:
    """Print overall assessment of the correlation relationship.

    Args:
        summary: From correlation_summary.
        tail_dep: From compute_tail_dependence.
        regime_df: From detect_correlation_regimes.
    """
    print("\n  Overall Assessment:")
    print("  " + "-" * 55)

    pearson = summary["pearson"]
    if abs(pearson) > 0.8:
        strength = "STRONG"
    elif abs(pearson) > 0.5:
        strength = "MODERATE"
    elif abs(pearson) > 0.3:
        strength = "WEAK"
    else:
        strength = "NEGLIGIBLE"

    direction = "positive" if pearson > 0 else "negative"
    print(f"  Correlation strength: {strength} {direction} (r={pearson:.3f})")
    print(f"  95% CI: [{summary['ci_95_low']:.3f}, {summary['ci_95_high']:.3f}]")

    # Spearman vs Pearson divergence
    diff = abs(summary["spearman"] - summary["pearson"])
    if diff > 0.1:
        print(f"  Spearman-Pearson gap: {diff:.3f} — suggests non-linear relationship")

    # Tail dependence assessment
    lower = tail_dep["lower_tail_dependence"]
    if lower > pearson and pearson > 0:
        print(f"  WARNING: Tail dependence ({lower:.2f}) > normal correlation ({pearson:.2f})")
        print(f"  Assets crash together MORE than normal correlation suggests.")
    elif lower > 0.5:
        print(f"  High tail dependence ({lower:.2f}) — significant crash co-movement.")

    # Regime stability
    regime_counts = regime_df["regime"].value_counts()
    pct_abnormal = 1.0 - regime_counts.get("normal", 0) / len(regime_df)
    if pct_abnormal > 0.2:
        print(f"  Unstable correlation: {pct_abnormal:.0%} of time in abnormal regime.")
    else:
        print(f"  Correlation regime is relatively stable ({pct_abnormal:.0%} abnormal).")

    # Current regime
    current_regime = regime_df["regime"].iloc[-1]
    current_z = regime_df["zscore"].iloc[-1]
    if current_regime != "normal":
        print(f"  ALERT: Currently in {current_regime} regime (z={current_z:.2f}).")

    print("\n  Note: This is analysis output, not financial advice.")


# ── Main ────────────────────────────────────────────────────────────
def main() -> None:
    """Run rolling correlation analysis."""
    parser = argparse.ArgumentParser(
        description="Rolling correlation with regime detection and tail dependence"
    )
    parser.add_argument(
        "--demo", action="store_true",
        help="Use synthetic demo data with correlation regime changes"
    )
    parser.add_argument(
        "--coins", type=str, default="bitcoin,ethereum",
        help="Two comma-separated CoinGecko IDs (default: bitcoin,ethereum)"
    )
    parser.add_argument(
        "--days", type=int, default=365,
        help="Days of history (default: 365)"
    )
    args = parser.parse_args()

    print("=" * 60)
    print("  ROLLING CORRELATION ANALYSIS")
    print("=" * 60)

    if args.demo:
        print("\n  Mode: DEMO (synthetic data with regime changes)")
        returns_a, returns_b, name_a, name_b = generate_demo_data()
    else:
        coins = args.coins.split(",")
        if len(coins) != 2:
            print("  Error: provide exactly 2 coin IDs (e.g., --coins bitcoin,ethereum)")
            sys.exit(1)

        name_a, name_b = coins[0], coins[1]
        print(f"\n  Fetching {args.days}d data for {name_a} and {name_b}...")

        with httpx.Client() as client:
            prices_a = fetch_prices(name_a, days=args.days, client=client)
            time.sleep(REQUEST_DELAY)
            prices_b = fetch_prices(name_b, days=args.days, client=client)

        if prices_a is None or prices_b is None:
            print("  Error: could not fetch price data. Try --demo mode.")
            sys.exit(1)

        prices = pd.DataFrame({"a": prices_a, "b": prices_b}).dropna()
        if len(prices) < 30:
            print(f"  Error: only {len(prices)} data points — need at least 30.")
            sys.exit(1)

        returns_a = prices["a"].pct_change().dropna()
        returns_b = prices["b"].pct_change().dropna()
        # Align
        common_idx = returns_a.index.intersection(returns_b.index)
        returns_a = returns_a.loc[common_idx]
        returns_b = returns_b.loc[common_idx]

    print(f"\n  Pair: {name_a} / {name_b}")
    print(f"  Observations: {len(returns_a)}")
    print(f"  Date range: {returns_a.index[0].strftime('%Y-%m-%d')} to "
          f"{returns_a.index[-1].strftime('%Y-%m-%d')}")

    # Full-sample correlation
    summary = correlation_summary(returns_a, returns_b, name_a, name_b)
    print("\n  Full-Sample Correlation:")
    print("  " + "-" * 45)
    print(f"  Pearson:    {summary['pearson']:>+.4f}")
    print(f"  Spearman:   {summary['spearman']:>+.4f}")
    print(f"  Kendall:    {summary['kendall']:>+.4f}")
    print(f"  t-statistic: {summary['t_statistic']:.2f}")
    print(f"  95% CI:     [{summary['ci_95_low']:.3f}, {summary['ci_95_high']:.3f}]")

    # Rolling correlation
    rolling_df = compute_rolling_correlations(returns_a, returns_b)
    ewma_corr = ewma_correlation(returns_a, returns_b, span=60)
    print_rolling_summary(rolling_df, ewma_corr)

    # Regime detection (using medium window correlation)
    medium_corr = rolling_df["corr_medium"].dropna()
    if len(medium_corr) > 30:
        regime_df = detect_correlation_regimes(medium_corr)
        print_regime_timeline(regime_df)

        # Regime changes
        changes = find_regime_changes(regime_df)
        if changes:
            print("  Regime Change Events:")
            print("  " + "-" * 45)
            for date_str, from_r, to_r in changes:
                print(f"  {date_str}: {from_r} -> {to_r}")
            print()
    else:
        regime_df = pd.DataFrame({"regime": ["normal"], "zscore": [0.0]})
        print("  Insufficient data for regime detection (need >90 observations).\n")

    # Tail dependence
    tail_dep = compute_tail_dependence(returns_a, returns_b)
    print("\n  Tail Dependence Analysis:")
    print("  " + "-" * 50)
    print(f"  Lower tail (crash together):  {tail_dep['lower_tail_dependence']:.3f}")
    print(f"  Upper tail (rally together):  {tail_dep['upper_tail_dependence']:.3f}")
    print(f"  Joint crash events:           {tail_dep['joint_crash_events']}")
    print(f"  Joint rally events:           {tail_dep['joint_rally_events']}")
    print(f"  Extreme threshold:            {tail_dep['extreme_threshold_pct']:.1f}th percentile")

    # Comparison: tail dependence vs normal correlation
    if tail_dep["lower_tail_dependence"] > abs(summary["pearson"]):
        print("  >> Tail dependence EXCEEDS normal correlation — crashes are more correlated")
        print("     than average co-movement suggests. Diversification benefit overstated.")
    else:
        print("  >> Tail dependence within expected range given normal correlation.")

    # Overall assessment
    print_correlation_assessment(summary, tail_dep, regime_df)
    print("=" * 60)


if __name__ == "__main__":
    main()