Multivariate bubble tests

devtools::load_all("exuber", quiet = TRUE)

cat("=== 1. Formula check: coexplosive_stat() vs independent brute-force loop ===\n")
set.seed(1)
Tn <- 100
x <- rnorm(Tn)
y <- 2 + 0.5 * x + rnorm(Tn)
lag <- 2L

res <- exuber:::coexplosive_stat(y, x, lag)

# Independent brute-force recomputation: manual loop, manual OLS via lm(),
# manual cumulative sum loop (not vectorized the same way as the package code)
lo <- max(lag, 0L) + 1L
hi <- Tn + min(lag, 0L)
yy <- y[lo:hi]
xx <- x[(lo:hi) - lag]
fit_lm <- lm(yy ~ xx)
e_brute <- residuals(fit_lm)
n_brute <- length(e_brute)
sigma2_brute <- sum(e_brute^2) / n_brute
S_manual_sum <- 0
running <- 0
for (t in seq_along(e_brute)) {
  running <- running + e_brute[t]
  S_manual_sum <- S_manual_sum + running^2
}
S_brute <- S_manual_sum / (sigma2_brute * n_brute^2)

cat(sprintf("package S=%.8f  brute S=%.8f  diff=%.2e\n", res$S, S_brute, abs(res$S - S_brute)))
stopifnot(abs(res$S - S_brute) < 1e-8)
cat("OK: exact match\n\n")

cat("=== 2. Size under H0, homoskedastic errors (should be near nominal 5%) ===\n")
run_h0_homo <- function(seed) {
  set.seed(seed)
  Tn <- 150
  # x contains a genuine explosive episode (Model 1: unit root -> explosive)
  Te <- 90
  ex <- c(cumsum(rnorm(Te)), NA)
  ex <- ex[1:Te]
  expl <- ex[Te] * 1.05^(1:(Tn - Te)) + cumsum(rnorm(Tn - Te, sd = 0.3))
  x <- c(ex, expl)
  # y co-explosive with x at lag 0: y = a + b*x + iid noise
  y <- 1 + 0.8 * x + rnorm(Tn, sd = 1)
  out <- exuber:::cobubble_test(y, x, lag = 0L, nboot = 199L, seed = 1)
  out$reject
}
rejections <- sapply(1:100, run_h0_homo)
cat(sprintf("Empirical size (homoskedastic, 100 reps): %.3f (nominal 0.05)\n\n", mean(rejections)))

cat("=== 3. Size under H0, heteroskedastic errors (tests the wild bootstrap's whole point) ===\n")
run_h0_hetero <- function(seed) {
  set.seed(seed)
  Tn <- 150
  Te <- 90
  ex <- cumsum(rnorm(Te))
  expl <- ex[Te] * 1.05^(1:(Tn - Te)) + cumsum(rnorm(Tn - Te, sd = 0.3))
  x <- c(ex, expl)
  # heteroskedastic errors: sd jumps from 1 to 4 partway through
  sd_pattern <- c(rep(1, Tn %/% 2), rep(4, Tn - Tn %/% 2))
  y <- 1 + 0.8 * x + rnorm(Tn, sd = sd_pattern)
  out <- exuber:::cobubble_test(y, x, lag = 0L, nboot = 199L, seed = 1)
  out$reject
}
rejections_hetero <- sapply(1:100, run_h0_hetero)
cat(sprintf("Empirical size (heteroskedastic, 100 reps): %.3f (nominal 0.05)\n\n", mean(rejections_hetero)))

cat("=== 4. Power under H1: y and x are NOT co-explosive (independent explosive episodes) ===\n")
run_h1 <- function(seed) {
  set.seed(seed)
  Tn <- 150
  Te <- 90
  ex <- cumsum(rnorm(Te))
  expl_x <- ex[Te] * 1.05^(1:(Tn - Te)) + cumsum(rnorm(Tn - Te, sd = 0.3))
  x <- c(ex, expl_x)
  # y has its OWN independent explosive episode (not driven by x at all)
  ey <- cumsum(rnorm(Te))
  expl_y <- ey[Te] * 1.05^(1:(Tn - Te)) + cumsum(rnorm(Tn - Te, sd = 0.3))
  y <- ey_full <- c(ey, expl_y)
  out <- exuber:::cobubble_test(y, x, lag = 0L, nboot = 199L, seed = 1)
  out$reject
}
rejections_h1 <- sapply(1:60, run_h1)
cat(sprintf("Empirical power (H1, independent explosive episodes, 60 reps): %.3f (should be high)\n\n", mean(rejections_h1)))

cat("=== 5. Lag recovery: true lag = 3, does coexplosive_select_lag() find it? ===\n")
run_lag_recovery <- function(seed) {
  set.seed(seed)
  Tn <- 200
  Te <- 120
  true_lag <- 3L
  ex <- cumsum(rnorm(Te))
  expl <- ex[Te] * 1.06^(1:(Tn - Te)) + cumsum(rnorm(Tn - Te, sd = 0.3))
  x <- c(ex, expl)
  # y co-explosive with x lagged by true_lag: y_t = a + b*x_{t-true_lag} + noise
  # build y such that y[t] depends on x[t - true_lag]
  y <- rep(NA_real_, Tn)
  for (t in (true_lag + 1):Tn) y[t] <- 1 + 0.8 * x[t - true_lag] + rnorm(1, sd = 0.5)
  y[1:true_lag] <- x[1:true_lag] + rnorm(true_lag, sd = 0.5)
  est_lag <- exuber:::coexplosive_select_lag(y, x, lags = -6:6)
  est_lag
}
lags_est <- sapply(1:20, run_lag_recovery)
cat("Estimated lags across 20 seeds (true = 3):\n")
print(table(lags_est))

"""Validation of cobubble_test() -- Evripidou, Harvey, Leybourne & Sollis
(2022)'s co-explosive behaviour test. Same folder/base name as the R
script it cross-checks: radf_cobubble_validation.R. See
docs/multivariate.md, "Co-bubble test", for the full write-up.

No RNG-bit-exact comparison against R is attempted (numpy's Generator vs
R's RNG -- see sim.py's module docstring for why that's a project-wide
non-goal); this reproduces the SAME five checks the R script runs, each
against its own independently-generated data, and (like the R script)
mostly cares about direction/magnitude (empirical size near nominal,
power high, lag recovered) rather than exact numbers. Run standalone from
the exuber-project/ root:

  uv run --project pyexuber python docs/replication/multivariate/radf_cobubble_validation.py

Not imported by pyexuber's own pytest suite for the same repo-boundary
reason as radf_common_validation.py (see that file's docstring); these
checks are re-implemented directly in pyexuber/tests/test_multivariate.py.
"""

import numpy as np

from exuber.cobubble_test import _coexplosive_select_lag, _coexplosive_stat, cobubble_test


def check_formula_exact(seed: int = 1) -> None:
    """coexplosive_stat() vs. an independently written brute-force
    computation (separate lstsq call, manual loop-based cumulative sum
    instead of vectorized cumsum) -- mirrors the R script's check 1."""
    rng = np.random.default_rng(seed)
    tn = 100
    x = rng.normal(size=tn)
    y = 2 + 0.5 * x + rng.normal(size=tn)
    lag = 2

    S, _resid, _sigma2, _n = _coexplosive_stat(y, x, lag)

    lo, hi = max(lag, 0), tn + min(lag, 0)
    yy, xx = y[lo:hi], x[lo - lag : hi - lag]
    design = np.column_stack([np.ones(len(yy)), xx])
    beta, *_ = np.linalg.lstsq(design, yy, rcond=None)
    e_brute = yy - design @ beta
    n_brute = len(e_brute)
    sigma2_brute = np.sum(e_brute**2) / n_brute
    running = 0.0
    s_manual_sum = 0.0
    for e in e_brute:
        running += e
        s_manual_sum += running**2
    s_brute = s_manual_sum / (sigma2_brute * n_brute**2)

    diff = abs(S - s_brute)
    print(f"package S={S:.8f}  brute S={s_brute:.8f}  diff={diff:.2e}")
    assert diff < 1e-8
    print("OK: exact match\n")


def _build_dgp(seed: int, hetero: bool) -> tuple[np.ndarray, np.ndarray]:
    rng = np.random.default_rng(seed)
    tn, te = 150, 90
    ex = np.cumsum(rng.normal(size=te))
    expl = ex[-1] * 1.05 ** np.arange(1, tn - te + 1) + np.cumsum(
        rng.normal(scale=0.3, size=tn - te)
    )
    x = np.concatenate([ex, expl])
    if hetero:
        sd_pattern = np.concatenate([np.ones(tn // 2), np.full(tn - tn // 2, 4.0)])
        y = 1 + 0.8 * x + rng.normal(size=tn) * sd_pattern
    else:
        y = 1 + 0.8 * x + rng.normal(size=tn)
    return y, x


def check_size(seed_base: int, hetero: bool, reps: int = 100) -> float:
    rejections = []
    for i in range(reps):
        y, x = _build_dgp(seed_base + i, hetero)
        out = cobubble_test(y, x, lag=0, nboot=199, seed=1)
        rejections.append(out.reject)
    size = float(np.mean(rejections))
    label = "heteroskedastic" if hetero else "homoskedastic"
    print(f"Empirical size ({label}, {reps} reps): {size:.3f} (nominal 0.05)")
    return size


def check_power(seed_base: int = 3000, reps: int = 60) -> float:
    rejections = []
    for i in range(reps):
        rng = np.random.default_rng(seed_base + i)
        tn, te = 150, 90
        ex = np.cumsum(rng.normal(size=te))
        expl_x = ex[-1] * 1.05 ** np.arange(1, tn - te + 1) + np.cumsum(
            rng.normal(scale=0.3, size=tn - te)
        )
        x = np.concatenate([ex, expl_x])
        ey = np.cumsum(rng.normal(size=te))
        expl_y = ey[-1] * 1.05 ** np.arange(1, tn - te + 1) + np.cumsum(
            rng.normal(scale=0.3, size=tn - te)
        )
        y = np.concatenate([ey, expl_y])
        out = cobubble_test(y, x, lag=0, nboot=199, seed=1)
        rejections.append(out.reject)
    power = float(np.mean(rejections))
    print(f"Empirical power (H1, independent explosive episodes, {reps} reps): {power:.3f}")
    return power


def check_lag_recovery(true_lag: int = 3, seeds: int = 20) -> list[int]:
    estimated = []
    for seed in range(1, seeds + 1):
        rng = np.random.default_rng(seed)
        tn, te = 200, 120
        ex = np.cumsum(rng.normal(size=te))
        expl = ex[-1] * 1.06 ** np.arange(1, tn - te + 1) + np.cumsum(
            rng.normal(scale=0.3, size=tn - te)
        )
        x = np.concatenate([ex, expl])
        y = np.full(tn, np.nan)
        for t in range(true_lag, tn):
            y[t] = 1 + 0.8 * x[t - true_lag] + rng.normal(scale=0.5)
        y[:true_lag] = x[:true_lag] + rng.normal(scale=0.5, size=true_lag)
        estimated.append(_coexplosive_select_lag(y, x, range(-6, 7)))
    print(f"Estimated lags across {seeds} seeds (true={true_lag}): {estimated}")
    return estimated


if __name__ == "__main__":
    print("=== 1. Formula check ===")
    check_formula_exact()
    print("=== 2. Size under H0, homoskedastic errors ===")
    check_size(seed_base=1, hetero=False)
    print("\n=== 3. Size under H0, heteroskedastic errors ===")
    check_size(seed_base=1, hetero=True)
    print("\n=== 4. Power under H1 ===")
    check_power()
    print("\n=== 5. Lag recovery (true lag = 3) ===")
    check_lag_recovery()

# Replication script for radf_common() and radf_common_cv() (common-bubble
# detection via PCA + PSY, Chen, Phillips & Shi 2023). See
# docs/multivariate.md, "Common-bubble detection via PCA + PSY".
#
# Theorem 4.3 states that the null limit of the PSY statistic on the first
# principal component equals the univariate GSADF null, whatever the panel
# width N. At practical N the null quantile grows with N instead. This is why
# radf_common_cv() simulates an N-dependent null and radf_mc_cv() is not used.
Sys.setenv(NOT_CRAN = "true")
options(exuber.parallel = FALSE, exuber.show_progress = FALSE)
devtools::load_all("exuber", quiet = TRUE)

cat("=== radf_common()'s null quantiles grow with panel width N ===\n")
cat("(sharpest possible null: N independent random walks, no true common factor)\n")
cat("T = 250 throughout, seed = 271828\n\n")

cv_univariate <- radf_mc_cv(250, nrep = 2000, seed = 1)
cat("radf_mc_cv(250)'s 95% gsadf_cv (the univariate benchmark):",
  round(cv_univariate$gsadf_cv["95%"], 3), "\n\n"
)

set.seed(271828)
results <- list()
for (N in c(6, 20, 50, 100)) {
  gsadf_null <- replicate(400, {
    panel <- replicate(N, cumsum(rnorm(250)))
    radf_common(panel)$gsadf
  })
  q <- quantile(gsadf_null, c(0.9, 0.95, 0.99))
  results[[as.character(N)]] <- q
  cat(sprintf(
    "N=%3d: 90%%=%.3f 95%%=%.3f 99%%=%.3f  (diff from radf_mc_cv 95%%: %+.2f)\n",
    N, q[1], q[2], q[3], q[2] - cv_univariate$gsadf_cv["95%"]
  ))
}

cat("\nPublished (doc): N=6 -> (2.420,2.662,3.057); N=20 -> (3.225,3.451,4.144);\n")
cat("                 N=50 -> (4.057,4.365,4.921); N=100 -> (4.880,5.203,5.836)\n\n")

cat("=== radf_common_cv() reproduces this N-dependence directly ===\n")
cv_small <- radf_common_cv(n = 250, N = 4, minw = NULL, nrep = 500, seed = 1)
cv_large <- radf_common_cv(n = 250, N = 30, minw = NULL, nrep = 500, seed = 1)
cat("radf_common_cv(N=4)  gsadf_cv:", round(cv_small$gsadf_cv, 3), "\n")
cat("radf_common_cv(N=30) gsadf_cv:", round(cv_large$gsadf_cv, 3), "\n")
cat("N=30's null quantile is higher than N=4's:",
  all(cv_large$gsadf_cv > cv_small$gsadf_cv), "\n\n"
)

cat("=== Full test-common.R suite ===\n")
testthat::test_file(
  "exuber/tests/testthat/test-common.R",
  reporter = "summary"
)

"""Validation of radf_common()/radf_common_cv() -- Chen, Phillips & Shi
(2023)'s common-bubble detection via PCA + PSY. Same folder/base name as
the R script it cross-checks: radf_common_validation.R. See
docs/multivariate.md, "Common-bubble detection via PCA + PSY", for the
full narrative: Theorem
4.3's claim that the PSY-on-PC1 statistic's null is asymptotically
identical to plain univariate GSADF's does NOT hold at practical panel
widths N -- the true null quantile *grows* with N. That is why
radf_common_cv() (not radf_mc_cv()) must be used for critical values, and
why this script reproduces the SAME directional finding rather than
hiding it, per this project's instruction to keep the honest caveat in
the Python port too.

Reference numbers in the comments below (not asserted bit-for-bit --
cv.py's own module docstring already establishes numpy's Generator vs R's
RNG can't be matched bit-for-bit for any stochastic simulation in this
project) were produced by:

  Rscript -e '
    Sys.setenv(NOT_CRAN = "true")
    options(exuber.parallel = FALSE, exuber.show_progress = FALSE)
    devtools::load_all("exuber", quiet = TRUE)
    cv4  <- radf_common_cv(n = 100, N = 4,  nrep = 300, seed = 42)
    cv20 <- radf_common_cv(n = 100, N = 20, nrep = 300, seed = 42)
    cat("N=4  gsadf_cv:", cv4$gsadf_cv, "\n")
    cat("N=20 gsadf_cv:", cv20$gsadf_cv, "\n")
  '
  (run from exuber-project/ root)

  N=4  gsadf_cv: 2.184956 2.393796 3.348111   (90/95/99%)
  N=20 gsadf_cv: 3.281803 3.478511 3.968891   (90/95/99%)

  -- confirms the N-dependence: N=20's 95% (3.48) is well above N=4's
  (2.39), the same direction as docs/multivariate.md's own published
  table (N=6: 2.66, N=20: 3.45-3.48, N=50: 4.33-4.37, N=100: 5.20).

Two checks, run standalone from the exuber-project/ root:
  uv run --project pyexuber python docs/replication/multivariate/radf_common_validation.py
  1. PCA formula-exact check (pure numpy, no _core/radf() needed): the
     module's internal _pca() helper against an independent
     eigendecomposition of the sample covariance matrix.
  2. N-dependence sanity check (needs the compiled _core extension --
     see pyexuber/CLAUDE.md for why it can't be built on every machine):
     radf_common_cv() at two panel widths, checking the 95% gsadf_cv
     grows with N, matching the R finding above.

Not imported by pyexuber's own pytest suite: pyexuber is developed and
CI-tested in its own repo (kvasilopoulos/pyexuber), which doesn't check
out this umbrella repo's docs/ -- so pyexuber/tests/test_multivariate.py
re-implements these same two checks directly instead of importing this
file across that repo boundary. Keep the two in sync by hand if either
changes.
"""

import numpy as np

from exuber.radf_common import _pca, radf_common_cv


def check_pca_formula_exact(seed: int = 0) -> None:
    rng = np.random.default_rng(seed)
    n, nc, r = 60, 5, 2
    x = np.cumsum(rng.normal(size=(n, nc)), axis=0)

    loadings, scores, explained = _pca(x, r)

    # Independent recomputation via eigendecomposition of the sample
    # covariance matrix -- a different numerical path than SVD.
    xc = x - x.mean(axis=0)
    cov = (xc.T @ xc) / (n - 1)
    eigvals, eigvecs = np.linalg.eigh(cov)
    order = np.argsort(eigvals)[::-1]
    eigvals, eigvecs = eigvals[order], eigvecs[:, order]

    # Eigenvector sign is arbitrary -- align each column's sign before
    # comparing, then require agreement to near machine precision.
    for j in range(r):
        sign = np.sign(np.dot(loadings[:, j], eigvecs[:, j])) or 1.0
        diff = float(np.max(np.abs(loadings[:, j] - sign * eigvecs[:, j])))
        assert diff < 1e-10, f"loadings column {j}: diff={diff:.2e}"

    expected_ratio = eigvals[:r] / np.sum(eigvals)
    assert np.allclose(explained, expected_ratio, atol=1e-10)
    assert np.allclose(scores, xc @ loadings, atol=1e-10)

    print(f"PCA formula-exact check OK (max loadings diff < 1e-10, r={r})")


def check_cv_grows_with_n(seed: int = 1) -> None:
    cv_small = radf_common_cv(n=60, N=4, nrep=150, seed=seed)
    cv_large = radf_common_cv(n=60, N=20, nrep=150, seed=seed)

    small_95 = cv_small.gsadf_cv[1]
    large_95 = cv_large.gsadf_cv[1]
    print(f"N=4  95% gsadf_cv: {small_95:.3f}")
    print(f"N=20 95% gsadf_cv: {large_95:.3f}")
    assert large_95 > small_95, (
        "expected N=20's null quantile to exceed N=4's (see R reference numbers "
        "in this file's module docstring) -- Theorem 4.3's N-independence claim "
        "does not hold at practical panel widths"
    )
    print("N-dependence check OK: N=20's 95% gsadf_cv > N=4's")


if __name__ == "__main__":
    check_pca_formula_exact()
    check_cv_grows_with_n()

# Validation of contagion_reg(), the bubble contagion regression of
# Greenaway-McGrevy & Phillips (2016): the fixed-window AR(1) sequence, the
# single-delay Nadaraya-Watson regression and the leave-one-out bandwidth.
# See docs/multivariate.md, "Contagion regression".
#
# No published numeric table exists to validate against, because the source
# paper reports its results as Figures 7-8. Each closed-form piece is checked
# against a brute-force computation instead, followed by a directional check
# of the estimator.
#
# Run from the exuber-project/ root.

devtools::load_all("exuber", quiet = TRUE)

cat("=== 1. Fixed-window AR(1) coefficient sequence (eq. 1) vs. brute-force lm() ===\n")
set.seed(1)
n <- 150
S <- 50
core <- cumsum(rnorm(n))
beta_core <- exuber:::contagion_fixed_window_beta(core, S)
for (t_check in c(60, 80, 100, 130, 150)) {
  win <- (t_check - S + 1):t_check
  fit <- lm(core[win[-1]] ~ core[win[-length(win)]])
  cf <- beta_core[as.character(t_check)]
  cat(sprintf(
    "t=%d closed-form=%.8f lm()=%.8f |diff|=%.2e\n",
    t_check, cf, coef(fit)[2], abs(cf - coef(fit)[2])
  ))
}

cat("\n=== 2. Nadaraya-Watson ratio (eq. 6) vs. manual weighted-least-squares ===\n")
y <- 0.5 * core + cumsum(rnorm(n, sd = 0.5))
bc <- exuber:::contagion_fixed_window_beta(core, S)
bj <- exuber:::contagion_fixed_window_beta(y, S)
d <- 2
r_test <- 0.5
h_test <- 0.2
fast <- exuber:::contagion_nw_delta2(bc, bj, n, r_test, h_test, d)
s <- as.integer(names(bj))
core_shift <- (bc - mean(bc))[as.character(s - d)]
valid <- !is.na(core_shift)
s2 <- s[valid]
bjc <- (bj - mean(bj))[valid]
csh <- core_shift[valid]
w <- dnorm((s2 / n - r_test) / h_test) / h_test
manual <- sum(w * bjc * csh) / sum(w * csh^2)
cat(sprintf("fast=%.10f manual=%.10f |diff|=%.2e\n", fast, manual, abs(fast - manual)))

cat("\n=== 3. LOOCV SSE (eq. 7) vs. manual leave-one-out double loop ===\n")
fast_sse <- exuber:::contagion_loocv_sse(h_test, bc, bj, n, d)
m <- length(s2)
manual_sse <- 0
for (i in seq_len(m)) {
  r_i <- s2[i] / m
  num <- 0
  den <- 0
  for (p in seq_len(m)) {
    if (p == i) next
    wp <- dnorm((s2[p] / n - r_i) / h_test) / h_test
    num <- num + wp * bjc[p] * csh[p]
    den <- den + wp * csh[p]^2
  }
  pred <- (num / den) * csh[i]
  manual_sse <- manual_sse + (bjc[i] - pred)^2
}
cat(sprintf("fast=%.10f manual=%.10f |diff|=%.2e\n", fast_sse, manual_sse, abs(fast_sse - manual_sse)))

cat("\n=== 4. Bandwidth CV: interior optimum, SSE no worse than either H_T endpoint ===\n")
h_opt <- exuber:::contagion_bandwidth_cv(bc, bj, n, d)
H_T <- c(m^(-1 / 2), m^(-1 / 10))
cat("H_T:", H_T, " h_opt:", h_opt, "\n")
cat(
  "SSE at h_opt:", exuber:::contagion_loocv_sse(h_opt, bc, bj, n, d),
  " SSE at H_T[1]:", exuber:::contagion_loocv_sse(H_T[1], bc, bj, n, d),
  " SSE at H_T[2]:", exuber:::contagion_loocv_sse(H_T[2], bc, bj, n, d), "\n"
)

cat("\n=== 5. Directional sensible-behavior check: planted contagion vs. independent series ===\n")
set.seed(42)
nrep <- 15
planted_range <- indep_range <- numeric(nrep)
for (i in seq_len(nrep)) {
  set.seed(1000 + i)
  core_i <- cumsum(rnorm(n))
  y_planted <- numeric(n)
  y_planted[1:10] <- rnorm(10)
  for (t in 11:n) {
    local_rho <- 0.5 + 0.4 * tanh((core_i[max(t - 3, 1)] - core_i[max(t - 13, 1)]) / 5)
    y_planted[t] <- local_rho * y_planted[t - 1] + rnorm(1)
  }
  y_indep <- cumsum(rnorm(n))

  out_planted <- contagion_reg(y_planted, core_i, S = S, d = 3, h = 0.3)
  out_indep <- contagion_reg(y_indep, core_i, S = S, d = 3, h = 0.3)
  planted_range[i] <- diff(range(out_planted$delta2))
  indep_range[i] <- diff(range(out_indep$delta2))
}
cat("mean range(delta2), planted contagion: ", mean(planted_range), "\n")
cat("mean range(delta2), independent series:", mean(indep_range), "\n")

cat("\ndone\n")

"""Validation of contagion_reg() -- Greenaway-McGrevy & Phillips (2016)'s
bubble contagion regression, minimum-viable subset (fixed-window AR(1)
sequence, single-delay Nadaraya-Watson regression, LOOCV bandwidth). Same
folder/base name as the R script it cross-checks: radf_contagion_validation.R.
See docs/multivariate.md, "Contagion regression", for the full write-up
(the window-width convention and the matrix orientation of the LOOCV SSE
helper are both checked directly here).

No published numeric table exists to validate against -- the source
paper's own results are Figures 7-8, not tabulated numbers (same
situation the R script's own header documents). Validated instead via
brute-force cross-checks of each closed-form piece (to near machine
precision, no RNG involved -- these are NOT stochastic checks, unlike
radf_common_validation.py/radf_cobubble_validation.py), plus a
directional sensible-behavior check. Run standalone from the
exuber-project/ root:

  uv run --project pyexuber python docs/replication/multivariate/radf_contagion_validation.py

Not imported by pyexuber's own pytest suite for the same repo-boundary
reason as radf_common_validation.py (see that file's docstring); these
checks are re-implemented directly in pyexuber/tests/test_multivariate.py.
"""

import numpy as np

from exuber.contagion_reg import (
    _contagion_bandwidth_cv,
    _contagion_fixed_window_beta,
    _contagion_loocv_sse,
    _contagion_nw_delta2,
    contagion_reg,
)


def check_fixed_window_beta(seed: int = 1) -> None:
    """eq. 1 vs. brute-force lstsq at several window-end dates."""
    rng = np.random.default_rng(seed)
    n, S = 150, 50
    core = np.cumsum(rng.normal(size=n))
    t_core, beta_core = _contagion_fixed_window_beta(core, S)

    pos = {int(t): i for i, t in enumerate(t_core)}
    for t_check in (60, 80, 100, 130, 150):
        win = core[t_check - S : t_check]  # S levels, python 0-indexed
        design = np.column_stack([np.ones(S - 1), win[:-1]])
        beta_lstsq, *_ = np.linalg.lstsq(design, win[1:], rcond=None)
        cf = beta_core[pos[t_check]]
        diff = abs(cf - beta_lstsq[1])
        print(f"t={t_check} closed-form={cf:.8f} lstsq={beta_lstsq[1]:.8f} |diff|={diff:.2e}")
        assert diff < 1e-8


def check_nw_ratio(seed: int = 1) -> None:
    """eq. 6 vs. a manual weighted-least-squares ratio."""
    rng = np.random.default_rng(seed)
    n, S = 150, 50
    core = np.cumsum(rng.normal(size=n))
    y = 0.5 * core + np.cumsum(rng.normal(scale=0.5, size=n))
    t_core, beta_core = _contagion_fixed_window_beta(core, S)
    t_j, beta_j = _contagion_fixed_window_beta(y, S)
    d, r_test, h_test = 2, np.array([0.5]), 0.2

    fast = _contagion_nw_delta2(t_core, beta_core, t_j, beta_j, n, r_test, h_test, d)[0]

    bcore_c = beta_core - beta_core.mean()
    bj_c = beta_j - beta_j.mean()
    pos = {int(t): i for i, t in enumerate(t_core)}
    idx = np.array([pos.get(int(t) - d, -1) for t in t_j])
    valid = idx >= 0
    s2 = t_j[valid]
    bjc = bj_c[valid]
    csh = bcore_c[idx[valid]]
    w = np.exp(-0.5 * ((s2 / n - r_test[0]) / h_test) ** 2) / np.sqrt(2 * np.pi) / h_test
    manual = np.sum(w * bjc * csh) / np.sum(w * csh**2)

    diff = abs(fast - manual)
    print(f"fast={fast:.10f} manual={manual:.10f} |diff|={diff:.2e}")
    assert diff < 1e-10
    return t_core, beta_core, t_j, beta_j, bjc, csh, s2, n, d, h_test


def check_loocv_sse(t_core, beta_core, t_j, beta_j, bjc, csh, s2, n, d, h_test) -> None:
    """eq. 7 vs. a manual leave-one-out double loop (the correct
    orientation is K.T @ v, not K @ v, since the kernel weight matrix isn't
    symmetric)."""
    fast_sse = _contagion_loocv_sse(h_test, t_core, beta_core, t_j, beta_j, n, d)

    m = len(s2)
    manual_sse = 0.0
    for i in range(m):
        r_i = s2[i] / m
        num = den = 0.0
        for p in range(m):
            if p == i:
                continue
            wp = np.exp(-0.5 * ((s2[p] / n - r_i) / h_test) ** 2) / np.sqrt(2 * np.pi) / h_test
            num += wp * bjc[p] * csh[p]
            den += wp * csh[p] ** 2
        pred = (num / den) * csh[i]
        manual_sse += (bjc[i] - pred) ** 2

    diff = abs(fast_sse - manual_sse)
    print(f"fast={fast_sse:.10f} manual={manual_sse:.10f} |diff|={diff:.2e}")
    assert diff < 1e-8


def check_bandwidth_cv(t_core, beta_core, t_j, beta_j, n, d) -> None:
    """Interior optimum, SSE no worse than either H_T endpoint (eq. 7)."""
    h_opt = _contagion_bandwidth_cv(t_core, beta_core, t_j, beta_j, n, d)
    m = len(beta_j)
    H_T = (m ** (-1 / 2), m ** (-1 / 10))
    sse_opt = _contagion_loocv_sse(h_opt, t_core, beta_core, t_j, beta_j, n, d)
    sse_lo = _contagion_loocv_sse(H_T[0], t_core, beta_core, t_j, beta_j, n, d)
    sse_hi = _contagion_loocv_sse(H_T[1], t_core, beta_core, t_j, beta_j, n, d)
    print(f"H_T: {H_T}  h_opt: {h_opt}")
    print(f"SSE at h_opt: {sse_opt}  SSE at H_T[0]: {sse_lo}  SSE at H_T[1]: {sse_hi}")
    assert sse_opt <= sse_lo + 1e-8
    assert sse_opt <= sse_hi + 1e-8


def check_sensible_behavior(seed_base: int = 1000, nrep: int = 15) -> None:
    """Directional check (no ground truth to match): a satellite series
    whose local persistence genuinely tracks the core's own should show a
    visibly wider range of estimated delta_2(r) than an independent series."""
    n, S = 150, 50
    planted_range = np.empty(nrep)
    indep_range = np.empty(nrep)
    for i in range(nrep):
        rng = np.random.default_rng(seed_base + i)
        core_i = np.cumsum(rng.normal(size=n))
        y_planted = np.empty(n)
        y_planted[:10] = rng.normal(size=10)
        for t in range(10, n):
            local_rho = 0.5 + 0.4 * np.tanh(
                (core_i[max(t - 3, 0)] - core_i[max(t - 13, 0)]) / 5
            )
            y_planted[t] = local_rho * y_planted[t - 1] + rng.normal()
        y_indep = np.cumsum(rng.normal(size=n))

        out_planted = contagion_reg(y_planted, core_i, S=S, d=3, h=0.3)
        out_indep = contagion_reg(y_indep, core_i, S=S, d=3, h=0.3)
        planted_range[i] = np.ptp(out_planted.delta2)
        indep_range[i] = np.ptp(out_indep.delta2)

    print(f"mean range(delta2), planted contagion: {planted_range.mean():.4f}")
    print(f"mean range(delta2), independent series: {indep_range.mean():.4f}")


if __name__ == "__main__":
    print("=== 1. Fixed-window AR(1) coefficient sequence (eq. 1) ===")
    check_fixed_window_beta()
    print("\n=== 2. Nadaraya-Watson ratio (eq. 6) ===")
    ctx = check_nw_ratio()
    print("\n=== 3. LOOCV SSE (eq. 7) ===")
    check_loocv_sse(*ctx)
    print("\n=== 4. Bandwidth CV ===")
    check_bandwidth_cv(ctx[0], ctx[1], ctx[2], ctx[3], ctx[7], ctx[8])
    print("\n=== 5. Directional sensible-behavior check ===")
    check_sensible_behavior()
    print("\ndone")