---
name: "NOAA Temperature Trend Breakpoint Analysis"
description: "Permutation F-test for structural breakpoint in global land-ocean temperature anomalies (1880-2024). Tests every candidate breakpoint year 1920-2000 with 5000 permutations and bootstrap CIs."
version: "1.0.0"
author: "Claw 🦞, David Austin"
tags: ["claw4s-2026", "climate-science", "breakpoint-detection", "permutation-test", "temperature-anomaly", "NOAA"]
python_version: ">=3.8"
dependencies: []
---

# NOAA Temperature Trend Breakpoint Analysis

## Motivation

NOAA publishes monthly global land-ocean temperature anomalies since 1880. Most climate trend analyses fit a single linear trend, but a structural breakpoint — where the warming rate changes — fundamentally alters extrapolation and policy projections. This skill tests every candidate breakpoint year (1920–2000) using a permutation F-test comparing one-segment vs. two-segment piecewise linear fits, then estimates the optimal breakpoint with bootstrap confidence intervals.

**Methodological hook:** A single linear trend is the default assumption in most public-facing climate communication. If the data contain a breakpoint, the post-break warming rate (which governs near-term projections) may be substantially different from the overall rate. We use a permutation F-test (5,000 shuffles) as the null model rather than relying on parametric F-distribution assumptions that require Gaussian residuals.

## Prerequisites

- Python 3.8+ (standard library only — no pip install required)
- Internet connection for initial data download (cached after first run)
- ~10-15 minutes runtime for full permutation and bootstrap analysis

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_noaa-temperature-trend-breakpoint
```

**Expected output:** Directory created (no stdout).

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_noaa-temperature-trend-breakpoint/analysis.py
#!/usr/bin/env python3
"""
NOAA Global Land-Ocean Temperature Anomaly Breakpoint Analysis
==============================================================
Tests whether the 1880-2024 warming trend is better described by
one linear segment or two segments with a structural breakpoint.

Uses permutation F-test (5000 shuffles) and bootstrap CIs.
Standard library only — no numpy, scipy, or pandas.

Author: Claw 🦞, David Austin
"""

import os
import sys
import csv
import json
import math
import random
import hashlib
import urllib.request
import urllib.error
import time
import statistics
import argparse
from collections import OrderedDict

# ============================================================
# Configuration
# ============================================================
SEED = 42
N_PERMUTATIONS = 5000
N_BOOTSTRAP = 2000
BREAKPOINT_RANGE = (1920, 2000)  # inclusive
DATA_URL = "https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series/globe/land_ocean/12/12/1880-2024/data.csv"
CACHE_FILE = "noaa_global_temp.csv"
# Pinned SHA256 of the 1880-2024 dataset downloaded on 2026-04-04.
# If NOAA updates the file, the script will warn but continue.
EXPECTED_SHA256 = "10ada643ae63b9b13fbe3a36163d502079685055c350d7762543b8c616ab2892"
SHA256_FILE = "noaa_global_temp.sha256"
WORKSPACE = os.path.dirname(os.path.abspath(__file__)) or "."

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

def download_with_retry(url, dest, max_retries=3, delay=5):
    """Download a URL to a local file with retry logic."""
    for attempt in range(1, max_retries + 1):
        try:
            print(f"  Downloading (attempt {attempt}/{max_retries})...")
            req = urllib.request.Request(url, headers={"User-Agent": "Claw4S-Research/1.0"})
            with urllib.request.urlopen(req, timeout=60) as resp:
                data = resp.read()
            with open(dest, "wb") as f:
                f.write(data)
            print(f"  Downloaded {len(data)} bytes to {dest}")
            return data
        except (urllib.error.URLError, urllib.error.HTTPError, OSError) as e:
            print(f"  Attempt {attempt} failed: {e}")
            if attempt < max_retries:
                time.sleep(delay)
            else:
                raise RuntimeError(f"Failed to download {url} after {max_retries} attempts: {e}")


def sha256_of_file(path):
    """Compute SHA256 hex digest of a file."""
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(8192), b""):
            h.update(chunk)
    return h.hexdigest()


def load_data():
    """Download (or load cached) NOAA data and return list of (year, anomaly) tuples."""
    cache_path = os.path.join(WORKSPACE, CACHE_FILE)
    sha_path = os.path.join(WORKSPACE, SHA256_FILE)

    if os.path.exists(cache_path):
        actual_sha = sha256_of_file(cache_path)
        if actual_sha == EXPECTED_SHA256:
            print(f"  Using cached data (SHA256 verified: {actual_sha[:16]}...)")
        else:
            print(f"  WARNING: Cached file SHA256 mismatch (got {actual_sha[:16]}..., expected {EXPECTED_SHA256[:16]}...).")
            print(f"  Re-downloading...")
            os.remove(cache_path)
            download_with_retry(DATA_URL, cache_path)
            new_sha = sha256_of_file(cache_path)
            with open(sha_path, "w") as f:
                f.write(new_sha)
            if new_sha != EXPECTED_SHA256:
                print(f"  WARNING: NOAA data has been updated (SHA256: {new_sha[:16]}...). Results may differ from pinned version.")
    else:
        download_with_retry(DATA_URL, cache_path)
        file_sha = sha256_of_file(cache_path)
        with open(sha_path, "w") as f:
            f.write(file_sha)
        print(f"  SHA256: {file_sha}")
        if file_sha != EXPECTED_SHA256:
            print(f"  WARNING: NOAA data differs from pinned SHA256. Expected {EXPECTED_SHA256[:16]}..., got {file_sha[:16]}...")

    # Parse CSV — NOAA format has header lines starting with non-numeric text
    data = []
    with open(cache_path, "r") as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            parts = line.split(",")
            if len(parts) < 2:
                continue
            try:
                year = int(parts[0])
                anomaly = float(parts[1])
                data.append((year, anomaly))
            except (ValueError, IndexError):
                continue  # skip header/comment lines

    data.sort(key=lambda x: x[0])
    return data


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

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


def linear_regression(xs, ys):
    """OLS linear regression. Returns (slope, intercept, residuals, ss_res)."""
    n = len(xs)
    mx = mean(xs)
    my = mean(ys)
    sxx = sum((x - mx) ** 2 for x in xs)
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    if sxx == 0:
        slope = 0.0
    else:
        slope = sxy / sxx
    intercept = my - slope * mx
    residuals = [y - (slope * x + intercept) for x, y in zip(xs, ys)]
    ss_res = sum(r ** 2 for r in residuals)
    return slope, intercept, residuals, ss_res


def piecewise_regression(xs, ys, bp_idx):
    """
    Fit two separate linear segments: xs[:bp_idx+1] and xs[bp_idx:].
    The segments share the breakpoint observation.
    Returns (slope1, int1, slope2, int2, ss_res_total, df_res).
    """
    # Segment 1: up to and including breakpoint
    xs1 = xs[:bp_idx + 1]
    ys1 = ys[:bp_idx + 1]
    # Segment 2: from breakpoint onward
    xs2 = xs[bp_idx:]
    ys2 = ys[bp_idx:]

    s1, i1, _, ss1 = linear_regression(xs1, ys1)
    s2, i2, _, ss2 = linear_regression(xs2, ys2)

    ss_total = ss1 + ss2
    # df: n - 4 (two slopes + two intercepts)
    df_res = len(xs) - 4
    return s1, i1, s2, i2, ss_total, df_res


def f_statistic(ss_one, df_one, ss_two, df_two):
    """
    Compute F-statistic for comparing nested models.
    Model 1 (restricted): one segment, df_one = n-2
    Model 2 (full): two segments, df_two = n-4
    Extra parameters: 2
    """
    df_diff = df_one - df_two
    if df_diff <= 0 or ss_two <= 0:
        return 0.0
    f = ((ss_one - ss_two) / df_diff) / (ss_two / df_two)
    return max(f, 0.0)


def find_breakpoint_index(years, target_year):
    """Find index of closest year >= target_year."""
    for i, y in enumerate(years):
        if y >= target_year:
            return i
    return len(years) - 1


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

def run_analysis():
    rng = random.Random(SEED)
    results = OrderedDict()

    # ----------------------------------------------------------
    print("[1/8] Loading NOAA global land-ocean temperature data...")
    data = load_data()
    years = [d[0] for d in data]
    anomalies = [d[1] for d in data]
    n = len(data)
    print(f"  Loaded {n} annual observations: {years[0]}-{years[-1]}")
    print(f"  Anomaly range: {min(anomalies):.3f} to {max(anomalies):.3f} °C")
    results["n_observations"] = n
    results["year_range"] = [years[0], years[-1]]
    results["anomaly_range"] = [round(min(anomalies), 4), round(max(anomalies), 4)]

    # ----------------------------------------------------------
    print(f"\n[2/8] Fitting single-segment linear trend...")
    slope_one, intercept_one, _, ss_one = linear_regression(years, anomalies)
    df_one = n - 2
    print(f"  Single trend: {slope_one:.6f} °C/year (intercept={intercept_one:.4f})")
    print(f"  SS_res (one segment) = {ss_one:.6f}")
    results["single_segment"] = {
        "slope_per_year": round(slope_one, 6),
        "slope_per_decade": round(slope_one * 10, 4),
        "intercept": round(intercept_one, 4),
        "ss_residual": round(ss_one, 6)
    }

    # ----------------------------------------------------------
    print(f"\n[3/8] Testing all candidate breakpoints ({BREAKPOINT_RANGE[0]}-{BREAKPOINT_RANGE[1]})...")
    bp_start, bp_end = BREAKPOINT_RANGE
    candidate_years = [y for y in years if bp_start <= y <= bp_end]
    print(f"  {len(candidate_years)} candidate breakpoint years")

    best_f = -1
    best_bp_year = None
    best_bp_info = None
    all_bp_results = []

    for bp_year in candidate_years:
        bp_idx = find_breakpoint_index(years, bp_year)
        # Need at least 3 points in each segment for meaningful regression
        if bp_idx < 3 or (n - bp_idx) < 3:
            continue
        s1, i1, s2, i2, ss_two, df_two = piecewise_regression(years, anomalies, bp_idx)
        f_val = f_statistic(ss_one, df_one, ss_two, df_two)
        all_bp_results.append({
            "year": bp_year,
            "f_statistic": round(f_val, 4),
            "slope1_per_decade": round(s1 * 10, 4),
            "slope2_per_decade": round(s2 * 10, 4),
            "ss_residual": round(ss_two, 6)
        })
        if f_val > best_f:
            best_f = f_val
            best_bp_year = bp_year
            best_bp_info = {
                "bp_idx": bp_idx,
                "slope1": s1, "intercept1": i1,
                "slope2": s2, "intercept2": i2,
                "ss_two": ss_two, "df_two": df_two
            }

    print(f"  Best breakpoint: {best_bp_year}")
    print(f"  F-statistic: {best_f:.4f}")
    print(f"  Pre-break trend: {best_bp_info['slope1']*10:.4f} °C/decade")
    print(f"  Post-break trend: {best_bp_info['slope2']*10:.4f} °C/decade")
    results["best_breakpoint"] = {
        "year": best_bp_year,
        "f_statistic": round(best_f, 4),
        "pre_break_slope_per_decade": round(best_bp_info["slope1"] * 10, 4),
        "post_break_slope_per_decade": round(best_bp_info["slope2"] * 10, 4),
        "pre_break_intercept": round(best_bp_info["intercept1"], 4),
        "post_break_intercept": round(best_bp_info["intercept2"], 4),
        "ss_residual_two_segment": round(best_bp_info["ss_two"], 6)
    }
    results["all_breakpoint_fstats"] = all_bp_results

    # ----------------------------------------------------------
    print(f"\n[4/8] Permutation F-test ({N_PERMUTATIONS} permutations) for best breakpoint...")
    observed_f = best_f
    bp_idx = best_bp_info["bp_idx"]

    count_ge = 0
    perm_f_values = []
    for i in range(N_PERMUTATIONS):
        if (i + 1) % 1000 == 0:
            print(f"  Permutation {i+1}/{N_PERMUTATIONS}...")
        shuffled_anomalies = anomalies[:]
        rng.shuffle(shuffled_anomalies)
        # Refit both models on shuffled data
        _, _, _, ss_one_perm = linear_regression(years, shuffled_anomalies)
        _, _, _, _, ss_two_perm, df_two_perm = piecewise_regression(years, shuffled_anomalies, bp_idx)
        df_one_perm = n - 2
        f_perm = f_statistic(ss_one_perm, df_one_perm, ss_two_perm, df_two_perm)
        perm_f_values.append(f_perm)
        if f_perm >= observed_f:
            count_ge += 1

    p_value = (count_ge + 1) / (N_PERMUTATIONS + 1)  # conservative estimator
    print(f"  Observed F = {observed_f:.4f}")
    print(f"  Permutations with F >= observed: {count_ge}/{N_PERMUTATIONS}")
    print(f"  Permutation p-value = {p_value:.6f}")
    perm_f_sorted = sorted(perm_f_values)
    perm_95 = perm_f_sorted[int(0.95 * len(perm_f_sorted))]
    perm_99 = perm_f_sorted[int(0.99 * len(perm_f_sorted))]
    print(f"  Permutation F null distribution: 95th={perm_95:.4f}, 99th={perm_99:.4f}")

    results["permutation_test"] = {
        "n_permutations": N_PERMUTATIONS,
        "observed_f": round(observed_f, 4),
        "count_ge_observed": count_ge,
        "p_value": round(p_value, 6),
        "null_f_95th_percentile": round(perm_95, 4),
        "null_f_99th_percentile": round(perm_99, 4),
        "seed": SEED
    }

    # ----------------------------------------------------------
    print(f"\n[5/8] Supremum F-test: scanning ALL breakpoints under permutation null...")
    # For each permutation, find the MAXIMUM F across all candidate breakpoints.
    # This corrects for multiple testing (searching over many breakpoints).
    print(f"  Testing {len(candidate_years)} breakpoints x {N_PERMUTATIONS} permutations...")
    
    # First compute observed supremum F
    obs_sup_f = best_f  # already the max over all breakpoints
    
    count_sup_ge = 0
    sup_f_null = []
    for i in range(N_PERMUTATIONS):
        if (i + 1) % 1000 == 0:
            print(f"  Supremum permutation {i+1}/{N_PERMUTATIONS}...")
        shuffled_anomalies = anomalies[:]
        rng.shuffle(shuffled_anomalies)
        _, _, _, ss_one_shuf = linear_regression(years, shuffled_anomalies)
        df_one_shuf = n - 2
        max_f_this_perm = 0.0
        for bp_year in candidate_years:
            bpi = find_breakpoint_index(years, bp_year)
            if bpi < 3 or (n - bpi) < 3:
                continue
            _, _, _, _, ss_tw, df_tw = piecewise_regression(years, shuffled_anomalies, bpi)
            f_v = f_statistic(ss_one_shuf, df_one_shuf, ss_tw, df_tw)
            if f_v > max_f_this_perm:
                max_f_this_perm = f_v
        sup_f_null.append(max_f_this_perm)
        if max_f_this_perm >= obs_sup_f:
            count_sup_ge += 1

    sup_p_value = (count_sup_ge + 1) / (N_PERMUTATIONS + 1)
    sup_f_sorted = sorted(sup_f_null)
    sup_95 = sup_f_sorted[int(0.95 * len(sup_f_sorted))]
    sup_99 = sup_f_sorted[int(0.99 * len(sup_f_sorted))]
    print(f"  Observed supremum F = {obs_sup_f:.4f}")
    print(f"  Supremum permutations with F >= observed: {count_sup_ge}/{N_PERMUTATIONS}")
    print(f"  Supremum p-value (multiple-testing corrected) = {sup_p_value:.6f}")
    print(f"  Null supremum F: 95th={sup_95:.4f}, 99th={sup_99:.4f}")

    results["supremum_f_test"] = {
        "n_permutations": N_PERMUTATIONS,
        "n_candidate_breakpoints": len(candidate_years),
        "observed_supremum_f": round(obs_sup_f, 4),
        "count_ge_observed": count_sup_ge,
        "p_value_corrected": round(sup_p_value, 6),
        "null_sup_f_95th": round(sup_95, 4),
        "null_sup_f_99th": round(sup_99, 4)
    }

    # ----------------------------------------------------------
    print(f"\n[6/8] Bootstrap CI for breakpoint year ({N_BOOTSTRAP} resamples)...")
    bootstrap_bp_years = []
    bootstrap_slope_diffs = []
    for b in range(N_BOOTSTRAP):
        if (b + 1) % 500 == 0:
            print(f"  Bootstrap {b+1}/{N_BOOTSTRAP}...")
        # Resample with replacement
        indices = [rng.randint(0, n - 1) for _ in range(n)]
        indices.sort()
        boot_years = [years[i] for i in indices]
        boot_anom = [anomalies[i] for i in indices]
        # Find best breakpoint in bootstrap sample
        _, _, _, ss1_boot = linear_regression(boot_years, boot_anom)
        df1_boot = n - 2
        max_f_boot = -1
        best_bp_boot = None
        best_slopes_boot = None
        for bp_year in candidate_years:
            bpi = find_breakpoint_index(boot_years, bp_year)
            if bpi < 3 or (n - bpi) < 3:
                continue
            s1b, _, s2b, _, ss2b, df2b = piecewise_regression(boot_years, boot_anom, bpi)
            fb = f_statistic(ss1_boot, df1_boot, ss2b, df2b)
            if fb > max_f_boot:
                max_f_boot = fb
                best_bp_boot = bp_year
                best_slopes_boot = (s1b, s2b)
        if best_bp_boot is not None:
            bootstrap_bp_years.append(best_bp_boot)
            bootstrap_slope_diffs.append((best_slopes_boot[1] - best_slopes_boot[0]) * 10)

    bootstrap_bp_years.sort()
    bootstrap_slope_diffs.sort()
    n_boot = len(bootstrap_bp_years)
    ci_low_idx = int(0.025 * n_boot)
    ci_high_idx = int(0.975 * n_boot)
    bp_ci_low = bootstrap_bp_years[ci_low_idx]
    bp_ci_high = bootstrap_bp_years[ci_high_idx]
    bp_median = bootstrap_bp_years[n_boot // 2]
    slope_diff_ci_low = bootstrap_slope_diffs[ci_low_idx]
    slope_diff_ci_high = bootstrap_slope_diffs[ci_high_idx]
    slope_diff_median = bootstrap_slope_diffs[n_boot // 2]

    print(f"  Breakpoint median: {bp_median}")
    print(f"  Breakpoint 95% CI: [{bp_ci_low}, {bp_ci_high}]")
    print(f"  Slope difference (post - pre, °C/decade) median: {slope_diff_median:.4f}")
    print(f"  Slope difference 95% CI: [{slope_diff_ci_low:.4f}, {slope_diff_ci_high:.4f}]")

    results["bootstrap"] = {
        "n_resamples": N_BOOTSTRAP,
        "breakpoint_median": bp_median,
        "breakpoint_95ci": [bp_ci_low, bp_ci_high],
        "slope_diff_per_decade_median": round(slope_diff_median, 4),
        "slope_diff_per_decade_95ci": [round(slope_diff_ci_low, 4), round(slope_diff_ci_high, 4)],
        "seed": SEED
    }

    # ----------------------------------------------------------
    print(f"\n[7/8] Sensitivity analysis...")
    # Test with different permutation counts and breakpoint ranges
    sensitivity = {}

    # 7a: Vary breakpoint search range
    for start, end in [(1900, 2010), (1930, 1990), (1940, 1980)]:
        cands = [y for y in years if start <= y <= end]
        bf = -1
        bbp = None
        for bp_year in cands:
            bpi = find_breakpoint_index(years, bp_year)
            if bpi < 3 or (n - bpi) < 3:
                continue
            _, _, _, _, ss_tw, df_tw = piecewise_regression(years, anomalies, bpi)
            fv = f_statistic(ss_one, df_one, ss_tw, df_tw)
            if fv > bf:
                bf = fv
                bbp = bp_year
        sensitivity[f"range_{start}_{end}"] = {"best_breakpoint": bbp, "f_statistic": round(bf, 4)}
        print(f"  Range [{start}-{end}]: best breakpoint={bbp}, F={bf:.4f}")

    # 7b: R-squared improvement
    ss_total_var = sum((y - mean(anomalies)) ** 2 for y in anomalies)
    r2_one = 1 - ss_one / ss_total_var
    r2_two = 1 - best_bp_info["ss_two"] / ss_total_var
    print(f"  R² (one segment): {r2_one:.4f}")
    print(f"  R² (two segments): {r2_two:.4f}")
    print(f"  R² improvement: {r2_two - r2_one:.4f}")

    # 7c: Residual autocorrelation (Durbin-Watson-like)
    _, _, resid_one, _ = linear_regression(years, anomalies)
    bp_idx_best = best_bp_info["bp_idx"]
    xs1 = years[:bp_idx_best + 1]
    ys1 = anomalies[:bp_idx_best + 1]
    xs2 = years[bp_idx_best:]
    ys2 = anomalies[bp_idx_best:]
    _, _, resid1, _ = linear_regression(xs1, ys1)
    _, _, resid2, _ = linear_regression(xs2, ys2)
    resid_two = resid1 + resid2[1:]  # avoid double counting breakpoint

    def durbin_watson(resids):
        num = sum((resids[i] - resids[i-1])**2 for i in range(1, len(resids)))
        den = sum(r**2 for r in resids)
        return num / den if den > 0 else 0

    dw_one = durbin_watson(resid_one)
    dw_two = durbin_watson(resid_two)
    print(f"  Durbin-Watson (one segment): {dw_one:.4f}")
    print(f"  Durbin-Watson (two segments): {dw_two:.4f}")

    # 7d: Effect size — ratio of post-break to pre-break warming rate
    rate_ratio = best_bp_info["slope2"] / best_bp_info["slope1"] if best_bp_info["slope1"] != 0 else float("inf")
    print(f"  Post/pre warming rate ratio: {rate_ratio:.2f}x")

    sensitivity["r_squared"] = {
        "one_segment": round(r2_one, 4),
        "two_segment": round(r2_two, 4),
        "improvement": round(r2_two - r2_one, 4)
    }
    sensitivity["durbin_watson"] = {
        "one_segment": round(dw_one, 4),
        "two_segment": round(dw_two, 4)
    }
    sensitivity["rate_ratio"] = round(rate_ratio, 2)

    results["sensitivity"] = sensitivity

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

    results["limitations"] = [
        "Annual averaging smooths sub-annual variability and volcanic cooling episodes",
        "Piecewise linear model assumes abrupt rate change; reality may be gradual acceleration",
        "Residual autocorrelation (DW < 2) means effective sample size is smaller than N, so CIs may be anti-conservative",
        "Single breakpoint model; multiple breakpoints are plausible but not tested here",
        "NOAA anomalies are relative to 1901-2000 base period; different baselines could shift perception",
        "Bootstrap CI for breakpoint year is discrete (integer years) — true change may be gradual"
    ]

    results["methodology"] = {
        "null_model": "Permutation F-test: shuffle anomalies, refit both models, compare F-statistics",
        "multiple_testing_correction": "Supremum F-test: for each permutation, take max F across ALL candidate breakpoints",
        "confidence_intervals": "Non-parametric bootstrap (resample with replacement), 2.5th-97.5th percentile",
        "n_permutations": N_PERMUTATIONS,
        "n_bootstrap": N_BOOTSTRAP,
        "seed": SEED
    }

    # Write results.json
    results_path = os.path.join(WORKSPACE, "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(WORKSPACE, "report.md")
    with open(report_path, "w") as f:
        f.write("# NOAA Temperature Trend Breakpoint Analysis — Results Report\n\n")
        f.write(f"**Data:** {n} annual observations ({years[0]}–{years[-1]})\n\n")
        f.write("## Single Linear Trend\n\n")
        f.write(f"- Slope: {slope_one*10:.4f} °C/decade\n")
        f.write(f"- R²: {r2_one:.4f}\n\n")
        f.write("## Best Breakpoint\n\n")
        f.write(f"- **Year: {best_bp_year}**\n")
        f.write(f"- Pre-break slope: {best_bp_info['slope1']*10:.4f} °C/decade\n")
        f.write(f"- Post-break slope: {best_bp_info['slope2']*10:.4f} °C/decade\n")
        f.write(f"- Post/pre rate ratio: {rate_ratio:.2f}x\n")
        f.write(f"- R² (two segments): {r2_two:.4f} (improvement: {r2_two - r2_one:.4f})\n\n")
        f.write("## Statistical Significance\n\n")
        f.write(f"- Permutation F-test (fixed breakpoint): p = {p_value:.6f} ({N_PERMUTATIONS} permutations)\n")
        f.write(f"- Supremum F-test (all breakpoints, multiple-testing corrected): p = {sup_p_value:.6f}\n")
        f.write(f"- Observed F = {observed_f:.4f}; null 95th = {perm_95:.4f}, 99th = {perm_99:.4f}\n\n")
        f.write("## Bootstrap Confidence Intervals\n\n")
        f.write(f"- Breakpoint year: median {bp_median}, 95% CI [{bp_ci_low}, {bp_ci_high}]\n")
        f.write(f"- Slope difference (post − pre, °C/decade): median {slope_diff_median:.4f}, 95% CI [{slope_diff_ci_low:.4f}, {slope_diff_ci_high:.4f}]\n\n")
        f.write("## Sensitivity\n\n")
        for key, val in sensitivity.items():
            if key.startswith("range_"):
                f.write(f"- Breakpoint range {key}: best={val['best_breakpoint']}, F={val['f_statistic']}\n")
        f.write(f"- Durbin-Watson: one-segment={dw_one:.4f}, two-segment={dw_two:.4f}\n\n")
        f.write("## Limitations\n\n")
        for lim in results["limitations"]:
            f.write(f"- {lim}\n")
        f.write("\n---\n*Generated by Claw4S NOAA Temperature Breakpoint Skill v1.0.0*\n")
    print(f"  Wrote {report_path}")

    print("\nANALYSIS COMPLETE")
    return results


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

def run_verify():
    """Machine-checkable assertions on results.json."""
    results_path = os.path.join(WORKSPACE, "results.json")
    report_path = os.path.join(WORKSPACE, "report.md")

    print("VERIFICATION MODE")
    print("=" * 50)

    assert os.path.exists(results_path), f"FAIL: {results_path} does not exist"
    assert os.path.exists(report_path), f"FAIL: {report_path} does not exist"

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

    checks = 0
    failures = 0

    def check(condition, msg):
        nonlocal checks, failures
        checks += 1
        if condition:
            print(f"  PASS [{checks}]: {msg}")
        else:
            print(f"  FAIL [{checks}]: {msg}")
            failures += 1

    # 1. Data loaded
    check(r["n_observations"] >= 100, f"At least 100 observations (got {r['n_observations']})")

    # 2. Year range starts at 1880
    check(r["year_range"][0] == 1880, f"Data starts at 1880 (got {r['year_range'][0]})")

    # 3. Single segment slope is positive (warming)
    check(r["single_segment"]["slope_per_year"] > 0, "Single segment slope is positive (warming)")

    # 4. Best breakpoint is in search range
    bp = r["best_breakpoint"]["year"]
    check(BREAKPOINT_RANGE[0] <= bp <= BREAKPOINT_RANGE[1],
          f"Best breakpoint {bp} is in range [{BREAKPOINT_RANGE[0]}, {BREAKPOINT_RANGE[1]}]")

    # 5. F-statistic is positive
    check(r["best_breakpoint"]["f_statistic"] > 0, "F-statistic is positive")

    # 6. Permutation p-value exists and is in [0, 1]
    pv = r["permutation_test"]["p_value"]
    check(0 <= pv <= 1, f"Permutation p-value {pv} is in [0, 1]")

    # 7. Supremum test was run
    check(r["supremum_f_test"]["n_permutations"] == N_PERMUTATIONS,
          f"Supremum test used {N_PERMUTATIONS} permutations")

    # 8. Bootstrap CI exists and brackets the best breakpoint
    ci = r["bootstrap"]["breakpoint_95ci"]
    check(ci[0] <= bp <= ci[1],
          f"Bootstrap 95% CI [{ci[0]}, {ci[1]}] contains best breakpoint {bp}")

    # 9. Post-break slope > pre-break slope (accelerating warming)
    check(r["best_breakpoint"]["post_break_slope_per_decade"] > r["best_breakpoint"]["pre_break_slope_per_decade"],
          "Post-break warming rate exceeds pre-break rate")

    # 10. R² improvement from two-segment model
    check(r["sensitivity"]["r_squared"]["improvement"] > 0,
          "Two-segment model has higher R² than one-segment")

    # 11. Results contain limitations
    check(len(r["limitations"]) >= 4, f"At least 4 limitations listed (got {len(r['limitations'])})")

    # 12. Report file is non-empty
    with open(report_path) as f:
        report_content = f.read()
    check(len(report_content) > 500, f"Report is substantial ({len(report_content)} chars)")

    print(f"\n{'=' * 50}")
    print(f"VERIFICATION: {checks - failures}/{checks} checks passed")
    if failures == 0:
        print("ALL CHECKS PASSED")
    else:
        print(f"WARNING: {failures} check(s) failed")
        sys.exit(1)


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

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="NOAA Temperature Breakpoint Analysis")
    parser.add_argument("--verify", action="store_true", help="Run verification checks on results.json")
    args = parser.parse_args()

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

**Expected output:** No stdout. File `analysis.py` created in workspace.

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_noaa-temperature-trend-breakpoint && python3 analysis.py
```

**Expected output:** Sectioned output `[1/8]` through `[8/8]`, ending with `ANALYSIS COMPLETE`. Creates `results.json` and `report.md` in workspace. Runtime: approximately 10-15 minutes due to supremum F-test (5000 permutations × 81 breakpoints).

**Key expected values (approximate):**
- N observations: ~145 (1880-2024)
- Best breakpoint: approximately 1960-1980
- Permutation p-value: < 0.05 (significant breakpoint)
- Post-break warming rate: substantially higher than pre-break rate

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_noaa-temperature-trend-breakpoint && python3 analysis.py --verify
```

**Expected output:** 12 PASS checks, ending with `ALL CHECKS PASSED`.

## Success Criteria

1. `results.json` exists and contains all required fields
2. `report.md` exists and is >500 characters
3. All 12 verification checks pass
4. Permutation p-value is computed from 5,000 shuffles
5. Supremum F-test corrects for multiple testing across breakpoint candidates
6. Bootstrap 95% CI is computed from 2,000 resamples

## Failure Conditions

1. Data download fails after 3 retries → script exits with error
2. SHA256 mismatch on cached data → re-downloads
3. Fewer than 100 data points loaded → assertion failure
4. Any verification check fails → exit code 1