C08 · Forecasting atmospheric CO2 — Bayesian Modelling in Practice

The brief

A climate-policy analyst is preparing a briefing. Two numbers will go on the first slide:

When will the monthly mean CO2 concentration at Mauna Loa first exceed 450 ppm? Not a single date - a distribution she can quote as "most likely X, almost certainly by Y".
How fast is CO2 rising now, compared with the 1960s? In ppm per year, with uncertainty.

She has been burned before by forecasts with error bars that turned out to be decorative. So before she trusts your answer you must show that your model, fitted without the last ten years, would have forecast those ten years well - and that its intervals meant what they said.

This is the most famous time series in climate science and it looks easy: a smooth rise with a seasonal wiggle. The obvious GP model for it has a flaw that a ten-year hold-out barely exposes, and the model that fixes the flaw will not sample until you think hard about which component is allowed to explain what.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task3") reveals hints one level at a time: nudge → approach → code skeleton.
h.check("task5", rmse=...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done.
Assumes you have worked through E05 (GPs and HSGP).

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import preliz as pz
import pymc as pm
import pytensor.tensor as pt
from pymc.model.transform.optimization import freeze_dims_and_data

from pymc_challenges import Hints, data

RANDOM_SEED = 1958
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")

h = Hints("C08")
h.tasks()

C08 - Mauna Loa CO2: a forecast that knows what it does not know

task0 - Decompose it by eye (3 hints, 2 check(s))
task1 - Hold out ten years and set the bar (3 hints, 2 check(s))
task2 - The obvious GP (3 hints, 2 check(s))
task3 - A structural model, first attempt (3 hints)
task4 - Tame it with priors (3 hints, 2 check(s))
task5 - Hold-out evaluation (3 hints, 2 check(s))
task6 - Answer the analyst (3 hints, 3 check(s))

The data

One row per month. average is the monthly mean in ppm and decimal date is the middle of the month as a fractional year. deseasonalized is NOAA's own seasonal adjustment - you may look at it, but your model must work from average. NOAA files use negative numbers as missing-value codes; find out which columns are affected before you trust any of them.

data.describe("co2")
co2 = data.load("co2").rename(columns={"decimal date": "t", "average": "ppm"})
co2.tail()

co2: Monthly mean atmospheric CO2 (ppm) from March 1958 to the present.
  source: NOAA Global Monitoring Laboratory, Mauna Loa monthly mean CO2.
  url:    https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.csv

	year	month	t	ppm	deseasonalized	ndays	sdev	unc
817	2026	4	2026.2917	431.12	428.61	23	1.16	0.46
818	2026	5	2026.3750	432.34	429.06	19	0.64	0.28
819	2026	6	2026.4583	431.43	428.99	19	0.27	0.12
820	2026	7	2026.5417	429.13	428.78	21	0.61	0.26
821	2026	8	2026.6250	427.55	429.51	17	0.36	0.17

THRESHOLD = 450.0  # ppm
N_TEST = 120  # months held out: the last ten years

Task 0 · Decompose it by eye

Deliver

A check for missing-value codes. Is ppm affected?
A plot of the full series, and a close-up of a few years that shows the seasonal cycle.
The average seasonal profile (12 numbers) after removing the trend in a simple, model-free way. In which month does CO2 peak, and what is the peak-to-trough range? Is that range the same in the 1960s and the 2010s?
A rough growth rate (ppm/yr) per decade from annual means. This is the sanity check for everything your model says later.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

NOAA marks missing values with negative numbers (-1, -9.99, -0.99, sometimes -99.99). CO2 concentrations, day counts and standard deviations cannot be negative, so one comparison per column tells you where the codes are.

For the seasonal profile you need to remove the trend without a model. What is the average of any 12 consecutive months of a purely seasonal signal?

2 Approach

A centred 12-month moving average contains (almost) no seasonal signal, so ppm - rolling mean is "seasonal + irregular". Average it by calendar month with groupby(co2.month) to get the profile; repeat on the 1960s and the 2010s separately to see whether the cycle has changed.

For growth rates, annual means remove the seasonal cycle: groupby("year").mean().diff() is ppm per year, and you can average that by decade.

3 Skeleton

print((co2.drop(columns="t") < 0).sum())

detrended = co2.ppm - co2.ppm.rolling(12, center=True).mean()
profile = detrended.groupby(co2.month).mean()
peak_month = profile.___()                       # index of the maximum
seasonal_range = profile.max() - profile.___()

annual = co2.groupby("year").ppm.mean()
growth = annual.diff()
growth.groupby((growth.index // 10) * 10).___()

Self-checks in the notebook (h.check("task0", …)):

peak_month — calendar month (1-12) in which the average seasonal profile peaks
seasonal_range — peak-to-trough range of the average seasonal profile, a little over 6 ppm

Task 1 · Hold out ten years and set the bar

Split off the last N_TEST months as a test set. You may not touch it for fitting until Task 6.

A model is only interesting if it beats something trivial. Deliver the hold-out RMSE of two baselines:

seasonal naive: every future month equals the same calendar month in the last training year;
seasonal naive with drift: the same, plus k times the average annual increase over the last ten training years for a month that lies k years (1, 2, ...) after its source month. Measure the increase as the change in the 12-month mean between the last training year and the year ending ten years before it, divided by ten.

train, test = co2.iloc[:-N_TEST], co2.iloc[-N_TEST:]
print(f"train: {train.t.iloc[0]:.2f} - {train.t.iloc[-1]:.2f} ({len(train)} months)")
print(f"test:  {test.t.iloc[0]:.2f} - {test.t.iloc[-1]:.2f} ({len(test)} months)")

train: 1958.20 - 2016.62 (702 months)
test:  2016.71 - 2026.62 (120 months)

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Both baselines are a few lines of NumPy - no model. The only fiddly part is the indexing: test month number i (0-based) should copy training month -12 + i % 12, and it lies i // 12 + 1 years after that source month.

2 Approach

Seasonal naive: last_year[np.arange(N_TEST) % 12] where last_year is the final 12 training values.

Drift: the mean of the last 12 training values minus the mean of the 12 values ending ten years earlier ([-132:-120]), divided by 10. It should come out a little above 2 ppm/yr. Add drift * years_ahead to the seasonal-naive forecast.

3 Skeleton

y_train, y_test = train.ppm.to_numpy(), test.ppm.to_numpy()
last_year = y_train[-12:]
naive = last_year[np.arange(N_TEST) % 12]
drift = (y_train[-12:].mean() - y_train[___:___].mean()) / 10
years_ahead = np.arange(N_TEST) // 12 + 1
naive_drift = naive + ___
rmse = lambda f: np.sqrt(np.mean((y_test - f) ** 2))

Self-checks in the notebook (h.check("task1", …)):

rmse_naive — hold-out RMSE (ppm) of the seasonal-naive forecast, about 15.5
rmse_drift — hold-out RMSE (ppm) of seasonal naive with drift, about 2.1

Task 2 · The obvious GP

E05 showed that a GP with a stationary kernel is a wonderful smoother. Use one as a forecaster: a single stationary kernel plus Gaussian noise, fitted to the training data, nothing else. Use pm.gp.HSGP and make sure the approximation is valid over the whole range you will predict on.

Performance note (not a modelling hint). When the inputs of an HSGP are a pm.Data container, the sine basis is recomputed at every gradient evaluation. Sample inside with freeze_dims_and_data(model): (imported above) to turn the data into constants for sampling; keep using the original model with pm.set_data for prediction. In this notebook it is the difference between one minute and half an hour.

Deliver

The fit, with diagnostics. Compare the posterior lengthscale with the length of the record. What is the model trying to tell you?
A forecast plot over the hold-out period with a 90% predictive band - and the same forecast extended to 2060.
Hold-out RMSE and 90% coverage, next to the baselines.
A short paragraph: what does the forecast do far from the data, and why would any stationary kernel do the same, whatever its hyperparameters? What does the noise term sigma represent in this model?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Set the model up exactly as in E05: inputs in a pm.Data container, an HSGP prior, a Normal likelihood with shape=f.shape. Two things differ. You will predict outside the training range, so the boundary L must enclose the forecast dates too. And the series does not hover around zero - what does a zero-mean GP do with the mean of y, and where does that leave you far from the data?

2 Approach

Centre CO2 on the training mean and measure time in years from a convenient origin. Give pm.gp.hsgp_approx.approx_hsgp_hyperparams the x-range you will predict on (up to 2060) and the lengthscale range of your prior; it returns m and c, and L = c * half-range. Pass L=[...] explicitly to pm.gp.HSGP - if you pass c, L is derived from the training inputs only. The data determine the function tightly, so parametrization="centered" samples better (E05, section 6).

Sample inside with freeze_dims_and_data(model):, then predict with the original model: pm.set_data({"X": new_x}) and pm.sample_posterior_predictive(idata, var_names=["y", "f"], predictions=True).

Coverage: for each test month take the 5% and 95% quantiles of the predictive draws and count how often the observation falls between them.

For the "why": write down the GP prior at a time point 100 years away from all data. Does it depend on the data? On the hyperparameters?

3 Skeleton

m, c = pm.gp.hsgp_approx.approx_hsgp_hyperparams(
    x_range=[train.t.min() - T0, 2060 - T0], lengthscale_range=[___, ___], cov_func="matern52")
L = c * (2060 - train.t.min()) / 2

with pm.Model() as gp_model:
    X = pm.Data("X", (train.t.to_numpy() - T0)[:, None])
    ell = pm.InverseGamma("ell", ___, ___)
    eta = pm.HalfNormal("eta", ___)            # the rise is ~100 ppm
    sigma = pm.HalfNormal("sigma", ___)
    gp = pm.gp.HSGP(m=[m], L=[L], parametrization="centered",
                    cov_func=eta**2 * pm.gp.cov.Matern52(1, ls=ell))
    f = gp.prior("f", X=X)
    pm.Normal("y", f, sigma, observed=train.ppm.to_numpy() - PPM_MEAN, shape=f.shape)

with freeze_dims_and_data(gp_model):
    gp_idata = pm.sample(random_seed=RANDOM_SEED)

with gp_model:
    pm.set_data({"X": (t_new - T0)[:, None]})
    pred = pm.sample_posterior_predictive(gp_idata, var_names=["y", "f"], predictions=True)
draws = pred.predictions["y"].stack(sample=("chain", "draw")).to_numpy() + PPM_MEAN
lo, hi = np.quantile(draws, [0.05, 0.95], axis=___)

Self-checks in the notebook (h.check("task2", …)):

sigma — posterior mean of the noise sd in ppm - about 2, far more than measurement error. What is it absorbing?
rmse — hold-out RMSE (ppm) of the GP-only forecast mean; it should be WORSE than the drift baseline

Task 3 · A structural model, first attempt

Build an additive model for the training period in which every component has one job:

a trend that keeps rising when extrapolated, plus slow, smooth departures from it (decades) - nobody promised that the future follows a textbook curve exactly;
a seasonal cycle that repeats every year but whose size may change slowly;
medium-term irregular variation (a few years: El Niño, volcanic eruptions, recessions);
observation noise.

Work in sensible units (centre CO2, shift and scale time) so that priors can be read.

For this first attempt stay generic: give every GP component the same weakly informative priors for lengthscale and amplitude, as you would if you knew nothing about which component should explain what.

Deliver

The fit and its diagnostics - all of them.
Evidence (a table or a plot of the hyperparameter posterior) showing what the data cannot decide. Where do the divergences live?
One paragraph: why is this posterior hard? Think about what happens to the likelihood if one component hands some of its signal to another.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

"Additive" means mu = trend + slow departures + seasonal + irregular, each its own block of code. Decide for each block: what is the simplest form that can do this one job, and does it extrapolate the way I want? A polynomial extrapolates upwards; a zero-mean GP extrapolates to zero, which is exactly right for a departure from something.

When the fit goes wrong (it will), ask: if long gave 1 ppm of smooth signal to irregular, or to the quadratic, would the likelihood notice?

2 Approach

One workable set of blocks: a quadratic in scaled time; an HSGP with an ExpQuad kernel for slow departures; an HSGP with a Matern52 kernel for irregular variation; Fourier terms sin(2*pi*k*t), cos(2*pi*k*t) for k = 1..4 with coefficients that may drift linearly in time (that is what pm.gp.HSGPPeriodic builds internally, minus a constant column that would clash with your intercept); Normal noise.

Each HSGP needs its own m and L covering every date you will ever predict (2041). Short lengthscales need many basis functions: check how much prior variance your m keeps at the shortest lengthscale you consider, as in E05 section 6.

Generic priors: the same InverseGamma lengthscale prior (say 90% in 1-50 years) and the same HalfNormal amplitude prior for both GPs. Then look at az.summary (r_hat, ESS), count divergences, and scatter-plot pairs of hyperparameters with the divergent draws highlighted (idata.sample_stats["diverging"]).

3 Skeleton

with pm.Model() as model:
    X = pm.Data("X", (df.t.to_numpy() - T0)[:, None])
    x = X[:, 0]
    xs = x / 30
    b = pm.Normal("b", 0, [___, ___, ___], shape=3)
    quadratic = b[0] + b[1] * xs + b[2] * xs**2

    ell_long = pm.InverseGamma("ell_long", ___);  eta_long = pm.HalfNormal("eta_long", ___)
    ell_irr = pm.InverseGamma("ell_irr", ___);    eta_irr = pm.HalfNormal("eta_irr", ___)
    long = pm.gp.HSGP(m=[___], L=[___], cov_func=eta_long**2 * pm.gp.cov.ExpQuad(1, ls=ell_long)).prior("long", X=X)
    irregular = pm.gp.HSGP(m=[___], L=[___], cov_func=eta_irr**2 * pm.gp.cov.Matern52(1, ls=ell_irr)).prior("irregular", X=X)
    trend = pm.Deterministic("trend", quadratic + long)

    k = np.arange(1, 5)
    angle = 2 * np.pi * x[:, None] * k[None, :]
    fourier = pt.concatenate([pt.sin(angle), pt.cos(angle)], axis=1)
    beta_seas = pm.Normal("beta_seas", 0, ___, shape=8)
    beta_drift = pm.Normal("beta_drift", 0, ___, shape=8)
    seasonal = pm.Deterministic("seasonal", fourier @ beta_seas + xs * (fourier @ beta_drift))

    sigma = pm.HalfNormal("sigma", ___)
    mu = trend + seasonal + irregular
    pm.Normal("y", mu, sigma, observed=df.ppm.to_numpy() - PPM_MEAN, shape=mu.shape)

with freeze_dims_and_data(model):
    idata = pm.sample(random_seed=RANDOM_SEED)
div = idata.sample_stats["diverging"].stack(sample=("chain", "draw")).to_numpy()

Task 4 · Tame it with priors

Without removing a component, and without heroic sampler settings (a moderate target_accept is fine), make the decomposition identifiable through priors alone. Every prior should be one you can defend to the analyst in physical units (years, ppm).

Deliver

Your priors and the argument for each.
Clean diagnostics.
A plot of the fitted components over the training period.
Is the seasonal cycle growing? Give the peak-to-trough range at the start and at the end of the training period, with uncertainty.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Two different things are wrong and both are about priors. (1) Two components are allowed the same lengthscales, so they are interchangeable. (2) A HalfNormal amplitude prior puts its mode at "this component does not exist" - and when a component is switched off, what constrains its coefficients? Do you actually believe the slow departures might be exactly zero?

2 Approach

Give each GP a territory in years: the irregular component gets lengthscales of one to a few years (El Niño, eruptions); the long one gets "decades" - clearly longer than that, but clearly shorter than the ~60-year record, beyond which it is just another polynomial. Narrow priors are fine when you can argue for them: pz.maxent(pz.InverseGamma(), 1, 5, 0.9), or a pm.Gamma("ell_long", mu=25, sigma=5).

For the long amplitude use a prior with no mass at zero, e.g. pm.Gamma("eta_long", mu=1.5, sigma=0.5): a boundary-avoiding prior closes the funnel. A target_accept of 0.9-0.95 mops up the remaining few divergences.

For the seasonal range: the seasonal Deterministic has one value per training month and draw; pick the 12 rows of one year and take max - min per draw.

3 Skeleton

ell_long = pm.Gamma("ell_long", mu=___, sigma=___)       # decades
eta_long = pm.Gamma("eta_long", mu=___, sigma=___)       # ppm, away from zero
irr = pz.maxent(pz.InverseGamma(), lower=___, upper=___, mass=0.9, plot=False)
ell_irr = pm.InverseGamma("ell_irr", alpha=float(irr.alpha), beta=float(irr.beta))
eta_irr = pm.HalfNormal("eta_irr", ___)

with freeze_dims_and_data(model):
    idata = pm.sample(target_accept=0.95, random_seed=RANDOM_SEED)

seas = idata.posterior["seasonal"].stack(sample=("chain", "draw")).to_numpy()   # (months, draws)
rows = (train.year == 1960).to_numpy()
range_1960 = seas[rows].max(axis=0) - seas[rows].___(axis=0)

Self-checks in the notebook (h.check("task4", …)):

sigma — posterior mean noise sd of a monthly mean, about 0.3 ppm
ell_irr — posterior mean lengthscale of the irregular component, about 1.5 years

Task 5 · Would it have worked? Hold-out evaluation

Deliver

A forecast plot for the hold-out period with a 90% predictive band, against the truth.
A table with RMSE for both baselines, the stationary GP and the structural model, plus coverage of the 50% and 90% predictive intervals for the two Bayesian models.
An honest paragraph for the analyst: is the model accurate? Is it calibrated? How much can 120 months of hold-out data tell you about calibration, given how forecast errors behave over time?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

You already wrote everything you need in Task 2: predict at the test dates, summarise the predictive draws. The new part is judgement. Plot the forecast errors against time before you interpret a coverage number - do they look like 120 independent coin flips?

2 Approach

pm.set_data with the test dates, pm.sample_posterior_predictive(..., var_names=["y"], predictions=True). RMSE uses the predictive mean; coverage uses the quantiles of the predictive draws of y (which include the noise) - not of mu or trend.

Compute the lag-1 autocorrelation of the errors (np.corrcoef(e[:-1], e[1:])). With strongly autocorrelated errors, the effective number of independent checks of your interval is far below 120.

3 Skeleton

with structural_model:
    pm.set_data({"X": (test.t.to_numpy() - T0)[:, None]})
    pred = pm.sample_posterior_predictive(idata, var_names=["y"], predictions=True)
draws = pred.predictions["y"].stack(sample=("chain", "draw")).to_numpy() + PPM_MEAN  # (months, draws)
obs = test.ppm.to_numpy()
rmse = np.sqrt(np.mean((obs - draws.mean(axis=1)) ** 2))
lo, hi = np.quantile(draws, [___, ___], axis=1)
coverage90 = np.mean((obs >= lo) & (obs <= ___))

Self-checks in the notebook (h.check("task5", …)):

rmse — hold-out RMSE (ppm) of the structural model - it should clearly beat the drift baseline (2.06)
coverage90 — share of held-out months inside the central 90% predictive interval

Task 6 · Answer the analyst

Refit your structural model on all the data.

Deliver

The posterior predictive distribution of the first month in which the monthly mean exceeds 450 ppm: a plot, the median, and statements of the form "with probability p by the end of year Y". Use predictive draws of the observable monthly mean, not the trend.
The growth rate of the trend (no seasonal, no irregular component), averaged over the last ten years of data, against the same for January 1960 - January 1970, in ppm/yr, and their ratio - all with 94% intervals.
Three sentences for the slide, and one caveat the analyst must not leave out.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

"The first month above 450 ppm" is a property of one possible future. You have thousands of possible futures: each posterior predictive draw of y over the coming months is one. Compute the answer in each and look at the distribution. Why would doing this with the smooth trend instead of y give the wrong answer?

2 Approach

Rebuild the model on all of co2 and sample. Build a grid of future mid-month dates with the same convention as the data (year + (month - 0.5) / 12), e.g. from pd.period_range, out to 2040. Predict y and trend on observed + future dates in one call.

First exceedance per draw: (y_future > 450).argmax(axis=...) gives the index of the first True; check with .any() that every draw crosses before the grid ends.

Growth rates: the trend draws at two dates ten years apart, differenced and divided by the time between them - one number per draw, then summarise. Do the same for January 1960 to January 1970, and take the ratio draw by draw.

3 Skeleton

last = pd.Period(year=int(co2.year.iloc[-1]), month=int(co2.month.iloc[-1]), freq="M")
future = pd.period_range(last + 1, "2040-12", freq="M")
t_future = (future.year + (future.month - 0.5) / 12).to_numpy()
t_all = np.concatenate([co2.t.to_numpy(), t_future])

with final_model:
    pm.set_data({"X": (t_all - T0)[:, None]})
    pred = pm.sample_posterior_predictive(final_idata, var_names=["y", "trend"], predictions=True)
y = pred.predictions["y"].stack(sample=("chain", "draw")).to_numpy() + PPM_MEAN      # (dates, draws)
y_future = y[len(co2):]
first_idx = (y_future > THRESHOLD).___(axis=0)
first_t = t_future[first_idx]
np.median(first_t), pd.Series(future[first_idx].year).value_counts(normalize=True).sort_index().cumsum()

trend = pred.predictions["trend"].stack(sample=("chain", "draw")).to_numpy()
i_last = len(co2) - 1
growth_recent = (trend[i_last] - trend[i_last - 120]) / (co2.t.iloc[i_last] - co2.t.iloc[___])

Self-checks in the notebook (h.check("task6", …)):

first_exceed_median — median (over predictive draws) of the decimal date of the first month above 450 ppm, about 2033.3
growth_1960s — posterior mean growth rate of the trend, Jan 1960 - Jan 1970, about 0.9 ppm/yr
growth_recent — posterior mean growth rate of the trend over the last ten years of data, about 2.5 ppm/yr

Going further

Does the slow-departure component earn its place? Drop it, refit on the training period and repeat Task 5. What happens to the RMSE and to the coverage of the 90% interval, and why?
Rolling-origin evaluation. Refit with training data ending in 1986, 1996 and 2006 and forecast ten years each time. Does the 90% interval cover about 90% of the time across origins? Would you widen or narrow the prior on the amplitude of the slow departures?
Marginalise the coefficients. Given the hyperparameters, this model is linear and Gaussian in all ~130 coefficients, so they can be integrated out analytically, leaving NUTS a handful of hyperparameters. Implement the marginal likelihood with the Woodbury identity (cost $O(nm^2)$) and compare speed and results.
A physical covariate. The irregular component should correlate with El Niño. Find an ENSO index (e.g. NOAA's ONI), lag it by a few months, add it as a regressor and see how much of your irregular component it absorbs.

h.progress()

C08: 0/21 hints used, 13/13 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Decompose it by eye

#Task 1 · Hold out ten years and set the bar

#Task 2 · The obvious GP

#Task 3 · A structural model, first attempt

#Task 4 · Tame it with priors

#Task 5 · Would it have worked? Hold-out evaluation

#Task 6 · Answer the analyst

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Decompose it by eye

Task 1 · Hold out ten years and set the bar

Task 2 · The obvious GP

Task 3 · A structural model, first attempt

Task 4 · Tame it with priors

Task 5 · Would it have worked? Hold-out evaluation

Task 6 · Answer the analyst

Going further