---
name: "FAR 117 Fatigue DiD with Cargo and Commuter Placebos"
description: "Evaluates whether the January 2014 FAR 117 flight/duty/rest rule reduced fatigue-related ASRS reports in Part 121 passenger operations relative to exempt Part 121F cargo and Part 135 commuter controls, using difference-in-differences across four outcome transformations, a 5,000-permutation placebo test on random treatment dates, a passenger/cargo share-robustness sweep, and a power analysis under a −5% log-rate alternative."
version: "2.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "aviation-safety", "far-117", "difference-in-differences", "placebo-test", "asrs", "policy-evaluation"]
python_version: ">=3.8"
dependencies: []
data_source: "NASA Aviation Safety Reporting System (ASRS), https://asrs.arc.nasa.gov"
data_revision: "ASRS Program Briefing + Air Carrier Flight Crew Fatigue Report Set + Commuter Fatigue Report Set; cumulative counts through August/September 2024; provenance pinned by SHA256."
---

# FAR 117 Fatigue DiD with Cargo and Commuter Placebos

> **Use this skill when:** you need to test whether a dated policy change (here: January 4, 2014 FAR 117 flight/duty/rest rule) produced a detectable change in a voluntary-reporting time series (here: fatigue-coded ASRS reports for Part 121 passenger carriers), using a carve-out/exempt group as a negative control, with placebo-date permutation, bootstrap CIs, share-robustness, and a power calibration. The statistical machinery is domain-agnostic; only the `DOMAIN CONFIGURATION` block at the top of the script is aviation-specific.

## Research Question

Did the January 4, 2014 Federal Aviation Regulation Part 117 flight/duty/rest rule reduce the monthly rate of fatigue-coded NASA Aviation Safety Reporting System (ASRS) incident reports filed by Part 121 passenger crews, **net of** secular trends present in the exempt Part 121F all-cargo segment and the exempt Part 135 commuter segment?

We answer this via a two-control difference-in-differences design with (i) four low-count outcome transformations, (ii) five passenger/cargo-share specifications plus a time-varying share, (iii) a 5,000-permutation placebo-date null, and (iv) a power analysis under a realistic –5% log-rate alternative DGP.

## When to Use This Skill

Use this skill when you need to evaluate whether the January 4, 2014 Federal Aviation Regulation Part 117 (flight/duty/rest) rule reduced fatigue-related Aviation Safety Reporting System (ASRS) reports in Part 121 passenger carriers, while guarding against confounding from concurrent industry changes by comparing against the exempt Part 121F all-cargo segment **and** exempt Part 135 commuter segment as placebo (never-treated) groups.

## Prerequisites

- **Python**: 3.8 or later. **No pip install.** Standard library only (no numpy, pandas, scipy).
- **Operating system**: Linux / macOS / WSL. Windows paths should substitute `/tmp` with an equivalent writable dir.
- **Network access**: HTTPS egress to `asrs.arc.nasa.gov` for the first run to fetch three PDFs (≈1.3 MB total). Subsequent runs use the on-disk cache and require no network.
- **Disk space**: ≈5 MB for the workspace, PDFs, and artifacts.
- **Expected runtime**: ≈2–5 minutes on one modern CPU core (dominated by 5,000 placebo permutations + bootstrap + multi-share sensitivity).
- **Environment variables**: none required.
- **Write access**: workspace directory only (defaults to `/tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did/`).

### Preconditions (summary)
- Python 3.8+ standard library only; no third-party installs.
- Internet egress on first run to download three NASA ASRS PDF publications (≈1.3 MB total). Subsequent runs use the local cache.
- Runtime ≈2–5 minutes end-to-end on a single modern CPU.
- No write access required outside the workspace directory.

## Adaptation Guidance

This skill implements a **generic difference-in-differences (DiD) evaluator**. All domain-specific values live in the `DOMAIN CONFIGURATION` block at the top of `script.py`. The statistical machinery below that block is reusable across policies and datasets.

### What to change (domain-specific constants)

| Constant | Purpose | Example alternate value |
|---|---|---|
| `DATA_URLS` | Map of logical name → remote URL of input PDFs/data | swap for your own stable URLs |
| `EXPECTED_SHA256` | Pinned integrity hash per URL | recompute on first download |
| `TREATMENT_DATE` | ISO date when the policy takes effect | `"2021-01-01"` for a 2021 rule |
| `PRE_WINDOW_MONTHS` / `POST_WINDOW_MONTHS` | Window half-widths in months | 36 / 36 for a 6-year window |
| `TREATED_LABEL` / `CONTROL_LABEL` / `ALT_CONTROL_LABEL` | Display names for treated + two controls | domain-specific group names |
| `ANNUAL_INTAKE` | Year → total intake anchor for panel reconstruction | dict of year → count |
| `TREATED_SHARE` / `CONTROL_SHARE` / `ALT_CONTROL_SHARE` | Baseline operator-type shares of intake | domain-specific splits |
| `SHARE_SENSITIVITY` | List of (treated, control) split pairs to sweep | `[(0.6,0.4), ...]` |
| `PASS_SHARE_START` / `PASS_SHARE_END` | Endpoints of linear share drift for time-varying run | 0.9 → 0.7 |
| `SEASONAL_SHAPE` | 12-element monthly scaling factor (sums to 12) | sine wave, flat, etc. |
| `TRUE_LOG_EFFECT_ALT` | Effect size under alt DGP for power analysis | −0.10 for a stronger prior |
| `FATIGUE_KEYWORDS` | Narrative-flag keyword list for the PDF corpus | e.g. `["emission", "nox", ...]` for pollution |

### What to keep unchanged (general methodology)

The following functions are domain-agnostic and operate on an abstract panel of `{group, time, y, D, treated_group}` records:

- `two_way_fe_estimator(panel)` — Two-way FE DiD point estimate
- `cluster_robust_se(panel, beta)` — Cluster-robust SE at the time level
- `bootstrap_ci(panel, n_boot, block_len, rng)` — Block bootstrap CI
- `permutation_placebo_test(panel, n_perm, rng)` — Placebo-date permutation null
- `run_did_with_outcome(panel, outcome_key)` — DiD on a chosen outcome transformation
- Outcome transformations: `y`, `log1p(y)`, `asinh(y)`, `log(y + 0.5)`

### How to swap in a new dataset

1. **Inputs**: update `DATA_URLS` and `EXPECTED_SHA256` to point to your publication(s). The PDF text extractor is generic but you can swap `load_data()` for a CSV / JSON loader. Returned structure must be a list of `{group, time, y, y_log1p, y_asinh, y_log_half, D, treated_group}` dicts.
2. **Panel reconstruction**: if you have per-record data, replace `build_panel()` with a direct aggregator (groupby month + group). Keep the same output schema.
3. **Statistical run**: nothing in `run_analysis()` requires modification — it operates on the panel schema.
4. **Verification**: update plausibility bounds in `verify()` if your effect-size units differ (e.g., tighten `abs(beta) <= 2.0` for low-variance outcomes).

---

## Overview

Federal Aviation Regulation Part 117 — promulgated by the FAA in response to the Colgan Air 3407 accident and effective January 4, 2014 — tightened flight-time, duty-period, and rest requirements for Part 121 air-carrier pilots. Part 121F all-cargo carriers were explicitly **carved out** of the rule (the "cargo carve-out"), and Part 135 commuter operators were separately exempt, producing two natural comparison groups. A common narrative claims FAR 117 improved aviation safety. This skill tests that claim with a difference-in-differences (DiD) framework against **both** placebo control groups, across **four** outcome transformations (raw counts, log(1+y), asinh, log(y+0.5)), under **five** passenger/cargo share specifications, with a **5,000-iteration placebo-date permutation test**.

**Methodological hook.** Existing FAR 117 evaluations report linear DiD coefficients with clustered standard errors but do not (a) cross-check against a second placebo control, (b) systematically vary the low-count transformation to probe Jensen's-inequality bias in log(1+y), or (c) permute the treatment date to confirm that secular trend alone does not flag January 2014 as special. We add all three and a calibrated power analysis under a realistic −5% log-rate alternative.

**Data.** The skill downloads three NASA ASRS publications (Program Briefing PDF, Air Carrier FAR 121 Flight Crew Fatigue Report Set, Commuter Fatigue Report Set) from stable nasa.gov URLs, pins their contents by SHA256, and **reconstructs** a monthly Part 121 Passenger vs. Part 121F Cargo vs. Part 135 Commuter panel from published annual aggregates (ASRS yearly intake, fatigue-category share, seasonal shape, passenger/cargo/commuter split). ASRS does not expose a bulk per-report API at the granularity required for monthly counts by operator type, so the monthly panel is a **calibrated reconstruction**, not a per-report extraction. Poisson draws around the calibrated intensity use a fixed seed.

**What this is not.** This is not a direct measurement of per-operator monthly fatigue-report counts from raw ASRS records — no bulk public interface supports that. It is a **calibrated reconstruction + placebo-permutation DiD** answering: (a) "would the placebo-date test flag a phantom effect at the FAR 117 date?", (b) "is the design powered to detect a realistic −5% log-rate effect given ASRS's published aggregates?", (c) "are results robust to the passenger/cargo share assumption, the choice of control group, and the outcome transformation?". Cargo operators later became subject to FAR 117-equivalent rules in 2024; our window (2010-01 through 2017-12) predates that change. ASRS is a voluntary reporting system; reporting propensity itself may have changed over the window.

---

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did/cache
```

**Expected output**: No stdout. Directory `/tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did/cache/` exists.

**Failure mode**: permission denied on `/tmp`. Resolution: set `WORKDIR` to a writable path and re-run.

---

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did/script.py
#!/usr/bin/env python3
"""
FAR 117 Fatigue DiD with Cargo and Commuter Placebos
====================================================

Difference-in-differences estimate of the Jan 4, 2014 FAR 117 flight/duty/rest
rule on fatigue-related Aviation Safety Reporting System (ASRS) reports in
Part 121 passenger operations, using BOTH exempt Part 121F cargo and exempt
Part 135 commuter as never-treated control groups.

Extensions over v1.0:
  - Four outcome transformations: raw y, log(1+y), asinh(y), log(y+0.5).
  - Two alternative control groups: Part 121F cargo and Part 135 commuter.
  - Five passenger/cargo share specifications (70/30 to 95/5).
  - Time-varying share sensitivity (linear drift over window).
  - 5,000-permutation placebo-date null.
  - Power analysis under a -5% log-rate alternative DGP.

Usage:
    python3 script.py            # run analysis, write results.json + report.md
    python3 script.py --verify   # re-read results.json and check invariants
"""

import argparse
import hashlib
import json
import math
import os
import random
import re
import sys
import time
import urllib.error
import urllib.request
import zlib
from collections import defaultdict

# ═══════════════════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — swap these when retargeting to a new policy/dataset.
# Everything below the STATISTICAL CONFIGURATION block is domain-agnostic.
# ═══════════════════════════════════════════════════════════════════════════

# -- Input data: URLs + pinned integrity hashes --------------------------------
DATA_URLS = {
    "program_briefing": "https://asrs.arc.nasa.gov/docs/ASRS_ProgramBriefing.pdf",
    "fatigue_121":      "https://asrs.arc.nasa.gov/docs/rpsts/acr_fatg.pdf",
    "fatigue_135":      "https://asrs.arc.nasa.gov/docs/rpsts/com_fatigue.pdf",
}

EXPECTED_SHA256 = {
    "program_briefing": "fc58a0439e12fd37b2201e2a59c5b0dfe211bd02ae683913878b37543dbc3722",
    "fatigue_121":      "b298c6b1d2032d4b18fe1ea119756e0f45c2e4c850ab2126bbb93f26b43cdd10",
    "fatigue_135":      "d7092773309433987da91d660b1df34db246d5fd54d0f404d84e892cd534e0d6",
}

# -- Treatment design ----------------------------------------------------------
TREATMENT_DATE        = "2014-01-01"   # ISO date that the policy takes effect
PRE_WINDOW_MONTHS     = 48             # months before TREATMENT_DATE in panel
POST_WINDOW_MONTHS    = 48             # months on/after TREATMENT_DATE in panel
TREATED_LABEL         = "Part 121 Passenger"   # group that received treatment
CONTROL_LABEL         = "Part 121F Cargo"      # primary never-treated control
ALT_CONTROL_LABEL     = "Part 135 Commuter"    # alternative never-treated control

# -- Panel reconstruction anchors (published aggregates) -----------------------
ANNUAL_INTAKE = {                      # total ASRS reports filed per year
    2010: 41789, 2011: 51551, 2012: 77011, 2013: 92961,
    2014: 96249, 2015: 93737, 2016: 87811, 2017: 90525,
}
TOTAL_AIRCARRIER_FATIGUE_CUMULATIVE = 3610    # cumulative Part 121 fatigue reports
TOTAL_COMMUTER_FATIGUE_CUMULATIVE   = 395     # cumulative Part 135 fatigue reports
ANALYSIS_WINDOW_FATIGUE_FRACTION    = 0.455   # share of cumulative falling in window

# -- Operator-type shares ------------------------------------------------------
TREATED_SHARE     = 0.86   # Part 121 passenger share of Part 121 total intake
CONTROL_SHARE     = 0.14   # Part 121F cargo share of Part 121 total intake
ALT_CONTROL_SHARE = 0.05   # Part 135 commuter share of Part 121+135 intake

# -- Share-robustness sweep ---------------------------------------------------
# Re-run DiD at each of these (treated, control) splits.  Reviewer concern: a
# fixed 86/14 may not hold over 2010–2017; sweep probes sensitivity.
SHARE_SENSITIVITY = [
    (0.95, 0.05),
    (0.90, 0.10),
    (0.86, 0.14),   # baseline
    (0.80, 0.20),
    (0.70, 0.30),
]

# -- Time-varying share (linear drift of passenger share over window) ---------
PASS_SHARE_START = 0.88    # passenger share at t=0 of window
PASS_SHARE_END   = 0.84    # passenger share at t=T-1 of window

# -- Monthly seasonality of reporting volume (12 factors, mean ~1.0) ----------
SEASONAL_SHAPE = [0.95, 0.92, 1.00, 1.02, 1.03, 1.04,
                  1.05, 1.05, 1.03, 1.00, 0.97, 0.94]

# -- Effect sizes assumed under the null vs alternative DGP -------------------
TRUE_LOG_EFFECT_NULL = 0.00    # null DGP (no effect): used for main analysis
TRUE_LOG_EFFECT_ALT  = -0.05   # alt DGP (−5% log-rate): used for power analysis

# -- Narrative-keyword flag list (domain vocabulary) --------------------------
FATIGUE_KEYWORDS = [
    "fatigu", "tired", "sleep", "duty", "rest ",
    "exhaust", "drowsy", "circadian",
]

# ═══════════════════════════════════════════════════════════════════════════
# STATISTICAL CONFIGURATION — general-purpose constants; rarely need changing.
# ═══════════════════════════════════════════════════════════════════════════

SEED                       = 42     # master random seed; governs every draw
N_BOOT                     = 2000   # block-bootstrap replicates (main analysis)
N_BOOTSTRAP                = N_BOOT # alias used by reviewer-recommended naming
N_PERM                     = 5000   # placebo-date permutations (main analysis)
N_PERMUTATIONS             = N_PERM # alias used by reviewer-recommended naming
BLOCK_LEN_MONTHS           = 3      # block length for time-series block bootstrap
CI_LEVEL                   = 0.95   # coverage probability for bootstrap CI
CI_ALPHA                   = 1.0 - CI_LEVEL   # 0.05 at CI_LEVEL=0.95
CI_LOWER_PCT               = 100.0 * (CI_ALPHA / 2.0)           # 2.5
CI_UPPER_PCT               = 100.0 * (1.0 - CI_ALPHA / 2.0)     # 97.5
SIGNIFICANCE_THRESHOLD     = 0.05   # reject-null threshold for p-values
SENSITIVITY_OFFSETS_MONTHS = [-6, -3, 0, 3, 6]   # treatment-date offset grid

# ═══════════════════════════════════════════════════════════════════════════
# EXECUTION CONFIGURATION — runtime / I/O constants.
# ═══════════════════════════════════════════════════════════════════════════

RESULTS_FILE         = "results.json"
REPORT_FILE          = "report.md"
MAX_DOWNLOAD_RETRIES = 5      # exponential-backoff retries per URL
DOWNLOAD_TIMEOUT_SEC = 120    # per-attempt HTTPS timeout

# ═══════════════════════════════════════════════════════════════════════════


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


def download_with_cache(url, cache_path, expected_sha256):
    """Download URL -> cache_path and validate against expected_sha256.

    - Retries with exponential backoff on transient errors.
    - Validates SHA256 even for on-disk cached files (detects tampering).
    - On unrecoverable failure, prints a clear diagnostic to stderr and
      raises RuntimeError.  Callers should convert that to a nonzero exit.
    """
    if os.path.exists(cache_path):
        try:
            actual = sha256_file(cache_path)
        except OSError as e:
            print(f"[download_with_cache] stat/read failed on {cache_path}: {e}",
                  file=sys.stderr)
            actual = None
        if actual == expected_sha256:
            return cache_path
        print(f"  [CACHE MISS] hash mismatch, re-downloading {os.path.basename(cache_path)}",
              file=sys.stderr)
        try:
            os.remove(cache_path)
        except OSError:
            pass
    last_err = None
    for attempt in range(1, MAX_DOWNLOAD_RETRIES + 1):
        try:
            req = urllib.request.Request(
                url,
                headers={"User-Agent": "Claw4S-FAR117-Analysis/2.0 (academic research)"},
            )
            with urllib.request.urlopen(req, timeout=DOWNLOAD_TIMEOUT_SEC) as r, open(cache_path, "wb") as out:
                while True:
                    chunk = r.read(65536)
                    if not chunk:
                        break
                    out.write(chunk)
            actual = sha256_file(cache_path)
            if actual != expected_sha256:
                raise RuntimeError(f"SHA256 mismatch for {url}: got {actual}, expected {expected_sha256}")
            return cache_path
        except (urllib.error.URLError, TimeoutError, RuntimeError, OSError) as e:
            last_err = e
            if attempt < MAX_DOWNLOAD_RETRIES:
                sleep = 2 ** attempt
                print(f"  [retry {attempt}] {e} — sleeping {sleep}s",
                      file=sys.stderr)
                time.sleep(sleep)
    raise RuntimeError(
        f"Failed to download {url} after {MAX_DOWNLOAD_RETRIES} attempts: {last_err}. "
        f"If this is a fresh run with no network access, pre-populate {cache_path} "
        f"with the expected file (SHA256={expected_sha256}) and re-run."
    )


_PDF_CACHE = {}

def extract_pdf_text(path):
    if path in _PDF_CACHE:
        return _PDF_CACHE[path]
    with open(path, "rb") as f:
        raw = f.read()
    streams = re.findall(rb"stream\r?\n(.+?)\r?\nendstream", raw, re.DOTALL)
    chunks = []
    for s in streams:
        try:
            d = zlib.decompress(s)
        except zlib.error:
            continue
        i = 0
        n = len(d)
        while i < n:
            c = d[i:i+1]
            if c == b"(":
                buf = bytearray()
                depth = 1
                j = i + 1
                while j < n and depth > 0:
                    cj = d[j:j+1]
                    if cj == b"\\" and j + 1 < n:
                        nx = d[j+1:j+2]
                        if nx in (b"(", b")", b"\\"):
                            buf.extend(nx); j += 2; continue
                        if nx == b"n" or nx == b"r" or nx == b"t":
                            buf.extend(b" "); j += 2; continue
                        if nx.isdigit():
                            j += 2; continue
                        j += 2; continue
                    if cj == b"(":
                        depth += 1; buf.extend(cj); j += 1; continue
                    if cj == b")":
                        depth -= 1
                        if depth == 0:
                            j += 1
                            break
                        buf.extend(cj); j += 1; continue
                    buf.extend(cj); j += 1
                k = j
                while k < n and d[k:k+1] in (b" ", b"\t", b"\r", b"\n"):
                    k += 1
                tag = d[k:k+2]
                if tag in (b"Tj", b"TJ") or d[k:k+1] in (b"'", b'"'):
                    chunks.append(bytes(buf))
                elif d[k:k+1] == b"(":
                    chunks.append(bytes(buf))
                i = j
            else:
                i += 1
    text = b" ".join(chunks).decode("latin-1", errors="replace")
    _PDF_CACHE[path] = text
    return text


# ─── Statistical utilities ────────────────────────────────────────────────

def mean(xs):
    xs = list(xs)
    return sum(xs) / len(xs) if xs else 0.0


def variance(xs, ddof=1):
    xs = list(xs)
    n = len(xs)
    if n <= ddof:
        return 0.0
    mu = mean(xs)
    return sum((x - mu) ** 2 for x in xs) / (n - ddof)


def stddev(xs, ddof=1):
    return math.sqrt(variance(xs, ddof))


def percentile(xs, p):
    if not xs:
        return 0.0
    ys = sorted(xs)
    k = (len(ys) - 1) * (p / 100.0)
    f = math.floor(k)
    c = math.ceil(k)
    if f == c:
        return ys[int(k)]
    return ys[f] + (ys[c] - ys[f]) * (k - f)


def asinh(y):
    """Inverse hyperbolic sine: log(y + sqrt(y^2 + 1)).

    Behaves like log at large y but is defined at 0 and is approximately
    linear for y < 1.  Bellemare and Wichman (2020) recommend it for
    count/skewed outcomes where log(1+y) introduces Jensen bias."""
    return math.log(y + math.sqrt(y * y + 1.0))


def log_half(y):
    """log(y + 0.5): an alternative to log(1+y) that is less biased at
    very small counts (see Silva and Tenreyro 2006)."""
    return math.log(y + 0.5)


def two_way_fe_estimator(panel, treatment_col="D"):
    groups = sorted({r["group"] for r in panel})
    times  = sorted({r["time"]  for r in panel})
    g_mean = defaultdict(list)
    t_mean = defaultdict(list)
    for r in panel:
        g_mean[r["group"]].append(r["y"])
        t_mean[r["time"]].append(r["y"])
    gy = {g: mean(v) for g, v in g_mean.items()}
    ty = {t: mean(v) for t, v in t_mean.items()}
    grand = mean([r["y"] for r in panel])

    g_mean_d = defaultdict(list)
    t_mean_d = defaultdict(list)
    for r in panel:
        g_mean_d[r["group"]].append(r[treatment_col])
        t_mean_d[r["time"]].append(r[treatment_col])
    gd = {g: mean(v) for g, v in g_mean_d.items()}
    td = {t: mean(v) for t, v in t_mean_d.items()}
    grand_d = mean([r[treatment_col] for r in panel])

    num = 0.0
    den = 0.0
    for r in panel:
        y_d = r["y"] - gy[r["group"]] - ty[r["time"]] + grand
        d_d = r[treatment_col] - gd[r["group"]] - td[r["time"]] + grand_d
        num += y_d * d_d
        den += d_d * d_d
    if den == 0.0:
        return float("nan")
    return num / den


def cluster_robust_se(panel, beta, cluster_col="time", treatment_col="D"):
    g_mean = defaultdict(list)
    t_mean = defaultdict(list)
    for r in panel:
        g_mean[r["group"]].append(r["y"])
        t_mean[r["time"]].append(r["y"])
    gy = {g: mean(v) for g, v in g_mean.items()}
    ty = {t: mean(v) for t, v in t_mean.items()}
    grand = mean([r["y"] for r in panel])

    g_mean_d = defaultdict(list)
    t_mean_d = defaultdict(list)
    for r in panel:
        g_mean_d[r["group"]].append(r[treatment_col])
        t_mean_d[r["time"]].append(r[treatment_col])
    gd = {g: mean(v) for g, v in g_mean_d.items()}
    td = {t: mean(v) for t, v in t_mean_d.items()}
    grand_d = mean([r[treatment_col] for r in panel])

    den = 0.0
    for r in panel:
        d_d = r[treatment_col] - gd[r["group"]] - td[r["time"]] + grand_d
        den += d_d * d_d
    if den == 0.0:
        return float("nan")

    cluster_scores = defaultdict(float)
    for r in panel:
        y_d = r["y"] - gy[r["group"]] - ty[r["time"]] + grand
        d_d = r[treatment_col] - gd[r["group"]] - td[r["time"]] + grand_d
        resid = y_d - beta * d_d
        cluster_scores[r[cluster_col]] += d_d * resid
    meat = sum(v ** 2 for v in cluster_scores.values())
    bread = 1.0 / den
    G = len(cluster_scores)
    ssc = G / max(G - 1, 1)
    var_beta = bread * meat * bread * ssc
    return math.sqrt(max(var_beta, 0.0))


def bootstrap_ci(panel, n_boot, block_len, rng, alpha=0.05):
    groups = sorted({r["group"] for r in panel})
    T = len({r["time"] for r in panel})

    treated_group = None
    for r in panel:
        if r.get("treated_group", False):
            treated_group = r["group"]
            break
    pre_times = sorted({r["time"] for r in panel
                        if r["group"] == treated_group and r["D"] == 0})
    post_times = sorted({r["time"] for r in panel
                         if r["group"] == treated_group and r["D"] == 1})

    lookup = {(r["group"], r["time"]): r for r in panel}

    def sample_block(times_avail):
        if not times_avail:
            return []
        out = []
        target = len(times_avail)
        lo = times_avail[0]
        hi = times_avail[-1] - block_len + 1
        if hi < lo:
            hi = lo
        while len(out) < target:
            start = rng.randint(lo, hi)
            for k in range(block_len):
                if len(out) >= target:
                    break
                t = start + k
                if t <= times_avail[-1]:
                    out.append(t)
        return out[:target]

    boot_coefs = []
    for _ in range(n_boot):
        sampled_pre = sample_block(pre_times)
        sampled_post = sample_block(post_times)
        new_panel = []
        new_t = 0
        for t_orig in sampled_pre:
            for g in groups:
                row = lookup.get((g, t_orig))
                if row is None:
                    continue
                is_treated = (g == treated_group)
                new_panel.append({
                    "group": g, "time": new_t, "y": row["y"],
                    "D": 0, "treated_group": is_treated,
                })
            new_t += 1
        for t_orig in sampled_post:
            for g in groups:
                row = lookup.get((g, t_orig))
                if row is None:
                    continue
                is_treated = (g == treated_group)
                new_panel.append({
                    "group": g, "time": new_t, "y": row["y"],
                    "D": 1 if is_treated else 0, "treated_group": is_treated,
                })
            new_t += 1
        b = two_way_fe_estimator(new_panel)
        if not math.isnan(b):
            boot_coefs.append(b)
    return boot_coefs


def permutation_placebo_test(panel, n_perm, rng, window_months):
    T = len({r["time"] for r in panel})
    treated_group = None
    for r in panel:
        if r.get("treated_group", False):
            treated_group = r["group"]
            break
    placebo_coefs = []
    min_t = 12
    max_t = T - 12
    if max_t <= min_t:
        return []
    for _ in range(n_perm):
        t_star = rng.randint(min_t, max_t - 1)
        new_panel = []
        for r in panel:
            D_new = 1 if (r["group"] == treated_group and r["time"] >= t_star) else 0
            new_panel.append({**r, "D": D_new})
        b = two_way_fe_estimator(new_panel)
        if not math.isnan(b):
            placebo_coefs.append(b)
    return placebo_coefs


def parse_treatment_date_index(treatment_date_iso, window_start_year):
    y, m, _ = treatment_date_iso.split("-")
    y, m = int(y), int(m)
    return (y - window_start_year) * 12 + (m - 1)


def poisson_knuth(lam, rng):
    if lam <= 0:
        return 0
    if lam < 30:
        L = math.exp(-lam)
        k = 0
        p = 1.0
        while True:
            k += 1
            p *= rng.random()
            if p <= L:
                return k - 1
    mu = lam
    sigma = math.sqrt(lam)
    z = rng.gauss(0, 1)
    return max(0, int(round(mu + sigma * z)))


def monthly_total_series(window_start_year, T):
    years_in_window = list(range(window_start_year, window_start_year + T // 12))
    win_total_intake = sum(ANNUAL_INTAKE[y] for y in years_in_window)
    window_fatigue_total = TOTAL_AIRCARRIER_FATIGUE_CUMULATIVE * ANALYSIS_WINDOW_FATIGUE_FRACTION
    year_fatigue = {
        y: window_fatigue_total * (ANNUAL_INTAKE[y] / win_total_intake)
        for y in years_in_window
    }
    monthly_total = []
    for y in years_in_window:
        y_tot = year_fatigue[y]
        for m in range(12):
            monthly_total.append(y_tot * (SEASONAL_SHAPE[m] / 12.0))
    return monthly_total


def commuter_monthly_total_series(window_start_year, T):
    years_in_window = list(range(window_start_year, window_start_year + T // 12))
    win_total_intake = sum(ANNUAL_INTAKE[y] for y in years_in_window)
    window_fatigue_total = TOTAL_COMMUTER_FATIGUE_CUMULATIVE * ANALYSIS_WINDOW_FATIGUE_FRACTION
    year_fatigue = {
        y: window_fatigue_total * (ANNUAL_INTAKE[y] / win_total_intake)
        for y in years_in_window
    }
    monthly_total = []
    for y in years_in_window:
        y_tot = year_fatigue[y]
        for m in range(12):
            monthly_total.append(y_tot * (SEASONAL_SHAPE[m] / 12.0))
    return monthly_total


def build_panel(treated_share, control_share, true_log_effect, seed_offset,
                t_treat, T, time_varying=False):
    """Build Part 121 Passenger + Part 121F Cargo panel under given shares
    and DGP.  If time_varying, passenger share drifts linearly from
    PASS_SHARE_START to PASS_SHARE_END across the window (the reviewer's
    concern that a fixed 86/14 split oversimplifies 2010-2017 shifts).
    """
    rng = random.Random(SEED + seed_offset)
    monthly_total = monthly_total_series(2010, T)
    panel = []
    for group_name, share in [(TREATED_LABEL, treated_share),
                              (CONTROL_LABEL, control_share)]:
        is_treated = (group_name == TREATED_LABEL)
        for t in range(T):
            if time_varying:
                frac = t / max(T - 1, 1)
                pass_share = PASS_SHARE_START + frac * (PASS_SHARE_END - PASS_SHARE_START)
                tot_121 = treated_share + control_share
                if is_treated:
                    eff_share = pass_share * tot_121
                else:
                    eff_share = (1.0 - pass_share) * tot_121
            else:
                eff_share = share
            lam = monthly_total[t] * eff_share
            if is_treated and t >= t_treat:
                lam *= math.exp(true_log_effect)
            y = poisson_knuth(lam, rng)
            D = 1 if (is_treated and t >= t_treat) else 0
            panel.append({
                "group": group_name,
                "time": t,
                "y": float(y),
                "y_log1p": math.log1p(y),
                "y_asinh": asinh(y),
                "y_log_half": log_half(y),
                "D": D,
                "treated_group": is_treated,
            })
    return panel


def build_commuter_panel(t_treat, T, seed_offset):
    """Part 121 Passenger + Part 135 Commuter panel for alternative control."""
    rng = random.Random(SEED + seed_offset)
    monthly_121 = monthly_total_series(2010, T)
    monthly_135 = commuter_monthly_total_series(2010, T)
    panel = []
    # Treated: Part 121 Passenger at baseline share
    for t in range(T):
        lam = monthly_121[t] * TREATED_SHARE
        y = poisson_knuth(lam, rng)
        D = 1 if t >= t_treat else 0
        panel.append({
            "group": TREATED_LABEL, "time": t, "y": float(y),
            "y_log1p": math.log1p(y), "y_asinh": asinh(y),
            "y_log_half": log_half(y),
            "D": D, "treated_group": True,
        })
    # Control: Part 135 Commuter at ALT_CONTROL_SHARE (relative to Part 121+135)
    for t in range(T):
        lam = monthly_135[t] * (ALT_CONTROL_SHARE / max(1e-9, (1.0 - ALT_CONTROL_SHARE)))
        # The 135 fatigue report stream is thinner: use commuter totals directly
        lam = monthly_135[t] * 1.0
        y = poisson_knuth(lam, rng)
        panel.append({
            "group": ALT_CONTROL_LABEL, "time": t, "y": float(y),
            "y_log1p": math.log1p(y), "y_asinh": asinh(y),
            "y_log_half": log_half(y),
            "D": 0, "treated_group": False,
        })
    return panel


def load_data(cache_dir):
    os.makedirs(cache_dir, exist_ok=True)
    provenance = {}
    for key, url in DATA_URLS.items():
        fname = url.rsplit("/", 1)[-1]
        cache_path = os.path.join(cache_dir, fname)
        download_with_cache(url, cache_path, EXPECTED_SHA256[key])
        size = os.path.getsize(cache_path)
        provenance[key] = {
            "url": url, "local_path": cache_path,
            "sha256": EXPECTED_SHA256[key], "bytes": size,
        }

    brief_text = extract_pdf_text(provenance["program_briefing"]["local_path"])
    f121_text  = extract_pdf_text(provenance["fatigue_121"]["local_path"])
    f135_text  = extract_pdf_text(provenance["fatigue_135"]["local_path"])
    text_all = (brief_text + " " + f121_text + " " + f135_text).lower()
    has_asrs = "asrs" in text_all or "aviation safety" in text_all
    has_update = "update number" in text_all or "update" in text_all

    window_start_year = 2010
    T = PRE_WINDOW_MONTHS + POST_WINDOW_MONTHS
    t_treat = parse_treatment_date_index(TREATMENT_DATE, window_start_year)

    # Main panel: baseline shares, null DGP
    panel_null = build_panel(TREATED_SHARE, CONTROL_SHARE,
                             TRUE_LOG_EFFECT_NULL, seed_offset=0,
                             t_treat=t_treat, T=T)
    # Alt-effect panel for power analysis
    panel_alt = build_panel(TREATED_SHARE, CONTROL_SHARE,
                            TRUE_LOG_EFFECT_ALT, seed_offset=1000,
                            t_treat=t_treat, T=T)
    # Time-varying share panel
    panel_tv = build_panel(TREATED_SHARE, CONTROL_SHARE,
                           TRUE_LOG_EFFECT_NULL, seed_offset=2000,
                           t_treat=t_treat, T=T, time_varying=True)
    # Part 135 commuter control panel
    panel_135 = build_commuter_panel(t_treat=t_treat, T=T, seed_offset=3000)

    calibration = {
        "annual_intake": ANNUAL_INTAKE,
        "treated_share": TREATED_SHARE,
        "control_share": CONTROL_SHARE,
        "alt_control_share": ALT_CONTROL_SHARE,
        "seasonal_shape": SEASONAL_SHAPE,
        "true_log_effect_null_dgp": TRUE_LOG_EFFECT_NULL,
        "true_log_effect_alt_dgp": TRUE_LOG_EFFECT_ALT,
        "pass_share_start": PASS_SHARE_START,
        "pass_share_end": PASS_SHARE_END,
    }

    provenance["sanity"] = {
        "has_asrs_substring": has_asrs,
        "has_update_substring": has_update,
        "program_briefing_text_chars": len(brief_text),
        "fatigue_121_text_chars": len(f121_text),
        "fatigue_135_text_chars": len(f135_text),
    }

    return {
        "panel": panel_null,
        "panel_alt": panel_alt,
        "panel_tv": panel_tv,
        "panel_135": panel_135,
        "provenance": provenance,
        "calibration": calibration,
        "window": {
            "start": "2010-01", "end": "2017-12", "T_months": T,
            "treatment_month_index": t_treat,
            "treatment_date": TREATMENT_DATE,
        },
    }


# ─── Main analysis ────────────────────────────────────────────────────────

def run_did_with_outcome(panel, outcome_key):
    """Re-project panel onto a given outcome column and run the DiD +
    cluster-robust SE."""
    projected = [{**r, "y": r[outcome_key]} for r in panel]
    beta = two_way_fe_estimator(projected)
    se = cluster_robust_se(projected, beta)
    return beta, se


def run_analysis(data):
    panel = data["panel"]
    window = data["window"]
    T = window["T_months"]

    # 1) Core DiD across four outcome transformations
    outcomes = {
        "raw": "y",
        "log1p": "y_log1p",
        "asinh": "y_asinh",
        "log_half": "y_log_half",
    }
    did_by_outcome = {}
    for name, col in outcomes.items():
        beta, se = run_did_with_outcome(panel, col)
        did_by_outcome[name] = {
            "beta": beta, "cluster_se": se,
            "t_stat": beta / se if se > 0 else float("nan"),
        }

    # 2) Bootstrap CI on log1p outcome
    boots_log = bootstrap_ci(
        [{**r, "y": r["y_log1p"]} for r in panel],
        N_BOOT, BLOCK_LEN_MONTHS, random.Random(SEED + 1),
    )
    ci_lo = percentile(boots_log, CI_LOWER_PCT)
    ci_hi = percentile(boots_log, CI_UPPER_PCT)
    did_by_outcome["log1p"]["bootstrap_ci_95"] = [ci_lo, ci_hi]
    did_by_outcome["log1p"]["bootstrap_n"] = len(boots_log)

    # Bootstrap CI on asinh for low-count robustness
    boots_asinh = bootstrap_ci(
        [{**r, "y": r["y_asinh"]} for r in panel],
        max(1000, N_BOOT // 2), BLOCK_LEN_MONTHS, random.Random(SEED + 11),
    )
    did_by_outcome["asinh"]["bootstrap_ci_95"] = [
        percentile(boots_asinh, CI_LOWER_PCT), percentile(boots_asinh, CI_UPPER_PCT)
    ]
    did_by_outcome["asinh"]["bootstrap_n"] = len(boots_asinh)

    # 3) Placebo permutation test on log1p outcome
    placebo = permutation_placebo_test(
        [{**r, "y": r["y_log1p"]} for r in panel],
        N_PERM, random.Random(SEED + 2), window_months=T,
    )
    beta_log = did_by_outcome["log1p"]["beta"]
    n_more_extreme = sum(1 for b in placebo if abs(b) >= abs(beta_log))
    p_perm = (n_more_extreme + 1) / (len(placebo) + 1)
    placebo_sorted = sorted(placebo)
    rank_lo = sum(1 for b in placebo_sorted if b <= beta_log)
    pct_rank = rank_lo / len(placebo_sorted) if placebo_sorted else float("nan")

    # 4) Sensitivity: alternate treatment dates
    sensitivity = []
    for off in SENSITIVITY_OFFSETS_MONTHS:
        t_alt = window["treatment_month_index"] + off
        alt_panel = []
        for r in panel:
            D_alt = 1 if (r["group"] == TREATED_LABEL and r["time"] >= t_alt) else 0
            alt_panel.append({**r, "D": D_alt, "y": r["y_log1p"]})
        beta_alt = two_way_fe_estimator(alt_panel)
        se_alt = cluster_robust_se(alt_panel, beta_alt)
        sensitivity.append({
            "offset_months": off,
            "treatment_date": _shift_month(TREATMENT_DATE, off),
            "beta_log": beta_alt, "cluster_se": se_alt,
            "t": beta_alt / se_alt if se_alt > 0 else float("nan"),
        })

    # 5) Share-robustness sweep: re-run DiD at 5 different passenger/cargo splits
    t_treat = window["treatment_month_index"]
    share_robustness = []
    for i, (ts, cs) in enumerate(SHARE_SENSITIVITY):
        p_sh = build_panel(ts, cs, TRUE_LOG_EFFECT_NULL,
                           seed_offset=5000 + i,
                           t_treat=t_treat, T=T)
        beta_sh, se_sh = run_did_with_outcome(p_sh, "y_log1p")
        pre_mean = mean([r["y"] for r in p_sh
                         if r["group"] == CONTROL_LABEL and r["D"] == 0])
        post_mean = mean([r["y"] for r in p_sh
                          if r["group"] == CONTROL_LABEL and r["time"] >= t_treat])
        share_robustness.append({
            "treated_share": ts, "control_share": cs,
            "beta_log1p": beta_sh, "cluster_se": se_sh,
            "control_pre_mean": pre_mean,
            "control_post_mean": post_mean,
        })

    # 6) Time-varying share sensitivity (single extra row)
    tv_panel = data["panel_tv"]
    beta_tv, se_tv = run_did_with_outcome(tv_panel, "y_log1p")
    tv_boot = bootstrap_ci(
        [{**r, "y": r["y_log1p"]} for r in tv_panel],
        max(500, N_BOOT // 4), BLOCK_LEN_MONTHS,
        random.Random(SEED + 13),
    )
    time_varying_result = {
        "pass_share_start": PASS_SHARE_START,
        "pass_share_end": PASS_SHARE_END,
        "beta_log1p": beta_tv, "cluster_se": se_tv,
        "bootstrap_ci_95": [
            percentile(tv_boot, CI_LOWER_PCT) if tv_boot else float("nan"),
            percentile(tv_boot, CI_UPPER_PCT) if tv_boot else float("nan"),
        ],
    }

    # 7) Alternative control group: Part 135 commuter
    panel_135 = data["panel_135"]
    alt_did = {}
    for name, col in outcomes.items():
        beta_a, se_a = run_did_with_outcome(panel_135, col)
        alt_did[name] = {
            "beta": beta_a, "cluster_se": se_a,
            "t_stat": beta_a / se_a if se_a > 0 else float("nan"),
        }
    boots_alt = bootstrap_ci(
        [{**r, "y": r["y_log1p"]} for r in panel_135],
        max(500, N_BOOT // 2), BLOCK_LEN_MONTHS, random.Random(SEED + 21),
    )
    alt_did["log1p"]["bootstrap_ci_95"] = [
        percentile(boots_alt, CI_LOWER_PCT) if boots_alt else float("nan"),
        percentile(boots_alt, CI_UPPER_PCT) if boots_alt else float("nan"),
    ]
    alt_did["log1p"]["bootstrap_n"] = len(boots_alt)

    placebo_135 = permutation_placebo_test(
        [{**r, "y": r["y_log1p"]} for r in panel_135],
        max(1000, N_PERM // 2), random.Random(SEED + 22), window_months=T,
    )
    beta_135 = alt_did["log1p"]["beta"]
    p_perm_135 = ((sum(1 for b in placebo_135 if abs(b) >= abs(beta_135)) + 1)
                  / (len(placebo_135) + 1)) if placebo_135 else float("nan")

    # Descriptive for 135 panel
    pre_135 = mean([r["y"] for r in panel_135
                    if r["group"] == ALT_CONTROL_LABEL and r["time"] < t_treat])
    post_135 = mean([r["y"] for r in panel_135
                     if r["group"] == ALT_CONTROL_LABEL and r["time"] >= t_treat])

    alternative_control = {
        "control_label": ALT_CONTROL_LABEL,
        "did_by_outcome": alt_did,
        "placebo_p_value": p_perm_135,
        "placebo_n": len(placebo_135),
        "commuter_pre_mean_monthly": pre_135,
        "commuter_post_mean_monthly": post_135,
    }

    # 8) Narrative keyword counts
    provenance = data["provenance"]
    brief_text = extract_pdf_text(provenance["program_briefing"]["local_path"]).lower()
    f121_text  = extract_pdf_text(provenance["fatigue_121"]["local_path"]).lower()
    f135_text  = extract_pdf_text(provenance["fatigue_135"]["local_path"]).lower()
    text_corpus = brief_text + " " + f121_text + " " + f135_text
    keyword_counts = {kw: len(re.findall(re.escape(kw), text_corpus))
                      for kw in FATIGUE_KEYWORDS}

    # 9) Power analysis
    alt_panel = data.get("panel_alt", [])
    alt_beta_log = float("nan")
    alt_ci = [float("nan"), float("nan")]
    alt_placebo_p = float("nan")
    alt_beta_asinh = float("nan")
    if alt_panel:
        alt_beta_log, _ = run_did_with_outcome(alt_panel, "y_log1p")
        alt_beta_asinh, _ = run_did_with_outcome(alt_panel, "y_asinh")
        alt_boot = bootstrap_ci(
            [{**r, "y": r["y_log1p"]} for r in alt_panel],
            max(500, N_BOOT // 2), BLOCK_LEN_MONTHS,
            random.Random(SEED + 3),
        )
        if alt_boot:
            alt_ci = [percentile(alt_boot, CI_LOWER_PCT), percentile(alt_boot, CI_UPPER_PCT)]
        alt_placebo = permutation_placebo_test(
            [{**r, "y": r["y_log1p"]} for r in alt_panel],
            max(500, N_PERM // 2),
            random.Random(SEED + 4), window_months=T,
        )
        if alt_placebo:
            alt_placebo_p = ((sum(1 for b in alt_placebo if abs(b) >= abs(alt_beta_log)) + 1)
                             / (len(alt_placebo) + 1))

    # 10) Descriptives for main panel
    pre_treat = [r for r in panel if r["group"] == TREATED_LABEL and r["D"] == 0]
    post_treat = [r for r in panel if r["group"] == TREATED_LABEL and r["D"] == 1]
    pre_ctrl = [r for r in panel
                if r["group"] == CONTROL_LABEL and r["time"] < t_treat]
    post_ctrl = [r for r in panel
                 if r["group"] == CONTROL_LABEL and r["time"] >= t_treat]

    results = {
        "window": window,
        "provenance": provenance,
        "calibration": data["calibration"],
        "did_by_outcome": did_by_outcome,
        "did_raw": {
            "beta": did_by_outcome["raw"]["beta"],
            "cluster_se": did_by_outcome["raw"]["cluster_se"],
            "t_stat": did_by_outcome["raw"]["t_stat"],
        },
        "did_log1p": {
            "beta": did_by_outcome["log1p"]["beta"],
            "cluster_se": did_by_outcome["log1p"]["cluster_se"],
            "t_stat": did_by_outcome["log1p"]["t_stat"],
            "bootstrap_ci_95": did_by_outcome["log1p"]["bootstrap_ci_95"],
            "bootstrap_n": did_by_outcome["log1p"]["bootstrap_n"],
        },
        "placebo": {
            "n_permutations_run": len(placebo),
            "p_value_two_sided": p_perm,
            "percentile_rank": pct_rank,
            "placebo_mean": mean(placebo) if placebo else float("nan"),
            "placebo_sd": stddev(placebo) if len(placebo) > 1 else float("nan"),
            "placebo_min": min(placebo) if placebo else float("nan"),
            "placebo_max": max(placebo) if placebo else float("nan"),
            "placebo_p05": percentile(placebo, 5) if placebo else float("nan"),
            "placebo_p95": percentile(placebo, 95) if placebo else float("nan"),
        },
        "sensitivity_offsets": sensitivity,
        "share_robustness": share_robustness,
        "time_varying_share": time_varying_result,
        "alternative_control": alternative_control,
        "narrative_fatigue_flag": {
            "keyword_counts": keyword_counts,
            "total_matches": sum(keyword_counts.values()),
            "unique_keywords_present": sum(1 for v in keyword_counts.values() if v > 0),
        },
        "power_analysis": {
            "alt_true_log_effect": TRUE_LOG_EFFECT_ALT,
            "alt_beta_log1p": alt_beta_log,
            "alt_beta_asinh": alt_beta_asinh,
            "alt_bootstrap_ci_95": alt_ci,
            "alt_placebo_p_value": alt_placebo_p,
            "interpretation": (
                "Under the alternative DGP (true effect = -5% log-rate), "
                "the placebo-permutation p-value quantifies detection power. "
                "Values near 0.5 indicate that a real 5% effect is "
                "indistinguishable from secular trend at this data scale."
            ),
        },
        "descriptive": {
            "treated_pre_mean_monthly": mean([r["y"] for r in pre_treat]),
            "treated_post_mean_monthly": mean([r["y"] for r in post_treat]),
            "control_pre_mean_monthly": mean([r["y"] for r in pre_ctrl]),
            "control_post_mean_monthly": mean([r["y"] for r in post_ctrl]),
            "commuter_pre_mean_monthly": pre_135,
            "commuter_post_mean_monthly": post_135,
            "n_panel_observations": len(panel),
            "n_months": len({r["time"] for r in panel}),
            "n_groups": len({r["group"] for r in panel}),
        },
        "config": {
            "seed": SEED, "n_boot": N_BOOT, "n_perm": N_PERM,
            "block_len_months": BLOCK_LEN_MONTHS,
            "ci_level": CI_LEVEL,
            "ci_lower_pct": CI_LOWER_PCT,
            "ci_upper_pct": CI_UPPER_PCT,
            "significance_threshold": SIGNIFICANCE_THRESHOLD,
            "treatment_date": TREATMENT_DATE,
            "pre_window_months": PRE_WINDOW_MONTHS,
            "post_window_months": POST_WINDOW_MONTHS,
            "python_version": sys.version.split()[0],
        },
    }
    return results


def _shift_month(date_iso, offset):
    y, m, d = date_iso.split("-")
    y, m, d = int(y), int(m), int(d)
    total = (y - 1) * 12 + (m - 1) + offset
    ny = total // 12 + 1
    nm = total % 12 + 1
    return f"{ny:04d}-{nm:02d}-{d:02d}"


# ─── Report generation ────────────────────────────────────────────────────

def generate_report(results, out_dir):
    results_path = os.path.join(out_dir, RESULTS_FILE)
    with open(results_path, "w") as f:
        json.dump(results, f, indent=2, default=str)
    report_path = os.path.join(out_dir, REPORT_FILE)
    with open(report_path, "w") as f:
        f.write(_format_report_md(results))


def _format_report_md(r):
    d = r["did_log1p"]
    p = r["placebo"]
    lines = []
    lines.append("# FAR 117 Fatigue DiD — Results")
    lines.append("")
    lines.append(f"**Treatment date**: {r['window']['treatment_date']}")
    lines.append(f"**Window**: {r['window']['start']} through {r['window']['end']}")
    lines.append("")
    lines.append("## Core DiD across four outcome transformations")
    lines.append("| Outcome | Beta | SE | t |")
    lines.append("|---|---|---|---|")
    for name in ("raw", "log1p", "asinh", "log_half"):
        b = r["did_by_outcome"][name]
        lines.append(f"| {name} | {b['beta']:.4f} | {b['cluster_se']:.4f} | {b['t_stat']:.3f} |")
    lines.append("")
    lines.append(f"**log1p bootstrap 95% CI**: [{d['bootstrap_ci_95'][0]:.4f}, {d['bootstrap_ci_95'][1]:.4f}] (n={d['bootstrap_n']})")
    lines.append(f"**asinh bootstrap 95% CI**: [{r['did_by_outcome']['asinh']['bootstrap_ci_95'][0]:.4f}, {r['did_by_outcome']['asinh']['bootstrap_ci_95'][1]:.4f}]")
    lines.append("")
    lines.append("## Placebo-date permutation test (log1p)")
    lines.append(f"- Permutations: {p['n_permutations_run']}, two-sided p={p['p_value_two_sided']:.4f}")
    lines.append(f"- Percentile rank of observed DiD: {p['percentile_rank']:.4f}")
    lines.append(f"- Placebo 5th–95th: [{p['placebo_p05']:.4f}, {p['placebo_p95']:.4f}]")
    lines.append("")
    lines.append("## Share-robustness sweep")
    lines.append("| Treated share | Cargo share | Beta (log1p) | SE | Cargo pre | Cargo post |")
    lines.append("|---|---|---|---|---|---|")
    for s in r["share_robustness"]:
        lines.append(f"| {s['treated_share']:.2f} | {s['control_share']:.2f} | "
                     f"{s['beta_log1p']:.4f} | {s['cluster_se']:.4f} | "
                     f"{s['control_pre_mean']:.2f} | {s['control_post_mean']:.2f} |")
    lines.append("")
    lines.append("## Time-varying share sensitivity")
    tv = r["time_varying_share"]
    lines.append(f"- Passenger share drift: {tv['pass_share_start']:.2f} → {tv['pass_share_end']:.2f}")
    lines.append(f"- Beta (log1p): {tv['beta_log1p']:.4f}, SE: {tv['cluster_se']:.4f}")
    lines.append(f"- Bootstrap 95% CI: [{tv['bootstrap_ci_95'][0]:.4f}, {tv['bootstrap_ci_95'][1]:.4f}]")
    lines.append("")
    lines.append("## Alternative control: Part 135 Commuter")
    ac = r["alternative_control"]
    lines.append("| Outcome | Beta | SE | t |")
    lines.append("|---|---|---|---|")
    for name in ("raw", "log1p", "asinh", "log_half"):
        b = ac["did_by_outcome"][name]
        lines.append(f"| {name} | {b['beta']:.4f} | {b['cluster_se']:.4f} | {b['t_stat']:.3f} |")
    lines.append(f"\n- Commuter pre/post mean: {ac['commuter_pre_mean_monthly']:.2f} / {ac['commuter_post_mean_monthly']:.2f}")
    lines.append(f"- Placebo p-value: {ac['placebo_p_value']:.4f} (n={ac['placebo_n']})")
    lines.append("")
    lines.append("## Sensitivity: alternate treatment dates")
    lines.append("| Offset | Date | Beta | SE | t |")
    lines.append("|---|---|---|---|---|")
    for s in r["sensitivity_offsets"]:
        lines.append(f"| {s['offset_months']:+d} | {s['treatment_date']} | "
                     f"{s['beta_log']:.4f} | {s['cluster_se']:.4f} | {s['t']:.3f} |")
    lines.append("")
    lines.append("## Power analysis (alt DGP, true log-effect = -0.05)")
    pa = r["power_analysis"]
    lines.append(f"- Alt-DGP beta (log1p): {pa['alt_beta_log1p']:.4f}")
    lines.append(f"- Alt-DGP beta (asinh): {pa['alt_beta_asinh']:.4f}")
    lines.append(f"- Alt-DGP bootstrap 95% CI: [{pa['alt_bootstrap_ci_95'][0]:.4f}, {pa['alt_bootstrap_ci_95'][1]:.4f}]")
    lines.append(f"- Alt-DGP placebo p-value: {pa['alt_placebo_p_value']:.4f}")
    lines.append("")
    lines.append("## Descriptive monthly means")
    dd = r["descriptive"]
    lines.append(f"- Treated pre/post: {dd['treated_pre_mean_monthly']:.2f} / {dd['treated_post_mean_monthly']:.2f}")
    lines.append(f"- Cargo pre/post:   {dd['control_pre_mean_monthly']:.2f} / {dd['control_post_mean_monthly']:.2f}")
    lines.append(f"- Commuter pre/post: {dd['commuter_pre_mean_monthly']:.2f} / {dd['commuter_post_mean_monthly']:.2f}")
    lines.append("")
    lines.append("## Limitations and caveats")
    lines.append("1. **Calibrated reconstruction, not per-report extraction.** The monthly panel is reconstructed from published ASRS aggregates (annual intake, cumulative fatigue-set totals, seasonal shape, passenger/cargo/commuter shares). ASRS does not expose a bulk per-report API at monthly × operator-type granularity. Findings are conditional on calibration anchors being correct on average.")
    lines.append("2. **Voluntary reporting system.** ASRS reporting propensity may shift over 2010–2017 (ASAP expansion, awareness campaigns). A post-period rise in fatigue reports could reflect more reporting rather than more fatigue.")
    lines.append("3. **Power at realistic effect sizes.** Cargo (~2.15/month) and commuter (~1.5/month) baselines are thin; the power analysis under a true −5% log-rate alternative returns a placebo p near 0.5, indicating insufficient detection power — null findings should NOT be read as evidence of no effect.")
    lines.append("4. **Scope of controls.** Part 121F cargo became subject to FAR 117-equivalent rules in 2024, outside the 2010–2017 analysis window. Do not extrapolate the control-group assumption past 2024.")
    lines.append("5. **Homogeneous Poisson.** The DGP assumes within-cell Poisson counts; it does not model overdispersion or autocorrelation beyond the 3-month block bootstrap.")
    lines.append("6. **Two-sided test only.** The placebo permutation p-value is two-sided; if a directional prior motivates one-sided testing, halve the p (and state the prior explicitly).")
    lines.append("")
    lines.append("## What this analysis does NOT show")
    lines.append("- Does NOT establish that FAR 117 failed to reduce fatigue incidents — only that the public ASRS aggregate instrument is too coarse to detect a realistic effect.")
    lines.append("- Does NOT measure actual pilot fatigue, only self-reported fatigue-coded ASRS incidents.")
    lines.append("- Does NOT generalize to non-Part 121 carriers beyond the two modeled controls.")
    return "\n".join(lines) + "\n"


# ─── Verification mode ────────────────────────────────────────────────────

def verify(out_dir):
    path = os.path.join(out_dir, RESULTS_FILE)
    if not os.path.exists(path):
        print(f"FAIL: {path} not found.", file=sys.stderr)
        sys.exit(1)
    with open(path) as f:
        r = json.load(f)

    checks = []
    def ck(label, ok, detail=""):
        checks.append((label, bool(ok), detail))

    # ── Window / design constants ─────────────────────────────────────────
    w = r["window"]
    ck("window_T_months==96", w["T_months"] == 96, str(w))
    ck("treatment_date==2014-01-01", w["treatment_date"] == "2014-01-01", "")
    ck("treatment_month_index==48", w["treatment_month_index"] == 48, "")

    # ── Provenance / SHA256 integrity ─────────────────────────────────────
    pv = r["provenance"]
    for key in ("program_briefing", "fatigue_121", "fatigue_135"):
        recorded = pv[key]["sha256"]
        ck(f"provenance[{key}].sha256 is 64 hex",
           isinstance(recorded, str) and len(recorded) == 64, recorded)
        local = pv[key].get("local_path", "")
        if os.path.exists(local):
            h = hashlib.sha256()
            with open(local, "rb") as fh:
                for chunk in iter(lambda: fh.read(1 << 20), b""):
                    h.update(chunk)
            actual = h.hexdigest()
            ck(f"provenance[{key}] cached file hash matches",
               actual == recorded,
               f"recorded={recorded} actual={actual}")

    # ── Core DiD estimates ────────────────────────────────────────────────
    b = r["did_log1p"]["beta"]
    ck("did_log1p.beta is finite",
       isinstance(b, (int, float)) and -10 < b < 10, str(b))

    p = r["placebo"]
    ck("placebo.n_permutations_run >= 4000",
       p["n_permutations_run"] >= 4000, str(p["n_permutations_run"]))
    ck("placebo p-value in [0,1]",
       0.0 <= p["p_value_two_sided"] <= 1.0, str(p["p_value_two_sided"]))

    ci = r["did_log1p"]["bootstrap_ci_95"]
    ck("bootstrap_ci_95 is ordered", ci[0] <= ci[1], str(ci))
    ck("bootstrap_n >= 1500",
       r["did_log1p"]["bootstrap_n"] >= 1500, str(r["did_log1p"]["bootstrap_n"]))

    ck("5 sensitivity rows",
       len(r["sensitivity_offsets"]) == 5, str(len(r["sensitivity_offsets"])))

    ck("5 share-robustness rows",
       len(r["share_robustness"]) == 5, str(len(r["share_robustness"])))

    ac = r["alternative_control"]
    ck("alternative control is Part 135 Commuter",
       ac["control_label"] == "Part 135 Commuter", str(ac["control_label"]))
    ck("alternative control has 4 outcome rows",
       len(ac["did_by_outcome"]) == 4, str(len(ac["did_by_outcome"])))
    ck("alt placebo p in [0,1]",
       0.0 <= ac["placebo_p_value"] <= 1.0, str(ac["placebo_p_value"]))

    tv = r["time_varying_share"]
    ck("time-varying share beta finite",
       isinstance(tv["beta_log1p"], (int, float)) and -10 < tv["beta_log1p"] < 10,
       str(tv["beta_log1p"]))

    ck("did_by_outcome has 4 rows",
       len(r["did_by_outcome"]) == 4, str(len(r["did_by_outcome"])))

    dd = r["descriptive"]
    ck("treated_pre_mean > 0",
       dd["treated_pre_mean_monthly"] > 0, str(dd["treated_pre_mean_monthly"]))
    ck("n_panel_observations==192",
       dd["n_panel_observations"] == 192, str(dd["n_panel_observations"]))
    ck("n_months==96", dd["n_months"] == 96, str(dd["n_months"]))
    ck("n_groups==2", dd["n_groups"] == 2, str(dd["n_groups"]))

    nfl = r["narrative_fatigue_flag"]
    ck("at least one fatigue keyword matched",
       nfl["unique_keywords_present"] >= 1, str(nfl))

    pa = r["power_analysis"]
    ck("power_analysis.alt_true_log_effect == -0.05",
       abs(pa["alt_true_log_effect"] - (-0.05)) < 1e-9,
       str(pa["alt_true_log_effect"]))
    ck("power_analysis.alt_placebo_p in [0,1]",
       0.0 <= pa["alt_placebo_p_value"] <= 1.0,
       str(pa["alt_placebo_p_value"]))

    se = r["did_log1p"]["cluster_se"]
    ck("cluster_se > 0", se > 0, str(se))

    ck("bootstrap CI near point estimate",
       (ci[0] - 0.5) <= b <= (ci[1] + 0.5), f"b={b}, ci={ci}")

    # ── NEW: effect-size plausibility (prevents absurd DiD outputs) ──────
    ck("|did_log1p.beta| <= 2.0 (plausibility)",
       abs(b) <= 2.0, f"beta={b}")
    asinh_beta = r["did_by_outcome"]["asinh"]["beta"]
    asinh_se = r["did_by_outcome"]["asinh"]["cluster_se"]
    cohen_d_proxy = abs(asinh_beta) / asinh_se if asinh_se > 1e-9 else float("inf")
    ck("Cohen's d proxy (|asinh beta| / SE) < 5",
       cohen_d_proxy < 5.0, f"d_proxy={cohen_d_proxy:.3f}")

    # ── NEW: CI width sanity (rules out degenerate bootstrap) ────────────
    ci_width = ci[1] - ci[0]
    ck("bootstrap CI width > 0",
       ci_width > 0.0, f"width={ci_width}")
    ck("bootstrap CI width > 1% of max(|beta|, 1e-6) (not degenerate)",
       ci_width > 0.01 * max(abs(b), 1e-6), f"width={ci_width}, beta={b}")

    # ── NEW: placebo distribution must be centered near zero ─────────────
    # (falsification: if this fails, the panel itself has a latent effect
    # that makes the DiD identification suspect)
    ck("placebo distribution centered near 0 (|mean| <= 0.15)",
       abs(p["placebo_mean"]) <= 0.15, f"placebo_mean={p['placebo_mean']}")
    ck("placebo 5th–95th range is strictly positive",
       p["placebo_p95"] > p["placebo_p05"],
       f"p05={p['placebo_p05']} p95={p['placebo_p95']}")
    ck("placebo sd > 0 (null distribution not degenerate)",
       p["placebo_sd"] > 0.0, f"sd={p['placebo_sd']}")

    # ── NEW: sensitivity sweep coherence (findings hold under offsets) ───
    sens_betas = [s["beta_log"] for s in r["sensitivity_offsets"]]
    sens_finite = all(isinstance(x, (int, float)) and abs(x) < 10.0 for x in sens_betas)
    ck("all sensitivity betas finite and |b| < 10",
       sens_finite, str(sens_betas))
    # The observed DiD must itself appear as the "offset=0" row:
    zero_offset_rows = [s for s in r["sensitivity_offsets"] if s["offset_months"] == 0]
    ck("sensitivity_offsets includes offset_months=0",
       len(zero_offset_rows) == 1, str(len(zero_offset_rows)))

    # ── NEW: share-robustness coherence ──────────────────────────────────
    share_betas = [s["beta_log1p"] for s in r["share_robustness"]]
    ck("all share-robustness betas finite and |b| < 5",
       all(isinstance(x, (int, float)) and abs(x) < 5.0 for x in share_betas),
       str(share_betas))
    baseline_rows = [s for s in r["share_robustness"]
                     if abs(s["treated_share"] - 0.86) < 1e-9]
    ck("share_robustness includes baseline 0.86/0.14 row",
       len(baseline_rows) == 1, str(len(baseline_rows)))

    # ── NEW: alt-DGP power analysis direction ────────────────────────────
    # Under true_log_effect = -0.05, the alt-DGP DiD should be negative or
    # at least below the null-DGP DiD in expectation.  (Machine-checkable
    # negative/falsification check: alt-DGP beta should be LESS than null-
    # DGP beta + 0.5 log-units.)
    null_beta = r["did_log1p"]["beta"]
    alt_beta = pa["alt_beta_log1p"]
    ck("alt-DGP beta < null-DGP beta + 0.5 (effect direction sanity)",
       alt_beta < null_beta + 0.5,
       f"alt={alt_beta}, null={null_beta}")

    # ── NEW: config/seed bookkeeping ─────────────────────────────────────
    cfg = r["config"]
    ck("config.seed == 42", cfg["seed"] == 42, str(cfg["seed"]))
    ck("config.n_perm >= 5000", cfg["n_perm"] >= 5000, str(cfg["n_perm"]))
    ck("config.n_boot >= 2000", cfg["n_boot"] >= 2000, str(cfg["n_boot"]))

    # ── NEW: alt control (Part 135) coherence ────────────────────────────
    ck("alt-control commuter_pre_mean_monthly > 0",
       ac["commuter_pre_mean_monthly"] > 0.0,
       str(ac["commuter_pre_mean_monthly"]))
    ck("alt-control commuter_post_mean_monthly > 0",
       ac["commuter_post_mean_monthly"] > 0.0,
       str(ac["commuter_post_mean_monthly"]))

    # ── NEW: descriptive non-degeneracy ──────────────────────────────────
    ck("control_pre_mean_monthly > 0",
       dd["control_pre_mean_monthly"] > 0.0,
       str(dd["control_pre_mean_monthly"]))

    # ── NEW: calibration bookkeeping ─────────────────────────────────────
    cal = r["calibration"]
    ck("calibration.treated_share + control_share == 1.0",
       abs((cal["treated_share"] + cal["control_share"]) - 1.0) < 1e-9,
       f"sum={cal['treated_share'] + cal['control_share']}")

    # ── NEW: config bookkeeping (CI level, significance threshold) ───────
    ck("config.ci_level ~= 0.95",
       abs(cfg.get("ci_level", 0.95) - 0.95) < 1e-9, str(cfg.get("ci_level")))
    ck("config.significance_threshold in (0, 0.5)",
       0.0 < cfg.get("significance_threshold", 0.05) < 0.5,
       str(cfg.get("significance_threshold")))
    ck("config.block_len_months == 3",
       cfg.get("block_len_months") == 3, str(cfg.get("block_len_months")))

    # ── NEW: seasonal shape invariants (mean ~= 1.0, 12 months) ──────────
    ss = cal.get("seasonal_shape", [])
    ck("seasonal_shape has 12 entries", len(ss) == 12, str(len(ss)))
    if ss:
        ss_mean = sum(ss) / len(ss)
        ck("seasonal_shape mean within [0.95, 1.05]",
           0.95 <= ss_mean <= 1.05, f"mean={ss_mean:.4f}")

    # ── NEW: bootstrap CI coverage sanity on log1p outcome ───────────────
    # The bootstrap mean should be reasonably close to the point estimate.
    log1p_cis_midpoint = (ci[0] + ci[1]) / 2.0
    ck("bootstrap CI midpoint within 1.0 log-unit of point estimate",
       abs(log1p_cis_midpoint - b) < 1.0,
       f"mid={log1p_cis_midpoint:.3f}, beta={b:.3f}")

    # ── NEW: share-robustness betas span a non-trivial range ─────────────
    ck("share_robustness betas vary (not all identical)",
       (max(share_betas) - min(share_betas)) > 1e-6 if share_betas else False,
       f"range=[{min(share_betas):.4f},{max(share_betas):.4f}]")

    fail = [c for c in checks if not c[1]]
    ok = [c for c in checks if c[1]]
    print("VERIFICATION:")
    for label, pass_, detail in checks:
        flag = "PASS" if pass_ else "FAIL"
        print(f"  [{flag}] {label}" + (f"  ({detail})" if not pass_ else ""))
    print(f"\n{len(ok)}/{len(checks)} checks passed.")
    if fail:
        print(f"FAILURES: {len(fail)}", file=sys.stderr)
        sys.exit(1)
    print("ALL CHECKS PASSED")


def main():
    ap = argparse.ArgumentParser(description="FAR 117 Fatigue DiD with Cargo and Commuter Placebos")
    ap.add_argument("--verify", action="store_true",
                    help="Re-read results.json and check invariants, then exit.")
    args = ap.parse_args()

    here = os.path.dirname(os.path.abspath(__file__))
    out_dir = here
    cache_dir = os.path.join(here, "cache")
    try:
        os.makedirs(cache_dir, exist_ok=True)
    except OSError as e:
        print(f"ERROR: could not create cache dir {cache_dir}: {e}",
              file=sys.stderr)
        sys.exit(2)

    if args.verify:
        verify(out_dir)
        return

    print("[1/7] Downloading + verifying ASRS publications...")
    try:
        data = load_data(cache_dir)
    except RuntimeError as e:
        print(f"ERROR: data load failed: {e}", file=sys.stderr)
        print("If offline, pre-populate the cache with the 3 expected PDFs.",
              file=sys.stderr)
        sys.exit(2)
    except Exception as e:
        print(f"UNEXPECTED ERROR during data load: {type(e).__name__}: {e}",
              file=sys.stderr)
        sys.exit(3)
    print(f"       OK — 3 PDFs cached at {cache_dir}")

    print("[2/7] Building monthly panels (main, alt-DGP, time-varying, Part 135)...")
    n_panel = len(data["panel"])
    print(f"       OK — main panel: {n_panel} obs")

    print("[3/7] Core DiD across 4 outcome transformations + cluster SE + bootstrap...")
    results = run_analysis(data)
    for name in ("raw", "log1p", "asinh", "log_half"):
        b = results["did_by_outcome"][name]
        print(f"       {name}: beta={b['beta']:.4f}, SE={b['cluster_se']:.4f}")

    print(f"[4/7] Placebo-date permutation test ({results['placebo']['n_permutations_run']} perms)...")
    p = results["placebo"]
    print(f"       OK — two-sided p={p['p_value_two_sided']:.4f}")

    print("[5/7] Share-robustness sweep + time-varying share + alternative Part 135 control...")
    for s in results["share_robustness"]:
        print(f"       split={s['treated_share']:.2f}/{s['control_share']:.2f}: "
              f"beta={s['beta_log1p']:.4f}")
    ac = results["alternative_control"]
    print(f"       Part 135 control: beta(log1p)={ac['did_by_outcome']['log1p']['beta']:.4f}, "
          f"placebo p={ac['placebo_p_value']:.4f}")

    print("[6/7] Sensitivity + narrative-keyword flag + power analysis...")
    pa = results["power_analysis"]
    print(f"       OK — alt-DGP placebo p={pa['alt_placebo_p_value']:.4f}")

    print("[7/7] Writing results.json and report.md...")
    generate_report(results, out_dir)
    print(f"       OK — {os.path.join(out_dir, RESULTS_FILE)}")
    print(f"       OK — {os.path.join(out_dir, REPORT_FILE)}")

    print("ANALYSIS COMPLETE")


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

**Expected output**: No stdout (file written silently).

**Success check**: `ls /tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did/script.py` shows a file of ≈30+ kilobytes.

---

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did && python3 script.py
```

**Expected output** (7 sections ending with `ANALYSIS COMPLETE`):

```
[1/7] Downloading + verifying ASRS publications...
[2/7] Building monthly panels ...
[3/7] Core DiD across 4 outcome transformations ...
[4/7] Placebo-date permutation test (5000 perms)...
[5/7] Share-robustness sweep + time-varying share + alternative Part 135 control...
[6/7] Sensitivity + narrative-keyword flag + power analysis...
[7/7] Writing results.json and report.md...
ANALYSIS COMPLETE
```

**Success criteria**: runs in ≤ 5 minutes, writes `results.json` and `report.md`, ends with `ANALYSIS COMPLETE`.

---

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did && python3 script.py --verify
```

**Expected output**: all checks print `[PASS]`, followed by `N/N checks passed.` (with N ≥ 40), followed by the literal line `ALL CHECKS PASSED`. Exit code is `0`.

**Failure mode**: any `[FAIL]` line, followed by exit 1 and a `FAILURES: k` line on stderr. Inspect the detail, fix the upstream cause, re-run step 3.

---

## Expected Artifacts

```
/tmp/claw4s_auto_far-117-2014-and-fatigue-coded-asrs-incidents-did/
├── script.py
├── results.json
├── report.md
└── cache/
    ├── ASRS_ProgramBriefing.pdf
    ├── acr_fatg.pdf
    └── com_fatigue.pdf
```

## Reproducibility Contract

- `SEED = 42` governs every stochastic step (Poisson draws, bootstrap, permutations, share-robustness panels, alt-DGP panel, Part 135 panel).
- Downloaded PDFs are integrity-checked by SHA256; drift causes abort.
- Analysis uses only the Python 3.8+ standard library.
- Runtime is ≤ 5 minutes on a modern CPU with warm cache.
- Repeated runs must produce byte-identical `results.json` given the same cached PDFs and Python version (verified by re-running and diffing).

## Success Criteria

A run is considered successful when **all** of the following hold (each is machine-checkable in `--verify` mode):

1. `python3 script.py` exits 0 and stdout ends with the literal line `ANALYSIS COMPLETE`.
2. `results.json` and `report.md` are written to the workspace directory.
3. `python3 script.py --verify` exits 0, prints `[PASS]` for every check, and ends with the literal line `ALL CHECKS PASSED`.
4. `results.json` contains `window.T_months == 96`, `treatment_date == "2014-01-01"`, `treatment_month_index == 48`.
5. `results.json` contains bootstrap 95% confidence intervals with `ci[0] <= point_estimate <= ci[1]` (with a small slack).
6. `results.json` contains a placebo permutation test with `n_permutations_run >= 4000` and `0.0 <= p_value_two_sided <= 1.0`.
7. `results.json` reports four outcome transformations (raw, log1p, asinh, log_half), five share-robustness rows, five sensitivity-offset rows, one time-varying-share row, and one alternative-control (Part 135 commuter) block.
8. Effect-size plausibility: `|did_log1p.beta| <= 2.0` and the `asinh` standardized effect size (Cohen's d proxy) has `|d| < 5`.
9. CI width sanity: `bootstrap_ci[1] - bootstrap_ci[0]` is at least 1% of `max(|beta|, 1e-6)`.
10. Falsification/negative control: the placebo-date permutation distribution is centered near zero (mean within ±0.15 log-units of 0) and its 5th–95th percentile range is strictly positive.
11. All three PDFs in `cache/` match their expected SHA256 (pinned by `EXPECTED_SHA256`).

## Failure Conditions

Treat a run as **failed** (do not use its results) if any of the following hold:

1. `ANALYSIS COMPLETE` is not the last non-empty line of stdout.
2. `ALL CHECKS PASSED` is not the last non-empty line of the `--verify` stdout.
3. Any PDF SHA256 mismatches the pinned value (indicates upstream NASA publication update — re-pin and re-validate before trusting).
4. A downstream DiD coefficient is `NaN` or `inf` (indicates a degenerate panel: missing group or missing time index).
5. The placebo permutation returns fewer than 4,000 valid coefficients (too many NaNs = numerical instability).
6. Any SHA256 in `provenance` is not exactly 64 hex characters (a file was silently truncated or overwritten).
7. The panel has fewer than 192 rows (2 groups × 96 months) — indicates data corruption.
8. Network access is unavailable on first run and no cache exists — script fails explicitly with a descriptive error and exit code 2 (not 1).

### Known methodological caveats (what this study does **not** show)

- This is a **calibrated reconstruction** of the monthly Part 121 / Part 121F / Part 135 fatigue-report panel from NASA ASRS published aggregates (annual intake, cumulative fatigue-set totals, seasonal shape, passenger/cargo/commuter share). It is **not** a per-report extraction — ASRS does not expose a bulk public interface at that granularity. All findings are **conditional on** the reconstruction calibration being correct on average.
- ASRS is a **voluntary** reporting system. Reporting propensity may have shifted over 2010–2017 (e.g., ASAP program expansion, cabin-crew awareness campaigns), and a post-period bump in reports could reflect more reporting rather than more fatigue.
- A null DiD coefficient could reflect insufficient statistical power rather than true equivalence. The power analysis under a –5% log-rate alternative quantifies this directly.
- Part 121F cargo operators became subject to FAR 117-equivalent rules in 2024, **after** our 2010–2017 analysis window. Readers should not extrapolate the control-group assumption past 2024.
- The panel DGP uses a homogeneous Poisson within each (group, month) cell; we do not model overdispersion or autocorrelation beyond the block-bootstrap.