---
name: "WUI Structure Loss: Fire-Behaviour Change vs. Housing Expansion"
description: "Decomposes the 2000-2023 rise in California wildfire structure destruction into exposure-driven (housing-expansion) and behaviour-driven (loss-per-structure) components using a housing-fixed counterfactual. Applies block bootstrap confidence intervals over years, a year-label permutation null for the housing-expansion contribution, and sensitivity analyses over start year, mega-fire year exclusion, and alternative exposure denominators."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "wildfire", "wui", "california", "hazards", "exposure", "decomposition", "bootstrap", "permutation-test", "reproducibility"]
python_version: ">=3.8"
dependencies: []
data_source: "CAL FIRE Redbook annual fire statistics (https://www.fire.ca.gov/what-we-do/fire-resource-assessment-program/fire-perimeters); U.S. Census Bureau Housing Unit Estimates (https://www.census.gov/data/tables/time-series/demo/popest/2010s-state-detail.html, https://www.census.gov/data/tables/time-series/demo/popest/2020s-state-detail.html); California Department of Finance E-5 series (https://dof.ca.gov/forecasting/demographics/estimates/e-5-population-and-housing-estimates-for-cities-counties-and-the-state-2020-2024/)."
data_revision: "Structure-loss counts are CAL FIRE Redbook annual totals for all jurisdictions (direct protection area + state responsibility area). Housing unit counts are US Census Bureau Population Estimates Program (PEP) state-level housing unit estimates for California, vintage 2023 release; values are SHA256-locked in the analysis script."
---

# WUI Structure Loss: Fire-Behaviour Change vs. Housing Expansion

## Research Question

**Is the 2000-2023 rise in California wildfire structure destruction primarily a
"fire-behaviour" signal (each structure now faces a higher per-unit hazard) or a
"housing-expansion" signal (more structures are simply exposed)?**

We decompose the log-change in annual structure loss (L) between an early window
(2000-2009) and a late window (2014-2023) into

    Δ log(L) = Δ log(H) + Δ log(L / H)

where H is the annual count of California housing units (the exposure denominator)
and L/H is the per-housing-unit loss rate.  The **housing-expansion share** is
`Δ log(H) / Δ log(L)`; the **loss-rate share** is `Δ log(L/H) / Δ log(L)`.

A housing-expansion share near 1 means "the rise is mostly because there are more
houses"; a loss-rate share near 1 means "per-house hazard has risen".  We further
compute a housing-fixed counterfactual loss series
`L_t × (H_early / H_t)` that holds exposure at the early-window geometric mean,
and use a circular block bootstrap + year-label permutation null to place
frequentist uncertainty on each share.

## When to Use This Skill

Use this skill when you need to test whether the widely-reported rise in wildfire structure destruction in California (and in WUI regions generally) is driven primarily by changing fire behaviour (hotter, more intense, faster-spreading fires) or by expansion of the housing stock that is exposed to fire — i.e., when you need to separate the "fire side" from the "exposure side" of a rising hazard curve before attributing it to climate or to land-use policy.

### Preconditions
- **Python version:** 3.8+ (standard library only — no numpy, scipy, pandas, requests, or third-party statistical packages).
- **Network:** None required. Primary input data are pinned tables (CAL FIRE Redbook annual totals; US Census PEP housing unit estimates for California) embedded as Python literals and SHA256-verified at load time. The script will optionally fetch a Census-provided CSV as an external provenance anchor if `CLAW_FETCH_CENSUS=1` is set in the environment; this is disabled by default to keep the skill offline-deterministic.
- **Runtime:** Fresh run ~30-90 seconds (10,000 block-bootstrap resamples and 10,000 year-label permutations on a 24-year series). `--verify` completes in <10 seconds.
- **Disk:** <5 MB (the analysis writes only `results.json`, `report.md`, and a small provenance JSON).
- **CPU:** Single-threaded; all statistical computation is pure Python.

## Adaptation Guidance

The statistical core of this skill (log-ratio Kaya-style exposure decomposition, block bootstrap over years, year-label permutation null for the exposure contribution, and sensitivity sweep over window and outlier-exclusion choices) is domain-agnostic and can be re-pointed at any hazard-loss-vs-exposure decomposition question (flood damages vs. coastal housing stock, hailstorm insurance claims vs. vehicles registered, heat-wave deaths vs. population >65, etc.).

- **Change the jurisdiction, hazard, or time window:** Edit the `DOMAIN CONFIGURATION` block in `analysis.py`. `STRUCTURE_LOSS_BY_YEAR` holds the annual loss count; `HOUSING_UNITS_BY_YEAR` holds the annual exposure count. Replace both tables with any paired year→count dictionaries covering the same years. `START_YEAR` and `END_YEAR` constrain the analysis window.
- **Change the exposure denominator:** Replace `HOUSING_UNITS_BY_YEAR` with any exposure proxy (WUI-specific housing counts, population, insured value). The decomposition is index-neutral. Alternative denominators can be supplied via `ALT_EXPOSURE_BY_YEAR` and are tested in the sensitivity section automatically.
- **Change the sensitivity set:** `MEGA_FIRE_YEARS` lists years to exclude in a leave-out sensitivity analysis. `SENSITIVITY_WINDOWS` gives alternative `(start, end)` ranges used to check whether the finding is robust to endpoint choice.
- **Change the null model strength:** `N_BOOTSTRAP` and `N_PERMUTATION` control Monte-Carlo resolution; `BLOCK_LENGTH` controls the circular block bootstrap block size (set to 1 for IID bootstrap on annual data).
- **Swap the loss variable:** Any annual count or dollar figure works. Dollar figures inflate interpretation but don't change the decomposition math because logs are taken.
- **What stays the same:** `fisher_log_decomposition()` (exact additive decomposition of Δlog(L) into Δlog(H) + Δlog(L/H)), `block_bootstrap_ci()` (circular block bootstrap over years), `permutation_null_housing_share()` (year-label shuffle), `welch_t()`, `sens_slope()`, and the `--verify` assertion harness are all general-purpose.
- **Replace the data entirely:** Put any paired series in `STRUCTURE_LOSS_BY_YEAR` and `HOUSING_UNITS_BY_YEAR`. Update `STRUCTURE_LOSS_SHA256` and `HOUSING_UNITS_SHA256` to the hashes of the new canonical strings (the script prints the expected hash if the constants are set to the empty string).

## Overview

**The claim under test.** Annual counts of structures destroyed by wildfire in California rose from a few hundred in the early 2000s to five figures in the late 2010s, and press coverage routinely treats this as evidence of an intensifying wildfire regime driven by climate change. At the same time, the California housing stock grew, and an increasing fraction of that stock is situated in the wildland-urban interface (WUI). Before attributing rising losses to fire-behaviour change, one must hold the *exposure denominator* — the number of structures at risk — fixed.

**The confound.** Structure loss is a joint outcome of hazard (fire behaviour) and exposure (structures at risk). If exposure grows, loss can rise mechanically even if per-structure hazard is flat. Unadjusted loss counts therefore conflate two independent drivers. A correct test of "are fires worse?" requires decomposing Δ(loss) into Δ(exposure) and Δ(loss-per-exposure).

**What this skill does.** Using CAL FIRE Redbook annual structure-destroyed totals for California 2000-2023 and US Census Bureau Population Estimates Program (PEP) state-level housing unit counts for California, it computes the log-additive Kaya-style decomposition

    Δ log(L) = Δ log(H) + Δ log(L / H)

year-over-year and over the full window. The *housing-expansion share* is the fraction of the total log-change in L attributable to Δ log(H); the *loss-rate share* is the complement. A housing-fixed counterfactual loss series is computed as `L_t × (H_2000 / H_t)` and compared to the observed series.

**Null models.** (a) *Circular block bootstrap* over the year index (block length 3), 10,000 resamples — 95% percentile CIs on the housing-expansion share, the loss-rate share, and the housing-fixed counterfactual's rise. (b) *Year-label permutation* — shuffle the assignment of housing-unit counts to years 10,000 times and recompute the housing-expansion share under each permutation; report the fraction of permutations in which the permuted share equals or exceeds the observed share in absolute value. This tests whether the observed housing contribution is distinguishable from a purely random pairing of years.

**Methodological hook.** Exposure-adjusted loss. Most popular discussion of "are fires getting worse?" compares raw loss counts. This skill makes the exposure denominator explicit, uses a housing-fixed counterfactual to quantify how much of the observed rise survives holding housing constant, and benchmarks the housing contribution against a year-label permutation null. We additionally test robustness to the inclusion of mega-fire years (2017, 2018, 2020, 2021), which concentrate most of the loss, and to the choice of exposure denominator (total housing vs. California Department of Finance single-family unit counts).

**What this is not.** Not a causal identification of climate's role in California fire behaviour; the residual "loss-rate" component bundles climate, fuel, suppression doctrine, ignition patterns, and the changing spatial overlap of housing with fire perimeters. Not a WUI-raster analysis (SILVIS WUI change layers are not parseable in pure-Python stdlib). Not a structure-by-structure exposure estimate; it uses state-level housing unit totals as the exposure proxy.

## Limitations, Assumptions, and Failure Modes

These caveats are also emitted into `report.md` by the script, so they travel with the artefact.

1. **Data limitation — aggregation.** CAL FIRE Redbook annual totals conflate state-responsibility-area and direct-protection-area losses; they do not stratify by WUI class, structure type, or jurisdiction. The analysis therefore cannot distinguish exurban/WUI losses from backcountry losses, and a single mega-fire can dominate an annual count.
2. **Methodological assumption — exposure proxy.** Exposure is proxied by state-level housing-unit totals, not WUI-specific unit counts. The Kaya-style decomposition is exact with respect to this proxy, but it cannot reveal intra-state spatial shifts of housing into higher-hazard zones (e.g., the Sierra Nevada foothills). The `ALT_EXPOSURE_BY_YEAR` sensitivity partially addresses the choice-of-denominator question but not the spatial-overlap question.
3. **Methodological assumption — geometric means and logs.** The primary decomposition uses geometric means of annual counts within decade-long windows. This is exact for an identity in log space but implicitly down-weights extreme years; a reader who wants an arithmetic-mean-based decomposition must restate the identity (which is no longer additive in logs).
4. **Methodological assumption — independence.** Circular block bootstrap (block length 3) is our concession to within-window serial correlation. Measured lag-1 autocorrelation is reported in the output. If the true dependence structure is long-memory, block-bootstrap CIs will be too narrow; the reported `lag1_log_loss_full_window` and `lag1_log_rate_full_window` diagnose this.
5. **Scenario that breaks the analysis.** Zero-loss years, or years where `STRUCTURE_LOSS_BY_YEAR` is missing a value inside the decade windows, cause the log decomposition to error; the script halts with a clear `ValueError` (exit 3) rather than silently imputing.
6. **Temporal scope.** The 2000–2023 window is short (24 years). A longer series could revise both the loss-rate and the housing-expansion attributions. Sensitivity over alternative start/end years and an outlier-excluded (mega-fire-excluded) window are reported as robustness checks; the finding survives both, but one cannot rule out a regime shift outside the sample.
7. **What the results do NOT show.** They do NOT identify climate's causal role in fire behaviour, do NOT attribute loss to any specific fire or ignition source, and do NOT control for the spatial overlap of housing with fire perimeters. They also do NOT measure dollar value of loss, insurance payouts, or displacement; the loss variable is a structure count.
8. **Deterministic but not bit-for-bit across Python versions.** Results are seeded (`SEED=42`) and the same Python patch version reproduces `results.json` byte-for-byte; cross-version `math.log` differences in the last ulp can in principle perturb the 12th decimal, but all published and `--verify`-asserted quantities are robust to this.

## Controls and Comparators

The script evaluates the observed housing-expansion share against four distinct nulls / comparators, so that the finding is not a single-test artefact:

- **(C1) Circular block bootstrap** (`N_BOOTSTRAP=10000`, block length 3): 95% percentile CIs on `share_H`, `share_R`, `delta_logL`, and `delta_logH`, resampled within early and late windows.
- **(C2) Year-label permutation null** (`N_PERMUTATION=10000`): housing values are shuffled across the union of early + late years; the resulting empirical distribution of `share_H` under the null of exchangeable year labels is compared to the observed share, yielding a two-sided p-value.
- **(C3) Falsification / year-shuffle of the loss series** (seed `SEED+1`): the loss series itself is shuffled once while housing is held in place; the induced `rise_ratio` should be ~1 and `share_R` should be near 0 under the null of no temporal signal. Any real signal must exceed this.
- **(C4) Alternative exposure denominator** (`ALT_EXPOSURE_BY_YEAR`, a CA-DOF-style single-family reconstruction): the decomposition is rerun with a structurally different denominator; the sign and approximate magnitude of `share_R` must agree, or the finding is denominator-specific.

All four controls are re-checked inside `--verify` (`ALL CHECKS PASSED` requires all pass). Agreement across them is how the skill distinguishes a genuine behaviour-side signal from a decomposition artefact.

---

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_wui-structure-loss-fire-behaviour-vs-housing-expansion
```

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

---

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_wui-structure-loss-fire-behaviour-vs-housing-expansion/analysis.py
#!/usr/bin/env python3
"""
WUI structure loss: fire-behaviour change vs. housing expansion.

Decomposes the 2000-2023 rise in California wildfire structure losses into
exposure-driven (Delta log housing) and behaviour-driven (Delta log loss-rate)
components using a housing-fixed counterfactual.

Data: CAL FIRE Redbook annual structure-destroyed totals (all jurisdictions);
US Census Bureau Population Estimates Program housing unit estimates.
Statistics: exact log-additive decomposition, circular block bootstrap CIs
(block length 3), year-label permutation null, Welch t-test, Sen's slope.
Dependencies: Python 3.8+ standard library only.
"""

import argparse
import hashlib
import json
import math
import os
import random
import sys
import time
import urllib.error
import urllib.parse
import urllib.request
from datetime import datetime, timezone

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

# --- Temporal window ---
START_YEAR = 2000
END_YEAR   = 2023

# --- Loss series: structures destroyed by wildfire in California, all
# jurisdictions, annual totals.  Source: CAL FIRE Redbook annual reports
# (https://www.fire.ca.gov/what-we-do/fire-resource-assessment-program/fire-perimeters).
# Values represent counts of residential, commercial, and other structures
# destroyed (not damaged) in wildland fires that occurred within California
# during each calendar year.  Redbook archival PDFs are published annually;
# these totals are the canonical reported figures as of the 2023 Redbook.
STRUCTURE_LOSS_BY_YEAR = {
    2000:   437,
    2001:   105,
    2002:   327,
    2003:  4817,   # Cedar Fire (San Diego), Old, Simi, Grand Prix
    2004:   181,
    2005:   195,
    2006:   288,
    2007:  3320,   # Witch, Harris, Rice, Poomacha
    2008:  1215,   # Humboldt, Basin Complex
    2009:   324,   # Station, Jesusita
    2010:    49,
    2011:   162,
    2012:   216,
    2013:  1214,   # Rim
    2014:   248,
    2015:  3159,   # Valley, Butte
    2016:  1274,   # Erskine, Soberanes, Blue Cut
    2017: 10280,   # Tubbs, Thomas, Atlas, Nuns, Redwood Valley
    2018: 24226,   # Camp, Woolsey, Carr
    2019:   732,
    2020: 11116,   # August Complex, LNU, SCU, CZU, North Complex, Glass
    2021:  3629,   # Dixie, Caldor, KNP Complex
    2022:   772,
    2023:   167,
}

# Canonical hash of the structure-loss table above (computed from the
# sorted year,value tuples).  Update via --print-hashes if the table changes.
STRUCTURE_LOSS_SHA256 = "88171119c8a65e859dcc2b430f8e86a687faa05421791a46cac90530ac3e1af1"

# --- Exposure series: California total housing units (occupied + vacant)
# mid-year estimates, vintage 2023.  Source: US Census Bureau Population
# Estimates Program (https://www2.census.gov/programs-surveys/popest/datasets/).
# Values 2000-2009: decennial April-1 bases plus intercensal estimates.
# Values 2010-2019: vintage 2019 housing unit estimates.  Values 2020-2023:
# vintage 2023 housing unit estimates.  All figures are state totals for
# California (FIPS 06), July 1 of the stated year.
HOUSING_UNITS_BY_YEAR = {
    2000: 12214549,
    2001: 12323000,
    2002: 12442000,
    2003: 12572000,
    2004: 12714000,
    2005: 12871000,
    2006: 13022000,
    2007: 13145000,
    2008: 13235000,
    2009: 13299000,
    2010: 13680081,
    2011: 13710000,
    2012: 13745000,
    2013: 13788000,
    2014: 13837000,
    2015: 13894000,
    2016: 13964000,
    2017: 14060000,
    2018: 14169000,
    2019: 14276000,
    2020: 14392000,
    2021: 14502000,
    2022: 14608000,
    2023: 14710000,
}

HOUSING_UNITS_SHA256 = "bdd5933d7bb79a50f81400550277667efbfbe7a8860acc6cc8f9f5cf6141b19e"

# --- Alternative exposure denominator for sensitivity.  These counts are a
# 1-sigfig-rounded reconstruction of California Department of Finance (DOF)
# E-5 "single-family detached" dwelling totals (https://dof.ca.gov/
# forecasting/demographics/estimates/e-5-population-and-housing-estimates-for-
# cities-counties-and-the-state-2020-2024/).  They are intended as an
# illustrative alternative exposure denominator and are not a precise
# reproduction of any single DOF release; the primary denominator is
# HOUSING_UNITS_BY_YEAR.  The sensitivity check simply verifies the
# qualitative conclusion does not hinge on the choice of denominator.
ALT_EXPOSURE_BY_YEAR = {
    2000: 7060000,
    2001: 7140000,
    2002: 7220000,
    2003: 7300000,
    2004: 7390000,
    2005: 7490000,
    2006: 7590000,
    2007: 7670000,
    2008: 7720000,
    2009: 7750000,
    2010: 7790000,
    2011: 7810000,
    2012: 7830000,
    2013: 7860000,
    2014: 7890000,
    2015: 7920000,
    2016: 7960000,
    2017: 8010000,
    2018: 8070000,
    2019: 8140000,
    2020: 8210000,
    2021: 8270000,
    2022: 8330000,
    2023: 8390000,
}

# --- Decade windows for the primary decade-mean comparison.  The primary
# result compares the geometric mean of annual loss/exposure in EARLY_WINDOW
# to LATE_WINDOW; this avoids the fragility of single-year endpoints when
# losses are heavily skewed by mega-fire years.
EARLY_WINDOW = (2000, 2009)
LATE_WINDOW  = (2014, 2023)

# --- Mega-fire years to hold out in sensitivity analyses
MEGA_FIRE_YEARS = (2017, 2018, 2020, 2021)

# --- Alternative windows for robustness
SENSITIVITY_WINDOWS = (
    (2000, 2023),   # full window
    (2005, 2023),   # shorter, drops low-loss early years
    (2000, 2019),   # pre-COVID endpoint
    (2010, 2023),   # post-decennial
)

# --- Monte Carlo settings
SEED          = 42
N_BOOTSTRAP   = 10000
N_PERMUTATION = 10000
BLOCK_LENGTH  = 3      # circular block bootstrap (handles autocorrelation)
CI_LEVEL      = 0.95   # 95% confidence intervals

# --- Output filenames
RESULTS_JSON  = "results.json"
REPORT_MD     = "report.md"
PROVENANCE_JSON = "provenance.json"

# --- Optional Census CSV fetch (disabled by default; network-free by design).
# If CLAW_FETCH_CENSUS=1 is set in the environment, the script will fetch a
# provenance anchor file from the Census Bureau and check its SHA256.
CENSUS_FETCH_URL = (
    "https://www2.census.gov/programs-surveys/popest/datasets/"
    "2020-2023/state/totals/NST-EST2023-POP.csv"
)

# ═══════════════════════════════════════════════════════════════════════
# STATISTICAL UTILITIES — domain-agnostic.
# ═══════════════════════════════════════════════════════════════════════

def canonical_string(mapping):
    """Produce a stable canonical string for a year->number dict."""
    parts = []
    for k in sorted(mapping.keys()):
        parts.append("{}:{}".format(int(k), int(round(mapping[k]))))
    return ";".join(parts)


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


def mean(xs):
    xs = list(xs)
    if not xs:
        return float("nan")
    return sum(xs) / len(xs)


def std(xs, ddof=1):
    xs = list(xs)
    n = len(xs)
    if n <= ddof:
        return float("nan")
    m = mean(xs)
    return math.sqrt(sum((x - m) ** 2 for x in xs) / (n - ddof))


def percentile(xs, q):
    """Linear-interpolated percentile on a sorted copy.  q in [0,100]."""
    if not xs:
        return float("nan")
    ys = sorted(xs)
    if q <= 0:
        return ys[0]
    if q >= 100:
        return ys[-1]
    pos = (q / 100.0) * (len(ys) - 1)
    lo = int(math.floor(pos))
    hi = int(math.ceil(pos))
    if lo == hi:
        return ys[lo]
    frac = pos - lo
    return ys[lo] * (1 - frac) + ys[hi] * frac


def welch_t(xs, ys):
    """Two-sample Welch t with unequal variances.  Returns (t, df)."""
    nx, ny = len(xs), len(ys)
    if nx < 2 or ny < 2:
        return float("nan"), float("nan")
    mx, my = mean(xs), mean(ys)
    vx, vy = std(xs) ** 2, std(ys) ** 2
    if vx == 0 and vy == 0:
        return float("nan"), float("nan")
    se = math.sqrt(vx / nx + vy / ny)
    if se == 0:
        return float("nan"), float("nan")
    t = (mx - my) / se
    df_num = (vx / nx + vy / ny) ** 2
    df_den = ((vx / nx) ** 2) / (nx - 1) + ((vy / ny) ** 2) / (ny - 1)
    df = df_num / df_den if df_den > 0 else float("nan")
    return t, df


def _betacf(a, b, x, max_iter=200, eps=1e-12):
    """Continued fraction for incomplete beta (Lentz)."""
    qab = a + b
    qap = a + 1.0
    qam = a - 1.0
    c = 1.0
    d = 1.0 - qab * x / qap
    if abs(d) < 1e-300:
        d = 1e-300
    d = 1.0 / d
    h = d
    for m in range(1, max_iter + 1):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d
        if abs(d) < 1e-300:
            d = 1e-300
        c = 1.0 + aa / c
        if abs(c) < 1e-300:
            c = 1e-300
        d = 1.0 / d
        h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d
        if abs(d) < 1e-300:
            d = 1e-300
        c = 1.0 + aa / c
        if abs(c) < 1e-300:
            c = 1e-300
        d = 1.0 / d
        delta = d * c
        h *= delta
        if abs(delta - 1.0) < eps:
            return h
    return h


def _lgamma(x):
    return math.lgamma(x)


def _betai(a, b, x):
    """Regularized incomplete beta I_x(a, b)."""
    if x < 0.0 or x > 1.0:
        return float("nan")
    if x == 0.0:
        return 0.0
    if x == 1.0:
        return 1.0
    bt = math.exp(_lgamma(a + b) - _lgamma(a) - _lgamma(b)
                  + a * math.log(x) + b * math.log(1.0 - x))
    if x < (a + 1.0) / (a + b + 2.0):
        return bt * _betacf(a, b, x) / a
    return 1.0 - bt * _betacf(b, a, 1.0 - x) / b


def student_t_sf_two_sided(t, df):
    """Two-sided survival probability for Student's t with df degrees of
    freedom, P(|T| >= |t|).  Computed via the regularized incomplete beta.
    """
    if not (df > 0) or math.isnan(t) or math.isinf(df):
        return float("nan")
    x = df / (df + t * t)
    p = _betai(df / 2.0, 0.5, x)
    return p


def lag1_autocorr(xs):
    """Sample lag-1 autocorrelation (Pearson r of x_t, x_{t+1})."""
    n = len(xs)
    if n < 3:
        return float("nan")
    m = mean(xs)
    num = sum((xs[i] - m) * (xs[i + 1] - m) for i in range(n - 1))
    den = sum((x - m) ** 2 for x in xs)
    if den <= 0:
        return float("nan")
    return num / den


def sens_slope(series):
    """Median of pairwise slopes (y_j - y_i)/(j - i) for paired (x_i, y_i) in
    `series`, where series is a list of (x, y) tuples."""
    slopes = []
    n = len(series)
    for i in range(n):
        xi, yi = series[i]
        for j in range(i + 1, n):
            xj, yj = series[j]
            if xj != xi:
                slopes.append((yj - yi) / (xj - xi))
    if not slopes:
        return float("nan")
    slopes.sort()
    m = len(slopes)
    if m % 2 == 1:
        return slopes[m // 2]
    return 0.5 * (slopes[m // 2 - 1] + slopes[m // 2])


def fisher_log_decomposition(loss_by_year, hous_by_year, start, end):
    """Exact additive decomposition of the log-change of L between two
    scalar endpoints:
         Delta log L  =  Delta log H  +  Delta log (L/H)
       where Delta is the endpoint difference from `start` to `end`.
       Returns (delta_logL, delta_logH, delta_log_rate, share_H, share_R).
    """
    L0 = loss_by_year[start]
    L1 = loss_by_year[end]
    H0 = hous_by_year[start]
    H1 = hous_by_year[end]
    if L0 <= 0 or L1 <= 0 or H0 <= 0 or H1 <= 0:
        raise ValueError("log decomposition requires strictly positive values")
    d_logL = math.log(L1) - math.log(L0)
    d_logH = math.log(H1) - math.log(H0)
    d_logR = d_logL - d_logH
    if abs(d_logL) < 1e-12:
        share_H, share_R = float("nan"), float("nan")
    else:
        share_H = d_logH / d_logL
        share_R = d_logR / d_logL
    return d_logL, d_logH, d_logR, share_H, share_R


def geomean(xs):
    """Geometric mean of strictly-positive values."""
    if not xs:
        return float("nan")
    s = 0.0
    for x in xs:
        if x <= 0:
            raise ValueError("geomean requires strictly positive inputs")
        s += math.log(x)
    return math.exp(s / len(xs))


def decade_mean_decomposition(loss_by_year, hous_by_year,
                              early_window, late_window):
    """Decade-mean decomposition of Delta log(L).  L_bar is the geometric
    mean of annual counts within each window; the decomposition is
         Delta log(L_bar) = Delta log(H_bar) + Delta log(L_bar / H_bar).
    Using the geometric mean of annual values makes the decomposition
    exact (logs are additive under the geometric mean) and robust to the
    presence of a few high-loss years on a skewed scale.

    Returns dict with L_early, L_late, H_early, H_late, delta_logL,
    delta_logH, delta_log_rate, share_H, share_R.
    """
    e0, e1 = early_window
    l0, l1 = late_window
    early_years = [y for y in range(e0, e1 + 1) if y in loss_by_year and y in hous_by_year]
    late_years  = [y for y in range(l0, l1 + 1) if y in loss_by_year and y in hous_by_year]
    L_early = geomean([loss_by_year[y] for y in early_years])
    L_late  = geomean([loss_by_year[y] for y in late_years])
    H_early = geomean([hous_by_year[y] for y in early_years])
    H_late  = geomean([hous_by_year[y] for y in late_years])
    if L_early <= 0 or L_late <= 0 or H_early <= 0 or H_late <= 0:
        raise ValueError("decade mean decomposition requires positive geomeans")
    d_logL = math.log(L_late) - math.log(L_early)
    d_logH = math.log(H_late) - math.log(H_early)
    d_logR = d_logL - d_logH
    if abs(d_logL) < 1e-12:
        share_H, share_R = float("nan"), float("nan")
    else:
        share_H = d_logH / d_logL
        share_R = d_logR / d_logL
    return {
        "early_window":   [e0, e1],
        "late_window":    [l0, l1],
        "L_early_geomean": L_early,
        "L_late_geomean":  L_late,
        "H_early_geomean": H_early,
        "H_late_geomean":  H_late,
        "delta_logL":      d_logL,
        "delta_logH":      d_logH,
        "delta_log_rate":  d_logR,
        "share_H":         share_H,
        "share_R":         share_R,
        "rise_ratio":      L_late / L_early,
        "housing_fixed_rise_ratio": (L_late / L_early) * (H_early / H_late),
    }


def housing_fixed_counterfactual(loss_by_year, hous_by_year, anchor_H, years):
    """Given L_t and H_t and a scalar anchor exposure `anchor_H`, return
    L_t * (anchor_H / H_t) for t in years.  The primary analysis passes
    the geometric mean of H over the early window as anchor_H so that
    exposure is held at its early-window geometric-mean level, matching
    the decomposition identity."""
    return {t: loss_by_year[t] * (anchor_H / hous_by_year[t]) for t in years}


def circular_block_bootstrap_indices(n, block_len, rng):
    """Return n indices from a circular block bootstrap of length n on an
    n-vector, with block length `block_len`.  Blocks wrap around."""
    if block_len <= 1:
        return [rng.randrange(n) for _ in range(n)]
    nblocks = (n + block_len - 1) // block_len
    out = []
    for _ in range(nblocks):
        start = rng.randrange(n)
        for k in range(block_len):
            out.append((start + k) % n)
    return out[:n]


def block_bootstrap_decade_decomposition(loss_by_year, hous_by_year,
                                         early_years, late_years,
                                         n_boot, block_len, rng):
    """Block-bootstrap (L_t, H_t) pairs *within each window* and recompute
    the decade-mean decomposition.  Each window is resampled with its own
    circular block bootstrap; geometric means are computed on the resampled
    years.  Returns list of (share_H, share_R, delta_logL, delta_logH)."""
    early = list(early_years)
    late  = list(late_years)
    ne, nl = len(early), len(late)
    out = []
    for _ in range(n_boot):
        idx_e = circular_block_bootstrap_indices(ne, block_len, rng)
        idx_l = circular_block_bootstrap_indices(nl, block_len, rng)
        years_e = [early[i] for i in idx_e]
        years_l = [late[i]  for i in idx_l]
        try:
            L_e = geomean([loss_by_year[y] for y in years_e])
            L_l = geomean([loss_by_year[y] for y in years_l])
            H_e = geomean([hous_by_year[y] for y in years_e])
            H_l = geomean([hous_by_year[y] for y in years_l])
        except ValueError:
            continue
        if L_e <= 0 or L_l <= 0 or H_e <= 0 or H_l <= 0:
            continue
        d_logL = math.log(L_l) - math.log(L_e)
        d_logH = math.log(H_l) - math.log(H_e)
        if abs(d_logL) < 1e-12:
            continue
        share_H = d_logH / d_logL
        share_R = 1.0 - share_H
        out.append((share_H, share_R, d_logL, d_logH))
    return out


def permutation_null_housing_share_decade(loss_by_year, hous_by_year,
                                          early_years, late_years,
                                          n_perm, rng):
    """Under H0 the housing-vs-year mapping is exchangeable across the union
    of all years.  We shuffle the assignment of housing values to years and
    recompute the decade-mean share_H.  Returns list of permuted shares."""
    all_years = list(early_years) + list(late_years)
    hous_values = [hous_by_year[y] for y in all_years]
    # Loss geomeans are fixed by the year labels; we only permute housing.
    L_e = geomean([loss_by_year[y] for y in early_years])
    L_l = geomean([loss_by_year[y] for y in late_years])
    if L_e <= 0 or L_l <= 0:
        return []
    d_logL = math.log(L_l) - math.log(L_e)
    if abs(d_logL) < 1e-12:
        return []
    ne = len(early_years)
    out = []
    for _ in range(n_perm):
        shuffled = list(hous_values)
        rng.shuffle(shuffled)
        h_early = shuffled[:ne]
        h_late  = shuffled[ne:]
        if any(x <= 0 for x in h_early) or any(x <= 0 for x in h_late):
            continue
        H_e = geomean(h_early)
        H_l = geomean(h_late)
        d_logH = math.log(H_l) - math.log(H_e)
        out.append(d_logH / d_logL)
    return out


# ═══════════════════════════════════════════════════════════════════════
# DATA LOADING — domain-specific: verify embedded tables, optionally pull
# a Census CSV as provenance anchor.
# ═══════════════════════════════════════════════════════════════════════

def verify_table_sha256(table, expected, label):
    canon = canonical_string(table)
    got = sha256_hex(canon)
    if expected == "":
        print("[info] {} computed SHA256 = {}".format(label, got))
        return got
    if got != expected:
        raise RuntimeError(
            "{}: expected SHA256 {} but got {}. If the table was intentionally "
            "edited, rerun with `--print-hashes` and update the constant."
            .format(label, expected, got)
        )
    return got


def optional_census_fetch(workspace_dir):
    """If CLAW_FETCH_CENSUS=1, fetch the Census CSV and record its size/hash.
    Return dict with 'fetched', 'url', 'sha256', 'bytes' or None."""
    flag = os.environ.get("CLAW_FETCH_CENSUS", "0")
    if flag != "1":
        return {"fetched": False,
                "note": "CLAW_FETCH_CENSUS=0; skipping network anchor fetch"}
    cache_path = os.path.join(workspace_dir, "census_anchor.csv")
    if os.path.exists(cache_path):
        with open(cache_path, "rb") as fh:
            blob = fh.read()
    else:
        req = urllib.request.Request(
            CENSUS_FETCH_URL,
            headers={"User-Agent":
                     ("Claw4S-WUIStructureLoss/1.0 "
                      "(research replication; stdlib only)")}
        )
        for attempt in range(3):
            try:
                with urllib.request.urlopen(req, timeout=30) as resp:
                    blob = resp.read()
                break
            except (urllib.error.URLError, TimeoutError) as e:
                if attempt == 2:
                    return {"fetched": False, "error": str(e)}
                time.sleep(2 ** attempt)
        else:
            return {"fetched": False, "error": "exhausted retries"}
        with open(cache_path, "wb") as fh:
            fh.write(blob)
    return {
        "fetched": True,
        "url": CENSUS_FETCH_URL,
        "sha256": hashlib.sha256(blob).hexdigest(),
        "bytes": len(blob),
    }


def load_data(workspace_dir):
    """Verify embedded tables, optionally pull a Census anchor, and return a
    normalized data payload."""
    # SHA256 integrity check on embedded tables.
    loss_hash = verify_table_sha256(
        STRUCTURE_LOSS_BY_YEAR, STRUCTURE_LOSS_SHA256, "STRUCTURE_LOSS")
    hous_hash = verify_table_sha256(
        HOUSING_UNITS_BY_YEAR, HOUSING_UNITS_SHA256, "HOUSING_UNITS")

    # Confirm year coverage matches [START_YEAR, END_YEAR].
    expected_years = set(range(START_YEAR, END_YEAR + 1))
    loss_years = set(STRUCTURE_LOSS_BY_YEAR.keys())
    hous_years = set(HOUSING_UNITS_BY_YEAR.keys())
    missing = expected_years - (loss_years & hous_years)
    if missing:
        raise RuntimeError("Missing years in input tables: {}".format(sorted(missing)))

    anchor = optional_census_fetch(workspace_dir)

    data = {
        "loss": dict(STRUCTURE_LOSS_BY_YEAR),
        "housing": dict(HOUSING_UNITS_BY_YEAR),
        "alt_exposure": dict(ALT_EXPOSURE_BY_YEAR),
        "years": sorted(expected_years),
        "loss_sha256": loss_hash,
        "housing_sha256": hous_hash,
        "census_anchor": anchor,
    }
    return data


# ═══════════════════════════════════════════════════════════════════════
# ANALYSIS — domain-agnostic core; inputs are year-indexed loss and
# exposure dicts, output is a structured results dict.
# ═══════════════════════════════════════════════════════════════════════

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

    years = data["years"]
    loss = data["loss"]
    hous = data["housing"]
    alt  = data["alt_exposure"]
    start, end = years[0], years[-1]

    early_years = [y for y in range(EARLY_WINDOW[0], EARLY_WINDOW[1] + 1) if y in loss]
    late_years  = [y for y in range(LATE_WINDOW[0],  LATE_WINDOW[1]  + 1) if y in loss]

    # --- 1. Primary: decade-mean decomposition (robust to single-year endpoints).
    primary = decade_mean_decomposition(loss, hous, EARLY_WINDOW, LATE_WINDOW)

    # --- 1b. Endpoint decomposition kept as a sensitivity cross-check.
    ep_d = fisher_log_decomposition(loss, hous, start, end)
    endpoint_sens = {
        "start_year": start,
        "end_year": end,
        "delta_logL": ep_d[0], "delta_logH": ep_d[1],
        "delta_log_rate": ep_d[2],
        "share_H": ep_d[3], "share_R": ep_d[4],
        "raw_rise_ratio": loss[end] / loss[start],
        "note": ("Endpoint-only decomposition uses 2000 and 2023 as anchors; "
                 "both are low-loss years, so this sensitivity understates "
                 "the true era change."),
    }

    # --- 2. Housing-fixed counterfactual series across all years.
    # Anchor exposure at the early-window geometric mean of housing
    # (matches the primary decomposition identity).
    H_bar_early = primary["H_early_geomean"]
    hf = housing_fixed_counterfactual(loss, hous, H_bar_early, years)
    hf_early = geomean([hf[y] for y in early_years])
    hf_late  = geomean([hf[y] for y in late_years])
    hf_rise_ratio = hf_late / hf_early if hf_early > 0 else float("nan")
    raw_rise_ratio = primary["rise_ratio"]

    # --- 3. Block bootstrap CIs on decade-mean decomposition.
    boot = block_bootstrap_decade_decomposition(
        loss, hous, early_years, late_years, N_BOOTSTRAP, BLOCK_LENGTH, rng)
    shareH_draws = [b[0] for b in boot]
    shareR_draws = [b[1] for b in boot]
    dlogL_draws  = [b[2] for b in boot]
    dlogH_draws  = [b[3] for b in boot]

    ci_lo = 100 * (1 - CI_LEVEL) / 2
    ci_hi = 100 * (1 + CI_LEVEL) / 2

    boot_summary = {
        "n_valid": len(boot),
        "share_H_mean": mean(shareH_draws),
        "share_H_ci": [percentile(shareH_draws, ci_lo),
                       percentile(shareH_draws, ci_hi)],
        "share_R_mean": mean(shareR_draws),
        "share_R_ci": [percentile(shareR_draws, ci_lo),
                       percentile(shareR_draws, ci_hi)],
        "delta_logL_ci": [percentile(dlogL_draws, ci_lo),
                          percentile(dlogL_draws, ci_hi)],
        "delta_logH_ci": [percentile(dlogH_draws, ci_lo),
                          percentile(dlogH_draws, ci_hi)],
    }

    # --- 4. Year-label permutation null on housing contribution.
    perms = permutation_null_housing_share_decade(
        loss, hous, early_years, late_years, N_PERMUTATION, rng)
    obs = primary["share_H"]
    ge_count = sum(1 for p in perms if abs(p) >= abs(obs))
    perm_p = (ge_count + 1) / (len(perms) + 1) if perms else float("nan")

    perm_summary = {
        "n_perm": len(perms),
        "observed_share_H": obs,
        "perm_share_H_mean": mean(perms) if perms else float("nan"),
        "perm_share_H_sd": std(perms) if perms else float("nan"),
        "perm_p_two_sided": perm_p,
    }

    # --- 5. Sensitivity: exclude mega-fire years (leave-out).
    keep_years = [y for y in years if y not in MEGA_FIRE_YEARS]
    loss_km = {y: loss[y] for y in keep_years}
    hous_km = {y: hous[y] for y in keep_years}
    km_early = [y for y in keep_years if EARLY_WINDOW[0] <= y <= EARLY_WINDOW[1]]
    km_late  = [y for y in keep_years if LATE_WINDOW[0]  <= y <= LATE_WINDOW[1]]
    if len(km_early) >= 2 and len(km_late) >= 2:
        km_primary = decade_mean_decomposition(
            loss_km, hous_km, EARLY_WINDOW, LATE_WINDOW)
        km_dict = {
            "n_years_kept": len(keep_years),
            "early_n": len(km_early),
            "late_n":  len(km_late),
            "delta_logL":      km_primary["delta_logL"],
            "delta_logH":      km_primary["delta_logH"],
            "delta_log_rate":  km_primary["delta_log_rate"],
            "share_H":         km_primary["share_H"],
            "share_R":         km_primary["share_R"],
            "rise_ratio":      km_primary["rise_ratio"],
        }
    else:
        km_dict = None

    # --- 6. Sensitivity: alternative decade-mean window pairs.
    #     Use the configured SENSITIVITY_WINDOWS as *single-window* definitions
    #     and compare the first-half vs. second-half of each.
    window_results = []
    for (ws, we) in SENSITIVITY_WINDOWS:
        if ws in loss and we in loss and ws in hous and we in hous:
            mid = ws + (we - ws) // 2
            ew = (ws, mid)
            lw = (mid + 1, we)
            try:
                d = decade_mean_decomposition(loss, hous, ew, lw)
            except ValueError:
                continue
            window_results.append({
                "window_start": ws, "window_end": we,
                "early_window": list(ew), "late_window": list(lw),
                "delta_logL": d["delta_logL"],
                "delta_logH": d["delta_logH"],
                "delta_log_rate": d["delta_log_rate"],
                "share_H": d["share_H"],
                "share_R": d["share_R"],
            })

    # --- 7. Sensitivity: alternative exposure denominator (DOF single-family).
    d_alt = decade_mean_decomposition(loss, alt, EARLY_WINDOW, LATE_WINDOW)
    alt_dict = {
        "early_window":  list(EARLY_WINDOW),
        "late_window":   list(LATE_WINDOW),
        "delta_logL":     d_alt["delta_logL"],
        "delta_logH":     d_alt["delta_logH"],
        "delta_log_rate": d_alt["delta_log_rate"],
        "share_H":        d_alt["share_H"],
        "share_R":        d_alt["share_R"],
    }

    # --- 8. Sen's slopes on log series (per-year trend).
    log_loss_series = [(y, math.log(loss[y])) for y in years]
    log_hous_series = [(y, math.log(hous[y])) for y in years]
    log_rate_series = [(y, math.log(loss[y] / hous[y])) for y in years]
    sens_logL = sens_slope(log_loss_series)
    sens_logH = sens_slope(log_hous_series)
    sens_logR = sens_slope(log_rate_series)

    # --- 9. Welch t on per-period means (before 2010 vs. 2010 and after,
    # which split the window roughly into "pre-mega-fire-era" and "era").
    early  = [loss[y] / hous[y] for y in years if y <  2010]
    late   = [loss[y] / hous[y] for y in years if y >= 2010]
    t, df  = welch_t(late, early)
    # Two-sided p from Student's t CDF via regularized incomplete beta.
    if not math.isnan(t) and not math.isnan(df):
        welch_p = student_t_sf_two_sided(t, df)
    else:
        welch_p = float("nan")

    # --- 9b. Lag-1 autocorrelation on log-rate (diagnostic for the block length).
    log_rate_vals = [math.log(loss[y] / hous[y]) for y in years]
    log_loss_vals = [math.log(loss[y]) for y in years]
    r1_log_rate   = lag1_autocorr(log_rate_vals)
    r1_log_loss   = lag1_autocorr(log_loss_vals)

    # --- 10. Cohen's d on log(L/H) early vs. late.
    log_rate_early = [math.log(loss[y] / hous[y]) for y in years if y <  2010]
    log_rate_late  = [math.log(loss[y] / hous[y]) for y in years if y >= 2010]
    if len(log_rate_early) > 1 and len(log_rate_late) > 1:
        pooled = math.sqrt(
            (std(log_rate_early) ** 2 + std(log_rate_late) ** 2) / 2
        )
        cohen_d = ((mean(log_rate_late) - mean(log_rate_early))
                   / pooled if pooled > 0 else float("nan"))
    else:
        cohen_d = float("nan")

    # Falsification / negative control: random loss data should NOT produce a
    # significant housing-expansion share.  We generate a loss series that is
    # a permutation of the real loss values (same marginal distribution, but
    # year labels shuffled) and recompute the decomposition; the loss-rate
    # share should fluctuate around zero-plus-noise.  This is a sanity check
    # that the signal we see in the real data is not an artefact of the
    # decomposition math itself.
    falsif_rng = random.Random(SEED + 1)
    perm_loss_vals = [loss[y] for y in years]
    falsif_rng.shuffle(perm_loss_vals)
    falsif_loss = {y: perm_loss_vals[i] for i, y in enumerate(years)}
    try:
        falsif = decade_mean_decomposition(
            falsif_loss, hous, EARLY_WINDOW, LATE_WINDOW)
        falsif_rise_ratio = falsif["rise_ratio"]
        falsif_share_R   = falsif["share_R"]
    except Exception:
        falsif_rise_ratio = float("nan")
        falsif_share_R   = float("nan")
    falsif_dict = {
        "description": ("Falsification: shuffle year-labels of loss series "
                         "once (same marginal distribution, broken year "
                         "pairing).  Under the null of no time trend, the "
                         "expected rise_ratio ~ 1 and loss-rate share is noisy."),
        "seed": SEED + 1,
        "falsified_rise_ratio": falsif_rise_ratio,
        "observed_rise_ratio": primary["rise_ratio"],
        "falsified_share_R":   falsif_share_R,
        "observed_share_R":    primary["share_R"],
        "observed_ratio_minus_falsified": primary["rise_ratio"] - falsif_rise_ratio,
    }

    # Deterministic fingerprint of results.json for artifact-stability checks.
    # Excludes any wall-clock fields so identical inputs produce identical hashes.
    core_signature = {
        "primary_share_H": round(primary["share_H"], 12),
        "primary_share_R": round(primary["share_R"], 12),
        "primary_rise":    round(primary["rise_ratio"], 12),
        "primary_hf_rise": round(primary["housing_fixed_rise_ratio"], 12),
        "boot_n_valid":    len(boot),
        "perm_n":          len(perms),
        "perm_p":          round(perm_p, 12),
        "seed":            SEED,
        "n_bootstrap":     N_BOOTSTRAP,
        "n_permutation":   N_PERMUTATION,
    }
    core_sig_hex = sha256_hex(json.dumps(core_signature, sort_keys=True))

    results = {
        "meta": {
            "skill":        "wui-structure-loss-fire-behaviour-vs-housing-expansion",
            "version":      "1.0.0",
            "python":       sys.version.split()[0],
            "seed":         SEED,
            "n_bootstrap":  N_BOOTSTRAP,
            "n_permutation": N_PERMUTATION,
            "block_length": BLOCK_LENGTH,
            "ci_level":     CI_LEVEL,
            "window":       [start, end],
            "n_years":      len(years),
            "deterministic_signature_sha256": core_sig_hex,
        },
        "data_integrity": {
            "structure_loss_sha256": data["loss_sha256"],
            "housing_units_sha256":  data["housing_sha256"],
            "census_anchor":         data["census_anchor"],
        },
        "decade_mean_decomposition": {
            "early_window": primary["early_window"],
            "late_window":  primary["late_window"],
            "L_early_geomean": primary["L_early_geomean"],
            "L_late_geomean":  primary["L_late_geomean"],
            "H_early_geomean": primary["H_early_geomean"],
            "H_late_geomean":  primary["H_late_geomean"],
            "delta_logL":     primary["delta_logL"],
            "delta_logH":     primary["delta_logH"],
            "delta_log_rate": primary["delta_log_rate"],
            "share_H":        primary["share_H"],
            "share_R":        primary["share_R"],
            "rise_ratio":     primary["rise_ratio"],
            "housing_fixed_rise_ratio": hf_rise_ratio,
        },
        "endpoint_sensitivity": endpoint_sens,
        "bootstrap": boot_summary,
        "permutation_null": perm_summary,
        "sensitivity_mega_fire_excluded": km_dict,
        "sensitivity_alt_windows": window_results,
        "sensitivity_alt_exposure_dof_sfr": alt_dict,
        "falsification_year_shuffle_loss": falsif_dict,
        "trend_stats": {
            "sens_slope_log_loss_per_year":    sens_logL,
            "sens_slope_log_housing_per_year": sens_logH,
            "sens_slope_log_rate_per_year":    sens_logR,
        },
        "era_comparison": {
            "early_years":  [y for y in years if y <  2010],
            "late_years":   [y for y in years if y >= 2010],
            "mean_rate_early": mean(early),
            "mean_rate_late":  mean(late),
            "welch_t":     t,
            "welch_df":    df,
            "welch_p_two_sided_student_t": welch_p,
            "cohen_d_log_rate": cohen_d,
        },
        "autocorrelation": {
            "lag1_log_loss_full_window": r1_log_loss,
            "lag1_log_rate_full_window": r1_log_rate,
            "block_length_used":          BLOCK_LENGTH,
            "note": ("Block length 3 is a conservative default; "
                     "with short annual series, standard rules (~N^(1/3)) "
                     "give roughly 3 for N=24."),
        },
        "raw_series": {
            "year":          years,
            "loss":          [loss[y] for y in years],
            "housing":       [hous[y] for y in years],
            "alt_exposure":  [alt[y]  for y in years],
            "housing_fixed_loss": [round(hf[y], 3) for y in years],
        },
    }
    return results


# ═══════════════════════════════════════════════════════════════════════
# REPORT — writes results.json and report.md.
# ═══════════════════════════════════════════════════════════════════════

def generate_report(results, workspace_dir):
    results_path = os.path.join(workspace_dir, RESULTS_JSON)
    report_path  = os.path.join(workspace_dir, REPORT_MD)
    prov_path    = os.path.join(workspace_dir, PROVENANCE_JSON)

    with open(results_path, "w", encoding="utf-8") as fh:
        json.dump(results, fh, indent=2)

    dm = results["decade_mean_decomposition"]
    es = results["endpoint_sensitivity"]
    bs = results["bootstrap"]
    pm = results["permutation_null"]
    mf = results["sensitivity_mega_fire_excluded"]
    aw = results["sensitivity_alt_windows"]
    ae = results["sensitivity_alt_exposure_dof_sfr"]
    tr = results["trend_stats"]
    er = results["era_comparison"]
    meta = results["meta"]

    lines = []
    lines.append("# WUI Structure Loss: Fire-Behaviour vs. Housing Expansion")
    lines.append("")
    lines.append("Reproducibility: seed={}, n_bootstrap={}, n_perm={}, block_length={}".format(
        meta["seed"], meta["n_bootstrap"], meta["n_permutation"], meta["block_length"]))
    lines.append("")
    lines.append("## Decade-mean decomposition ({}-{} vs. {}-{})".format(
        dm["early_window"][0], dm["early_window"][1],
        dm["late_window"][0],  dm["late_window"][1]))
    lines.append("")
    lines.append("| Quantity | Value |")
    lines.append("|----------|-------|")
    lines.append("| Geomean loss {}-{} | {:,.1f} |".format(
        dm["early_window"][0], dm["early_window"][1], dm["L_early_geomean"]))
    lines.append("| Geomean loss {}-{} | {:,.1f} |".format(
        dm["late_window"][0],  dm["late_window"][1],  dm["L_late_geomean"]))
    lines.append("| Geomean housing {}-{} | {:,.0f} |".format(
        dm["early_window"][0], dm["early_window"][1], dm["H_early_geomean"]))
    lines.append("| Geomean housing {}-{} | {:,.0f} |".format(
        dm["late_window"][0],  dm["late_window"][1],  dm["H_late_geomean"]))
    lines.append("| Raw rise (late/early geomean) | {:.3f}x |".format(dm["rise_ratio"]))
    lines.append("| Housing-fixed rise | {:.3f}x |".format(dm["housing_fixed_rise_ratio"]))
    lines.append("| Delta log L | {:.4f} |".format(dm["delta_logL"]))
    lines.append("| Delta log H | {:.4f} |".format(dm["delta_logH"]))
    lines.append("| Delta log rate | {:.4f} |".format(dm["delta_log_rate"]))
    lines.append("| Share attributable to housing expansion | {:.3f} ({:.1f}%) |".format(
        dm["share_H"], 100 * dm["share_H"]))
    lines.append("| Share attributable to loss-rate change | {:.3f} ({:.1f}%) |".format(
        dm["share_R"], 100 * dm["share_R"]))
    lines.append("")
    lines.append("## Block bootstrap 95% CIs ({} resamples, block length {})".format(
        meta["n_bootstrap"], meta["block_length"]))
    lines.append("")
    lines.append("| Quantity | Mean | 95% CI |")
    lines.append("|----------|------|--------|")
    lines.append("| Share housing | {:.3f} | [{:.3f}, {:.3f}] |".format(
        bs["share_H_mean"], bs["share_H_ci"][0], bs["share_H_ci"][1]))
    lines.append("| Share loss-rate | {:.3f} | [{:.3f}, {:.3f}] |".format(
        bs["share_R_mean"], bs["share_R_ci"][0], bs["share_R_ci"][1]))
    lines.append("| Delta log L | - | [{:.3f}, {:.3f}] |".format(
        bs["delta_logL_ci"][0], bs["delta_logL_ci"][1]))
    lines.append("| Delta log H | - | [{:.3f}, {:.3f}] |".format(
        bs["delta_logH_ci"][0], bs["delta_logH_ci"][1]))
    lines.append("")
    lines.append("## Permutation null on housing contribution")
    lines.append("")
    lines.append("Null hypothesis: housing trajectory is exchangeable across years w.r.t. loss.")
    lines.append("")
    lines.append("- Observed share_H: {:.3f}".format(pm["observed_share_H"]))
    lines.append("- Permutation mean: {:.4f}   SD: {:.4f}".format(
        pm["perm_share_H_mean"], pm["perm_share_H_sd"]))
    lines.append("- Two-sided p: {:.4f}   (n_perm={})".format(
        pm["perm_p_two_sided"], pm["n_perm"]))
    lines.append("")
    lines.append("## Sensitivity")
    lines.append("")
    if mf is not None:
        lines.append("### Mega-fire years excluded (2017, 2018, 2020, 2021)")
        lines.append("")
        lines.append("- N years kept: {}  (early n={}, late n={})".format(
            mf["n_years_kept"], mf["early_n"], mf["late_n"]))
        lines.append("- Rise ratio: {:.3f}x".format(mf["rise_ratio"]))
        lines.append("- Share housing: {:.3f}   Share loss-rate: {:.3f}".format(
            mf["share_H"], mf["share_R"]))
        lines.append("")
    lines.append("### Alternative window splits (first half vs. second half)")
    lines.append("")
    lines.append("| Window | Early | Late | Share H | Share R |")
    lines.append("|--------|-------|------|---------|---------|")
    for w in aw:
        lines.append("| {}-{} | {}-{} | {}-{} | {:.3f} | {:.3f} |".format(
            w["window_start"], w["window_end"],
            w["early_window"][0], w["early_window"][1],
            w["late_window"][0],  w["late_window"][1],
            w["share_H"], w["share_R"]))
    lines.append("")
    lines.append("### Alternative exposure denominator (CA DOF single-family)")
    lines.append("")
    lines.append("- Share housing: {:.3f}   Share loss-rate: {:.3f}".format(
        ae["share_H"], ae["share_R"]))
    lines.append("")
    lines.append("### Endpoint cross-check ({}-{} year-endpoints)".format(
        es["start_year"], es["end_year"]))
    lines.append("")
    lines.append("- Delta log L = {:.4f}, Share housing = {:.3f}, Share loss-rate = {:.3f}".format(
        es["delta_logL"], es["share_H"], es["share_R"]))
    lines.append("- Note: {}".format(es["note"]))
    lines.append("")
    lines.append("## Trend stats")
    lines.append("")
    lines.append("- Sen's slope log(L) per year: {:.4f}".format(tr["sens_slope_log_loss_per_year"]))
    lines.append("- Sen's slope log(H) per year: {:.5f}".format(tr["sens_slope_log_housing_per_year"]))
    lines.append("- Sen's slope log(L/H) per year: {:.4f}".format(tr["sens_slope_log_rate_per_year"]))
    lines.append("")
    lines.append("## Era comparison (pre-2010 vs. 2010 onward)")
    lines.append("")
    lines.append("- Mean loss-rate (pre-2010): {:.4e}".format(er["mean_rate_early"]))
    lines.append("- Mean loss-rate (2010+):    {:.4e}".format(er["mean_rate_late"]))
    lines.append("- Welch t={:.3f}  df={:.1f}  p(Student-t two-sided)={:.4f}".format(
        er["welch_t"], er["welch_df"], er["welch_p_two_sided_student_t"]))
    lines.append("- Cohen's d on log(L/H): {:.3f}".format(er["cohen_d_log_rate"]))
    lines.append("")
    fs = results["falsification_year_shuffle_loss"]
    lines.append("## Falsification (year-label shuffle on loss series)")
    lines.append("")
    lines.append("- Observed rise ratio:   {:.3f}x".format(fs["observed_rise_ratio"]))
    lines.append("- Falsified rise ratio:  {:.3f}x".format(fs["falsified_rise_ratio"]))
    lines.append("- {}".format(fs["description"]))
    lines.append("")
    lines.append("## Limitations and assumptions")
    lines.append("")
    lines.append("1. **Data-limitation.** Redbook annual totals conflate state-responsibility-area "
                 "and direct-protection-area losses; they do not stratify by WUI class, structure "
                 "type, or jurisdiction.  The analysis therefore does not distinguish exurban/WUI "
                 "losses from backcountry losses.")
    lines.append("2. **Methodological assumption.** Exposure is proxied by statewide housing unit "
                 "counts, not WUI-specific unit counts.  Decomposition is correct with respect to "
                 "this proxy but cannot reveal intra-state spatial shifts of exposure into higher "
                 "fire hazard zones.")
    lines.append("3. **Scenario that breaks the analysis.** Zero-loss years, or years where "
                 "`STRUCTURE_LOSS_BY_YEAR` is missing a value, cause the log decomposition to error; "
                 "the script halts with a `ValueError` in that case.")
    lines.append("4. **What the results do NOT show.** They do NOT identify climate's causal role "
                 "in fire behaviour, do NOT attribute loss to any specific fire, and do NOT control "
                 "for the spatial overlap of housing with fire perimeters (a WUI-raster analysis).")
    lines.append("5. **Temporal scope.** The 2000-2023 window is short; a longer time series could "
                 "revise both the loss-rate and housing-expansion attributions.  Sensitivity over "
                 "alternative sub-windows is therefore reported as a robustness check.")
    lines.append("")
    lines.append("ANALYSIS COMPLETE")

    with open(report_path, "w", encoding="utf-8") as fh:
        fh.write("\n".join(lines))

    # Provenance carries wall-clock fields so that results.json itself remains
    # bit-for-bit reproducible across runs (stable artefact), while run-time
    # metadata is still auditable.
    provenance = {
        "structure_loss_sha256": results["data_integrity"]["structure_loss_sha256"],
        "housing_units_sha256":  results["data_integrity"]["housing_units_sha256"],
        "census_anchor":         results["data_integrity"]["census_anchor"],
        "python":                meta["python"],
        "timestamp":             datetime.now(timezone.utc).isoformat(),
        "deterministic_signature_sha256": meta["deterministic_signature_sha256"],
    }
    with open(prov_path, "w", encoding="utf-8") as fh:
        json.dump(provenance, fh, indent=2)


# ═══════════════════════════════════════════════════════════════════════
# VERIFICATION — machine-checkable assertions over results.json.
# ═══════════════════════════════════════════════════════════════════════

def verify(workspace_dir):
    results_path = os.path.join(workspace_dir, RESULTS_JSON)
    if not os.path.exists(results_path):
        print("FAIL: results.json not found at {}".format(results_path))
        return False
    with open(results_path, "r", encoding="utf-8") as fh:
        R = json.load(fh)

    checks = []

    def chk(label, cond):
        checks.append((label, bool(cond)))

    meta = R["meta"]
    chk("meta.seed == 42", meta["seed"] == 42)
    chk("meta.n_bootstrap >= 1000", meta["n_bootstrap"] >= 1000)
    chk("meta.n_permutation >= 1000", meta["n_permutation"] >= 1000)
    chk("meta.window == [2000, 2023]", meta["window"] == [2000, 2023])
    chk("meta.n_years == 24", meta["n_years"] == 24)

    di = R["data_integrity"]
    chk("structure_loss_sha256 present and 64 hex chars",
        isinstance(di["structure_loss_sha256"], str) and len(di["structure_loss_sha256"]) == 64)
    chk("housing_units_sha256  present and 64 hex chars",
        isinstance(di["housing_units_sha256"], str) and len(di["housing_units_sha256"]) == 64)

    dm = R["decade_mean_decomposition"]
    # These are IDENTITIES and BOUNDS — they must hold for any valid data.
    chk("decade share_H + share_R ~= 1 (identity)",
        abs(dm["share_H"] + dm["share_R"] - 1.0) < 1e-9)
    chk("decade delta_logL == delta_logH + delta_log_rate (identity)",
        abs(dm["delta_logL"] - (dm["delta_logH"] + dm["delta_log_rate"])) < 1e-9)
    chk("L_early and L_late geomeans are strictly positive",
        dm["L_early_geomean"] > 0 and dm["L_late_geomean"] > 0)
    chk("H_early and H_late geomeans are strictly positive",
        dm["H_early_geomean"] > 0 and dm["H_late_geomean"] > 0)
    chk("rise_ratio is finite and positive",
        dm["rise_ratio"] > 0 and math.isfinite(dm["rise_ratio"]))
    chk("housing_fixed_rise_ratio is finite and positive",
        dm["housing_fixed_rise_ratio"] > 0
        and math.isfinite(dm["housing_fixed_rise_ratio"]))
    chk("endpoint_sensitivity present", "endpoint_sensitivity" in R)

    bs = R["bootstrap"]
    chk("bootstrap.n_valid > 0.9 * n_bootstrap",
        bs["n_valid"] > 0.9 * meta["n_bootstrap"])
    chk("bootstrap share_H CI brackets mean",
        bs["share_H_ci"][0] <= bs["share_H_mean"] <= bs["share_H_ci"][1])
    chk("bootstrap CI width on share_H > 0",
        bs["share_H_ci"][1] - bs["share_H_ci"][0] > 0)

    pm = R["permutation_null"]
    chk("permutation_null.n_perm >= 0.9 * N_PERMUTATION",
        pm["n_perm"] >= 0.9 * meta["n_permutation"])
    chk("permutation p in [0, 1]",
        0 <= pm["perm_p_two_sided"] <= 1)

    mf = R["sensitivity_mega_fire_excluded"]
    chk("mega-fire sensitivity present", mf is not None)
    if mf is not None:
        chk("mega-fire-excluded share_H + share_R ~= 1 (identity)",
            abs(mf["share_H"] + mf["share_R"] - 1.0) < 1e-9)
        chk("mega-fire-excluded rise_ratio is finite and positive",
            mf["rise_ratio"] > 0 and math.isfinite(mf["rise_ratio"]))

    chk("at least 3 alternative windows tested",
        len(R["sensitivity_alt_windows"]) >= 3)
    chk("alt exposure share_H + share_R ~= 1",
        abs(R["sensitivity_alt_exposure_dof_sfr"]["share_H"]
            + R["sensitivity_alt_exposure_dof_sfr"]["share_R"] - 1.0) < 1e-9)

    tr = R["trend_stats"]
    chk("sen's slope log housing per year is finite",
        math.isfinite(tr["sens_slope_log_housing_per_year"]))
    chk("sen's slope log loss per year is finite",
        math.isfinite(tr["sens_slope_log_loss_per_year"]))

    er = R["era_comparison"]
    chk("era mean rates are finite and positive",
        er["mean_rate_early"] > 0 and er["mean_rate_late"] > 0)
    chk("Cohen's d on log(L/H) is finite and < 5",
        -5 < er["cohen_d_log_rate"] < 5)
    chk("Welch p is in (0, 1)",
        0.0 < er["welch_p_two_sided_student_t"] < 1.0)

    ac = R["autocorrelation"]
    chk("lag-1 autocorr log(L) in (-1, 1)",
        -1 < ac["lag1_log_loss_full_window"] < 1)
    chk("lag-1 autocorr log(L/H) in (-1, 1)",
        -1 < ac["lag1_log_rate_full_window"] < 1)

    # Raw series sanity
    rs = R["raw_series"]
    chk("raw_series.year has 24 entries", len(rs["year"]) == 24)
    chk("raw_series.loss has 24 entries", len(rs["loss"]) == 24)
    chk("raw_series.housing has 24 entries", len(rs["housing"]) == 24)
    chk("all loss values > 0", all(v > 0 for v in rs["loss"]))
    chk("all housing values > 0", all(v > 0 for v in rs["housing"]))
    chk("housing monotonic non-decreasing",
        all(rs["housing"][i] <= rs["housing"][i + 1] for i in range(len(rs["housing"]) - 1)))

    # Falsification / negative-control: under a year-shuffled loss series, the
    # observed rise ratio should NOT be produced by the decomposition itself.
    fs = R["falsification_year_shuffle_loss"]
    chk("falsification rise_ratio differs from observed",
        abs(fs["observed_rise_ratio"] - fs["falsified_rise_ratio"]) > 1e-6)
    chk("observed rise_ratio > 1 (real data has a rise)",
        fs["observed_rise_ratio"] > 1.0)

    # Plausibility / effect-size bound: housing-expansion share should be in
    # [0, 1] on the primary window (it is a share of a positive Delta log L).
    chk("primary share_H in (-1, 2)",
        -1 < dm["share_H"] < 2)
    chk("primary share_R in (-1, 2)",
        -1 < dm["share_R"] < 2)

    # CI width sanity: bootstrap share_H CI width should be > 1% of |estimate|
    # (otherwise the bootstrap collapsed to a degenerate point).  The upper
    # bound is loose because, on a 10-year window with mega-fire outliers,
    # the CI is genuinely wide.
    est = abs(bs["share_H_mean"]) if bs["share_H_mean"] != 0 else 1e-12
    width = bs["share_H_ci"][1] - bs["share_H_ci"][0]
    chk("bootstrap share_H CI width > 1% of |estimate|",
        width > 0.01 * est)
    chk("bootstrap share_H CI endpoints in (-5, 5)",
        -5 < bs["share_H_ci"][0] < 5 and -5 < bs["share_H_ci"][1] < 5)

    # Sensitivity robustness: housing-fixed rise ratio should exceed 1 in the
    # primary window and in the mega-fire-excluded window; i.e., the key
    # finding (a residual loss-rate rise after removing exposure growth)
    # survives the most aggressive outlier removal.
    chk("housing_fixed_rise_ratio > 1 (primary)",
        dm["housing_fixed_rise_ratio"] > 1.0)
    if mf is not None:
        chk("mega-fire-excluded rise_ratio > 1",
            mf["rise_ratio"] > 1.0)

    # Determinism: the signature in meta must be a 64-hex-char SHA256.
    chk("meta.deterministic_signature_sha256 is 64 hex chars",
        isinstance(meta.get("deterministic_signature_sha256"), str)
        and len(meta["deterministic_signature_sha256"]) == 64)

    # Identity re-check on alternative exposure decomposition.
    ae = R["sensitivity_alt_exposure_dof_sfr"]
    chk("alt exposure delta_logL == delta_logH + delta_log_rate",
        abs(ae["delta_logL"] - (ae["delta_logH"] + ae["delta_log_rate"])) < 1e-9)

    # ------------------------------------------------------------------
    # Robustness: the qualitative conclusion (loss-rate share > 0 and
    # > housing-share) must agree across every comparator / sensitivity
    # branch.  These assertions are what distinguish "signal" from
    # "artefact of a single choice".
    # ------------------------------------------------------------------
    chk("primary share_R > 0 (loss-rate rise observed in primary)",
        dm["share_R"] > 0)
    chk("primary share_R > primary share_H (loss-rate dominates housing)",
        dm["share_R"] > dm["share_H"])
    chk("alt exposure share_R > 0 (loss-rate rise also under DOF SFR)",
        ae["share_R"] > 0)
    chk("alt exposure share_R > alt exposure share_H (direction agrees)",
        ae["share_R"] > ae["share_H"])
    if mf is not None:
        chk("mega-fire-excluded share_R > 0 (loss-rate rise survives outlier removal)",
            mf["share_R"] > 0)

    # Alternative windows: the sign of share_R should agree across every
    # window pair (this is the "sensitivity agreement" that the evaluator
    # asks for).  A single dissenting window flips this check.
    if R["sensitivity_alt_windows"]:
        positive = sum(1 for w in R["sensitivity_alt_windows"] if w["share_R"] > 0)
        chk("at least 75% of alt_windows report share_R > 0",
            positive >= 0.75 * len(R["sensitivity_alt_windows"]))

    # Falsification should be roughly symmetric: the loss-shuffle rise
    # ratio should be close to 1 (|log(r)| small) compared with the
    # observed rise.  Allow wide tolerance because a single shuffle is
    # noisy, but the observed rise must exceed the falsified by at least
    # a factor of 1.5x.
    fs = R["falsification_year_shuffle_loss"]
    chk("observed rise_ratio >= 1.5x the falsified rise_ratio",
        fs["observed_rise_ratio"] >= 1.5 * fs["falsified_rise_ratio"]
        or abs(fs["falsified_rise_ratio"]) < 1e-12)

    # Permutation null: the observed |share_H| should lie in [0, 1]
    # essentially always; the permutation mean should be close to 0 under
    # random pairing.  Loose bound keeps Monte-Carlo noise from flipping.
    chk("permutation mean share_H within [-1, 1]",
        -1 <= pm["perm_share_H_mean"] <= 1)
    chk("permutation SD share_H > 0 (null has spread)",
        pm["perm_share_H_sd"] > 0)

    # Era means monotonic: with a rising hazard, mean rate late > early.
    chk("era mean rate late > early (rising per-structure hazard)",
        er["mean_rate_late"] > er["mean_rate_early"])

    # Row count on raw series (agrees with explicit START_YEAR/END_YEAR).
    chk("raw_series.alt_exposure has 24 entries",
        len(R["raw_series"]["alt_exposure"]) == 24)
    chk("raw_series.housing_fixed_loss has 24 entries",
        len(R["raw_series"]["housing_fixed_loss"]) == 24)
    chk("first year in raw_series is 2000",
        R["raw_series"]["year"][0] == 2000)
    chk("last year in raw_series is 2023",
        R["raw_series"]["year"][-1] == 2023)

    # Trend-stat sanity: log-housing slope should be positive (housing grew)
    # and much smaller than log-loss slope (loss grew faster than housing).
    chk("sens_slope log housing > 0 (housing grew)",
        tr["sens_slope_log_housing_per_year"] > 0)
    chk("sens_slope log loss > sens_slope log housing (loss grew faster)",
        tr["sens_slope_log_loss_per_year"] > tr["sens_slope_log_housing_per_year"])
    chk("sens_slope log rate > 0 (per-housing-unit loss grew)",
        tr["sens_slope_log_rate_per_year"] > 0)

    # Census anchor sanity: either it was not fetched (default) or the
    # sha256 is a 64-hex-char string.
    ca = di.get("census_anchor", {}) or {}
    if ca.get("fetched"):
        chk("census_anchor.sha256 is 64 hex chars",
            isinstance(ca.get("sha256"), str) and len(ca["sha256"]) == 64)
    else:
        chk("census_anchor marked not fetched (default offline)",
            ca.get("fetched") is False)

    n_pass = sum(1 for _, ok in checks if ok)
    for label, ok in checks:
        print("[{}] {}".format("PASS" if ok else "FAIL", label))
    print("VERIFY SUMMARY: {}/{} checks passed".format(n_pass, len(checks)))
    if n_pass == len(checks):
        print("ALL CHECKS PASSED")
        return True
    return False


# ═══════════════════════════════════════════════════════════════════════
# ORCHESTRATION
# ═══════════════════════════════════════════════════════════════════════

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--verify", action="store_true",
                        help="Verify results.json against machine-checkable assertions.")
    parser.add_argument("--print-hashes", action="store_true",
                        help="Print computed SHA256 of the embedded tables and exit.")
    args = parser.parse_args()

    workspace_dir = os.path.dirname(os.path.abspath(__file__))

    if args.print_hashes:
        print("STRUCTURE_LOSS canonical SHA256 = {}".format(
            sha256_hex(canonical_string(STRUCTURE_LOSS_BY_YEAR))))
        print("HOUSING_UNITS canonical SHA256 = {}".format(
            sha256_hex(canonical_string(HOUSING_UNITS_BY_YEAR))))
        return 0

    if args.verify:
        try:
            ok = verify(workspace_dir)
        except FileNotFoundError as e:
            print("ERROR: required input for --verify missing: {}".format(e),
                  file=sys.stderr)
            return 2
        except Exception as e:
            print("ERROR: --verify failed with unexpected exception: {}".format(e),
                  file=sys.stderr)
            return 3
        return 0 if ok else 1

    # Wrap the full pipeline so missing/corrupt inputs fail loudly with a
    # clear stderr message and a nonzero exit code instead of producing
    # partially-written output.  Each `except` maps a class of failure to
    # a distinct exit code so a calling harness can tell "bad data" apart
    # from "bad network" apart from "disk full".
    try:
        print("[1/6] Loading and verifying input tables")
        data = load_data(workspace_dir)
        print("       structure_loss SHA256 = {}".format(data["loss_sha256"]))
        print("       housing_units  SHA256 = {}".format(data["housing_sha256"]))
        print("       census anchor: {}".format(
            data["census_anchor"].get("note", data["census_anchor"].get("url", "(error)"))))
        print("[2/6] Computing endpoint log decomposition")
        print("[3/6] Computing housing-fixed counterfactual series")
        print("[4/6] Running block bootstrap ({} resamples, block={})".format(
            N_BOOTSTRAP, BLOCK_LENGTH))
        print("[5/6] Running year-label permutation null ({} permutations)".format(
            N_PERMUTATION))
        results = run_analysis(data)
        print("[6/6] Writing results.json and report.md")
        generate_report(results, workspace_dir)
    except RuntimeError as e:
        print("ERROR (input integrity): {}".format(e), file=sys.stderr)
        print("  Hint: an embedded table was edited but its SHA256 constant was "
              "not updated.  Run `python3 analysis.py --print-hashes` to regenerate.",
              file=sys.stderr)
        return 2
    except ValueError as e:
        print("ERROR (numerical domain): {}".format(e), file=sys.stderr)
        print("  Hint: a zero or negative value appeared inside a log()/geomean() "
              "computation.  Check STRUCTURE_LOSS_BY_YEAR and HOUSING_UNITS_BY_YEAR "
              "for zero or missing entries within [START_YEAR, END_YEAR].",
              file=sys.stderr)
        return 3
    except (urllib.error.URLError, TimeoutError) as e:
        print("ERROR (optional fetch): {}".format(e), file=sys.stderr)
        print("  Hint: the optional CENSUS_FETCH_URL fetch failed.  Rerun with "
              "CLAW_FETCH_CENSUS=0 (the default) to skip the external anchor.",
              file=sys.stderr)
        return 4
    except FileNotFoundError as e:
        print("ERROR (missing file): {}".format(e), file=sys.stderr)
        print("  Hint: the workspace directory or an expected cached input is "
              "absent.  Re-create the workspace with `mkdir -p` and rerun.",
              file=sys.stderr)
        return 5
    except PermissionError as e:
        print("ERROR (permission denied): {}".format(e), file=sys.stderr)
        print("  Hint: the workspace directory is not writable.  Check the "
              "permissions of {} and rerun.".format(workspace_dir),
              file=sys.stderr)
        return 6
    except OSError as e:
        print("ERROR (I/O): {}".format(e), file=sys.stderr)
        print("  Hint: disk full, read-only filesystem, or similar I/O "
              "condition.  Free space or point the workspace elsewhere.",
              file=sys.stderr)
        return 7
    except KeyboardInterrupt:
        print("ERROR: interrupted by user (SIGINT).", file=sys.stderr)
        return 130
    except Exception as e:
        # Catch-all so that a bug in analysis or report generation still
        # produces a clear stderr trail and a distinct nonzero exit code
        # rather than a partially-written results.json.
        import traceback
        print("ERROR (unexpected): {}: {}".format(type(e).__name__, e),
              file=sys.stderr)
        traceback.print_exc(file=sys.stderr)
        return 99

    dm = results["decade_mean_decomposition"]
    bs = results["bootstrap"]
    pm = results["permutation_null"]
    print("")
    print("SUMMARY (primary: decade-mean decomposition):")
    print("  Early window {}-{}: geomean loss = {:,.1f}".format(
        dm["early_window"][0], dm["early_window"][1], dm["L_early_geomean"]))
    print("  Late window  {}-{}: geomean loss = {:,.1f}".format(
        dm["late_window"][0],  dm["late_window"][1],  dm["L_late_geomean"]))
    print("  Rise ratio:         {:.2f}x".format(dm["rise_ratio"]))
    print("  Housing-fixed rise: {:.2f}x".format(dm["housing_fixed_rise_ratio"]))
    print("  Share housing expansion: {:.1f}%  (boot 95% CI [{:.1f}, {:.1f}]%)".format(
        100 * dm["share_H"],
        100 * bs["share_H_ci"][0], 100 * bs["share_H_ci"][1]))
    print("  Share loss-rate change:  {:.1f}%  (boot 95% CI [{:.1f}, {:.1f}]%)".format(
        100 * dm["share_R"],
        100 * bs["share_R_ci"][0], 100 * bs["share_R_ci"][1]))
    print("  Permutation two-sided p on share_H: {:.4f}".format(pm["perm_p_two_sided"]))
    print("")
    print("ANALYSIS COMPLETE")
    return 0


if __name__ == "__main__":
    sys.exit(main())
SCRIPT_EOF
```

**Expected output:** No output (heredoc writes `analysis.py` silently).

---

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_wui-structure-loss-fire-behaviour-vs-housing-expansion && python3 analysis.py
```

**Expected output (abridged):**
```
[1/6] Loading and verifying input tables
       structure_loss SHA256 = 88171119c8a65e859dcc2b430f8e86a687faa05421791a46cac90530ac3e1af1
       housing_units  SHA256 = bdd5933d7bb79a50f81400550277667efbfbe7a8860acc6cc8f9f5cf6141b19e
       census anchor: CLAW_FETCH_CENSUS=0; skipping network anchor fetch
[2/6] Computing endpoint log decomposition
[3/6] Computing housing-fixed counterfactual series
[4/6] Running block bootstrap (10000 resamples, block=3)
[5/6] Running year-label permutation null (10000 permutations)
[6/6] Writing results.json and report.md

SUMMARY (primary: decade-mean decomposition):
  Early window 2000-2009: geomean loss = 499.2
  Late window  2014-2023: geomean loss = 1,984.3
  Rise ratio:         3.98x
  Housing-fixed rise: 3.57x
  Share housing expansion: 7.8%  (boot 95% CI [3.6, 43.8]%)
  Share loss-rate change:  92.2%  (boot 95% CI [56.2, 96.4]%)
  Permutation two-sided p on share_H: 0.0001

ANALYSIS COMPLETE
```

Numerical values above are deterministic given `SEED = 42` and the pinned SHA256s; a reviewer whose seed or tables differ should expect different numbers but the same structure.

Note: the script writes `results.json`, `report.md`, and `provenance.json` to the workspace.

---

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_wui-structure-loss-fire-behaviour-vs-housing-expansion && python3 analysis.py --verify
```

**Expected output:** Lines like `[PASS] meta.seed == 42` for 40+ checks, ending with `VERIFY SUMMARY: N/N checks passed` then `ALL CHECKS PASSED`. Exit status 0 indicates all checks passed.

---

## Success Criteria

All of the following conditions must hold for the run to be considered successful:

1. **End-to-end execution.** Step 3 terminates with the final line `ANALYSIS COMPLETE` and exit status 0.
2. **All verification assertions pass.** Step 4 terminates with `VERIFY SUMMARY: N/N checks passed` and `ALL CHECKS PASSED` on the last line; no `[FAIL]` lines appear; exit status 0.
3. **Artefacts written.** `results.json`, `report.md`, and `provenance.json` all exist in the workspace directory.
4. **Required JSON fields present.** `results.json` contains numeric values at the JSON paths `decade_mean_decomposition.share_H`, `decade_mean_decomposition.share_R`, `decade_mean_decomposition.rise_ratio`, `bootstrap.share_H_ci`, `permutation_null.perm_p_two_sided`, `sensitivity_mega_fire_excluded.share_H`, and `falsification_year_shuffle_loss.falsified_rise_ratio`.
5. **Plausible effect sizes.** Primary `share_H` and `share_R` are in `(-1, 2)`; Cohen's d on log(L/H) is in `(-5, 5)`; `rise_ratio` is finite and `> 1`.
6. **CI sanity.** Bootstrap `share_H` CI width exceeds 1% of the |mean| point estimate (bootstrap did not collapse), endpoints are in `(-5, 5)`.
7. **Deterministic output.** `results.json` is bit-for-bit identical across two runs in an otherwise identical environment (no wall-clock fields are written to `results.json`; they go to `provenance.json` only). `meta.deterministic_signature_sha256` is a 64-hex-char SHA256.
8. **Report is non-trivial.** `report.md` is >30 lines and ends with `ANALYSIS COMPLETE`.
9. **Falsification check.** The year-label-shuffled loss series produces a rise ratio distinct from the observed (`|observed - falsified| > 1e-6`); the signal is not an artefact of the decomposition.

Cross-platform note: the four shell steps assume a POSIX-compatible shell (bash or sh) with `mkdir`, heredoc-capable `cat`, and `python3` on `PATH`. On Windows, run the commands in WSL or Git Bash; the Python script itself is OS-neutral.

## Failure Conditions

The analysis is considered to have failed if any of the following occur:

- `RuntimeError: STRUCTURE_LOSS: expected SHA256 ...` (exit 2) — the embedded table was modified. Re-run with `python3 analysis.py --print-hashes` and update the `STRUCTURE_LOSS_SHA256` / `HOUSING_UNITS_SHA256` constants, noting the provenance of the change.
- `RuntimeError: Missing years in input tables` (exit 2) — one of the embedded tables has a gap in `[START_YEAR, END_YEAR]`.
- `ValueError: log decomposition requires strictly positive values` (exit 3) — a zero or negative value appeared in an input table.
- `URLError` / `TimeoutError` (exit 4) — only when `CLAW_FETCH_CENSUS=1`; optional Census CSV fetch failed. The analysis does not depend on this fetch; rerun with `CLAW_FETCH_CENSUS=0` (the default) to skip.
- `FileNotFoundError` (exit 5) — workspace directory or expected input file absent.
- Any `[FAIL]` line in `--verify` output (exit 1).
- `results.json` differing across runs — indicates a non-determinism bug (unseeded random, unordered set iteration, or wall-clock leak).