Forecasting the Keeling Curve

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

September 29, 2026

Let’s explore the Keeling curve a little further, as described in Lecture 7. There we proposed modeling the monthly carbon dioxide concentration at Mauna Loa by \[ y\approx\theta_0+\theta_1t+\theta_2t^2+\theta_3\cos(2\pi t)+\theta_4\sin(2\pi t), \] with $t$ in years since 1958, and asked what concentration we should expect in May 2030. Let’s fit that model with least squares!

The data

NOAA publishes the monthly means as a CSV file, with a header of comments marked by #. The column decimal date is the middle of each month in years (March 1958 is $1958.2027$), and average is the monthly mean concentration in parts per million (ppm).

import datetime

import numpy as np
import plotly.graph_objects as go
import polars as pl
from plotly.subplots import make_subplots
from sklearn.datasets import fetch_file

path = fetch_file("https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.csv")
df = pl.read_csv(path, comment_prefix="#")

year = df["decimal date"].to_numpy()       # calendar time, mid-month
y = df["average"].to_numpy()               # monthly mean CO2, ppm
t = year - 1958                            # years since 1958
n = len(y)

first, last = df.row(0, named=True), df.row(-1, named=True)
print(f"{n} months, from {first['year']}-{first['month']:02d} to {last['year']}-{last['month']:02d}")
print(f"concentration from {y.min():.2f} to {y.max():.2f} ppm")

822 months, from 1958-03 to 2026-08
concentration from 312.42 to 432.34 ppm

Figure style

# The same muted palette as the other examples, which reads correctly against
# both the light and the dark version of this site.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
BLUE, ORANGE, GREEN, PURPLE = "#2a78d6", "#eb6834", "#1baf7a", "#9b5de5"


def style(fig, height=420):
    """Apply the shared layout to a plot with axes."""
    fig.update_layout(
        template="simple_white",
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        font=dict(color=INK, size=13),
        margin=dict(l=70, r=30, t=30, b=55),
        height=height,
        hovermode="x unified",
        legend=dict(orientation="h", y=-0.18, x=0),
    )
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    return fig

Plot the data

fig = go.Figure(go.Scatter(x=year, y=y, mode="lines", name="monthly mean",
                           line=dict(color=INK, width=1.5)))
style(fig).update_layout(showlegend=False)
fig.update_xaxes(title_text="year")
fig.update_yaxes(title_text="CO₂ (ppm)")
fig.show()

Figure 1: The Keeling curve: monthly mean CO₂ at Mauna Loa. Zoom in on any few years to see the seasonal cycle.

Fitting the model

Each month contributes one row $(1,t_i,t_i^2,\cos(2\pi t_i),\sin(2\pi t_i))$ to the design matrix $\bfX\in\Rbb^{n\times5}$. Since $t$ and the calendar year differ by the integer $1958$, the cosine and sine columns are the same whichever we use; only the polynomial columns care. We write the design matrix for a general polynomial degree, since we will want to vary it later.

def design(t, degree=2):
    """Design matrix with columns 1, t, ..., t^degree, cos(2 pi t), sin(2 pi t)."""
    powers = [t ** k for k in range(degree + 1)]
    return np.column_stack(powers + [np.cos(2 * np.pi * t), np.sin(2 * np.pi * t)])


X = design(t)
theta, rss, rank, sing_vals = np.linalg.lstsq(X, y, rcond=None)

print(f"X is {X.shape[0]} x {X.shape[1]} with rank {rank}")
for name, value in zip(["1", "t", "t^2", "cos(2 pi t)", "sin(2 pi t)"], theta):
    print(f"theta for {name:12s} = {value:9.4f}")
print(f"\nroot mean squared residual: {np.sqrt(rss[0] / n):.3f} ppm")

X is 822 x 5 with rank 5
theta for 1            =  314.5689
theta for t            =    0.7419
theta for t^2          =    0.0135
theta for cos(2 pi t)  =   -0.9887
theta for sin(2 pi t)  =    2.6589

root mean squared residual: 0.991 ppm

The matrix has full column rank, so by the set of least squares solutions the coefficients are unique, and we can read them off.

The trend $\theta_0+\theta_1t+\theta_2t^2$ has growth rate $\theta_1+2\theta_2t$ ppm per year.
The seasonal part can be rewritten as a single shifted cosine, $\theta_3\cos(2\pi t)+\theta_4\sin(2\pi t)=A\cos\big(2\pi(t-\tau)\big)$, with amplitude $A=\sqrt{\theta_3^2+\theta_4^2}$ and peak at the fraction $\tau$ of the year where $(\cos 2\pi\tau,\sin 2\pi\tau)=(\theta_3,\theta_4)/A$.

for when in [1958, 1990, 2026]:
    print(f"growth rate in {when}: {theta[1] + 2 * theta[2] * (when - 1958):.2f} ppm/year")

A = np.hypot(theta[3], theta[4])
tau = (np.arctan2(theta[4], theta[3]) / (2 * np.pi)) % 1
peak = datetime.date(2001, 1, 1) + datetime.timedelta(days=tau * 365.25)
print(f"\nseasonal amplitude {A:.2f} ppm (peak to trough {2 * A:.2f} ppm), peaking around {peak:%B %d}")

growth rate in 1958: 0.74 ppm/year
growth rate in 1990: 1.60 ppm/year
growth rate in 2026: 2.58 ppm/year

seasonal amplitude 2.84 ppm (peak to trough 5.67 ppm), peaking around April 23

So the rate of increase has more than tripled since 1958, and every year the concentration swings by about $5.7$ ppm, highest in spring and lowest in autumn. (The Northern Hemisphere’s plants take up carbon dioxide over the summer growing season and release it again over the winter.)

Figure 2 shows the fitted values $\widehat{\bfy}=\bfX\widehat{\bftheta}$ together with the residuals $\widehat{\bfr}=\bfy-\widehat{\bfy}$.

Plot the fit and the residuals

t_may_2030 = 2030 + 4.5 / 12 - 1958                   # mid-May, as in the data
y_may_2030 = (design(np.array([t_may_2030])) @ theta)[0]
year_grid = np.linspace(year[0], 2031, 3000)
fit_grid = design(year_grid - 1958) @ theta
residuals = y - X @ theta

fig = make_subplots(rows=2, cols=1, shared_xaxes=True, row_heights=[0.68, 0.32],
                    vertical_spacing=0.05)
fig.add_trace(go.Scatter(x=year, y=y, mode="lines", name="data",
                         line=dict(color=INK, width=1.5)), row=1, col=1)
fig.add_trace(go.Scatter(x=year_grid, y=fit_grid, mode="lines", name="least squares fit",
                         line=dict(color=BLUE, width=1.5)), row=1, col=1)
fig.add_trace(go.Scatter(x=[2030 + 4.5 / 12], y=[y_may_2030], mode="markers",
                         name="May 2030 forecast",
                         marker=dict(color=ORANGE, size=10, line=dict(color="white", width=1))),
              row=1, col=1)
fig.add_trace(go.Scatter(x=year, y=residuals, mode="lines", name="residual",
                         showlegend=False, line=dict(color=INK, width=1)), row=2, col=1)
fig.add_hline(y=0, line=dict(color=BLUE, width=1), row=2, col=1)
style(fig, height=600)
fig.update_xaxes(title_text="year", row=2, col=1)
fig.update_yaxes(title_text="CO₂ (ppm)", row=1, col=1)
fig.update_yaxes(title_text="residual (ppm)", row=2, col=1)
fig.show()

Figure 2: Top: the data and the least squares fit, extended to the end of 2030, with the May 2030 forecast marked. Bottom: the residuals.

Five numbers describe nearly seven decades of data to within about a part per million! The residuals do roughly drift as time goes on. I’m told that this is partly due to El Nino, but I’m not a domain expert and one should avoid overly interpreting a model whose data they don’t understand.

Measure time from 1958, not from year 0

Using the calendar year instead of time from 1958 can actually cause some numerical problems. With the calendar year itself, the columns $1$, $t$, and $t^2$ are nearly parallel (every $t_i$ is between $1958$ and $2027$), and the design matrix is badly ill-conditioned.

for label, time in [("calendar year", year), ("years since 1958", t)]:
    Z = design(time)
    print(f"{label:17s} kappa(X) = {np.linalg.cond(Z):.1e}   kappa(X^T X) = {np.linalg.cond(Z.T @ Z):.1e}")

calendar year     kappa(X) = 4.5e+10   kappa(X^T X) = 1.7e+21
years since 1958  kappa(X) = 6.4e+03   kappa(X^T X) = 4.1e+07

With calendar years, $\kappa(\bfX^\top\bfX)$ is far beyond $10^{16}$, so solving the normal equations can lose every significant digit we have. The fitted curve is the same either way in exact arithmetic, since both design matrices have the same column space.

Is it any good at forecasting?

A small residual on the data we fit is not the same thing as a good forecast. To check, pretend that it is the start of 2006: fit on the data before 2006 only, forecast the twenty years since, and compare against what actually happened. To see why the lecture chose a quadratic trend, we try every polynomial degree from $1$ to $5$ (keeping the same seasonal columns).

def rmse(r):
    return np.sqrt(np.mean(r ** 2))


cutoff = 2006
train = year < cutoff
test = ~train
forecasts = {}

print(f"{'degree':>6s}{'train RMSE':>12s}{'test RMSE':>12s}")
for degree in range(1, 6):
    Z = design(t, degree)
    theta_train, *_ = np.linalg.lstsq(Z[train], y[train], rcond=None)
    forecasts[degree] = design(year_grid - 1958, degree) @ theta_train
    print(f"{degree:6d}{rmse(y[train] - Z[train] @ theta_train):12.2f}"
          f"{rmse(y[test] - Z[test] @ theta_train):12.2f}")

degree  train RMSE   test RMSE
     1        2.23       14.88
     2        0.95        2.15
     3        0.88        8.58
     4        0.82        9.19
     5        0.76       60.09

Plot the backtest

fig = go.Figure()
fig.add_trace(go.Scatter(x=year, y=y, mode="lines", name="data",
                         line=dict(color=INK, width=1.5)))
for degree, color in [(1, ORANGE), (2, BLUE), (3, GREEN), (4, PURPLE)]:
    fig.add_trace(go.Scatter(x=year_grid, y=forecasts[degree], mode="lines",
                             name=f"degree {degree}",
                             line=dict(color=color, width=2.5 if degree == 2 else 1.5)))
fig.add_vline(x=cutoff, line=dict(color=INK, width=1, dash="dash"))
style(fig, height=480)
fig.update_xaxes(title_text="year", range=[1990, 2031])
fig.update_yaxes(title_text="CO₂ (ppm)", range=[345, 450])
fig.show()

Figure 3: Fitting on the data before 2006 (left of the dashed line) and forecasting forward. Higher-degree trends fit the past slightly better and the future much worse. Note, I cut the degree 4 curve off because it starts exploding off into the distance.

This is a central tension of prediction. Every extra column enlarges the column space, so the training error can only go down with the degree. That is, least squares is projecting onto a bigger subspace. But a high-degree polynomial spends its extra freedom bending to fit wiggles in the past, and outside the data those bends run off to infinity. This can be referred to as a generalization gap. A straight line is too rigid in the other direction and cannot keep up with the acceleration.

Warning

Beware! I think we’ll have time to discuss this later, but modern models often do not need to be carefully designed to balance generalization and model size. They just need to be really fucking big. This is visible through the double descent phenomena, see https://en.wikipedia.org/wiki/Double_descent.

--- title: "Forecasting the Keeling Curve" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-09-29 number-sections: false jupyter: python3 format: html: toc: true toc-depth: 2 code-line-numbers: true code-tools: true out-width: "100%" column: page execute: echo: true draft: false --- Let's explore the Keeling curve a little further, as described in [Lecture 7](../../lectures/lecture-07.qmd#linear-regression). There we proposed modeling the monthly carbon dioxide concentration at Mauna Loa by $$ y\approx\theta_0+\theta_1t+\theta_2t^2+\theta_3\cos(2\pi t)+\theta_4\sin(2\pi t), $$ with $t$ in years since 1958, and asked what concentration we should expect in May 2030. Let's fit that model with least squares! ## The data NOAA publishes the monthly means as a CSV file, with a header of comments marked by `#`. The column `decimal date` is the middle of each month in years (March 1958 is $1958.2027$), and `average` is the monthly mean concentration in parts per million (ppm). ```{python} import datetime import numpy as np import plotly.graph_objects as go import polars as pl from plotly.subplots import make_subplots from sklearn.datasets import fetch_file path = fetch_file("https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.csv") df = pl.read_csv(path, comment_prefix="#") year = df["decimal date"].to_numpy() # calendar time, mid-month y = df["average"].to_numpy() # monthly mean CO2, ppm t = year - 1958 # years since 1958 n = len(y) first, last = df.row(0, named=True), df.row(-1, named=True) print(f"{n} months, from {first['year']}-{first['month']:02d} to {last['year']}-{last['month']:02d}") print(f"concentration from {y.min():.2f} to {y.max():.2f} ppm") ``` ```{python} #| code-fold: true #| code-summary: "Figure style" # The same muted palette as the other examples, which reads correctly against # both the light and the dark version of this site. INK = "#7c7c78" GRID = "rgba(128, 128, 128, 0.28)" BLUE, ORANGE, GREEN, PURPLE = "#2a78d6", "#eb6834", "#1baf7a", "#9b5de5" def style(fig, height=420): """Apply the shared layout to a plot with axes.""" fig.update_layout( template="simple_white", paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)", font=dict(color=INK, size=13), margin=dict(l=70, r=30, t=30, b=55), height=height, hovermode="x unified", legend=dict(orientation="h", y=-0.18, x=0), ) fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) return fig ``` ```{python} #| label: fig-keeling #| fig-cap: "The Keeling curve: monthly mean CO₂ at Mauna Loa. Zoom in on any few years to see the seasonal cycle." #| fig-alt: "A line plot of monthly CO2 concentration from 1958 to 2026, rising from about 315 ppm to about 430 ppm, getting steeper over time, with a small yearly oscillation." #| code-fold: true #| code-summary: "Plot the data" fig = go.Figure(go.Scatter(x=year, y=y, mode="lines", name="monthly mean", line=dict(color=INK, width=1.5))) style(fig).update_layout(showlegend=False) fig.update_xaxes(title_text="year") fig.update_yaxes(title_text="CO₂ (ppm)") fig.show() ``` ## Fitting the model Each month contributes one row $(1,t_i,t_i^2,\cos(2\pi t_i),\sin(2\pi t_i))$ to the design matrix $\bfX\in\Rbb^{n\times5}$. Since $t$ and the calendar year differ by the integer $1958$, the cosine and sine columns are the same whichever we use; only the polynomial columns care. We write the design matrix for a general polynomial degree, since we will want to vary it later. ```{python} def design(t, degree=2): """Design matrix with columns 1, t, ..., t^degree, cos(2 pi t), sin(2 pi t).""" powers = [t ** k for k in range(degree + 1)] return np.column_stack(powers + [np.cos(2 * np.pi * t), np.sin(2 * np.pi * t)]) X = design(t) theta, rss, rank, sing_vals = np.linalg.lstsq(X, y, rcond=None) print(f"X is {X.shape[0]} x {X.shape[1]} with rank {rank}") for name, value in zip(["1", "t", "t^2", "cos(2 pi t)", "sin(2 pi t)"], theta): print(f"theta for {name:12s} = {value:9.4f}") print(f"\nroot mean squared residual: {np.sqrt(rss[0] / n):.3f} ppm") ``` The matrix has full column rank, so by [the set of least squares solutions](../../lectures/lecture-07.qmd#thm-least-squares-uniqueness) the coefficients are unique, and we can read them off. - The *trend* $\theta_0+\theta_1t+\theta_2t^2$ has growth rate $\theta_1+2\theta_2t$ ppm per year. - The *seasonal part* can be rewritten as a single shifted cosine, $\theta_3\cos(2\pi t)+\theta_4\sin(2\pi t)=A\cos\big(2\pi(t-\tau)\big)$, with amplitude $A=\sqrt{\theta_3^2+\theta_4^2}$ and peak at the fraction $\tau$ of the year where $(\cos 2\pi\tau,\sin 2\pi\tau)=(\theta_3,\theta_4)/A$. ```{python} for when in [1958, 1990, 2026]: print(f"growth rate in {when}: {theta[1] + 2 * theta[2] * (when - 1958):.2f} ppm/year") A = np.hypot(theta[3], theta[4]) tau = (np.arctan2(theta[4], theta[3]) / (2 * np.pi)) % 1 peak = datetime.date(2001, 1, 1) + datetime.timedelta(days=tau * 365.25) print(f"\nseasonal amplitude {A:.2f} ppm (peak to trough {2 * A:.2f} ppm), peaking around {peak:%B %d}") ``` So the rate of increase has more than tripled since 1958, and every year the concentration swings by about $5.7$ ppm, highest in spring and lowest in autumn. (The Northern Hemisphere's plants take up carbon dioxide over the summer growing season and release it again over the winter.) @fig-fit shows the fitted values $\widehat{\bfy}=\bfX\widehat{\bftheta}$ together with the residuals $\widehat{\bfr}=\bfy-\widehat{\bfy}$. ```{python} #| label: fig-fit #| fig-cap: "Top: the data and the least squares fit, extended to the end of 2030, with the May 2030 forecast marked. Bottom: the residuals." #| fig-alt: "Top panel: the data with the fitted curve on top, nearly indistinguishable, and the fitted curve continuing past 2026 to 2031 with a marked point at May 2030. Bottom panel: residuals within about plus or minus 3 ppm, with slow multi-year waves rather than pure noise." #| code-fold: true #| code-summary: "Plot the fit and the residuals" t_may_2030 = 2030 + 4.5 / 12 - 1958 # mid-May, as in the data y_may_2030 = (design(np.array([t_may_2030])) @ theta)[0] year_grid = np.linspace(year[0], 2031, 3000) fit_grid = design(year_grid - 1958) @ theta residuals = y - X @ theta fig = make_subplots(rows=2, cols=1, shared_xaxes=True, row_heights=[0.68, 0.32], vertical_spacing=0.05) fig.add_trace(go.Scatter(x=year, y=y, mode="lines", name="data", line=dict(color=INK, width=1.5)), row=1, col=1) fig.add_trace(go.Scatter(x=year_grid, y=fit_grid, mode="lines", name="least squares fit", line=dict(color=BLUE, width=1.5)), row=1, col=1) fig.add_trace(go.Scatter(x=[2030 + 4.5 / 12], y=[y_may_2030], mode="markers", name="May 2030 forecast", marker=dict(color=ORANGE, size=10, line=dict(color="white", width=1))), row=1, col=1) fig.add_trace(go.Scatter(x=year, y=residuals, mode="lines", name="residual", showlegend=False, line=dict(color=INK, width=1)), row=2, col=1) fig.add_hline(y=0, line=dict(color=BLUE, width=1), row=2, col=1) style(fig, height=600) fig.update_xaxes(title_text="year", row=2, col=1) fig.update_yaxes(title_text="CO₂ (ppm)", row=1, col=1) fig.update_yaxes(title_text="residual (ppm)", row=2, col=1) fig.show() ``` Five numbers describe nearly seven decades of data to within about a part per million! The residuals do roughly drift as time goes on. I'm told that this is partly due to El Nino, but I'm not a domain expert and one should avoid overly interpreting a model whose data they don't understand. ::: {.callout-warning} ## Measure time from 1958, not from year 0 Using the calendar year instead of time from 1958 can actually cause some numerical problems. With the calendar year itself, the columns $1$, $t$, and $t^2$ are nearly parallel (every $t_i$ is between $1958$ and $2027$), and the design matrix is badly [ill-conditioned](../../lectures/lecture-07.qmd#the-full-rank-case-the-hat-matrix). ```{python} for label, time in [("calendar year", year), ("years since 1958", t)]: Z = design(time) print(f"{label:17s} kappa(X) = {np.linalg.cond(Z):.1e} kappa(X^T X) = {np.linalg.cond(Z.T @ Z):.1e}") ``` With calendar years, $\kappa(\bfX^\top\bfX)$ is far beyond $10^{16}$, so solving the normal equations can lose every significant digit we have. The fitted curve is the same either way in exact arithmetic, since both design matrices have the same column space. ::: ## Is it any good at forecasting? A small residual on the data we fit is not the same thing as a good forecast. To check, pretend that it is the start of 2006: fit on the data before 2006 only, forecast the twenty years since, and compare against what actually happened. To see why the lecture chose a *quadratic* trend, we try every polynomial degree from $1$ to $5$ (keeping the same seasonal columns). ```{python} def rmse(r): return np.sqrt(np.mean(r ** 2)) cutoff = 2006 train = year < cutoff test = ~train forecasts = {} print(f"{'degree':>6s}{'train RMSE':>12s}{'test RMSE':>12s}") for degree in range(1, 6): Z = design(t, degree) theta_train, *_ = np.linalg.lstsq(Z[train], y[train], rcond=None) forecasts[degree] = design(year_grid - 1958, degree) @ theta_train print(f"{degree:6d}{rmse(y[train] - Z[train] @ theta_train):12.2f}" f"{rmse(y[test] - Z[test] @ theta_train):12.2f}") ``` ```{python} #| label: fig-backtest #| fig-cap: "Fitting on the data before 2006 (left of the dashed line) and forecasting forward. Higher-degree trends fit the past slightly better and the future much worse. Note, I cut the degree 4 curve off because it starts exploding off into the distance." #| fig-alt: "The data with four forecast curves fitted on data before 2006. After 2006, the degree 1 curve falls well below the data, the degree 2 curve tracks the data closely, and the degree 3 and 4 curves bend away from the data." #| code-fold: true #| code-summary: "Plot the backtest" fig = go.Figure() fig.add_trace(go.Scatter(x=year, y=y, mode="lines", name="data", line=dict(color=INK, width=1.5))) for degree, color in [(1, ORANGE), (2, BLUE), (3, GREEN), (4, PURPLE)]: fig.add_trace(go.Scatter(x=year_grid, y=forecasts[degree], mode="lines", name=f"degree {degree}", line=dict(color=color, width=2.5 if degree == 2 else 1.5))) fig.add_vline(x=cutoff, line=dict(color=INK, width=1, dash="dash")) style(fig, height=480) fig.update_xaxes(title_text="year", range=[1990, 2031]) fig.update_yaxes(title_text="CO₂ (ppm)", range=[345, 450]) fig.show() ``` This is a central tension of prediction. Every extra column enlarges the column space, so the training error can only go down with the degree. That is, least squares is projecting onto a bigger subspace. But a high-degree polynomial spends its extra freedom bending to fit wiggles in the past, and outside the data those bends run off to infinity. This can be referred to as a generalization gap. A straight line is too rigid in the other direction and cannot keep up with the acceleration. :::{.callout-warning} Beware! I think we'll have time to discuss this later, but modern models often do not need to be carefully designed to balance generalization and model size. They just need to be *really fucking big*. This is visible through the double descent phenomena, see <https://en.wikipedia.org/wiki/Double_descent>. :::