# Skill: Granger Causality — Sunspot Number → Global Temperature Anomaly

## Overview
This skill tests whether the international sunspot number (SILSO v2.0) Granger-causes global land-ocean temperature anomalies (NOAA NCEI) over the period 1850–2024. The analysis includes cross-correlation, bidirectional Granger tests at multiple lag orders, partial correlation on AR residuals, a permutation test, sliding-window stability analysis, and detrended robustness checks.

## Data Requirements
1. **Sunspot data**: SILSO v2.0 monthly sunspot number CSV (semicolon-delimited; columns: year;month;decimal_year;ssn;...). Download from https://www.sidc.be/SILSO/INFO/snmtotcsv.php
2. **Temperature data**: NOAA NCEI global land-ocean monthly temperature anomaly. Download from https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series/globe/land_ocean/1/0/1850-2024/data.csv

## Dependencies
- Pure Python standard library only (math, random, csv, os)
- random.seed(42) for reproducibility

## Analysis Script

```python
#!/usr/bin/env python3
"""
Granger Causality Analysis: Solar Activity → Global Temperature Anomaly
Pure Python stdlib implementation.
"""
import math
import random
import csv
import os

random.seed(42)

# ─────────────────────────────────────────────────────────
# 1.  LINEAR ALGEBRA HELPERS  (small-matrix OLS)
# ─────────────────────────────────────────────────────────

def mat_zeros(r, c):
    return [[0.0]*c for _ in range(r)]

def mat_mul(A, B):
    rA, cA = len(A), len(A[0])
    rB, cB = len(B), len(B[0])
    assert cA == rB
    C = mat_zeros(rA, cB)
    for i in range(rA):
        for k in range(cA):
            aik = A[i][k]
            for j in range(cB):
                C[i][j] += aik * B[k][j]
    return C

def mat_T(A):
    r, c = len(A), len(A[0])
    return [[A[i][j] for i in range(r)] for j in range(c)]

def mat_inv(M):
    """Gauss-Jordan inverse for small matrices."""
    n = len(M)
    aug = [M[i][:] + [1.0 if j == i else 0.0 for j in range(n)] for i in range(n)]
    for col in range(n):
        # partial pivot
        mx = abs(aug[col][col]); mi = col
        for row in range(col+1, n):
            if abs(aug[row][col]) > mx:
                mx = abs(aug[row][col]); mi = row
        aug[col], aug[mi] = aug[mi], aug[col]
        piv = aug[col][col]
        if abs(piv) < 1e-15:
            raise ValueError("Singular matrix")
        for j in range(2*n):
            aug[col][j] /= piv
        for row in range(n):
            if row == col:
                continue
            f = aug[row][col]
            for j in range(2*n):
                aug[row][j] -= f * aug[col][j]
    return [aug[i][n:] for i in range(n)]

def ols_fit(X, y):
    """OLS via normal equations, optimized for speed. X is list-of-rows, y is list. Returns (beta, RSS)."""
    n = len(y)
    k = len(X[0])
    # Build X'X and X'y directly (avoids full transpose + matmul)
    XtX = [[0.0]*k for _ in range(k)]
    Xty = [0.0]*k
    for i in range(n):
        xi = X[i]
        yi = y[i]
        for a in range(k):
            xia = xi[a]
            Xty[a] += xia * yi
            for b in range(a, k):
                v = xia * xi[b]
                XtX[a][b] += v
                if a != b:
                    XtX[b][a] += v
    # Solve via Gauss-Jordan
    XtX_inv = mat_inv(XtX)
    # beta = XtX_inv @ Xty
    beta = [0.0]*k
    for i in range(k):
        s = 0.0
        for j in range(k):
            s += XtX_inv[i][j] * Xty[j]
        beta[i] = s
    # RSS
    rss = 0.0
    for i in range(n):
        yhat = 0.0
        xi = X[i]
        for j in range(k):
            yhat += xi[j] * beta[j]
        rss += (y[i] - yhat)**2
    return beta, rss

# ─────────────────────────────────────────────────────────
# 2.  STATISTICAL HELPERS
# ─────────────────────────────────────────────────────────

def _log_gamma(x):
    """Lanczos approximation for log Gamma."""
    if x <= 0:
        return float('inf')
    coefs = [76.18009172947146, -86.50532032941677, 24.01409824083091,
             -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5]
    tmp = x + 5.5
    ser = 1.000000000190015
    for j, c in enumerate(coefs):
        ser += c / (x + 1 + j)
    return (x + 0.5) * math.log(tmp) - tmp + math.log(2.5066282746310005 * ser / x)

def _beta_cf(a, b, x, max_iter=200, eps=1e-12):
    """Continued-fraction evaluation for incomplete beta (Lentz)."""
    qab = a + b
    qap = a + 1.0
    qam = a - 1.0
    c = 1.0
    d = 1.0 - qab * x / qap
    if abs(d) < 1e-30: d = 1e-30
    d = 1.0 / d
    h = d
    for m in range(1, max_iter + 1):
        m2 = 2 * m
        # even step
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d
        if abs(d) < 1e-30: d = 1e-30
        c = 1.0 + aa / c
        if abs(c) < 1e-30: c = 1e-30
        d = 1.0 / d
        h *= d * c
        # odd step
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d
        if abs(d) < 1e-30: d = 1e-30
        c = 1.0 + aa / c
        if abs(c) < 1e-30: c = 1e-30
        d = 1.0 / d
        delta = d * c
        h *= delta
        if abs(delta - 1.0) < eps:
            return h
    return h

def _betainc(a, b, x):
    """Regularised incomplete beta I_x(a,b)."""
    if x < 0.0 or x > 1.0:
        return 0.0
    if x == 0.0 or x == 1.0:
        return x
    lbeta = _log_gamma(a + b) - _log_gamma(a) - _log_gamma(b)
    front = math.exp(lbeta + a * math.log(x) + b * math.log(1.0 - x))
    if x < (a + 1.0) / (a + b + 2.0):
        return front * _beta_cf(a, b, x) / a
    else:
        return 1.0 - front * _beta_cf(b, a, 1.0 - x) / b

def f_cdf(f_val, d1, d2):
    """CDF of F(d1, d2) distribution at f_val."""
    if f_val <= 0:
        return 0.0
    x = d1 * f_val / (d1 * f_val + d2)
    return _betainc(d1 / 2.0, d2 / 2.0, x)

def f_pvalue(f_val, d1, d2):
    """Upper-tail p-value: P(F > f_val)."""
    return 1.0 - f_cdf(f_val, d1, d2)

# ─────────────────────────────────────────────────────────
# 3.  DATA LOADING
# ─────────────────────────────────────────────────────────

BASE = os.path.dirname(os.path.abspath(__file__))

def load_sunspot(path):
    """Return dict {(year,month): ssn}."""
    data = {}
    with open(path, 'r') as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith('#'):
                continue
            parts = line.split(';')
            if len(parts) < 4:
                continue
            try:
                yr = int(parts[0].strip())
                mo = int(parts[1].strip())
                ssn = float(parts[3].strip())
            except (ValueError, IndexError):
                continue
            if ssn < 0:
                ssn = 0.0  # missing marker
            data[(yr, mo)] = ssn
    return data

def load_temperature(path):
    """Return dict {(year,month): anomaly}."""
    data = {}
    with open(path, 'r') as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith('#') or line.startswith('Date'):
                continue
            parts = line.split(',')
            if len(parts) < 2:
                continue
            try:
                date_str = parts[0].strip()
                anom = float(parts[1].strip())
                yr = int(date_str[:4])
                mo = int(date_str[4:])
            except (ValueError, IndexError):
                continue
            data[(yr, mo)] = anom
    return data

def load_gistemp(path):
    """Load NASA GISTEMP CSV: rows = years, columns = months."""
    data = {}
    with open(path, 'r') as f:
        lines = f.readlines()
    header_idx = None
    for i, line in enumerate(lines):
        if line.strip().startswith('Year'):
            header_idx = i
            break
    if header_idx is None:
        return data
    for line in lines[header_idx+1:]:
        parts = line.strip().split(',')
        if len(parts) < 13:
            continue
        try:
            yr = int(parts[0])
        except ValueError:
            continue
        for mo in range(1, 13):
            val = parts[mo].strip()
            if val == '***' or val == '':
                continue
            try:
                data[(yr, mo)] = float(val)
            except ValueError:
                continue
    return data

# Load data
sunspot_path = os.path.join(os.path.dirname(BASE), 'sunspot_cusum', 'monthly_sunspot.csv')
if not os.path.exists(sunspot_path):
    sunspot_path = os.path.join(BASE, '..', 'sunspot_cusum', 'monthly_sunspot.csv')

ssn_data = load_sunspot(sunspot_path)

# Try NOAA monthly first, fall back to GISTEMP
noaa_path = os.path.join(BASE, 'temperature_data.csv')
gistemp_path = os.path.join(BASE, 'gistemp.csv')

temp_data = load_temperature(noaa_path)
temp_source = "NOAA NCEI"
if len(temp_data) < 100:
    temp_data = load_gistemp(gistemp_path)
    temp_source = "NASA GISTEMP v4"

# Align on common months
common_keys = sorted(set(ssn_data.keys()) & set(temp_data.keys()))
# Filter to <= 2024-12
common_keys = [(y,m) for y,m in common_keys if y <= 2024]

ssn_series = [ssn_data[k] for k in common_keys]
temp_series = [temp_data[k] for k in common_keys]

print("=" * 70)
print("GRANGER CAUSALITY ANALYSIS: SUNSPOT NUMBER → GLOBAL TEMPERATURE")
print("=" * 70)
print(f"\nTemperature source: {temp_source}")
print(f"Sunspot source: SILSO v2.0 (Royal Observatory of Belgium)")
print(f"Common period: {common_keys[0][0]}-{common_keys[0][1]:02d} to "
      f"{common_keys[-1][0]}-{common_keys[-1][1]:02d}")
print(f"N observations: {len(common_keys)}")

# ─────────────────────────────────────────────────────────
# 4.  DESCRIPTIVE STATISTICS
# ─────────────────────────────────────────────────────────

def mean(x):
    return sum(x) / len(x)

def std(x):
    m = mean(x)
    return math.sqrt(sum((xi - m)**2 for xi in x) / (len(x) - 1))

def pearson_r(x, y):
    n = len(x)
    mx, my = mean(x), mean(y)
    sx, sy = std(x), std(y)
    if sx == 0 or sy == 0:
        return 0.0
    return sum((x[i]-mx)*(y[i]-my) for i in range(n)) / ((n-1)*sx*sy)

print(f"\n--- Descriptive Statistics ---")
print(f"Sunspot number: mean={mean(ssn_series):.2f}, std={std(ssn_series):.2f}, "
      f"min={min(ssn_series):.1f}, max={max(ssn_series):.1f}")
print(f"Temperature anomaly: mean={mean(temp_series):.4f}, std={std(temp_series):.4f}, "
      f"min={min(temp_series):.2f}, max={max(temp_series):.2f}")
print(f"Contemporaneous Pearson r = {pearson_r(ssn_series, temp_series):.4f}")

# ─────────────────────────────────────────────────────────
# 5.  CROSS-CORRELATION FUNCTION
# ─────────────────────────────────────────────────────────

print(f"\n{'='*70}")
print("CROSS-CORRELATION FUNCTION (lags -240 to +240 months)")
print("Positive lag k: sunspot leads temperature by k months")
print("="*70)

N = len(ssn_series)
mx_s = mean(ssn_series)
mx_t = mean(temp_series)
sd_s = std(ssn_series)
sd_t = std(temp_series)

ccf = {}
max_lag = 240
for lag in range(-max_lag, max_lag + 1):
    s = 0.0
    cnt = 0
    for i in range(N):
        j = i - lag  # sunspot at time j, temp at time i
        if 0 <= j < N:
            s += (ssn_series[j] - mx_s) * (temp_series[i] - mx_t)
            cnt += 1
    if cnt > 0:
        ccf[lag] = s / (cnt * sd_s * sd_t)

# Report key lags
print(f"\n{'Lag (months)':>14} {'CCF':>10}")
print("-" * 26)
key_lags = list(range(-240, 241, 12))
for lag in key_lags:
    if lag in ccf:
        print(f"{lag:>14d} {ccf[lag]:>10.4f}")

# Find peak positive and negative
pos_peak_lag = max(range(-max_lag, max_lag+1), key=lambda k: ccf.get(k, -999))
neg_peak_lag = min(range(-max_lag, max_lag+1), key=lambda k: ccf.get(k, 999))
print(f"\nPeak positive CCF: lag={pos_peak_lag}, r={ccf[pos_peak_lag]:.4f}")
print(f"Peak negative CCF: lag={neg_peak_lag}, r={ccf[neg_peak_lag]:.4f}")

# Also report every 6 months near peak
print(f"\nDetailed CCF around peak (every 6 months, lag -60 to +60):")
print(f"{'Lag':>8} {'CCF':>10}")
print("-" * 20)
for lag in range(-60, 61, 6):
    if lag in ccf:
        print(f"{lag:>8d} {ccf[lag]:>10.4f}")

# ─────────────────────────────────────────────────────────
# 6.  GRANGER CAUSALITY TEST
# ─────────────────────────────────────────────────────────

def granger_test(y_series, x_series, p, label=""):
    """
    Test if x Granger-causes y at lag order p.
    Restricted: y_t = c + sum(a_i * y_{t-i})
    Unrestricted: y_t = c + sum(a_i * y_{t-i}) + sum(b_j * x_{t-j})
    Returns (F, p_value, RSS_r, RSS_u, n_obs)
    """
    n = len(y_series)
    if p >= n // 2:
        return None
    
    # Build matrices
    y_dep = []
    X_r = []  # restricted
    X_u = []  # unrestricted
    
    for t in range(p, n):
        y_dep.append(y_series[t])
        # Restricted: intercept + p lags of y
        row_r = [1.0]
        for i in range(1, p+1):
            row_r.append(y_series[t-i])
        X_r.append(row_r)
        # Unrestricted: intercept + p lags of y + p lags of x
        row_u = row_r[:] 
        for j in range(1, p+1):
            row_u.append(x_series[t-j])
        X_u.append(row_u)
    
    n_obs = len(y_dep)
    k_r = p + 1
    k_u = 2*p + 1
    
    try:
        _, rss_r = ols_fit(X_r, y_dep)
        _, rss_u = ols_fit(X_u, y_dep)
    except (ValueError, ZeroDivisionError):
        return None
    
    dof_num = p  # number of restrictions
    dof_den = n_obs - k_u
    
    if dof_den <= 0 or rss_u <= 0:
        return None
    
    F = ((rss_r - rss_u) / dof_num) / (rss_u / dof_den)
    pval = f_pvalue(F, dof_num, dof_den)
    
    return (F, pval, rss_r, rss_u, n_obs)

lag_orders = [1, 6, 12, 24, 36, 60, 120, 132]

print(f"\n{'='*70}")
print("GRANGER CAUSALITY TEST: Sunspot → Temperature")
print("H0: Sunspot number does NOT Granger-cause temperature anomaly")
print("="*70)
print(f"\n{'p':>6} {'F-stat':>12} {'p-value':>12} {'RSS_r':>12} {'RSS_u':>12} {'n_obs':>8} {'Signif':>8}")
print("-" * 72)

granger_results_st = {}
for p in lag_orders:
    res = granger_test(temp_series, ssn_series, p, "SSN→Temp")
    if res:
        F, pval, rss_r, rss_u, n_obs = res
        sig = "***" if pval < 0.001 else "**" if pval < 0.01 else "*" if pval < 0.05 else ""
        print(f"{p:>6d} {F:>12.4f} {pval:>12.6f} {rss_r:>12.4f} {rss_u:>12.4f} {n_obs:>8d} {sig:>8}")
        granger_results_st[p] = (F, pval)

# ─────────────────────────────────────────────────────────
# 7.  BIDIRECTIONAL TEST: Temperature → Sunspot
# ─────────────────────────────────────────────────────────

print(f"\n{'='*70}")
print("BIDIRECTIONAL TEST: Temperature → Sunspot")
print("H0: Temperature anomaly does NOT Granger-cause sunspot number")
print("="*70)
print(f"\n{'p':>6} {'F-stat':>12} {'p-value':>12} {'RSS_r':>12} {'RSS_u':>12} {'n_obs':>8} {'Signif':>8}")
print("-" * 72)

granger_results_ts = {}
for p in lag_orders:
    res = granger_test(ssn_series, temp_series, p, "Temp→SSN")
    if res:
        F, pval, rss_r, rss_u, n_obs = res
        sig = "***" if pval < 0.001 else "**" if pval < 0.01 else "*" if pval < 0.05 else ""
        print(f"{p:>6d} {F:>12.4f} {pval:>12.6f} {rss_r:>12.4f} {rss_u:>12.4f} {n_obs:>8d} {sig:>8}")
        granger_results_ts[p] = (F, pval)

# ─────────────────────────────────────────────────────────
# 8.  PARTIAL CORRELATION (AR residuals)
# ─────────────────────────────────────────────────────────

print(f"\n{'='*70}")
print("PARTIAL CORRELATION: SSN ↔ Temp after removing AR(p) components")
print("="*70)

for ar_p in [12, 24, 60, 132]:
    n = len(temp_series)
    if ar_p >= n // 2:
        continue
    
    # Fit AR(p) to temperature
    y_t = []; X_t = []
    for t in range(ar_p, n):
        y_t.append(temp_series[t])
        row = [1.0] + [temp_series[t-i] for i in range(1, ar_p+1)]
        X_t.append(row)
    try:
        beta_t, _ = ols_fit(X_t, y_t)
        resid_t = []
        for idx in range(len(y_t)):
            yhat = sum(X_t[idx][j]*beta_t[j] for j in range(len(beta_t)))
            resid_t.append(y_t[idx] - yhat)
    except:
        continue
    
    # Fit AR(p) to sunspot
    y_s = []; X_s = []
    for t in range(ar_p, n):
        y_s.append(ssn_series[t])
        row = [1.0] + [ssn_series[t-i] for i in range(1, ar_p+1)]
        X_s.append(row)
    try:
        beta_s, _ = ols_fit(X_s, y_s)
        resid_s = []
        for idx in range(len(y_s)):
            yhat = sum(X_s[idx][j]*beta_s[j] for j in range(len(beta_s)))
            resid_s.append(y_s[idx] - yhat)
    except:
        continue
    
    r_partial = pearson_r(resid_t, resid_s)
    n_eff = len(resid_t)
    # t-test for correlation
    if abs(r_partial) < 1.0:
        t_stat = r_partial * math.sqrt((n_eff - 2) / (1 - r_partial**2))
    else:
        t_stat = float('inf')
    print(f"AR({ar_p:>3d}) residuals: r_partial = {r_partial:>8.4f}, t = {t_stat:>8.3f}, n = {n_eff}")

# ─────────────────────────────────────────────────────────
# 9.  PERMUTATION TEST (for key lag orders)
# ─────────────────────────────────────────────────────────

print(f"\n{'='*70}")
print("PERMUTATION TEST: Sunspot → Temperature (500 shuffles)")
print("="*70)

perm_lags = [12]
n_perm = 500

def granger_f_only(y_series, x_series, p):
    """Fast Granger F-stat only (no p-value computation)."""
    n = len(y_series)
    if p >= n // 2:
        return None
    y_dep = []
    X_r = []
    X_u = []
    for t in range(p, n):
        y_dep.append(y_series[t])
        row_r = [1.0]
        for i in range(1, p+1):
            row_r.append(y_series[t-i])
        X_r.append(row_r)
        row_u = row_r[:]
        for j in range(1, p+1):
            row_u.append(x_series[t-j])
        X_u.append(row_u)
    n_obs = len(y_dep)
    k_u = 2*p + 1
    try:
        _, rss_r = ols_fit(X_r, y_dep)
        _, rss_u = ols_fit(X_u, y_dep)
    except:
        return None
    dof_den = n_obs - k_u
    if dof_den <= 0 or rss_u <= 0:
        return None
    F = ((rss_r - rss_u) / p) / (rss_u / dof_den)
    return F

for p in perm_lags:
    F_actual = granger_f_only(temp_series, ssn_series, p)
    if F_actual is None:
        continue
    
    count_exceed = 0
    ssn_shuffled = ssn_series[:]
    for perm_i in range(n_perm):
        random.shuffle(ssn_shuffled)
        F_perm = granger_f_only(temp_series, ssn_shuffled, p)
        if F_perm is not None and F_perm >= F_actual:
            count_exceed += 1
    
    perm_pval = (count_exceed + 1) / (n_perm + 1)
    print(f"Lag p={p:>3d}: F_actual={F_actual:.4f}, "
          f"perm p-value={perm_pval:.4f} ({count_exceed}/{n_perm} exceeded)")

# ─────────────────────────────────────────────────────────
# 10.  SLIDING WINDOW GRANGER TEST
# ─────────────────────────────────────────────────────────

print(f"\n{'='*70}")
print("SLIDING WINDOW GRANGER TEST (50-year windows, step=5 years)")
print("Testing: Sunspot → Temperature, lag order p=12")
print("="*70)

window_months = 50 * 12  # 600 months
step_months = 5 * 12     # 60 months
test_p = 12

print(f"\n{'Window':>20} {'F-stat':>10} {'p-value':>10} {'Signif':>8}")
print("-" * 52)

sliding_results = []
start = 0
while start + window_months <= N:
    end = start + window_months
    y_win = temp_series[start:end]
    x_win = ssn_series[start:end]
    
    yr_start = common_keys[start][0]
    yr_end = common_keys[end-1][0]
    
    res = granger_test(y_win, x_win, test_p)
    if res:
        F, pval, _, _, n_obs = res
        sig = "***" if pval < 0.001 else "**" if pval < 0.01 else "*" if pval < 0.05 else ""
        print(f"{yr_start:>8d}-{yr_end:<8d} {F:>10.4f} {pval:>10.6f} {sig:>8}")
        sliding_results.append((yr_start, yr_end, F, pval))
    
    start += step_months

# Also do p=60 (5-year) sliding window
print(f"\n{'='*70}")
print("SLIDING WINDOW GRANGER TEST (50-year windows, step=5 years)")
print("Testing: Sunspot → Temperature, lag order p=60")
print("="*70)

test_p2 = 60
print(f"\n{'Window':>20} {'F-stat':>10} {'p-value':>10} {'Signif':>8}")
print("-" * 52)

start = 0
while start + window_months <= N:
    end = start + window_months
    y_win = temp_series[start:end]
    x_win = ssn_series[start:end]
    
    yr_start = common_keys[start][0]
    yr_end = common_keys[end-1][0]
    
    res = granger_test(y_win, x_win, test_p2)
    if res:
        F, pval, _, _, n_obs = res
        sig = "***" if pval < 0.001 else "**" if pval < 0.01 else "*" if pval < 0.05 else ""
        print(f"{yr_start:>8d}-{yr_end:<8d} {F:>10.4f} {pval:>10.6f} {sig:>8}")
    
    start += step_months

# ─────────────────────────────────────────────────────────
# 11.  DETRENDED ANALYSIS
# ─────────────────────────────────────────────────────────

print(f"\n{'='*70}")
print("DETRENDED ANALYSIS: Removing linear trend from temperature")
print("="*70)

# Fit linear trend to temperature
X_trend = [[1.0, i] for i in range(N)]
beta_trend, _ = ols_fit(X_trend, temp_series)
temp_detrended = [temp_series[i] - (beta_trend[0] + beta_trend[1]*i) for i in range(N)]

print(f"Linear trend: {beta_trend[1]*120:.4f} °C/decade")
print(f"Pearson r (SSN vs detrended temp) = {pearson_r(ssn_series, temp_detrended):.4f}")

print(f"\nGranger test on DETRENDED temperature:")
print(f"{'p':>6} {'F-stat':>12} {'p-value':>12} {'Signif':>8}")
print("-" * 42)

for p in lag_orders:
    res = granger_test(temp_detrended, ssn_series, p)
    if res:
        F, pval, _, _, n_obs = res
        sig = "***" if pval < 0.001 else "**" if pval < 0.01 else "*" if pval < 0.05 else ""
        print(f"{p:>6d} {F:>12.4f} {pval:>12.6f} {sig:>8}")

# ─────────────────────────────────────────────────────────
# 12.  SUMMARY
# ─────────────────────────────────────────────────────────

print(f"\n{'='*70}")
print("SUMMARY OF FINDINGS")
print("="*70)

print(f"""
1. DATA: {len(common_keys)} monthly observations, {common_keys[0][0]}-{common_keys[-1][0]}
   Sunspot source: SILSO v2.0
   Temperature source: {temp_source}

2. CROSS-CORRELATION:
   Peak positive CCF at lag {pos_peak_lag} months: r = {ccf[pos_peak_lag]:.4f}
   Peak negative CCF at lag {neg_peak_lag} months: r = {ccf[neg_peak_lag]:.4f}
   Contemporaneous (lag 0): r = {ccf[0]:.4f}

3. GRANGER CAUSALITY (Sunspot → Temperature):""")

for p in lag_orders:
    if p in granger_results_st:
        F, pval = granger_results_st[p]
        sig = "YES" if pval < 0.05 else "NO"
        print(f"   p={p:>3d}: F={F:.4f}, p={pval:.6f} → Significant at 5%: {sig}")

print(f"""
4. GRANGER CAUSALITY (Temperature → Sunspot):""")

for p in lag_orders:
    if p in granger_results_ts:
        F, pval = granger_results_ts[p]
        sig = "YES" if pval < 0.05 else "NO"
        print(f"   p={p:>3d}: F={F:.4f}, p={pval:.6f} → Significant at 5%: {sig}")

print(f"""
5. INTERPRETATION:
   - The Granger causality results are reported above with exact p-values.
   - The bidirectional test serves as a falsification check.
   - Sliding window analysis reveals temporal stability of the relationship.
   - Detrended analysis controls for the long-run warming trend.
""")

print("Analysis complete.")
```

## Expected Output Summary
- **N = 2,100** monthly observations (1850-01 to 2024-12)
- **Contemporaneous Pearson r = −0.003** (essentially zero)
- **Granger causality (SSN → Temp)**: No significance at any lag order (smallest p = 0.169 at p=24)
- **Granger causality (Temp → SSN)**: No significance except marginal at p=6 (p=0.048, does not survive Bonferroni correction)
- **Partial correlations**: All |r| < 0.013
- **Permutation test** (p=12): p-value = 0.513, confirming null
- **Sliding windows**: Only 1 of 26 windows significant at p=12 (consistent with chance)
- **Detrended**: All null

## Key Methods Implemented
1. **OLS via normal equations** with Gauss-Jordan matrix inversion (partial pivoting)
2. **F-distribution p-values** via regularized incomplete beta function (Lanczos log-gamma + Lentz continued fraction)
3. **Permutation inference** as non-parametric validation
4. **Cross-correlation function** at ±240 month lags
5. **Sliding-window temporal stability** analysis