# Operational S2S extreme-heat-day forecast for Vietnam — results report

*Paper-ready synthesis. Every number is traceable to a run and a file (see the Evidence pack
at the end). Regenerate with the stage scripts listed in `CLAUDE.md`.*

## 1. Objective and forecast product

We produce a calibrated **subseasonal-to-seasonal (S2S) probability of extreme heat** over
Vietnam out to 45 days. The event is the **extreme-heat day (EHD)**, defined per grid cell and
calendar day as

> a day is an EHD ⇔ `Tx ≥ min(Tx90p, 35 °C)`,

where `Tx90p` is the calendar-day 90th percentile of observed daily maximum temperature over
the 2002–2023 baseline, estimated from a centred 15-day calendar-day window. The `min(·)` rule
is fixed by the user; the `max(Tx90p, 35 °C)` variant is retained only as a diagnostic and is
never substituted into the product.

The forecast target is aggregated into **six fixed, non-overlapping lead-time blocks** (days
since initialisation): **W1** 0–10, **W2** 11–17, **W3** 18–24, **W4** 25–31, **W5** 32–38,
**W6** 39–45. The **primary** product is `P(≥1 EHD in the block)`; the **secondary** product is
`P(heat-wave = ≥3 consecutive EHDs)` with block-boundary straddling handled on the continuous
0–45 day series and attributed to every block a run overlaps. The secondary target is **verified
under the same LOYO protocol** (§4.2); it is **not yet wired into the operational NetCDF** producer
(a documented next step, §6) — the delivered file below carries the primary product only.

The delivered field is a **CF-compliant NetCDF per initialisation** carrying the primary
probability `prob_ehd(block, lat, lon)` and a binary operational flag `ehd_flag(block, lat, lon)`
(probability ≥ the region's decision threshold), produced by `scripts/10_operational_product.py`.

## 2. Data and its consequences for method

- **Control-only reforecast (1 member).** The S2S archive contains only control (`_cf_`) files;
  there is no ensemble. Probabilities are therefore learned by **single-member MOS**
  (distributional/classification calibration of the deterministic member), **not** by an
  ensemble exceedance fraction. Ensemble diagnostics (rank histogram, spread–error) are
  undefined and reported as N/A.
- **Forecast already on the 0.25° observation grid** (41×65). No downscaling stage; the upstream
  interpolation is a documented limitation.
- **Forecast Tx is 6-hourly `mx2t6`** → daily Tx = daily max of the four 6-hourly values; day 0 =
  first 24 h after init, giving lead days 0–45.
- **Calendar mismatch:** gridded obs is `365_day` (29 Feb removed); forecast is Gregorian. We
  match on `(month, day)` and drop 29 Feb from forecast-derived Tx.
- **Target vs validation:** the **gridded** Tx is the calibration/verification target (spatially
  complete, aligned to the forecast grid); the **200 station** records are independent point
  validation. Stations run ≈ +1.5 to +2.2 °C warmer than their enclosing 0.25° cell, so
  station-scale EHD frequency is 2–5× the gridded frequency — a bound on point-scale risk, not a
  target (see `docs/OBS_RECONCILIATION.md`).

## 3. Method

**Features (per init × land cell × block), forecast-only:** block max / mean / min of forecast
daily Tx; fraction of block days with forecast `Tx ≥ threshold`; how far the hottest forecast day
clears the block-mean threshold; the block ordinal (lead dependence); the init-date seasonal
harmonics `sin/cos(2π·doy/365)`; and grid latitude/longitude. (`src/ehd/features.py`.)

**Calibration families:** two references — `raw` (forecast block-exceedance fraction, the
uncalibrated baseline) and `climat` (training base rate) — and two trainable single-member MOS
classifiers, standardised **logistic regression** and **histogram gradient boosting (HGB)**. Each
trainable family is additionally **isotonic-recalibrated**. The production/headline family is
**`hgb_recal`**.

**Decision threshold:** a per-**region × block** probability threshold that maximises Peirce's
skill score, fit on recalibrated training predictions and applied to make the binary flag.

**Cross-validation (leakage rule).** Skill is estimated by **leave-one-year-out (LOYO)** over the
20 reforecast years 2003–2022. For a held-out year *y*, none of the following sees *y*: (1) the
observed climatology and Tx90p threshold, (2) the calibration model, (3) the isotonic
recalibration, and (4) the decision threshold. Each stage records its fit-years; the validator
(§7) re-checks this. The production model (§6) is a separate **all-years refit** for deployment —
its honest skill is the LOYO estimate below, never the in-sample fit.

## 4. Headline cross-validated skill (primary target, `hgb_recal`)

Pooled over all 20 held-out folds (~1.12 M out-of-sample forecasts per block), scored against a
per-fold base-rate climatology, with a 95 % **whole-year block bootstrap** (resamples entire
years, preserving serial correlation):

| block | days | base rate | BSS vs clim [95 % CI] | AUC | BSS vs raw | significant |
|-------|------|-----------|-----------------------|-----|-----------|-------------|
| W1 | 0–10  | 0.431 | **0.356** [0.332, 0.383] | 0.845 | 0.517 | yes |
| W2 | 11–17 | 0.331 | **0.118** [0.092, 0.144] | 0.710 | 0.310 | yes |
| W3 | 18–24 | 0.331 | **0.037** [0.009, 0.063] | 0.630 | 0.288 | yes (marginal) |
| W4 | 25–31 | 0.331 | 0.014 [−0.012, 0.039] | 0.603 | 0.281 | no |
| W5 | 32–38 | 0.331 | 0.004 [−0.021, 0.027] | 0.591 | 0.280 | no |
| W6 | 39–45 | 0.331 | −0.006 [−0.028, 0.018] | 0.579 | 0.278 | no |

**Headline statement.** *The calibrated product has significant Brier skill over climatology out
to W3 (day 24); at and beyond W4 (day 25+) it is statistically indistinguishable from
climatology.* AUC stays > 0.5 at all leads (weak discrimination persists) but calibrated
probability skill (BSS) reaches zero by W4. **Calibration always beats the raw forecast**
(BSS_raw +0.28…+0.52) because the raw block-exceedance probability is badly miscalibrated (flat
reliability; `results/figures/cv_reliability.png`). This honest short-horizon result is preferred
over a decorated one.

Source: `results/tables/cv_scores_by_block.csv`, `cv_bootstrap_ci.csv`;
figure `results/figures/cv_skill_by_lead.png`.

### 4.1 Binary contingency inventory (headline family, at the region decision threshold)

The full ten-score binary inventory the scoring spec requires, pooled over folds for `hgb_recal`
at the operational per-region × block decision threshold (the flag actually issued):

| block | base | POD | POFD | FAR | CSI | ETS | HSS | Peirce | EDI | SEDI | freq bias |
|-------|------|-----|------|-----|-----|-----|-----|--------|-----|------|-----------|
| W1 | 0.431 | 0.761 | 0.237 | 0.291 | 0.580 | 0.351 | 0.520 | 0.524 | 0.682 | 0.682 | 1.07 |
| W2 | 0.331 | 0.644 | 0.345 | 0.520 | 0.379 | 0.160 | 0.275 | 0.299 | 0.415 | 0.417 | 1.34 |
| W3 | 0.331 | 0.583 | 0.402 | 0.582 | 0.322 | 0.090 | 0.165 | 0.181 | 0.257 | 0.259 | 1.39 |
| W4 | 0.331 | 0.557 | 0.414 | 0.600 | 0.303 | 0.070 | 0.130 | 0.143 | 0.203 | 0.205 | 1.39 |
| W5 | 0.331 | 0.549 | 0.428 | 0.612 | 0.294 | 0.058 | 0.109 | 0.121 | 0.171 | 0.173 | 1.41 |
| W6 | 0.331 | 0.521 | 0.415 | 0.617 | 0.283 | 0.051 | 0.097 | 0.106 | 0.148 | 0.152 | 1.36 |

Discrimination (Peirce/EDI/SEDI) and the event-focused CSI/ETS/HSS all decay monotonically with
lead and the frequency bias grows past 1 (the long-lead over-forecasting noted in §6), matching the
BSS collapse above. The per-**region** × block version of the same ten scores is in
`cv_scores_by_region_block.csv`. Source: `results/tables/cv_binary_scores_by_block.csv`,
`cv_scores_by_region_block.csv`; producer `08_verify.py`.

### 4.2 Secondary target — heat-wave (≥3 consecutive EHDs)

The secondary product is verified with the **same leakage-safe LOYO protocol** as the primary
(same `hgb_recal` family, per-fold Tx90p climatology, isotonic recalibration, and per-region × block
decision threshold; runner `scripts/06b_run_loyo_secondary.py`, scored by `08b_verify_secondary.py`).
The heat-wave event is rarer than the primary event (block base rate ≈ 0.11–0.16 vs 0.33–0.43).

| block | base rate | BSS vs clim [95 % CI] | AUC | SEDI | freq bias | significant |
|-------|-----------|-----------------------|-----|------|-----------|-------------|
| W1 | 0.158 | **0.269** [0.240, 0.300] | 0.849 | 0.691 | 2.00 | yes |
| W2 | 0.136 | **0.092** [0.058, 0.128] | 0.727 | 0.453 | 2.63 | yes |
| W3 | 0.136 | **0.033** [0.002, 0.059] | 0.649 | 0.302 | 2.59 | yes (marginal) |
| W4 | 0.136 | 0.011 [−0.023, 0.041] | 0.614 | 0.223 | 2.57 | no |
| W5 | 0.135 | 0.010 [−0.021, 0.037] | 0.612 | 0.221 | 2.66 | no |
| W6 | 0.112 | −0.005 [−0.033, 0.022] | 0.594 | 0.188 | 2.82 | no |

**Interpretation.** The heat-wave product has the **same skill horizon as the primary** — significant
Brier skill over climatology out to W3 (marginal), indistinguishable from climatology at W4+ — so the
"trust W1–W3, treat W4–W6 as climatology" rule carries over. Two honest caveats: (i) discrimination is
weaker than the primary at matched lead (rarer event, fewer positives to learn from), and (ii) the
binary heat-wave **flag over-forecasts strongly** (frequency bias 2.0–2.8): the Peirce-optimal
threshold on a rare event buys hit rate at the cost of many false alarms, so the flag should be read
as a *screening* alert alongside the probability, not a calibrated yes/no. Source:
`results/tables/cv_scores_secondary_by_block.csv`, `cv_bootstrap_secondary_ci.csv`,
`cv_binary_secondary_by_block.csv`; producers `06b_run_loyo_secondary.py`, `08b_verify_secondary.py`.

## 5. Stratified skill

### 5.1 Spatial / regional gradient (per-cell FDR field significance)
Per-cell Brier skill vs climatology with Benjamini–Hochberg FDR (q = 0.05) over a whole-year
bootstrap: significant cells = **539 / 504 / 192 / 116 / 68 / 30** of 541 for W1…W6. Regionally
(BSS by region × block), the tropical south retains skill much longer than the north:

- **R6 Tay Nguyen / R7 Nam Bo (south):** BSS ≈ 0.18–0.19 at W2, ≈ 0.10–0.12 at W3, and still
  positive at extended range — ≈ +0.07–0.08 at W4, decaying to ≈ +0.04 by W6.
- **R1–R3 (north):** skill gone by W3 and slightly negative W3–W6.

Interpretation: tropical heat is more persistent and less synoptically driven than the north.
**Operationally the skilful horizon is lead- *and* region-dependent: north ≈ W2, south ≈ W3+.**
Source: `outputs/ehd_skill_grid.nc`, `results/tables/cv_scores_by_region_block.csv`,
`regional_partition.csv`; figures `region_map.png`, `cv_skill_map_W1.png`, `regional_partition.png`.

### 5.2 Init-season stratification (`hgb_recal`, BSS vs climatology)

| block | DJF | MAM | JJA | SON |
|-------|-----|-----|-----|-----|
| W1 | 0.420 | 0.407 | 0.263 | 0.343 |
| W2 | 0.117 | 0.177 | 0.076 | 0.104 |
| W3 | 0.012 | 0.075 | 0.029 | 0.031 |
| W4 | −0.010 | 0.051 | −0.004 | 0.019 |
| W5 | −0.014 | 0.039 | −0.005 | −0.003 |
| W6 | −0.021 | 0.025 | −0.011 | −0.016 |

**MAM (the pre-monsoon peak-heat season) is the most predictable at long lead** — the only season
with positive BSS through W6 — whereas DJF and JJA lose skill by W3. This matters operationally:
the product is most trustworthy at extended range exactly in the season when extreme heat is most
consequential. Source: `results/tables/cv_scores_by_season_block.csv`; figure `cv_strata.png`.

### 5.3 Base-rate severity stratification (`hgb_recal`, BSS vs climatology)
Land cells split into Low/Moderate/High terciles of their climatological EHD frequency within
each block:

| block | Low | Moderate | High |
|-------|-----|----------|------|
| W1 | 0.317 | 0.357 | 0.392 |
| W2 | 0.125 | 0.105 | 0.124 |
| W3 | 0.047 | 0.021 | 0.043 |
| W4 | 0.027 | −0.002 | 0.017 |
| W5 | 0.015 | −0.010 | 0.008 |
| W6 | 0.007 | −0.021 | −0.002 |

At W1 skill rises with base rate (frequent-EHD cells are best forecast); by W4–W6 all terciles
are near zero. Source: `results/tables/cv_scores_by_severity_block.csv`.

## 6. Operational product and real-time path

`scripts/10_operational_product.py train` fits the **production model on all 20 years**: the
full-period Tx90p threshold field, one HGB classifier + isotonic recalibration per block, and the
per-region × block Peirce-optimal decision thresholds. Artefacts: `outputs/production_model/`
(`ehd_production_model.joblib`, `tx90p_threshold.nc`, `manifest.json`).

`... produce --year YYYY` scores every init of a forecast file through the saved model and writes
one CF-NetCDF per init to `outputs/operational/ehd_product_<YYYYMMDD>.nc` with
`prob_ehd(block, lat, lon)` and `ehd_flag(block, lat, lon)`.

**Real-time note.** The model consumes single-member (control) features. A real-time run must feed
the same single control member built identically (daily Tx = max of four 6-hourly `mx2t6`, drop
29 Feb, match on month/day). If a real-time *ensemble* is issued, reduce it to the control/mean
first. Verify the operational model version matches the reforecast the calibration was trained on;
document any mismatch. The **binary flag over-forecasts at long lead** (frequency bias 1.2–1.6 at
W4–W6), so it must be read *with* the verification scores in §4–5, not alone.

## 7. Independent validation

`scripts/11_independent_validator.py` re-reads the raw LOYO predictions and recomputes the
headline per-block Brier, BSS vs climatology, AUC and the bootstrap CI **from scratch in plain
numpy, without importing `ehd.verify` or reusing stage-08 code**, then asserts agreement with the
published tables within tolerance and re-checks fold coverage/leakage. Output:
`results/tables/validator_headline_rederivation.csv` (PASS/FAIL printed to stdout).

## 8. Scoring inventory
Calibrated forecasts are scored against **both** references (climatology and the raw forecast).
Probabilistic: Brier / BSS (+ decomposition via reliability), reliability curves, ROC/AUC,
sharpness. Binary (at the region decision threshold): the full ten-score inventory — POD, POFD, FAR, CSI,
ETS, HSS, Peirce, EDI, SEDI, frequency bias — reported by block (§4.1) and by region × block
(`cv_scores_by_region_block.csv`), and for the v2 zoo/ensembles (§10) and the vote ensemble (§11).
Continuous MOS metrics (MAE/RMSE/bias/CRPS) are implemented in `src/ehd/verify.py` for the
Tx side. Significance: whole-year block bootstrap; field significance across cells via FDR.
Rank histogram / spread–error are N/A (single member).

## 9. Limitations
1. **No ensemble** — probabilities come from single-member MOS, not ensemble spread.
2. **Upstream interpolation** to 0.25° is baked into the archive and cannot be undone here.
3. **Gridded-vs-station warm offset** — point-scale EHD risk exceeds the gridded product.
4. **Skill horizon is short** — W1–W3 pooled, region- and season-dependent; W4–W6 climatology-
   equivalent on the pooled national average.
5. **Long-lead binary flag over-forecasts** (frequency bias > 1) — use with the scores.
6. **Region shapefile `7,NB_region.shp` is corrupt**; maps use the gridded `reg` mask + coastlines.
7. **Heat-wave (secondary) flag over-forecasts** (frequency bias 2.0–2.8, §4.2) and the product is
   verified but not yet in the operational NetCDF — wiring it into `10_operational_product.py` is
   the documented next step.

## 10. v2 enhancements — severity predictand, expanded zoo, cascade, ensemble

The v1 result above is the load-bearing deliverable and is unchanged. v2 (user request, 2026-09-14)
extends the *calibration and verification* only; it keeps every v1 hard constraint (the six fixed
blocks, `EHD = min(Tx90p, 35 °C)`, the leakage rule) and the same event (`yEHDE = ≥1 EHD in block`,
identical to the v1 primary target), so v1 and v2 are directly comparable.

**What changed.**
1. **Severity predictand.** Regressors learn block-max observed severity
   `sev = max_{d∈block}(Tx_obs,d − Thr_d)`, `Thr = min(Tx90p, 35)`; by construction `sev ≥ 0 ⇔ yEHDE`.
   Predicted `sev` is mapped to `P(yEHDE)` by an isotonic map fit on TRAIN only. Forecast severity is
   also added as a feature.
2. **Expanded model zoo (all single-member MOS, LOYO-nested):** classifiers `logistic, hgb, rf, mlp`
   (torch ANN); severity regressors `ridge_sev, hgb_sev, mlp_sev`; and a **fully-convolutional CNN**
   plus a **deep torch ANN** (`cnn`, `ann_torch`) added in Phase B — nine members total.
3. **Two-level cascade.** LV1 (`hgb`) screens samples it calls near-certain non-events
   (`P(yEHDE) ≤ 1 − cutoff`, cutoff swept 0.90→0.95→0.99) and passes their LV1 probability through;
   LV2 zoo models are trained on the surviving harder samples. The cutoff sweep moves scores only in
   the third decimal — **the easy-negative screen is ~skill-neutral** here.
4. **Whole-dataset (2003–2022 LOYO-pooled) binary verification** with the full contingency inventory
   (POD, POFD, FAR, CSI, ETS, HSS, Peirce, EDI, SEDI, frequency bias) per block × model × cutoff.

**Ensemble vs. members (headline of v2).** To compare fairly, **every** member — including the
CNN and deep ANN — is scored on the *identical* leakage-safe pipeline as the ensembles (per-held-out-
year, per-region Peirce decision threshold; `scripts/14b_single_member_scores.py`). On that common
footing the best single member is the **fully-convolutional CNN at every block**. The stacked
ensemble (`ens_stack`) is the best forecast at **short lead (W1–W2)**, adding a clear increment over
the CNN; but at **long lead (W3–W6) the CNN matches or slightly exceeds `ens_stack`** — it wins on
BSS and AUC at W3–W6 and on SEDI at W3 and W6. This is consistent with §11, where the CNN is the
strongest single member at extended range. So stacking helps most where the signal is strongest
(short lead); at long lead no combination beats the single CNN. (`ens_mean` is marginally best at W6.)

| block | best single member = CNN — SEDI / BSS / AUC | `ens_stack` — SEDI / BSS / AUC | winner |
|-------|---------------------------------------------|--------------------------------|--------|
| W1 | 0.707 / 0.385 / 0.859 | **0.715 / 0.391 / 0.864** | ens_stack |
| W2 | 0.445 / 0.136 / 0.722 | **0.454 / 0.139 / 0.725** | ens_stack |
| W3 | **0.292 / 0.063 / 0.647** | 0.291 / 0.059 / 0.642 | ~tie (CNN on BSS/AUC) |
| W4 | 0.245 / **0.046 / 0.622** | **0.249** / 0.041 / 0.616 | split (CNN BSS/AUC, ens SEDI) |
| W5 | 0.190 / **0.033 / 0.606** | **0.203** / 0.029 / 0.596 | split (CNN BSS/AUC, ens SEDI) |
| W6 | **0.181 / 0.027 / 0.593** | 0.170 / 0.022 / 0.581 | CNN |

(Bold = better value for that metric at that lead. CNN scores from `v2_single_member_scores_by_block.csv`;
`ens_stack` from `v2_scores_with_ensemble.csv` — both on the stage-14 pipeline.)

**Interpretation.** v2 confirms the v1 horizon story — strong skill at W1, meaningful through W3,
fading to near-climatology by W4–W6. Stacking adds a clear increment at short lead
(W1 SEDI 0.707→0.715, BSS 0.385→0.391) but does **not** beat the single CNN at long lead, where the
CNN is the most skilful individual model. Per-cell binary skill maps per lead week
(`outputs/v2_skill_maps.nc`) reproduce the north-vs-south gradient of §5.1, and a self-contained
`results/v2_report.html` bundles the tables, skill-by-lead curve, and maps. Source:
`results/tables/v2_scores_with_ensemble.csv`, `v2_single_member_scores_by_block.csv`,
`v2_best_single_member_by_block.csv`, `v2_best_model_by_block.csv`,
`v2_cv_scores_by_model_block.csv`; figures `results/figures/v2_ensemble_skill_by_lead.png`,
`v2_skillmap_*.png`.

## 11. Vote-fraction ensemble and regional event forecast (prompt §6)

A second, model-agnostic combination and an operational regional product, built entirely on the
existing leakage-safe LOYO predictions (`scripts/17_vote_ensemble_regional.py`).

**Members and vote.** The pool is **23 AI members** — 7 sklearn families
(logistic, hgb, rf, mlp, ridge_sev, hgb_sev, mlp_sev) × 3 cascade cutoffs (1.0, 0.90, 0.99) plus
the CNN and deep-ANN (torch) members. Each member casts a binary yEHDE vote using **its own
leakage-safe per-region×block decision threshold**, and the probability is the vote fraction
`P(yEHDE) = m/n` (n = 23). The deterministic forecast is `EHDE = 1 where P > P₀`.

**§6.1 Member comparison.** Ranked on the stated priority (high POD, low POFD/FAR, esp. W3–W6),
the **CNN is the strongest single member at long lead** (W3: POD 0.62, POFD 0.41, FAR 0.57,
SEDI 0.30 — best of the 23); the linear members (logistic, ridge_sev) are weakest. Leave-one-
member-out shows the vote is **highly redundant**: dropping any one member changes ensemble Peirce
by ≲0.002 (CNN the largest positive contributor). Source
`results/tables/v6_member_scores_by_block.csv`, `v6_member_contribution.csv`;
figure `v6_member_pod_pofd_far.png`.

**§6.2–6.3 Vote fraction and optimal P₀.** The Peirce-optimal P₀ sits near a **simple majority
vote** (P₀ ≈ 0.5 at W1/W5/W6, ≈ 0.46 at W2–W4). At the per-fold P₀ (leakage-safe):

| block | POD | POFD | FAR | Peirce | SEDI | P₀ |
|-------|-----|------|-----|--------|------|-----|
| W1 | 0.767 | 0.249 | 0.300 | 0.517 | 0.674 | 0.50 |
| W2 | 0.680 | 0.377 | 0.529 | 0.303 | 0.423 | 0.46 |
| W3 | 0.607 | 0.414 | 0.579 | 0.193 | 0.274 | 0.46 |
| W4 | 0.589 | 0.425 | 0.593 | 0.164 | 0.235 | 0.48 |
| W5 | 0.556 | 0.423 | 0.606 | 0.133 | 0.191 | 0.50 |
| W6 | 0.549 | 0.427 | 0.612 | 0.122 | 0.175 | 0.50 |

The unweighted vote roughly matches the **best single member** and is a little below the trained
`ens_stack` (§10) at short lead (W1 SEDI 0.674 vs 0.715) — expected, since the vote does no
fitting. It is the most transparent product. Source `results/tables/v6_vote_scores_by_block.csv`,
`v6_p0_sweep.csv`; figure `v6_p0_sweep.png`; per-cell maps `outputs/v6_vote_maps.nc`,
`results/figures/v6_vote_skillmap.png`.

**§6.4 Regional event forecast.** For region `R_j`, week `W_k` and each init, forecast coverage
`P = M_j/N_j` (fraction of the region's cells forecast yEHDE) is thresholded at `P_thres` to alert
the region; the same M/N rule on the observed field defines the regional truth. **Aggregating to
region×week is markedly more skilful than per-cell** (W1 Peirce 0.65 / SEDI 0.80, still positive at
W6), because cell-scale noise averages out. `P_thres` chosen per week by Peirce is unstable
(W5 → 0.10); the fixed-0.5 majority reference (`*_at0p5` columns) is steadier and recommended for
operations. Source `results/tables/v6_regional_by_week.csv`, `v6_regional_event_table.csv`;
figure `v6_regional_event_heatmap.png`.

## 12. Stage 2 — extended predictors, feature selection, and calibration recovery

Stage 1 (§1–11) learns `P(EHD)` from forecast-Tx statistics only. Stage 2 asks whether **additional
predictors** — land-surface, large-scale circulation, and teleconnections — extend the skilful
horizon, using a principled, leakage-safe **feature-selection-first** design.

**Candidate pool (28 predictors, four mechanisms).** Beyond the 10 Stage-1 forecast-backbone
features: **land-surface** (soil moisture `sm20/sm100`, deep soil temp `st100`) and **local
thermodynamics** (`skt`, `cape`, `tcw`, `tcc`, `d2m`) block-aggregated on the model grid
(`scripts/20a_build_ext_grid.py`); **large-scale circulation** indices (Z500 over a Vietnam/SCS box,
a subtropical box, their gradient, and box-mean T850) from the wider ERA domain; **teleconnections**
(ENSO Niño-3.4 for the init month; MJO OMI PC1/PC2/amplitude at init) and **persistence** (antecedent
15-day observed EHD-frequency and mean severity). All are leakage-safe: gridded/large-scale fields are
pure forecast (threshold-independent), teleconnections are observed init states, and persistence uses
the fold threshold on observations strictly before init.

**Screening (leakage-safe, nested).** For each lead block we rank candidates by **out-of-sample
permutation importance** — the drop in held-out-year AUC when a predictor is shuffled, aggregated over
the 20 LOYO folds (`scripts/20c_feature_screen.py`) — cross-checked with univariate mutual
information. Selection keeps the backbone plus any candidate above a noise floor on mean OOS
importance; stability (std) is reported, **not** used as a gate (episodic teleconnections are
high-variance yet dominant). The result is physically coherent and recovers the textbook S2S
transition **out-of-sample**:

- **Short lead (W1–W2):** forecast-Tx statistics and box-mean **T850** dominate.
- **Extended range (W3–W6):** the forecast's own Tx signal is exhausted and slow-varying predictors
  take over — **ENSO Niño-3.4 is the single most important predictor at W3–W6**, followed by MJO,
  Z500/T850, and antecedent persistence.
- **Honest null:** soil moisture/temperature never clear the noise floor (only CAPE at W1) — in this
  control reforecast, land-surface memory adds little once ENSO/circulation/persistence are present.

Per-block selected extras: W1 +6 (T850, persistence, Z500-grad, MJO, CAPE); W2–W4 +9 (the full ENSO
+ MJO + Z500 + T850 + persistence suite); W5 +7; W6 +4 (ENSO, Z500, persistence). Source
`results/tables/feature_importance_by_block.csv`, `feature_selection_by_block.csv`; figure
`results/figures/feature_importance_heatmap.png`.

**Naive addition miscalibrates.** A matched HGB retrain (backbone vs backbone+selected, nested LOYO,
`scripts/20d_retrain_selected.py`/`20e_verify_v3.py`) showed the extended predictors **raise
discrimination** at W3–W6 (AUC +0.02…+0.05, SEDI +0.05…+0.12) but **degrade calibrated BSS** — a
sharper-but-less-reliable model, because the isotonic map had been fit in-sample.

**Calibration recovery (the Stage-2 model).** Refitting calibration **out-of-fold** — an inner
K-fold over training years produces OOF probabilities on which we fit (i) isotonic recalibration and
(ii) a logistic **stack** of {backbone, selected} — then refitting the HGBs on the full training years
(`scripts/20f_recover_calibration.py`/`20g_verify_v3cal.py`) converts the discrimination gain into
**calibrated** skill. The stack beats the backbone at **every** lead (whole-year bootstrap on the gain
excludes zero everywhere):

| block | base BSS | sel BSS | **stack BSS** | stack−base gain [95 % CI] | AUC base→stack | SEDI base→stack |
|-------|----------|---------|---------------|---------------------------|----------------|-----------------|
| W1 | 0.348 | 0.357 | **0.377** | +0.029 [0.025, 0.034] | 0.842→0.856 | 0.674→0.700 |
| W2 | 0.115 | 0.123 | **0.142** | +0.027 [0.014, 0.042] | 0.707→0.730 | 0.411→0.460 |
| W3 | 0.039 | 0.061 | **0.072** | +0.033 [0.015, 0.049] | 0.622→0.658 | 0.242→0.307 |
| W4 | 0.022 | 0.045 | **0.051** | +0.029 [0.014, 0.043] | 0.582→0.639 | 0.173→0.258 |
| W5 | 0.015 | 0.050 | **0.053** | +0.038 [0.020, 0.054] | 0.570→0.639 | 0.143→0.274 |
| W6 | 0.011 | 0.031 | **0.034** | +0.023 [0.011, 0.036] | 0.535→0.608 | 0.101→0.215 |

**Horizon extension.** With an absolute whole-year bootstrap of BSS vs climatology, the Stage-2 stack
is **significant at every block through W6** (day 45): W4 0.051 [0.028, 0.069], W5 0.053 [0.030,
0.070], W6 0.034 [0.017, 0.048] — whereas the Stage-1 headline lost significance at W4. Two separable
effects contribute: the OOF calibration lifts even the backbone marginally (base becomes significant
at W4–W5), and the selected extended predictors add substantially on top (the significant stack−base
gains above). Together they push significant extreme-heat probability skill from **W3 to W6**, with the
long-lead skill physically sourced from ENSO/MJO/circulation/persistence. The Stage-2 flag also
over-forecasts less at long lead than Stage 1 (W4 frequency bias 1.49 vs base 1.78). Source
`results/tables/v3cal_scores_by_block.csv`, `v3cal_bss_gain_bootstrap.csv`,
`v3cal_absolute_bss_ci.txt`, `v3cal_binary_by_block.csv`; figure `results/figures/v3cal_skill_by_lead.png`.

**Caveats.** Stage-2 models are trained under a 200k-row cap for feasibility (matched across base and
selected, so gains are fair); the horizon extension combines a calibration-methodology gain with a
predictor gain, both reported separately; and, as throughout, the reforecast is control-only so no
ensemble-spread predictors exist. These do not change any Stage-1 hard constraint.

## 12b. Stage 3 — hybrid model zoo, spatial deep learning, and the diversity stack

Stage 2 (§12) established the **lead-dependent extended predictors** but exercised them through a
single calibrated HGB. Stage 3 follows the hybrid-design recommendation (error-diverse, calibrated
probability stack over a lead-aware zoo) and closes the gap the earlier zoo left open: it lets the
**full model zoo consume the lead-dependent predictor slice** and adds two model classes v2 lacked —
**gradient-boosting diversity** (XGBoost, LightGBM alongside sklearn HGB) and **true spatial deep
learning** (a residual U-Net over the Vietnam grid, alongside the plain CNN). Every member is trained
LOYO-nested with out-of-fold isotonic calibration and per-region×block Peirce thresholds — the leakage
rule is unchanged, and so are all v1 hard constraints (six fixed blocks, `EHD = min(Tx90p,35)`).

**The 13-member calibrated stack.** Members: `logistic, hgb, rf, mlp, xgb, lgbm` classifiers;
`ridge_sev, hgb_sev, mlp_sev, xgb_sev, lgbm_sev` severity regressors; and two spatial DL members
`cnn, resunet`. A LOYO-safe logistic stack (L2; an L1-sparse variant selects members automatically)
combines their out-of-sample probabilities per block. The stack has **positive Brier skill over
climatology at every lead, and equals or beats the best single member throughout:**

| block | base rate | best single member | **ens_stack BSS** | ens_mean BSS | ens_stack SEDI | ens_stack Peirce |
|-------|-----------|--------------------|-------------------|--------------|----------------|------------------|
| W1 | 0.431 | resunet 0.421 | **0.431** | 0.383 | 0.747 | 0.589 |
| W2 | 0.331 | xgb/cnn 0.140 | **0.157** | 0.140 | 0.487 | 0.353 |
| W3 | 0.331 | logistic/resunet 0.078 | **0.087** | 0.078 | 0.356 | 0.253 |
| W4 | 0.331 | logistic 0.066 | **0.071** | 0.062 | 0.320 | 0.226 |
| W5 | 0.331 | logistic 0.058 | **0.062** | 0.058 | 0.323 | 0.228 |
| W6 | 0.331 | logistic 0.053 | **0.058** | 0.053 | 0.296 | 0.208 |

Source `results/tables/v3_stack_scores_by_block.csv`, `v3_stack_diversity_by_block.csv`,
`v3zoo_scores_by_model_block.csv`, `v3zoo_best_model_by_block.csv`; figure
`results/figures/v3_stack_by_lead.png`.

**Two physically coherent findings — which member wins depends on lead.**
- **Short lead (W1–W2):** the new members win. Gradient boosting leads among the tabular families
  (W1 lgbm BSS 0.372, hgb 0.367, xgb 0.349 vs logistic 0.249), and the **spatial DL is best of all**
  — resunet BSS 0.421, cnn 0.409 at W1 — because near the initialisation there is a strong,
  spatially-structured forecast-Tx signal for a convolutional model to exploit.
- **Extended range (W3–W6):** plain **logistic regression on the lead-dependent large-scale
  predictors (ENSO/Z500/T850/persistence) is the single best family** (BSS W3 0.076, W4 0.066,
  W5 0.058, W6 0.053) — the boosting trees overfit the weak teleconnection signal, exactly as the
  hybrid design predicted. The stack recovers the logistic long-lead skill and adds a little on top.

**Diversity diagnostics confirm the stack is more than an average.** Per-member error-correlation and
leave-one-member-out ΔBSS (`v3_stack_diversity_by_block.csv`) show the spatial DL is the most
**error-diverse** member at short lead (W1 resunet mean |err-corr| 0.87 — the lowest in the pool — vs
~0.95 for the trees) and the only member with a clearly positive leave-one-out contribution
(W1 ΔBSS +0.020 for resunet, +0.004 for cnn). At long lead the tabular members are highly correlated
(|err-corr| ~0.98) and near-redundant; logistic carries the signal and the L1-sparse stack collapses
to it, so `ens_l1 ≡ ens_stack` at every block. This is why the stack's short-lead gain is large
(W1 0.383→0.431) and its long-lead gain modest but consistent.

**Horizon extension over Stage 2.** The Stage-3 stack extends the skilful horizon further than the
Stage-2 calibrated model at long lead: **W4 0.051→0.071, W5 0.053→0.062, W6 0.034→0.058** (BSS vs
climatology), while matching or exceeding it at short lead (W1 0.377→0.431). Significant probability
skill therefore holds through **W6 (day 45)**, now with a short-lead boost sourced from spatial DL and
a long-lead floor sourced from the large-scale predictors.

**Tested and rejected.** Two deeper DL families were built and evaluated but **excluded from
production**: a feature-tokenizer **Transformer** and a deep torch **ANN** are BSS-negative at long
lead (overfitting the weak extended-range signal, matching the design note's own Transformer caveat).
They remain in the codebase (`src/ehd/gridmodels.py`) for reproducibility but are not stacked.

**Caveats.** (i) Members run under a training-row cap for feasibility, matched across families so
comparisons are fair. (ii) **Honest DL limitation:** the residual U-Net operates on Vietnam-grid
channel images; the large-scale predictors enter only as broadcast box-index channels, **not** as true
spatial fields over the extended domain — feeding gridded Z500/T850/SST as spatial channels needs a
new full-field cache. **This was subsequently built and tested in §12c (Stage 4): the true spatial
field adds no skill**, so the limitation turns out not to cost skill. (iii) Control-only reforecast throughout,
so no ensemble-spread predictors exist. None of this changes a v1 hard constraint. New code:
`src/ehd/gridmodels.py` (`TorchResUNet`, `TorchFTTransformer`), extended `src/ehd/calibration.py`
(xgb/lgbm families), `config/config.yaml` `predictors:` block, `scripts/21_run_loyo_v3zoo.py`,
`scripts/22b_run_dl_v3.py`, `scripts/23_v3_diversity_stack.py`, `scripts/25_verify_v3zoo.py`.

## 12c. Stage 4 — true full-field spatial deep learning (the honest limitation, tested)

§12b closed with an honest caveat (ii): the Stage-3 residual U-Net saw the large-scale predictors
only as **broadcast box-index channels**, never as **true spatial fields** over the extended domain.
Stage 4 removes that limitation directly and tests whether it was costing skill. A **comprehensive
full-field cache** was built (`scripts/27_build_fullfield_cache.py`): block-mean of the raw gridded
S2S fields over the six fixed lead blocks, ocean retained, row-aligned to the forecast init order —
pure reductions, leakage-free. **55 channels**: extended-domain (71×61, 1°) Z500/T850/U200/V200 and
Vietnam-grid (41×65, 0.25°) surface + multi-level T/H/W/Q/U/V (`outputs/fullfield_cache/ff_*.npz`,
6.2 GB / 20 yr). A **two-branch fully-convolutional net** (`src/ehd/ffnet.py` `FullFieldResUNet`)
consumes it: a residual U-Net over the Vietnam grid (61 channels = 51 full-field + 10 rasterised
Stage-1 forecast features) plus a small CNN over the extended domain that pools to a context vector
tiled back onto the VN grid. Trained LOYO-nested with masked BCE, out-of-fold isotonic calibration
and per-region×block thresholds — leakage rule and all v1 hard constraints unchanged
(`scripts/28_run_loyo_ffnet.py`, 20-fold LOYO → `outputs/loyo_pred_v4ff/`).

**Verdict — a clean negative result: the true spatial field adds essentially no skill.** Scored
out-of-sample against the Stage-3 box-index resunet baseline, the deltas are ≈ 0 or slightly
**negative** at every lead:

| block | ffresunet BSS | box-index resunet BSS | Δ BSS | ffresunet SEDI | box-index SEDI | Δ SEDI |
|-------|---------------|-----------------------|-------|----------------|----------------|--------|
| W1 | 0.392 | 0.421 | **−0.028** | 0.723 | 0.740 | −0.017 |
| W2 | 0.139 | 0.139 | −0.000 | 0.456 | 0.477 | −0.021 |
| W3 | 0.079 | 0.078 | +0.001 | 0.356 | 0.342 | +0.013 |
| W4 | 0.050 | 0.052 | −0.001 | 0.290 | 0.303 | −0.013 |
| W5 | 0.042 | 0.046 | −0.005 | 0.263 | 0.307 | −0.043 |
| W6 | 0.044 | 0.052 | −0.007 | 0.264 | 0.290 | −0.026 |

Source `results/tables/v4ff_scores_by_model_block.csv`, `v4_vs_v3_resunet.csv`; figure
`results/figures/v4ff_skill_by_lead.png`.

**Interpretation.** For this predictand — P(≥1 EHD in block) over Vietnam — the discriminative signal
lives in the **local, grid-cell-scale forecast value**, not in the large-scale spatial *morphology* of
the extended domain. Once the box index (which cell/region a sample belongs to) and the local forecast
Tx are supplied, resolving the actual Z500/T850 patterns as spatial channels adds no separable
information the calibrated members were missing; the extra capacity even costs a little at W1 and the
long leads (mild overfitting of the weak extended-range signal, consistent with §12b). **The §12b
"honest limitation" was therefore not costing skill.** The full-field ResU-Net is **not merged into the
production diversity stack** — there is no skill lift to justify its cost — but the cache and net remain
in the codebase as the definitive test of the spatial-morphology hypothesis. New code:
`scripts/27_build_fullfield_cache.py`, `src/ehd/ffnet.py`, `scripts/28_run_loyo_ffnet.py`,
`config/config.yaml` `fullfield:` block.

## 13. Evidence pack — number → source
| Claim | File(s) | Producing script |
|-------|---------|------------------|
| Headline BSS/AUC per block | `results/tables/cv_scores_by_block.csv` | `08_verify.py` |
| Headline binary inventory per block (10 scores) | `results/tables/cv_binary_scores_by_block.csv` | `08_verify.py` |
| Secondary (heat-wave) LOYO scores + CI + binary inventory | `results/tables/cv_scores_secondary_by_block.csv`, `cv_bootstrap_secondary_ci.csv`, `cv_binary_secondary_by_block.csv` | `06b_run_loyo_secondary.py`, `08b_verify_secondary.py` |
| Secondary out-of-sample predictions (raw) | `outputs/loyo_pred_secondary/pred_hw_*.npz` | `06b_run_loyo_secondary.py` |
| Bootstrap CIs / significance | `results/tables/cv_bootstrap_ci.csv` | `08_verify.py` |
| Region × block scores | `results/tables/cv_scores_by_region_block.csv` | `08_verify.py` |
| Regional partition + names | `results/tables/regional_partition.csv`, `region_centroids.csv` | `09_regional_products.py`, `01_*` |
| FDR field significance / skill grid | `outputs/ehd_skill_grid.nc` | `09_regional_products.py` |
| Season stratification | `results/tables/cv_scores_by_season_block.csv` | `08_verify.py` |
| Severity stratification | `results/tables/cv_scores_by_severity_block.csv` | `08_verify.py` |
| Reliability / skill-by-lead / strata figures | `results/figures/cv_reliability.png`, `cv_skill_by_lead.png`, `cv_strata.png` | `08_verify.py` |
| Obs reconciliation | `results/tables/obs_reconciliation_*.csv`, `docs/OBS_RECONCILIATION.md` | `02_obs_reconciliation.py` |
| Production model + per-init product | `outputs/production_model/`, `outputs/operational/` | `10_operational_product.py` |
| Independent re-derivation | `results/tables/validator_headline_rederivation.csv` | `11_independent_validator.py` |
| Out-of-sample predictions (raw) | `outputs/loyo_pred/pred_*.npz` | `06_run_loyo.py` |
| v2 per-model binary scores | `results/tables/v2_cv_scores_by_model_block.csv`, `v2_best_model_by_block.csv` | `12_run_loyo_v2.py`, `13_verify_v2.py` |
| v2 ensemble scores | `results/tables/v2_scores_with_ensemble.csv` | `14_ensemble_v2.py` |
| v2 single-member scores on the ensemble pipeline (incl. CNN/ANN) | `results/tables/v2_single_member_scores_by_block.csv`, `v2_best_single_member_by_block.csv` | `14b_single_member_scores.py` |
| v2 per-cell skill maps | `outputs/v2_skill_maps.nc`, `results/figures/v2_skillmap_*.png` | `15_skill_maps_v2.py` |
| v2 CNN + deep ANN predictions | `outputs/loyo_pred_v2_cnn/pred_v2cnn_*.npz` | `12b_run_cnn_ann.py` |
| v2 self-contained HTML report | `results/v2_report.html` | `16_report_v2_html.py` |
| v2 out-of-sample predictions (raw) | `outputs/loyo_pred_v2/pred_v2_*.npz` | `12_run_loyo_v2.py` |
| §6 member comparison + LOMO contribution | `results/tables/v6_member_scores_by_block.csv`, `v6_member_contribution.csv` | `17_vote_ensemble_regional.py` |
| §6 vote ensemble + P0 sweep | `results/tables/v6_vote_scores_by_block.csv`, `v6_p0_sweep.csv` | `17_vote_ensemble_regional.py` |
| §6 regional event table | `results/tables/v6_regional_by_week.csv`, `v6_regional_event_table.csv` | `17_vote_ensemble_regional.py` |
| §6 vote-ensemble skill maps | `outputs/v6_vote_maps.nc`, `results/figures/v6_*.png` | `17_vote_ensemble_regional.py` |
| Stage-2 gridded candidate predictors (cache) | `outputs/ext_predictors/extgrid_*.npz` | `20a_build_ext_grid.py` |
| Stage-2 candidate design matrix (cache) | `outputs/candidates/cand_*.npz` | `20b_assemble_candidates.py` |
| Stage-2 feature importance + selection | `results/tables/feature_importance_by_block.csv`, `feature_selection_by_block.csv`; fig `feature_importance_heatmap.png` | `20c_feature_screen.py` |
| Stage-2 selected-retrain (naive) scores | `results/tables/v3_scores_by_block.csv`, `v3_bss_gain_bootstrap.csv` | `20d_retrain_selected.py`, `20e_verify_v3.py` |
| Stage-2 calibration-recovery scores | `results/tables/v3cal_scores_by_block.csv`, `v3cal_bss_gain_bootstrap.csv`, `v3cal_absolute_bss_ci.txt`, `v3cal_binary_by_block.csv`; fig `v3cal_skill_by_lead.png` | `20f_recover_calibration.py`, `20g_verify_v3cal.py` |
| Stage-2 out-of-sample predictions (raw) | `outputs/loyo_pred_v3sel/`, `outputs/loyo_pred_v3cal/` | `20d`, `20f` |
| Stage-3 per-model zoo binary scores + best-by-block | `results/tables/v3zoo_scores_by_model_block.csv`, `v3zoo_best_model_by_block.csv` | `21_run_loyo_v3zoo.py`, `25_verify_v3zoo.py` |
| Stage-3 13-member stack scores + diversity diagnostics | `results/tables/v3_stack_scores_by_block.csv`, `v3_stack_diversity_by_block.csv`; fig `v3_stack_by_lead.png` | `23_v3_diversity_stack.py` |
| Stage-3 out-of-sample predictions (raw): zoo + spatial DL | `outputs/loyo_pred_v3zoo/pred_v3zoo_*.npz`, `outputs/loyo_pred_v3dl/pred_v3dl_*.npz` | `21_run_loyo_v3zoo.py`, `22b_run_dl_v3.py` |
| Stage-4 full-field cache (55 channels, 20 yr) | `outputs/fullfield_cache/ff_*.npz` | `27_build_fullfield_cache.py` |
| Stage-4 true-spatial ffresunet scores + head-to-head vs box-index | `results/tables/v4ff_scores_by_model_block.csv`, `v4ff_best_model_by_block.csv`, `v4_vs_v3_resunet.csv`; fig `v4ff_skill_by_lead.png` | `28_run_loyo_ffnet.py`, `25_verify_v3zoo.py` |
| Stage-4 out-of-sample predictions (raw) | `outputs/loyo_pred_v4ff/pred_v4ff_*.npz` | `28_run_loyo_ffnet.py` |