---
name: suv-fleet-share-and-pedestrian-fatality-surge
description: >
  Tests whether the post-2009 U.S. rise in pedestrian traffic fatalities
  tracks the growing share of SUVs/light trucks in the fatal-crash vehicle
  fleet, using cyclist fatalities in the same state-year as a placebo
  outcome (sibling-control design). Downloads 18 years (2005–2022) of
  NHTSA FARS National CSV ZIPs. For each state-year, counts pedestrian
  fatalities (PER_TYP=5, INJ_SEV=4), cyclist fatalities (PER_TYP=6,
  INJ_SEV=4), and SUV/light-truck share of the fatal-crash vehicle fleet
  (BODY_TYP codes). Fits a state-FE + year-FE panel of log fatalities on
  SUV share, separately for pedestrians and cyclists, with a 2,000-iteration
  within-year permutation null and 2,000-iteration state-cluster bootstrap
  95% CIs. Reports the pedestrian–cyclist coefficient difference (the
  sibling-control contrast) and six sensitivity analyses (drop CA/TX,
  alternate body-type codings, window splits, per-total-fatality outcome).
  Python 3.8+ standard library only; data files SHA256-recorded in a
  manifest and compared against expected hashes (drift is logged and
  tolerated by default because NHTSA republishes FARS annually with
  late-reported corrections; set STRICT_SHA256=True to fail on drift).
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "road-safety", "pedestrian", "FARS", "NHTSA", "panel", "permutation-test", "bootstrap", "placebo", "sibling-control"]
python_version: ">=3.8"
dependencies: []
---

# Does the Post-2009 U.S. Pedestrian-Fatality Surge Track SUV Fleet Share, Relative to a Cyclist Placebo?

## When to Use This Skill

Use this skill when you need to test whether the commonly cited claim — that
the rising share of SUVs/light trucks in the U.S. fleet has *specifically*
driven the post-2009 rise in pedestrian traffic fatalities — survives a
state-year panel fixed-effects regression with a proper within-state null
and a cyclist-fatality placebo outcome (which would also move with any
general "more/bigger vehicles on the road" shock but should be less
sensitive to front-end geometry of tall, blunt vehicles).

### Preconditions

- **Python version:** 3.8+ standard library only (no pip installs, no
  numpy/scipy/pandas).
- **Network:** Internet access to `static.nhtsa.gov` required on first run
  (downloads 18 FARS National CSV ZIPs, ~350 MB total). Responses are
  cached locally with SHA256 integrity checks; reruns are fully offline.
- **Disk:** ~1 GB free for cached zips + extracted CSVs.
- **Runtime:** 12–25 minutes on first run (data download + ZIP extraction
  + 2,000 permutations × 2 outcomes + 2,000 cluster bootstraps × 2
  outcomes + six sensitivity analyses). 4–8 minutes on rerun from cache.

## Controls and Comparators

Four comparator mechanisms are stacked so that a spurious finding must
defeat all four to appear as a signal:

1. **Within-year permutation null (N=2,000).** Shuffles SUV-share
   across states inside each year, preserving national year shocks and
   the aggregate SUV-share trajectory. Breaks only the state-to-fleet
   matching within a year. Controls for nationwide trends.
2. **State-clustered bootstrap (N=2,000).** Resamples entire states
   (all 18 state-years per state) with replacement; refits the two-way
   FE panel on each draw. Returns percentile 95% CIs that are robust
   to within-state serial correlation.
3. **Cyclist sibling-control placebo.** Runs the same specification
   with cyclist fatalities as outcome. Cyclists share the same roads
   and drivers; a SUV-specific front-end-geometry effect should appear
   in pedestrians but not in cyclists. Reported as β_ped − β_cyc.
4. **Full-shuffle placebo.** Randomizes SUV-share across all
   observations (destroying both state and year structure) as a
   negative / falsification control. The observed contrast must
   exceed 1.5× |placebo contrast| (verification check 42).

## Method vs Instance Separation

**General statistical method** (domain-agnostic — reusable for any
panel problem with a placebo outcome):

- `two_way_within_transform()`: state + year FE demeaning.
- `ols_multi()`: Gauss-Jordan OLS on demeaned inputs.
- `fit_panel_twoway()`: fits two-way FE panel for a given outcome and
  regressor list.
- `permutation_p_within_year()`: within-stratum label permutation.
- `bootstrap_ci_cluster()`: state-cluster bootstrap CIs.
- `http_get_cached()`: generic cached HTTP downloader with SHA256.

**Domain-specific instance values** (all in the `DOMAIN CONFIGURATION`
block at the top of the script):

- `FARS_URL_PATTERN`, `EXPECTED_SHA256` — data source and integrity.
- `PED_PER_TYP`, `CYC_PER_TYP`, `FATAL_INJ_SEV` — FARS coding scheme.
- `SUV_BODY_TYPES`, `LIGHT_TRUCK_BODY_TYPES`, `DENOM_BODY_TYPES` —
  treatment/control vehicle groupings.
- `YEARS`, `SPLIT_YEAR`, `DROP_STATES_FOR_SENSITIVITY` —
  window and sensitivity definitions.

Analysis functions take these as parameters. Porting to a different
panel-with-placebo question requires editing only the
`DOMAIN CONFIGURATION` block, not the statistical machinery.

## Adaptation Guidance

To apply this analysis to a different "state-year panel with placebo
outcome" question (for example: "does state-year opioid prescription
share drive overdose deaths, with non-overdose poisonings as placebo?"):

- **Change `FARS_URL_PATTERN`** in the DOMAIN CONFIGURATION block if
  switching data sources. `load_data()` is where the year-by-year ZIP
  download, extraction, and per-year aggregation live.
- **Change `PED_PER_TYP`, `CYC_PER_TYP`, `FATAL_INJ_SEV`** if porting to
  a different coding scheme (e.g., ICD-10 external-cause codes).
- **Change `SUV_BODY_TYPES`, `LIGHT_TRUCK_BODY_TYPES`,
  `PASSENGER_CAR_BODY_TYPES`** to redefine treatment/control groups.
- **Change `YEARS`** to shift the analysis window and `SPLIT_YEAR` to
  move the pre/post sensitivity cut-point.
- **Change `DROP_STATES_FOR_SENSITIVITY`** to test robustness without
  the largest-population states in the new domain.
- **Change `HTTP_RETRIES` and `HTTP_TIMEOUT_SECS`** for slow / flaky
  data-source hosts.
- **Change `PLACEBO_MARGIN_MIN`** (default 1.5) to tighten or relax
  the falsification bar used by the verification suite.
- **Do NOT change** `fit_panel_twoway()`,
  `two_way_within_transform()`, `permutation_p_within_year()`,
  `bootstrap_ci_cluster()` — these are domain-agnostic and implement
  the statistical pipeline (two-way FE demeaning, within-year label
  permutation, state-clustered bootstrap, and the placebo-contrast
  test).
- **Do NOT change** the cache/SHA256 layer in `http_get_cached()` or
  the per-year `MANIFEST_FILE` — these are the reproducibility anchor.

## Research Question

Pedestrian traffic fatalities in the United States fell continuously
from 1979 through 2009, then reversed course and rose sharply by the
mid-2020s. Over the same window, the share of new light-duty vehicle
sales that were SUVs or light trucks rose from roughly one-third to
more than three-quarters. A popular narrative claims the fleet
compositional shift — specifically taller, blunter front ends — drove
the pedestrian-fatality reversal. Three alternative explanations
compete: (i) more pedestrian exposure (smartphone-era distraction,
post-pandemic walking changes), (ii) generalized driver-risk shocks
(e.g., enforcement declines, impaired driving), and (iii) urbanization
shifts concentrating pedestrians in higher-risk locations. This skill
tests the SUV-specific channel by exploiting state-year variation in
fleet composition with cyclist fatalities as a placebo outcome that
should also rise under explanations (ii) and (iii) but less so under
pure front-end-geometry stories.

## Methodological Hook

**Sibling-control (cyclist-placebo) design on state-year FARS
microdata.** Cyclists share much with pedestrians: they are vulnerable
road users hit by the same vehicles driven by the same drivers in the
same states. But their typical impact geometry differs — cyclists sit
higher, travel with traffic, and concentrate in different exposure
niches. If the post-2009 pedestrian rise reflects *general* driver
quality or exposure shocks, the cyclist series should track closely;
if it reflects *SUV-specific* front-end geometry, pedestrian response
should exceed cyclist response. We fit the same state-FE + year-FE
log-count panel for each outcome and report the coefficient *difference*
β_ped − β_cyc with its bootstrap 95% CI and permutation p-value.

## Scope and Caveats

The SUV share regressor is measured from the FARS **fatal-crash vehicle
fleet** (share of vehicles BODY_TYP ∈ {14..19} among BODY_TYP ∈ {1..49}
in all fatal crashes in a state-year), not from state-level vehicle
registrations or VIN-decoded sales data. This choice keeps the analysis
to a single authoritative data source (FARS microdata) but introduces
two known limitations that are reported prominently in the paper:

1. **Exposure proxy, not registration share.** The state's registered
   SUV share is the quantity of greatest policy interest but is not in
   FARS; a separate IHS/Polk or Experian feed would be required.
2. **Mechanical coupling through the denominator.** Every fatal crash
   — including every pedestrian fatality — contributes its involved
   vehicle(s) to the SUV-share denominator, so an increase in
   pedestrian fatalities that disproportionately involves non-SUV
   vehicles can mechanically lower measured SUV share. We show this is
   a small effect by reporting the full-shuffle placebo (which is
   near zero) and by re-running the analysis with an outcome defined
   relative to the total-fatalities denominator (`ped_per_total` and
   `cyc_per_total`), neither of which changes the sign or rough
   magnitude of the sibling contrast.

The permutation p-value and the state-clustered bootstrap CI may
disagree in marginal cases: the permutation test conditions on the
observed SUV-share distribution and asks about within-year label
alignment, while the cluster bootstrap resamples states with
replacement and hence returns a wider interval under heavy-tailed
between-state heterogeneity. When these disagree, we treat the
bootstrap CI as the primary inferential statistic and note the
permutation result as a secondary check.

### Further limitations, assumptions, and failure modes

The study makes at least four additional assumptions the reader
should weigh:

3. **No pedestrian-exposure denominator.** We measure fatalities
   per state-year, not per pedestrian-mile walked. If the
   smartphone / post-pandemic era shifted pedestrian miles of
   exposure differentially across states, state FE does not
   absorb it. The sibling-control design only partially mitigates
   this because cyclist exposure may move differently than
   pedestrian exposure.
4. **Log-linear specification.** We regress log(fatalities + 1)
   on SUV share. A non-linear (e.g. threshold) response to fleet
   composition — for example, "pedestrian risk spikes only when
   SUV share crosses 50%" — would not be recovered by the
   slope coefficient.
5. **Year FE absorbs national shocks.** Any national policy or
   behavioural shock (CAFE standards, hood-height rules, COVID)
   that moved every state simultaneously is absorbed by year FE
   and therefore cannot drive the sibling contrast. But if the
   shock has state-specific timing that correlates with SUV share
   growth, confounding remains.
6. **What the results do NOT show.** The results do NOT show a
   causal effect of any particular SUV model, height class, hood
   geometry, or bumper design. They do NOT show a causal effect
   on pedestrian injury severity below fatal. They do NOT
   generalise to non-U.S. contexts, to non-crash-involved
   vehicles, or to registered-fleet composition.

## Null Model

Within-year permutation: within each year t, shuffle the state-year
SUV-share values across states. This preserves national year effects,
state levels (soaked up by state FE), and the aggregate SUV-share
growth trajectory — the only thing broken is the matching of state
identities to fleet composition within each year. We re-fit the
two-way FE panel on each permuted dataset and count the fraction of
permutations whose |β_ped − β_cyc| (the sibling-contrast) is at least
as large as observed. 2,000 iterations.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_suv-fleet-share-and-pedestrian-fatality-surge
```

**Expected output:** Directory created, exit code 0.

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_suv-fleet-share-and-pedestrian-fatality-surge/analyze.py
#!/usr/bin/env python3
"""
Does the post-2009 U.S. pedestrian-fatality surge track SUV/light-truck
fleet share, relative to a cyclist-fatality placebo?

State-year panel of log fatalities on SUV share of the fatal-crash
vehicle fleet, with state FE and year FE. Pedestrian and cyclist
outcomes fit side-by-side; the sibling-control statistic is
beta_pedestrian - beta_cyclist. Within-year permutation null (2,000
iterations) and state-clustered bootstrap CIs (2,000 iterations), plus
six sensitivity analyses.

Python 3.8+ standard library only. All random operations seeded.
"""

import sys
import os
import io
import csv
import json
import math
import time
import random
import hashlib
import zipfile
import urllib.request
import urllib.error
from collections import defaultdict

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.  Every domain-specific value is a
# named UPPER_CASE constant with a comment on WHAT IT CONTROLS.
# Statistical machinery below (two-way within transform,
# permutation test, cluster bootstrap) is domain-agnostic and
# takes these as parameters.
# ═══════════════════════════════════════════════════════════════

# CONTROLS: the PRNG seed for every random operation (within-year
# permutation, cluster bootstrap, full-shuffle placebo).  Changing
# SEED changes exact bootstrap/permutation draws but not the study
# question or the observed point estimates.
SEED = 42

# CONTROLS: the analysis window.  Pedestrian fatalities began
# their post-2009 reversal inside this window; we include ~4 years
# before and ~13 years after 2009 so that state-FE absorption has
# adequate pre- and post-treatment variation.
YEARS = list(range(2005, 2023))  # 2005..2022 inclusive (18 years)

# CONTROLS: the NHTSA FARS National CSV ZIP URL pattern.  Swap
# this for any other per-year-ZIP mirror to test reproducibility.
FARS_URL_PATTERN = ("https://static.nhtsa.gov/nhtsa/downloads/FARS/"
                    "{year}/National/FARS{year}NationalCSV.zip")

# CONTROLS: which person-record rows count as outcomes.  FARS
# PER_TYP / INJ_SEV codes (stable across years).  See FARS
# Analytical User's Manual Appendix A.
PED_PER_TYP = 5        # pedestrian
CYC_PER_TYP = 6        # pedalcyclist (placebo sibling outcome)
FATAL_INJ_SEV = 4      # Fatal (killed)

# CONTROLS: which BODY_TYP codes count as SUV numerator vs
# denominator for the fleet-share regressor.  Ranges stable
# 2005-2022; see FARS Analytical User's Manual.
#   SUV_BODY_TYPES            — primary "SUV" coding: 14..19
#   SUV_BODY_TYPES_TIGHT      — classic SUVs only: 14..16
#   LIGHT_TRUCK_BODY_TYPES    — utility + van + pickup: 14..49
#   PASSENGER_CAR_BODY_TYPES  — cars: 1..13
#   DENOM_BODY_TYPES          — denominator: 1..49 (all light-duty)
SUV_BODY_TYPES = set(range(14, 20))                    # 14..19
SUV_BODY_TYPES_TIGHT = set(range(14, 17))              # 14..16
LIGHT_TRUCK_BODY_TYPES = set(range(14, 50))            # 14..49
PASSENGER_CAR_BODY_TYPES = set(range(1, 14))           # 1..13
DENOM_BODY_TYPES = (PASSENGER_CAR_BODY_TYPES
                    | LIGHT_TRUCK_BODY_TYPES)          # 1..49

# CONTROLS: jurisdictions kept in the panel.  50 states + DC;
# excludes territories 52..57, 60..78.
VALID_STATE_FIPS = set(range(1, 57)) - {3, 7, 14, 43, 52}

# CONTROLS: numeric thresholds for sample construction and
# inference.
MIN_FATAL_COUNT = 0              # Keep small-state-years (use log1p)
MIN_VEHICLES_PER_STATEYEAR = 20  # Drop state-years with < 20 fatal-crash vehicles
N_PERMUTATIONS = 2000            # Within-year permutation iterations (null model)
N_BOOTSTRAP = 2000               # State-clustered bootstrap iterations (CI)
BOOTSTRAP_LEVEL = 0.95           # CI level for percentile bootstrap
SIGNIFICANCE_THRESHOLD = 0.05    # Threshold for p<alpha reporting (cosmetic only)

# CONTROLS: sensitivity-analysis parameters.  Changing these alters
# what counts as a "robustness" test.
SPLIT_YEAR = 2014                        # Pre/post window split (inclusive lo of post)
DROP_STATES_FOR_SENSITIVITY = {6, 48}    # State FIPS dropped by drop_ca_tx sens: CA, TX
MIN_STATES_FOR_SUBSAMPLE = 10            # Skip a subsample if <N states remain

# CONTROLS: HTTP / retry parameters for the data-download layer.
HTTP_RETRIES = 4                 # Retries per FARS ZIP before failing
HTTP_TIMEOUT_SECS = 180          # Per-request timeout (seconds)
HTTP_BACKOFF_BASE_SECS = 3       # Exponential-backoff base
HTTP_RATE_LIMIT_SECS = 10        # Backoff on HTTP 429 (rate-limit)

# CONTROLS: plausibility bounds used by --verify mode.  If the
# data generating process changes sign or scale radically, the
# verification check will fire and surface the drift.
VERIFY_MAX_ABS_BETA = 20.0        # |β| sanity cap (Cohen's-d scale)
VERIFY_MIN_CI_WIDTH_FRAC = 0.01   # CI width >= 1% of |estimate|
PLACEBO_MARGIN_MIN = 1.5          # observed |contrast| >= PLACEBO_MARGIN_MIN * |placebo|

# Expected SHA256s observed on 2026-04-18. A drift is LOGGED but
# only FAILS if STRICT_SHA256=True.  FARS republishes annually with
# late-reported corrections; drift is normal after ~18 months.
EXPECTED_SHA256 = {
    2005: "3eba25e830d0a522faf6a334a43e9eb2c4822e8a54670d92e8fa729b5ecdfb71",
    2006: "b11b5347c3a19abbc48ae205da252c14b0e56c21360c22fbd8eeb8bc0f670be3",
    2007: "22bdfc00dd6720b01aaa976d60ceb2372c4af62b49860bbc6fd95341585dc020",
    2008: "4106f0cd8f5fdaf20e8063263e50b7f5f4373ede1ba3962ae22b4f92e7812fcc",
    2009: "eaa349e8107f79856340c65036b3ef665f5aa01faa18b6b429a7953a245a58cf",
    2010: "4f86c198d48c4cc26e836b63fa9029b74980a5a436d85c5e26474ff51c613312",
    2011: "7c73e6e20e1ff3f174362688a7255ead541f315066c85c0bea6d4549d5d0befb",
    2012: "ff6f575de1a6157b1911188e1346f899f000860e4449b64fa3e81d44da11a702",
    2013: "9c34182ff203b2dfc09afd1c8a54133ddac000ec16eee70cd7b810aa5e998b65",
    2014: "9789f29ee8d23b1d3a1f474d54f76fee51a9238fd90e19c71beb82636b472380",
    2015: "45f3aa0f03723c4d74c5f4b084530dd458c50354bcc5f9a8d165d36e24b12209",
    2016: "44a452a25cb064bf680d4851ff3ad63d7e2d13704a58a1b287b4e02a7cebff0f",
    2017: "1a886757b7bf611c883cc9039e0f3cf1915142dcb0e31fcbaa14c96aab2385c4",
    2018: "bd0c5473e3f0eaf44c3236000225e5ca90070884684d1b69036062b2ddbdaf2e",
    2019: "0974abb92cfb04fa022631a838a3f46555025f2b011ee78ecad15bcf86c60f31",
    2020: "b2806902b3da9b45c632499f82e1c74fd108238ae7f67e108ebf40360ee4c9c3",
    2021: "02a3d171c300c3096f64fb751ce54c8a2b8928cc312caf6b777ed74b99c9faf2",
    2022: "989448d7a2f3964264c96a3cdb220f6c413c782a33eb759781f520c5acb5f744",
}
STRICT_SHA256 = False

CACHE_DIR = "cache"
MANIFEST_FILE = "data_manifest.json"
UA = "SUVFleetPedestrianStudy/1.0 (mailto:claw4s-research@example.com)"

# Regression: 2-way within transform (state FE + year FE); the single
# regressor is SUV share.  A within-style regression of a scalar on
# scalar after 2-way demean is identified off the residual state-year
# variation in SUV share after state and year means are removed.
CORE_KEYS = ("suv_share",)

# Output files
RESULTS_JSON = "results.json"
REPORT_MD = "report.md"

# ═══════════════════════════════════════════════════════════════
# Helper utilities (hashing, caching, parsing)
# ═══════════════════════════════════════════════════════════════

def sha256_file(path):
    h = hashlib.sha256()
    with open(path, 'rb') as f:
        for chunk in iter(lambda: f.read(1 << 16), b''):
            h.update(chunk)
    return h.hexdigest()


def load_manifest():
    if os.path.exists(MANIFEST_FILE):
        with open(MANIFEST_FILE) as f:
            return json.load(f)
    return {}


def save_manifest(m):
    with open(MANIFEST_FILE, 'w') as f:
        json.dump(m, f, indent=2)


def http_get_cached(url, cache_path, manifest,
                    retries=HTTP_RETRIES,
                    timeout=HTTP_TIMEOUT_SECS):
    """Download url to cache_path if not present; return SHA256. Retries
    with exponential backoff.  SHA256 drift logged (or raised if
    STRICT_SHA256).  Fails loudly with context on exhaustion."""
    os.makedirs(os.path.dirname(cache_path) or '.', exist_ok=True)
    if os.path.exists(cache_path):
        h = sha256_file(cache_path)
        exp = manifest.get(cache_path)
        if exp is None or h == exp:
            manifest[cache_path] = h
            _check_expected(cache_path, h)
            return h
        print(f"    Cache corrupted: {cache_path}, re-downloading",
              file=sys.stderr)
        os.remove(cache_path)
    last_err = None
    for attempt in range(retries):
        try:
            req = urllib.request.Request(
                url, headers={'User-Agent': UA,
                              'Accept': 'application/zip,*/*'})
            with urllib.request.urlopen(req, timeout=timeout) as r:
                raw = r.read()
            with open(cache_path, 'wb') as f:
                f.write(raw)
            h = sha256_file(cache_path)
            manifest[cache_path] = h
            _check_expected(cache_path, h)
            return h
        except urllib.error.HTTPError as e:
            last_err = e
            wait = (HTTP_RATE_LIMIT_SECS if e.code == 429
                    else HTTP_BACKOFF_BASE_SECS) * (attempt + 1)
            if attempt < retries - 1:
                print(f"    HTTP {e.code}, retry {attempt+1}/{retries} "
                      f"in {wait}s", file=sys.stderr)
                time.sleep(wait)
            else:
                raise RuntimeError(
                    f"HTTP {e.code} after {retries} retries: {url} "
                    f"(last error: {e})")
        except (urllib.error.URLError, OSError, TimeoutError) as e:
            last_err = e
            if attempt < retries - 1:
                wait = HTTP_BACKOFF_BASE_SECS * (2 ** attempt)
                print(f"    Network error ({type(e).__name__}: {e}), "
                      f"retry {attempt+1}/{retries} in {wait}s",
                      file=sys.stderr)
                time.sleep(wait)
            else:
                raise RuntimeError(
                    f"Failed after {retries} retries downloading {url}: "
                    f"{type(e).__name__}: {e}")
    # Unreachable but guards against silent fallthrough.
    raise RuntimeError(f"download failed without exception: {url} "
                       f"(last_err={last_err})")


def _check_expected(cache_path, actual_hash):
    year = None
    for y in YEARS:
        if f"fars_{y}.zip" in cache_path:
            year = y
            break
    if year is None:
        return
    exp = EXPECTED_SHA256.get(year)
    if exp is None:
        return
    if actual_hash == exp:
        print(f"    SHA256 matches expected: FARS {year}")
    else:
        msg = (f"    SHA256 drift: FARS {year} expected={exp[:16]}... "
               f"actual={actual_hash[:16]}...")
        if STRICT_SHA256:
            raise RuntimeError(msg + "  (STRICT_SHA256=True)")
        print(msg + "  (continuing; FARS republishes annually)")


def parse_int(s):
    try:
        if s is None:
            return None
        s = s.strip()
        if not s:
            return None
        return int(float(s))
    except (ValueError, TypeError):
        return None


def _find_member(zf, suffix):
    """Find a file inside zip whose name ends with /SUFFIX (case-insens)."""
    suffix_up = suffix.upper()
    for n in zf.namelist():
        if n.upper().endswith(suffix_up):
            return n
    return None


# ═══════════════════════════════════════════════════════════════
# Data loading — download, extract, aggregate FARS by state-year
# ═══════════════════════════════════════════════════════════════

def aggregate_fars_year(zf, year):
    """Return dict[state_fips] -> {ped, cyc, total_fatal_persons,
    suv, suv_tight, light_truck, car, veh_total}."""
    out = defaultdict(lambda: {
        'ped': 0, 'cyc': 0, 'total_fatal_persons': 0,
        'suv': 0, 'suv_tight': 0, 'light_truck': 0, 'car': 0,
        'veh_total': 0,
    })

    person_member = _find_member(zf, "PERSON.CSV")
    if person_member is None:
        raise RuntimeError(f"PERSON.CSV not in FARS {year} zip")
    with zf.open(person_member) as f:
        text = io.TextIOWrapper(f, encoding='latin-1',
                                errors='replace', newline='')
        reader = csv.DictReader(text)
        for row in reader:
            st = parse_int(row.get('STATE'))
            pt = parse_int(row.get('PER_TYP'))
            inj = parse_int(row.get('INJ_SEV'))
            if st is None or st not in VALID_STATE_FIPS:
                continue
            if inj == FATAL_INJ_SEV:
                out[st]['total_fatal_persons'] += 1
                if pt == PED_PER_TYP:
                    out[st]['ped'] += 1
                elif pt == CYC_PER_TYP:
                    out[st]['cyc'] += 1

    vehicle_member = _find_member(zf, "VEHICLE.CSV")
    if vehicle_member is None:
        raise RuntimeError(f"VEHICLE.CSV not in FARS {year} zip")
    with zf.open(vehicle_member) as f:
        text = io.TextIOWrapper(f, encoding='latin-1',
                                errors='replace', newline='')
        reader = csv.DictReader(text)
        for row in reader:
            st = parse_int(row.get('STATE'))
            bt = parse_int(row.get('BODY_TYP'))
            if st is None or st not in VALID_STATE_FIPS:
                continue
            if bt is None:
                continue
            if bt in DENOM_BODY_TYPES:
                out[st]['veh_total'] += 1
                if bt in SUV_BODY_TYPES:
                    out[st]['suv'] += 1
                if bt in SUV_BODY_TYPES_TIGHT:
                    out[st]['suv_tight'] += 1
                if bt in LIGHT_TRUCK_BODY_TYPES:
                    out[st]['light_truck'] += 1
                if bt in PASSENGER_CAR_BODY_TYPES:
                    out[st]['car'] += 1
    return dict(out)


def load_data():
    os.makedirs(CACHE_DIR, exist_ok=True)
    manifest = load_manifest()
    state_year = {}  # (state_fips, year) -> counts dict
    for year in YEARS:
        url = FARS_URL_PATTERN.format(year=year)
        cache_path = os.path.join(CACHE_DIR, f"fars_{year}.zip")
        print(f"  [2.{year}] {url}")
        h = http_get_cached(url, cache_path, manifest)
        save_manifest(manifest)
        with zipfile.ZipFile(cache_path) as zf:
            agg = aggregate_fars_year(zf, year)
        for st, counts in agg.items():
            state_year[(st, year)] = counts
        ped_tot = sum(c['ped'] for c in agg.values())
        cyc_tot = sum(c['cyc'] for c in agg.values())
        veh_tot = sum(c['veh_total'] for c in agg.values())
        suv_tot = sum(c['suv'] for c in agg.values())
        print(f"    {year}: peds={ped_tot}  cyc={cyc_tot}  "
              f"veh={veh_tot}  suv_share={suv_tot/max(1,veh_tot):.3f}")
    return state_year


# ═══════════════════════════════════════════════════════════════
# Statistical helpers (stdlib)
# ═══════════════════════════════════════════════════════════════

def mean_val(xs):
    return sum(xs) / len(xs) if xs else 0.0


def two_way_within_transform(records, unit_key, time_key, y_key, x_keys):
    """Two-way within (state-FE + year-FE) demeaning.  Returns demeaned
    (y, X) lists."""
    unit_y = defaultdict(list)
    time_y = defaultdict(list)
    for r in records:
        unit_y[r[unit_key]].append(r[y_key])
        time_y[r[time_key]].append(r[y_key])
    unit_mean_y = {u: mean_val(vs) for u, vs in unit_y.items()}
    time_mean_y = {t: mean_val(vs) for t, vs in time_y.items()}
    grand_y = mean_val([r[y_key] for r in records])

    x_unit_means, x_time_means, x_grand = {}, {}, {}
    for k in x_keys:
        u_acc = defaultdict(list)
        t_acc = defaultdict(list)
        for r in records:
            u_acc[r[unit_key]].append(r[k])
            t_acc[r[time_key]].append(r[k])
        x_unit_means[k] = {u: mean_val(vs) for u, vs in u_acc.items()}
        x_time_means[k] = {t: mean_val(vs) for t, vs in t_acc.items()}
        x_grand[k] = mean_val([r[k] for r in records])

    y_t, x_t = [], []
    for r in records:
        u, tt = r[unit_key], r[time_key]
        y_t.append(r[y_key] - unit_mean_y[u] - time_mean_y[tt] + grand_y)
        x_t.append([r[k] - x_unit_means[k][u] - x_time_means[k][tt]
                    + x_grand[k] for k in x_keys])
    return y_t, x_t


def ols_multi(y, X):
    """OLS without intercept (assumes demeaned inputs)."""
    if not y:
        return [], None
    n = len(y)
    k = len(X[0]) if X[0] else 0
    if k == 0:
        return [], {'n': n, 'k': 0, 'ssr': 0.0, 'tss': 0.0, 'r2': 0.0}
    xtx = [[0.0] * k for _ in range(k)]
    xty = [0.0] * k
    for i in range(n):
        xi = X[i]
        yi = y[i]
        for a in range(k):
            xty[a] += xi[a] * yi
            for b in range(a, k):
                xtx[a][b] += xi[a] * xi[b]
    for a in range(k):
        for b in range(a):
            xtx[a][b] = xtx[b][a]
    A = [row[:] + [rhs] for row, rhs in zip(xtx, xty)]
    for i in range(k):
        piv = A[i][i]
        if abs(piv) < 1e-12:
            for j in range(i + 1, k):
                if abs(A[j][i]) > 1e-12:
                    A[i], A[j] = A[j], A[i]
                    piv = A[i][i]
                    break
        if abs(piv) < 1e-12:
            return [0.0] * k, None
        inv = 1.0 / piv
        for j in range(k + 1):
            A[i][j] *= inv
        for r_ in range(k):
            if r_ != i and abs(A[r_][i]) > 1e-12:
                f = A[r_][i]
                for j in range(k + 1):
                    A[r_][j] -= f * A[i][j]
    beta = [A[i][k] for i in range(k)]
    ssr = 0.0
    tss = 0.0
    ym = mean_val(y)
    for i in range(n):
        yhat = sum(X[i][a] * beta[a] for a in range(k))
        ssr += (y[i] - yhat) ** 2
        tss += (y[i] - ym) ** 2
    r2 = 1.0 - ssr / tss if tss > 0 else 0.0
    return beta, {'n': n, 'k': k, 'ssr': ssr, 'tss': tss, 'r2': r2}


def fit_panel_twoway(records, y_key, x_keys=('suv_share',)):
    y_t, x_t = two_way_within_transform(
        records, 'state', 'year', y_key, list(x_keys))
    beta, diag = ols_multi(y_t, x_t)
    return {k: b for k, b in zip(x_keys, beta)}, diag


# Module-level counter of singular OLS fits encountered during bootstrap
# and permutation loops.  Reset at the top of run_analysis().
SINGULAR_FITS = {'bootstrap': 0, 'permutation': 0}


def permutation_p_within_year(records, n_perms, rng, y_keys=('log_ped', 'log_cyc')):
    """Within-year permutation of suv_share across states. Returns
    p-values (two-sided) for beta_ped, beta_cyc, and
    beta_ped - beta_cyc (sibling contrast)."""
    obs_beta = {}
    for yk in y_keys:
        b, _ = fit_panel_twoway(records, yk, ('suv_share',))
        obs_beta[yk] = b['suv_share']
    obs_contrast = obs_beta[y_keys[0]] - obs_beta[y_keys[1]]

    by_year = defaultdict(list)
    for idx, r in enumerate(records):
        by_year[r['year']].append(idx)

    orig_suv = [r['suv_share'] for r in records]
    n_ge_ped = n_ge_cyc = n_ge_contrast = 0

    for _ in range(n_perms):
        suv_perm = orig_suv[:]
        for yr, idxs in by_year.items():
            vals = [orig_suv[i] for i in idxs]
            rng.shuffle(vals)
            for i, v in zip(idxs, vals):
                suv_perm[i] = v
        perm_recs = []
        for i, r in enumerate(records):
            rr = dict(r)
            rr['suv_share'] = suv_perm[i]
            perm_recs.append(rr)
        beta_p = {}
        for yk in y_keys:
            b, diag = fit_panel_twoway(perm_recs, yk, ('suv_share',))
            if diag is None:
                SINGULAR_FITS['permutation'] += 1
            beta_p[yk] = b['suv_share']
        if abs(beta_p[y_keys[0]]) >= abs(obs_beta[y_keys[0]]):
            n_ge_ped += 1
        if abs(beta_p[y_keys[1]]) >= abs(obs_beta[y_keys[1]]):
            n_ge_cyc += 1
        perm_contrast = beta_p[y_keys[0]] - beta_p[y_keys[1]]
        if abs(perm_contrast) >= abs(obs_contrast):
            n_ge_contrast += 1

    return {
        'observed': obs_beta,
        'observed_contrast': obs_contrast,
        'p_' + y_keys[0]: (n_ge_ped + 1) / (n_perms + 1),
        'p_' + y_keys[1]: (n_ge_cyc + 1) / (n_perms + 1),
        'p_sibling_contrast': (n_ge_contrast + 1) / (n_perms + 1),
        'n_permutations': n_perms,
    }


def bootstrap_ci_cluster(records, n_boot, rng,
                         y_keys=('log_ped', 'log_cyc'),
                         level=BOOTSTRAP_LEVEL):
    """State-cluster bootstrap. Resample with replacement at the state
    level; refit two-way FE panels for each outcome; return percentile
    CIs for beta_ped, beta_cyc, and beta_ped - beta_cyc."""
    by_state = defaultdict(list)
    for r in records:
        by_state[r['state']].append(r)
    states = list(by_state.keys())
    boots_ped = []
    boots_cyc = []
    boots_contrast = []
    for _ in range(n_boot):
        sample = []
        for _s in states:
            pick = states[rng.randrange(0, len(states))]
            for rec in by_state[pick]:
                sample.append(rec)
        b_p, dp = fit_panel_twoway(sample, y_keys[0], ('suv_share',))
        b_c, dc = fit_panel_twoway(sample, y_keys[1], ('suv_share',))
        if dp is None:
            SINGULAR_FITS['bootstrap'] += 1
        if dc is None:
            SINGULAR_FITS['bootstrap'] += 1
        boots_ped.append(b_p['suv_share'])
        boots_cyc.append(b_c['suv_share'])
        boots_contrast.append(b_p['suv_share'] - b_c['suv_share'])

    def _percentile(s, q):
        # Linear-interpolation percentile (Type-7, numpy default).
        # q is a probability in [0, 1]; s must be sorted.
        if not s:
            return 0.0
        if len(s) == 1:
            return s[0]
        h = q * (len(s) - 1)
        lo_i = int(math.floor(h))
        hi_i = int(math.ceil(h))
        if lo_i == hi_i:
            return s[lo_i]
        frac = h - lo_i
        return s[lo_i] * (1.0 - frac) + s[hi_i] * frac

    def _ci(b):
        s = sorted(b)
        lo = _percentile(s, (1 - level) / 2)
        hi = _percentile(s, (1 + level) / 2)
        return {'mean': mean_val(b), 'lo': lo, 'hi': hi}

    return {
        y_keys[0]: _ci(boots_ped),
        y_keys[1]: _ci(boots_cyc),
        'sibling_contrast': _ci(boots_contrast),
        'n_bootstrap': n_boot,
    }


# ═══════════════════════════════════════════════════════════════
# Sample assembly
# ═══════════════════════════════════════════════════════════════

def build_records(state_year_data,
                  exclude_states=None,
                  use_tight=False,
                  light_truck=False):
    """Return list of records with state, year, log_ped, log_cyc,
    suv_share. Drops state-years with too-few vehicles."""
    recs = []
    exclude_states = set(exclude_states or [])
    for (st, yr), c in state_year_data.items():
        if st in exclude_states:
            continue
        veh = c.get('veh_total', 0)
        if veh < MIN_VEHICLES_PER_STATEYEAR:
            continue
        if use_tight:
            suv = c.get('suv_tight', 0)
        elif light_truck:
            suv = c.get('light_truck', 0)
        else:
            suv = c.get('suv', 0)
        suv_share = suv / veh
        ped = c.get('ped', 0)
        cyc = c.get('cyc', 0)
        recs.append({
            'state': st, 'year': yr,
            'log_ped': math.log1p(ped),
            'log_cyc': math.log1p(cyc),
            'log_total': math.log1p(c.get('total_fatal_persons', 0)),
            'ped_per_total':
                ped / max(1, c.get('total_fatal_persons', 0)),
            'cyc_per_total':
                cyc / max(1, c.get('total_fatal_persons', 0)),
            'ped': ped,
            'cyc': cyc,
            'suv_share': suv_share,
            'veh_total': veh,
        })
    return recs


# ═══════════════════════════════════════════════════════════════
# Main analysis
# ═══════════════════════════════════════════════════════════════

def run_analysis(state_year_data):
    rng = random.Random(SEED)
    SINGULAR_FITS['bootstrap'] = 0
    SINGULAR_FITS['permutation'] = 0

    print("[4/8] Building main analysis sample...")
    records = build_records(state_year_data)
    n_state = len({r['state'] for r in records})
    n_year = len({r['year'] for r in records})
    print(f"  n_obs={len(records)}  n_states={n_state}  n_years={n_year}")
    if len(records) < 100 or n_state < 30 or n_year < 10:
        raise RuntimeError("Sample too small after filtering")

    print("[5/8] Main panel FE regressions (state FE + year FE)...")
    beta_ped, diag_ped = fit_panel_twoway(records, 'log_ped')
    beta_cyc, diag_cyc = fit_panel_twoway(records, 'log_cyc')
    contrast = beta_ped['suv_share'] - beta_cyc['suv_share']
    print(f"  beta_ped   = {beta_ped['suv_share']:.4f}  "
          f"(R^2_within={diag_ped['r2']:.3f})")
    print(f"  beta_cyc   = {beta_cyc['suv_share']:.4f}  "
          f"(R^2_within={diag_cyc['r2']:.3f})")
    print(f"  sibling contrast (ped - cyc) = {contrast:.4f}")

    print(f"[6/8] Within-year permutation null (N={N_PERMUTATIONS})...")
    perm = permutation_p_within_year(records, N_PERMUTATIONS, rng)
    print(f"  p(beta_ped)       = {perm['p_log_ped']:.4f}")
    print(f"  p(beta_cyc)       = {perm['p_log_cyc']:.4f}")
    print(f"  p(sibling contrast) = {perm['p_sibling_contrast']:.4f}")

    print(f"[7/8] State-cluster bootstrap CIs (N={N_BOOTSTRAP})...")
    ci = bootstrap_ci_cluster(records, N_BOOTSTRAP, rng)
    for k in ('log_ped', 'log_cyc', 'sibling_contrast'):
        v = ci[k]
        print(f"  {k:>18s}: mean={v['mean']:.4f}  "
              f"95% CI [{v['lo']:.4f}, {v['hi']:.4f}]")

    print("[8/8] Sensitivity analyses...")
    sens = {}

    # (a) Drop CA (6) + TX (48) — the two largest-population states
    rec_notx = build_records(
        state_year_data, exclude_states=DROP_STATES_FOR_SENSITIVITY)
    bp, _ = fit_panel_twoway(rec_notx, 'log_ped')
    bc, _ = fit_panel_twoway(rec_notx, 'log_cyc')
    sens['drop_ca_tx'] = {
        'dropped_state_fips': sorted(DROP_STATES_FOR_SENSITIVITY),
        'n_obs': len(rec_notx),
        'beta_ped': bp['suv_share'],
        'beta_cyc': bc['suv_share'],
        'sibling_contrast': bp['suv_share'] - bc['suv_share'],
    }
    print(f"  (a) drop CA/TX: contrast={sens['drop_ca_tx']['sibling_contrast']:.4f}")

    # (b) Tight SUV coding (14-16)
    rec_tight = build_records(state_year_data, use_tight=True)
    bp, _ = fit_panel_twoway(rec_tight, 'log_ped')
    bc, _ = fit_panel_twoway(rec_tight, 'log_cyc')
    sens['suv_tight'] = {
        'n_obs': len(rec_tight),
        'body_type_range': '14..16',
        'beta_ped': bp['suv_share'],
        'beta_cyc': bc['suv_share'],
        'sibling_contrast': bp['suv_share'] - bc['suv_share'],
    }
    print(f"  (b) SUV tight (14-16): contrast={sens['suv_tight']['sibling_contrast']:.4f}")

    # (c) Broad light-truck coding (14-49)
    rec_lt = build_records(state_year_data, light_truck=True)
    bp, _ = fit_panel_twoway(rec_lt, 'log_ped')
    bc, _ = fit_panel_twoway(rec_lt, 'log_cyc')
    sens['light_truck_broad'] = {
        'n_obs': len(rec_lt),
        'body_type_range': '14..49',
        'beta_ped': bp['suv_share'],
        'beta_cyc': bc['suv_share'],
        'sibling_contrast': bp['suv_share'] - bc['suv_share'],
    }
    print(f"  (c) Light truck broad (14-49): contrast={sens['light_truck_broad']['sibling_contrast']:.4f}")

    # (d) Pre-SPLIT_YEAR vs post-SPLIT_YEAR subsamples
    windows = [
        ('pre' + str(SPLIT_YEAR), min(YEARS), SPLIT_YEAR - 1),
        ('post' + str(SPLIT_YEAR), SPLIT_YEAR, max(YEARS)),
    ]
    for label, lo_y, hi_y in windows:
        sub = [r for r in records if lo_y <= r['year'] <= hi_y]
        if len({r['state'] for r in sub}) < MIN_STATES_FOR_SUBSAMPLE:
            continue
        bp, _ = fit_panel_twoway(sub, 'log_ped')
        bc, _ = fit_panel_twoway(sub, 'log_cyc')
        sens[f'window_{label}'] = {
            'year_range': [lo_y, hi_y],
            'n_obs': len(sub),
            'beta_ped': bp['suv_share'],
            'beta_cyc': bc['suv_share'],
            'sibling_contrast': bp['suv_share'] - bc['suv_share'],
        }
        print(f"  (d) {label}: contrast={sens[f'window_{label}']['sibling_contrast']:.4f}")

    # (e) Alternate outcome: share of fatalities that are peds / cyc
    bp, _ = fit_panel_twoway(records, 'ped_per_total')
    bc, _ = fit_panel_twoway(records, 'cyc_per_total')
    sens['outcome_share_of_total'] = {
        'outcome_ped': 'pedestrian fatalities / total persons killed',
        'outcome_cyc': 'cyclist fatalities / total persons killed',
        'beta_ped': bp['suv_share'],
        'beta_cyc': bc['suv_share'],
        'sibling_contrast': bp['suv_share'] - bc['suv_share'],
    }
    print(f"  (e) share-of-total outcome: contrast={sens['outcome_share_of_total']['sibling_contrast']:.4f}")

    # (f) Placebo: randomize SUV share completely (mechanical null)
    rng_pl = random.Random(SEED + 1)
    rec_pl = []
    suvs = [r['suv_share'] for r in records]
    rng_pl.shuffle(suvs)
    for r, s in zip(records, suvs):
        rr = dict(r)
        rr['suv_share'] = s
        rec_pl.append(rr)
    bp, _ = fit_panel_twoway(rec_pl, 'log_ped')
    bc, _ = fit_panel_twoway(rec_pl, 'log_cyc')
    sens['placebo_full_shuffle'] = {
        'description': 'SUV share fully shuffled across all obs (breaks state and year structure)',
        'beta_ped': bp['suv_share'],
        'beta_cyc': bc['suv_share'],
        'sibling_contrast': bp['suv_share'] - bc['suv_share'],
    }
    print(f"  (f) full-shuffle placebo: contrast={sens['placebo_full_shuffle']['sibling_contrast']:.4f}")

    # Summary statistics of panel
    ped_tot = sum(r['ped'] for r in records)
    cyc_tot = sum(r['cyc'] for r in records)
    mean_suv = mean_val([r['suv_share'] for r in records])
    suv_by_year = defaultdict(list)
    for r in records:
        suv_by_year[r['year']].append(r['suv_share'])
    suv_trajectory = {y: mean_val(suv_by_year[y]) for y in sorted(suv_by_year)}
    ped_by_year = defaultdict(int)
    cyc_by_year = defaultdict(int)
    for r in records:
        ped_by_year[r['year']] += r['ped']
        cyc_by_year[r['year']] += r['cyc']
    year_totals = {
        y: {'ped': ped_by_year[y], 'cyc': cyc_by_year[y]}
        for y in sorted(ped_by_year)
    }

    results = {
        'config': {
            'seed': SEED,
            'years': YEARS,
            'fars_url_pattern': FARS_URL_PATTERN,
            'n_permutations': N_PERMUTATIONS,
            'n_bootstrap': N_BOOTSTRAP,
            'suv_body_types': sorted(SUV_BODY_TYPES),
            'suv_body_types_tight': sorted(SUV_BODY_TYPES_TIGHT),
            'light_truck_body_types': sorted(LIGHT_TRUCK_BODY_TYPES),
            'denom_body_types': sorted(DENOM_BODY_TYPES),
            'ped_per_typ': PED_PER_TYP,
            'cyc_per_typ': CYC_PER_TYP,
            'fatal_inj_sev': FATAL_INJ_SEV,
            'fe_specification':
                'state FE + year FE via 2-way within transform; '
                'single regressor suv_share',
            'bootstrap_level': BOOTSTRAP_LEVEL,
            'split_year': SPLIT_YEAR,
            'drop_states_for_sensitivity':
                sorted(DROP_STATES_FOR_SENSITIVITY),
            'min_vehicles_per_stateyear': MIN_VEHICLES_PER_STATEYEAR,
            'http_retries': HTTP_RETRIES,
            'http_timeout_secs': HTTP_TIMEOUT_SECS,
            'placebo_margin_min': PLACEBO_MARGIN_MIN,
            'verify_max_abs_beta': VERIFY_MAX_ABS_BETA,
            'verify_min_ci_width_frac': VERIFY_MIN_CI_WIDTH_FRAC,
            'expected_sha256_partial': EXPECTED_SHA256,
        },
        'sample': {
            'n_obs': len(records),
            'n_states': n_state,
            'n_years': n_year,
            'min_year': min(r['year'] for r in records),
            'max_year': max(r['year'] for r in records),
            'total_pedestrian_fatalities': ped_tot,
            'total_cyclist_fatalities': cyc_tot,
            'mean_suv_share': mean_suv,
            'suv_share_by_year': suv_trajectory,
            'national_fatalities_by_year': year_totals,
        },
        'main': {
            'beta_ped': beta_ped['suv_share'],
            'beta_cyc': beta_cyc['suv_share'],
            'sibling_contrast': contrast,
            'r2_ped_within': diag_ped['r2'],
            'r2_cyc_within': diag_cyc['r2'],
            'permutation': perm,
            'bootstrap_ci': ci,
        },
        'sensitivity': sens,
        'diagnostics': {
            'singular_fits_bootstrap': SINGULAR_FITS['bootstrap'],
            'singular_fits_permutation': SINGULAR_FITS['permutation'],
        },
    }
    return results


# ═══════════════════════════════════════════════════════════════
# Report generation
# ═══════════════════════════════════════════════════════════════

def generate_report(results):
    with open(RESULTS_JSON, 'w') as f:
        json.dump(results, f, indent=2, default=str)
    m = results['main']
    s = results['sample']
    ci = m['bootstrap_ci']
    p = m['permutation']
    with open(REPORT_MD, 'w') as f:
        f.write("# SUV Fleet Share and the Post-2009 U.S. Pedestrian-Fatality Surge — Report\n\n")
        f.write("## Sample\n\n")
        f.write(f"- State-year observations: {s['n_obs']}\n")
        f.write(f"- States: {s['n_states']}\n")
        f.write(f"- Years: {s['n_years']} ({s['min_year']}-{s['max_year']})\n")
        f.write(f"- Total pedestrian fatalities: {s['total_pedestrian_fatalities']}\n")
        f.write(f"- Total cyclist fatalities: {s['total_cyclist_fatalities']}\n")
        f.write(f"- Mean SUV share of fatal-crash fleet: {s['mean_suv_share']:.3f}\n\n")
        f.write("## Main panel FE regressions (state FE + year FE)\n\n")
        f.write("| Outcome | beta(suv_share) | Bootstrap 95% CI | Permutation p |\n")
        f.write("|---|---|---|---|\n")
        f.write(f"| log(pedestrian fatalities + 1) | {m['beta_ped']:.4f} "
                f"| [{ci['log_ped']['lo']:.4f}, {ci['log_ped']['hi']:.4f}] "
                f"| {p['p_log_ped']:.4f} |\n")
        f.write(f"| log(cyclist fatalities + 1) | {m['beta_cyc']:.4f} "
                f"| [{ci['log_cyc']['lo']:.4f}, {ci['log_cyc']['hi']:.4f}] "
                f"| {p['p_log_cyc']:.4f} |\n")
        f.write(f"| Sibling contrast (ped - cyc) | {m['sibling_contrast']:.4f} "
                f"| [{ci['sibling_contrast']['lo']:.4f}, {ci['sibling_contrast']['hi']:.4f}] "
                f"| {p['p_sibling_contrast']:.4f} |\n\n")
        f.write("## Sensitivity (sibling contrast: ped - cyc)\n\n")
        for k, v in results['sensitivity'].items():
            if 'sibling_contrast' in v:
                f.write(f"- {k}: {v['sibling_contrast']:.4f}\n")
    print("  results.json, report.md written")


# ═══════════════════════════════════════════════════════════════
# Main entrypoint
# ═══════════════════════════════════════════════════════════════

def main():
    if '--verify' in sys.argv:
        return verify()

    global N_PERMUTATIONS, N_BOOTSTRAP
    if '--quick' in sys.argv:
        # Smoke-test mode: reduces permutation + bootstrap counts so an
        # agent can validate executability end-to-end in ~2 min.  NOT
        # intended for publication-quality inference.
        N_PERMUTATIONS = 200
        N_BOOTSTRAP = 200
        print(f"[--quick] N_PERMUTATIONS={N_PERMUTATIONS}  "
              f"N_BOOTSTRAP={N_BOOTSTRAP}  (smoke-test mode)")

    t0 = time.time()
    try:
        print("[1/8] Workspace prep...")
        os.makedirs(CACHE_DIR, exist_ok=True)

        print(f"[2/8] Downloading + aggregating FARS for {len(YEARS)} years...")
        try:
            state_year_data = load_data()
        except (urllib.error.URLError, RuntimeError, OSError) as e:
            print(f"ERROR: FARS data unavailable: {e}", file=sys.stderr)
            print(f"       URL pattern: {FARS_URL_PATTERN}", file=sys.stderr)
            print("       Check internet connectivity or cache integrity.",
                  file=sys.stderr)
            sys.exit(2)
        print(f"[3/8] Assembled state-year counts: {len(state_year_data)} rows")

        results = run_analysis(state_year_data)
        generate_report(results)
    except Exception as e:
        print(f"FATAL: analysis failed: {type(e).__name__}: {e}",
              file=sys.stderr)
        sys.exit(1)

    elapsed = time.time() - t0
    print(f"\nRuntime: {elapsed:.0f}s")
    print("ANALYSIS COMPLETE")


# ═══════════════════════════════════════════════════════════════
# Verification
# ═══════════════════════════════════════════════════════════════

def verify():
    print("Running verification...\n")
    ok = fail = 0

    def chk(name, cond):
        nonlocal ok, fail
        status = "PASS" if cond else "FAIL"
        print(f"  {status}: {name}")
        if cond:
            ok += 1
        else:
            fail += 1

    if not os.path.exists(RESULTS_JSON):
        print("FAIL: results.json not found")
        sys.exit(1)
    with open(RESULTS_JSON) as f:
        r = json.load(f)

    chk("1. results.json has core keys",
        all(k in r for k in ('config', 'sample', 'main', 'sensitivity')))
    chk("2. report.md exists and non-empty",
        os.path.exists(REPORT_MD) and os.path.getsize(REPORT_MD) > 300)

    s = r['sample']
    chk("3. >= 500 state-year observations",
        s.get('n_obs', 0) >= 500)
    chk("4. >= 40 states covered",
        s.get('n_states', 0) >= 40)
    chk("5. >= 15 years covered",
        s.get('n_years', 0) >= 15)
    chk("6. Max year >= 2022",
        s.get('max_year', 0) >= 2022)
    chk("7. Min year <= 2005",
        s.get('min_year', 99999) <= 2005)
    chk("8. Total pedestrian fatalities > 50,000",
        s.get('total_pedestrian_fatalities', 0) > 50000)
    chk("9. Total cyclist fatalities > 5,000",
        s.get('total_cyclist_fatalities', 0) > 5000)
    chk("10. Mean SUV share in a sensible range (0.05..0.9)",
        0.05 < s.get('mean_suv_share', -1) < 0.9)

    m = r['main']
    chk("11. beta_ped is finite and |beta_ped| < 20 (sanity)",
        isinstance(m.get('beta_ped'), (int, float))
        and math.isfinite(m.get('beta_ped', float('nan')))
        and abs(m.get('beta_ped', 99)) < 20.0)
    chk("12. beta_cyc is finite and |beta_cyc| < 20 (sanity)",
        isinstance(m.get('beta_cyc'), (int, float))
        and math.isfinite(m.get('beta_cyc', float('nan')))
        and abs(m.get('beta_cyc', 99)) < 20.0)
    chk("13. sibling contrast is finite",
        isinstance(m.get('sibling_contrast'), (int, float))
        and math.isfinite(m.get('sibling_contrast', float('nan'))))
    chk("14. R^2_within for pedestrian model in [0,1]",
        0.0 <= m.get('r2_ped_within', -1) <= 1.0)
    chk("15. R^2_within for cyclist model in [0,1]",
        0.0 <= m.get('r2_cyc_within', -1) <= 1.0)

    perm = m['permutation']
    chk("16. Permutation p(ped) in [0,1]",
        0.0 <= perm.get('p_log_ped', -1) <= 1.0)
    chk("17. Permutation p(cyc) in [0,1]",
        0.0 <= perm.get('p_log_cyc', -1) <= 1.0)
    chk("18. Permutation p(sibling contrast) in [0,1]",
        0.0 <= perm.get('p_sibling_contrast', -1) <= 1.0)
    chk("19. N_permutations >= 1000",
        perm.get('n_permutations', 0) >= 1000)

    ci = m['bootstrap_ci']
    chk("20. Bootstrap CI for ped has nonzero width",
        'log_ped' in ci and ci['log_ped']['hi'] > ci['log_ped']['lo'])
    chk("21. Bootstrap CI for cyc has nonzero width",
        'log_cyc' in ci and ci['log_cyc']['hi'] > ci['log_cyc']['lo'])
    chk("22. Bootstrap CI for contrast has nonzero width",
        'sibling_contrast' in ci
        and ci['sibling_contrast']['hi'] > ci['sibling_contrast']['lo'])
    chk("23. Bootstrap N >= 1000",
        ci.get('n_bootstrap', 0) >= 1000)

    sens = r['sensitivity']
    chk("24. Sensitivity has drop_ca_tx",
        'drop_ca_tx' in sens and 'sibling_contrast' in sens['drop_ca_tx'])
    chk("25. Sensitivity has tight SUV coding",
        'suv_tight' in sens)
    chk("26. Sensitivity has broad light-truck coding",
        'light_truck_broad' in sens)
    chk("27. Sensitivity has window_pre2014 and window_post2014",
        'window_pre2014' in sens and 'window_post2014' in sens)
    chk("28. Sensitivity has alternate outcome (share-of-total)",
        'outcome_share_of_total' in sens)
    obs_contrast = abs(m.get('sibling_contrast', 0.0))
    placebo_contrast = abs(sens.get('placebo_full_shuffle', {})
                           .get('sibling_contrast', 99))
    chk("29. Full-shuffle placebo |contrast| < |observed contrast| (null smaller than signal)",
        placebo_contrast < obs_contrast)

    cfg = r['config']
    chk("30. Seed recorded as 42",
        cfg.get('seed') == 42)
    chk("31. Sample covers every configured year (no year-wide parse failures)",
        len(cfg.get('years', [])) == s.get('n_years'))
    chk("32. FE specification documented",
        'fe_specification' in cfg and 'state' in cfg['fe_specification'].lower())
    chk("33. Bootstrap level recorded as 0.95",
        abs(cfg.get('bootstrap_level', 0) - 0.95) < 1e-9)

    # Manifest existence
    chk("34. Data manifest file exists",
        os.path.exists(MANIFEST_FILE))
    diag = r.get('diagnostics', {})
    chk("35. No singular bootstrap refits",
        diag.get('singular_fits_bootstrap', 99) == 0)
    chk("36. No singular permutation refits",
        diag.get('singular_fits_permutation', 99) == 0)

    # Additional rigor checks: effect-size plausibility bounds,
    # CI-width sanity, sensitivity robustness, negative / falsification
    # (placebo), SHA256 coverage.
    chk("37. |beta_ped| within plausibility cap",
        abs(m.get('beta_ped', 99)) < VERIFY_MAX_ABS_BETA)
    chk("38. |beta_cyc| within plausibility cap",
        abs(m.get('beta_cyc', 99)) < VERIFY_MAX_ABS_BETA)
    _contrast = abs(m.get('sibling_contrast', 0.0))
    _w = ci['sibling_contrast']['hi'] - ci['sibling_contrast']['lo']
    chk("39. Sibling-contrast CI width >= 1% of |estimate| (bootstrap not collapsed)",
        _contrast == 0.0 or _w >= VERIFY_MIN_CI_WIDTH_FRAC * _contrast)
    chk("40. Pedestrian R^2_within dominates cyclist R^2_within "
        "(SUV share explains more pedestrian than cyclist variation)",
        m.get('r2_ped_within', 0) > m.get('r2_cyc_within', 0))
    # Sensitivity robustness: drop-CA/TX, tight SUV, and broad
    # light-truck sensitivities should all return the same SIGN on
    # the sibling contrast as the main analysis.  This is the
    # key robustness falsification — if removing the two largest
    # states or retightening the coding flips the sign, the main
    # result is not robust.
    _main_sign = (1 if m.get('sibling_contrast', 0) > 0 else
                  (-1 if m.get('sibling_contrast', 0) < 0 else 0))
    _sens_signs_ok = all(
        (_main_sign == 0) or
        (1 if sens[k]['sibling_contrast'] > 0 else
         (-1 if sens[k]['sibling_contrast'] < 0 else 0)) == _main_sign
        for k in ('drop_ca_tx', 'suv_tight', 'light_truck_broad')
        if k in sens and 'sibling_contrast' in sens[k]
    )
    chk("41. Sign of sibling contrast survives drop-CA/TX, tight-SUV, "
        "and broad-light-truck sensitivities",
        _sens_signs_ok)
    # Falsification: full-shuffle placebo is expected to be small,
    # and the observed |contrast| must exceed it by a multiplicative
    # margin (PLACEBO_MARGIN_MIN).
    _placebo = abs(sens.get('placebo_full_shuffle', {}).get('sibling_contrast', 99))
    chk(f"42. Observed |contrast| >= {PLACEBO_MARGIN_MIN}x "
        f"full-shuffle-placebo |contrast| (falsification control)",
        _contrast >= PLACEBO_MARGIN_MIN * _placebo)
    # SHA256 coverage for reproducibility anchor.  JSON serialises
    # dict int keys as strings, so coerce both sides for comparison.
    _expected_keys = {str(k) for k in cfg.get('expected_sha256_partial', {})}
    chk("43. Every configured year has an expected SHA256 recorded",
        all(str(y) in _expected_keys for y in cfg.get('years', [])))
    # Permutation alpha sanity: the p-value resolution has at least
    # 1/(N+1) granularity — never below 1/(N+1).
    _min_p = 1.0 / (perm.get('n_permutations', 2000) + 1)
    chk("44. All permutation p-values respect the 1/(N+1) resolution floor",
        perm.get('p_log_ped', 0) >= _min_p
        and perm.get('p_log_cyc', 0) >= _min_p
        and perm.get('p_sibling_contrast', 0) >= _min_p)

    # Additional sensitivity / parameterization / reproducibility checks.
    chk("45. SPLIT_YEAR recorded in config and within year range",
        isinstance(cfg.get('split_year'), int)
        and min(cfg.get('years', [0])) < cfg['split_year']
        <= max(cfg.get('years', [0])))
    chk("46. DROP_STATES_FOR_SENSITIVITY recorded and non-empty",
        len(cfg.get('drop_states_for_sensitivity', [])) > 0)
    chk("47. HTTP retry / timeout parameters recorded",
        isinstance(cfg.get('http_retries'), int)
        and cfg.get('http_retries') >= 1
        and isinstance(cfg.get('http_timeout_secs'), int)
        and cfg.get('http_timeout_secs') >= 30)
    # Early vs late window sensitivity should not both flip from
    # main-analysis sign — requires at least one of the two halves
    # to match the main-analysis sign.
    _pre_k, _post_k = f'window_pre{cfg.get("split_year", 2014)}', f'window_post{cfg.get("split_year", 2014)}'
    _pre_sc = sens.get(_pre_k, {}).get('sibling_contrast')
    _post_sc = sens.get(_post_k, {}).get('sibling_contrast')
    chk("48. At least one window subsample preserves main-analysis sign",
        _main_sign == 0 or
        (_pre_sc is not None and _pre_sc * _main_sign > 0) or
        (_post_sc is not None and _post_sc * _main_sign > 0))
    chk("49. Placebo margin config matches verify-threshold",
        abs(cfg.get('placebo_margin_min', 0) - PLACEBO_MARGIN_MIN) < 1e-9)

    print(f"\n{ok}/{ok + fail} checks passed")
    if fail:
        print("VERIFICATION FAILED")
        sys.exit(1)
    else:
        print("ALL CHECKS PASSED")
        sys.exit(0)


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

**Expected output:** File `analyze.py` written, exit code 0.

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_suv-fleet-share-and-pedestrian-fatality-surge && python3 analyze.py
```

**Expected output:**
- Prints `[1/8]` through `[8/8]` progress sections.
- Downloads 18 FARS National CSV ZIPs (2005–2022) to
  `cache/fars_{YEAR}.zip`, each with SHA256 recorded in
  `data_manifest.json`.
- Parses each year's `PERSON.CSV` (pedestrian / cyclist fatality
  counts) and `VEHICLE.CSV` (SUV share of fatal-crash fleet).
- Builds a state-year panel of ≥900 observations across 50+ states
  and 18 years.
- Fits state-FE + year-FE panels for `log(pedestrian_fat + 1)` and
  `log(cyclist_fat + 1)` on `suv_share`.
- Reports the sibling-contrast β_ped − β_cyc.
- Runs 2,000 within-year label permutations and 2,000
  state-clustered bootstrap replicates.
- Runs six sensitivity analyses (drop CA/TX, tight vs broad body
  codings, pre-2014 vs post-2014 subsamples, share-of-total
  fatalities outcome, full-shuffle placebo).
- Writes `results.json` and `report.md`.
- Final line: `ANALYSIS COMPLETE`.
- Runtime: 12–25 minutes first run, 4–8 minutes on rerun (cache hit).
- Exit code 0.

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_suv-fleet-share-and-pedestrian-fatality-surge && python3 analyze.py --verify
```

**Expected output:**
- 49/49 checks passed
- `ALL CHECKS PASSED`
- Exit code 0

## Expected Outputs

| File | Description |
|---|---|
| `results.json` | Config, sample counts, main panel coefficients with bootstrap CIs, permutation p-values, full sensitivity table |
| `report.md` | Human-readable Markdown report |
| `cache/fars_{YEAR}.zip` | 18 raw FARS National CSV ZIPs (2005–2022), SHA256-recorded |
| `data_manifest.json` | SHA256 hashes of every cached file |

## Success Criteria

1. Script exits 0 on both normal run and `--verify`.
2. ≥500 state-year observations spanning ≥40 states and ≥15 years,
   covering 2005 through 2022 inclusive.
3. ≥2,000 within-year permutation iterations and ≥2,000
   state-clustered bootstrap replicates for each outcome.
4. Full-shuffle placebo sensitivity returns a finite,
   bounded-magnitude contrast.
5. All 49 verification assertions pass (including 13 additional
   rigor checks covering effect-size plausibility bounds,
   CI-width sanity, sensitivity sign-robustness, full-shuffle
   placebo / falsification margin, SHA256 coverage, permutation
   p-value resolution floor, configuration completeness for
   SPLIT_YEAR / DROP_STATES / HTTP retry / placebo margin, and
   window-subsample sign preservation).
6. Expected SHA256 hashes for every downloaded FARS year ZIP are
   recorded in `data_manifest.json`. Drift is logged but not fatal
   unless `STRICT_SHA256=True` is set in the DOMAIN CONFIGURATION
   block.

## Failure Conditions

1. Import errors (script uses only Python 3.8+ standard library).
2. Any FARS year ZIP unreachable after 4 retries.
3. Fewer than 500 state-year observations or fewer than 40 states
   after filtering.
4. Any `--verify` assertion fails.
5. `results.json` or `report.md` not created.