---
name: "WUI-expansion vs. fire-behavior decomposition of US wildfire structure losses"
description: "Decomposes the observed 1999-2023 rise in US wildfire structure destruction into an exposure component (WUI housing-unit growth, SILVIS 1990/2000/2010 census series) and a hazard component (losses per WUI housing unit, derived from NIFC annual structure-loss totals). Runs a log-linear shift-share decomposition, a constant-1990-WUI counterfactual, a paired-series permutation test that breaks the exposure-hazard pairing while preserving each series's own temporal structure, a moving-block bootstrap for CIs, and six sensitivity checks (window, mega-fire exclusion, WUI-extrapolation assumption, interface-vs-intermix, residence-only vs total-structures, log-linear vs additive decomposition)."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "wildfire", "wui", "structure-loss", "attribution", "shift-share", "silvis", "nifc", "counterfactual"]
python_version: ">=3.8"
dependencies: []
data_source: "NIFC Wildland Fire Summary and Statistics Annual Reports 1999-2023 (https://www.nifc.gov/fire-information/statistics) for annual structure-loss totals. SILVIS WUI Maps v2 / Radeloff et al. 2018 PNAS (http://silvis.forest.wisc.edu/data/wui-change/) for 1990/2000/2010 housing-unit counts within the US wildland-urban interface."
data_revision: "NIFC year-end figures as published in the 'National Interagency Coordination Center Wildland Fire Summary and Statistics Annual Report' series (2023 edition is the most recent at the time of skill authoring). SILVIS values are the 1990/2000/2010 decennial census overlays reported in Radeloff et al., PNAS 115(13), 2018."
---

# WUI-expansion vs. fire-behavior decomposition of US wildfire structure losses

## Trigger / When to Use This Skill

**Trigger:** Use this skill when you need to test whether an apparent trend in time-series loss data is a genuine hazard signal or a reporting/exposure artifact, using a negative-control (frozen-exposure counterfactual) design paired with a shift-share log decomposition.

**In one sentence.** Given a paired (loss_t, exposure_t) series, this skill quantifies what fraction of log(loss) trend is mechanical exposure growth vs. intensifying per-exposure hazard rate, and tests significance against both a marginal-shuffle null and an autocorrelation-preserving circular-shift null.

**Specifically use this skill when:**

- You have an annual count of losses (fatalities, damaged structures, insurance claims, strikes, outages) that rose over time.
- You have an independently-measured annual or census-year series for the *exposed population* that also changed over time.
- You suspect — or want to rule out — the possibility that "more things to damage" is a non-trivial share of the trend, before attributing the rise to a hazard-driver narrative (climate change, worsening fire weather, ageing infrastructure, etc).

**Do NOT use this skill when:** your loss series and exposure series are not independently measured (e.g. both derived from the same table); your exposure measurement is annual but lumpy (e.g. single census values only at endpoints); your loss counts are dominated by a single outlier year that cannot be identified *a priori*; or you want causal attribution of the residual `β_R` to a specific physical driver (this skill measures the non-exposure residual but does not attribute it).

The skill is designed specifically for the common claim "worsening fire weather and climate change have driven the recent surge in wildfire structure destruction," where the US WUI housing-unit count simultaneously grew by roughly 41% over the period. More generally it fits any loss-count series that factors multiplicatively as `losses = exposure × hazard_rate` with an independently-measured exposure series.

## Research Question

**Primary hypothesis (H1).** Of the log-linear trend in annual US wildfire structure losses 1999-2023, the share explained by WUI housing-unit growth is non-zero and non-trivial.

**Null (H0).** The exposure share of the trend is indistinguishable from what would be produced by random year-pairing between the loss series and the exposure series, once we control for each series's own autocorrelation.

**Formal decomposition.** Let `L_t` be annual structures destroyed, `H_t` WUI housing-unit count, and `R_t ≡ L_t / H_t` the hazard rate. Then `ln(L_t) = ln(H_t) + ln(R_t)`. Fitting separate OLS trends on year gives `β_L = β_H + β_R`. Define the **exposure share of the trend** as `β_H / β_L`.

**Falsifiability.**
- If `β_H / β_L` is statistically indistinguishable from 0 under a circular-shift null that preserves exposure autocorrelation (p > 0.05), H1 is rejected and WUI growth does not meaningfully explain the trend.
- If the constant-1990-WUI counterfactual's avoided-loss share is close to 0, H1 is rejected.
- If the sensitivity suite shows the exposure share flipping sign or swinging by more than 30 percentage points across reasonable window/extrapolation/loss-column choices, the primary result is not robust enough to support H1.

**Method.** Log-linear shift-share OLS + Laspeyres additive cross-check, constant-baseline counterfactual, moving-block bootstrap CIs (block=3, 5000 iter), paired-series marginal-shuffle null (2000 iter, seeded), exact circular-shift null (n-1 shifts), six sensitivity checks (analysis window, mega-fire exclusion, WUI post-2010 extrapolation rule, interface-vs-total WUI, residences-only vs total-structures, leave-one-year-out).

### Preconditions

- **Python version:** 3.8 or newer. **Standard library only.** No numpy, scipy, pandas, or any pip install.
- **Network needs:** Network access is **optional**. The NIFC annual totals and the three SILVIS WUI census-year snapshots are small and are embedded in the script with pinned SHA-256 hashes; lightweight reachability probes against the NIFC and SILVIS source URLs are attempted once and failure is non-fatal.
- **Approximate runtime:** 25-60 seconds end-to-end on a single CPU (includes 5,000-iteration moving-block bootstrap, 2,000-iteration paired-series permutation test, leave-one-out and window sensitivity sweeps, and verification).
- **Disk:** < 20 KB cached data, < 400 KB output.

## Adaptation Guidance

This skill implements a **generic shift-share attribution for a multiplicatively-factorable loss series**. The statistical machinery in `run_analysis` is domain-agnostic; only the `# DOMAIN CONFIGURATION` block in the script references wildfires, housing, and WUI specifics.

What to change in `DOMAIN CONFIGURATION` when adapting:

- `LOSS_DATA_URL`, `LOSS_EMBEDDED_CSV`, `LOSS_EMBEDDED_SHA256` — swap to your annual loss series.
- `EXPOSURE_CENSUS_YEARS`, `EXPOSURE_CENSUS_VALUES`, `EXPOSURE_EMBEDDED_SHA256` — swap to your exposure census snapshots.
- `BASELINE_YEAR` — the year whose exposure level defines the "constant-exposure counterfactual" (1990 here).
- `ANALYSIS_START_YEAR`, `ANALYSIS_END_YEAR` — analysis window.
- `MEGA_EVENT_YEARS` — years to drop in the outlier-robustness check.
- `SENSITIVITY_EXTRAPOLATION_MODES` — list of alternative extrapolation assumptions beyond 2010.
- `DOMAIN_LABEL`, `LOSS_LABEL`, `EXPOSURE_LABEL` — cosmetic, used in `generate_report`.

What stays the same (do not edit):

- `ols_fit`, `matrix_inverse` — hand-rolled OLS used for ≤2 predictors.
- `interp_exposure` — piecewise-linear interpolation between census years with a pluggable post-last-census extrapolation rule.
- `shift_share_log_decomposition` — β_L = β_H + β_R identity plus share computation.
- `counterfactual_total_loss` — computes Σ H_base · R_t and the implied avoided-loss share.
- `paired_permutation_test` — permutes the year pairing between the exposure and hazard-rate series while leaving each series's own autocorrelation intact.
- `moving_block_bootstrap` — autocorrelation-robust CI for the attribution share.
- `run_analysis` — orchestrates all of the above generically.

Example domain swaps: coastal-flood building losses vs. coastal housing stock (NOAA SLR viewer + ACS); heat-related mortality vs. population over 65 in heat-exposed counties (CDC WONDER + Census); aviation bird-strike damage vs. commercial flight operations (FAA Wildlife Strike DB + BTS T-100). Each reduces to "how much of the trend is the exposure base vs. the rate of damage per exposed unit?"

### Concrete domain-swap worked example

To port this skill to **coastal-flood building losses vs. US coastal-county housing stock (2000-2023)**:

1. Replace `LOSS_EMBEDDED_CSV` with year, buildings_damaged, homes_only from NOAA National Flood Insurance Program annual reports. Recompute `LOSS_EMBEDDED_SHA256` with `hashlib.sha256(csv.encode()).hexdigest()`.
2. Replace `EXPOSURE_CENSUS_YEARS = [2000, 2010, 2020]` and `EXPOSURE_CENSUS_VALUES = [coastal_housing_2000, ..., coastal_housing_2020]` from ACS coastal-county roll-ups; recompute `EXPOSURE_EMBEDDED_SHA256`.
3. Set `ANALYSIS_START_YEAR = 2000`, `ANALYSIS_END_YEAR = 2023`, `BASELINE_YEAR = 2000`.
4. Set `MEGA_EVENT_YEARS = [2005, 2017]` (Katrina, Harvey-Irma-Maria).
5. Set `DOMAIN_LABEL = "US coastal flood building losses"`, `LOSS_LABEL = "buildings damaged"`, `EXPOSURE_LABEL = "coastal housing units"`.
6. Leave `SENSITIVITY_EXTRAPOLATION_MODES` alone unless census-year gap changes; then adjust `linear_2000_2010` keys accordingly.
7. Everything under `# STATISTICAL HELPERS` and `run_analysis` is unchanged.
8. Update `verify()` assertion [14] (WUI-growth threshold) to a threshold appropriate to coastal-housing growth (≈1.15x rather than 1.2x).

## Success Criteria

The analysis is considered to have run correctly when ALL of the following hold. They are auto-checked by `analyze.py --verify` (30 assertions) and by structural checks on `results.json`.

1. **Execution completes.** `analyze.py` exits 0 and the final stdout line is literally `ANALYSIS COMPLETE`.
2. **Verification completes.** `analyze.py --verify` exits 0 and stdout contains `VERIFICATION PASSED`.
3. **Data integrity.** Both SHA-256 hashes in `meta` equal the pinned `LOSS_EMBEDDED_SHA256` and `EXPOSURE_EMBEDDED_SHA256` (not `BOOTSTRAP`).
4. **Shift-share identity holds.** `|β_L − (β_H + β_R)| < 1e-8` and `exposure_share + hazard_share = 1.0` to machine precision.
5. **Coverage.** `n_years ≥ 20`, window spans 1999-2023 (or adapted equivalent).
6. **Bootstrap depth.** ≥4000 valid bootstrap samples and ≥1800 valid marginal-shuffle samples; exact circular-shift null has ≥20 non-zero shifts.
7. **Plausibility.** Exposure share ∈ [−2, 2]; avoided-loss share ∈ (−1, 1); lag-1 residual ACF ∈ [−1, 1]; WUI growth multiplier > 1.2x.
8. **Sensitivity coverage.** All 6 sensitivity checks ran (window, mega-exclude, extrapolation, interface-only, residences-only, leave-one-out) with non-empty output; `flat_2010` exposure share ≤ `linear_2000_2010` exposure share (monotonicity).

## Failure Conditions

This analysis is **not** trustworthy as produced if any of the following are true. Each corresponds to a specific check the user should perform before citing the numbers:

1. **Post-2010 WUI is extrapolated, not observed.** The SILVIS WUI census snapshots are 1990/2000/2010 only. Values for 2011-2023 come from linear extrapolation of the 2000-2010 slope and may overstate or understate true WUI growth. Because of this construction, the `β_H` OLS trend against years has R² ≈ 0.998 — this is an artifact of the interpolation, not a real statistical coincidence. The `wui_extrapolation_mode` sensitivity shows how much the exposure share changes under each alternative. *Check:* the `flat_2010` row of `sensitivity.wui_extrapolation_mode` in `results.json` sets a strict lower bound on the exposure share; if that lower bound is negative or near-zero, the attribution is not robust.
2. **Reported hazard-rate trend is not causal.** `β_R` captures *any* non-exposure-linked change: climate, fuel-load, ignition density, defensible-space adoption, suppression budgets, building-code stringency, or reporting-completeness changes at NIFC. The skill quantifies how much of `β_L` is NOT mechanically due to `H_t` — it does not attribute the residue to a physical driver. Any claim "the residual is caused by climate change" is out of scope.
3. **NIFC totals are operational, not audited.** NIFC's "structures destroyed" series is an operational year-end tally. It mixes residences, commercial buildings, outbuildings and minor structures; the residence/non-residence split changed reporting-standard at least once during the window. Single-year values are noisy by ±10-20%. *Check:* the `residences_only_losses` sensitivity isolates the more consistently-reported subset and should move the shares by less than ~15 points if reporting drift is not dominating the trend.
4. **The shift-share identity is a mechanical identity, not evidence of causation.** `β_L = β_H + β_R` is an algebraic fact, not a test of whether WUI-growth *caused* loss-growth. The null tests (marginal-shuffle + circular-shift) are the actual significance tests for whether the cross-series pairing is non-random. If both null p-values are > 0.05, the paired pattern is not statistically distinguishable from a random pairing under the corresponding null, and the decomposition should be reported as descriptive only.
5. **Three-point exposure series.** With only 3 directly-measured SILVIS census values, the OLS fit to log(H_t) against year is near-perfectly linear by construction — the reported `r2_log_exposure_trend` close to 1.0 is an artifact of the interpolation scheme, NOT evidence of a real mechanistic relationship. Do not report this R² as a finding.
6. **Out-of-scope.** This analysis does NOT: (a) attribute `β_R` to any physical driver; (b) separate interface-WUI structure losses from intermix-WUI; (c) use fire-perimeter microdata (MTBS) or parcel-level structure age; (d) compute dollar losses (counts only); (e) infer causality about specific policy interventions.

## Overview

**Research question.** NIFC records show US wildfire structure losses rose from roughly 1,000 structures/year in the early 2000s to a multi-year mean near 10,000 structures/year in 2017-2023. The dominant popular narrative attributes this rise to climate change and worsening fire weather. The SILVIS Wildland-Urban Interface (WUI) dataset shows the housing-unit population that physically sits inside the WUI grew from 30.8 million in 1990 to 43.4 million in 2010, and (by linear extrapolation) to ~49-51 million by 2020. The question this skill quantifies is: **of the observed rise in annual structure losses 1999-2023, what fraction is explained by the mechanical growth of the exposed housing-unit base in the WUI, and what fraction is explained by an increased hazard rate (losses per exposed housing unit)?**

**Approach — shift-share.** Let `L_t` be NIFC structures destroyed in year `t`, `H_t` the WUI housing-unit count in year `t`, and `R_t ≡ L_t / H_t` the hazard rate (losses per WUI housing unit). Then exactly

    ln(L_t) = ln(H_t) + ln(R_t)

Fitting separate OLS trends on year to each of these three log-series yields slopes that satisfy `β_L = β_H + β_R` by construction. We define the **exposure share of the trend** as `β_H / β_L` and the **hazard share** as `β_R / β_L`. These sum to 1 exactly.

**Approach — counterfactual.** We additionally compute a dollar-identifiable counterfactual: **what would 1999-2023 cumulative structure losses have been if WUI housing had stayed frozen at the 1990 level?**

    L_t^CF = H_1990 · R_t
    avoided-loss share = 1 − Σ_t L_t^CF / Σ_t L_t

**Null model.** Two complementary nulls are reported.

*Marginal-shuffle (Fisher-randomization) null.* Shuffle the year-index of the exposure series iid over 2,000 iterations. This preserves the marginal distribution of `H_t` but deliberately destroys its time order. It is a test of whether the observed cross-series pairing is non-random.

*Circular-shift null (autocorrelation-preserving).* Enumerate all `n − 1` nonzero circular shifts `H_t → H_{(t+k) mod n}` and recompute the exposure share at each shift. A circular shift preserves the autocorrelation structure and trend shape of the exposure series exactly; only the cross-series year-pairing is broken (with a wrap-around discontinuity). Because shift enumeration is exact, the minimum attainable two-sided p is `1 / (n − 1) ≈ 0.042` at `n = 25`; this discreteness is reported honestly.

**Rigor.** 5,000-iteration moving-block bootstrap (block length 3) for all trend slopes, share statistics, and counterfactual avoided-loss fraction; 2,000-iteration marginal-shuffle null plus exact `n−1`-shift circular-shift null; six sensitivity analyses (analysis window, mega-fire-year exclusion, WUI post-2010 extrapolation assumption, interface-only WUI, residence-only losses, leave-one-year-out).

**What this is not.** This is an **exposure-vs-rate attribution**, not a causal decomposition of the rate. It does not attribute `β_R` to any specific physical driver (climate, fuel accumulation, suppression policy, ignition patterns). It does not use fire-perimeter microdata (MTBS) or building-level structure-age data. The 1990-2020 WUI series combines three directly-measured SILVIS census snapshots (1990, 2000, 2010) with a linear extrapolation for 2011-2023; the main analysis uses this default extrapolation and the sensitivity suite reports how shares move under alternatives.

---

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_wui-expansion-explains-more-structure-loss-variance-than-fir/cache
```

**Expected output:** No output (directory created silently).

**Success condition:** `/tmp/claw4s_auto_wui-expansion-explains-more-structure-loss-variance-than-fir/cache` exists and is writable.

---

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_wui-expansion-explains-more-structure-loss-variance-than-fir/analyze.py
#!/usr/bin/env python3
"""
WUI-expansion vs. fire-behavior decomposition of US wildfire structure losses.

Fits a log-linear shift-share decomposition to NIFC annual structures-destroyed
1999-2023 paired with SILVIS WUI housing-unit counts (1990/2000/2010 census
snapshots linearly interpolated/extrapolated). Decomposes the total trend into
an exposure component (β_H / β_L) and a hazard component (β_R / β_L). Runs a
constant-1990-WUI counterfactual, a paired-series permutation test null, a
moving-block bootstrap for CIs, and a full sensitivity suite. Standard library
only (Python 3.8+).
"""

import argparse
import hashlib
import json
import math
import os
import random
import sys
import urllib.request
import urllib.error


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

DOMAIN_LABEL = "US wildfire structure losses"
LOSS_LABEL = "structures destroyed"
EXPOSURE_LABEL = "WUI housing units"

LOSS_DATA_URL = "https://www.nifc.gov/fire-information/statistics"
EXPOSURE_DATA_URL = "http://silvis.forest.wisc.edu/data/wui-change/"

# NIFC annual total structures destroyed (residences + minor structures +
# commercial), 1999-2023. Source: NIFC "National Interagency Coordination
# Center Wildland Fire Summary and Statistics Annual Report" series.
# Column: year, structures, residences_only
# residences_only columns are NIFC residences-destroyed (subset of structures)
# used in the residence-only sensitivity check. Where NIFC published only
# a total, residences_only is reported as the NIFC "residences" subtotal.
LOSS_EMBEDDED_CSV = """year,structures,residences_only
1999,748,394
2000,861,484
2001,731,421
2002,2381,835
2003,5781,3710
2004,1241,547
2005,967,433
2006,1150,562
2007,2936,2356
2008,2611,1497
2009,1414,598
2010,1069,450
2011,5484,2904
2012,4244,2216
2013,1151,556
2014,1576,1108
2015,4636,2620
2016,4313,1928
2017,12306,6973
2018,25790,18137
2019,963,536
2020,17904,10488
2021,5972,3577
2022,2717,1420
2023,4369,2221
"""

# SHA-256 of LOSS_EMBEDDED_CSV as UTF-8 bytes.
LOSS_EMBEDDED_SHA256 = "3e491c7329063cf0a4e63155d5a68cd4657e5f652c3ab56546407f4d206ea804"
# SHA-256 of canonical JSON serialization of the SILVIS census object.
EXPOSURE_EMBEDDED_SHA256 = "eba054d411731c1bada92dd3fb5eb78445bf9ee440cc1da433d11ceba4135e92"

# SILVIS WUI census-year housing-unit counts (1990/2000/2010).
# Source: Radeloff et al. 2018 PNAS Table 1 and SI, "Rapid growth of the US
# wildland-urban interface raises wildfire risk," PNAS 115(13):3314-3319.
# All figures are in total WUI (interface + intermix).
EXPOSURE_CENSUS_YEARS = [1990, 2000, 2010]
EXPOSURE_CENSUS_VALUES = [30826000, 37296000, 43438000]  # housing units
# SILVIS WUI area (km^2) in the same snapshots, used as a cross-check only.
EXPOSURE_CENSUS_AREA_KM2 = [581167, 672000, 770437]

# Interface-only WUI housing units (sensitivity). Interface = WUI adjacent to
# wildland; intermix = WUI interspersed with wildland. From Radeloff et al. 2018
# Table 1. Approximate split 30% interface, 70% intermix.
EXPOSURE_INTERFACE_ONLY = [9300000, 11400000, 13200000]

# Post-2010 extrapolation for the default analysis and sensitivity alternatives.
# "linear_2000_2010": slope = (H_2010 - H_2000) / 10, extrapolated forward.
# "linear_1990_2010": slope = (H_2010 - H_1990) / 20, extrapolated forward.
# "flat_2010":        H_t = H_2010 for t > 2010 (conservative lower bound).
SENSITIVITY_EXTRAPOLATION_MODES = ["linear_2000_2010", "linear_1990_2010", "flat_2010"]
DEFAULT_EXTRAPOLATION_MODE = "linear_2000_2010"

ANALYSIS_START_YEAR = 1999
ANALYSIS_END_YEAR = 2023
BASELINE_YEAR = 1990  # counterfactual "frozen WUI" base year
MEGA_EVENT_YEARS = [2017, 2018, 2020]  # dropped in mega-fire sensitivity
ALTERNATIVE_START_YEARS = [2005, 2010]  # shorter-window sensitivities

RANDOM_SEED = 42
BOOTSTRAP_ITERATIONS = 5000
BLOCK_BOOTSTRAP_LENGTH = 3
PERMUTATION_ITERATIONS = 2000

# Workspace can be overridden via WUI_WORKSPACE env var for portability; the
# default is the auto-run path this skill instructs the agent to create.
WORKSPACE = os.environ.get(
    "WUI_WORKSPACE",
    "/tmp/claw4s_auto_wui-expansion-explains-more-structure-loss-variance-than-fir",
)
CACHE_DIR = os.path.join(WORKSPACE, "cache")
RESULTS_JSON = os.path.join(WORKSPACE, "results.json")
REPORT_MD = os.path.join(WORKSPACE, "report.md")


# ═══════════════════════════════════════════════════════════════
# STATISTICAL HELPERS (domain-agnostic)
# ═══════════════════════════════════════════════════════════════


def matrix_inverse(A):
    """Invert an n-by-n matrix via Gauss-Jordan elimination. Pure Python."""
    n = len(A)
    M = [row[:] + [1.0 if i == j else 0.0 for j in range(n)] for i, row in enumerate(A)]
    for i in range(n):
        pivot_row = i
        max_abs = abs(M[i][i])
        for r in range(i + 1, n):
            if abs(M[r][i]) > max_abs:
                max_abs = abs(M[r][i])
                pivot_row = r
        if max_abs < 1e-12:
            raise ValueError("Singular matrix in matrix_inverse")
        M[i], M[pivot_row] = M[pivot_row], M[i]
        pivot = M[i][i]
        M[i] = [v / pivot for v in M[i]]
        for r in range(n):
            if r != i and M[r][i] != 0.0:
                factor = M[r][i]
                M[r] = [M[r][c] - factor * M[i][c] for c in range(2 * n)]
    return [row[n:] for row in M]


def matmul(A, B):
    n, m = len(A), len(A[0])
    p = len(B[0])
    return [[sum(A[i][k] * B[k][j] for k in range(m)) for j in range(p)] for i in range(n)]


def matvec(A, v):
    return [sum(a * b for a, b in zip(row, v)) for row in A]


def transpose(A):
    return [list(col) for col in zip(*A)]


def ols_fit(X, y):
    """OLS: returns (beta, residuals, se_beta, r_squared, yhat)."""
    n, p = len(X), len(X[0])
    Xt = transpose(X)
    XtX = matmul(Xt, X)
    XtX_inv = matrix_inverse(XtX)
    Xty = matvec(Xt, y)
    beta = matvec(XtX_inv, Xty)
    yhat = [sum(X[i][k] * beta[k] for k in range(p)) for i in range(n)]
    resid = [y[i] - yhat[i] for i in range(n)]
    dof = n - p
    rss = sum(r * r for r in resid)
    sigma2 = rss / dof if dof > 0 else float("nan")
    se_beta = [math.sqrt(max(0.0, sigma2 * XtX_inv[k][k])) for k in range(p)]
    y_mean = sum(y) / n
    tss = sum((yi - y_mean) ** 2 for yi in y)
    r2 = 1.0 - rss / tss if tss > 0 else float("nan")
    return beta, resid, se_beta, r2, yhat


def trend_design(years):
    """Design matrix for y = β0 + β1·t (simple linear trend)."""
    return [[1.0, float(t)] for t in years]


def fit_trend(years, y):
    """Fit y ~ 1 + t and return (intercept, slope, r_squared, se_slope)."""
    X = trend_design(years)
    beta, _, se, r2, _ = ols_fit(X, y)
    return beta[0], beta[1], r2, se[1]


def pct(lst, p):
    """Linear-interpolation percentile (R-7 / numpy default)."""
    lst_sorted = sorted(lst)
    n = len(lst_sorted)
    if n == 0:
        return float("nan")
    if n == 1:
        return lst_sorted[0]
    idx = p * (n - 1)
    lo = int(idx)
    hi = min(lo + 1, n - 1)
    frac = idx - lo
    return lst_sorted[lo] * (1 - frac) + lst_sorted[hi] * frac


def pearson_r(a, b):
    n = len(a)
    ma = sum(a) / n
    mb = sum(b) / n
    num = sum((ai - ma) * (bi - mb) for ai, bi in zip(a, b))
    denom = math.sqrt(sum((ai - ma) ** 2 for ai in a) * sum((bi - mb) ** 2 for bi in b))
    return num / denom if denom > 0 else float("nan")


def interp_exposure(year, census_years, census_values, extrapolation_mode):
    """Piecewise-linear interpolation between census years, with a chosen
    extrapolation rule for years beyond the last census year."""
    ys = census_years
    vs = census_values
    if year <= ys[0]:
        slope = (vs[1] - vs[0]) / (ys[1] - ys[0])
        return vs[0] + slope * (year - ys[0])
    for i in range(len(ys) - 1):
        if ys[i] <= year <= ys[i + 1]:
            slope = (vs[i + 1] - vs[i]) / (ys[i + 1] - ys[i])
            return vs[i] + slope * (year - ys[i])
    if extrapolation_mode == "linear_2000_2010":
        slope = (vs[-1] - vs[-2]) / (ys[-1] - ys[-2])
        return vs[-1] + slope * (year - ys[-1])
    elif extrapolation_mode == "linear_1990_2010":
        slope = (vs[-1] - vs[0]) / (ys[-1] - ys[0])
        return vs[-1] + slope * (year - ys[-1])
    elif extrapolation_mode == "flat_2010":
        return vs[-1]
    else:
        raise ValueError(f"Unknown extrapolation_mode: {extrapolation_mode}")


def shift_share_log_decomposition(years, loss, exposure):
    """Fit ln(L), ln(H), ln(R) vs year, return β_L, β_H, β_R and shares.

    By construction β_L = β_H + β_R. The exposure share is β_H / β_L and
    the hazard share is β_R / β_L; they sum to 1 exactly.
    """
    ln_L = [math.log(l) for l in loss]
    ln_H = [math.log(h) for h in exposure]
    ln_R = [math.log(l / h) for l, h in zip(loss, exposure)]
    _, b_L, r2_L, se_L = fit_trend(years, ln_L)
    _, b_H, r2_H, se_H = fit_trend(years, ln_H)
    _, b_R, r2_R, se_R = fit_trend(years, ln_R)
    if abs(b_L) < 1e-12:
        exposure_share = float("nan")
        hazard_share = float("nan")
    else:
        exposure_share = b_H / b_L
        hazard_share = b_R / b_L
    return {
        "beta_log_loss_per_year": b_L,
        "beta_log_exposure_per_year": b_H,
        "beta_log_hazard_per_year": b_R,
        "se_beta_log_loss": se_L,
        "se_beta_log_exposure": se_H,
        "se_beta_log_hazard": se_R,
        "r2_log_loss_trend": r2_L,
        "r2_log_exposure_trend": r2_H,
        "r2_log_hazard_trend": r2_R,
        "exposure_share_of_trend": exposure_share,
        "hazard_share_of_trend": hazard_share,
        "identity_residual_abs": abs(b_L - (b_H + b_R)),
    }


def counterfactual_total_loss(years, loss, exposure, baseline_exposure):
    """Compute Σ L_t^CF = Σ H_base · R_t and the avoided-loss share.

    The avoided-loss share is (Σ L_t − Σ L_t^CF) / Σ L_t; i.e. the fraction of
    cumulative observed losses that would not have occurred if exposure had
    stayed at the baseline level.
    """
    R = [l / h for l, h in zip(loss, exposure)]
    L_cf = [baseline_exposure * r for r in R]
    sum_obs = sum(loss)
    sum_cf = sum(L_cf)
    avoided = sum_obs - sum_cf
    share = avoided / sum_obs if sum_obs > 0 else float("nan")
    return {
        "baseline_exposure": baseline_exposure,
        "baseline_year": BASELINE_YEAR,
        "sum_observed": sum_obs,
        "sum_counterfactual": sum_cf,
        "avoided_losses_total": avoided,
        "avoided_loss_share": share,
        "per_year_counterfactual": [
            {"year": y, "L_observed": l, "L_counterfactual": lc, "R": r}
            for y, l, lc, r in zip(years, loss, L_cf, R)
        ],
    }


def moving_block_bootstrap(years, loss, exposure, baseline_exposure,
                           block_len, rng, iters):
    """Moving-block bootstrap on paired (loss, exposure) series.

    For each iteration: resample length-n by concatenating randomly-chosen
    blocks of contiguous (year, loss, exposure) triples. Refit the
    shift-share and counterfactual and record the four summary statistics.
    """
    n = len(years)
    max_start = n - block_len
    if max_start < 0:
        return None
    starts = list(range(0, max_start + 1))

    exposure_shares = []
    hazard_shares = []
    avoided_shares = []
    beta_L_samples = []
    beta_H_samples = []
    beta_R_samples = []

    for _ in range(iters):
        sample_yrs = []
        sample_loss = []
        sample_exp = []
        while len(sample_yrs) < n:
            s = rng.choice(starts)
            for k in range(block_len):
                if len(sample_yrs) >= n:
                    break
                sample_yrs.append(years[s + k])
                sample_loss.append(loss[s + k])
                sample_exp.append(exposure[s + k])
        try:
            ss = shift_share_log_decomposition(sample_yrs, sample_loss, sample_exp)
            exposure_shares.append(ss["exposure_share_of_trend"])
            hazard_shares.append(ss["hazard_share_of_trend"])
            beta_L_samples.append(ss["beta_log_loss_per_year"])
            beta_H_samples.append(ss["beta_log_exposure_per_year"])
            beta_R_samples.append(ss["beta_log_hazard_per_year"])
            cf = counterfactual_total_loss(
                sample_yrs, sample_loss, sample_exp, baseline_exposure
            )
            avoided_shares.append(cf["avoided_loss_share"])
        except (ValueError, ZeroDivisionError):
            continue

    return {
        "block_length": block_len,
        "iterations": len(exposure_shares),
        "exposure_share_ci": [pct(exposure_shares, 0.025), pct(exposure_shares, 0.975)],
        "hazard_share_ci": [pct(hazard_shares, 0.025), pct(hazard_shares, 0.975)],
        "avoided_loss_share_ci": [pct(avoided_shares, 0.025), pct(avoided_shares, 0.975)],
        "beta_log_loss_ci": [pct(beta_L_samples, 0.025), pct(beta_L_samples, 0.975)],
        "beta_log_exposure_ci": [pct(beta_H_samples, 0.025), pct(beta_H_samples, 0.975)],
        "beta_log_hazard_ci": [pct(beta_R_samples, 0.025), pct(beta_R_samples, 0.975)],
    }


def marginal_shuffle_null(years, loss, exposure, rng, iters):
    """Marginal-shuffle (Fisher randomization) null for the exposure share.

    Under H0, the cross-series year-pairing between exposure and loss is
    arbitrary. We iid-shuffle the exposure series's year-index while holding
    (year, loss) fixed and recompute the exposure share. This preserves the
    *marginal distribution* of exposure but deliberately destroys its time
    order; it is a test of whether the observed cross-series pairing is
    non-random, not a test that respects exposure's autocorrelation.

    Two-sided p = fraction of iterations with |exposure_share_perm| ≥
    |exposure_share_observed|.
    """
    observed = shift_share_log_decomposition(years, loss, exposure)
    observed_share = observed["exposure_share_of_trend"]
    observed_abs = abs(observed_share)

    null_shares = []
    n = len(years)
    indices = list(range(n))
    for _ in range(iters):
        shuffled = indices[:]
        rng.shuffle(shuffled)
        permuted_exposure = [exposure[i] for i in shuffled]
        try:
            ss = shift_share_log_decomposition(years, loss, permuted_exposure)
            null_shares.append(ss["exposure_share_of_trend"])
        except (ValueError, ZeroDivisionError):
            continue

    if not null_shares:
        return None
    n_ge = sum(1 for s in null_shares if abs(s) >= observed_abs)
    p_value = n_ge / len(null_shares)
    return {
        "iterations": len(null_shares),
        "observed_exposure_share": observed_share,
        "null_exposure_share_mean": sum(null_shares) / len(null_shares),
        "null_exposure_share_sd": math.sqrt(
            sum((s - sum(null_shares) / len(null_shares)) ** 2 for s in null_shares)
            / (len(null_shares) - 1)
        ) if len(null_shares) > 1 else float("nan"),
        "null_exposure_share_p025": pct(null_shares, 0.025),
        "null_exposure_share_p975": pct(null_shares, 0.975),
        "p_two_sided": p_value,
    }


def circular_shift_null(years, loss, exposure):
    """Circular-shift null for the exposure share (autocorrelation-preserving).

    Under H0, we enumerate all n-1 nonzero circular shifts of the exposure
    series relative to the (year, loss) series and recompute the exposure
    share for each shift. A circular shift preserves exposure's own
    autocorrelation and trend shape perfectly; it breaks only the
    cross-series year-pairing, introducing a wrap-around discontinuity at
    the shift boundary.

    This is a stronger null than iid-shuffling because the shifted series
    still has a trend, still has its autocorrelation, and only the pairing
    to loss is misaligned. The exact enumeration yields a coarse p-value
    grid (minimum attainable p = 1/(n-1)).

    Two-sided p = fraction of shifts with |exposure_share_shifted| ≥
    |exposure_share_observed|.
    """
    n = len(years)
    observed = shift_share_log_decomposition(years, loss, exposure)
    observed_share = observed["exposure_share_of_trend"]
    observed_abs = abs(observed_share)

    shares = []
    for k in range(1, n):
        shifted = exposure[k:] + exposure[:k]
        try:
            ss = shift_share_log_decomposition(years, loss, shifted)
            shares.append(ss["exposure_share_of_trend"])
        except (ValueError, ZeroDivisionError):
            continue

    if not shares:
        return None
    n_ge = sum(1 for s in shares if abs(s) >= observed_abs)
    p_value = n_ge / len(shares)
    return {
        "shifts_evaluated": len(shares),
        "observed_exposure_share": observed_share,
        "null_exposure_share_mean": sum(shares) / len(shares),
        "null_exposure_share_p025": pct(shares, 0.025),
        "null_exposure_share_p975": pct(shares, 0.975),
        "p_two_sided": p_value,
        "minimum_attainable_p": 1.0 / len(shares),
    }


def residual_acf_lag1(residuals):
    n = len(residuals)
    mu = sum(residuals) / n
    num = sum((residuals[i] - mu) * (residuals[i - 1] - mu) for i in range(1, n))
    den = sum((r - mu) ** 2 for r in residuals)
    return num / den if den > 0 else float("nan")


# ═══════════════════════════════════════════════════════════════
# DATA LOADING
# ═══════════════════════════════════════════════════════════════


def _sha256_of_string(s):
    return hashlib.sha256(s.encode("utf-8")).hexdigest()


def verify_and_cache_data():
    """Write embedded data to cache, verify SHA-256, probe source URLs."""
    os.makedirs(CACHE_DIR, exist_ok=True)
    loss_path = os.path.join(CACHE_DIR, "nifc_structures.csv")
    with open(loss_path, "wb") as f:
        f.write(LOSS_EMBEDDED_CSV.encode("utf-8"))
    loss_sha = _sha256_of_string(LOSS_EMBEDDED_CSV)

    exposure_obj = {
        "years": EXPOSURE_CENSUS_YEARS,
        "housing_units": EXPOSURE_CENSUS_VALUES,
        "area_km2": EXPOSURE_CENSUS_AREA_KM2,
        "interface_only": EXPOSURE_INTERFACE_ONLY,
    }
    exposure_json = json.dumps(exposure_obj, sort_keys=True)
    exposure_path = os.path.join(CACHE_DIR, "silvis_wui_census.json")
    with open(exposure_path, "w") as f:
        f.write(exposure_json)
    exposure_sha = _sha256_of_string(exposure_json)

    sha_marker = os.path.join(CACHE_DIR, ".data_sha256.txt")
    with open(sha_marker, "w") as f:
        f.write(f"loss_csv  {loss_sha}\n")
        f.write(f"exposure  {exposure_sha}\n")

    if LOSS_EMBEDDED_SHA256 == "BOOTSTRAP":
        print(f"  [note] LOSS_EMBEDDED_SHA256 is BOOTSTRAP; computed: {loss_sha}")
        print(f"  [note] EXPOSURE_EMBEDDED_SHA256 computed: {exposure_sha}")
    else:
        if loss_sha != LOSS_EMBEDDED_SHA256:
            raise ValueError(
                f"Loss data integrity check failed. "
                f"Expected {LOSS_EMBEDDED_SHA256}, got {loss_sha}."
            )
        if exposure_sha != EXPOSURE_EMBEDDED_SHA256:
            raise ValueError(
                f"Exposure data integrity check failed. "
                f"Expected {EXPOSURE_EMBEDDED_SHA256}, got {exposure_sha}."
            )

    def _reach(url):
        try:
            req = urllib.request.Request(url, headers={"User-Agent": "claw4s-skill/1.0"})
            with urllib.request.urlopen(req, timeout=3) as resp:
                return resp.status == 200
        except (urllib.error.URLError, OSError, TimeoutError):
            return False

    return {
        "loss_path": loss_path,
        "exposure_path": exposure_path,
        "loss_sha256": loss_sha,
        "exposure_sha256": exposure_sha,
        "nifc_url_reachable": _reach(LOSS_DATA_URL),
        "silvis_url_reachable": _reach(EXPOSURE_DATA_URL),
    }


def load_data():
    """Load NIFC + SILVIS data into a clean structure."""
    meta = verify_and_cache_data()
    rows = []
    with open(meta["loss_path"], "r") as f:
        header = f.readline().strip().split(",")
        for line in f:
            if not line.strip():
                continue
            parts = line.strip().split(",")
            rows.append({
                "year": int(parts[0]),
                "structures": int(parts[1]),
                "residences_only": int(parts[2]),
            })
    rows.sort(key=lambda r: r["year"])
    with open(meta["exposure_path"], "r") as f:
        exposure_obj = json.loads(f.read())
    return {
        "loss_rows": rows,
        "exposure": exposure_obj,
        "loss_sha256": meta["loss_sha256"],
        "exposure_sha256": meta["exposure_sha256"],
        "nifc_url_reachable": meta["nifc_url_reachable"],
        "silvis_url_reachable": meta["silvis_url_reachable"],
    }


# ═══════════════════════════════════════════════════════════════
# ANALYSIS
# ═══════════════════════════════════════════════════════════════


def build_series(data, start_year, end_year, loss_column="structures",
                 exposure_variant="total", extrapolation_mode=DEFAULT_EXTRAPOLATION_MODE):
    """Assemble year-aligned (year, loss, exposure) series."""
    rows = [r for r in data["loss_rows"] if start_year <= r["year"] <= end_year]
    rows.sort(key=lambda r: r["year"])
    years = [r["year"] for r in rows]
    loss = [float(r[loss_column]) for r in rows]
    census_years = data["exposure"]["years"]
    if exposure_variant == "total":
        census_values = data["exposure"]["housing_units"]
    elif exposure_variant == "interface_only":
        census_values = data["exposure"]["interface_only"]
    else:
        raise ValueError(f"Unknown exposure_variant: {exposure_variant}")
    exposure = [interp_exposure(y, census_years, census_values, extrapolation_mode)
                for y in years]
    return years, loss, exposure


def run_decomposition_on(years, loss, exposure, baseline_exposure, rng,
                         do_bootstrap=True, do_permutation=True):
    """Run the full decomposition (main + null + CIs) on a single series."""
    ss = shift_share_log_decomposition(years, loss, exposure)
    cf = counterfactual_total_loss(years, loss, exposure, baseline_exposure)

    bb = None
    if do_bootstrap:
        bb = moving_block_bootstrap(
            years, loss, exposure, baseline_exposure,
            BLOCK_BOOTSTRAP_LENGTH, rng, BOOTSTRAP_ITERATIONS,
        )
    perm_marginal = None
    perm_circular = None
    if do_permutation:
        perm_marginal = marginal_shuffle_null(years, loss, exposure, rng, PERMUTATION_ITERATIONS)
        perm_circular = circular_shift_null(years, loss, exposure)

    # Residuals of the ln(L) trend -> lag-1 ACF diagnostic
    ln_L = [math.log(l) for l in loss]
    X = trend_design(years)
    _, resid, _, _, _ = ols_fit(X, ln_L)
    acf1 = residual_acf_lag1(resid)

    return {
        "n": len(years),
        "year_range": [years[0], years[-1]],
        "shift_share": ss,
        "counterfactual": cf,
        "moving_block_bootstrap": bb,
        "marginal_shuffle_null": perm_marginal,
        "circular_shift_null": perm_circular,
        "residual_acf_lag1_log_loss": acf1,
        "pearson_corr_logL_logH": pearson_r(
            [math.log(l) for l in loss], [math.log(h) for h in exposure]
        ),
    }


def additive_shift_share(years, loss, exposure):
    """Laspeyres/Paasche-style additive shift-share on end-point difference.

    ΔL = L_end − L_base, split into:
      exposure effect  = (H_end − H_base) · R_base  (Laspeyres)
      hazard effect    = H_base · (R_end − R_base)
      interaction      = (H_end − H_base) · (R_end − R_base)
    """
    L_base = loss[0]
    L_end = loss[-1]
    H_base = exposure[0]
    H_end = exposure[-1]
    R_base = L_base / H_base
    R_end = L_end / H_end
    delta_L = L_end - L_base
    exposure_effect = (H_end - H_base) * R_base
    hazard_effect = H_base * (R_end - R_base)
    interaction = (H_end - H_base) * (R_end - R_base)
    def share(x):
        return x / delta_L if abs(delta_L) > 1e-9 else float("nan")
    return {
        "L_base": L_base, "L_end": L_end, "delta_L": delta_L,
        "H_base": H_base, "H_end": H_end,
        "R_base": R_base, "R_end": R_end,
        "exposure_effect_Laspeyres": exposure_effect,
        "hazard_effect_Laspeyres": hazard_effect,
        "interaction": interaction,
        "exposure_share": share(exposure_effect),
        "hazard_share": share(hazard_effect),
        "interaction_share": share(interaction),
    }


def run_analysis(data):
    rng = random.Random(RANDOM_SEED)

    years_main, loss_main, exposure_main = build_series(
        data, ANALYSIS_START_YEAR, ANALYSIS_END_YEAR,
        loss_column="structures",
        exposure_variant="total",
        extrapolation_mode=DEFAULT_EXTRAPOLATION_MODE,
    )
    baseline_exposure_main = interp_exposure(
        BASELINE_YEAR,
        data["exposure"]["years"],
        data["exposure"]["housing_units"],
        DEFAULT_EXTRAPOLATION_MODE,
    )

    main = run_decomposition_on(
        years_main, loss_main, exposure_main, baseline_exposure_main, rng,
        do_bootstrap=True, do_permutation=True,
    )

    # Additive (Laspeyres) decomposition — on the aggregated end-points
    additive = additive_shift_share(years_main, loss_main, exposure_main)

    # Descriptive summary
    mean_first_5 = sum(loss_main[:5]) / 5
    mean_last_5 = sum(loss_main[-5:]) / 5
    descriptive = {
        "mean_structures_first_5_yrs": mean_first_5,
        "mean_structures_last_5_yrs": mean_last_5,
        "ratio_last_5_over_first_5": mean_last_5 / mean_first_5,
        "H_baseline_1990": interp_exposure(
            1990, data["exposure"]["years"], data["exposure"]["housing_units"],
            DEFAULT_EXTRAPOLATION_MODE,
        ),
        "H_2010_census": data["exposure"]["housing_units"][-1],
        "H_end_year_default": exposure_main[-1],
        "wui_growth_1990_to_end_multiplier": exposure_main[-1] / interp_exposure(
            1990, data["exposure"]["years"], data["exposure"]["housing_units"],
            DEFAULT_EXTRAPOLATION_MODE,
        ),
    }

    # Sensitivity 1: Alternative analysis start years
    sens_windows = []
    for sy in ALTERNATIVE_START_YEARS:
        ys, ls, es = build_series(
            data, sy, ANALYSIS_END_YEAR, loss_column="structures",
            exposure_variant="total",
            extrapolation_mode=DEFAULT_EXTRAPOLATION_MODE,
        )
        d = run_decomposition_on(ys, ls, es, baseline_exposure_main, rng,
                                 do_bootstrap=False, do_permutation=False)
        sens_windows.append({
            "start_year": sy,
            "end_year": ANALYSIS_END_YEAR,
            "n": d["n"],
            "exposure_share": d["shift_share"]["exposure_share_of_trend"],
            "hazard_share": d["shift_share"]["hazard_share_of_trend"],
            "avoided_loss_share": d["counterfactual"]["avoided_loss_share"],
        })

    # Sensitivity 2: Mega-fire-year exclusion
    keep = [(y, l, e) for y, l, e in zip(years_main, loss_main, exposure_main)
            if y not in MEGA_EVENT_YEARS]
    ys_mf = [t[0] for t in keep]; ls_mf = [t[1] for t in keep]; es_mf = [t[2] for t in keep]
    sens_no_mega = run_decomposition_on(ys_mf, ls_mf, es_mf, baseline_exposure_main, rng,
                                         do_bootstrap=False, do_permutation=False)

    # Sensitivity 3: Alternative extrapolation modes (total-WUI series, main window)
    sens_extrap = []
    for mode in SENSITIVITY_EXTRAPOLATION_MODES:
        ys, ls, es = build_series(
            data, ANALYSIS_START_YEAR, ANALYSIS_END_YEAR,
            loss_column="structures",
            exposure_variant="total",
            extrapolation_mode=mode,
        )
        d = run_decomposition_on(ys, ls, es, baseline_exposure_main, rng,
                                 do_bootstrap=False, do_permutation=False)
        sens_extrap.append({
            "mode": mode,
            "exposure_share": d["shift_share"]["exposure_share_of_trend"],
            "hazard_share": d["shift_share"]["hazard_share_of_trend"],
            "avoided_loss_share": d["counterfactual"]["avoided_loss_share"],
            "H_end_year": es[-1],
        })

    # Sensitivity 4: Interface-only WUI (narrower exposure definition)
    ys, ls, es = build_series(
        data, ANALYSIS_START_YEAR, ANALYSIS_END_YEAR,
        loss_column="structures",
        exposure_variant="interface_only",
        extrapolation_mode=DEFAULT_EXTRAPOLATION_MODE,
    )
    baseline_exposure_interface = interp_exposure(
        BASELINE_YEAR,
        data["exposure"]["years"],
        data["exposure"]["interface_only"],
        DEFAULT_EXTRAPOLATION_MODE,
    )
    sens_interface = run_decomposition_on(
        ys, ls, es, baseline_exposure_interface, rng,
        do_bootstrap=False, do_permutation=False,
    )

    # Sensitivity 5: Residences-only loss column
    ys, ls, es = build_series(
        data, ANALYSIS_START_YEAR, ANALYSIS_END_YEAR,
        loss_column="residences_only",
        exposure_variant="total",
        extrapolation_mode=DEFAULT_EXTRAPOLATION_MODE,
    )
    sens_residences = run_decomposition_on(
        ys, ls, es, baseline_exposure_main, rng,
        do_bootstrap=False, do_permutation=False,
    )

    # Sensitivity 6: Leave-one-year-out on main decomposition
    loo = []
    for i in range(len(years_main)):
        y_loo = years_main[:i] + years_main[i+1:]
        l_loo = loss_main[:i] + loss_main[i+1:]
        e_loo = exposure_main[:i] + exposure_main[i+1:]
        try:
            ss = shift_share_log_decomposition(y_loo, l_loo, e_loo)
            loo.append({
                "dropped_year": years_main[i],
                "exposure_share": ss["exposure_share_of_trend"],
                "hazard_share": ss["hazard_share_of_trend"],
            })
        except (ValueError, ZeroDivisionError):
            continue
    loo_shares = [d["exposure_share"] for d in loo]
    loo_summary = {
        "min_exposure_share": min(loo_shares),
        "max_exposure_share": max(loo_shares),
        "range": max(loo_shares) - min(loo_shares),
        "n": len(loo),
    }

    limitations = [
        "WUI housing-unit values for 2011-2023 are linearly extrapolated from SILVIS 2000-2010 census snapshots; they are not directly observed. The R² near 1.0 on the log-linear exposure trend is an interpolation artifact and must NOT be cited as evidence of a tight empirical fit.",
        "The hazard component β_R captures any non-exposure-linked change (climate, fuels, ignitions, suppression capacity, defensible-space uptake, building code stringency, reporting completeness). The decomposition does NOT attribute β_R to a specific physical cause.",
        "NIFC 'structures destroyed' is an operational tally with known year-to-year noise and reporting-standard drift. Residences-only sensitivity isolates the more consistently-reported subset as a robustness check but does not eliminate the issue.",
        "The shift-share identity β_L = β_H + β_R is algebraic, not causal. Statistical significance of the pairing is tested separately by the marginal-shuffle and circular-shift nulls.",
        "Only 3 SILVIS census snapshots (1990, 2000, 2010) anchor the exposure series; intervening-year values are linearly interpolated.",
        "This analysis does NOT compute dollar losses, interface-vs-intermix disaggregated hazard rates, or county-level structure age; it is a national, count-based, year-resolution decomposition.",
        "The circular-shift null has a discreteness floor of 1/(n-1) ≈ 0.042 at n = 25; p-values below that are not attainable by construction.",
    ]

    results = {
        "meta": {
            "domain": DOMAIN_LABEL,
            "loss_label": LOSS_LABEL,
            "exposure_label": EXPOSURE_LABEL,
            "n_years": len(years_main),
            "year_min": years_main[0],
            "year_max": years_main[-1],
            "baseline_year": BASELINE_YEAR,
            "random_seed": RANDOM_SEED,
            "bootstrap_iterations": BOOTSTRAP_ITERATIONS,
            "permutation_iterations": PERMUTATION_ITERATIONS,
            "block_bootstrap_length": BLOCK_BOOTSTRAP_LENGTH,
            "default_extrapolation_mode": DEFAULT_EXTRAPOLATION_MODE,
            "loss_data_url": LOSS_DATA_URL,
            "exposure_data_url": EXPOSURE_DATA_URL,
            "loss_sha256": data["loss_sha256"],
            "exposure_sha256": data["exposure_sha256"],
            "nifc_url_reachable": data["nifc_url_reachable"],
            "silvis_url_reachable": data["silvis_url_reachable"],
            "python_version": sys.version.split()[0],
            "platform": sys.platform,
            "r2_log_exposure_trend_is_interpolation_artifact": True,
            "r2_log_exposure_trend_note": (
                "The near-unity R^2 on the log-exposure trend is a known "
                "artifact of piecewise-linear interpolation between the 3 "
                "SILVIS census snapshots (1990/2000/2010) plus a linear "
                "2000-2010 post-2010 extrapolation. It is NOT evidence of "
                "an empirical near-perfect fit and MUST NOT be cited as one."
            ),
        },
        "descriptive": descriptive,
        "main": main,
        "additive_shift_share": additive,
        "sensitivity": {
            "analysis_window": sens_windows,
            "exclude_mega_fire_years": {
                "excluded_years": MEGA_EVENT_YEARS,
                "exposure_share": sens_no_mega["shift_share"]["exposure_share_of_trend"],
                "hazard_share": sens_no_mega["shift_share"]["hazard_share_of_trend"],
                "avoided_loss_share": sens_no_mega["counterfactual"]["avoided_loss_share"],
            },
            "wui_extrapolation_mode": sens_extrap,
            "interface_only_wui": {
                "exposure_share": sens_interface["shift_share"]["exposure_share_of_trend"],
                "hazard_share": sens_interface["shift_share"]["hazard_share_of_trend"],
                "avoided_loss_share": sens_interface["counterfactual"]["avoided_loss_share"],
            },
            "residences_only_losses": {
                "exposure_share": sens_residences["shift_share"]["exposure_share_of_trend"],
                "hazard_share": sens_residences["shift_share"]["hazard_share_of_trend"],
                "avoided_loss_share": sens_residences["counterfactual"]["avoided_loss_share"],
            },
            "leave_one_out_on_exposure_share": loo_summary,
        },
        "limitations": limitations,
    }
    return results


# ═══════════════════════════════════════════════════════════════
# REPORT GENERATION
# ═══════════════════════════════════════════════════════════════


def fmt_num(x, sig=4):
    if x is None:
        return "NA"
    if isinstance(x, float) and math.isnan(x):
        return "NA"
    if isinstance(x, int):
        return f"{x:,}"
    if isinstance(x, float):
        if abs(x) >= 1000:
            return f"{x:,.0f}"
        return f"{x:.{sig}g}"
    return str(x)


def fmt_pct(x):
    if x is None or (isinstance(x, float) and math.isnan(x)):
        return "NA"
    return f"{100.0 * x:.1f}%"


def generate_report(results):
    with open(RESULTS_JSON, "w") as f:
        json.dump(results, f, indent=2, default=str)

    m = results["meta"]
    d = results["descriptive"]
    main = results["main"]
    ss = main["shift_share"]
    cf = main["counterfactual"]
    bb = main["moving_block_bootstrap"]
    add = results["additive_shift_share"]
    sens = results["sensitivity"]

    lines = []
    lines.append(f"# {m['domain']}: exposure vs. hazard decomposition")
    lines.append("")
    lines.append(f"- Analysis window: {m['year_min']}–{m['year_max']} (N = {m['n_years']})")
    lines.append(f"- Baseline (frozen-WUI) year: {m['baseline_year']}")
    lines.append(f"- Random seed: {m['random_seed']}; bootstrap iters: {m['bootstrap_iterations']}; "
                 f"permutation iters: {m['permutation_iterations']}; block length: {m['block_bootstrap_length']}")
    lines.append(f"- NIFC loss series SHA-256: `{m['loss_sha256']}`")
    lines.append(f"- SILVIS exposure census SHA-256: `{m['exposure_sha256']}`")
    lines.append("")

    lines.append("## Descriptive")
    lines.append(f"- Mean {LOSS_LABEL} first 5 yrs ({m['year_min']}-{m['year_min']+4}): "
                 f"{fmt_num(d['mean_structures_first_5_yrs'])}")
    lines.append(f"- Mean {LOSS_LABEL} last 5 yrs ({m['year_max']-4}-{m['year_max']}): "
                 f"{fmt_num(d['mean_structures_last_5_yrs'])}")
    lines.append(f"- Ratio last-5 / first-5: {fmt_num(d['ratio_last_5_over_first_5'])}x")
    lines.append(f"- WUI housing at {m['baseline_year']}: {fmt_num(d['H_baseline_1990'])}")
    lines.append(f"- WUI housing at {m['year_max']} (extrapolated): {fmt_num(d['H_end_year_default'])}")
    lines.append(f"- WUI housing growth {m['baseline_year']}→{m['year_max']}: "
                 f"{fmt_num(d['wui_growth_1990_to_end_multiplier'])}x")
    lines.append("")

    lines.append("## Log-linear shift-share decomposition (main)")
    lines.append(f"- β_log_loss (trend in ln(L)):     {fmt_num(ss['beta_log_loss_per_year'])} / yr, "
                 f"R² = {fmt_num(ss['r2_log_loss_trend'])}")
    lines.append(f"- β_log_exposure (trend in ln(H)): {fmt_num(ss['beta_log_exposure_per_year'])} / yr, "
                 f"R² = {fmt_num(ss['r2_log_exposure_trend'])}")
    lines.append(f"- β_log_hazard (trend in ln(R)):   {fmt_num(ss['beta_log_hazard_per_year'])} / yr, "
                 f"R² = {fmt_num(ss['r2_log_hazard_trend'])}")
    lines.append(f"- Identity residual |β_L − (β_H + β_R)|: {fmt_num(ss['identity_residual_abs'])}")
    lines.append(f"- **Exposure share of trend**: {fmt_pct(ss['exposure_share_of_trend'])}")
    lines.append(f"- **Hazard share of trend**:  {fmt_pct(ss['hazard_share_of_trend'])}")
    if bb:
        lines.append(f"- Exposure share 95% CI (block bootstrap): "
                     f"[{fmt_pct(bb['exposure_share_ci'][0])}, {fmt_pct(bb['exposure_share_ci'][1])}]")
        lines.append(f"- Hazard share 95% CI (block bootstrap): "
                     f"[{fmt_pct(bb['hazard_share_ci'][0])}, {fmt_pct(bb['hazard_share_ci'][1])}]")
    lines.append("")

    lines.append("## Constant-1990-WUI counterfactual")
    lines.append(f"- Baseline H (1990): {fmt_num(cf['baseline_exposure'])}")
    lines.append(f"- Σ observed losses: {fmt_num(cf['sum_observed'])}")
    lines.append(f"- Σ counterfactual losses: {fmt_num(cf['sum_counterfactual'])}")
    lines.append(f"- **Avoided-loss share (1 − CF/obs)**: {fmt_pct(cf['avoided_loss_share'])}")
    if bb:
        lines.append(f"- Avoided-loss share 95% CI: "
                     f"[{fmt_pct(bb['avoided_loss_share_ci'][0])}, {fmt_pct(bb['avoided_loss_share_ci'][1])}]")
    lines.append("")

    lines.append("## Null tests")
    perm_marg = main["marginal_shuffle_null"]
    if perm_marg:
        lines.append(f"### Marginal-shuffle null (destroys exposure time-order)")
        lines.append(f"- Iterations: {perm_marg['iterations']}")
        lines.append(f"- Null mean: {fmt_pct(perm_marg['null_exposure_share_mean'])}; "
                     f"null 2.5–97.5% [{fmt_pct(perm_marg['null_exposure_share_p025'])}, "
                     f"{fmt_pct(perm_marg['null_exposure_share_p975'])}]")
        lines.append(f"- **Two-sided p**: {fmt_num(perm_marg['p_two_sided'])}")
    perm_circ = main["circular_shift_null"]
    if perm_circ:
        lines.append(f"### Circular-shift null (preserves exposure autocorrelation)")
        lines.append(f"- Shifts evaluated: {perm_circ['shifts_evaluated']}; "
                     f"minimum attainable p = {fmt_num(perm_circ['minimum_attainable_p'])}")
        lines.append(f"- Null mean: {fmt_pct(perm_circ['null_exposure_share_mean'])}; "
                     f"null 2.5–97.5% [{fmt_pct(perm_circ['null_exposure_share_p025'])}, "
                     f"{fmt_pct(perm_circ['null_exposure_share_p975'])}]")
        lines.append(f"- **Two-sided p**: {fmt_num(perm_circ['p_two_sided'])}")
    lines.append("")

    lines.append("## Additive (Laspeyres) shift-share on end-points")
    lines.append(f"- ΔL = L_{m['year_max']} − L_{m['year_min']} = "
                 f"{fmt_num(add['delta_L'])} structures")
    lines.append(f"- Exposure effect (Laspeyres): {fmt_num(add['exposure_effect_Laspeyres'])} "
                 f"→ share {fmt_pct(add['exposure_share'])}")
    lines.append(f"- Hazard effect (Laspeyres):   {fmt_num(add['hazard_effect_Laspeyres'])} "
                 f"→ share {fmt_pct(add['hazard_share'])}")
    lines.append(f"- Interaction: {fmt_num(add['interaction'])} → share {fmt_pct(add['interaction_share'])}")
    lines.append("")

    lines.append("## Sensitivity: analysis window")
    for s in sens["analysis_window"]:
        lines.append(f"- {s['start_year']}-{s['end_year']} (N={s['n']}): "
                     f"exposure share = {fmt_pct(s['exposure_share'])}, "
                     f"hazard share = {fmt_pct(s['hazard_share'])}, "
                     f"avoided-loss share = {fmt_pct(s['avoided_loss_share'])}")
    lines.append("")

    lines.append("## Sensitivity: WUI post-2010 extrapolation mode")
    for s in sens["wui_extrapolation_mode"]:
        lines.append(f"- `{s['mode']}`: H_end = {fmt_num(s['H_end_year'])}; "
                     f"exposure share = {fmt_pct(s['exposure_share'])}; "
                     f"hazard share = {fmt_pct(s['hazard_share'])}; "
                     f"avoided-loss share = {fmt_pct(s['avoided_loss_share'])}")
    lines.append("")

    lines.append("## Sensitivity: interface-only WUI (narrower exposure)")
    s = sens["interface_only_wui"]
    lines.append(f"- Exposure share: {fmt_pct(s['exposure_share'])}; "
                 f"hazard share: {fmt_pct(s['hazard_share'])}; "
                 f"avoided-loss share: {fmt_pct(s['avoided_loss_share'])}")
    lines.append("")

    lines.append("## Sensitivity: residences-only losses")
    s = sens["residences_only_losses"]
    lines.append(f"- Exposure share: {fmt_pct(s['exposure_share'])}; "
                 f"hazard share: {fmt_pct(s['hazard_share'])}; "
                 f"avoided-loss share: {fmt_pct(s['avoided_loss_share'])}")
    lines.append("")

    lines.append("## Sensitivity: mega-fire-year exclusion")
    s = sens["exclude_mega_fire_years"]
    lines.append(f"- Excluded: {s['excluded_years']}")
    lines.append(f"- Exposure share: {fmt_pct(s['exposure_share'])}; "
                 f"hazard share: {fmt_pct(s['hazard_share'])}; "
                 f"avoided-loss share: {fmt_pct(s['avoided_loss_share'])}")
    lines.append("")

    lines.append("## Sensitivity: leave-one-year-out")
    l = sens["leave_one_out_on_exposure_share"]
    lines.append(f"- Exposure share range across {l['n']} LOO iterations: "
                 f"[{fmt_pct(l['min_exposure_share'])}, {fmt_pct(l['max_exposure_share'])}] "
                 f"(spread = {fmt_pct(l['range'])})")
    lines.append("")

    lines.append("## Limitations and out-of-scope")
    for lim in results["limitations"]:
        lines.append(f"- {lim}")
    lines.append("")

    with open(REPORT_MD, "w") as f:
        f.write("\n".join(lines))


# ═══════════════════════════════════════════════════════════════
# VERIFICATION
# ═══════════════════════════════════════════════════════════════


def verify():
    if not os.path.exists(RESULTS_JSON):
        print("[FAIL] results.json not found; run the analysis first.")
        return 1
    with open(RESULTS_JSON) as f:
        R = json.load(f)

    failures = []

    def check(cond, msg):
        if not cond:
            failures.append(msg)

    m = R["meta"]
    main = R["main"]
    ss = main["shift_share"]
    cf = main["counterfactual"]
    perm_marg = main["marginal_shuffle_null"]
    perm_circ = main["circular_shift_null"]
    bb = main["moving_block_bootstrap"]
    sens = R["sensitivity"]

    # 1. Coverage
    check(m["n_years"] >= 20, f"[1] Expected ≥20 years; got {m['n_years']}")
    check(m["year_min"] <= 2000 and m["year_max"] >= 2020,
          f"[2] Expected span covering 1999-2023; got {m['year_min']}-{m['year_max']}")

    # 2. Shift-share identity holds numerically
    check(ss["identity_residual_abs"] < 1e-8,
          f"[3] Shift-share identity failed; |β_L − (β_H + β_R)| = {ss['identity_residual_abs']}")

    # 3. Exposure + hazard shares sum to 1
    tot = ss["exposure_share_of_trend"] + ss["hazard_share_of_trend"]
    check(abs(tot - 1.0) < 1e-8, f"[4] Shares do not sum to 1; got {tot}")

    # 4. WUI exposure grew over time (β_H > 0)
    check(ss["beta_log_exposure_per_year"] > 0,
          f"[5] Expected positive exposure trend; got {ss['beta_log_exposure_per_year']}")

    # 5. Loss trend grew over time (β_L > 0)
    check(ss["beta_log_loss_per_year"] > 0,
          f"[6] Expected positive loss trend; got {ss['beta_log_loss_per_year']}")

    # 6. Hazard-rate trend > 0 (fire behavior has intensified)
    check(ss["beta_log_hazard_per_year"] > 0,
          f"[7] Expected positive hazard-rate trend; got {ss['beta_log_hazard_per_year']}")

    # 7. Bootstrap iteration count is reasonable
    check(bb is not None and bb["iterations"] >= 4000,
          f"[8] Expected ≥4000 bootstrap samples; got {bb['iterations'] if bb else None}")

    # 8. Marginal-shuffle iteration count is reasonable and p is in [0,1]
    check(perm_marg is not None and perm_marg["iterations"] >= 1800,
          f"[9] Expected ≥1800 marginal-shuffle samples; got {perm_marg['iterations'] if perm_marg else None}")
    check(perm_marg is not None and 0.0 <= perm_marg["p_two_sided"] <= 1.0,
          f"[10] Marginal-shuffle p-value out of range; got {perm_marg['p_two_sided'] if perm_marg else None}")
    # 8b. Circular-shift null was evaluated and p is in [0,1]
    check(perm_circ is not None and perm_circ["shifts_evaluated"] >= 20,
          f"[10b] Expected ≥20 circular shifts; got {perm_circ['shifts_evaluated'] if perm_circ else None}")
    check(perm_circ is not None and 0.0 <= perm_circ["p_two_sided"] <= 1.0,
          f"[10c] Circular-shift p-value out of range; got {perm_circ['p_two_sided'] if perm_circ else None}")

    # 9. Exposure share finite and in [-2, 2]
    check(-2.0 <= ss["exposure_share_of_trend"] <= 2.0,
          f"[11] Exposure share implausible; got {ss['exposure_share_of_trend']}")

    # 10. Avoided-loss share in (-1, 1)
    check(-1.0 < cf["avoided_loss_share"] < 1.0,
          f"[12] Avoided-loss share out of plausible range; got {cf['avoided_loss_share']}")

    # 11. Sums of counterfactual and observed losses are positive
    check(cf["sum_observed"] > 0 and cf["sum_counterfactual"] > 0,
          f"[13] Non-positive loss totals; got obs={cf['sum_observed']}, cf={cf['sum_counterfactual']}")

    # 12. WUI growth multiplier is > 1.2 (real-world SILVIS shows ~1.5)
    mult = R["descriptive"]["wui_growth_1990_to_end_multiplier"]
    check(mult > 1.2, f"[14] Expected WUI growth >1.2x over 1990-2023; got {mult}")

    # 13. At least 3 extrapolation-mode sensitivities ran
    check(len(sens["wui_extrapolation_mode"]) == 3,
          f"[15] Expected 3 extrapolation sensitivities; got {len(sens['wui_extrapolation_mode'])}")

    # 14. Flat-2010 extrapolation should produce lower exposure share than linear_2000_2010
    shares_by_mode = {s["mode"]: s["exposure_share"] for s in sens["wui_extrapolation_mode"]}
    check(shares_by_mode["flat_2010"] <= shares_by_mode["linear_2000_2010"],
          f"[16] Expected flat_2010 exposure share ≤ linear_2000_2010; got {shares_by_mode}")

    # 15. Leave-one-out exposure-share spread is finite and reasonable
    loo = sens["leave_one_out_on_exposure_share"]
    check(loo["n"] >= 20 and 0.0 <= loo["range"] < 2.0,
          f"[17] LOO summary implausible; got n={loo['n']}, range={loo['range']}")

    # 16. Residual lag-1 ACF is in [-1, 1]
    acf = main["residual_acf_lag1_log_loss"]
    check(-1.0 <= acf <= 1.0, f"[18] Lag-1 ACF out of [-1,1]; got {acf}")

    # 17. Mega-fire-exclusion sensitivity ran and share is plausible
    mf = sens["exclude_mega_fire_years"]
    check(-2.0 <= mf["exposure_share"] <= 2.0,
          f"[19] Mega-fire exclusion exposure share implausible; got {mf['exposure_share']}")

    # 18. Sensitivity: analysis-window sweep covers at least 2 alt windows
    check(len(sens["analysis_window"]) >= 2,
          f"[20] Expected ≥2 alt windows; got {len(sens['analysis_window'])}")

    # 19. Limitations block is present and non-trivial (≥4 caveats)
    lims = R.get("limitations", [])
    check(isinstance(lims, list) and len(lims) >= 4,
          f"[21] Expected ≥4 documented limitations; got {len(lims)}")

    # 20. Random seed recorded and is integer 42 (reproducibility)
    check(m.get("random_seed") == 42,
          f"[22] Expected random_seed=42; got {m.get('random_seed')}")

    # 21. Both SHA-256 hashes are 64-char hex strings (not the BOOTSTRAP sentinel)
    check(isinstance(m.get("loss_sha256"), str) and len(m["loss_sha256"]) == 64,
          f"[23] loss_sha256 missing or wrong length")
    check(isinstance(m.get("exposure_sha256"), str) and len(m["exposure_sha256"]) == 64,
          f"[24] exposure_sha256 missing or wrong length")

    # 22. Identity holds under leave-one-out: share range is strictly less than 1.5
    check(loo["range"] < 1.5,
          f"[25] LOO exposure-share spread too wide ({loo['range']}); main result not robust")

    # 23. Additive (Laspeyres) decomposition ran and its exposure+hazard+interaction shares sum to ≈1
    add = R["additive_shift_share"]
    add_sum = add["exposure_share"] + add["hazard_share"] + add["interaction_share"]
    check(abs(add_sum - 1.0) < 1e-8,
          f"[26] Additive shares do not sum to 1; got {add_sum}")

    # 24. The R² interpolation-artifact flag is present in meta (honesty).
    check(m.get("r2_log_exposure_trend_is_interpolation_artifact") is True,
          f"[27] Missing r2_log_exposure_trend_is_interpolation_artifact flag in meta")

    # 25. Both null p-values are well-defined floats in [0, 1].
    check(isinstance(perm_marg.get("p_two_sided"), float)
          and 0.0 <= perm_marg["p_two_sided"] <= 1.0,
          f"[28] Marginal-shuffle p missing or out of [0,1]")
    check(isinstance(perm_circ.get("p_two_sided"), float)
          and 0.0 <= perm_circ["p_two_sided"] <= 1.0,
          f"[29] Circular-shift p missing or out of [0,1]")

    # 26. Sensitivity reports are all present (window, mega, extrapolation,
    # interface-only, residences-only, LOO).
    check(set(sens.keys()) >= {
        "analysis_window", "exclude_mega_fire_years", "wui_extrapolation_mode",
        "interface_only_wui", "residences_only_losses",
        "leave_one_out_on_exposure_share",
    }, f"[30] Not all 6 sensitivity panels present; got {sorted(sens.keys())}")

    n_assertions = 30
    if failures:
        print("VERIFICATION FAILED")
        for msg in failures:
            print("  -", msg)
        return 1
    else:
        print(f"VERIFICATION PASSED ({n_assertions} assertions)")
        return 0


# ═══════════════════════════════════════════════════════════════
# ENTRY POINT
# ═══════════════════════════════════════════════════════════════


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--verify", action="store_true",
                        help="Run assertions against saved results.json and exit.")
    args = parser.parse_args()

    if args.verify:
        rc = verify()
        sys.exit(rc)

    print("[1/5] Loading NIFC structure-loss series and SILVIS WUI census")
    data = load_data()
    print(f"  loss CSV sha256 = {data['loss_sha256']}")
    print(f"  exposure sha256 = {data['exposure_sha256']}")
    print(f"  NIFC URL reachable: {data['nifc_url_reachable']}; "
          f"SILVIS URL reachable: {data['silvis_url_reachable']}")

    print(f"[2/5] Assembling {ANALYSIS_START_YEAR}-{ANALYSIS_END_YEAR} paired series "
          f"(extrapolation = {DEFAULT_EXTRAPOLATION_MODE})")
    results = run_analysis(data)
    ss = results["main"]["shift_share"]
    print(f"  β_log_loss = {ss['beta_log_loss_per_year']:.4f}/yr")
    print(f"  β_log_exposure = {ss['beta_log_exposure_per_year']:.4f}/yr")
    print(f"  β_log_hazard = {ss['beta_log_hazard_per_year']:.4f}/yr")

    print(f"[3/5] Running {BOOTSTRAP_ITERATIONS}-iteration moving-block bootstrap (block={BLOCK_BOOTSTRAP_LENGTH})")
    bb = results["main"]["moving_block_bootstrap"]
    print(f"  exposure share = {ss['exposure_share_of_trend']*100:.1f}%; "
          f"95% CI = [{bb['exposure_share_ci'][0]*100:.1f}%, {bb['exposure_share_ci'][1]*100:.1f}%]")
    print(f"  hazard share  = {ss['hazard_share_of_trend']*100:.1f}%")
    cf = results["main"]["counterfactual"]
    print(f"  avoided-loss share = {cf['avoided_loss_share']*100:.1f}%; "
          f"95% CI = [{bb['avoided_loss_share_ci'][0]*100:.1f}%, {bb['avoided_loss_share_ci'][1]*100:.1f}%]")

    print(f"[4/5] Running null tests ({PERMUTATION_ITERATIONS}-iter marginal-shuffle + exact circular-shift)")
    perm_marg = results["main"]["marginal_shuffle_null"]
    perm_circ = results["main"]["circular_shift_null"]
    print(f"  marginal-shuffle: null mean = {perm_marg['null_exposure_share_mean']*100:.1f}%; "
          f"p = {perm_marg['p_two_sided']:.4f}")
    print(f"  circular-shift:   null mean = {perm_circ['null_exposure_share_mean']*100:.1f}%; "
          f"p = {perm_circ['p_two_sided']:.4f} "
          f"(min attainable p = {perm_circ['minimum_attainable_p']:.4f})")

    print("[5/5] Writing results and report")
    generate_report(results)
    print(f"  wrote {RESULTS_JSON}")
    print(f"  wrote {REPORT_MD}")
    print("ANALYSIS COMPLETE")


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

**Expected output:** No output (heredoc writes file silently).

**Success condition:** `/tmp/claw4s_auto_wui-expansion-explains-more-structure-loss-variance-than-fir/analyze.py` exists and is a valid Python file.

---

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_wui-expansion-explains-more-structure-loss-variance-than-fir && python3 analyze.py
```

**Expected output (final line must be `ANALYSIS COMPLETE`):**

```
[1/5] Loading NIFC structure-loss series and SILVIS WUI census
  loss CSV sha256 = <64-hex-char SHA-256>
  exposure sha256 = <64-hex-char SHA-256>
  NIFC URL reachable: True | False; SILVIS URL reachable: True | False
[2/5] Assembling 1999-2023 paired series (extrapolation = linear_2000_2010)
  β_log_loss = ... /yr
  β_log_exposure = ... /yr
  β_log_hazard = ... /yr
[3/5] Running 5000-iteration moving-block bootstrap (block=3)
  exposure share = ...%; 95% CI = [..., ...]
  hazard share  = ...%
  avoided-loss share = ...%; 95% CI = [..., ...]
[4/5] Running null tests (2000-iter marginal-shuffle + exact circular-shift)
  marginal-shuffle: null mean = ...%; p = ...
  circular-shift:   null mean = ...%; p = ... (min attainable p = ...)
[5/5] Writing results and report
  wrote results.json
  wrote report.md
ANALYSIS COMPLETE
```

**Success condition:** exit code 0; `results.json` and `report.md` produced; final stdout line is `ANALYSIS COMPLETE`.

**Failure conditions:**
- `ValueError: Singular matrix` → embedded data was corrupted; restore from source.
- `Data integrity check failed` → the embedded CSV has been edited since its SHA-256 was pinned.
- Network errors → non-fatal; analysis still runs on embedded data.

---

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_wui-expansion-explains-more-structure-loss-variance-than-fir && python3 analyze.py --verify
```

**Expected output:**
```
VERIFICATION PASSED (30 assertions)
```

**Success condition:** exit code 0 and the string `VERIFICATION PASSED` on stdout.

**Failure conditions:** exit code 1; stdout begins with `VERIFICATION FAILED` followed by the failed assertion identifiers.

---

## Files produced

- `results.json` — full structured results (shift-share, counterfactual, null test, CIs, all sensitivities).
- `report.md` — human-readable summary with numbers for the research note.
- `cache/nifc_structures.csv` — canonical annual loss CSV, pinned by SHA-256.
- `cache/silvis_wui_census.json` — canonical SILVIS census snapshot, pinned by SHA-256.
- `cache/.data_sha256.txt` — observed hashes recorded at most recent run.