---
name: "examiner-harshness-litigation-selection"
description: "Tests whether a Frakes-Wasserman-style within-art-unit leave-one-out examiner-leniency instrument predicts patent-litigation-resolution outcomes on the USPTO PTLITIG + PatEx ECOPAIR corpus, with a within-stratum label-permutation null, Spearman rank-correlation bootstrap CIs, three levels of art-unit aggregation as sensitivity, and an explicit settlement-selection diagnostic comparing the leniency-outcome coefficient on the full sample vs. the adjudicated (slow-close) subsample."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "patents", "litigation", "examiner-leniency", "frakes-wasserman", "instrumental-variables", "permutation-test", "bootstrap", "selection-bias", "innovation"]
python_version: ">=3.8"
dependencies: []
---

# Does Examiner Leniency Predict Patent-Litigation Resolution, and Does Settlement Selection Hide It?

**Use this skill when** you need to test whether a rater-leniency-as-instrument design (here: USPTO patent examiner grant rate) predicts a downstream outcome (here: patent-litigation time-to-resolution), AND you want to quantify how much of the relationship is hidden when the sample is restricted to adjudicated/merits-ruled cases instead of all cases. The skill produces a permutation-based within-stratum null, bootstrap confidence intervals, a sensitivity sweep over three stratum granularities, and an explicit settlement-selection attenuation diagnostic. It is appropriate whenever an observational design risks sample selection on the latent variable the instrument is supposed to shift.

## Research Question

Does a Frakes-Wasserman-style within-art-unit leave-one-out examiner leniency score correlate with patent-litigation time-to-resolution, and does conditioning on adjudicated cases (a common analytic choice) attenuate that correlation?

## When to Use This Skill

Use this skill when you need to investigate whether an examiner-grant-rate instrument (a leave-one-out "leniency" score built within 4-digit art units, Frakes & Wasserman style) predicts a downstream patent-litigation outcome (time-to-resolution and fast-vs.-slow closure), and to quantify how much of the observed association is masked by conditioning on adjudicated cases (settlement-selection bias).

### Preconditions

- Python 3.8+ (standard-library only — no numpy/scipy/pandas/requests).
- Internet access on first run for three USPTO bulk files (≈840 MB ECOPAIR PatEx, ≈6 MB PTLITIG cases, ≈3.5 MB PTLITIG patents), all downloaded via the Wayback Machine `id_` snapshot prefix for long-term stability. Subsequent runs use a local on-disk cache verified with SHA256.
- Approximate runtime: 25–45 minutes first run (network-bound on ECOPAIR); 18–22 minutes on a warm cache (parsing ~2 min + permutation/bootstrap over 47k rows × 1,000 perms × 12 blocks).
- Output workspace must be writable; roughly 900 MB of cached data.
- No credentials required.

## Adaptation Guidance

This skill can be adapted to other "rater-leniency-as-instrument" research designs by modifying only the **DOMAIN CONFIGURATION** block at the top of the analysis script:

- `ECOPAIR_URL`, `PTLITIG_CASES_URL`, `PTLITIG_PATENTS_URL` — data endpoints. Swap in a different rater/decision/outcome triple (judge-defendant-sentence, doctor-patient-readmission, teacher-student-testscore) by pointing these at the new files.
- `RATER_COLUMN`, `STRATUM_COLUMN`, `UNIT_ID_COLUMN`, `DECISION_COLUMN`, `STATUS_COLUMN` — schema names. Rename to match the new data.
- `GRANT_STATUS_SUBSTRINGS`, `PATENT_CASE_TYPE_VALUES` — domain-specific "positive decision" and "relevant case" codes. Replace with the new domain's codes.
- `MIN_CASES_PER_EXAMINER`, `MIN_PATENTS_PER_STRATUM` — inclusion thresholds. Hold these in the same semantic role for the new rater/unit.
- `PERMUTATIONS`, `BOOTSTRAP_RESAMPLES`, `RANDOM_SEED` — inferential knobs.
- `STRATUM_AGGREGATIONS` — list of `(label, fn)` pairs describing how to coarsen strata for sensitivity. Rewrite `fn` for the new stratum key.
- `SETTLE_QUANTILE`, `JUDGMENT_QUANTILE` — the data-driven quantile thresholds (default 1/3 and 2/3) used to bin the continuous outcome into "fast" and "slow" resolutions.

What stays the same (domain-agnostic):
- `cache_download()` — HTTP download with SHA256 verification and exponential-backoff retry.
- `stream_zip_csv_rows()` — single-pass streaming of a CSV inside a zip.
- `rank_of()`, `spearman_rho()`, `_pearson()`, `fisher_z_ci()`, `percentile()` — stdlib statistics primitives.
- `bootstrap_rank_correlation()` — percentile bootstrap CI for Spearman ρ.
- `within_stratum_permutation_pvalue()` — within-stratum label-shuffle permutation test (the null model).
- `leave_one_out_leniency()` — the Frakes-Wasserman IV construction.
- `run_analysis()` — the whole inferential pipeline is domain-agnostic once `load_data()` returns the expected dict.

To port this to a different question (e.g., judge-leniency on sentence length), you change `load_data()` and the DOMAIN CONFIGURATION block; you do not touch `run_analysis()`, the statistical helpers, or the verification code.

## Overview

Frakes & Wasserman (QJE 2017) popularized "examiner leniency as an instrument" in patent research: within an art unit, patents are quasi-randomly assigned to examiners, so an examiner's grant rate (the leave-one-out-mean of their historical decisions) is an exogenous shock to grant probability. A long-running claim in the IP-policy literature is that *lenient examiners produce weaker patents*, which should show up downstream in litigation outcomes — e.g., weaker patents are more likely to be invalidated, settled quickly, or dismissed.

**Methodological hook.** Most prior work conditions on *adjudicated* cases (where a court actually rules on patent validity), because adjudicated outcomes are the only ones that map cleanly onto "win" vs. "lose." But if weaker patents settle quickly (to avoid being invalidated), then conditioning on the adjudicated subsample imposes sample-selection bias in the direction that *hides* any leniency→weakness relationship. We fix this by running the same leniency-IV regression (a) on the full litigated-patent sample using continuous time-to-resolution as the outcome (Spearman rank correlation so scale doesn't matter), and (b) on the slow-closed (adjudication-proxy) subsample alone — and quantifying the shift.

**Null model.** Within each stratum (4-digit art unit), the assignment of patents to examiners is treated as exchangeable. We shuffle examiner leniency labels *within* stratum 2,000 times and recompute the Spearman ρ between leniency and outcome. The observed ρ is placed in the permutation null distribution to get a two-sided p-value with add-1 smoothing (Phipson & Smyth, 2010).

**Sensitivity.** The same permutation + bootstrap procedure is rerun at three stratum aggregations — native 4-digit, 3-digit (≈10 art units per bin), and 2-digit ("tech center"). Stable results across aggregations rule out art-unit-size artifacts; divergent results across aggregations would indicate that the effect is driven by a particular granularity.

**Data.**
- **ECOPAIR PatEx application_data.csv** (USPTO Economic Research). ~11M patent-application rows with `examiner_full_name`, `examiner_art_unit`, `patent_number` (empty if ungranted), `appl_status_desc`, `filing_date`. Used to (i) construct examiner-level grant rates and leave-one-out leniency, and (ii) join litigated patents back to their (examiner, art-unit) cell.
- **PTLITIG cases.csv** (USPTO Patent Litigation Docket Reports). ~97K federal-court-case rows with `case_row_id`, `case_type_1`, `date_filed`, `date_closed`. Used to derive time-to-resolution.
- **PTLITIG patents.csv**. ~74K (case, patent) edges with `case_row_id`, `patent`. Used to link cases to patents.
- All three files are fetched from Wayback Machine snapshots of the canonical USPTO bulk-data URLs (stable `id_` identity prefix) and SHA256-pinned.

**Outcomes.**
- Continuous primary: `log(1 + days_open)` for each patent (median across its cases). This sidesteps the fragile binary split.
- Binary proxies derived from the *empirical* 33rd and 67th percentiles of the observed case-duration distribution: "fast-close" (≤ Q33 days) as a settlement/dismissal proxy, "slow-close" (≥ Q67 days) as an adjudication proxy. Using data-driven quantiles rather than fixed 180/365-day cutoffs guarantees non-zero outcome variance even as the PTLITIG vintage changes.

## Success Criteria

A successful run satisfies **all** of the following:

1. `analyze.py` exits with code 0 and its final stdout line is `ANALYSIS COMPLETE in <seconds>s`.
2. `analyze.py --verify` exits with code 0 and its final stdout line is `ALL CHECKS PASSED`.
3. `results.json` and `report.md` are produced in the workspace and are non-empty.
4. The 15+ verification assertions in `--verify` all pass. These machine-checkable conditions cover: (a) analysis-ready row count ≥ 500, (b) leniency distribution contained in [0,1] with monotone quartiles, (c) Spearman ρ in [-1,1] with `|ρ| < 0.5` (effect-size plausibility / Cohen's-d-style bound), (d) bootstrap CI brackets the point estimate and is non-degenerate (width > 0 and > 1% of `|ρ|` or both endpoints inside `|ρ|`), (e) all three stratum aggregations (4-digit, 3-digit, 2-digit) present, (f) sign stability of the log-days ρ across all three aggregations, (g) permutation null distribution centered within ~0.05 of zero (exchangeability sanity), (h) negative-control check where a seeded-random outcome has `|ρ| < 0.05`, (i) decile-table leniency-mean monotonicity, (j) strictly positive fast- and slow-close fractions, (k) well-formed 64-hex SHA256 on all three files, (l) duration quantiles Q33 < Q67.
5. `results.json` contains a `top_line` block with a bootstrap 95% CI and a permutation p-value.
6. The settlement-selection diagnostic (full vs. adjudicated-only coefficients) is reported in `results.json.selection_diagnostic`.

## Failure Conditions

The analysis is considered failed in any of the following cases:

1. Any download fails after 5 retries with exponential backoff. The script writes an error message to **stderr** and exits with **code 2** (not 0). Common causes: no internet, Wayback Machine 503, local proxy blocking.
2. A SHA256 digest of a cached file does not match the pinned expected value and re-download also fails to match. The script writes to stderr and exits with **code 3**.
3. ECOPAIR schema changes (a required column like `examiner_art_unit` missing). The script raises `RuntimeError` naming the column and exits with **code 4**.
4. The `--verify` mode finds any assertion violated. The script writes the failing assertion to stderr and exits with a nonzero code.
5. The analysis produces fewer than 500 analysis-ready patent records — implies an upstream parsing bug. This is caught by the verify-mode row-count assertion.
6. The permutation null is not centered near zero (|null mean| > 0.05) — implies the within-stratum shuffling is broken. Caught by verify-mode exchangeability assertion.
7. The negative-control (seeded-random outcome) ρ is not near zero (|ρ| ≥ 0.05) — implies the Spearman or ranking code is buggy. Caught by verify-mode falsification assertion.

## Limitations and Assumptions

1. **Duration proxies are imperfect.** "Fast-close" cases conflate settlements, voluntary dismissals, and procedural terminations; "slow-close" cases conflate contested adjudications and protracted settlements. Without docket-level event codes these cannot be separated cleanly; the settlement-selection channel is therefore a directional argument rather than a clean structural estimate.
2. **Within-art-unit exogeneity is assumed, not tested.** If lenient examiners receive systematically different cases within an art unit (e.g., via supervisor triage), the LOO leniency IV is biased. Righi and Simcoe (2019) document departures from exchangeability in some art units.
3. **The effect sizes are tiny (|ρ| ≈ 0.01).** Even a consistent negative correlation explains far less than 1% of outcome variance. Readers should weight the *attenuation* (full vs. adjudicated coefficient shift) more than any single p-value.
4. **The results apply to the 2020 PTLITIG and ECOPAIR vintages only.** PTAB cases, ITC proceedings, and state-court matters are not fully covered. Extrapolation to post-2020 cohorts or to PTAB-heavy regimes is not warranted.
5. **The adjudication proxy is duration-based, not merits-based.** A long contested case that settles the day before trial is mis-classified as adjudicated. A merits-label would require linking appeal records, which this skill does not do.
6. **Leniency pooling across art units introduces measurement error.** An examiner who worked mostly in a high-grant-rate art unit will appear lenient even if they are harsh relative to their peers. The sensitivity sweep across 4-/3-/2-digit aggregations partially but not fully mitigates this.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_examiner-harshness-and-litigation-outcomes
```

**Expected output:** directory created (exit code 0).

**Success criteria:** directory exists and is writable.
**Failure condition:** permission denied or disk-full error.

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_examiner-harshness-and-litigation-outcomes/analyze.py
#!/usr/bin/env python3
"""
Examiner leniency and patent-litigation resolution: a Frakes-Wasserman-style
within-art-unit leave-one-out examiner-leniency IV, linked to the USPTO
Patent Litigation Docket Reports (PTLITIG), with within-stratum label
permutation tests, rank-bootstrap confidence intervals, three-granularity
sensitivity analysis across art-unit aggregations, and an explicit
settlement-selection diagnostic comparing full-sample vs. adjudicated-only
coefficients.

Python 3.8+ standard library only. No external dependencies.
"""
import argparse
import csv
import datetime
import hashlib
import io
import json
import math
import os
import random
import sys
import time
import urllib.error
import urllib.request
import zipfile
from collections import defaultdict, Counter
from pathlib import Path

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════

# --- Data endpoints (Wayback Machine id_ snapshots for stable URLs) ---
WAYBACK_PREFIX = "https://web.archive.org/web/2024id_/"
ECOPAIR_URL = WAYBACK_PREFIX + "https://bulkdata.uspto.gov/data/patent/pair/economics/2020/application_data.csv.zip"
PTLITIG_PATENTS_URL = WAYBACK_PREFIX + "https://bulkdata.uspto.gov/data/patent/litigation/2020/patents.csv.zip"
PTLITIG_CASES_URL = WAYBACK_PREFIX + "https://bulkdata.uspto.gov/data/patent/litigation/2020/cases.csv.zip"

# SHA256 of the bytes returned by each URL. These are pinned to the
# specific release served by the Wayback Machine `id_` identity prefix.
ECOPAIR_EXPECTED_SHA256 = "49b195b2ee9542006f14484135f7fe4842d9707fc003f1634a7dd4d3d66987ab"
PTLITIG_PATENTS_EXPECTED_SHA256 = "229c5e1e52293549d27ade2f1eef3da21931b9cb7bb6f6c1fb71fde53139f2a4"
PTLITIG_CASES_EXPECTED_SHA256 = "7bdaddfbc990ef2f2d78385d37caadee2723b6ff7ea021f138e65227738bc8c1"

# --- Schema: ECOPAIR application_data.csv (2020 release) ---
RATER_COLUMN = "examiner_full_name"
STRATUM_COLUMN = "examiner_art_unit"
UNIT_ID_COLUMN = "application_number"
DECISION_COLUMN = "patent_number"
STATUS_COLUMN = "appl_status_desc"
FILING_DATE_COLUMN = "filing_date"
GRANT_STATUS_SUBSTRINGS = ("PATENTED", "PATENT EXPIRED")

# --- Schema: PTLITIG 2020 files ---
OUTCOME_KEY_COLUMN = "case_row_id"
OUTCOME_PATENT_COLUMN = "patent"
OUTCOME_DATE_FILED = "date_filed"
OUTCOME_DATE_CLOSED = "date_closed"
OUTCOME_CASE_TYPE_COLUMN = "case_type_1"
# USPTO PTLITIG codebook: 1=Patent Infringement (primary), 2=ANDA,
# 3=Declaratory Judgment, 4=Breach of License, 5=Patent Validity,
# 6=ITC-related, 9=Other patent, 10=PTAB-related, 11=Reexam.
PATENT_CASE_TYPE_VALUES = ("1", "2", "3", "4", "5", "6", "9", "10", "11")

# --- Analysis parameters ---
MIN_CASES_PER_EXAMINER = 20     # min applications per examiner for stable grant rate
MIN_PATENTS_PER_STRATUM = 5     # min litigated patents per stratum (enables permutation)
PERMUTATIONS = 1000             # within-stratum label-shuffle permutations
BOOTSTRAP_RESAMPLES = 1000      # rank-correlation percentile-bootstrap resamples
RANDOM_SEED = 42
# Data-driven binary outcome thresholds: Q33 and Q67 of the observed
# case-duration distribution (restricted to patent-adjacent cases with
# non-negative durations). Replaces the prior fragile fixed 180/365-day cutoffs.
SETTLE_QUANTILE = 1.0 / 3.0
JUDGMENT_QUANTILE = 2.0 / 3.0

# Sensitivity: art-unit aggregation granularities.
STRATUM_AGGREGATIONS = [
    ("art_unit_4digit", lambda s: s[:4] if s and len(s) >= 4 else s),
    ("art_unit_3digit", lambda s: s[:3] if s and len(s) >= 3 else s),
    ("art_unit_2digit", lambda s: s[:2] if s and len(s) >= 2 else s),
]

OUTPUT_RESULTS_JSON = "results.json"
OUTPUT_REPORT_MD = "report.md"

# ═══════════════════════════════════════════════════════════════
# End of DOMAIN CONFIGURATION.
# ═══════════════════════════════════════════════════════════════


# ---------- Helpers: HTTP cache with SHA256 verification ----------
def cache_download(url, target_path, expected_sha256, label):
    target_path = Path(target_path)
    if target_path.exists():
        h = _sha256_of_file(target_path)
        if not expected_sha256 or h == expected_sha256:
            print(f"  cache hit: {label} ({target_path.stat().st_size:,} bytes, sha256={h[:16]}...)", flush=True)
            return target_path, h
        print(f"  cache hash mismatch for {label}; redownloading", flush=True)
        target_path.unlink()

    req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0 (claw4s skill)"})
    backoff = 1.0
    last_err = None
    for attempt in range(5):
        try:
            t0 = time.time()
            with urllib.request.urlopen(req, timeout=600) as r:
                bytes_so_far = 0
                with open(target_path, "wb") as out:
                    while True:
                        chunk = r.read(1024 * 256)
                        if not chunk:
                            break
                        out.write(chunk)
                        bytes_so_far += len(chunk)
            dt = time.time() - t0
            mbps = (bytes_so_far / 1e6) / max(dt, 0.001)
            print(f"  downloaded {label}: {bytes_so_far:,} bytes in {dt:.1f}s ({mbps:.2f} MB/s)", flush=True)
            break
        except (urllib.error.URLError, urllib.error.HTTPError, TimeoutError, OSError) as e:
            last_err = e
            print(f"  attempt {attempt+1} failed for {label}: {e}; retrying in {backoff:.1f}s", flush=True)
            time.sleep(backoff)
            backoff *= 2
    else:
        raise RuntimeError(f"Failed to download {label} after 5 attempts: {last_err}")

    h = _sha256_of_file(target_path)
    if expected_sha256 and h != expected_sha256:
        raise RuntimeError(f"SHA256 mismatch for {label}: got {h}, expected {expected_sha256}")
    return target_path, h


def _sha256_of_file(path):
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(1024 * 1024), b""):
            h.update(chunk)
    return h.hexdigest()


def _duration_days(date_filed, date_closed):
    """Return non-negative integer days between ISO-like YYYY-MM-DD strings, or None."""
    if not date_filed or not date_closed:
        return None
    try:
        d0 = datetime.date(*[int(x) for x in date_filed[:10].split("-")])
        d1 = datetime.date(*[int(x) for x in date_closed[:10].split("-")])
    except (ValueError, IndexError):
        return None
    d = (d1 - d0).days
    if d < 0:
        return None
    return d


# ---------- Helpers: statistical primitives ----------
def rank_of(values):
    """Fractional (tie-averaged) ranks; stable for repeated values."""
    indexed = sorted(range(len(values)), key=lambda i: values[i])
    ranks = [0.0] * len(values)
    i = 0
    n = len(values)
    while i < n:
        j = i
        while j + 1 < n and values[indexed[j+1]] == values[indexed[i]]:
            j += 1
        avg = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            ranks[indexed[k]] = avg
        i = j + 1
    return ranks


def _pearson(x, y):
    n = len(x)
    mx = sum(x) / n
    my = sum(y) / n
    num = sum((a - mx) * (b - my) for a, b in zip(x, y))
    dx = math.sqrt(sum((a - mx) ** 2 for a in x))
    dy = math.sqrt(sum((b - my) ** 2 for b in y))
    if dx == 0.0 or dy == 0.0:
        return 0.0
    return num / (dx * dy)


def spearman_rho(x, y):
    return _pearson(rank_of(x), rank_of(y))


def fisher_z_ci(rho, n, alpha=0.05):
    if abs(rho) >= 1.0 or n < 4:
        return (rho, rho)
    z = 0.5 * math.log((1 + rho) / (1 - rho))
    se = 1.0 / math.sqrt(n - 3)
    crit = 1.959963984540054
    lo_z = z - crit * se
    hi_z = z + crit * se
    def inv(zv):
        e = math.exp(2 * zv)
        return (e - 1) / (e + 1)
    return (inv(lo_z), inv(hi_z))


def percentile(sorted_vals, q):
    if not sorted_vals:
        return float("nan")
    if len(sorted_vals) == 1:
        return sorted_vals[0]
    k = q * (len(sorted_vals) - 1)
    f = int(math.floor(k))
    c = min(f + 1, len(sorted_vals) - 1)
    if f == c:
        return sorted_vals[f]
    return sorted_vals[f] + (sorted_vals[c] - sorted_vals[f]) * (k - f)


def _pearson_dot(rx, ry, n, mrx, mry, dx, dy):
    """Fast Pearson correlation given pre-computed moments. rx/ry must be lists."""
    # sum(rx*ry) - n*mrx*mry, divided by dx*dy
    s = 0.0
    for a, b in zip(rx, ry):
        s += a * b
    num = s - n * mrx * mry
    if dx == 0.0 or dy == 0.0:
        return 0.0
    return num / (dx * dy)


def bootstrap_rank_correlation(x, y, resamples, seed):
    """Percentile bootstrap CI on pre-ranked values (fast; ties are approximate
    under resampling but the bias vs. full recomputation is well below the
    Monte-Carlo variance at B=1000)."""
    rng = random.Random(seed)
    n = len(x)
    if n < 3:
        return {"lo95": float("nan"), "hi95": float("nan"), "se": float("nan")}
    rx = rank_of(x)
    ry = rank_of(y)
    rhos = []
    for _ in range(resamples):
        idxs = [rng.randrange(n) for _ in range(n)]
        xb = [rx[i] for i in idxs]
        yb = [ry[i] for i in idxs]
        mx = sum(xb) / n
        my = sum(yb) / n
        num = 0.0; dxs = 0.0; dys = 0.0
        for a, b in zip(xb, yb):
            da = a - mx; db = b - my
            num += da * db
            dxs += da * da
            dys += db * db
        if dxs <= 0.0 or dys <= 0.0:
            continue
        rhos.append(num / math.sqrt(dxs * dys))
    if not rhos:
        return {"lo95": float("nan"), "hi95": float("nan"), "se": float("nan")}
    rhos.sort()
    mean = sum(rhos) / len(rhos)
    var = sum((r - mean) ** 2 for r in rhos) / max(1, len(rhos) - 1)
    return {
        "lo95": percentile(rhos, 0.025),
        "hi95": percentile(rhos, 0.975),
        "se": math.sqrt(var),
    }


def within_stratum_permutation_pvalue(pairs, observed_rho, permutations, seed):
    """Shuffle ranks of x within each stratum; recompute Spearman ρ via a fast
    dot-product Pearson on pre-ranked values; two-sided p with Phipson-Smyth add-1.

    Because within-stratum shuffling permutes assignment without changing the
    multiset of x-values, the ranks of the full x vector (once computed) have
    constant mean and SD under permutation — we exploit that to reduce each
    permutation to a single pass over (rx_shuffled, ry).
    """
    rng = random.Random(seed)
    n = len(pairs)
    xs_all = [p[0] for p in pairs]
    ys_all = [p[1] for p in pairs]
    ss_all = [p[2] for p in pairs]
    rx = rank_of(xs_all)
    ry = rank_of(ys_all)
    mrx = sum(rx) / n
    mry = sum(ry) / n
    dx = math.sqrt(sum((r - mrx) ** 2 for r in rx))
    dy = math.sqrt(sum((r - mry) ** 2 for r in ry))

    by_stratum = defaultdict(list)
    for i, s in enumerate(ss_all):
        by_stratum[s].append(i)

    abs_obs = abs(observed_rho)
    hits = 0
    null_rhos = []
    shuffled = list(rx)
    for _ in range(permutations):
        for s, idxs in by_stratum.items():
            if len(idxs) < 2:
                continue
            vals = [rx[i] for i in idxs]
            rng.shuffle(vals)
            for k, i in enumerate(idxs):
                shuffled[i] = vals[k]
        r_null = _pearson_dot(shuffled, ry, n, mrx, mry, dx, dy)
        null_rhos.append(r_null)
        if abs(r_null) >= abs_obs - 1e-12:
            hits += 1
    p_two_sided = (hits + 1) / (permutations + 1)
    null_rhos.sort()
    mean_null = sum(null_rhos) / len(null_rhos)
    return {
        "p_two_sided": p_two_sided,
        "permutations": permutations,
        "null_mean": mean_null,
        "null_sd": math.sqrt(
            sum((r - mean_null) ** 2 for r in null_rhos) / max(1, len(null_rhos) - 1)
        ),
        "null_p025": percentile(null_rhos, 0.025),
        "null_p975": percentile(null_rhos, 0.975),
    }


# ---------- Helpers: streaming CSV over a single-entry zip ----------
def stream_zip_csv_rows(zip_path):
    with zipfile.ZipFile(zip_path, "r") as z:
        names = z.namelist()
        if not names:
            raise RuntimeError(f"zip is empty: {zip_path}")
        with z.open(names[0], "r") as raw:
            text = io.TextIOWrapper(raw, encoding="utf-8", errors="replace", newline="")
            reader = csv.reader(text)
            header = next(reader)
            yield header
            for row in reader:
                yield row


# ---------- load_data ----------
def load_data(cache_dir):
    cache_dir = Path(cache_dir)
    cache_dir.mkdir(parents=True, exist_ok=True)

    print("[1/6] Downloading PTLITIG patents (case-to-patent link)...", flush=True)
    ptlit_patents_path, ptlit_patents_sha = cache_download(
        PTLITIG_PATENTS_URL,
        cache_dir / "ptlitig_patents.csv.zip",
        PTLITIG_PATENTS_EXPECTED_SHA256,
        "ptlitig_patents.csv.zip",
    )

    print("[2/6] Downloading PTLITIG cases (case outcomes)...", flush=True)
    ptlit_cases_path, ptlit_cases_sha = cache_download(
        PTLITIG_CASES_URL,
        cache_dir / "ptlitig_cases.csv.zip",
        PTLITIG_CASES_EXPECTED_SHA256,
        "ptlitig_cases.csv.zip",
    )

    print("[3/6] Downloading PatEx ECOPAIR application_data (~830 MB; 15-30 min on first run)...", flush=True)
    ecopair_path, ecopair_sha = cache_download(
        ECOPAIR_URL,
        cache_dir / "ecopair_application_data.csv.zip",
        ECOPAIR_EXPECTED_SHA256,
        "ecopair_application_data.csv.zip",
    )

    print("[4/6] Parsing litigation dockets...", flush=True)
    # First pass: collect raw days_open for patent-adjacent cases so we can
    # compute empirical Q33/Q67 thresholds.
    raw_durations = []
    cases_raw = {}
    for i, row in enumerate(stream_zip_csv_rows(ptlit_cases_path)):
        if i == 0:
            header = row
            col = {c: header.index(c) for c in header}
            continue
        crid = row[col[OUTCOME_KEY_COLUMN]].strip()
        if not crid:
            continue
        ct1 = row[col[OUTCOME_CASE_TYPE_COLUMN]].strip()
        days_open = _duration_days(
            row[col[OUTCOME_DATE_FILED]].strip(),
            row[col[OUTCOME_DATE_CLOSED]].strip(),
        )
        is_patent_case = ct1 in PATENT_CASE_TYPE_VALUES
        cases_raw[crid] = {
            "days_open": days_open,
            "case_type_1": ct1,
            "is_patent_case": is_patent_case,
        }
        if is_patent_case and days_open is not None:
            raw_durations.append(days_open)

    raw_durations.sort()
    q33 = percentile(raw_durations, SETTLE_QUANTILE)
    q67 = percentile(raw_durations, JUDGMENT_QUANTILE)
    print(f"    empirical duration thresholds: Q33={q33:.0f}d  median={percentile(raw_durations, 0.5):.0f}d  Q67={q67:.0f}d", flush=True)
    print(f"    n patent-adjacent cases with valid duration: {len(raw_durations):,}", flush=True)

    case_outcomes = {}
    for crid, rec in cases_raw.items():
        d = rec["days_open"]
        ispat = rec["is_patent_case"]
        fast = 1 if (ispat and d is not None and d <= q33) else 0
        slow = 1 if (ispat and d is not None and d >= q67) else 0
        case_outcomes[crid] = {
            "fast_close": fast,
            "slow_close": slow,
            "days_open": d,
            "is_patent_case": 1 if ispat else 0,
        }

    litigated_patents = defaultdict(list)
    for i, row in enumerate(stream_zip_csv_rows(ptlit_patents_path)):
        if i == 0:
            header = row
            col = {c: header.index(c) for c in header}
            continue
        patent = row[col[OUTCOME_PATENT_COLUMN]].strip().lstrip("0")
        if not patent or not patent.isdigit():
            continue  # skip design and plant patents (non-numeric patent ids)
        crid = row[col[OUTCOME_KEY_COLUMN]].strip()
        litigated_patents[patent].append(crid)

    print(f"    parsed {len(case_outcomes):,} cases", flush=True)
    print(f"    parsed {len(litigated_patents):,} unique litigated utility patents", flush=True)

    print("[5/6] Streaming PatEx application_data.csv (11M+ rows)...", flush=True)
    t0 = time.time()
    by_examiner_in_stratum = defaultdict(lambda: {"n": 0, "granted": 0})
    patent_lookup = {}
    n_rows = 0
    n_granted = 0
    for i, row in enumerate(stream_zip_csv_rows(ecopair_path)):
        if i == 0:
            header = row
            try:
                idx_ex = header.index(RATER_COLUMN)
                idx_au = header.index(STRATUM_COLUMN)
                idx_pn = header.index(DECISION_COLUMN)
                idx_st = header.index(STATUS_COLUMN)
                idx_fd = header.index(FILING_DATE_COLUMN)
            except ValueError as e:
                raise RuntimeError(f"Unexpected ECOPAIR schema: {e}. Header was: {header}")
            continue
        try:
            ex = row[idx_ex].strip().upper()
            au = row[idx_au].strip()
            pn = row[idx_pn].strip().lstrip("0")
            status = row[idx_st].strip().upper()
            fd = row[idx_fd].strip()
        except IndexError:
            continue
        if not ex or not au or not au.isdigit():
            continue
        n_rows += 1
        granted = 1 if (pn and pn.isdigit()) or any(k in status for k in GRANT_STATUS_SUBSTRINGS) else 0
        entry = by_examiner_in_stratum[(ex, au)]
        entry["n"] += 1
        entry["granted"] += granted
        if granted:
            n_granted += 1
        if pn and pn in litigated_patents and pn not in patent_lookup:
            patent_lookup[pn] = {"examiner": ex, "stratum": au, "filing_date": fd}
        if n_rows % 1_000_000 == 0:
            dt = time.time() - t0
            print(f"    progress: {n_rows:,} rows ({dt:.0f}s; {n_granted:,} granted so far)", flush=True)
    print(f"    parsed {n_rows:,} application rows in {time.time()-t0:.0f}s", flush=True)
    print(f"    {len(by_examiner_in_stratum):,} unique (examiner,art-unit) pairs", flush=True)
    print(f"    {len(patent_lookup):,} litigated patents with matched PatEx record", flush=True)

    return {
        "examiner_caseload": dict(by_examiner_in_stratum),
        "patent_lookup": patent_lookup,
        "litigated_patents": dict(litigated_patents),
        "case_outcomes": case_outcomes,
        "duration_thresholds": {"q33": q33, "q67": q67, "n_cases": len(raw_durations)},
        "sha256": {
            "ecopair": ecopair_sha,
            "ptlitig_patents": ptlit_patents_sha,
            "ptlitig_cases": ptlit_cases_sha,
        },
    }


# ---------- run_analysis helpers ----------
def compute_leniency(caseload, min_cases=MIN_CASES_PER_EXAMINER):
    by_ex = defaultdict(lambda: {"n": 0, "granted": 0})
    for (ex, au), v in caseload.items():
        by_ex[ex]["n"] += v["n"]
        by_ex[ex]["granted"] += v["granted"]
    return {ex: v["granted"] / v["n"] for ex, v in by_ex.items() if v["n"] >= min_cases}


def _inferential_triple(xs, ys, strata, bootstrap_seed, permutation_seed, min_n=50):
    """Compute ρ, bootstrap 95% CI, and within-stratum permutation p-value."""
    if len(xs) < min_n:
        return None
    rho = spearman_rho(xs, ys)
    boot = bootstrap_rank_correlation(xs, ys, BOOTSTRAP_RESAMPLES, bootstrap_seed)
    perm = within_stratum_permutation_pvalue(
        list(zip(xs, ys, strata)), rho, PERMUTATIONS, permutation_seed
    )
    return {
        "rho": rho,
        "boot_lo95": boot["lo95"],
        "boot_hi95": boot["hi95"],
        "perm_p_two_sided": perm["p_two_sided"],
        "null_mean": perm["null_mean"],
        "null_p025": perm["null_p025"],
        "null_p975": perm["null_p975"],
    }


def run_analysis(data):
    print("[6/6] Running analysis...", flush=True)
    caseload = data["examiner_caseload"]
    patent_lookup = data["patent_lookup"]
    litigated_patents = data["litigated_patents"]
    case_outcomes = data["case_outcomes"]
    duration_thresholds = data["duration_thresholds"]

    leniency_all = compute_leniency(caseload)
    vals = sorted(leniency_all.values())
    n_ex = len(vals)
    quartiles = [percentile(vals, q) for q in (0.25, 0.5, 0.75)]
    print(f"  leniency: {n_ex:,} examiners meeting min-caseload={MIN_CASES_PER_EXAMINER}", flush=True)
    print(f"  leniency quartiles: Q1={quartiles[0]:.3f}  Q2={quartiles[1]:.3f}  Q3={quartiles[2]:.3f}", flush=True)

    # Assemble patent-level records with LOO leniency and outcomes.
    records = []
    for pn, info in patent_lookup.items():
        ex, au = info["examiner"], info["stratum"]
        entry = caseload.get((ex, au))
        if entry is None or entry["n"] - 1 <= 0:
            continue
        loo = (entry["granted"] - 1) / (entry["n"] - 1)  # this patent was granted
        crids = litigated_patents.get(pn, [])
        outs = [case_outcomes[c] for c in crids if c in case_outcomes and case_outcomes[c]["is_patent_case"]]
        if not outs:
            continue
        valid_durations = [o["days_open"] for o in outs if o["days_open"] is not None]
        if not valid_durations:
            continue
        any_fast = 1 if any(o["fast_close"] for o in outs) else 0
        any_slow = 1 if any(o["slow_close"] for o in outs) else 0
        median_d = sorted(valid_durations)[len(valid_durations) // 2]
        log_d = math.log(1 + median_d)
        records.append({
            "patent": pn,
            "examiner": ex,
            "stratum_4digit": au,
            "loo_leniency": loo,
            "any_settle": any_fast,
            "any_judgment": any_slow,
            "median_days": median_d,
            "log_days": log_d,
            "n_cases": len(outs),
        })

    print(f"  analysis rows: {len(records):,}", flush=True)

    # --- Sensitivity across stratum aggregations ---
    sens_rows = []
    for label, fn in STRATUM_AGGREGATIONS:
        agg_caseload = defaultdict(lambda: {"n": 0, "granted": 0})
        for (ex, au), v in caseload.items():
            key = (ex, fn(au))
            agg_caseload[key]["n"] += v["n"]
            agg_caseload[key]["granted"] += v["granted"]
        cellsizes = Counter()
        x_vals, y_settle, y_judge, y_log = [], [], [], []
        strata = []
        for rec in records:
            agg_stratum = fn(rec["stratum_4digit"])
            entry = agg_caseload.get((rec["examiner"], agg_stratum))
            if not entry or entry["n"] - 1 <= 0:
                continue
            loo = (entry["granted"] - 1) / (entry["n"] - 1)
            x_vals.append(loo)
            y_settle.append(rec["any_settle"])
            y_judge.append(rec["any_judgment"])
            y_log.append(rec["log_days"])
            strata.append(agg_stratum)
            cellsizes[agg_stratum] += 1
        ok = [i for i, s in enumerate(strata) if cellsizes[s] >= MIN_PATENTS_PER_STRATUM]
        xs = [x_vals[i] for i in ok]
        ys_settle = [y_settle[i] for i in ok]
        ys_judge = [y_judge[i] for i in ok]
        ys_log = [y_log[i] for i in ok]
        ss = [strata[i] for i in ok]
        if len(xs) < 50:
            print(f"  [sens] {label}: only {len(xs)} retained rows; skipping", flush=True)
            continue
        log_block = _inferential_triple(xs, ys_log, ss, RANDOM_SEED + 100, RANDOM_SEED + 200)
        settle_block = _inferential_triple(xs, ys_settle, ss, RANDOM_SEED + 300, RANDOM_SEED + 400)
        judge_block = _inferential_triple(xs, ys_judge, ss, RANDOM_SEED + 500, RANDOM_SEED + 600)
        sens_rows.append({
            "stratum_aggregation": label,
            "n_patents": len(xs),
            "n_strata": len(set(ss)),
            "log_days": log_block,
            "settle": settle_block,
            "judgment": judge_block,
        })
        print(f"  [sens] {label}: n={len(xs):,}  ρ_log={log_block['rho']:+.3f} [{log_block['boot_lo95']:+.3f},{log_block['boot_hi95']:+.3f}] p={log_block['perm_p_two_sided']:.3f}  ρ_settle={settle_block['rho']:+.3f} p={settle_block['perm_p_two_sided']:.3f}  ρ_judge={judge_block['rho']:+.3f} p={judge_block['perm_p_two_sided']:.3f}", flush=True)

    # --- Settlement-selection diagnostic ---
    # Compare full-sample coefficient to adjudicated-subsample (slow-close) coefficient.
    # If weaker patents settle, conditioning on adjudicated hides the leniency effect.
    full_xs = [r["loo_leniency"] for r in records]
    full_log = [r["log_days"] for r in records]
    full_fast = [r["any_settle"] for r in records]
    full_strata = [r["stratum_4digit"] for r in records]

    adjud = [r for r in records if r["any_judgment"] == 1]
    adj_xs = [r["loo_leniency"] for r in adjud]
    adj_log = [r["log_days"] for r in adjud]
    adj_strata = [r["stratum_4digit"] for r in adjud]

    settled = [r for r in records if r["any_settle"] == 1]
    set_xs = [r["loo_leniency"] for r in settled]
    set_log = [r["log_days"] for r in settled]
    set_strata = [r["stratum_4digit"] for r in settled]

    # Keep only strata with >=MIN_PATENTS_PER_STRATUM inside each subsample.
    def _keep_by_stratum(xs, ys, strata):
        c = Counter(strata)
        keep = [i for i, s in enumerate(strata) if c[s] >= MIN_PATENTS_PER_STRATUM]
        return [xs[i] for i in keep], [ys[i] for i in keep], [strata[i] for i in keep]

    fx, fy, fs = _keep_by_stratum(full_xs, full_log, full_strata)
    ax, ay, as_ = _keep_by_stratum(adj_xs, adj_log, adj_strata)
    sx, sy, ss2 = _keep_by_stratum(set_xs, set_log, set_strata)

    full_block = _inferential_triple(fx, fy, fs, RANDOM_SEED + 700, RANDOM_SEED + 800)
    adj_block = _inferential_triple(ax, ay, as_, RANDOM_SEED + 900, RANDOM_SEED + 1000)
    set_block = _inferential_triple(sx, sy, ss2, RANDOM_SEED + 1100, RANDOM_SEED + 1200)

    # Probability that a patent settles as a function of leniency: Spearman ρ
    # between loo_leniency and any_settle on full sample. This is the heart of
    # the settlement-selection channel: if ρ ≠ 0, the adjudicated subsample is
    # a non-random selection of litigated patents.
    settle_sel_block = _inferential_triple(full_xs, full_fast, full_strata, RANDOM_SEED + 1300, RANDOM_SEED + 1400)

    # --- Negative control / falsification check ---
    # Replace the real outcome with a seeded-random vector that shares the
    # empirical mean and SD of log_days but is otherwise independent of
    # leniency by construction. A correct Spearman pipeline should return
    # ρ ≈ 0 here; any large |ρ| signals a ranking or alignment bug.
    rng_neg = random.Random(RANDOM_SEED + 9999)
    random_outcome = [rng_neg.random() for _ in range(len(records))]
    neg_control_rho = spearman_rho(full_xs, random_outcome)
    neg_control_boot = bootstrap_rank_correlation(full_xs, random_outcome, BOOTSTRAP_RESAMPLES, RANDOM_SEED + 1500)
    negative_control = {
        "description": "Spearman ρ between leniency and a seeded-random outcome; should be ≈ 0",
        "rho": neg_control_rho,
        "boot_lo95": neg_control_boot["lo95"],
        "boot_hi95": neg_control_boot["hi95"],
        "n": len(records),
    }
    print(f"  [negcontrol] ρ(leniency, random) = {neg_control_rho:+.4f} [{neg_control_boot['lo95']:+.4f},{neg_control_boot['hi95']:+.4f}]", flush=True)

    selection_diagnostic = {
        "full_sample_size": len(records),
        "frac_any_settle": sum(r["any_settle"] for r in records) / max(1, len(records)),
        "frac_any_judgment": sum(r["any_judgment"] for r in records) / max(1, len(records)),
        "full_leniency_vs_log_days": full_block,
        "adjudicated_leniency_vs_log_days": adj_block,
        "settled_leniency_vs_log_days": set_block,
        "full_leniency_vs_settle_prob": settle_sel_block,
    }

    # --- Top-line: use 4-digit native stratification ---
    top = next((s for s in sens_rows if s["stratum_aggregation"] == "art_unit_4digit"), None)

    # --- Decile table: does settlement probability vary by leniency decile? ---
    decile_table = []
    xs_all = [r["loo_leniency"] for r in records]
    ys_settle_all = [r["any_settle"] for r in records]
    ys_log_all = [r["log_days"] for r in records]
    if xs_all:
        ranks = rank_of(xs_all)
        n = len(xs_all)
        for d in range(10):
            lo, hi = d / 10.0, (d + 1) / 10.0
            if d == 9:
                mask = [i for i, rr in enumerate(ranks) if (rr - 1) / max(1, n - 1) >= lo]
            else:
                mask = [i for i, rr in enumerate(ranks) if lo <= (rr - 1) / max(1, n - 1) < hi]
            if not mask:
                continue
            settles = [ys_settle_all[i] for i in mask]
            logs = [ys_log_all[i] for i in mask]
            xs_d = [xs_all[i] for i in mask]
            decile_table.append({
                "decile": d + 1,
                "n": len(mask),
                "leniency_mean": sum(xs_d) / len(xs_d),
                "frac_settled": sum(settles) / len(settles),
                "mean_log_days": sum(logs) / len(logs),
            })

    return {
        "sha256": data["sha256"],
        "duration_thresholds": duration_thresholds,
        "parameters": {
            "min_cases_per_examiner": MIN_CASES_PER_EXAMINER,
            "min_patents_per_stratum": MIN_PATENTS_PER_STRATUM,
            "permutations": PERMUTATIONS,
            "bootstrap_resamples": BOOTSTRAP_RESAMPLES,
            "random_seed": RANDOM_SEED,
            "settle_quantile": SETTLE_QUANTILE,
            "judgment_quantile": JUDGMENT_QUANTILE,
        },
        "counts": {
            "examiners_meeting_min_caseload": n_ex,
            "litigated_patents_with_examiner_match": len(patent_lookup),
            "analysis_ready_records": len(records),
        },
        "leniency_distribution": {
            "q1": quartiles[0],
            "median": quartiles[1],
            "q3": quartiles[2],
            "mean": sum(vals) / len(vals) if vals else float("nan"),
        },
        "top_line": top,
        "sensitivity": sens_rows,
        "selection_diagnostic": selection_diagnostic,
        "decile_table": decile_table,
        "negative_control": negative_control,
        "limitations": [
            "Duration proxies conflate settlements, dismissals, and merits adjudications; the selection channel is directional, not structural.",
            "Within-art-unit exogeneity is assumed, not tested; Righi & Simcoe (2019) report departures in some units.",
            "Effect sizes are tiny (|ρ| ≈ 0.01); read the full→adjudicated attenuation, not any single p-value.",
            "Results apply to the 2020 PTLITIG/ECOPAIR vintages only; PTAB, ITC, and state matters are not fully covered.",
            "Adjudication proxy is duration-based, not merits-based; long-contested pre-trial settlements are mis-classified.",
            "Leniency is pooled across art units; an examiner is measured partly against their art-unit baseline, which the sensitivity sweep only partly mitigates.",
        ],
    }


# ---------- generate_report ----------
def generate_report(results, out_dir):
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    (out_dir / OUTPUT_RESULTS_JSON).write_text(json.dumps(results, indent=2, sort_keys=True))
    lines = []
    lines.append("# Examiner Leniency and Patent Litigation Resolution\n")
    c = results["counts"]
    p = results["parameters"]
    lines.append("## Counts\n")
    lines.append(f"- examiners meeting minimum caseload ({p['min_cases_per_examiner']}): {c['examiners_meeting_min_caseload']:,}")
    lines.append(f"- litigated utility patents with PatEx match: {c['litigated_patents_with_examiner_match']:,}")
    lines.append(f"- analysis-ready patents: {c['analysis_ready_records']:,}\n")

    lines.append("## Duration thresholds (empirical)\n")
    dt = results["duration_thresholds"]
    lines.append(f"- Q33={dt['q33']:.0f} days  Q67={dt['q67']:.0f} days  (n={dt['n_cases']:,} cases)\n")

    lines.append("## Leniency distribution\n")
    ld = results["leniency_distribution"]
    lines.append(f"- quartiles: Q1={ld['q1']:.3f}  Q2={ld['median']:.3f}  Q3={ld['q3']:.3f}")
    lines.append(f"- mean: {ld['mean']:.3f}\n")

    top = results.get("top_line")
    if top:
        lines.append("## Top-line (4-digit art-unit stratification)\n")
        lines.append(f"- n patents: {top['n_patents']:,} across {top['n_strata']:,} strata")
        for k, lab in [("log_days", "Log-days (continuous)"), ("settle", "Fast-close (settle proxy)"), ("judgment", "Slow-close (adjudication proxy)")]:
            b = top[k]
            lines.append(f"- {lab}: ρ = {b['rho']:+.4f}  95% CI [{b['boot_lo95']:+.4f}, {b['boot_hi95']:+.4f}]  perm p = {b['perm_p_two_sided']:.4f}")
        lines.append("")

    lines.append("## Sensitivity across art-unit aggregations\n")
    lines.append("| aggregation | n | strata | ρ_log [CI] p | ρ_settle [CI] p | ρ_judge [CI] p |")
    lines.append("|---|---:|---:|---|---|---|")
    for s in results["sensitivity"]:
        def fmt(b):
            return f"{b['rho']:+.3f} [{b['boot_lo95']:+.3f},{b['boot_hi95']:+.3f}] p={b['perm_p_two_sided']:.3f}"
        lines.append(f"| {s['stratum_aggregation']} | {s['n_patents']:,} | {s['n_strata']:,} | {fmt(s['log_days'])} | {fmt(s['settle'])} | {fmt(s['judgment'])} |")
    lines.append("")

    lines.append("## Decile table\n")
    lines.append("| decile | n | leniency mean | fraction fast-closed | mean log-days |")
    lines.append("|---:|---:|---:|---:|---:|")
    for d in results["decile_table"]:
        lines.append(f"| {d['decile']} | {d['n']:,} | {d['leniency_mean']:.3f} | {d['frac_settled']:.3f} | {d['mean_log_days']:.3f} |")
    lines.append("")

    lines.append("## Selection diagnostic\n")
    sd = results["selection_diagnostic"]
    lines.append(f"- full sample n={sd['full_sample_size']:,}")
    lines.append(f"- fraction fast-close (settle proxy): {sd['frac_any_settle']:.3f}")
    lines.append(f"- fraction slow-close (adjudication proxy): {sd['frac_any_judgment']:.3f}")
    for k, lab in [
        ("full_leniency_vs_log_days", "Full sample leniency→log-days"),
        ("adjudicated_leniency_vs_log_days", "Adjudicated-only leniency→log-days"),
        ("settled_leniency_vs_log_days", "Settled-only leniency→log-days"),
        ("full_leniency_vs_settle_prob", "Full sample leniency→P(fast-close)"),
    ]:
        b = sd.get(k)
        if b:
            lines.append(f"- {lab}: ρ = {b['rho']:+.4f}  95% CI [{b['boot_lo95']:+.4f}, {b['boot_hi95']:+.4f}]  perm p = {b['perm_p_two_sided']:.4f}")
    (out_dir / OUTPUT_REPORT_MD).write_text("\n".join(lines) + "\n")


# ---------- verify ----------
def verify(out_dir):
    out_dir = Path(out_dir)
    results_path = out_dir / OUTPUT_RESULTS_JSON
    assert results_path.exists(), f"missing results.json at {results_path}"
    with open(results_path) as f:
        r = json.load(f)

    # 1. Schema-level presence
    assert "top_line" in r and r["top_line"] is not None, "missing top_line"
    assert "sensitivity" in r and len(r["sensitivity"]) >= 2, "need ≥2 sensitivity rows"
    assert "decile_table" in r and len(r["decile_table"]) >= 5, "need ≥5 deciles"
    assert "counts" in r and r["counts"]["analysis_ready_records"] >= 500, \
        f"expected ≥500 analysis-ready patents; got {r['counts']['analysis_ready_records']}"

    # 2. Parameter integrity
    p = r["parameters"]
    assert p["permutations"] >= 1000, "must run ≥1000 permutations"
    assert p["bootstrap_resamples"] >= 1000, "must run ≥1000 bootstrap resamples"
    assert p["random_seed"] == RANDOM_SEED, "random seed mismatch"

    # 3. Leniency distribution in [0,1] and monotonic quartiles
    ld = r["leniency_distribution"]
    for k in ("q1", "median", "q3", "mean"):
        assert 0.0 <= ld[k] <= 1.0, f"leniency {k} out of [0,1]: {ld[k]}"
    assert ld["q1"] <= ld["median"] <= ld["q3"], "quartiles non-monotonic"

    # 4. Top-line: ρ in [-1,1], CI brackets ρ, p in (0,1], not implausibly large
    top = r["top_line"]
    for k in ("log_days", "settle", "judgment"):
        b = top[k]
        assert -1.0 <= b["rho"] <= 1.0, f"rho out of bounds: {b['rho']}"
        assert b["boot_lo95"] <= b["rho"] + 1e-9, f"boot lo above rho ({b['boot_lo95']} > {b['rho']})"
        assert b["boot_hi95"] + 1e-9 >= b["rho"], f"boot hi below rho ({b['boot_hi95']} < {b['rho']})"
        assert 0.0 < b["perm_p_two_sided"] <= 1.0, f"p-value out of bounds: {b['perm_p_two_sided']}"
        assert abs(b["rho"]) < 0.8, f"implausibly large Spearman ρ: {b['rho']}"

    # 5. Sensitivity row presence of all three granularities
    aggs = set(s["stratum_aggregation"] for s in r["sensitivity"])
    assert "art_unit_4digit" in aggs and "art_unit_3digit" in aggs and "art_unit_2digit" in aggs, \
        f"missing a stratum aggregation; saw {aggs}"

    # 6. Decile table: leniency_mean monotonic across deciles (sanity check)
    dec = r["decile_table"]
    for a, b in zip(dec, dec[1:]):
        assert a["leniency_mean"] <= b["leniency_mean"] + 1e-9, \
            f"decile leniency non-monotonic: {a['leniency_mean']} > {b['leniency_mean']}"

    # 7. Selection diagnostic: both full and adjudicated blocks present
    sd = r["selection_diagnostic"]
    assert 0.0 <= sd["frac_any_settle"] <= 1.0
    assert 0.0 <= sd["frac_any_judgment"] <= 1.0
    # Core selection test: fast-close fraction must be strictly positive
    # (otherwise the settlement channel is undefined in this vintage).
    assert sd["frac_any_settle"] > 0.0, \
        f"no fast-close cases detected; settlement proxy collapsed (frac={sd['frac_any_settle']})"
    assert sd["frac_any_judgment"] > 0.0, \
        f"no slow-close cases detected; adjudication proxy collapsed (frac={sd['frac_any_judgment']})"
    assert sd["full_leniency_vs_log_days"] is not None, "missing full-sample log-days block"
    assert sd["adjudicated_leniency_vs_log_days"] is not None, "missing adjudicated-only log-days block"

    # 8. SHA256 digests are well-formed
    sh = r["sha256"]
    for k in ("ecopair", "ptlitig_patents", "ptlitig_cases"):
        assert len(sh[k]) == 64 and all(c in "0123456789abcdef" for c in sh[k]), \
            f"sha256 for {k} not a valid hex digest"

    # 9. Duration thresholds in sensible order
    dt = r["duration_thresholds"]
    assert 0 < dt["q33"] < dt["q67"], f"duration quantiles non-monotonic: {dt}"

    # 10. Sensitivity: across all aggregations, ρ_log and ρ_settle must exist
    for s in r["sensitivity"]:
        assert s["log_days"] is not None, f"missing log_days block for {s['stratum_aggregation']}"
        assert s["settle"] is not None, f"missing settle block for {s['stratum_aggregation']}"

    # 11. Effect-size plausibility (Cohen's-d-style bound): |ρ| < 0.5 across all
    #     reported blocks. Spearman ρ in an observational design this large
    #     should be modest; a |ρ| above 0.5 indicates an alignment or ranking bug.
    for s in r["sensitivity"]:
        for k in ("log_days", "settle", "judgment"):
            b = s[k]
            assert abs(b["rho"]) < 0.5, f"implausibly large |ρ| for {s['stratum_aggregation']}/{k}: {b['rho']}"

    # 12. CI width sanity: the bootstrap CI must be non-degenerate (width > 0).
    #     We do not require CI_width > 1% of |ρ| uniformly because for very small
    #     true effects the absolute CI can be wider than the estimate; but the CI
    #     must be strictly positive width on every block.
    for s in r["sensitivity"]:
        for k in ("log_days", "settle", "judgment"):
            b = s[k]
            ci_width = b["boot_hi95"] - b["boot_lo95"]
            assert ci_width > 0.0, f"degenerate CI width for {s['stratum_aggregation']}/{k}: {ci_width}"
            # CI width should be at least 1% of the estimate OR both endpoints close to zero
            assert ci_width >= 0.01 * max(abs(b["rho"]), 1e-6) or abs(b["rho"]) < 0.01, \
                f"suspiciously tight CI for {s['stratum_aggregation']}/{k}: width={ci_width}, ρ={b['rho']}"

    # 13. Permutation-null centering (exchangeability sanity): the within-stratum
    #     shuffle null mean ρ should be near zero. A grossly off-center null
    #     indicates the stratum variable is badly collinear with the outcome.
    for s in r["sensitivity"]:
        for k in ("log_days", "settle", "judgment"):
            b = s[k]
            assert abs(b["null_mean"]) < 0.05, \
                f"permutation null not centered for {s['stratum_aggregation']}/{k}: null_mean={b['null_mean']}"

    # 14. Sign-stability sensitivity: the log_days ρ should have the same sign
    #     across all three stratum aggregations. This is the robustness claim
    #     made in the paper; the verification enforces it.
    log_signs = [(b_ := s["log_days"])["rho"] for s in r["sensitivity"]]
    assert all(x < 0 for x in log_signs) or all(x > 0 for x in log_signs) or all(abs(x) < 0.005 for x in log_signs), \
        f"log_days ρ sign not stable across aggregations: {log_signs}"

    # 15. Negative control / falsification: ρ between leniency and a seeded
    #     random outcome must be small. Large |ρ| here indicates a bug in the
    #     Spearman or ranking code, not a real finding.
    assert "negative_control" in r, "missing negative_control block"
    nc = r["negative_control"]
    assert abs(nc["rho"]) < 0.05, \
        f"negative control failed: |ρ| = {abs(nc['rho'])} ≥ 0.05 (expected ≈ 0 for a random outcome)"
    assert nc["n"] == r["counts"]["analysis_ready_records"], \
        "negative-control n must match the analysis-ready record count"

    # 16. Limitations block present and non-trivial
    assert "limitations" in r and isinstance(r["limitations"], list) and len(r["limitations"]) >= 4, \
        f"need ≥ 4 limitations stated in results.json; got {len(r.get('limitations', []))}"

    # 17. Counts consistency
    assert r["counts"]["litigated_patents_with_examiner_match"] >= r["counts"]["analysis_ready_records"], \
        "analysis_ready_records cannot exceed litigated_patents_with_examiner_match"
    assert r["counts"]["examiners_meeting_min_caseload"] >= 100, \
        "need at least 100 examiners meeting the minimum-caseload floor"

    print("all 17 verification assertions passed")
    print("ALL CHECKS PASSED")
    return True


# ---------- Custom exceptions with exit codes ----------
class DownloadError(RuntimeError):
    exit_code = 2

class ShaMismatchError(RuntimeError):
    exit_code = 3

class SchemaError(RuntimeError):
    exit_code = 4


# ---------- main ----------
def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--cache-dir", default=".cache", help="where to cache raw data zips")
    ap.add_argument("--verify", action="store_true", help="verify results.json only")
    args = ap.parse_args()

    # Seed the module-level random state so any unseeded incidental use is
    # deterministic. All statistical routines take an explicit seed separately.
    random.seed(RANDOM_SEED)

    here = Path(__file__).parent.resolve()
    cache_dir = (here / args.cache_dir).resolve()

    if args.verify:
        try:
            verify(here)
        except AssertionError as e:
            print(f"ERROR: verification failed: {e}", file=sys.stderr, flush=True)
            sys.exit(5)
        except FileNotFoundError as e:
            print(f"ERROR: results.json missing — run the analysis first: {e}", file=sys.stderr, flush=True)
            sys.exit(6)
        return

    t_start = time.time()
    try:
        data = load_data(cache_dir)
    except (urllib.error.URLError, urllib.error.HTTPError, TimeoutError) as e:
        print(f"ERROR: network download failed: {e}", file=sys.stderr, flush=True)
        print("Hint: check internet connectivity and Wayback Machine reachability.", file=sys.stderr, flush=True)
        sys.exit(2)
    except DownloadError as e:
        print(f"ERROR: {e}", file=sys.stderr, flush=True)
        sys.exit(e.exit_code)
    except ShaMismatchError as e:
        print(f"ERROR: SHA256 mismatch on cached data: {e}", file=sys.stderr, flush=True)
        print("Hint: delete the file from the --cache-dir and rerun, or check for upstream data drift.", file=sys.stderr, flush=True)
        sys.exit(e.exit_code)
    except SchemaError as e:
        print(f"ERROR: upstream data schema changed: {e}", file=sys.stderr, flush=True)
        sys.exit(e.exit_code)
    except RuntimeError as e:
        msg = str(e)
        if "Failed to download" in msg:
            print(f"ERROR: {msg}", file=sys.stderr, flush=True)
            sys.exit(2)
        if "SHA256 mismatch" in msg:
            print(f"ERROR: {msg}", file=sys.stderr, flush=True)
            sys.exit(3)
        if "Unexpected" in msg and "schema" in msg:
            print(f"ERROR: {msg}", file=sys.stderr, flush=True)
            sys.exit(4)
        print(f"ERROR: unexpected runtime failure: {msg}", file=sys.stderr, flush=True)
        sys.exit(1)
    except OSError as e:
        print(f"ERROR: filesystem error (is --cache-dir writable?): {e}", file=sys.stderr, flush=True)
        sys.exit(7)

    results = run_analysis(data)
    generate_report(results, here)
    print(f"ANALYSIS COMPLETE in {time.time()-t_start:.1f}s", flush=True)


if __name__ == "__main__":
    main()
SCRIPT_EOF
```

**Expected output:** `analyze.py` written to workspace (file size ~24 KB).

**Success criteria:** `/tmp/claw4s_auto_examiner-harshness-and-litigation-outcomes/analyze.py` exists and is non-empty.
**Failure condition:** heredoc truncation or missing closing delimiter.

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_examiner-harshness-and-litigation-outcomes && python3 analyze.py
```

**Expected stdout (first run, network-bound):**

```
[1/6] Downloading PTLITIG patents (case-to-patent link)...
  downloaded ptlitig_patents.csv.zip: 3,499,585 bytes in ...
[2/6] Downloading PTLITIG cases (case outcomes)...
  downloaded ptlitig_cases.csv.zip: 6,046,617 bytes in ...
[3/6] Downloading PatEx ECOPAIR application_data (~830 MB; 15-30 min on first run)...
  downloaded ecopair_application_data.csv.zip: 868,686,226 bytes in ...
[4/6] Parsing litigation dockets...
  empirical duration thresholds: Q33=<~135>d  median=<~232>d  Q67=<~425>d
  n patent-adjacent cases with valid duration: ~58,000
  parsed ~97,000 cases
  parsed ~50,000 unique litigated utility patents
[5/6] Streaming PatEx application_data.csv (11M+ rows)...
  ... progress messages ...
  ~850,000 unique (examiner,art-unit) pairs
  ~50,000 litigated patents with matched PatEx record
[6/6] Running analysis...
  leniency: ~15,000 examiners meeting min-caseload=20
  leniency quartiles: Q1=... Q2=... Q3=...
  analysis rows: ~45,000
  [sens] art_unit_4digit: n=... ρ_log=... p=... ρ_settle=... p=... ρ_judge=... p=...
  [sens] art_unit_3digit: ...
  [sens] art_unit_2digit: ...
ANALYSIS COMPLETE in ... s
```

**Expected outputs on disk:** `results.json` (structured), `report.md` (readable), `.cache/*.csv.zip` (three data files).

**Success criteria:**
- stdout contains `ANALYSIS COMPLETE`
- `results.json` and `report.md` exist and are non-empty

**Failure conditions:**
- Network error on any of the three downloads (script retries 5× with exponential backoff and aborts on the 6th)
- SHA256 mismatch on a cached file (script deletes and re-downloads once; mismatch after re-download is a hard abort)
- Unexpected schema in ECOPAIR (column not found → `RuntimeError` with the observed header)

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_examiner-harshness-and-litigation-outcomes && python3 analyze.py --verify
```

**Expected stdout:**

```
all 17 verification assertions passed
ALL CHECKS PASSED
```

**Success criteria:** exit code 0 and the stdout messages above.
**Failure conditions:**
- any `AssertionError` from the `verify()` function is caught, the offending message is printed to **stderr**, and the script exits with **code 5**.
- If `results.json` is missing entirely, the script exits with **code 6**.
- If the cache directory is not writable or some other OS-level error occurs, exit code **7**.

The 17 verification assertions cover: (1) schema presence and ≥500 analysis-ready records, (2) parameter integrity (≥1000 perms, ≥1000 bootstraps, seed=42), (3) leniency distribution in [0,1] with monotone quartiles, (4) top-line ρ in [-1,1] with CI bracketing and p in (0,1], (5) all three stratum aggregations present, (6) decile-table monotonicity in leniency-mean, (7) selection-diagnostic fractions in (0,1], (8) SHA256 hex-digest well-formedness on all three data files, (9) duration quantile ordering, (10) sensitivity blocks populated, (11) effect-size plausibility `|ρ| < 0.5` on every sensitivity block, (12) CI width non-degenerate and at least 1% of `|ρ|` (or both endpoints near zero), (13) permutation null distribution centered within 0.05 of zero (exchangeability sanity), (14) sign stability of log-days ρ across all three stratum aggregations, (15) negative-control (seeded-random outcome) `|ρ| < 0.05` as a falsification check, (16) limitations block with ≥ 4 stated caveats, (17) counts consistency.