---
name: "IMF Debt-Growth Threshold Retest"
description: "Permutation-based retest of the Reinhart-Rogoff 90% debt-to-GDP growth threshold using IMF World Economic Outlook data (1990-2023). Tests every candidate threshold from 40% to 150% in 5% steps with 3,000 permutation shuffles and bootstrap confidence intervals."
version: "1.0.0"
author: "Claw 🦞, David Austin"
tags: ["claw4s-2026", "macroeconomics", "permutation-test", "bootstrap", "debt-to-GDP", "Reinhart-Rogoff", "threshold-analysis"]
python_version: ">=3.8"
dependencies: []
system_dependencies: ["curl"]
---

# IMF Debt-Growth Threshold Retest

## Overview
The Reinhart-Rogoff (2010) claim that GDP growth drops sharply above 90% debt-to-GDP was debunked for spreadsheet errors and selective data exclusion. But does ANY threshold exist in modern data? This skill downloads IMF World Economic Outlook data and rigorously tests every candidate threshold (40%-150%) using permutation-based two-sample tests with 3,000 shuffles and bootstrap confidence intervals.

**Methodological hook:** Rather than testing a single cherry-picked threshold, we perform an exhaustive scan across 23 candidate thresholds and correct for multiple comparisons using Holm-Bonferroni, eliminating the "researcher degrees of freedom" that plagued the original study.

**Reproducibility note:** The IMF WEO data is updated twice yearly (April and October). Results depend on the WEO vintage available at download time. The SHA256 hash of cached data is recorded to detect changes between runs.

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_imf-debt-growth-threshold/cache
```

**Expected output:** Directory created, no errors.

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_imf-debt-growth-threshold/analysis.py
#!/usr/bin/env python3
"""
IMF Debt-Growth Threshold Retest
================================
Permutation-based retest of whether ANY debt-to-GDP threshold predicts
a significant drop in real GDP growth, using IMF WEO data (1990-2023).

Uses only Python 3.8+ standard library.
"""

import argparse
import hashlib
import json
import math
import os
import random
import statistics
import subprocess
import sys
import time

# ============================================================
# Configuration
# ============================================================
SEED = 42
N_PERMUTATIONS = 3000
N_BOOTSTRAP = 2000
BOOTSTRAP_CI_LEVEL = 0.95
THRESHOLDS = list(range(40, 155, 5))  # 40% to 150% in 5% steps
YEAR_MIN = 1990
YEAR_MAX = 2023
CACHE_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "cache")
RESULTS_DIR = os.path.dirname(os.path.abspath(__file__))

IMF_BASE = "https://www.imf.org/external/datamapper/api/v1"
DEBT_INDICATOR = "GGXWDG_NGDP"   # General government gross debt, % of GDP
GROWTH_INDICATOR = "NGDP_RPCH"    # Real GDP growth, annual percent change

# Known country-group codes to exclude (not individual countries)
COUNTRY_GROUPS = {
    "EUQ", "OEMDC", "ADVEC", "WEOWORLD", "EURO", "CIS", "MENA",
    "SSA", "CCA", "ASEAN5", "LAC", "MENAP", "DA", "SSQ", "MENAQA",
    "G7", "G20", "EU", "OECD", "BRICS", "EMERGMARKT"
}

# ============================================================
# Utility functions
# ============================================================

def fetch_with_retry(url, max_retries=3, timeout=60):
    """Fetch URL with retry logic using curl (IMF API blocks Python urllib)."""
    for attempt in range(max_retries):
        try:
            result = subprocess.run(
                ["curl", "-s", "-f", "--max-time", str(timeout),
                 "-H", "Accept: application/json", url],
                capture_output=True, timeout=timeout + 10
            )
            if result.returncode == 0 and result.stdout:
                return result.stdout
            err_msg = result.stderr.decode("utf-8", errors="replace").strip()
            raise RuntimeError(f"curl failed (rc={result.returncode}): {err_msg}")
        except (subprocess.TimeoutExpired, RuntimeError, OSError) as e:
            if attempt == max_retries - 1:
                raise
            wait = 2 ** attempt
            print(f"  Retry {attempt+1}/{max_retries} after {wait}s: {e}")
            time.sleep(wait)


def download_and_cache(indicator, cache_name):
    """Download IMF data for an indicator, cache locally, verify SHA256."""
    cache_path = os.path.join(CACHE_DIR, f"{cache_name}.json")
    sha_path = os.path.join(CACHE_DIR, f"{cache_name}.sha256")

    if os.path.exists(cache_path) and os.path.exists(sha_path):
        with open(cache_path, "rb") as f:
            data = f.read()
        actual_sha = hashlib.sha256(data).hexdigest()
        with open(sha_path, "r") as f:
            expected_sha = f.read().strip()
        if actual_sha == expected_sha:
            print(f"  Using cached {cache_name} (SHA256 verified: {actual_sha[:16]}...)")
            return json.loads(data)
        else:
            print(f"  Cache corrupted for {cache_name}, re-downloading...")

    url = f"{IMF_BASE}/{indicator}"
    print(f"  Downloading {indicator} from {url} ...")
    raw = fetch_with_retry(url)
    actual_sha = hashlib.sha256(raw).hexdigest()

    with open(cache_path, "wb") as f:
        f.write(raw)
    with open(sha_path, "w") as f:
        f.write(actual_sha)

    print(f"  Cached {cache_name} ({len(raw)} bytes, SHA256: {actual_sha[:16]}...)")
    return json.loads(raw)


def parse_imf_data(raw_json, indicator):
    """Parse IMF JSON into {country_code: {year_int: float_value}}."""
    values = raw_json.get("values", {}).get(indicator, {})
    result = {}
    for country, year_data in values.items():
        if country in COUNTRY_GROUPS:
            continue
        parsed = {}
        for year_str, val in year_data.items():
            try:
                yr = int(year_str)
                v = float(val)
                if YEAR_MIN <= yr <= YEAR_MAX and math.isfinite(v):
                    parsed[yr] = v
            except (ValueError, TypeError):
                continue
        if parsed:
            result[country] = parsed
    return result


def build_paired_dataset(debt_data, growth_data):
    """Build list of (debt_pct, gdp_growth) observations where both exist."""
    pairs = []
    countries_used = set()
    for country in debt_data:
        if country not in growth_data:
            continue
        for year in debt_data[country]:
            if year in growth_data[country]:
                pairs.append((debt_data[country][year], growth_data[country][year], country, year))
                countries_used.add(country)
    return pairs, countries_used


# ============================================================
# Statistical functions (stdlib only)
# ============================================================

def mean(values):
    """Compute arithmetic mean."""
    if not values:
        return float('nan')
    return sum(values) / len(values)


def std_dev(values):
    """Compute sample standard deviation."""
    if len(values) < 2:
        return float('nan')
    return statistics.stdev(values)


def cohens_d(group1, group2):
    """Compute Cohen's d effect size."""
    n1, n2 = len(group1), len(group2)
    if n1 < 2 or n2 < 2:
        return float('nan')
    m1, m2 = mean(group1), mean(group2)
    s1, s2 = std_dev(group1), std_dev(group2)
    pooled_s = math.sqrt(((n1-1)*s1**2 + (n2-1)*s2**2) / (n1+n2-2))
    if pooled_s == 0:
        return float('nan')
    return (m1 - m2) / pooled_s


def welch_t_stat(group1, group2):
    """Compute Welch's t-statistic."""
    n1, n2 = len(group1), len(group2)
    if n1 < 2 or n2 < 2:
        return float('nan'), float('nan')
    m1, m2 = mean(group1), mean(group2)
    s1, s2 = std_dev(group1), std_dev(group2)
    se = math.sqrt(s1**2/n1 + s2**2/n2)
    if se == 0:
        return float('nan'), float('nan')
    t = (m1 - m2) / se
    # Welch-Satterthwaite degrees of freedom
    num = (s1**2/n1 + s2**2/n2)**2
    den = (s1**2/n1)**2/(n1-1) + (s2**2/n2)**2/(n2-1)
    if den == 0:
        return t, float('nan')
    df = num / den
    return t, df


def permutation_test(below, above, n_perm, rng):
    """
    Two-sample permutation test for difference in means.
    H0: no difference in GDP growth between below/above threshold groups.
    Returns observed diff, p-value, and null distribution.
    """
    observed_diff = mean(below) - mean(above)
    combined = below + above
    n_below = len(below)
    count_extreme = 0
    null_diffs = []

    for _ in range(n_perm):
        rng.shuffle(combined)
        perm_below = combined[:n_below]
        perm_above = combined[n_below:]
        perm_diff = mean(perm_below) - mean(perm_above)
        null_diffs.append(perm_diff)
        if abs(perm_diff) >= abs(observed_diff):
            count_extreme += 1

    p_value = (count_extreme + 1) / (n_perm + 1)  # Conservative estimator
    return observed_diff, p_value, null_diffs


def bootstrap_ci(values, stat_fn, n_boot, ci_level, rng):
    """Bootstrap confidence interval for a statistic."""
    boot_stats = []
    n = len(values)
    for _ in range(n_boot):
        sample = [values[rng.randint(0, n-1)] for _ in range(n)]
        boot_stats.append(stat_fn(sample))
    boot_stats.sort()
    alpha = 1 - ci_level
    lo_idx = int(math.floor(alpha/2 * n_boot))
    hi_idx = int(math.ceil((1 - alpha/2) * n_boot)) - 1
    return boot_stats[lo_idx], boot_stats[hi_idx], boot_stats


def bootstrap_diff_ci(group1, group2, n_boot, ci_level, rng):
    """Bootstrap CI for difference in means between two groups."""
    diffs = []
    n1, n2 = len(group1), len(group2)
    for _ in range(n_boot):
        s1 = [group1[rng.randint(0, n1-1)] for _ in range(n1)]
        s2 = [group2[rng.randint(0, n2-1)] for _ in range(n2)]
        diffs.append(mean(s1) - mean(s2))
    diffs.sort()
    alpha = 1 - ci_level
    lo_idx = int(math.floor(alpha/2 * n_boot))
    hi_idx = int(math.ceil((1 - alpha/2) * n_boot)) - 1
    return diffs[lo_idx], diffs[hi_idx]


def holm_bonferroni(p_values_with_labels):
    """
    Holm-Bonferroni correction for multiple comparisons.
    Input: list of (label, p_value)
    Returns: list of (label, original_p, adjusted_p, significant_at_005)
    """
    m = len(p_values_with_labels)
    sorted_pvs = sorted(p_values_with_labels, key=lambda x: x[1])
    results = []
    max_adj_p = 0
    for i, (label, p) in enumerate(sorted_pvs):
        adj_p = min(p * (m - i), 1.0)
        adj_p = max(adj_p, max_adj_p)  # Enforce monotonicity
        max_adj_p = adj_p
        results.append((label, p, adj_p, adj_p < 0.05))
    # Sort back to original order by label
    results.sort(key=lambda x: x[0])
    return results


def spearman_rank(x_vals, y_vals):
    """Compute Spearman rank correlation."""
    n = len(x_vals)
    if n < 3:
        return float('nan')

    def rank_data(data):
        indexed = sorted(enumerate(data), key=lambda p: p[1])
        ranks = [0.0] * n
        i = 0
        while i < n:
            j = i
            while j < n - 1 and indexed[j+1][1] == indexed[j][1]:
                j += 1
            avg_rank = (i + j) / 2.0 + 1
            for k in range(i, j+1):
                ranks[indexed[k][0]] = avg_rank
            i = j + 1
        return ranks

    rx = rank_data(x_vals)
    ry = rank_data(y_vals)
    mean_rx = mean(rx)
    mean_ry = mean(ry)
    num = sum((a - mean_rx) * (b - mean_ry) for a, b in zip(rx, ry))
    den_x = math.sqrt(sum((a - mean_rx)**2 for a in rx))
    den_y = math.sqrt(sum((b - mean_ry)**2 for b in ry))
    if den_x == 0 or den_y == 0:
        return float('nan')
    return num / (den_x * den_y)


def fisher_z_ci(rho, n, ci_level=0.95):
    """Fisher z-transform confidence interval for correlation."""
    if n < 4 or abs(rho) >= 1:
        return float('nan'), float('nan')
    z = 0.5 * math.log((1 + rho) / (1 - rho))
    se = 1.0 / math.sqrt(n - 3)
    # z-critical for 95% CI
    z_crit = 1.96 if ci_level == 0.95 else 2.576
    lo_z = z - z_crit * se
    hi_z = z + z_crit * se
    lo = (math.exp(2*lo_z) - 1) / (math.exp(2*lo_z) + 1)
    hi = (math.exp(2*hi_z) - 1) / (math.exp(2*hi_z) + 1)
    return lo, hi


# ============================================================
# Sensitivity analyses
# ============================================================

def run_sensitivity_time_periods(pairs, rng):
    """Test if results hold across different time periods."""
    periods = [
        ("1990-2007 (pre-crisis)", 1990, 2007),
        ("2008-2014 (crisis era)", 2008, 2014),
        ("2015-2023 (post-crisis)", 2015, 2023),
    ]
    results = {}
    for pname, ymin, ymax in periods:
        sub = [(d, g) for d, g, c, y in pairs if ymin <= y <= ymax]
        if len(sub) < 20:
            results[pname] = {"n": len(sub), "note": "insufficient data"}
            continue
        debt_vals = [d for d, g in sub]
        growth_vals = [g for d, g in sub]
        rho = spearman_rank(debt_vals, growth_vals)
        # Test 90% threshold specifically
        below = [g for d, g in sub if d < 90]
        above = [g for d, g in sub if d >= 90]
        if len(below) >= 5 and len(above) >= 5:
            diff = mean(below) - mean(above)
            d_effect = cohens_d(below, above)
        else:
            diff = float('nan')
            d_effect = float('nan')
        results[pname] = {
            "n": len(sub),
            "spearman_rho": round(rho, 4),
            "mean_diff_at_90": round(diff, 3),
            "cohens_d_at_90": round(d_effect, 3)
        }
    return results


def run_sensitivity_exclude_outliers(pairs, rng):
    """Test excluding extreme GDP growth values (>|15%|)."""
    filtered = [(d, g, c, y) for d, g, c, y in pairs if abs(g) <= 15]
    n_excluded = len(pairs) - len(filtered)
    below = [g for d, g, c, y in filtered if d < 90]
    above = [g for d, g, c, y in filtered if d >= 90]
    if len(below) >= 5 and len(above) >= 5:
        diff, p, _ = permutation_test(below, above, 1000, rng)
        d_effect = cohens_d(below, above)
    else:
        diff, p, d_effect = float('nan'), float('nan'), float('nan')
    return {
        "n_original": len(pairs),
        "n_after_filter": len(filtered),
        "n_excluded": n_excluded,
        "mean_diff_at_90": round(diff, 3),
        "permutation_p": round(p, 4),
        "cohens_d": round(d_effect, 3)
    }


def run_sensitivity_income_groups(pairs, debt_data, growth_data, rng):
    """
    Split countries by median debt level and test within each group.
    This checks if results are driven by developing vs developed country confounding.
    """
    # Compute median debt per country
    country_med_debt = {}
    for country in debt_data:
        vals = list(debt_data[country].values())
        if vals:
            vals_sorted = sorted(vals)
            mid = len(vals_sorted) // 2
            country_med_debt[country] = vals_sorted[mid]

    if not country_med_debt:
        return {"note": "insufficient data"}

    overall_median = sorted(country_med_debt.values())[len(country_med_debt)//2]

    low_debt_countries = {c for c, m in country_med_debt.items() if m < overall_median}
    high_debt_countries = {c for c, m in country_med_debt.items() if m >= overall_median}

    results = {}
    for group_name, group_set in [("low-debt countries", low_debt_countries),
                                   ("high-debt countries", high_debt_countries)]:
        sub = [(d, g) for d, g, c, y in pairs if c in group_set]
        below = [g for d, g in sub if d < 90]
        above = [g for d, g in sub if d >= 90]
        if len(below) >= 5 and len(above) >= 5:
            diff = mean(below) - mean(above)
            d_effect = cohens_d(below, above)
        else:
            diff, d_effect = float('nan'), float('nan')
        results[group_name] = {
            "n_obs": len(sub),
            "n_countries": len(group_set),
            "mean_diff_at_90": round(diff, 3),
            "cohens_d_at_90": round(d_effect, 3)
        }
    results["median_debt_split_point"] = round(overall_median, 1)
    return results


def run_sensitivity_cluster_permutation(pairs, rng, threshold=90, n_perm=1000):
    """
    Cluster-level permutation test at a given threshold.
    Shuffles country LABELS (not individual observations) to respect
    within-country autocorrelation. This is a more conservative test.
    """
    # Group observations by country
    country_obs = {}
    for d, g, c, y in pairs:
        country_obs.setdefault(c, []).append((d, g))

    # Compute observed statistic: diff in mean growth below vs above threshold
    below = [g for d, g, c, y in pairs if d < threshold]
    above = [g for d, g, c, y in pairs if d >= threshold]
    if len(below) < 10 or len(above) < 10:
        return {"note": "insufficient data", "threshold": threshold}
    observed_diff = mean(below) - mean(above)

    # For each country, compute fraction of obs above threshold
    # Then shuffle country assignments to build null distribution
    countries = list(country_obs.keys())
    count_extreme = 0
    null_diffs = []

    for _ in range(n_perm):
        # Shuffle which country's data goes where
        shuffled_countries = countries[:]
        rng.shuffle(shuffled_countries)
        # Reassign: pair each country's obs with another country's debt values
        # Simpler approach: shuffle the country labels on the growth values
        perm_below = []
        perm_above = []
        for orig_c, shuf_c in zip(countries, shuffled_countries):
            orig_obs = country_obs[orig_c]  # (debt, growth) from original country
            shuf_obs = country_obs[shuf_c]  # we take growth from shuffled country
            # Match by index (both lists may differ in length, use min)
            for i in range(min(len(orig_obs), len(shuf_obs))):
                debt_val = orig_obs[i][0]  # Keep original debt
                growth_val = shuf_obs[i][1]  # Shuffled growth
                if debt_val < threshold:
                    perm_below.append(growth_val)
                else:
                    perm_above.append(growth_val)

        if perm_below and perm_above:
            perm_diff = mean(perm_below) - mean(perm_above)
            null_diffs.append(perm_diff)
            if abs(perm_diff) >= abs(observed_diff):
                count_extreme += 1

    p_value = (count_extreme + 1) / (len(null_diffs) + 1)
    return {
        "threshold": threshold,
        "observed_diff": round(observed_diff, 3),
        "cluster_permutation_p": round(p_value, 4),
        "n_countries": len(countries),
        "n_permutations": n_perm,
        "n_below": len(below),
        "n_above": len(above),
    }


# ============================================================
# Main analysis
# ============================================================

def run_analysis():
    print("=" * 70)
    print("IMF DEBT-GROWTH THRESHOLD RETEST")
    print("Permutation-based analysis with Holm-Bonferroni correction")
    print("=" * 70)

    os.makedirs(CACHE_DIR, exist_ok=True)
    rng = random.Random(SEED)

    # ----------------------------------------------------------
    print(f"\n[1/8] Downloading IMF data...")
    debt_raw = download_and_cache(DEBT_INDICATOR, "debt_gdp")
    growth_raw = download_and_cache(GROWTH_INDICATOR, "gdp_growth")

    # ----------------------------------------------------------
    print(f"\n[2/8] Parsing and merging data...")
    debt_data = parse_imf_data(debt_raw, DEBT_INDICATOR)
    growth_data = parse_imf_data(growth_raw, GROWTH_INDICATOR)
    pairs, countries_used = build_paired_dataset(debt_data, growth_data)
    print(f"  Observations (country-years): {len(pairs)}")
    print(f"  Countries with paired data: {len(countries_used)}")
    print(f"  Year range: {YEAR_MIN}-{YEAR_MAX}")

    if len(pairs) < 100:
        print("ERROR: Insufficient paired data points. Aborting.")
        sys.exit(1)

    debt_vals = [d for d, g, c, y in pairs]
    growth_vals = [g for d, g, c, y in pairs]

    # ----------------------------------------------------------
    print(f"\n[3/8] Overall correlation analysis...")
    rho = spearman_rank(debt_vals, growth_vals)
    rho_ci_lo, rho_ci_hi = fisher_z_ci(rho, len(pairs))
    print(f"  Spearman rho: {rho:.4f}")
    print(f"  95% CI (Fisher z): [{rho_ci_lo:.4f}, {rho_ci_hi:.4f}]")

    # Bootstrap CI for Spearman
    def spearman_boot(indices_unused):
        # This is called with a resampled list
        return spearman_rank(indices_unused, [growth_vals[i] for i in range(len(indices_unused))])

    # Direct bootstrap of Spearman correlation
    rho_boot_lo, rho_boot_hi, _ = bootstrap_ci(
        list(range(len(pairs))),
        lambda idxs: spearman_rank([debt_vals[i] for i in idxs], [growth_vals[i] for i in idxs]),
        N_BOOTSTRAP, BOOTSTRAP_CI_LEVEL, rng
    )
    print(f"  95% CI (Bootstrap, {N_BOOTSTRAP} resamples): [{rho_boot_lo:.4f}, {rho_boot_hi:.4f}]")

    # ----------------------------------------------------------
    print(f"\n[4/8] Threshold scan: permutation tests at {len(THRESHOLDS)} thresholds...")
    print(f"  Thresholds: {THRESHOLDS[0]}% to {THRESHOLDS[-1]}% (step 5%)")
    print(f"  Permutations per threshold: {N_PERMUTATIONS}")

    threshold_results = []
    p_values_for_correction = []

    for thresh in THRESHOLDS:
        below = [g for d, g, c, y in pairs if d < thresh]
        above = [g for d, g, c, y in pairs if d >= thresh]

        if len(below) < 10 or len(above) < 10:
            threshold_results.append({
                "threshold": thresh,
                "note": "insufficient observations in one group",
                "n_below": len(below),
                "n_above": len(above)
            })
            continue

        # Permutation test
        obs_diff, perm_p, null_dist = permutation_test(below, above, N_PERMUTATIONS, rng)

        # Effect size
        d_effect = cohens_d(below, above)

        # Welch t-test
        t_stat, t_df = welch_t_stat(below, above)

        # Bootstrap CI for difference in means
        diff_ci_lo, diff_ci_hi = bootstrap_diff_ci(below, above, N_BOOTSTRAP, BOOTSTRAP_CI_LEVEL, rng)

        result = {
            "threshold": thresh,
            "n_below": len(below),
            "n_above": len(above),
            "mean_below": round(mean(below), 3),
            "mean_above": round(mean(above), 3),
            "observed_diff": round(obs_diff, 3),
            "permutation_p": round(perm_p, 4),
            "cohens_d": round(d_effect, 3),
            "welch_t": round(t_stat, 3),
            "welch_df": round(t_df, 1),
            "bootstrap_ci_lo": round(diff_ci_lo, 3),
            "bootstrap_ci_hi": round(diff_ci_hi, 3),
        }
        threshold_results.append(result)
        p_values_for_correction.append((thresh, perm_p))

        status = "*" if perm_p < 0.05 else " "
        print(f"  {status} {thresh:>3d}%: diff={obs_diff:+.2f}pp, p={perm_p:.4f}, d={d_effect:+.3f}, "
              f"95%CI=[{diff_ci_lo:+.2f}, {diff_ci_hi:+.2f}], n=({len(below)},{len(above)})")

    # ----------------------------------------------------------
    print(f"\n[5/8] Multiple comparison correction (Holm-Bonferroni)...")
    corrected = holm_bonferroni(p_values_for_correction)
    n_sig_raw = sum(1 for _, p in p_values_for_correction if p < 0.05)
    n_sig_corrected = sum(1 for _, _, adj_p, sig in corrected if sig)
    print(f"  Thresholds significant at p<0.05 (uncorrected): {n_sig_raw}/{len(p_values_for_correction)}")
    print(f"  Thresholds significant at p<0.05 (Holm-Bonferroni): {n_sig_corrected}/{len(p_values_for_correction)}")

    # Add corrected p-values to threshold results
    corrected_dict = {label: (orig_p, adj_p, sig) for label, orig_p, adj_p, sig in corrected}
    for tr in threshold_results:
        thresh = tr["threshold"]
        if thresh in corrected_dict:
            tr["adjusted_p"] = round(corrected_dict[thresh][1], 4)
            tr["significant_after_correction"] = corrected_dict[thresh][2]

    # Find threshold with largest effect size (not lowest p, since all hit floor)
    valid_results = [r for r in threshold_results if "cohens_d" in r]
    if valid_results:
        best = max(valid_results, key=lambda r: abs(r["cohens_d"]))
        print(f"\n  Largest effect size at threshold: {best['threshold']}%")
        print(f"    Cohen's d = {best['cohens_d']:+.3f}, Diff = {best['observed_diff']:+.2f} pp")
        print(f"    Uncorrected p = {best['permutation_p']:.4f}, Adjusted p = {best.get('adjusted_p', 'N/A')}")

    # ----------------------------------------------------------
    print(f"\n[6/8] Sensitivity analysis: time periods...")
    sensitivity_time = run_sensitivity_time_periods(pairs, rng)
    for period, res in sensitivity_time.items():
        if "note" in res:
            print(f"  {period}: {res['note']} (n={res['n']})")
        else:
            print(f"  {period}: n={res['n']}, rho={res['spearman_rho']:.3f}, "
                  f"diff@90={res['mean_diff_at_90']:+.2f}pp, d={res['cohens_d_at_90']:+.3f}")

    print(f"\n  Sensitivity: excluding outliers (|growth| > 15%)...")
    sensitivity_outliers = run_sensitivity_exclude_outliers(pairs, rng)
    print(f"    Excluded: {sensitivity_outliers['n_excluded']} obs")
    print(f"    Diff@90 = {sensitivity_outliers['mean_diff_at_90']:+.2f}pp, "
          f"p = {sensitivity_outliers['permutation_p']:.4f}, d = {sensitivity_outliers['cohens_d']:+.3f}")

    print(f"\n  Sensitivity: income groups (by median debt level)...")
    sensitivity_income = run_sensitivity_income_groups(pairs, debt_data, growth_data, rng)
    for group, res in sensitivity_income.items():
        if isinstance(res, dict) and "n_obs" in res:
            print(f"    {group}: n={res['n_obs']}, diff@90={res['mean_diff_at_90']:+.2f}pp, "
                  f"d={res['cohens_d_at_90']:+.3f}")
        elif group == "median_debt_split_point":
            print(f"    Split point: {res}% debt/GDP")

    print(f"\n  Sensitivity: cluster-level permutation (shuffling countries, not obs)...")
    cluster_results = {}
    for ct in [60, 90, 120]:
        cr = run_sensitivity_cluster_permutation(pairs, rng, threshold=ct, n_perm=1000)
        cluster_results[ct] = cr
        if "cluster_permutation_p" in cr:
            print(f"    Threshold {ct}%: diff={cr['observed_diff']:+.2f}pp, "
                  f"cluster p={cr['cluster_permutation_p']:.4f} (n_countries={cr['n_countries']})")

    # ----------------------------------------------------------
    print(f"\n[7/8] Descriptive statistics by debt quartile...")
    debt_sorted = sorted(debt_vals)
    q25 = debt_sorted[len(debt_sorted)//4]
    q50 = debt_sorted[len(debt_sorted)//2]
    q75 = debt_sorted[3*len(debt_sorted)//4]

    quartiles = [
        (f"Q1 (<{q25:.0f}%)", [g for d, g, c, y in pairs if d < q25]),
        (f"Q2 ({q25:.0f}-{q50:.0f}%)", [g for d, g, c, y in pairs if q25 <= d < q50]),
        (f"Q3 ({q50:.0f}-{q75:.0f}%)", [g for d, g, c, y in pairs if q50 <= d < q75]),
        (f"Q4 (>={q75:.0f}%)", [g for d, g, c, y in pairs if d >= q75]),
    ]

    for qname, qvals in quartiles:
        if qvals:
            m = mean(qvals)
            s = std_dev(qvals) if len(qvals) > 1 else 0
            lo, hi, _ = bootstrap_ci(qvals, mean, N_BOOTSTRAP, BOOTSTRAP_CI_LEVEL, rng)
            print(f"  {qname}: mean={m:.2f}%, sd={s:.2f}, 95%CI=[{lo:.2f}, {hi:.2f}], n={len(qvals)}")

    # ----------------------------------------------------------
    print(f"\n[8/8] Writing results...")

    results = {
        "metadata": {
            "analysis": "IMF Debt-Growth Threshold Retest",
            "version": "1.0.0",
            "author": "Claw, David Austin",
            "seed": SEED,
            "n_permutations": N_PERMUTATIONS,
            "n_bootstrap": N_BOOTSTRAP,
            "year_range": [YEAR_MIN, YEAR_MAX],
            "thresholds_tested": THRESHOLDS,
            "data_source": "IMF World Economic Outlook via datamapper API",
            "debt_indicator": DEBT_INDICATOR,
            "growth_indicator": GROWTH_INDICATOR,
        },
        "data_summary": {
            "n_observations": len(pairs),
            "n_countries": len(countries_used),
            "debt_mean": round(mean(debt_vals), 2),
            "debt_median": round(sorted(debt_vals)[len(debt_vals)//2], 2),
            "debt_min": round(min(debt_vals), 2),
            "debt_max": round(max(debt_vals), 2),
            "growth_mean": round(mean(growth_vals), 2),
            "growth_median": round(sorted(growth_vals)[len(growth_vals)//2], 2),
            "growth_min": round(min(growth_vals), 2),
            "growth_max": round(max(growth_vals), 2),
        },
        "overall_correlation": {
            "spearman_rho": round(rho, 4),
            "fisher_z_ci_95": [round(rho_ci_lo, 4), round(rho_ci_hi, 4)],
            "bootstrap_ci_95": [round(rho_boot_lo, 4), round(rho_boot_hi, 4)],
        },
        "threshold_scan": threshold_results,
        "multiple_comparison": {
            "method": "Holm-Bonferroni",
            "n_tests": len(p_values_for_correction),
            "n_significant_uncorrected": n_sig_raw,
            "n_significant_corrected": n_sig_corrected,
        },
        "sensitivity": {
            "time_periods": sensitivity_time,
            "exclude_outliers": sensitivity_outliers,
            "income_groups": sensitivity_income,
            "cluster_permutation": cluster_results,
        },
        "quartile_analysis": {
            qname: {
                "mean_growth": round(mean(qvals), 3),
                "sd": round(std_dev(qvals), 3) if len(qvals) > 1 else None,
                "n": len(qvals)
            }
            for qname, qvals in quartiles if qvals
        }
    }

    results_path = os.path.join(RESULTS_DIR, "results.json")
    with open(results_path, "w") as f:
        json.dump(results, f, indent=2)
    print(f"  Wrote {results_path}")

    # Write report.md
    report_path = os.path.join(RESULTS_DIR, "report.md")
    with open(report_path, "w") as f:
        f.write("# IMF Debt-Growth Threshold Retest: Results Report\n\n")
        f.write(f"**Data:** {len(pairs)} country-year observations from {len(countries_used)} countries ({YEAR_MIN}-{YEAR_MAX})\n\n")
        f.write(f"**Source:** IMF World Economic Outlook ({DEBT_INDICATOR}, {GROWTH_INDICATOR})\n\n")

        f.write("## Overall Correlation\n\n")
        f.write(f"Spearman rho = {rho:.4f}, 95% Fisher CI [{rho_ci_lo:.4f}, {rho_ci_hi:.4f}], "
                f"Bootstrap CI [{rho_boot_lo:.4f}, {rho_boot_hi:.4f}]\n\n")

        f.write("## Threshold Scan Results\n\n")
        f.write("| Threshold | n_below | n_above | Diff (pp) | Perm p | Adj p | Cohen's d | 95% CI |\n")
        f.write("|-----------|---------|---------|-----------|--------|-------|-----------|--------|\n")
        for tr in threshold_results:
            if "permutation_p" in tr:
                sig_mark = "**" if tr.get("significant_after_correction", False) else ""
                f.write(f"| {sig_mark}{tr['threshold']}%{sig_mark} | {tr['n_below']} | {tr['n_above']} | "
                        f"{tr['observed_diff']:+.2f} | {tr['permutation_p']:.4f} | "
                        f"{tr.get('adjusted_p', 'N/A')} | {tr['cohens_d']:+.3f} | "
                        f"[{tr['bootstrap_ci_lo']:+.2f}, {tr['bootstrap_ci_hi']:+.2f}] |\n")

        f.write("\n## Multiple Comparison Correction\n\n")
        f.write(f"- Method: Holm-Bonferroni\n")
        f.write(f"- Tests: {len(p_values_for_correction)}\n")
        f.write(f"- Significant (uncorrected p<0.05): {n_sig_raw}\n")
        f.write(f"- Significant (corrected p<0.05): {n_sig_corrected}\n\n")

        f.write("## Sensitivity Analysis\n\n")
        f.write("### Time Periods\n\n")
        for period, res in sensitivity_time.items():
            if "note" not in res:
                f.write(f"- {period}: n={res['n']}, rho={res['spearman_rho']}, "
                        f"diff@90={res['mean_diff_at_90']:+.2f}pp, d={res['cohens_d_at_90']:+.3f}\n")

        f.write("\n### Outlier Exclusion\n\n")
        f.write(f"Excluding |growth|>15%: n={sensitivity_outliers['n_after_filter']}, "
                f"diff@90={sensitivity_outliers['mean_diff_at_90']:+.2f}pp, "
                f"p={sensitivity_outliers['permutation_p']:.4f}\n\n")

        f.write("### Country Income Groups\n\n")
        for group, res in sensitivity_income.items():
            if isinstance(res, dict) and "n_obs" in res:
                f.write(f"- {group}: n={res['n_obs']}, diff@90={res['mean_diff_at_90']:+.2f}pp, "
                        f"d={res['cohens_d_at_90']:+.3f}\n")

    print(f"  Wrote {report_path}")

    print("\n" + "=" * 70)
    print("ANALYSIS COMPLETE")
    print("=" * 70)

    return results


# ============================================================
# Verification mode
# ============================================================

def verify():
    """Machine-checkable assertions on results.json."""
    results_path = os.path.join(RESULTS_DIR, "results.json")
    print(f"Verifying {results_path}...")

    with open(results_path) as f:
        results = json.load(f)

    checks = []

    # Check 1: Sufficient data
    n_obs = results["data_summary"]["n_observations"]
    checks.append(("n_observations >= 500", n_obs >= 500, f"n={n_obs}"))

    # Check 2: Sufficient countries
    n_countries = results["data_summary"]["n_countries"]
    checks.append(("n_countries >= 50", n_countries >= 50, f"n={n_countries}"))

    # Check 3: All thresholds tested
    n_thresh = len(results["threshold_scan"])
    checks.append(("all 23 thresholds tested", n_thresh == 23, f"n={n_thresh}"))

    # Check 4: Correlation is in valid range
    rho = results["overall_correlation"]["spearman_rho"]
    checks.append(("spearman_rho in [-1, 1]", -1 <= rho <= 1, f"rho={rho}"))

    # Check 5: CI is properly ordered
    ci = results["overall_correlation"]["fisher_z_ci_95"]
    checks.append(("fisher_CI_lo < fisher_CI_hi", ci[0] < ci[1], f"CI={ci}"))

    # Check 6: Bootstrap CI is properly ordered
    bci = results["overall_correlation"]["bootstrap_ci_95"]
    checks.append(("bootstrap_CI_lo < bootstrap_CI_hi", bci[0] < bci[1], f"CI={bci}"))

    # Check 7: Multiple comparison correction present
    mc = results["multiple_comparison"]
    checks.append(("holm_bonferroni applied", mc["method"] == "Holm-Bonferroni", f"method={mc['method']}"))

    # Check 8: corrected <= uncorrected significant
    checks.append(("corrected_sig <= uncorrected_sig",
                    mc["n_significant_corrected"] <= mc["n_significant_uncorrected"],
                    f"{mc['n_significant_corrected']} <= {mc['n_significant_uncorrected']}"))

    # Check 9: Sensitivity analyses present
    checks.append(("time_period sensitivity present",
                    len(results["sensitivity"]["time_periods"]) >= 2,
                    f"n={len(results['sensitivity']['time_periods'])}"))

    # Check 10: Outlier sensitivity present
    checks.append(("outlier sensitivity present",
                    results["sensitivity"]["exclude_outliers"]["n_after_filter"] > 0,
                    f"n={results['sensitivity']['exclude_outliers']['n_after_filter']}"))

    # Check 11: Income group sensitivity present
    checks.append(("income_group sensitivity present",
                    "median_debt_split_point" in results["sensitivity"]["income_groups"],
                    str(list(results["sensitivity"]["income_groups"].keys()))))

    # Check 12: Cluster permutation sensitivity present
    checks.append(("cluster_permutation sensitivity present",
                    "cluster_permutation" in results["sensitivity"],
                    str(list(results["sensitivity"].keys()))))

    # Check 13: results.json has all required sections
    required_sections = ["metadata", "data_summary", "overall_correlation",
                         "threshold_scan", "multiple_comparison", "sensitivity", "quartile_analysis"]
    all_present = all(s in results for s in required_sections)
    checks.append(("all result sections present", all_present, str(list(results.keys()))))

    # Check 14: Permutation count matches config
    first_valid = next((r for r in results["threshold_scan"] if "permutation_p" in r), None)
    checks.append(("permutation_count matches config",
                    results["metadata"]["n_permutations"] == 3000,
                    f"n_perm={results['metadata']['n_permutations']}"))

    # Check 14: p-values in valid range
    all_p_valid = all(
        0 <= r.get("permutation_p", 0) <= 1
        for r in results["threshold_scan"]
        if "permutation_p" in r
    )
    checks.append(("all p-values in [0,1]", all_p_valid, ""))

    # Check 15: report.md exists
    report_exists = os.path.exists(os.path.join(RESULTS_DIR, "report.md"))
    checks.append(("report.md exists", report_exists, ""))

    print("\nVerification Results:")
    print("-" * 60)
    all_pass = True
    for name, passed, detail in checks:
        status = "PASS" if passed else "FAIL"
        if not passed:
            all_pass = False
        print(f"  [{status}] {name}: {detail}")

    print("-" * 60)
    if all_pass:
        print(f"ALL {len(checks)} CHECKS PASSED")
    else:
        n_fail = sum(1 for _, p, _ in checks if not p)
        print(f"FAILED: {n_fail}/{len(checks)} checks failed")
        sys.exit(1)


# ============================================================
# Entry point
# ============================================================

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="IMF Debt-Growth Threshold Retest")
    parser.add_argument("--verify", action="store_true", help="Run verification checks on results.json")
    args = parser.parse_args()

    if args.verify:
        verify()
    else:
        run_analysis()
SCRIPT_EOF
```

**Expected output:** Script file created at `/tmp/claw4s_auto_imf-debt-growth-threshold/analysis.py`, no errors.

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_imf-debt-growth-threshold && python3 analysis.py
```

**Expected output:**
- Sections `[1/8]` through `[8/8]` printed to stdout
- Ends with `ANALYSIS COMPLETE`
- Creates `results.json` and `report.md` in the workspace
- Runtime: 5-15 minutes (permutation tests are computationally intensive)

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_imf-debt-growth-threshold && python3 analysis.py --verify
```

**Expected output:**
- All 15 verification checks show `[PASS]`
- Ends with `ALL 16 CHECKS PASSED`

## Success Criteria

- [ ] `results.json` exists with all sections: metadata, data_summary, overall_correlation, threshold_scan, multiple_comparison, sensitivity, quartile_analysis
- [ ] `report.md` exists with threshold scan table and sensitivity results
- [ ] All 15 verification checks pass
- [ ] Data cached locally with SHA256 verification
- [ ] All random operations seeded (seed=42)
- [ ] No external dependencies beyond Python 3.8+ stdlib

## Failure Conditions

- IMF API unreachable → retry with exponential backoff (3 attempts)
- Fewer than 500 paired observations → abort with error message
- Any verification check fails → report which check and why
- Script crashes → check stderr for traceback