How raw model output becomes a number worth scoring

The methodology page explains how forecasts are scored here. This page explains how the forecasts our own group makes are built: the statistical layer that Weather Trader runs on top of the raw models, and the layer Forecast Scores runs on top of the verification record. Both are listed method by method, with the equation as implemented, the published work it comes from, and an honest status.

• Live issued and verified in public
• Shadow running and recorded, not yet driving a published number
• Research fitted and scored offline, has not passed its gate

The chain, end to end

Post-processing is not one model, it is a pipeline, and each stage fixes a different defect. The order matters: a stage that corrects a bias the next stage would have corrected anyway buys nothing, and a weighting stage asked to compensate for an uncorrected bias spends its weights on the wrong problem.

1. The target: calibration is defined by what you settle on

Before any statistics, the question "correct against what?" has to be answered exactly, because the same forecast is biased differently against two different truths. We carry two.

• The official record — the NWS daily climate report value (via RCC-ACIS) on the station's local standard day. This is what this site scores and what the Kalshi daily-temperature contracts settle on.
• The report record — the maximum (or minimum) of every routine and special METAR over the local calendar day, in whole °F. This is what the other venues pay on.

Measured over 8,048 station-days, the official high runs 0.75 °F above the METAR maximum (the difference is 0 on 36% of days, +1 on 52%, +2 on 10%) and the official low sits 0.55 °F below the METAR minimum. The official value catches short peaks the hourly reports miss. So the gap is modelled, not assumed: a per-station probability mass function over the integer difference, estimated walk-forward and shrunk toward the pooled estimate, and a distribution is moved from one basis to the other by convolving with it. Live

The second half of the target is which forecast was available. Every fit and every score here uses the last issue of each source published at or before the decision time — publication time, not nominal cycle time, read from the file's own timestamp. A model whose grids land six hours late is scored as it was actually available. Getting this wrong is the most common way a post-processing result turns out to be unreproducible in operations.

2. Decaying-average bias correction

The cheapest useful post-processing there is, and still the backbone of NOAA's own ensemble products. For each key — station, element, source, cycle, lead — keep one running number and no history:

The 1/(n+1) branch makes the early estimate a plain running mean, so a new station or a new model version is not corrected by noise; the n/(n+k) shrinkage does the same at application time. Half-lives differ by what is being corrected, and they are fitted, not chosen: 30 days for MOS regression statistics, 10 days for a station's residual term, 7 days for daily highs and 21 for daily lows in the station superensemble. Lows need a longer memory because their errors are driven by surface decoupling that recurs on a slower rhythm.

This is the decaying-average scheme NCEP has used for the global ensemble since 2006 (Cui et al. 2012), the same idea as the bias step inside the National Blend of Models (Hamill et al. 2017), and in the literature it is a Kalman filter with a fixed gain. Live — on every source we hold, and it is also the comparator every more elaborate method has to beat. "Raw plus the station's decaying bias" is a surprisingly strong forecast, and a method that cannot beat it is not earning its complexity.

3. Updatable Model Output Statistics

Bias correction cannot tell you that today's high will be held down by a marine layer the model mixed out, or that a clear calm night will decouple. That needs regression on the model's own fields — Model Output Statistics, in the line that starts with Glahn & Lowry (1972). We run it on a dozen models (HRRR, RRFS, ECMWF IFS, ECMWF AIFS and its ensemble, GFS, ICON, GEM, UK Met Office, WeatherNext 2 and 3), at about 2,400 CONUS ASOS/AWOS stations and, for AIFS, 2,264 more airports worldwide. Three of them — HRRR-MOS, RRFS-MOS and WN2-MOS — are scored on this site against everything else.

Form

Each equation predicts the model's residual, not the observation:

Fitting the residual means that heavy shrinkage degrades gracefully to "raw plus an intercept" rather than to nonsense. One equation exists per (element, region, cycle hour, projection); daily elements use a day slot in place of the projection.

Regions, and borrowing strength

Stations are grouped into 30 k-means regions on latitude, longitude·cos(lat) and elevation, with 1 km of elevation weighted like 3° of latitude, plus three Pacific marine-layer regions (within 30 km of the coast, below 300 m, split SoCal / central-north California / Oregon-Washington). Regional coefficients are shrunk toward the national equation for the same key; the national equation is shrunk toward zero. Penalties are expressed as pseudo-rows per standardised predictor in the national predictor scale — which matters more than it sounds: a predictor that barely varies inside one region (precipitation in a dry Southern California month) otherwise earns an unbounded coefficient, and in the prototype the first rain produced a forecast of 1,300 °F. At prediction time predictors are additionally clipped to ±4 standard deviations of the national key's training rows.

Updatable

Each key stores only decayed sufficient statistics — X'WX, X'Wy, y'Wy, n. Once a day they absorb the newly verified rows and the ridge is re-solved. Nothing is ever refitted from scratch and no history has to be kept, which is what makes a dozen models × 2,400 stations × every cycle and lead tractable on one box. This is the design of the Canadian UMOS system (Wilson & Vallée 2002, 2003), for the same reason they adopted it: when the underlying model changes, a stored cross-product matrix can be aged out, while an archive-based fit has to wait for a new archive. On a known model upgrade we multiply the regression statistics by 0.35 and the station terms by 0.5, and let the forgetting carry the new model's biases in over about a week.

Predictors

Deliberately restricted to fields every ensemble member also carries, so the same equations can run on members:

• Hourly: T, T(±3 h) − T, dew-point depression, u, v, speed, gust, total and low cloud, downward shortwave, 3-hour precipitation.
• Column: the full-mix temperature minus T2m — the temperature at the lowest of 925/850/700 hPa at least 50 hPa above the surface, brought down dry-adiabatically — boundary-layer height, 925 hPa wind, the 850−500 hPa thickness proxy, and the 850 hPa dew-point depression.
• Daily: the model's own day max/min from its hourly temperature, the hour of the extreme, peak-window cloud, solar, wind and dew-point depression, day precipitation, the day-to-day change, and a bounded 0–1 flag for "the model's extreme falls outside the diurnal window". The unbounded version of that flag extrapolated wildly.

No seasonal terms and no "model versus recent observations" predictor: both measurably hurt with less than a season of record, the latter especially on post-frontal days. Live

What it is worth

From the HRRR-MOS walk-forward (2024-07 to 2026-09, 809 issue days, 90% block-bootstrap intervals by issue date): MOS beats decaying-average bias correction at every lead and element, by 0.15–0.20 °F on daily highs, 0.09–0.32 °F on lows, 0.17–0.24 °F on hourly temperature, 0.27–0.47 °F on dew point. Against the National Blend on shared station-days it is 0.06–0.12 °F better on highs and 0.31–0.37 °F better on lows. On coastal California highs — the case the marine-layer regions exist for — it is 2.10 °F against bias-only's 2.28 °F, with the interval clear of zero.

And where it loses, which is the more useful half of the table: NBM's probability of precipitation is better than ours (Brier 0.0648 against 0.0691 for 6-hour PoP). Precipitation is the element where our post-processing is weakest, and §12 says what the literature would do about it.

Event probabilities

Sky-cover categories and PoP are fitted by REEP — linear regression on the 0/1 event indicator — then clipped to [0,1] and made monotone across thresholds, so P(≥0.10") can never exceed P(≥0.01"). REEP is crude next to logistic regression, and the fix is known (Wilks 2009; Messner et al. 2014). It is on the list. Live

4. Spread: EMOS / non-homogeneous Gaussian regression

A deterministic model has no spread at all, and an ensemble's raw spread is systematically too narrow. The standard repair is EMOS, also called non-homogeneous Gaussian regression (Gneiting et al. 2005): fit the predictive mean and the predictive variance separately, the variance as a function of the ensemble's own dispersion, by minimising the continuous ranked probability score.

The closed form is what makes it practical: the fit is a smooth four-parameter optimisation, not a likelihood over a mixture. Two variants run here.

• On a raw ensemble — the per-station EPS-EMOS in the station menu, and the GEFS trader, both in native model units so that temperature offsets and variance scaling fall out for free.
• On MOS members — each ensemble member is pushed through the same MOS equations, and the deterministic run supplies the centre:
  μ  = det + w·( mean(member MOS) − det )
  σ² = c·σ_key² + d·s²         σ_key = the equation's decayed residual SD
  with (w, c, d) fitted by minimum CRPS per element and day band, and (w) by Brier score for the REEP elements. Without members, σ = a·σ_key with the same fit. Quantiles come from the Gaussian and exceedances are continuity-corrected at whole degrees.

Live. The spread of a deterministic MOS is therefore the equation's own fitted residual spread, rescaled to match the realised spread-skill relationship — which is Veenhuis's (2013) second-moment calibration, arrived at from the same direction.

The failure mode worth naming: a stated σ nobody checked

Every σ above is stated. Whether it is right is a separate measurement, and the honest way to make it is the standardised residual z = (O − μ)/σ: if σ is right, z has standard deviation 1, and the share of outcomes inside the claimed 95% interval is 95%.

Twice now that check has caught something a CRPS table hid.

• The station consensus had z-variance 1.405 — its 2σ interval covered 91% of outcomes, not 95%. The cause was one level down: the EPS-EMOS member's own z-variance was 2.81, its stated σ roughly 40% too narrow, and inverse-variance weighting trusts a stated σ and so had made that member the heaviest in the blend. The fix measures √(mean z²) per source and lead and reads every σ through it, with a residual per-lead factor on the combination. Consensus CRPS improved at every lead (1.49→1.41 at day 1, 2.57→2.50 at day 7) and that member's weight fell from the largest to near nothing.
• The dynamic blend (§6) had z standard deviation drifting from 1.07 at day 0 to 1.89 at day 7 — a nominal 95% interval covering 90%, getting worse with every lead — because its σ was an unshrunk per-key error estimate with no predictability covariate. After rebuilding it as level × shape × scale, z standard deviation is 0.98–1.05 at every lead and the 95% interval covers 94.8%.

The shape term is EMOS again, with the weighted disagreement of the blend's inputs as the dispersion covariate instead of an ensemble variance. The tails are Student-t with the degrees of freedom read off the running kurtosis of z, which lands near 5 on real station data — a Gaussian predictive law was charging too little for the 1-in-30 bust. Shadow

5. Combining models: minimum variance on a shrunk covariance

Forecast combination is the oldest result in this field (Bates & Granger 1969): the average of two unbiased forecasts beats both unless they are perfectly correlated, and the optimal weights depend on the covariance of their errors, not on their individual skill. Our station superensemble makes every point source a member — NBM hourly and 6-hourly, NDFD, ECMWF HRES and EPS, AIFS, ICON, GEM, UK Met Office, GFS, GEFS, WeatherNext 2, HRRR, RRFS, LAMP, the classic GFS and NAM MOS, MOSMIX — and combines them in three stages, each fitted walk-forward.

Three details carry most of the value:

• The covariance, not just the variance. NBM hourly and NBM 6-hourly are nearly the same forecast; so are GFS, GFS MOS and LAMP. Minimum-variance weights make them share a weight instead of voting three times. Setting the off-diagonals to zero recovers inverse-MSE weighting, and setting S to the identity recovers the equal-weight average; both are available as options, and both are scored against the full form.
• Shrinkage, twice. A 20-member covariance estimated from a few weeks of decayed station history is badly conditioned, and inverting it unshrunk produces confident nonsense. Shrinking to a pooled estimate and then toward the diagonal is the Ledoit-Wolf prescription, and it is the difference between this working and not.
• Admission. A source enters only once its forgetting-weighted residual count passes a threshold at that moment type. A new model earns its way in walk-forward; it is never guessed into the price.

Stage 2 is a pooled ridge regression of the stage-1 error on issue-time regime features at leads 0–2; stage 3 is the CRPS-fitted spread of §4. Live, where it has earned the moment. On the Forecast Scores side the same logic appears as the skill-weighted MOS Consensus: the Lead Board's measured per-lead skill becomes the weights, and the consensus is archived as issued so that it is scored as a competitor rather than as a retrospective fit.

The inverse-MSE replica, and why we keep a deliberately simple version

Brown's (2026) working paper on the ForecastEx weather market found that a bias-corrected, equal-weight average of the public station guidance comes within 0.09 °F of the market's median the evening before — most of what the market does is that correction, performed by trading. We rebuild exactly that benchmark (each source debiased by its own trailing per-station residual, equal weights, a Gaussian whose width is the combination's own trailing error spread) and then try to beat it with inverse-MSE weights, on the same contracts, the same settlement rule and the same information cutoff. Keeping the naive version alive as a measured baseline is the only way to know whether the elaborate version is doing anything. Live

6. Dynamic performance weights

A separate line of work follows Abellan & Johnson's performance-based weighting: weight each model by the inverse of its recent root-mean-square error, recomputed from the last 30 days, per station, lead and element, with no calibration layer at all.

Our replay over 128 days, 1,192 stations and 6.8 million model-station-day pairs (828,000 scored blend rows) says something worth publishing plainly, because it is not what the headline of a weighting paper suggests:

• The dynamic weights are worth 2–3% against equal weights, which matches the paper's own finding.
• Bias correction is worth 16% on highs and 27% on lows — five to ten times as much. Every model's daily-high error at our stations shares the same cold sign, and a blend cannot average away a bias its members agree on.
• Static "optimal" weights fitted once are worse than equal weights, by up to 22%. Weights have to keep moving or not exist.
• The lineage-based family cap binds in practice, but the taxonomy is wrong: the measured error correlation of ECMWF's ensemble mean with its own high-resolution run is 0.97, GFS with GEFS is 0.77, and AIFS is the most independent member of the ECMWF family at 0.29. Families should be defined by measured error correlation, not by who runs the model.

Shadow. It is not published as a model on this site yet.

7. Mixtures: keeping the disagreement

Every method above returns one Gaussian. When ECMWF's high-resolution run says 97 °F and WeatherNext says 92 °F, a Gaussian reports 94 ± 2 and throws away the one thing that was interesting. Two mixture methods address it.

Bayesian model averaging (Raftery et al. 2005) keeps each source as its own kernel:

The weights are literally market shares — the posterior probability that each source is the best one. Missing members are handled the ensembleBMA way, renormalising each row's mixture over the members it actually has (Fraley, Raftery & Gneiting 2010), and it is scored by the closed-form mixture CRPS (Grimit et al. 2006) plus a PIT reliability index. Research, and the reason is instructive: BMA passed its gate, then stopped passing it (CRPS ratio 0.99–1.01) once the stated-σ fix in §4 landed. Most of BMA's apparent edge had been correcting one member's over-narrow σ. Fix the defect at its source and the fancy method has nothing left to earn.

Camps is the same idea for display rather than pricing: a Gaussian mixture fitted to the raw ensemble members by EM with a small C kernel, up to three components chosen by AIC, components below an 8% share or closer than 1.5 °F merged, and camp identities tracked across issues so that "the warm camp has grown from 12% to 37% over four runs" is a sentence the station card can say. Live

8. Whole degrees, brackets and ladders

Settlement is an integer in a bracket, so the last step of post-processing is discretisation — and it is not a rounding detail. A 2 °F bracket is about one standard deviation wide at day 1; the difference between integrating the predictive distribution properly and rounding the mean is the difference between a price and a guess.

Where the predictive law is Student-t the variate is scaled by √((ν−2)/ν) first so that the standard deviation remains σ. Going the other way — from a venue's price ladder to a distribution — the ladder is first made monotone by pooling adjacent violators (the PAV algorithm), then differenced into an implied pmf; its 50-cent crossing is the market's own deterministic forecast, which is how the market gets scored beside us on exactly the same footing. Live

9. The day in progress: a censored nowcast

A day-ahead distribution has no time-of-day model, so at 10 p.m. it is still spreading probability over highs the afternoon never reached. Half of a daily-temperature contract's life is spent in that state. The nowcast models the settled whole-degree extreme as

• L — the running observed bound: the METAR extreme so far, plus the 6-hour max/min groups inside the window on the official basis. Adding those groups makes the bound equal the official high on 96% of days.
• J — where the settled number sits relative to that bound. Exactly zero on the report basis, so the pmf is exactly zero past the floor; a walk-forward pmf over −2…+3 on the official basis.
• R — the rest-of-window extreme: a two-component Gaussian scale mixture around the newest hourly LAMP extreme, shifted by the observed METAR−LAMP deviation and LAMP's own expected further warming, with a spread that grows with hours left, run age and the size of that deviation.

It is fitted per local-hour bucket by maximum likelihood on the exact censored probability mass function — the likelihood of "the high was 81 and we already knew it was at least 79" is not the likelihood of a continuous observation — and separately per basis and variable. Against the day-ahead stack truncated at the same floor, using only runs published by that moment, it gains 0.38–0.76 in log score per station-hour; its PIT is flat and its 80% intervals cover 80–82%, where the stack's cover 67–70% for the low. It is admitted only in the hour buckets where the log-score gain's confidence interval clears zero and the ranked probability score does not get worse. The early-morning high stays on the day-ahead stack. Live

10. Combining with the market

The weather layer never looks at a price, which is what makes its skill measurable against the market. Measured, the best estimate of the outcome is a combination of the two, because each venue revises sluggishly. Three pooling methods are fitted per venue, variable and lead bucket, walk-forward:

The score they are all fitted to is the partition log score: minus the log probability we assigned to the cell of the venue's own ladder that the settled value fell in. Scoring on the integer grid instead would charge a bracket price for our resolution inside a bracket — the unit the market resolves on is the unit to optimise. A linear pool of 10% our nowcast into the market's brackets beats the market alone; at the touch, after fees, it still loses, so the trading rule does not take same-day contracts. Both halves of that sentence are the product. Shadow

11. Downscaling

Separately from station work, six global products (ECMWF IFS, ECMWF AIFS, the AIFS ensemble mean, GFS, GEM and ICON) are bias-corrected and downscaled onto NOAA's 2.5 km URMA analysis grid, with daily extremes rebuilt from the 6-hourly fields by Hermite interpolation where the native output is coarser than the diurnal cycle. The correction is learned per grid cell against the analysis — grid-native MOS, in the line of Glahn et al. (2009) — and the mean of two downscaled products is itself issued as a product, with nothing fitted on top of it. Live

12. The verification that drives all of it

None of the above can be trusted without a scoring layer that cannot be gamed by the fitting layer, so the scores are implemented as a separate module in the MET/METplus line that NOAA's own verification uses: continuous statistics (bias, MAE, RMSE, bias-corrected RMSE, correlation, anomaly correlation), the 2×2 contingency table and its scores (POD, FAR, CSI, Gilbert skill, Heidke, Peirce), and ensemble statistics — CRPS, the Talagrand rank histogram, the Brier score with its reliability-resolution-uncertainty decomposition (Murphy 1973), reliability-diagram bins and the spread-skill ratio.

Four rules matter more than the choice of score:

• Walk-forward, always. Every published row was forecast by equations fitted only to rows that had verified before its own publication time. There is a unit test whose name is the rule: no forecast is scored by equations that saw its outcome.
• Publication time, not cycle time. A source counts toward a forecast only if its run had actually landed; the file's own timestamp decides.
• Block-bootstrap intervals. Stations share a day's weather, so a thousand station-days from one week are not a thousand independent cases. Confidence intervals resample whole issue days. The leaderboard's "level with the leader" marker is the same idea: the gap between two sources is measured on the station-days they share, its standard error is computed with each day counting once, and a gap inside two standard errors is reported as a tie rather than as a rank.
• Admission gates, declared in advance. A method is promoted when a stated score beats a stated baseline on a stated window — and when it stops beating it, it is demoted. Section 7's BMA is the example.

Live

What the literature has that we do not use yet

This is the honest half of the page: methods that are published, benchmarked and in some cases operational elsewhere, that we have not adopted. Each entry says what it would fix here, which is the only reason to adopt one.

13. Nonparametric and machine-learning distributional regression

Every predictive distribution above is Gaussian (or a Gaussian mixture) with a variance that is linear in a dispersion covariate. The literature moved past that a decade ago:

• Quantile regression forests (Taillardat et al. 2016) — a random forest that returns the conditional distribution rather than a mean, with no distributional assumption and automatic interaction handling. This is the obvious next form for the daily high: our predictors are already assembled, and the non-linearity we most obviously miss (a marine layer either arrives or it does not) is exactly what a tree finds and a ridge cannot.
• Distributional regression networks (Rasp & Lerch 2018) — a neural network predicting the parameters of a distribution, trained directly on CRPS, with station embeddings that learn the local character instead of being given 33 hand-built regions. Their result is the one to beat: it outperformed EMOS and boosting while using fewer hand-engineered features.
• Bernstein quantile networks (Bremnes 2020) — a network emitting the coefficients of a Bernstein polynomial quantile function, which is monotone by construction. Monotonicity by construction is attractive here because we currently enforce it by clipping.
• Gradient-boosted distributional trees and the systematic comparisons (Schulz & Lerch 2022) — the benchmark that says how much of the gain is the method and how much is the features.

The honest counterweight: these methods need more training data than a 30-day decaying window, and the comparisons that favour them are run on multi-year reforecast archives. Our MOS products have 2–3 years for HRRR and seven weeks for a brand-new model. The right sequence is to fit on the deep archives (HRRR, ECMWF 0.25°) first, keep the updatable ridge where the record is short, and use EUPPBench (Demaeyer et al. 2023) as the external yardstick so that our comparison is not graded on our own homework. Planned

14. Distribution-free coverage: conformal prediction

Section 4 describes two cases of an interval that claimed 95% and delivered 90%, each found by inspection after the fact and repaired by a hand-built scale factor. Conformal prediction makes that a guarantee instead of a diagnostic: given exchangeable data, a conformalised interval has the nominal coverage by construction, with no distributional assumption at all. The directly applicable forms are conformalised quantile regression (Romano et al. 2019), which wraps any quantile model, and adaptive conformal inference (Gibbs & Candès 2021), which keeps coverage under distribution shift by adjusting the nominal level online — which is precisely our situation, since a model upgrade is a distribution shift.

Two caveats we should state rather than discover: weather data is not exchangeable (it is strongly serially dependent, and conditioning on the day fixes some but not all of that), and marginal coverage is not conditional coverage — an interval can cover 95% overall while covering 80% on the hard days, which is exactly the defect the σ-quintile reliability check in §4 was built to catch. Conformal methods would replace our hand-tuned κ with something provable; they would not remove the need for the reliability check. Planned

15. Isotonic distributional regression for the ladder

We already use the pool-adjacent-violators algorithm to make a price ladder monotone (§8). Isotonic distributional regression (Henzi, Ziegel & Gneiting 2021) is the same mathematics turned into a forecasting method: it estimates the full conditional distribution under an order constraint, is automatically calibrated, is simultaneously optimal for a broad class of loss functions, and has no tuning parameters. For "given our point forecast and today's dispersion, what is the distribution of the settled whole degree?" it is close to an ideal final calibration layer, and it is cheap. Planned

16. Precipitation, done properly

This is our largest acknowledged gap. Rain is scored on this site but barely post-processed: PoP is REEP, amounts are not calibrated at all, and the dynamic blend explicitly refuses to touch precipitation because a Gaussian blend of a zero-inflated quantity is meaningless. NBM beats us on PoP and we should expect it to. The literature is specific:

• Censored, shifted gamma distributions (Scheuerer & Hamill 2015) — one parametric family that handles the point mass at zero and the long right tail together, fitted by CRPS. This is the method inside NBM's precipitation guidance and the natural target.
• Extended logistic regression (Wilks 2009) and its heteroscedastic form (Messner et al. 2014) — a coherent set of threshold exceedance probabilities from one fit, which is what our clipped-and-monotonised REEP is a crude stand-in for.
• Generalised extreme value forms (Scheuerer 2014) for the amount distribution, and quantile mapping against the analysed climatology (the NBM approach, Hamill et al. 2017) where a long record exists.
• Neural and hybrid CSGD models for the conditional distribution where a deep reforecast archive is available.

We already run area-weighted precipitation verification on native grid cells over counties and states, with sufficient statistics that add across regions and time windows — so the scoring side for this is built and waiting. Planned

17. Coherent multivariate scenarios

Everything above calibrates one station, one element, one day at a time. The marginals are then right and the dependence is wrong — and several things we want to price are functions of the joint distribution: degree days over a week, a multi-day contract, a portfolio across cities, the probability that a cold spell holds for three days. Independent marginals will underestimate the tails of every one of those.

The standard repairs reorder calibrated marginals to inherit a realistic dependence structure: ensemble copula coupling (Schefzik, Thorarinsdottir & Gneiting 2013), which takes the rank order from the raw ensemble; the Schaake shuffle (Clark et al. 2004), which takes it from historical observations; and member-by-member post-processing (Van Schaeybroeck & Vannitsem 2015), which preserves it by construction. Scoring them requires multivariate proper scores — the variogram score (Scheuerer & Hamill 2015) or the energy score — because a marginal score cannot see a dependence error at all. This is the prerequisite for doing degree-day and energy products honestly rather than by assuming independence. Planned

18. Borrowing strength across stations

Our hierarchy is hand-built: 33 geographic regions, each shrunk toward a national equation. The literature does this by measured similarity instead of by geography. Similarity-based semi-local estimation (Lerch & Baran 2017) pools stations by the resemblance of their forecast-error characteristics and their climatology, rather than by distance; spatially adaptive EMOS (Scheuerer & Büermann 2014; Feldmann et al. 2015) interpolates the calibration parameters themselves across space. Both matter here for two concrete reasons: the thin keys in a 1,192-station network where some station-lead-cycle combinations have very few verified cases, and the gridded ambition of putting a calibrated card on thousands of places that have no settlement record of their own. Planned

19. One objective instead of a pipeline of heuristics

Stages 2 through 5 are fitted one after another, each by its own criterion, with shrinkage constants chosen by hand. That is defensible — each stage is interpretable and auditable, which matters when a number is being traded — but it is not optimal, and an external methodological review of our own work said so directly. The alternative is to fit the combination as a single constrained optimisation of one proper score:

with S chosen to match the product: CRPS for a continuous forecast, Brier or log score for a threshold contract, threshold-weighted CRPS (Gneiting & Ranjan 2011) when the money is in one tail, a multivariate score when the target is a joint quantity. Bayesian hierarchical stacking (Yao et al. 2022) is the version that keeps the shrinkage explicit and estimates it by rolling validation rather than by choosing it. A dependence penalty on correlated members replaces the hand-set family cap of §6. Planned

20. Analogs and reforecast training

We have an analog engine for the synoptic-scale question ("the last N times the flow looked like this"), using standardised Euclidean distance over issue-time features — a diagonal Mahalanobis metric, chosen because with a few hundred archived cases and 35 feature dimensions a full covariance cannot be estimated without overfitting it. What we do not have is the station-level analog literature: the analog ensemble (Delle Monache et al. 2011, 2013), which builds a predictive distribution at a point from the verifying observations of the most similar past forecasts at that same point, and reforecast-based calibration (Hamill & Whitaker 2006; Hamill et al. 2006; Hagedorn et al. 2008), which buys the long training record that every method in §13 wants. The analog ensemble is particularly attractive for rare events, where a parametric tail is an extrapolation and an analog is an observation. Planned

21. When the model arrives already calibrated

Post-processing exists because models are biased and underdispersed. The machine-learning forecast systems now in operation change that premise, and we score several of them (ECMWF AIFS and its ensemble, Google WeatherNext 2 and 3) alongside the physical models.

• Deterministic AI models (GraphCast, Lam et al. 2023; Pangu-Weather, Bi et al. 2023; AIFS, Lang et al. 2024) are trained to minimise a mean error, which makes them smooth — "blurry" in the literature's word — and a smooth field is flattered by some scores and penalised by others. Our own measurement is that AIFS has the most independent error of the ECMWF family (correlation 0.29 against its siblings' 0.97), which makes it valuable in a blend for a reason that has nothing to do with its own ranking.
• Generative and CRPS-trained ensembles (GenCast, Price et al. 2025; AIFS-CRPS, Lang et al. 2024; NeuralGCM, Kochkov et al. 2024) are trained directly on probabilistic scores, and arrive with calibration that a statistical layer was invented to supply. The open question is empirical and local: does a decaying-average station bias correction still improve a CRPS-trained ensemble at a specific airport? Our framework answers that question by construction, because the comparator is always "raw plus the station's bias" — and the answer may well be different for different models and different elements.
• Neural post-processing of fields (Grönquist et al. 2021) is the complementary direction: learn the correction on the grid, with spatial structure, rather than station by station.

Research

22. Verification upgrades we owe ourselves

• CORP reliability diagrams (Dimitriadis, Gneiting & Jordan 2021) — reliability curves estimated by isotonic regression instead of by arbitrary bins, with a consistent decomposition and uncertainty quantification. Our current reliability bins are chosen by hand, and a bin edge can change the story.
• Fair ensemble scores (Ferro 2014) — with 50 ECMWF members against 31 GEFS members, an ordinary finite-member CRPS rewards the larger ensemble for sampling density rather than for skill. Any head-to-head between ensembles of different size needs the fair version or a controlled member count.
• The spread-skill correction (Fortin et al. 2014) — the ratio of ensemble spread to error has a known finite-member factor; the naive ratio diagnoses underdispersion that is not there.
• Blocking by synoptic episode — our bootstrap resamples issue days, which handles stations sharing a day but not a ten-day blocking pattern contributing dozens of correlated cases. Resampling by episode or by week widens the intervals honestly.
• Threshold-weighted and quantile-weighted scores (Gneiting & Ranjan 2011) — CRPS is dominated by the centre of the distribution. When we repaired the σ of the dynamic blend, the tail behaviour improved dramatically and CRPS moved 1–2%, because CRPS was never the thing that was broken. A score that weights the region a contract pays in would have seen it.
• Structural breaks at model upgrades — version, resolution and ensemble size should be stored and conditioned on, not averaged across. We do decay old cases on a known upgrade; we do not yet estimate the offset between versions.

Planned

The order we intend to work in

Ranked by expected gain per unit of effort, with the gate each one must pass. A method is promoted when it passes on a stated window and demoted when it stops.

References

Given as author, year, title and journal. Follow the DOI where one is listed. Where a working paper or conference abstract is the only published form, it is labelled as such.

Last revised 21 October 2026. Statuses in this page are the state of the code, not an aspiration: a method labelled live is issuing numbers that are scored in public, and a method labelled research has either not passed its gate or has stopped passing it. Where a measured gain is quoted it comes from our own walk-forward record on the window stated beside it, and the window is usually short.