4.10 Time Series Forecasting

Forecasting means predicting future values of a time series from its past. Geoscience runs on forecasts: atmospheric CO₂, river discharge, glacier flow, groundwater levels, space weather. This notebook is a forecasting shootout. We take the canonical model families you have met in this book, point them all at the same two real data sets, evaluate them with the same metrics on the same temporal split, and let a comparison table decide.

The contestants, in order of increasing complexity:

1. Baselines: naive persistence and seasonal naive. Always first.
2. SARIMA: a classical statistical model (statsmodels).
3. Gradient boosting on lag features (lightgbm): classic ML from Chapter 3, adapted to time.
4. A small LSTM (torch): the recurrent architecture from notebook 4.4.
5. A small transformer encoder (torch): the attention pattern from notebook 4.4.

The data:

• Mauna Loa monthly CO₂ (NOAA Global Monitoring Laboratory): long, regular, strongly seasonal, with a smooth trend. The friendly case, and the primary benchmark.
• Jakobshavn Isbrae surface speed (Greenland): short, irregularly sampled, noisy. The hard case.

Learning goals

• Split time series data temporally and explain why random splits leak future information.
• Build lag features that use only the past.
• Evaluate forecasts with MSE, RMSE, MAE, MAPE, and MASE, computed in code.
• Turn a point forecaster into a probabilistic one with the pinball loss, then verify it: empirical coverage of the 90% interval, and CRPS.
• Measure how forecast error grows with lead time and locate the skill horizon.
• Quantify how much a model ranking moves across random seeds, and score rare events separately from the bulk.
• Judge when model complexity pays off and when a baseline wins.

The notebook ends with a class leaderboard in two tracks: a public CO₂ forecast that is a diagnostic only (Section 7 explains why it is trivially gameable, and why we say so in print), and a hidden synthetic series that carries the grading weight.

🖥️ Lecture slides — Session 26 (Wed Dec 2)

1. Data¶

1.1 Mauna Loa monthly CO₂¶

The Mauna Loa record is the longest continuous direct measurement of atmospheric CO₂, started by Charles Keeling in 1958. NOAA GML distributes monthly means in a text file: about 42 lines of # comments, one header line, then whitespace-delimited columns (year, month, decimal date, monthly average in ppm, deseasonalized value, number of days, standard deviation, uncertainty). We fetch it with pooch, which caches the file locally.

1.2 The leaderboard holdout: hands off the last 12 months¶

Before any modeling, we remove the final 12 months of the record: 2023-09 through 2024-08. These months are the class leaderboard holdout (Section 7) — a test set we will score exactly once. Nothing in the modeling sections touches them: no plot, no fit, no metric. We keep only their dates.

A smooth accelerating trend plus a sawtooth annual cycle (northern-hemisphere vegetation draws down CO₂ each summer). Highly predictable structure. Any model that cannot beat “repeat last year’s shape” here should embarrass its author.

1.3 Jakobshavn Isbrae surface speed¶

Jakobshavn Isbrae in west Greenland is one of the fastest glaciers on Earth. The file holds satellite-derived surface speed (m/yr) for a 10 x 10 pixel patch near the terminus, sampled irregularly (every 6 to 12 days) since 2015. Missing values are coded as -1.0; we mask them to NaN and average the valid pixels at each date into one mean-speed series.

Summer speed-ups ride on top of multi-year swings tied to terminus position and ocean forcing. Compared with CO₂ this series is short (about 430 samples over nine years), irregularly sampled, and its variability is not dominated by a clean repeating cycle.

2. Evaluation protocol¶

2.1 Temporal splits, identical for every model¶

A forecast model is only tested fairly on data from after everything it was trained on. A random train/test split of time series rows lets the model interpolate between neighbors of the test points, which inflates scores and predicts nothing. Every model below uses exactly these splits:

• CO₂: TEST = the final 48 months of the working series, 2019-09 through 2023-08. TRAIN = everything before, 1958-03 through 2019-08.
• Ice: TEST = the final 20% of dates. TRAIN = the first 80%.

All models forecast the full test window in one go (a 48-month horizon for CO₂), with no peeking at test values along the way.

2.2 Metrics¶

With forecast errors $e_t = y_t - \hat{y}_t$ over the $n$ test points:

• MSE $= \frac{1}{n}\sum e_t^2$. Penalizes large errors quadratically; units are squared, so it is awkward to read.
• RMSE $= \sqrt{\mathrm{MSE}}$. Back in data units; still dominated by the largest errors.
• MAE $= \frac{1}{n}\sum |e_t|$. Data units, treats all errors linearly. Easiest to interpret.
• MAPE $= \frac{100}{n}\sum |e_t / y_t|$. Unit-free percent, but it blows up near zero values and is asymmetric (over-forecasts and under-forecasts are penalized differently). Fine for CO₂ (~420 ppm) and ice speed, useless for series that cross zero.
• MASE (mean absolute scaled error): MAE divided by the MAE of the naive one-step forecast ($\hat{y}_t = y_{t-1}$) on the training data. MASE $< 1$ means the model beats naive persistence at its own one-step game; MASE $> 1$ means a “repeat the last value” rule was better. It is unit-free, defined at zero, and comparable across series, which is why the class leaderboard uses it.

One helper computes all five. Every number in the comparison tables comes out of this function. Never type metric values into a table by hand; the 2024 edition of this notebook did, and the table was wrong.

3. The shootout on CO₂¶

3.1 Baselines first, always¶

• Naive persistence: repeat the last training value across the whole 48-month horizon.
• Seasonal naive: repeat the last 12 training months, tiled four times. This copies last year’s seasonal shape but freezes the trend.

If a model cannot beat these, its complexity bought nothing.

3.2 SARIMA¶

SARIMA (seasonal ARIMA) models the differenced series with autoregressive and moving-average terms, plus seasonal counterparts at lag 12. We use the standard textbook order for Mauna Loa, $(1,1,1)\times(1,1,1)_{12}$: one regular difference for the trend, one seasonal difference for the annual cycle. We cap the optimizer iterations to keep the fit fast; this is a benchmark entry, not a tuning exercise.

3.3 Gradient boosting on lag features¶

Tree ensembles know nothing about time. We hand them time explicitly, as a table where each row predicts one month from features built only from earlier months:

• lags 1, 2, 3, 6, 12, 24 (the value $k$ months before the target),
• the rolling mean of the 12 months before the target (i.e. shifted by one step so the target itself is excluded),
• the calendar month (known in advance, so it is fair game).

One subtlety: a regression tree predicts averages of training targets, so it can never output a value outside the range it saw in training. CO₂ keeps rising past every training value, so a tree that predicts the level saturates near the training maximum and the forecast goes flat. The fix is to predict the one-month change and add it to the previous value; monthly changes are roughly stationary and sit well inside the training range. The features stay the same past-only lags.

The model predicts one step ahead. To cover the 48-month test window we forecast recursively: predict the next change, add it to the last value, append that prediction to the history, build features again, and so on. Errors compound, which is honest; the model gets no real observations from the test period.

3.4 A small LSTM¶

The deep models forecast directly: one forward pass maps a 36-month context window to all 48 horizon values at once (no recursion). Two practical points:

• Anchoring. CO₂ trends upward forever, so test windows sit at levels the network never saw in training. We subtract the last context value (“anchor”) from each window, train on the anchored shapes, and add the anchor back to the prediction. The network learns change relative to now, which is stationary enough to generalize.
• Scaling. Anchored windows are divided by one global standard deviation from the training windows, keeping inputs O(1).

The model is small and trains for a few epochs; raise both on your own machine if you want.

3.5 A small transformer encoder¶

Same windows, same anchoring, same direct multi-step head; only the sequence model changes. We reuse the pattern from notebook 4.4: project each scalar to a small embedding, add sinusoidal positional encoding (attention itself is order-blind), run two encoder layers, average over time, and map to the 48 horizon values.

3.6 The comparison, computed honestly¶

The table below is built programmatically from results_co2, the dictionary that register() filled from the actual predictions. No number in it was typed by hand.

What to look for in your table and plot:

• Naive persistence ignores both trend and season; its flat line is the floor everyone should beat, and its MASE is large because a 48-month horizon is much harder than the one-step game MASE is scaled by.
• Seasonal naive copies the annual shape but freezes the trend, so it drifts below the observations by roughly the trend times the horizon. Still, it beats several trained models. That is the recurring lesson of forecasting benchmarks.
• SARIMA encodes exactly the two things this series is made of, a differenced trend and a 12-month cycle, in a handful of parameters. On this kind of series it is very hard to beat.
• LightGBM, LSTM, transformer must learn trend and season from examples. In a typical run all three reach MASE below 1 and land between the seasonal baseline and SARIMA; recursive feedback (LightGBM) and limited training data (the deep pair) both cost accuracy on a 4-year horizon. Model capacity is not the bottleneck; problem structure is.

3.7 From point forecasts to intervals: the pinball loss, coverage, and CRPS¶

Every model so far answers “what will CO₂ be?” with one number per month. No serious user of a forecast settles for that: a reservoir operator, a carbon-budget analyst, and every operational weather center ask for the range of outcomes and how probable each is. The machinery is a small change to what we already have.

A model learns the quantile $q_\tau$ of the outcome when it is trained on the pinball loss (quantile loss)

which charges τ per unit of under-prediction and $1-\tau$ per unit of over-prediction; the minimizer of its expectation is exactly the τ-quantile. Predict $q_{0.05}$ and $q_{0.95}$ and you have a 90% prediction interval; $q_{0.5}$ is a median point forecast for free.

We give the LSTM from 3.4 a quantile head: instead of 48 numbers it emits $48 \times 19$ — one per horizon step per quantile level $0.05, 0.10, \dots, 0.95$ — trained on the pinball loss averaged over levels. Same windows, same anchoring, same epochs. That buys us a verifiable claim, and verification is the point:

• Empirical coverage: over the 48 test months, how often does the observation actually fall inside the 90% interval? A calibrated forecast covers about 90%. Less means the model is overconfident; much more means the intervals are too wide to inform any decision. This is the same calibration discipline as the reliability curves of Chapters 3 and 4.5, applied to intervals.
• One wrinkle: quantile outputs trained independently can cross ($\hat q_{0.40} > \hat q_{0.45}$ for some month). Sorting along the quantile axis at prediction time is the standard one-line repair, and we apply it below.

Say the number plainly: in this run the 90% interval catches 44 of 48 test months — 91.7% empirical coverage against a nominal 90%. That is a calibrated interval, and it was not guaranteed: quantile regression aims at calibration but nothing enforces it out of sample, which is why you compute coverage instead of asserting it. Two honest caveats. First, 48 consecutive months from a single forecast origin are strongly correlated, so this coverage estimate is itself noisy — a different test window could easily return 80% or 100%. Second, look at the printed widths: the interval grows from about 2.3 ppm at one month to about 6.3 ppm at four years. The model learned from the training windows that uncertainty compounds with lead time; nobody told it.

CRPS. Coverage checks one interval at one level. The continuous ranked probability score (CRPS) — the standard scoring rule of forecast verification, the number atmospheric ensemble systems report — scores the entire predictive distribution $F$ against the outcome $y$:

It has the units of the data (ppm here), and for a forecast collapsed to a single point it reduces exactly to the absolute error — so mean CRPS is directly comparable to a point forecast’s MAE, and it rewards both calibration and sharpness at once. The second identity is also the estimator: average the pinball loss over a dense quantile grid and double it. We already predict 19 quantiles, so the implementation is a few lines.

Read the three numbers together. The quantile LSTM’s mean CRPS (0.766 ppm) beats the MAE of its own median forecast (1.051 ppm) and of the point-forecast LSTM from 3.4 (0.886 ppm): spreading probability over the range where the outcome may fall earns real score even when the center is imperfect. That comparison is the entire argument for probabilistic forecasting, in one line of output. But SARIMA’s point forecast still wins at 0.281 ppm — a well-placed point beats an honestly wide distribution. Calibration does not excuse a biased or diffuse center; the verification community’s slogan is sharpness subject to calibration, and the natural next contestant would be SARIMA’s own forecast distribution (get_forecast(...).conf_int() gives its intervals — try scoring them).

3.8 Error vs lead time: the skill horizon¶

A single 48-month MAE flattens an entire dimension of forecast quality: how fast skill decays with lead time. Atmospheric verification never reports one number; it reports error as a function of lead, because a model that is excellent at 1 month and useless at 12 serves different decisions than one that is mediocre everywhere.

We measure it with a rolling-origin evaluation, the way operational systems are scored. The models stay frozen exactly as trained on data through 2019-08 — no retraining. The forecast origin then slides month by month through the test window; at each origin the model sees the observed history up to that month (an operational forecaster always has the latest observations) and issues forecasts for leads 1 through 12. Averaging the absolute error at each lead over all origins gives an error-vs-lead-time curve — with 37 origins per lead rather than the single sample per lead that one 48-month forecast provides.

The vocabulary to take away: the skill horizon is the lead time beyond which a model stops beating the reference forecast (persistence, here). Beyond it, the model adds nothing over “no change”. Weather models hit their skill horizon around two weeks; the question is where these models hit theirs on CO₂.

The curve and the printed check answer the skill-horizon question, and on CO₂ the answer is the opposite of weather. Persistence degrades steadily with lead — the trend and the seasonal cycle pull the series ever further from “no change” — while the trained models’ error grows far more slowly. LightGBM and the LSTM beat persistence at every lead from 1 to 12 months: their skill horizon lies beyond a year, past the longest lead we measured. Weather models lose to climatology in about two weeks because atmospheric chaos destroys the information in the initial state; monthly CO₂ keeps its skill because the structure that carries the forecast (trend plus annual cycle) never stops applying. The skill horizon is a property of the series as much as of the model.

The reference forecasts are still instructive. Seasonal naive beats persistence only at intermediate leads (3 through 9 in this run): at leads 1-2 last month’s value is simply very close to the truth, and near lead 12 persistence lands back in the same season, closing most of the gap. So even the choice of reference depends on lead time — one more reason a single-number summary of forecast skill hides more than it shows. For a series where the skill horizon arrives almost immediately, run this same rolling-origin evaluation on the Jakobshavn series of Section 4.

3.9 Is the ranking stable? Seed spread¶

Section 4 will warn that rankings from short series are unstable across random seeds. The same admission applies right here, so we quantify it instead of leaving it as a caveat. Every number in the 3.6 table came from one training run; a table entry is a sample from a distribution over seeds, and reporting a sample as if it were the mean is the same sin as reporting a point forecast without an interval.

The seed sensitivity differs by model class, and that is part of the lesson:

• naive, seasonal naive, SARIMA are deterministic: no seed, no spread.
• LightGBM takes a random_state, but as configured here (no row or feature subsampling) the algorithm is deterministic too — the seed is an inert knob until you enable bagging. We run 5 seeds to demonstrate that, not to estimate a spread.
• LSTM and transformer are stochastic twice over: weight initialization and batch shuffling (plus dropout, for the transformer). These are the entries whose table position you should distrust from a single run. They are also the slow entries, so we use 3 seeds each rather than 5 and say so.

The spread table sharpens the 3.6 ranking rather than overturning it — and shows exactly where the randomness lives:

• LightGBM posts the identical MASE for all five seeds (std 0.000), as promised: with no row or feature subsampling the algorithm is deterministic, and “we varied the seed” is evidence of robustness only when the algorithm actually consumes the seed.
• LSTM spans 0.808-0.856 (mean 0.827, std 0.021); the transformer spans 0.327-0.403 (mean 0.356, std 0.034) — a tenth of its own mean, from initialization, batch order, and dropout alone.
• The ranges do not overlap each other, SARIMA’s deterministic 0.256, or LightGBM’s 0.939, so on this run the shootout order survives the seeds. That is a conclusion we can now state with evidence instead of hope — and it is not a general law. On the ice series below, the gaps between models are of the same size as spreads like these, which is why its ranking should not be trusted from a single run. Report mean and spread; let the reader see whether the gaps beat the noise.

4. The ice series: a harder test¶

We run the ice series through the baselines, ARIMA, LightGBM, and the LSTM only. Fewer models on purpose: the series is short (about 350 training samples), irregularly sampled, and not dominated by a clean seasonal cycle, so a seasonal SARIMA has little to grab onto and a transformer has too little data to earn its parameters. Adding them would pad the table, not the insight.

4.1 Baselines and ARIMA¶

The sampling averages about 46 samples per year, so the “seasonal naive” tiles the last 46 training samples as an approximate last-year copy; with irregular sampling this alignment is only approximate, one more reason seasonal machinery helps less here. For the classical entry we use a plain ARIMA(1,1,1), no seasonal terms.

4.2 LightGBM and the LSTM, reused¶

Same helpers, new lag sets: lags up to 46 samples (about a year) and a short rolling mean. The LSTM uses a roughly one-year context (46 samples) and forecasts the whole test window directly, as before.

Note how much larger the MASE values are than for CO₂, across every model. The one-step naive game (the MASE denominator) is easy on a slowly varying speed series, but the multi-year test window contains changes in glacier behavior that no pattern in the training window announces. When the future is not written in the past, more model does not mean more forecast. Also expect the ranking among the trained models to be unstable here: with ~340 training samples, small choices (context length, lag set, random seed) move models up and down the table — we measured exactly this effect on the CO₂ shootout in Section 3.9, and it is worse on a series this short. Report that honestly rather than shopping for the settings that look best.

5. Rare events break bulk metrics¶

The defining data problem of the geosciences is a field that varies smoothly for years and then does something extreme for three days: floods on a discharge record, eruptions on a tremor record, offsets on a strain record. Every metric in this notebook so far averages over all test samples, and quiet samples outnumber event samples by roughly twelve to one — so a model can post an excellent MAE while missing every event that matters.

To show this cleanly we need ground truth about the events, so we build the series ourselves with mlgeo_synth.inject_rare_events: eight years of a daily seasonal background, punctuated by ~6 events per year lasting 2-6 days, with amplitudes drawn from a generalized Pareto distribution (tail="gpd", shape $\xi = 0.5$). That is a power-law tail — the extreme-value regime of real flood and surge records, where the variance is infinite and the largest event in any future window is routinely several times larger than the largest one yet observed. The generator returns an amplitude column: the true drawn size of each event, exactly the ground truth a bulk metric never sees.

The forecaster is deliberately given every advantage: one-step-ahead LightGBM (yesterday’s value is always in the feature set), trained on six years, tested on the last two. Then, instead of one MAE, we stratify the test error by the ground-truth amplitude.

A naive reading of the first row says the model is excellent: MAE 0.97 on a series whose seasonal swing alone is 6 units. The stratified rows tell the real story. Quiet days score 0.62 — the noise floor, as good as one-step forecasting gets. Small events (amplitude below 10) already cost 2.5. The seven samples in mid-size events cost about 10.6, the six samples in the largest events 18.3, and the single largest event — amplitude 50.2 — is missed by 50.8: the model predicted the seasonal background and the flood happened anyway. The right panel makes it systematic: peak error per event tracks the true amplitude along the error-equals-amplitude line. What the model has never seen, it misses in full.

Both causes are structural, not bugs. Event onsets are unpredictable from lag features — no yesterday-value tells you a spike starts tomorrow — and the GPD tail ($\xi = 0.5$) guarantees that the test window contains amplitudes beyond every training example, so the model is extrapolating precisely on the samples that matter most. Meanwhile the bulk metric hides all of it: event samples are 56 of 730 test days (about one in thirteen), so the overall MAE (0.97) sits within a factor of 1.6 of the quiet-day floor while the top stratum is thirty times worse. A leaderboard scored on overall MAE would congratulate this model.

The repair is in the metric, not the model: if events are what your forecast exists for — floods, eruptions, offsets — score them separately. Stratified tables like this one, event-window recall, or amplitude-weighted losses all work; averaging over the quiet days does not. This is the same fair-evaluation discipline as the hidden leaderboard track below: what you measure is what you optimize, so measure the thing that matters.

6. Exercises¶

Exercise 1. In the CO₂ LightGBM entry, remove lag_12 and lag_24 from the lag list and rerun. What happens to the MASE, and why?

Exercise 2. The deep models forecast all 48 months in one shot (direct multi-step); LightGBM predicts one month at a time and feeds predictions back in (recursive). Name one advantage of each strategy.

Exercise 3. Rerun Section 5 twice: once with tail="uniform" (the generator’s default: amplitudes uniform on 5-20) and once keeping tail="gpd" but raising rate_per_year to 24. Which change does more for the large-event stratum, and what does that tell you about the difference between rare and heavy-tailed?

7. The class leaderboard: two tracks¶

Now the 12 months we removed in Section 1.2 come back into play — but first, an honest disclosure about what a public leaderboard can and cannot measure.

Track A — CO₂ (diagnostic).

1. Produce a 12-month CO₂ forecast for 2023-09 through 2024-08, the held-out months, using any model trained on data up to 2023-08 only. On your honor — see the box above for why cheating here only defeats the diagnostic.
2. Save it to results/forecast_<uwnetid>.csv with exactly two columns: date (as YYYY-MM) and co2_ppm.
3. Submit it via a pull request to the course book repository, https://github.com/geo-smart/mlgeo-book. CI scores each PR and posts the standings.

Scoring is MASE against the held-out months, scaled by the naive one-step MAE of the training series (the same convention as forecast_metrics in Section 2.2 — the denominator comes from the pre-holdout record, not from the holdout itself).

Track B — the hidden synthetic series (graded). The file leaderboard/synth_forecast_history.csv in the course repository holds a daily displacement series from mlgeo_synth.gnss_series — secular trend, annual and semi-annual cycles, and colored noise, like the GNSS records of Chapter 2 — generated from a private seed and regenerated every year. Forecast the 90 days after the last history date and save results/forecast_hidden_<uwnetid>.csv with columns date (YYYY-MM-DD) and disp_mm, submitted in the same pull request. The truth series is never committed; the public leaderboard lists your submission as received, and the instructor’s run produces the scored table. Same metric: MASE, scaled by the history’s naive one-step error.

The cells below write an example file for each track using the recursive LightGBM model. Replace example with your UW NetID and the forecasts with your own.

Summary¶

• Baselines first: naive and seasonal-naive forecasts set the bar, and MASE measures every model against the naive one-step game.
• One temporal split, one metric helper, one programmatic table. Random splits and hand-typed numbers are how the 2024 edition went wrong; do not repeat either.
• On long, structured series (CO₂), a well-chosen classical model (SARIMA) is hard to beat, and deep models need care (anchoring, scaling) just to compete. On short, irregular series (Jakobshavn), every model struggles, and saying so is the correct result.
• Recursive vs direct multi-step forecasting is a real design choice; we used recursive for LightGBM and direct for the LSTM and transformer.
• A point forecast is the start, not the end. A quantile head trained on the pinball loss produced a 90% interval whose empirical coverage we measured instead of assumed, and CRPS — verification’s generalization of MAE — scored the whole distribution.
• Error is a function of lead time; plot it. The skill horizon — the lead at which a model stops beating persistence — is the operational summary of a forecast system.
• A ranking from one seed is one sample. The deterministic entries (baselines, SARIMA, this LightGBM configuration) have zero spread; the deep entries move enough across seeds to swap places. Report mean and spread, not a single draw.
• Bulk MAE can look excellent while every large rare event is missed: with heavy-tailed (GPD) amplitudes, the top event stratum carried errors thirty times the quiet-day MAE, and the largest event was missed almost in full. Stratify your scoring, or your metric will hide exactly the samples that matter.
• The leaderboard has two tracks because its public holdout is public data: the CO₂ score is a diagnostic (and trivially gameable — we say so in print), while the hidden synthetic series, regenerated yearly from a private seed, carries the grading weight.

Next: submit your two leaderboard forecasts, then continue to the final project assignment (4.20). For the sequence architectures themselves, look back at notebook 4.4; for the boosting machinery, Chapter 3.