Chapter 8
Time-dependent data

Lectures 15–16 · autocorrelation, autoregression, state-space models, Gaussian processes and honest forecasts

Every model so far has treated its observations as exchangeable: shuffle the rows of the data set and nothing is lost. Time series break that assumption. Today's temperature is close to yesterday's; this month's carbon dioxide concentration is last month's plus a little. Observations that arrive in order are relatives, not strangers, and a model that ignores the ordering makes two mistakes at once. It misses real structure (trends, cycles, persistence), and it reports far more certainty than the data contain, because it counts each of a run of similar observations as an independent piece of evidence.

This chapter models the ordering itself, using two of the most important time series in science: the monthly global surface temperature since 1880, and the monthly concentration of carbon dioxide measured on Mauna Loa since 1958. Both are publicly available and updated every month, and both are at the centre of live scientific arguments in which time-series methods decide the answer.

• Lecture 15 is about autocorrelation. We will ask how fast the world is warming and whether warming has accelerated, see how autocorrelated noise inflates the apparent precision of a trend, fit regressions with autoregressive errors in brms, and examine the "hiatus" of the 2000s as a lesson in choosing windows after seeing the data.
• Lecture 16 is about latent states and smooth functions. We will decompose the Keeling curve into a level, a growth rate and a seasonal cycle with a state-space model and the Kalman filter, ask whether the growth of carbon dioxide is still accelerating and whether the 2020 lockdowns left a mark, meet Gaussian processes as priors over whole functions, and check forecasts the only honest way: by predicting a future the model has not seen.

Lecture 15 — Autocorrelation

8.1 A warming world, month by month

NASA's Goddard Institute for Space Studies combines thousands of weather stations, ships and buoys into GISTEMP (Lenssen et al., 2019), a monthly estimate of the global average surface temperature since 1880. It is reported as an anomaly: the difference, in degrees Celsius, from the average for the same month in 1951–1980. Anomalies remove the seasonal cycle (global temperature has one, because the Northern Hemisphere has more land) and leave the change.

library(tidyverse)
library(brms)
library(posterior)
options(mc.cores = parallel::detectCores())

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

url <- "https://raw.githubusercontent.com/datasets/global-temp/main/data/monthly.csv"
gistemp <- read_csv(url, show_col_types = FALSE) |>
  filter(Source == "GISTEMP") |>
  transmute(year = as.integer(substr(Year, 1, 4)), month = as.integer(substr(Year, 6, 7)),
            temp = Mean) |>                                   # anomaly vs 1951-80, degrees C
  arrange(year, month)
range(gistemp$year)
annual <- gistemp |> group_by(year) |> summarise(temp = mean(temp), months = n()) |> filter(months == 12)
annual |> slice_max(temp, n = 5)                               # the five warmest years

[1] 1880 2025
# A tibble: 5 x 3
   year  temp months
  <int> <dbl>  <int>
1  2024  1.29     12
2  2025  1.19     12
3  2023  1.17     12
4  2016  1.01     12
5  2020  1.01     12

gistemp <- gistemp |> mutate(time = year + (month - 0.5) / 12)
ggplot(gistemp, aes(time, temp)) +
  annotate("rect", xmin = 1970, xmax = 2026, ymin = -Inf, ymax = Inf, fill = harbour[["grid"]], alpha = 0.5) +
  geom_line(colour = harbour[["muted"]], linewidth = 0.2, alpha = 0.7) +
  geom_point(data = annual, aes(year + 0.5, temp), colour = harbour[["rose"]], size = 1) +
  geom_hline(yintercept = 0, colour = harbour[["muted"]], linetype = "dashed") +
  labs(x = NULL, y = "temperature anomaly (°C)")

The three warmest years in the record are 2024, 2025 and 2023, followed by 2016 and 2020; 2024 was the warmest of all, about 1.3 °C above the 1951–1980 average. Since about 1970 the rise has been steady, and the question for this lecture is how to describe it. How fast is the world warming, with what uncertainty? And has the rate increased recently, as some climate scientists have argued since the record heat of 2023 and 2024? These are regression questions, but with a twist: the observations are a time series.

8.2 Autocorrelation: why time breaks independence

Fit a straight line to the monthly anomalies since 1970 and look at what is left over. If the residuals were independent, a warm month would be followed by a warm or cool month with equal probability. They are not: warm months come in runs, often lasting a year or more, because the climate system has memory. The oceans store heat, and the El Niño–Southern Oscillation shifts heat between the ocean and the atmosphere over periods of two to seven years. The autocorrelation function (ACF) measures this memory: the correlation between the residuals and themselves $k$ months earlier, for each lag $k$.

recent <- gistemp |> filter(year >= 1970, year <= 2025) |>
  mutate(t = (time - 2000) / 10)                           # decades since 2000
ols <- lm(temp ~ t, data = recent)
recent$resid <- resid(ols)
r_acf <- acf(recent$resid, lag.max = 24, plot = FALSE)
p1 <- ggplot(recent, aes(time, resid)) + geom_hline(yintercept = 0, colour = harbour[["muted"]]) +
  geom_line(colour = harbour[["sky"]], linewidth = 0.3) + labs(x = NULL, y = "residual (°C)")
p2 <- tibble(lag = r_acf$lag[-1], acf = r_acf$acf[-1]) |>
  ggplot(aes(lag, acf)) + geom_hline(yintercept = 0, colour = harbour[["muted"]]) +
  geom_hline(yintercept = c(-1, 1) * 1.96 / sqrt(nrow(recent)), linetype = "dashed", colour = harbour[["rose"]]) +
  geom_segment(aes(xend = lag, yend = 0), colour = harbour[["sky"]], linewidth = 1.2) +
  labs(x = "lag (months)", y = "autocorrelation")
patchwork::wrap_plots(p1, p2, widths = c(1.4, 1))

Why does this matter for the trend? Because the standard error of a regression slope assumes each of the 672 months is a fresh, independent draw of noise. With autocorrelated noise, consecutive months repeat much of the same information. For a residual process with lag-one autocorrelation $\rho$ (an AR(1) process, defined in the next section), the information in $n$ observations about a trend is roughly that of

independent ones. With $\rho =0.68$, that is $672\times 0.32/1.68\approx 130$: the 56 years of monthly data carry about as much information about the trend as 130 independent observations, not 672. The naive standard error of the slope is too small by a factor of about $\sqrt{672/130}\approx 2.3$. This is the same effective-sample-size idea that measured the information in correlated MCMC draws in Chapter 5, now applied to the data themselves.

At $\varphi =0$ both intervals behave. As $\varphi$ grows, the naive interval excludes zero far more often than 5% of the time: with $\varphi =0.7$, about 40% of pure-noise series show a "significant" trend, and with $\varphi =0.9$, two-thirds. The corrected interval does much better but is still too narrow (about 12% at $\varphi =0.7$), because it plugs in an autocorrelation estimated from only 60 points, which tends to be too low; a Bayesian model that estimates $\varphi$ and the trend together, and propagates the uncertainty in $\varphi$, does better still. Autocorrelated noise wanders, and a wandering series looks like a trend over any window shorter than its wandering. This is the single most common statistical error in the analysis of time series, and it is not confined to climate: it inflates the significance of trends in prices, crime rates, disease incidence and hospital waiting times whenever the residuals are treated as independent.

8.3 Regression with autoregressive errors

An autoregressive process of order $p$, AR($p$), lets each value depend on the $p$ before it. For order one,

The coefficient $\varphi$ is the memory. If $|\varphi |<1$ the process is stationary: it keeps returning to zero, a shock decays by a factor $\varphi$ every step, and the lag-$k$ autocorrelation is ${\varphi }^k$. At $\varphi =1$ it is a random walk, which never returns: every shock is permanent. For trend estimation we put the autoregression in the errors, not the mean: the temperature is a trend plus AR(1) noise,

In brms, the term ar(time = t, p = 1) adds AR(1) errors ordered by t. We fit the trend both ways.

pri <- c(prior(normal(0.5, 1), class = Intercept), prior(normal(0, 0.5), class = b),
         prior(exponential(5), class = sigma))
m_iid <- brm(temp ~ t, data = recent, prior = pri, seed = 8, refresh = 0, file = "fits/ch8_iid")
m_ar1 <- brm(temp ~ t + ar(time = t, p = 1), data = recent, prior = pri,
             seed = 8, refresh = 0, file = "fits/ch8_ar1")
rbind(independent = fixef(m_iid)["t", ], `AR(1) errors` = fixef(m_ar1)["t", ]) |> round(3)
posterior_summary(m_ar1, variable = c("ar[1]", "sigma")) |> round(3)

             Estimate Est.Error  Q2.5 Q97.5
independent     0.202     0.004 0.194 0.209
AR(1) errors    0.200     0.008 0.184 0.217
      Estimate Est.Error  Q2.5 Q97.5
ar[1]    0.681     0.029 0.624 0.737
sigma    0.112     0.003 0.106 0.118

Both models agree on the rate of warming, about 0.20 °C per decade since 1970. They disagree on its precision: the AR(1) model's posterior standard deviation for the slope is about twice as large, as the effective sample size predicted. With 56 years of data the trend is so strong that both intervals are far from zero. The difference matters when the signal is weaker, as it is in the next question.

8.4 Has warming accelerated?

After the record warmth of 2023 and 2024, James Hansen and colleagues (2025) argued that global warming has accelerated, from about 0.18 °C per decade between 1970 and 2010 to well above that since; others, applying formal change-point methods, concluded that an acceleration was not yet statistically detectable (Beaulieu et al., 2024). The question is important (it affects how soon the world passes 1.5 °C above pre-industrial temperatures) and it is a pure time-series problem. One simple way to pose it is a hinge model: a straight line whose slope may change in 2010,

where $\gamma$ is the increase in the warming rate after 2010. We fit it with independent errors and with AR(1) errors.

recent <- recent |> mutate(h = pmax(0, t - 1))           # decades since 2010, zero before
m_h_iid <- brm(temp ~ t + h, data = recent, prior = pri, seed = 8, refresh = 0, file = "fits/ch8_h_iid")
m_h_ar1 <- brm(temp ~ t + h + ar(time = t, p = 1), data = recent, prior = pri,
               seed = 8, refresh = 0, file = "fits/ch8_h_ar1")
summ <- function(m) { d <- as_draws_df(m)
  c(slope_before = mean(d$b_t), increase = mean(d$b_h), lo = quantile(d$b_h, .025),
    hi = quantile(d$b_h, .975), `Pr(increase > 0)` = mean(d$b_h > 0)) }
rbind(independent = summ(m_h_iid), `AR(1) errors` = summ(m_h_ar1)) |> round(3)

             slope_before increase lo.2.5% hi.97.5% Pr(increase > 0)
independent         0.171    0.155   0.119    0.193                1
AR(1) errors        0.169    0.155   0.075    0.234                1

Before 2010 the world warmed at about 0.17 °C per decade; the hinge model estimates that the rate increased by about 0.15 °C per decade after 2010, nearly doubling it. With independent errors the 95% interval for the increase is 0.12 to 0.19; with AR(1) errors it is about twice as wide, 0.08 to 0.23. In both cases essentially every posterior draw has a positive increase. On this model, the evidence for acceleration is strong.

But there is a catch, and it is not about autocorrelation. Why 2010? Because looking at the graph, that is roughly where the line seems to bend. A break point chosen after looking at the data will find a bend wherever the data happen to bend most, and the uncertainty in where the break is should be part of the answer. The widget lets you choose the window and the break year, and shows how much the conclusions move.

Choose 1998 to 2012 and the warming seems to have slowed, to about 0.13 °C per decade. That window starts at a record El Niño year and ends in a run of La Niña years, and for several years it was cited as evidence that global warming had stopped. Yet the corrected interval for its slope, about 0.05 to 0.21, excludes zero and includes the long-term rate of 0.20: fifteen years of autocorrelated monthly data could not distinguish a slowdown from an unchanged rate, and they clearly ruled out a stop. The same caution applies in reverse to the years since 2011, whose slope of about 0.39 °C per decade ends with the El Niño of 2023–2024. Move the break year and the estimated acceleration changes a great deal, from about 0.09 °C per decade with a break in 2000 to about 0.22 with a break in 2015, though it stays positive throughout. The honest summary is that the evidence for an increase in the rate is real, but its size and timing are uncertain, and the record years of 2023–2024 owe part of their warmth to El Niño.

Lecture 16 — Latent states and smooth functions

8.5 The Keeling curve

The second series is the most famous graph in climate science. Since March 1958, the concentration of carbon dioxide in the air has been measured continuously at the Mauna Loa Observatory in Hawaii, 3,400 metres up and far from local sources. The monthly averages show two things at once: a relentless rise, from about 315 parts per million (ppm) to over 425, and a regular annual cycle, as the vegetation of the Northern Hemisphere takes up carbon in summer and releases it in winter.

url_co2 <- "https://raw.githubusercontent.com/datasets/co2-ppm/main/data/co2-mm-mlo.csv"
co2 <- read_csv(url_co2, show_col_types = FALSE, skip = 1,
                col_names = c("ym", "time", "ppm", "deseasonalised", "days", "sd", "unc")) |>
  mutate(year = as.integer(substr(ym, 1, 4)), month = as.integer(substr(ym, 6, 7))) |>
  select(year, month, time, ppm)
c(months = nrow(co2), first = co2$ppm[1], last = tail(co2$ppm, 1))
co2 |> filter(month == 5, year %in% c(1960, 1990, 2020, 2025))   # May is the annual peak

months  first   last 
 822.0  315.7  427.6 
# A tibble: 4 x 4
   year month  time   ppm
  <int> <int> <dbl> <dbl>
1  1960     5 1960.  320.
2  1990     5 1990.  357.
3  2020     5 2020.  417.
4  2025     5 2025.  431.

8.6 State-space models and the Kalman filter

Regression with AR errors describes a series as a fixed function of time plus correlated noise. A state-space model takes a different view: there is an unobserved state that evolves over time, and each observation is a noisy view of the current state. For carbon dioxide, the state is the "true" smoothed concentration and its rate of growth, and the observation adds the seasonal cycle and month-to-month noise. Drawn as a causal graph, the states form a chain running forward in time, and each observation hangs off its own state.

The simplest state-space model is the local level model: the state ${\mu }_t$ is a random walk and the observation is the state plus noise,

Given the observations up to time $t-1$, our belief about ${\mu }_{t-1}$ is a normal distribution. The Kalman filter (Kalman, 1960) updates it in two steps each month. Predict: push the belief forward through the random walk, which keeps the mean and adds ${\sigma }_{\mu }^2$ to the variance. Update: combine the prediction with the new observation by the normal–normal rule behind the partial pooling of Chapter 6, a precision-weighted average. The posterior mean moves toward the observation by a fraction $K_t$, the Kalman gain, which is the share of the total variance that belongs to the prediction:

That is the whole algorithm: Bayes' rule, once per time step. Running it forward gives the filtered estimates, using only the past (what an observer would have known at each moment). A backward pass, the smoother, then revises every state using the whole series. The filter also delivers, as a by-product, the probability of each observation given the past, and their product is the likelihood of the hyperparameters ${\sigma }_{\mu }$ and ${\sigma }_y$ with the states integrated out. Here is the filter in ten lines of R, applied to annual global temperature.

kalman_level <- function(y, s_mu, s_y, m0 = y[1], P0 = 1) {
  n <- length(y); m <- P <- numeric(n); loglik <- 0; mp <- m0; Pp <- P0
  for (t in 1:n) {
    if (t > 1) { mp <- m[t - 1]; Pp <- P[t - 1] + s_mu^2 }      # predict
    K <- Pp / (Pp + s_y^2)                                       # Kalman gain
    loglik <- loglik + dnorm(y[t], mp, sqrt(Pp + s_y^2), log = TRUE)
    m[t] <- mp + K * (y[t] - mp); P[t] <- (1 - K) * Pp           # update
  }
  list(m = m, P = P, loglik = loglik)
}
fit_ml <- optim(c(log(0.05), log(0.1)), function(p)
  -kalman_level(annual$temp, exp(p[1]), exp(p[2]))$loglik)
c(s_mu = exp(fit_ml$par[1]), s_y = exp(fit_ml$par[2]))           # maximum likelihood

   s_mu     s_y 
0.06473 0.06854 

The maximum-likelihood estimates are ${\sigma }_{\mu }\approx 0.065$ °C (how much the underlying level moves in a year) and ${\sigma }_y\approx 0.069$ °C (how far a single year strays from it): the two are about equal, which in steady state gives a Kalman gain of about 0.6: each new year moves the estimated level about 60% of the way from the prediction toward the observation. The widget lets you explore them.

Make ${\sigma }_{\mu }$ large relative to ${\sigma }_y$ and the gain approaches one: the filter believes each new observation almost completely, and the level follows every wiggle. Make it small and the gain falls toward zero: the filter trusts its prediction and the level becomes a slow average. The likelihood is highest in between, and that compromise, chosen by the data, is the state-space model's version of the shrinkage of Chapter 6. Notice also that the filtered estimate lags behind the smoothed one during rapid warming: in real time, an observer cannot tell a new trend from a run of warm years until several years have passed.

8.7 Is carbon dioxide still accelerating?

For carbon dioxide we want three components: a smooth level, a slowly changing growth rate (the slope of the level, so that the level can rise faster and faster), and a seasonal cycle that may drift slowly. Each evolves by its own random walk, and the observation is their sum plus noise. This is a structural time-series model (Harvey, 1989); the growth rate is the quantity of interest. Stan can fit it efficiently because it has the Kalman filter built in: the distribution gaussian_dlm_obs computes the likelihood with every latent state integrated out, so the sampler explores only the four standard deviations. The states themselves are then recovered, for each posterior draw, by the Kalman smoother in the dlm package.

dlm_code <- "
data { int<lower=2> N; matrix[1, N] y; }
transformed data {
  matrix[6, 1] F = [[1], [0], [1], [0], [1], [0]];          // level + two seasonal harmonics
  matrix[6, 6] G = rep_matrix(0, 6, 6);
  real w1 = 2 * pi() / 12;  real w2 = 4 * pi() / 12;
  G[1, 1] = 1; G[1, 2] = 1; G[2, 2] = 1;                    // level grows by the slope
  G[3, 3] = cos(w1); G[3, 4] = sin(w1); G[4, 3] = -sin(w1); G[4, 4] = cos(w1);
  G[5, 5] = cos(w2); G[5, 6] = sin(w2); G[6, 5] = -sin(w2); G[6, 6] = cos(w2);
  vector[6] m0 = [315, 0.07, 0, 0, 0, 0]';
  matrix[6, 6] C0 = diag_matrix([100, 0.01, 10, 10, 10, 10]');
}
parameters { real<lower=0> s_level; real<lower=0> s_slope; real<lower=0> s_seas; real<lower=0> s_y; }
model {
  s_level ~ exponential(10); s_slope ~ exponential(100); s_seas ~ exponential(20); s_y ~ exponential(2);
  y ~ gaussian_dlm_obs(F, G, [s_y^2]',
        diag_matrix([s_level^2, s_slope^2, s_seas^2, s_seas^2, s_seas^2, s_seas^2]'), m0, C0);
}"
fit_dlm <- rstan::stan(model_code = dlm_code, data = list(N = nrow(co2), y = matrix(co2$ppm, nrow = 1)),
                       seed = 8, refresh = 0)
print(fit_dlm, digits = 4, probs = c(.05, .95))

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

            mean se_mean     sd       5%      95% n_eff   Rhat
s_level   0.1621  0.0002 0.0118   0.1432   0.1821  3581 0.9996
s_slope   0.0029  0.0000 0.0010   0.0017   0.0048  3232 1.0007
s_seas    0.0083  0.0000 0.0016   0.0060   0.0112  3655 0.9998
s_y       0.2157  0.0002 0.0092   0.2006   0.2312  3746 0.9994
lp__    477.3341  0.0306 1.4253 474.5831 478.9962  2172 1.0023

Samples were drawn using NUTS(diag_e) at Fri Oct 23 10:57:08 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

library(dlm)
draws_dlm <- as.data.frame(fit_dlm)
build <- function(p) dlmModPoly(2, dV = p$s_y^2, dW = c(p$s_level^2, p$s_slope^2),
                                m0 = c(315, 0.07), C0 = diag(c(100, 0.01))) +
                     dlmModTrig(s = 12, q = 2, dV = 0, dW = p$s_seas^2, m0 = rep(0, 4), C0 = diag(10, 4))
set.seed(8)
ids <- sample(nrow(draws_dlm), 200)
smooth_draws <- lapply(ids, function(i) dlmSmooth(co2$ppm, build(draws_dlm[i, ]))$s[-1, ])
growth <- sapply(smooth_draws, function(s) 12 * s[, 2])          # ppm per year, months x draws
level  <- sapply(smooth_draws, function(s) s[, 1])
season <- sapply(smooth_draws, function(s) s[, 3] + s[, 5])
gr <- co2 |> mutate(g = apply(growth, 1, median), g_lo = apply(growth, 1, quantile, .05),
                    g_hi = apply(growth, 1, quantile, .95), seas = apply(season, 1, median))
gr |> filter(month == 7, year %in% c(1960, 1970, 1980, 1990, 2000, 2010, 2015, 2019, 2020, 2025)) |>
  select(year, g, g_lo, g_hi) |> mutate(across(-year, ~ round(.x, 2)))

# A tibble: 10 x 4
    year     g  g_lo  g_hi
   <int> <dbl> <dbl> <dbl>
 1  1960  0.87  0.82  0.95
 2  1970  1.1   1.09  1.14
 3  1980  1.46  1.42  1.52
 4  1990  1.51  1.41  1.56
 5  2000  1.84  1.82  1.85
 6  2010  2.16  2.13  2.16
 7  2015  2.38  2.3   2.47
 8  2019  2.44  2.36  2.48
 9  2020  2.44  2.37  2.47
10  2025  2.46  2.4   2.47

lv <- co2 |> mutate(level = apply(level, 1, median))
p1 <- ggplot(lv, aes(time)) + geom_point(aes(y = ppm), colour = harbour[["muted"]], size = 0.3) +
  geom_line(aes(y = level), colour = harbour[["sky"]], linewidth = 0.8) + labs(x = NULL, y = "CO2 (ppm)")
p2 <- ggplot(gr, aes(time)) + geom_ribbon(aes(ymin = g_lo, ymax = g_hi), fill = harbour[["sky"]], alpha = 0.25) +
  geom_line(aes(y = g), colour = harbour[["sky"]], linewidth = 0.8) + labs(x = NULL, y = "growth (ppm/yr)")
p3 <- ggplot(gr, aes(time, seas)) + geom_line(colour = harbour[["moss"]], linewidth = 0.4) +
  labs(x = NULL, y = "seasonal (ppm)")
patchwork::wrap_plots(p1, p2, p3, ncol = 1, heights = c(1.3, 1, 0.8))

The sampler explores just four parameters, and all four are well estimated ($\widehat{R}=1.00$, effective sample sizes above 3,000). The growth rate tells a clear story. It rose from about 0.9 ppm per year around 1960 to 1.1 by 1970, 1.5 by 1980, 1.8 by 2000 and 2.2 by 2010, with a pause in the early 1990s (when the eruption of Mount Pinatubo cooled the climate and the economies of the former Soviet bloc collapsed). From 2015 to 2025 it barely changed, from about 2.4 to 2.5 ppm per year, with 90% intervals of about ±0.05 to 0.1. Carbon dioxide is still rising, by about 2.5 ppm a year, but on this model it is no longer rising faster each year, which matches the levelling-off of global fossil-fuel emissions over the same decade. Two cautions. The model has decided, through the estimated standard deviations, that year-to-year wobbles belong to the level rather than the growth rate, so the growth rate is smooth by construction. And single-year increases vary much more: the increase in 2024, boosted by El Niño and wildfires, was the largest on record. Whether that marks a renewed acceleration only further years will show.

What about 2020? Fossil fuel emissions fell by about 1.9 billion tonnes of carbon dioxide that year, roughly 5%. About 2.1 billion tonnes of carbon (7.8 billion tonnes of CO₂) raise the atmospheric concentration by 1 ppm, and only about half of each year's emissions stay in the air (the rest is taken up by oceans and plants), so the expected effect on the 2020 growth rate is about $1.9/7.8\times 0.5\approx 0.1$ to $0.2$ ppm. The year-to-year variation in the growth rate, driven by the response of tropical forests and oceans to El Niño and La Niña, is of the order of 0.5 ppm. A dip of 0.15 ppm in one year is buried in noise three times its size, and the model, which separates a slowly evolving growth rate from short-term variation, shows no break in 2020. Detecting it would have required years of reduced emissions. The 2020 non-detection was correctly predicted in advance by atmospheric scientists using exactly this signal-to-noise reasoning, and it is a good example of a power calculation for a time series.

8.8 Gaussian processes

Autoregression models a series through its own past; a state-space model builds it from evolving components. A Gaussian process (GP) takes a third view: it puts a prior directly on the whole unknown function $f(t)$. The prior says that the values of $f$ at any set of times are jointly normal, with correlations given by a kernel $k(t,t^{\prime })$ that encodes what kind of function we expect. The most common is the squared exponential kernel,

under which values at nearby times are strongly correlated and values more than a few length-scales $\ell$ apart are nearly independent; $\alpha$ is the typical size of the function's departures from its mean. Draws from this prior are smooth, wiggly curves, and conditioning on data gives a posterior over curves that pass near the data. Kernels can be added (a smooth trend plus a periodic seasonal component) and multiplied, which makes GPs a modular language for describing structure.

Three lessons come out of the widget. First, the marginal likelihood (the probability of the data with $f$ integrated out) is a built-in Occam's razor: very short length-scales fit every wiggle but are penalised for being able to fit anything, very long ones cannot follow the data, and the likelihood peaks in between. Second, the band is narrow where there are data and widens beyond them. Third, and most important for forecasting: with a stationary kernel such as the squared exponential, the forecast reverts to the mean once you are more than a length-scale beyond the data, because the prior says that distant values are unrelated to the observed ones and centred on the overall mean. A GP knows only what its kernel tells it. If you believe warming will continue, the kernel must say so, for example by adding a linear trend, and the forecast then depends on that assumption rather than discovering it.

In brms, gp(x) adds a GP term with a squared exponential kernel and estimates $\ell$ and $\alpha$ with priors rather than fixing them. Exact GPs need a Cholesky factorisation of an $n\times n$ matrix at every step of the sampler, which becomes slow even for a few hundred points, so we use the Hilbert-space approximate GP of Riutort-Mayol et al. (2023): gp(x, k = 25, c = 5/4) represents the GP with 25 basis functions, which is accurate for functions that are smooth relative to the range of the data and fits in seconds.

annual <- annual |> mutate(x = (year - 2000) / 50)          # half-centuries from 2000
m_gp <- brm(temp ~ gp(x, k = 25, c = 5/4), data = annual,
            prior = c(prior(normal(0, 1), class = Intercept),
                      prior(exponential(1), class = sdgp),
                      prior(exponential(5), class = sigma)),
            seed = 8, refresh = 0, control = list(adapt_delta = 0.99), file = "fits/ch8_gp")
posterior_summary(m_gp, variable = c("sdgp_gpx", "lscale_gpx", "sigma")) |> round(3)
future <- tibble(year = c(2025, 2030, 2035, 2045)) |> mutate(x = (year - 2000) / 50)
cbind(future["year"], round(fitted(m_gp, newdata = future, probs = c(.05, .95)), 2))

           Estimate Est.Error  Q2.5 Q97.5
sdgp_gpx      0.665     0.271 0.339 1.359
lscale_gpx    0.173     0.026 0.125 0.228
sigma         0.101     0.006 0.090 0.114
  year Estimate Est.Error   Q5  Q95
1 2025     1.14      0.05 1.05 1.23
2 2030     1.28      0.11 1.10 1.46
3 2035     1.38      0.19 1.08 1.70
4 2045     1.41      0.40 0.79 2.10

The estimated length-scale is about 0.17 half-centuries, or roughly 9 years (95% interval 6 to 11), the amplitude about 0.7 °C and the residual standard deviation about 0.10 °C. The forecast shows the behaviour the widget warned about. The GP continues the recent rise for about a length-scale: 1.28 °C in 2030 and 1.38 °C in 2035. Then it flattens, and its interval widens rapidly: by 2045 the forecast is 1.41 °C with a 90% interval from 0.8 to 2.1 °C. The flattening is not a prediction that warming will stop. It is the prior reasserting itself where the data run out, because the squared exponential kernel knows nothing about greenhouse gases.

8.9 Forecasting honestly

The only convincing test of a forecast is to make it before the outcome is known. With a time series we can simulate that: fit each model to the data up to some year, forecast the following years, and compare with what happened. This is leave-future-out validation (Bürkner, Gabry and Vehtari, 2020). Ordinary cross-validation, which leaves out random observations and predicts them from both their past and their future, overstates forecasting skill for time series, because interpolating between known neighbours is much easier than extrapolating.

train <- filter(annual, year <= 2015); test <- filter(annual, year > 2015)
m_gp_15 <- update(m_gp, newdata = train, seed = 8, refresh = 0)
m_lin_15 <- brm(temp ~ x, data = filter(train, year >= 1970),
                prior = c(prior(normal(0.5, 1), class = Intercept), prior(normal(0, 2), class = b),
                          prior(exponential(5), class = sigma)), seed = 8, refresh = 0)
score <- function(m) { set.seed(8); p <- posterior_predict(m, newdata = test)
  c(rmse = sqrt(mean((colMeans(p) - test$temp)^2)), bias = mean(colMeans(p) - test$temp),
    coverage_90 = mean(test$temp > apply(p, 2, quantile, .05) & test$temp < apply(p, 2, quantile, .95))) }
rbind(`GP, 1880-2015` = score(m_gp_15), `linear, 1970-2015` = score(m_lin_15)) |> round(3)

                   rmse   bias coverage_90
GP, 1880-2015     0.283 -0.246         0.4
linear, 1970-2015 0.193 -0.152         0.6

Fitted to the data up to 2015, the GP forecasts 2016–2025 too low by about 0.25 °C on average, and only 4 of the 10 years fall inside its 90% predictive intervals. The straight line fitted to 1970–2015 does better (an average error of −0.15 °C, 6 of 10 years covered), but it too is too low and too confident. Both models were surprised by the warmth of the last decade, which is itself evidence for the acceleration of Section 8.4 (and for the El Niño of 2023–2024). Neither model is a climate projection; projections come from physical climate models driven by emission scenarios. The point is the method: a forecast evaluated on a future it did not see reveals overconfidence that no in-sample check would.

8.10 Time series in practice

The three approaches of this chapter are less rivals than tools for different questions, and a careful analysis often uses more than one. Regression with AR or ARMA errors is the natural choice when the question is about the effect of predictors (a trend, a policy, a covariate) and the time structure is a nuisance to be accounted for. State-space models are the choice when the evolving components are themselves the object of interest (a growth rate, a seasonal pattern, an underlying level), when the data arrive in real time and estimates must be updated, or when there are missing observations (the Kalman filter simply skips the update step). Gaussian processes are the choice for flexible smoothing and interpolation with honest uncertainty, and for irregularly spaced times. All three are Bayesian models fitted with the same machinery as every other model in this book, and many standard GP kernels have exact state-space representations, so the distinctions are partly a matter of computation and convenience.

8.11 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. Beaulieu, C., Gallagher, C., Killick, R., Lund, R. and Shi, X. (2024). A recent surge in global warming is not detectable yet. Communications Earth & Environment, 5, 576.
2. Bürkner, P.-C., Gabry, J. and Vehtari, A. (2020). Approximate leave-future-out cross-validation for Bayesian time series models. Journal of Statistical Computation and Simulation, 90(14), 2499–2523.
3. Durbin, J. and Koopman, S. J. (2012). Time Series Analysis by State Space Methods (2nd ed.). Oxford University Press.
4. Foster, G. and Rahmstorf, S. (2011). Global temperature evolution 1979–2010. Environmental Research Letters, 6(4), 044022.
5. Granger, C. W. J. and Newbold, P. (1974). Spurious regressions in econometrics. Journal of Econometrics, 2(2), 111–120.
6. Graven, H. D. et al. (2013). Enhanced seasonal exchange of CO₂ by northern ecosystems since 1960. Science, 341(6150), 1085–1089.
7. Hansen, J. E. et al. (2025). Global warming has accelerated: are the United Nations and the public well-informed? Environment: Science and Policy for Sustainable Development, 67(1), 6–44.
8. Harvey, A. C. (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.
9. Hyndman, R. J. and Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts.
10. Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35–45.
11. Keeling, C. D. (1960). The concentration and isotopic abundances of carbon dioxide in the atmosphere. Tellus, 12(2), 200–203.
12. Lenssen, N. et al. (2019). Improvements in the GISTEMP uncertainty model. Journal of Geophysical Research: Atmospheres, 124(12), 6307–6326. Data via github.com/datasets/global-temp.
13. McElreath, R. (2020). Statistical Rethinking (2nd ed.). CRC Press.
14. Rahmstorf, S., Foster, G. and Cahill, N. (2017). Global temperature evolution: recent trends and some pitfalls. Environmental Research Letters, 12(5), 054001.
15. Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
16. Riutort-Mayol, G., Bürkner, P.-C., Andersen, M. R., Solin, A. and Vehtari, A. (2023). Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing, 33, 17.
17. Yule, G. U. (1926). Why do we sometimes get nonsense-correlations between time-series? Journal of the Royal Statistical Society, 89(1), 1–63.
18. Yule, G. U. (1927). On a method of investigating periodicities in disturbed series, with special reference to Wolfer's sunspot numbers. Philosophical Transactions of the Royal Society A, 226, 267–298.