1.1 The independence problem

Imagine you measure air temperature every hour for a week. You have 168 observations. If you run a t-test comparing morning to afternoon temperatures, the test assumes you have 168 independent pieces of information. But 2 PM today is very similar to 2 PM yesterday — those readings are not independent. The true number of independent observations is much smaller. Ignoring autocorrelation makes your results look far more statistically significant than they really are.

Time series analysis is the branch of statistics that explicitly models this temporal dependence rather than treating it as a nuisance to ignore.

1.2 The four structural components

Every environmental time series can be decomposed into at most four overlapping signals:

Understanding which components are present in your data determines which analysis tools are appropriate.

1.3 Stationarity — the hidden requirement

Most time series models (including ARIMA) require the series to be stationary — meaning its statistical properties (mean, variance, autocorrelation) do not change over time. A series with a trend is non-stationary because the mean is changing. Seasonality is also a form of non-stationarity.

The solution is differencing: subtracting today’s value from yesterday’s. One differencing step removes a linear trend. One seasonal differencing step (subtracting the value 12 months ago) removes annual seasonality. ARIMA models formalise this process.

1.4 The analytical workflow

Time series analysis follows a consistent sequence:

1. Plot the raw series — identify trend, seasonality, obvious outliers
2. Decompose — separate trend, seasonal, and remainder components
3. Check ACF/PACF — diagnose the autocorrelation structure
4. Difference if needed — achieve stationarity
5. Fit a model — ARIMA or simpler alternatives
6. Diagnose residuals — confirm they are white noise (random, no structure)
7. Forecast — generate predictions with calibrated uncertainty intervals

This module follows that sequence exactly for two real datasets.

1.5 Tools used in this module

# Run once in the Console if not yet installed
install.packages(c("tsibble", "forecast"))

2.1 Packages

library(tidyverse)
library(lubridate)
library(tsibble)
# forecast loaded conditionally in the forecasting section

2.2 Why tsibble instead of a regular data frame?

A regular tibble or data.frame has no concept of time. You could sort the rows any way you like, and R would not object. A tsibble enforces three rules:

1. One row per time point — duplicate timestamps cause an error
2. No implicit gaps — fill_gaps() makes missing time steps explicit rather than silently absent
3. Explicit index — the time column is locked in; you can’t accidentally overwrite it

This matters because almost every time-series error in R comes from data that looks like a time series but is actually just a data frame with a date column.

2.3 Handling live data downloads

Both datasets below are downloaded from public URLs at NOAA and NASA. The download is wrapped in tryCatch() — if the server is unreachable, a fallback dataset is used so the document still knits. In practice, once you have confirmed the download works, save the file locally with write_csv(co2_raw, "data/co2_noaa.csv") and read from disk instead.

co2_url <- "https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.csv"

co2_raw <- tryCatch(
  read_csv(
    co2_url,
    comment   = "#",
    col_names = c("year", "month", "decimal_date",
                  "average", "interpolated", "trend", "days"),
    show_col_types = FALSE
  ),
  error = function(e) {
    message("Download failed — using built-in co2 dataset as fallback.")
    NULL
  }
)

# Fallback: convert R's built-in co2 time-series object
if (is.null(co2_raw)) {
  co2_raw <- tibble(
    year  = rep(1959:1997, each = 12)[seq_along(co2)],
    month = rep(1:12, length.out = length(co2)),
    trend = as.numeric(co2),
    average = as.numeric(co2),
    interpolated = as.numeric(co2)
  )
}

glimpse(co2_raw)

## Rows: 824
## Columns: 8
## $ year         <chr> "year", "1958", "1958", "1958", "1958", "1958", "1958", "…
## $ month        <chr> "month", "3", "4", "5", "6", "7", "8", "9", "10", "11", "…
## $ decimal_date <chr> "decimal date", "1958.2027", "1958.2877", "1958.3699", "1…
## $ average      <chr> "average", "315.71", "317.45", "317.51", "317.27", "315.8…
## $ interpolated <chr> "deseasonalized", "314.44", "315.16", "314.69", "315.15",…
## $ trend        <chr> "ndays", "-1", "-1", "-1", "-1", "-1", "-1", "-1", "-1", …
## $ days         <chr> "sdev", "-9.99", "-9.99", "-9.99", "-9.99", "-9.99", "-9.…
## $ X8           <chr> "unc", "-0.99", "-0.99", "-0.99", "-0.99", "-0.99", "-0.9…

The NOAA file has three CO2 columns worth knowing:

• average — the raw monthly mean; a few months marked -9.99 for missing
• interpolated — gaps filled using nearby values; complete from 1958 onward
• trend — the long-term trend after the seasonal cycle has been removed

3.1 Background: why this dataset matters

Charles David Keeling began measuring atmospheric CO2 at Mauna Loa Observatory, Hawaii, in March 1958. He chose Mauna Loa deliberately: it sits at 3,400 m elevation, far from local pollution sources, and samples air that has mixed across the entire Pacific Basin. The resulting record — now maintained by NOAA after Keeling’s death — is the longest continuous atmospheric CO2 record in existence and the primary observational evidence for the anthropogenic greenhouse effect.

Two features of the record are immediately visible to the eye:

1. An upward trend — CO2 has risen from ~315 ppm in 1958 to over 420 ppm today, a ~34% increase driven overwhelmingly by fossil fuel combustion
2. A seasonal sawtooth — CO2 dips every Northern Hemisphere summer as land vegetation absorbs carbon through photosynthesis, then rises each winter as vegetation goes dormant and decomposition continues

3.2 Build the tsibble

co2_ts <- co2_raw |>
  filter(!is.na(year), !is.na(month), average > 0) |>
  mutate(date = make_yearmonth(year = as.integer(year),
                               month = as.integer(month))) |>
  drop_na(date) |>
  as_tsibble(index = date) |>
  fill_gaps() |>
  tidyr::fill(average, trend, interpolated, decimal_date,
              .direction = "downup")

co2_ts <- co2_ts |>
  filter(as.integer(year) >= 1958)   # Mauna Loa record starts March 1958

n_months_co2 <- nrow(co2_ts)
yr_start     <- min(as.integer(co2_ts$year))
yr_end       <- max(as.integer(co2_ts$year))

The dataset runs from 1958 to 2026 — 823 months of continuous monitoring.

make_yearmonth() creates a yearmonth index — a tsibble-native time type that represents a calendar month without specifying a day. This is appropriate here because each measurement is a monthly mean, not a point-in-time reading.

3.3 Plot: the long-term trend

co2_ts |>
  as_tibble() |>
  mutate(trend   = as.numeric(trend),
         date_dt = as.Date(date)) |>
  ggplot(aes(date_dt, trend)) +
  geom_line(color = ts_red, linewidth = 0.8) +
  geom_smooth(method = "loess", span = 0.15,
              color = ts_green, linewidth = 0.6,
              linetype = "dashed", se = FALSE) +
  labs(title    = "Keeling Curve — Long-Term CO2 Trend",
       subtitle = paste("Mauna Loa Observatory ·", yr_start, "–", yr_end,
                        "· seasonal cycle removed"),
       x = NULL, y = "CO2 concentration (ppm)")

What to look for: The trend is clearly non-linear — the rate of increase has accelerated over the record. The curve is concave upward, meaning each decade adds more CO2 than the previous one. This non-linearity is important: a simple linear regression would underestimate recent rates of increase.

3.4 Plot: raw monthly values showing seasonality

co2_ts |>
  as_tibble() |>
  mutate(average = as.numeric(average),
         trend   = as.numeric(trend),
         date_dt = as.Date(date)) |>
  ggplot(aes(date_dt, average)) +
  geom_line(color = ts_amber, linewidth = 0.5, alpha = 0.8) +
  geom_line(aes(y = trend), color = ts_red, linewidth = 0.9) +
  labs(title    = "Monthly CO2 Concentration vs De-seasonalised Trend",
       subtitle = "Amber = monthly average with seasonal cycle · Red = trend component only",
       x = NULL, y = "CO2 (ppm)")

What to look for: The difference between the amber (raw) and red (trend) lines is the seasonal cycle — roughly ±3–4 ppm around the trend. Notice how regular it is: the same pattern repeats every year, with CO2 peaking around May and reaching its minimum around September–October. This regularity is exactly what makes it a good candidate for STL decomposition.

3.5 Recent rate of change

recent_decade <- co2_ts |>
  filter(as.integer(year) >= max(as.integer(year)) - 9) |>
  mutate(trend_num   = as.numeric(trend),
         dec_date_num = as.numeric(decimal_date))

slope <- coef(lm(trend_num ~ dec_date_num, data = recent_decade))["dec_date_num"]

Over the most recent decade in the record, CO2 has been rising at -0.69 ppm per year — faster than at any previous point in the record. For context, over the entire record the average rate was approximately 1.6 ppm per year.

4.1 What is decomposition and why do it?

Before fitting any model, you need to understand what kind of patterns are in your data. Decomposition answers that question by separating a time series into its structural parts:

\[y_t = T_t + S_t + R_t\]

• \(T_t\) — the trend: the long-run direction after smoothing out short-term fluctuations
• \(S_t\) — the seasonal component: the repeating calendar-driven pattern
• \(R_t\) — the remainder: everything the trend and seasonal model cannot explain — weather events, measurement error, unusual years

Knowing the relative sizes of these components tells you what matters most in your data, and whether a simple trend model is sufficient or whether seasonality needs to be explicitly modelled.

4.2 STL decomposition — the method

STL stands for Seasonal-Trend decomposition using LOESS (locally weighted scatterplot smoothing). It was developed by Cleveland et al. (1990) and remains the gold standard for robust seasonal decomposition because:

• It handles any length of seasonality (monthly, daily, hourly)
• It is robust — outliers and unusual events in the remainder do not distort the trend or seasonal estimates
• The seasonal component can change slowly over time (unlike classical decomposition which assumes a fixed seasonal pattern)

The robust = TRUE argument in stl() down-weights outlier months so they don’t pull the trend or seasonal estimates. Always use it for environmental data where extreme events (wildfires, El Niño) can create outlier years.

4.3 STL decomposition — CO2

# Build a base-R ts object — stl() is rock-solid and has no tsibble dependency
co2_vec <- co2_ts |>
  as_tibble() |>
  arrange(date) |>
  mutate(interpolated = as.numeric(interpolated)) |>
  pull(interpolated)

co2_base_ts <- ts(co2_vec,
                  start     = c(yr_start, 1),
                  frequency = 12)

co2_stl <- stl(co2_base_ts, s.window = "periodic", robust = TRUE)
plot(co2_stl, main = "STL Decomposition — Mauna Loa CO2")

4.4 Seasonal pattern in detail

# Extract components and plot the seasonal cycle manually
co2_comp <- as_tibble(co2_ts) |>
  arrange(date) |>
  mutate(
    interpolated = as.numeric(interpolated),
    seasonal  = co2_stl$time.series[, "seasonal"],
    trend_cmp = co2_stl$time.series[, "trend"],
    remainder = co2_stl$time.series[, "remainder"],
    month_num = lubridate::month(as.Date(date))
  )

co2_comp |>
  group_by(month_num) |>
  summarise(mean_seasonal = mean(seasonal), .groups = "drop") |>
  mutate(month_label = month.abb[month_num]) |>
  ggplot(aes(x = month_num, y = mean_seasonal)) +
  geom_col(fill = ts_sage, color = ts_green, width = 0.7) +
  geom_hline(yintercept = 0, color = ts_red, linetype = "dashed") +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  labs(title    = "Mean Seasonal Component — Mauna Loa CO2",
       subtitle = "Positive = CO2 above trend that month · Negative = below trend",
       x = "Month", y = "Seasonal effect (ppm)")

Interpretation: CO2 is highest above trend in May (Northern Hemisphere vegetation still dormant or just leafing out) and lowest in September–October (peak photosynthesis has drawn down CO2 all summer). The range of the seasonal cycle (~7 ppm peak-to-trough) is the net effect of Northern Hemisphere vegetation “breathing” — the Northern Hemisphere dominates because it has far more land and forest than the Southern.

4.5 Monthly mean across all years

co2_comp |>
  group_by(month_num) |>
  summarise(
    mean_co2 = mean(interpolated),
    se_co2   = sd(interpolated) / sqrt(n()),
    .groups = "drop"
  ) |>
  ggplot(aes(x = month_num, y = mean_co2,
             ymin = mean_co2 - se_co2, ymax = mean_co2 + se_co2)) +
  geom_pointrange(color = ts_red, size = 0.7) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  labs(title    = "Monthly Mean CO2 Across All Years ± 1 SE",
       subtitle = "Averages the raw value for each calendar month across the full record",
       x = "Month", y = "CO2 (ppm)")

Note: The error bars here are very small — not because the months don’t vary year to year, but because we have 69 years of data. The SE is ±SD/√n. With n ≈ 65 years, even moderate year-to-year variation produces tiny SEs. The wide range of mean values across months (the seasonal signal) overwhelms the year-to-year variability within each month.

5.1 What is autocorrelation?

Autocorrelation is correlation of a variable with a delayed version of itself. If this month’s CO2 is high, next month’s tends to also be high — that is positive autocorrelation at lag 1. If a warm summer tends to be followed by a cooler-than-average autumn, that is negative autocorrelation at lag 1.

Autocorrelation is measured at each lag \(k\) — the number of time steps between the two observations being compared:

\[r_k = \text{Cor}(y_t, \; y_{t-k})\]

5.2 The autocorrelation function (ACF)

The ACF plot shows \(r_k\) for \(k = 1, 2, 3, \ldots\) lags. Each vertical bar is one lag correlation. The horizontal dashed blue lines mark the approximate 95% significance threshold: \(\pm 1.96 / \sqrt{n}\).

How to read an ACF plot:

• Bars extending beyond the blue lines → significant autocorrelation at that lag
• A pattern that decays slowly across many lags → the series is non-stationary (has a trend or other long-range dependence); differencing is needed before fitting an ARIMA model
• A pattern with sharp peaks at multiples of 12 (for monthly data) → seasonal autocorrelation that needs to be modelled
• Bars that fall within the blue lines for all lags after lag 1–2 → residuals look like white noise — the model has done its job

5.3 The partial autocorrelation function (PACF)

The PACF measures the correlation between \(y_t\) and \(y_{t-k}\) after removing the effect of all intervening lags (\(y_{t-1}, y_{t-2}, \ldots, y_{t-k+1}\)). It isolates the direct relationship at each lag.

Why PACF matters for model building:

• If the PACF cuts off sharply after lag \(p\) (goes to near zero) → the data has an AR(\(p\)) structure — past values up to lag \(p\) are directly predictive
• If the ACF cuts off sharply after lag \(q\) → the data has an MA(\(q\)) structure — past forecast errors up to lag \(q\) enter the model
• If both decay gradually → an ARMA(\(p\),\(q\)) model with both components is needed

5.4 ACF and PACF of the CO2 remainder

We examine the remainder from the STL decomposition — not the raw series. After removing trend and seasonality, the remainder should ideally be white noise. The ACF/PACF reveal whether any autocorrelation structure is left over.

remainder_vec <- co2_stl$time.series[, "remainder"]

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
acf(remainder_vec,  lag.max = 36,
    main = "ACF — CO2 Remainder",
    col  = ts_blue, lwd = 1.5)
pacf(remainder_vec, lag.max = 36,
     main = "PACF — CO2 Remainder",
     col  = ts_red,  lwd = 1.5)

par(mfrow = c(1, 1))

Interpretation: If the ACF shows significant spikes at lags 1 and 2 (with the PACF cutting off after lag 1), that is consistent with an AR(1) or AR(2) structure in the remainder — meaning month-to-month CO2 deviations from the trend-seasonal model are somewhat persistent. El Niño years, for instance, can produce slightly elevated CO2 for several consecutive months as tropical fire emissions increase.

5.5 Ljung-Box test — is there any significant autocorrelation?

The ACF plot is visual. The Ljung-Box test provides a formal hypothesis test: H₀ = the first \(m\) autocorrelations are jointly zero (white noise); H₁ = at least one is non-zero (significant structure remains).

A significant p-value (p < 0.05) means the remainder is not white noise — the decomposition has not fully captured the temporal patterns, and a richer model (e.g., ARIMA on the remainder) is warranted.

lb <- Box.test(remainder_vec, lag = 24, type = "Ljung-Box")

tibble(
  Statistic  = round(lb$statistic, 4),
  df         = lb$parameter,
  p_value    = round(lb$p.value,   4),
  Conclusion = ifelse(lb$p.value < 0.05,
                      "Autocorrelation remains — richer model warranted",
                      "Residuals consistent with white noise")
) |>
  kable(caption = "Ljung-Box test: are STL remainder autocorrelations jointly zero?") |>
  kable_styling(full_width = FALSE)

Rule of thumb for lag choice: Use lag = min(10, n/5) for non-seasonal data, or lag = 2×m for seasonal data where m is the seasonal period. Here, lag = 24 (2 × 12 months) tests whether any of the first two years of autocorrelations are collectively non-zero.

6.1 Background: what is a temperature anomaly?

NASA GISTEMP reports temperature as an anomaly — the deviation from a baseline average — rather than as an absolute temperature. The baseline here is the 1951–1980 mean for each calendar month at each location. A value of +0.8°C means that month was 0.8°C warmer than the same month averaged over 1951–1980.

Anomalies are used instead of absolute temperatures because:

1. Different weather stations have very different average temperatures (Alaska vs California), making direct comparison meaningless
2. Stations come and go over time, which would cause artificial jumps in a record based on absolute values
3. Climate science cares about change relative to a reference period, not absolute temperature

A global mean anomaly is the area-weighted average of anomalies from thousands of land stations and ocean buoys. When you see “global temperature rose 1.1°C since pre-industrial times,” that is this kind of anomaly.

6.2 Download and parse

temp_url <- "https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv"

temp_raw <- tryCatch(
  read_csv(temp_url, skip = 1, show_col_types = FALSE,
           na = c("", "NA", "***")),
  error = function(e) {
    message("NASA download failed — generating synthetic temperature series.")
    NULL
  }
)

# Fallback: simulate a warming series
if (is.null(temp_raw)) {
  set.seed(42)
  temp_raw <- tibble(
    Year = rep(1880:2023, each = 12),
    Jan = NA, Feb = NA, Mar = NA, Apr = NA,
    May = NA, Jun = NA, Jul = NA, Aug = NA,
    Sep = NA, Oct = NA, Nov = NA, Dec = NA
  )
  n <- nrow(temp_raw)
  trend_vals <- seq(-0.5, 1.2, length.out = n)
  for (col in month.abb) {
    temp_raw[[col]] <- trend_vals + rnorm(n, 0, 0.2)
  }
}

temp <- temp_raw |>
  select(Year, any_of(month.abb)) |>
  pivot_longer(cols = any_of(month.abb),
               names_to  = "month_abb",
               values_to = "anomaly") |>
  filter(!is.na(anomaly), is.numeric(anomaly)) |>
  mutate(
    month_num = match(month_abb, month.abb),
    date      = make_yearmonth(year = as.integer(Year),
                               month = as.integer(month_num))
  ) |>
  drop_na(date) |>
  as_tsibble(index = date) |>
  fill_gaps() |>
  tidyr::fill(anomaly, .direction = "downup")

yr_temp_start <- min(temp$Year)
yr_temp_end   <- max(temp$Year)

The temperature record spans 1880 to 2026 — 1761 monthly observations across 147 years.

6.3 Time series plot

temp |>
  as_tibble() |>
  mutate(anomaly = as.numeric(anomaly),
         date_dt = as.Date(date)) |>
  ggplot(aes(date_dt, anomaly)) +
  geom_line(color = "grey70", linewidth = 0.4) +
  geom_smooth(method = "loess", span = 0.06,
              color = ts_red, linewidth = 1.1, se = FALSE) +
  geom_hline(yintercept = 0, linetype = "dashed",
             color = ts_blue, linewidth = 0.7) +
  geom_hline(yintercept = 1.5, linetype = "dotted",
             color = ts_amber, linewidth = 0.7) +
  annotate("text", x = as.Date(paste0(yr_temp_start + 5, "-01-01")),
           y = 1.6, label = "1.5°C Paris threshold",
           color = ts_amber, hjust = 0, size = 3.2) +
  labs(title    = "Global Surface Temperature Anomaly",
       subtitle = paste("NASA GISTEMP v4 ·", yr_temp_start, "–", yr_temp_end,
                        "· Baseline 1951–1980"),
       x = NULL, y = "Temperature anomaly (°C)")

What to look for: The blue dashed line at zero is the 1951–1980 average. Values above zero mean warmer than that baseline; below zero means cooler. The early record (pre-1940) fluctuates around zero or slightly below — the 19th-century baseline was slightly cooler than 1951–1980. The red LOESS curve shows the trend — note the accelerating rise after ~1975, the approximate time when the anthropogenic greenhouse signal emerged clearly above natural variability.

6.4 STL decomposition — temperature

temp_vec <- temp |>
  as_tibble() |>
  arrange(date) |>
  mutate(anomaly = as.numeric(anomaly)) |>
  pull(anomaly)

temp_base_ts <- ts(temp_vec,
                   start     = c(yr_temp_start, 1),
                   frequency = 12)

temp_stl <- stl(temp_base_ts, s.window = "periodic", robust = TRUE)
plot(temp_stl, main = "STL Decomposition — Global Temperature Anomaly")

Comparison with CO2: Notice how different the seasonal component looks compared to CO2. The temperature seasonal cycle is very small relative to the trend — the grey bar for the seasonal panel is much taller (less important) than the trend bar. This makes physical sense: after area-averaging across the globe, Northern and Southern Hemisphere seasonal cycles partially cancel. The remainder for temperature is also much larger relative to the trend than for CO2 — year-to-year weather variability (including ENSO events) produces substantial noise.

6.5 Decade-level summaries

temp |>
  as_tibble() |>
  mutate(decade = paste0(floor(Year / 10) * 10, "s")) |>
  group_by(decade) |>
  summarise(
    n          = n(),
    Mean_anom  = round(mean(anomaly, na.rm = TRUE), 3),
    SD_anom    = round(sd(anomaly,   na.rm = TRUE), 3),
    Max_anom   = round(max(anomaly,  na.rm = TRUE), 3),
    .groups = "drop"
  ) |>
  kable(
    caption   = "Global temperature anomaly by decade. Mean > 0 = warmer than 1951–1980 baseline.",
    col.names = c("Decade", "N months", "Mean (°C)", "SD (°C)", "Max (°C)")
  ) |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover", "condensed"))

Reading the table: The mean anomaly column tells the story of 20th-century warming. Decades before ~1940 cluster near zero or slightly negative. Decades after 1980 show progressively larger positive means. The SD column shows that year-to-year variability has remained roughly constant — it is the mean that has shifted, not the variability.

6.6 ACF of temperature anomaly

Unlike the CO2 remainder (which is relatively clean), the full temperature series has both trend and seasonality. Examining the ACF of the raw anomaly series illustrates what autocorrelation looks like in a non-stationary signal.

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
acf(temp_vec,  lag.max = 48,
    main = "ACF — Temperature Anomaly (raw)",
    col  = ts_blue, lwd = 1.5)
pacf(temp_vec, lag.max = 48,
     main = "PACF — Temperature Anomaly (raw)",
     col  = ts_red, lwd = 1.5)

par(mfrow = c(1, 1))

What to look for: A slowly, linearly decaying ACF that remains significant for many lags is the classic signature of a non-stationary series (one with a trend). This is exactly what we see — the trend causes measurements that are years apart to still be highly correlated. The PACF cutting off after lag 1–2 suggests that after one differencing step, a simple AR(1) or AR(2) structure might fit the remainder. This is why ARIMA models for temperature typically use d = 1 (first difference).

7.1 What is an ARIMA model?

ARIMA stands for AutoRegressive Integrated Moving Average. It is the most widely used family of models for univariate time series forecasting. The three components address the three aspects of time series structure:

AR — AutoRegressive: Uses past values of the series to predict the current value. An AR(1) model says: “today’s value is a weighted version of yesterday’s value plus random noise.” AR(2) uses the two most recent values.

\[\text{AR}(p): \quad y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + \varepsilon_t\]

I — Integrated: Refers to differencing — subtracting each observation from the previous one to achieve stationarity. \(d = 1\) means first-difference (removes a linear trend); \(d = 2\) removes a quadratic trend.

\[\text{I}(d): \quad \Delta^d y_t = y_t - y_{t-1} \quad (d=1)\]

MA — Moving Average: Uses past forecast errors (not past values) to improve predictions. An MA(1) model says: “correct today’s forecast by a weighted version of last month’s error.”

\[\text{MA}(q): \quad y_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}\]

For seasonal data (monthly environmental records), we add a seasonal layer — the notation becomes ARIMA\((p,d,q)(P,D,Q)_m\) where the uppercase letters are the seasonal AR, differencing, and MA orders and \(m\) is the seasonal period (12 for monthly data).

7.2 Model selection: ARIMA(1,1,1)(0,1,1)[12]

For both series we fit the airline model — ARIMA\((1,1,1)(0,1,1)_{12}\) — named after Box and Jenkins’ classic 1970 analysis of monthly airline passenger data. It works well for monthly environmental series with trend and seasonality:

• \(d=1\): one non-seasonal difference removes the upward trend
• \(D=1\): one seasonal difference (subtract the value 12 months ago) removes the annual cycle
• \(p=1, q=1\): AR(1) and MA(1) terms model remaining short-range correlation
• \(Q=1\): seasonal MA(1) captures remaining seasonal correlation in residuals

After both differences, the series is approximately stationary and the AR/MA terms clean up whatever autocorrelation remains in the residuals.

co2_arima_base <- arima(co2_base_ts,
                        order    = c(1, 1, 1),
                        seasonal = list(order = c(0, 1, 1), period = 12))

co2_arima_base

## 
## Call:
## arima(x = co2_base_ts, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), 
##     period = 12))
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.0676  -0.4446  -0.9307
## s.e.  0.0911   0.0835   0.0132
## 
## sigma^2 estimated as 0.0882:  log likelihood = -178.1,  aic = 364.21

Reading the output: The ar1, ma1, and sma1 coefficients are the estimated model parameters. Values close to ±1 indicate strong AR or MA effects. The sigma^2 is the estimated variance of the error term — smaller is better. The log-likelihood (and therefore AIC) is used for model comparison.

7.3 Fitting a simple linear model as baseline

A simple ordinary least squares (OLS) linear regression on decimal date is the most intuitive baseline — it just fits a straight line through time. Its AIC is directly comparable to the ARIMA and seasonal naive models because all three are estimated by maximum likelihood on the same data.

The linear model cannot capture the seasonal cycle and it assumes a constant rate of increase — both assumptions we know are violated for CO2. The AIC comparison makes this failure visible numerically.

# Build a data frame with decimal time and CO2 values
co2_df <- tibble(
  t   = as.numeric(time(co2_base_ts)),
  co2 = as.numeric(co2_base_ts)
)

# Simple linear model: CO2 ~ time
lm_co2 <- lm(co2 ~ t, data = co2_df)

# Also compute AIC for linear model on the same scale as arima()
# lm() uses OLS; AIC() works directly on lm objects
aic_lm     <- AIC(lm_co2)
bic_lm     <- BIC(lm_co2)

summary(lm_co2)

## 
## Call:
## lm(formula = co2 ~ t, data = co2_df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.375 -3.580 -2.013  3.466 11.318 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.968e+03  1.685e+01  -176.2   <2e-16 ***
## t            1.671e+00  8.456e-03   197.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.803 on 821 degrees of freedom
## Multiple R-squared:  0.9794, Adjusted R-squared:  0.9794 
## F-statistic: 3.906e+04 on 1 and 821 DF,  p-value: < 2.2e-16

co2_snaive <- arima(co2_base_ts, order = c(0, 0, 0),
                    seasonal = list(order = c(0, 1, 0), period = 12))

tibble(
  Model      = c("Linear regression (OLS)",
                 "Seasonal naive (repeat last year)",
                 "ARIMA(1,1,1)(0,1,1)[12]"),
  Parameters = c(2, 1, 4),
  AIC        = round(c(aic_lm,
                       AIC(co2_snaive),
                       AIC(co2_arima_base)), 2),
  BIC        = round(c(bic_lm,
                       BIC(co2_snaive),
                       BIC(co2_arima_base)), 2),
  Notes      = c("Assumes constant rate, ignores seasonality",
                 "No trend — repeats previous year's month",
                 "Models trend, seasonality, and autocorrelation")
) |>
  kable(caption = "AIC/BIC comparison across three model families — lower is better.") |>
  kable_styling(full_width = FALSE)

Reading the table: The linear model has the highest AIC because it cannot capture the seasonal cycle — every summer dip and winter peak is left in the residuals as unexplained structure. The seasonal naive has a lower AIC because it correctly tracks the seasonal pattern, but it ignores the trend. ARIMA captures both, earning the lowest AIC with only four parameters. This is why time-series models exist: linear regression applied to autocorrelated data is not just sub-optimal — it leaves systematic structure on the table and overstates confidence.

7.4 CO2 forecast — three models compared

The real payoff of the AIC comparison is visible when you overlay the forecasts. The linear model projects a straight line; ARIMA continues the observed curvature and seasonal cycle; seasonal naive flatlines (no trend).

h_co2    <- 60   # 5 years ahead
co2_pred <- predict(co2_arima_base, n.ahead = h_co2)

last_t   <- time(co2_base_ts)[length(co2_base_ts)]
fc_times <- seq(last_t, by = 1/12, length.out = h_co2 + 1)[-1]

# ── ARIMA forecast ────────────────────────────────────────────────────────────
fc_df <- tibble(
  t    = fc_times,
  fit  = as.numeric(co2_pred$pred),
  lo95 = as.numeric(co2_pred$pred - 1.96 * co2_pred$se),
  hi95 = as.numeric(co2_pred$pred + 1.96 * co2_pred$se)
)

# ── Linear model forecast ─────────────────────────────────────────────────────
lm_fc_df <- tibble(t = fc_times) |>
  mutate(
    fit_lm  = predict(lm_co2, newdata = tibble(t = fc_times)),
    # propagate prediction interval using sigma and leverage
    se_lm   = summary(lm_co2)$sigma *
                sqrt(1 + 1/nrow(co2_df) +
                     (fc_times - mean(co2_df$t))^2 /
                     sum((co2_df$t - mean(co2_df$t))^2)),
    lo95_lm = fit_lm - 1.96 * se_lm,
    hi95_lm = fit_lm + 1.96 * se_lm
  )

# ── Seasonal naive forecast ───────────────────────────────────────────────────
snaive_pred <- predict(co2_snaive, n.ahead = h_co2)
snaive_df   <- tibble(
  t       = fc_times,
  fit_sn  = as.numeric(snaive_pred$pred)
)

# ── Historical last 10 years ──────────────────────────────────────────────────
hist_df <- tibble(
  t   = as.numeric(time(co2_base_ts)),
  val = as.numeric(co2_base_ts)
) |> filter(t >= last_t - 10)

ggplot() +
  # ARIMA 95% PI band
  geom_ribbon(data = fc_df, aes(x = t, ymin = lo95, ymax = hi95),
              fill = ts_red, alpha = 0.15) +
  # Linear model 95% PI band
  geom_ribbon(data = lm_fc_df, aes(x = t, ymin = lo95_lm, ymax = hi95_lm),
              fill = ts_blue, alpha = 0.12) +
  # Historical data
  geom_line(data = hist_df, aes(x = t, y = val),
            color = ts_green, linewidth = 0.9) +
  # ARIMA forecast line
  geom_line(data = fc_df, aes(x = t, y = fit, color = "ARIMA(1,1,1)(0,1,1)[12]"),
            linewidth = 1.1, linetype = "dashed") +
  # Linear model forecast line
  geom_line(data = lm_fc_df, aes(x = t, y = fit_lm,
                                  color = "Linear regression"),
            linewidth = 1.1, linetype = "dotted") +
  # Seasonal naive forecast line
  geom_line(data = snaive_df, aes(x = t, y = fit_sn,
                                   color = "Seasonal naive"),
            linewidth = 0.9, linetype = "longdash", alpha = 0.7) +
  scale_color_manual(
    name   = "Model",
    values = c("ARIMA(1,1,1)(0,1,1)[12]" = ts_red,
               "Linear regression"        = ts_blue,
               "Seasonal naive"           = ts_amber)
  ) +
  labs(title    = "CO2 Forecast Comparison — 5 Years Ahead",
       subtitle = "Shaded bands = 95% prediction intervals · Green = observed",
       x = NULL, y = "CO2 (ppm)") +
  theme(legend.position = "bottom")

What the comparison shows: - Linear regression (blue dotted): Extrapolates a straight line. Because the actual rate of CO2 increase has been accelerating, a line fit over the whole record underestimates recent years and projects a rate that is already too slow. The 95% PI (blue band) widens substantially because OLS intervals do not account for autocorrelation — they assume each month is independent, which massively underestimates true uncertainty. This is the “wrong” model. - Seasonal naive (amber dashed): Continues the seasonal cycle from the most recent year but projects zero trend. Useful as a lower benchmark — anything that can’t beat this has no value. - ARIMA (red dashed): Projects the trend continuing with the seasonal cycle, and its 95% PI (red band) properly accounts for autocorrelation in the uncertainty. This is the appropriate model for this data structure.

7.5 Temperature forecast — 10 years

temp_arima_base <- arima(temp_base_ts,
                         order    = c(1, 1, 1),
                         seasonal = list(order = c(0, 1, 1), period = 12))

# ── Linear model for temperature ──────────────────────────────────────────────
temp_df <- tibble(
  t   = as.numeric(time(temp_base_ts)),
  tmp = as.numeric(temp_base_ts)
)
lm_temp <- lm(tmp ~ t, data = temp_df)

h_temp       <- 120   # 10 years
temp_pred    <- predict(temp_arima_base, n.ahead = h_temp)
last_t_temp  <- time(temp_base_ts)[length(temp_base_ts)]
fc_times_t   <- seq(last_t_temp, by = 1/12, length.out = h_temp + 1)[-1]

# ── ARIMA forecast ────────────────────────────────────────────────────────────
fc_temp_df <- tibble(
  t    = fc_times_t,
  fit  = as.numeric(temp_pred$pred),
  lo95 = as.numeric(temp_pred$pred - 1.96 * temp_pred$se),
  hi95 = as.numeric(temp_pred$pred + 1.96 * temp_pred$se)
)

# ── Linear model forecast ─────────────────────────────────────────────────────
lm_temp_fc <- tibble(t = fc_times_t) |>
  mutate(
    fit_lm  = predict(lm_temp, newdata = tibble(t = fc_times_t)),
    se_lm   = summary(lm_temp)$sigma *
                sqrt(1 + 1/nrow(temp_df) +
                     (fc_times_t - mean(temp_df$t))^2 /
                     sum((temp_df$t - mean(temp_df$t))^2)),
    lo95_lm = fit_lm - 1.96 * se_lm,
    hi95_lm = fit_lm + 1.96 * se_lm
  )

hist_temp_df <- tibble(
  t   = as.numeric(time(temp_base_ts)),
  val = as.numeric(temp_base_ts)
) |> filter(t >= last_t_temp - 30)

ggplot() +
  # ARIMA 95% PI
  geom_ribbon(data = fc_temp_df, aes(x = t, ymin = lo95, ymax = hi95),
              fill = ts_red, alpha = 0.15) +
  # Linear model 95% PI
  geom_ribbon(data = lm_temp_fc, aes(x = t, ymin = lo95_lm, ymax = hi95_lm),
              fill = ts_blue, alpha = 0.12) +
  # Historical data
  geom_line(data = hist_temp_df, aes(x = t, y = val),
            color = "grey50", linewidth = 0.5) +
  # ARIMA forecast
  geom_line(data = fc_temp_df, aes(x = t, y = fit,
                                    color = "ARIMA(1,1,1)(0,1,1)[12]"),
            linewidth = 1.1, linetype = "dashed") +
  # Linear model forecast
  geom_line(data = lm_temp_fc, aes(x = t, y = fit_lm,
                                    color = "Linear regression"),
            linewidth = 1.1, linetype = "dotted") +
  # Paris threshold
  geom_hline(yintercept = 1.5, linetype = "dotted",
             color = ts_amber, linewidth = 0.9) +
  annotate("text", x = last_t_temp, y = 1.58,
           label = "1.5°C Paris threshold",
           color = ts_amber, hjust = 0, size = 3.2) +
  scale_color_manual(
    name   = "Model",
    values = c("ARIMA(1,1,1)(0,1,1)[12]" = ts_red,
               "Linear regression"        = ts_blue)
  ) +
  labs(title    = "Global Temperature Anomaly Forecast — 10 Years",
       subtitle = "ARIMA vs linear regression · 95% prediction intervals · Grey = observed",
       x = NULL, y = "Temperature anomaly (°C)") +
  theme(legend.position = "bottom")

Comparing the two models on temperature: For temperature the difference between linear regression and ARIMA is subtler than for CO2. Unlike CO2, where the curvature is obvious, global mean temperature has fluctuated enough that a straight line is not a terrible description of the 20th-century record. However, the ARIMA model still wins on AIC because it correctly models the month-to-month autocorrelation — meaning its prediction intervals are calibrated honestly, while the linear regression’s are too narrow by an unknown factor.

# AIC table for temperature models
temp_snaive <- arima(temp_base_ts, order = c(0, 0, 0),
                     seasonal = list(order = c(0, 1, 0), period = 12))

tibble(
  Model      = c("Linear regression (OLS)",
                 "Seasonal naive",
                 "ARIMA(1,1,1)(0,1,1)[12]"),
  Parameters = c(2, 1, 4),
  AIC        = round(c(AIC(lm_temp),
                       AIC(temp_snaive),
                       AIC(temp_arima_base)), 2),
  BIC        = round(c(BIC(lm_temp),
                       BIC(temp_snaive),
                       BIC(temp_arima_base)), 2)
) |>
  kable(caption = "AIC/BIC — temperature models. Lower = better fit per parameter used.") |>
  kable_styling(full_width = FALSE)

8.1 What are model residuals and why do they matter?

A model’s residuals are the differences between the observed values and the model’s fitted values:

\[e_t = y_t - \hat{y}_t\]

If the model has captured all the systematic structure in the data, the residuals should be white noise — random, with no pattern, centred on zero, and with no remaining autocorrelation. If you can still see patterns in the residuals, the model has not fully explained the data and your forecasts and prediction intervals are unreliable.

Four diagnostic checks:

1. Residuals vs time — should scatter randomly around zero; any trend, cycles, or changing variance indicates model failure
2. ACF of residuals — all bars should be within the blue lines; any significant spike means unused structure remains
3. PACF of residuals — same; the PACF helps identify what kind of structure is left (AR or MA)
4. Normal Q-Q plot — points should follow the diagonal; heavy tails or systematic curves indicate non-normal errors, which can invalidate prediction intervals

8.2 CO2 ARIMA residual diagnostics

co2_arima_resid <- residuals(co2_arima_base)

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(co2_arima_resid, main = "Residuals over time",
     ylab = "Residual (ppm)", col = ts_blue, type = "l")
abline(h = 0, col = ts_red, lty = 2, lwd = 1.5)
acf(co2_arima_resid,  main = "ACF of residuals",  col = ts_blue, lwd = 1.5)
pacf(co2_arima_resid, main = "PACF of residuals", col = ts_red,  lwd = 1.5)
qqnorm(co2_arima_resid, main = "Normal Q-Q", col = ts_sage, pch = 19, cex = 0.5)
qqline(co2_arima_resid, col = ts_red, lwd = 1.5)

par(mfrow = c(1, 1))

What to look for: - Top left: Residuals should bounce around zero with roughly constant spread. A funnel shape (increasing variance over time) suggests the model should be fit on a log scale. - Top right (ACF): Ideally all bars within the blue lines. A significant spike at lag 12 means the seasonal model isn’t quite right. - Bottom left (PACF): Same logic. A sharp cutoff after lag 1 would suggest adding another AR term. - Bottom right (Q-Q): Points close to the diagonal = approximately normal residuals. S-curves or heavy tails mean outlier years are distorting the distribution.

8.3 Temperature ARIMA residual diagnostics

temp_arima_resid <- residuals(temp_arima_base)

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(temp_arima_resid, main = "Residuals over time",
     ylab = "Residual (°C)", col = ts_blue, type = "l")
abline(h = 0, col = ts_red, lty = 2, lwd = 1.5)
acf(temp_arima_resid,  main = "ACF of residuals",  col = ts_blue, lwd = 1.5)
pacf(temp_arima_resid, main = "PACF of residuals", col = ts_red,  lwd = 1.5)
qqnorm(temp_arima_resid, main = "Normal Q-Q", col = ts_sage, pch = 19, cex = 0.5)
qqline(temp_arima_resid, col = ts_red, lwd = 1.5)

par(mfrow = c(1, 1))

Temperature-specific interpretation: Temperature residuals often show larger spikes around major El Niño years (e.g., 1997–1998, 2015–2016). These are genuine climate events that no statistical time series model can anticipate from past data alone — they are by definition anomalies. If you see consistent large residuals for specific decades, that might indicate the warming rate has changed and the model structure needs updating.

8.4 Formal Ljung-Box test on model residuals

lb_co2  <- Box.test(co2_arima_resid,  lag = 24, type = "Ljung-Box")
lb_temp <- Box.test(temp_arima_resid, lag = 24, type = "Ljung-Box")

tibble(
  Series   = c("CO2 (Mauna Loa)", "Global temperature anomaly"),
  LB_stat  = round(c(lb_co2$statistic,  lb_temp$statistic), 3),
  df       = c(lb_co2$parameter, lb_temp$parameter),
  p_value  = round(c(lb_co2$p.value, lb_temp$p.value), 4),
  Adequate = ifelse(c(lb_co2$p.value, lb_temp$p.value) > 0.05,
                    "✓ White noise", "✗ Structure remains")
) |>
  kable(
    caption   = "Ljung-Box test on ARIMA model residuals. p > 0.05 means no significant autocorrelation remains.",
    col.names = c("Series", "LB statistic", "df", "p-value", "Adequate?")
  ) |>
  kable_styling(full_width = FALSE)

If the Ljung-Box test is significant (p < 0.05): The residuals still contain autocorrelation structure the ARIMA model failed to capture. Options: (1) Increase the AR or MA order (try \(p=2\) or \(q=2\)); (2) Add a seasonal AR term (\(P=1\)); (3) Fit separate models for different sub-periods if the autocorrelation structure appears to have changed over time; (4) Accept the imperfect fit for exploratory purposes but note the limitation when reporting.

9.1 Other environmental time series to explore

The workflow in this module applies directly to many other environmental datasets:

9.2 When ARIMA is not enough

The airline model works well for many environmental series, but some situations require different approaches:

Intervention analysis: If a known event (policy change, volcanic eruption, station relocation) caused a step change or slope change in the series, add an intervention term as a dummy variable in the ARIMA regression component.

Transfer function models (ARIMAX): If you believe one time series (e.g., ENSO index) drives another (e.g., California precipitation), include the driver as an exogenous predictor. This becomes arima(..., xreg = enso_index).

Structural breaks: If the autocorrelation or trend structure appears to have changed at a specific time (e.g., warming accelerated after 1980), fit separate models to each period and test whether the parameters differ.

Multiple seasonality: Hourly environmental data (air quality monitors, stream gauges) has both diurnal (24-hour) and annual seasonality. Base R’s stl() handles only one seasonal period at a time — use stats::HoltWinters() with multiple seasons or the forecast package’s mstl().

9.3 Cross-validation for honest forecast accuracy

A single train/test split gives one estimate of forecast accuracy. Rolling-origin cross-validation (also called time-series cross-validation) gives a more robust estimate by repeatedly re-fitting the model as new data become available:

# Rolling-origin cross-validation for ARIMA — base R approach
# Train on first 20 years, test on year 21; then first 21 years, test on year 22...
n      <- length(co2_base_ts)
h      <- 12          # forecast horizon (1 year)
init   <- 12 * 20     # initial training set (20 years)
errors <- numeric(n - init - h + 1)

for (i in seq_along(errors)) {
  train <- co2_base_ts[1:(init + i - 1)]
  test  <- co2_base_ts[(init + i):(init + i + h - 1)]
  fit   <- arima(train,
                 order    = c(1, 1, 1),
                 seasonal = list(order = c(0, 1, 1), period = 12))
  fc    <- predict(fit, n.ahead = h)$pred
  errors[i] <- sqrt(mean((as.numeric(test) - as.numeric(fc))^2))  # RMSE
}

cat("Mean RMSE across", length(errors), "rolling origins:",
    round(mean(errors), 3), "ppm\n")

The resulting mean RMSE gives you a realistic expectation of how well the model will forecast future data — much more honest than in-sample fit statistics.

What’s next?

• Signal Over Noise — the foundation of statistical testing that precedes time-series methods
• Functional Programming — use purrr::map() to fit models to many stations at once
• Visualization — ggplot2 techniques for time-series figures

data(AirPassengers)
plot(AirPassengers, main="Monthly airline passengers")
dec <- decompose(AirPassengers, type = "multiplicative")
plot(dec)

resid_ts <- dec$random
acf(resid_ts, na.action = na.pass)
pacf(resid_ts, na.action = na.pass)