library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly)      # interactive figures and ALL 3-D graphics
library(forecast)
library(tseries)

How this chapter uses graphics

Every two-dimensional figure is drawn with ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=FALSE, echo=TRUE.

Every three-dimensional figure is drawn with plot_ly() and evaluated. Temporal dependence is inherently two-parameter, the autocorrelation is a function of lag and of the model coefficient, prediction intervals of horizon and persistence, covariance structures of lag and \(\rho\), and reading these from a fixed viewpoint discards the interaction that is the subject.

1 Learning objectives

After completing this chapter you will be able to:

State the stationarity and invertibility conditions for ARMA models in terms of the roots of the characteristic polynomials.
Interpret the ACF and PACF, with correct sampling bounds, and use them to identify candidate orders.
Test for a unit root using ADF and KPSS together, and read the four possible outcomes.
Explain why smoothing a series before fitting an ARIMA model fabricates autocorrelation and collapses forecast intervals.
Fit seasonal ARIMA models and check residuals with the Ljung–Box test.
Derive prediction intervals from the MA(\(\infty\)) representation and explain why they widen with horizon.
Evaluate forecasts against naive benchmarks using scale-free metrics (MASE, RMSSE) under rolling-origin cross-validation.
Fit linear and generalized linear mixed models, and choose a covariance structure.
Distinguish GEE from GLMM by target of inference, apply the sandwich variance, and use the attenuation relation between them.
Explain which longitudinal methods remain valid under MCAR, MAR, and MNAR dropout.

Estimated time: 14–18 hours including exercises. Prerequisites: Chapter 3 (matrix computing), Chapter 6 (temporal validation, §6.9), Chapter 9 (rolling-origin CV and forecast metrics, §9.17), and Chapter 10 (streaming and drift, §10.13).

2 PART I: TEMPORAL FOUNDATIONS

3 What makes temporal data different

Every method in Chapters 5–11 assumed observations were exchangeable. Time series and repeated measures violate that in two distinct ways, and the two require different machinery.

	Time series	Longitudinal / panel
Structure	One (or few) long sequences	Many short sequences, one per unit
\(n\) grows by	More time points	More units
Dependence	Within the sequence, across time	Within unit, across occasions
Asymptotics	\(T\to\infty\)	\(N\to\infty\), \(n_i\) fixed
Typical goal	Forecast the future	Estimate an effect, accounting for correlation
Core tools	ARIMA, state space, spectral	Mixed models, GEE, SEM

The distinction matters because the same word, “correlated observations” — points at opposite estimation problems. Forecasting exploits the dependence; longitudinal modeling treats it as a nuisance to be accounted for so that effect estimates and their standard errors are honest.

4 Stationarity

A process \(\{X_t\}\) is weakly (second-order) stationary if

\[ \begin{aligned} \mathbb E[X_t]&=\mu &&\text{(constant mean)}\\ \operatorname{Var}(X_t)&=\sigma^2<\infty &&\text{(constant variance)}\\ \operatorname{Cov}(X_t,X_{t+h})&=\gamma(h) &&\text{(covariance depends only on the lag \(h\))} \end{aligned} \]

Stationarity is what makes estimation possible at all: it is the assumption that lets a single realization stand in for the ensemble average, so that \(\hat\gamma(h)=\frac1T\sum_t(x_t-\bar x)(x_{t+h}-\bar x)\) estimates something fixed.

Common misconception: “stationary means the series has no trend.” A constant mean is only one of three requirements. A series can be trend-free and still non-stationary through changing variance (volatility clustering, common in financial and physiological signals) or through a changing autocorrelation structure (a regime shift).

Conversely, a series can look strongly patterned and be stationary: an AR(1) with \(\phi=0.95\) wanders in long excursions that resemble a trend, but its mean, variance, and autocovariance are all constant and it reverts. Distinguishing “trend-stationary” from “difference-stationary” is precisely what the unit-root tests of §12.5 are for.

4.1 Conditions on the polynomials

Write the backshift operator \(B x_t = x_{t-1}\). An ARMA(\(p,q\)) process is

\[\phi(B)\,X_t=\theta(B)\,\varepsilon_t, \qquad \phi(z)=1-\phi_1z-\cdots-\phi_pz^p, \quad \theta(z)=1+\theta_1z+\cdots+\theta_qz^q,\]

with \(\varepsilon_t\sim\mathrm{WN}(0,\sigma^2)\), white noise of variance \(\sigma^2\), a parameter to be estimated, not fixed at 1.

Stationarity and invertibility. \[\boxed{\;\text{stationary}\iff \text{all roots of }\phi(z)\text{ lie strictly \emph{outside} the unit circle}\;}\] \[\boxed{\;\text{invertible}\iff \text{all roots of }\theta(z)\text{ lie strictly \emph{outside} the unit circle}\;}\]

For AR(1), \(\phi(z)=1-\phi z\) has root \(1/\phi\), so stationarity is \(|\phi|<1\). A root on the unit circle (\(\phi=1\)) is a unit root, a random walk, which is not stationary and does not revert.

Invertibility matters because a non-invertible MA has an observationally equivalent invertible twin with the same autocovariance; requiring invertibility makes the parameterization unique and lets \(\varepsilon_t\) be recovered from the observed past.

set.seed(11)
T_s <- 400
sim_ar1 <- function(phi, n = T_s) as.numeric(arima.sim(list(ar = phi), n))

bind_rows(
  data.frame(t = 1:T_s, x = sim_ar1(0.5),  series = "phi = 0.50  (root 2.0)"),
  data.frame(t = 1:T_s, x = sim_ar1(0.95), series = "phi = 0.95  (root 1.05)"),
  data.frame(t = 1:T_s, x = cumsum(rnorm(T_s)), series = "phi = 1.00  (unit root)")) |>
  ggplot(aes(t, x)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_line(linewidth = 0.4, color = "steelblue") +
  facet_wrap(~ series, ncol = 1, scales = "free_y") +
  labs(title = "Approaching the unit circle",
       subtitle = "As the root approaches 1 the excursions lengthen; at the root the process stops reverting",
       x = "Time", y = NULL) +
  theme_dspa(10)

The middle panel is stationary despite looking trended, the root at 1.05 is outside the unit circle, so the process reverts, just slowly. The bottom panel has its root on the circle and does not revert at all. Telling these apart from a finite sample is genuinely hard, which is why unit-root tests have low power.

5 Autocorrelation

\[\rho(h)=\frac{\gamma(h)}{\gamma(0)},\qquad \hat\rho(h)=\frac{\sum_{t=1}^{T-h}(x_t-\bar x)(x_{t+h}-\bar x)}{\sum_{t=1}^{T}(x_t-\bar x)^2}\]

The partial autocorrelation \(\alpha(h)\) is the correlation between \(X_t\) and \(X_{t+h}\) after removing the linear effect of the intervening \(X_{t+1},\dots,X_{t+h-1}\).

The pair identifies model orders through a duality:

	ACF	PACF
AR(\(p\))	tails off (geometric / damped sine)	cuts off after lag \(p\)
MA(\(q\))	cuts off after lag \(q\)	tails off
ARMA(\(p,q\))	tails off	tails off

Why the PACF cuts off for AR. An AR(\(p\)) is exactly a regression of \(X_t\) on its own \(p\) previous values. Once those are conditioned on, nothing further back adds information, so \(\alpha(h)=0\) for \(h>p\) by construction.

The sampling bounds must come from the series being plotted. Under a white-noise null, \(\hat\rho(h)\approx N(0,1/T)\), giving the familiar \(\pm1.96/\sqrt T\). \(T\) is the length of the series whose ACF is shown, a bound computed from a different, longer or shorter series mis-scales every significance judgment on the plot.

For residual ACFs after fitting ARMA(\(p,q\)), even that bound is wrong: the low-lag residual autocorrelations have variance smaller than \(1/T\) because \(p+q\) parameters were fitted to remove them. The correct omnibus check is the Ljung–Box statistic with degrees of freedom reduced by \(p+q\) (§12.10).

# One helper, used everywhere. The bound is always computed from THIS series.
acf_frame <- function(x, lag.max = 36, type = c("correlation", "partial")) {
  type <- match.arg(type)
  a <- stats::acf(x, lag.max = lag.max, plot = FALSE, type = type,
                  na.action = na.pass)
  n <- sum(!is.na(x))
  data.frame(lag = as.numeric(a$lag), value = as.numeric(a$acf),
             bound = stats::qnorm(0.975) / sqrt(n),
             kind = ifelse(type == "correlation", "ACF", "PACF"))
}

plot_acf <- function(x, lag.max = 36, title = "") {
  d <- bind_rows(acf_frame(x, lag.max, "correlation") |> filter(lag > 0),
                 acf_frame(x, lag.max, "partial"))
  ggplot(d, aes(lag, value)) +
    geom_hline(yintercept = 0, color = "grey60") +
    geom_ribbon(aes(ymin = -bound, ymax = bound), fill = "grey88") +
    geom_segment(aes(xend = lag, yend = 0), linewidth = 0.7, color = "steelblue") +
    facet_wrap(~ kind, ncol = 1) +
    labs(title = title,
         subtitle = sprintf("Shaded band: +/- 1.96/sqrt(n) with n = %d, the length of THIS series",
                            sum(!is.na(x))),
         x = "Lag", y = NULL) +
    theme_dspa(10)
}

set.seed(21)
x_ar2 <- as.numeric(arima.sim(list(ar = c(0.6, 0.25)), 600))
x_ma2 <- as.numeric(arima.sim(list(ma = c(0.7, 0.5)), 600))

plot_acf(x_ar2, 24, "AR(2): PACF cuts off after lag 2, ACF tails off") /
  plot_acf(x_ma2, 24, "MA(2): ACF cuts off after lag 2, PACF tails off")

# --- Interactive equivalent ------------------------------------------------
d <- bind_rows(acf_frame(x_ar2, 24, "correlation") |> filter(lag > 0),
               acf_frame(x_ar2, 24, "partial"))
plot_ly(d, x = ~lag, y = ~value, color = ~kind, type = "bar") |>
  add_lines(x = ~range(lag), y = ~c(d$bound[1], d$bound[1]), name = "Upper bound",
            line = list(dash = "dash", color = "grey")) |>
  add_lines(x = ~range(lag), y = ~c(-d$bound[1], -d$bound[1]), name = "Lower bound",
            line = list(dash = "dash", color = "grey")) |>
  layout(title = "ACF and PACF", barmode = "group",
         xaxis = list(title = "Lag"), yaxis = list(title = "Correlation"))

The theoretical ACF of an AR(1) is \(\rho(h)=\phi^{|h|}\), a two-parameter family, and worth rotating:

lag_grid <- 0:30
phi_grid <- seq(-0.95, 0.95, length.out = 60)
Zacf <- outer(phi_grid, lag_grid, function(p, h) p^h)

plot_ly(x = lag_grid, y = phi_grid, z = Zacf, type = "surface",
        colorscale = "RdBu", colorbar = list(title = "rho(h)")) |>
  layout(title = "Theoretical AR(1) autocorrelation over lag and phi",
         scene = list(xaxis = list(title = "Lag h"),
                      yaxis = list(title = "AR coefficient phi"),
                      zaxis = list(title = "rho(h)")))

Rotate along the \(\phi\) axis. Near \(\phi=0\) the surface collapses to a spike at lag 0, white noise. As \(|\phi|\to1\) the decay flattens toward a ridge that never reaches zero, which is the geometry of a near-unit root. Negative \(\phi\) gives the alternating sheet.

6 The Wold decomposition

Why should ARMA models be the right family at all? Because of a representation theorem.

Wold decomposition. Every zero-mean weakly stationary process \(\{X_t\}\) admits a unique decomposition \[X_t=\underbrace{\sum_{j=0}^{\infty}\psi_j\varepsilon_{t-j}}_{\text{stochastic, }\psi_0=1,\ \sum\psi_j^2<\infty}+\underbrace{V_t}_{\text{deterministic}},\] where \(\{\varepsilon_t\}\) is white noise uncorrelated with \(V_t\).

Every stationary process is a linear filter of white noise. ARMA models are the family of rational such filters, \(\psi(B)=\theta(B)/\phi(B)\), and any \(\mathrm{MA}(\infty)\) with square-summable weights can be approximated arbitrarily well by a rational one with finitely many parameters. That is the justification for the whole enterprise.

The \(\psi\)-weights also give the forecast variance directly:

\[\boxed{\;\operatorname{Var}\big(X_{T+h}\mid \mathcal F_T\big)=\sigma^2\sum_{j=0}^{h-1}\psi_j^2\;}\]

which is why prediction intervals widen with horizon, and by exactly how much (§12.11).

# psi-weights of an ARMA, and the forecast variance they imply
psi_of <- function(ar = numeric(0), ma = numeric(0), n = 20)
  c(1, ARMAtoMA(ar = ar, ma = ma, lag.max = n - 1))

data.frame(
  h = 1:8,
  `AR(1) phi=0.5`  = round(cumsum(psi_of(ar = 0.5)[1:8]^2), 4),
  `AR(1) phi=0.9`  = round(cumsum(psi_of(ar = 0.9)[1:8]^2), 4),
  `MA(2)`          = round(cumsum(psi_of(ma = c(0.7, 0.5))[1:8]^2), 4),
  check.names = FALSE)

Read down the columns: the MA(2) forecast variance stops growing after \(h=3\), because an MA(\(q\)) is unpredictable beyond \(q\) steps and the interval saturates at the unconditional variance. The AR variances keep growing, and faster for larger \(\phi\).

7 Unit-root testing

Two tests with opposite nulls, and they are meant to be read together.

\[ \begin{aligned} \textbf{ADF (augmented Dickey--Fuller): }\quad &H_0:\ \text{unit root (non-stationary)}\\ \textbf{KPSS: }\quad &H_0:\ \text{(trend-)stationary} \end{aligned} \]

ADF	KPSS	Conclusion
Reject	Fail to reject	Stationary — both agree
Fail to reject	Reject	Unit root — both agree; difference the series
Fail to reject	Fail to reject	Inconclusive — not enough data to decide
Reject	Reject	Neither — suggests heteroscedasticity or a structural break

Common misconception: “the ADF test rejected, so the series is stationary.” Rejecting a unit-root null is evidence against a unit root, not evidence for stationarity, and the ADF test is well known to have low power against near-unit-root alternatives, so failing to reject frequently means only that \(T\) is small.

Running KPSS alongside converts a one-sided statement into a four-way decision, and the two “agree” cells are the only ones that support a confident conclusion.

set.seed(31)
series_list <- list(
  `white noise`          = rnorm(400),
  `AR(1), phi = 0.5`     = sim_ar1(0.5),
  `AR(1), phi = 0.95`    = sim_ar1(0.95),
  `random walk`          = cumsum(rnorm(400)),
  `RW, differenced`      = diff(cumsum(rnorm(400))))

do.call(rbind, lapply(names(series_list), function(nm) {
  x <- series_list[[nm]]
  adf <- suppressWarnings(tseries::adf.test(x, alternative = "stationary"))
  kps <- suppressWarnings(tseries::kpss.test(x, null = "Level"))
  data.frame(series = nm,
             adf_p = round(adf$p.value, 4),
             adf_rejects_unit_root = adf$p.value < 0.05,
             kpss_p = round(kps$p.value, 4),
             kpss_rejects_stationarity = kps$p.value < 0.05,
             ndiffs_suggested = forecast::ndiffs(x))
}))

forecast::ndiffs() automates the decision by applying the test repeatedly until it stops rejecting. Note the fourth row: the random walk needs one difference, and the fifth row confirms that one difference suffices.

8 PART II: ARIMA

9 The Beijing PM2.5 series

The Beijing PM2.5 dataset records hourly concentrations of particulate matter under 2.5 microns, expressed as an Air Quality Index, from 2008 to 2016, roughly 69,000 hourly observations.

pm_raw <- if (ONLINE) {
  dspa_read("https://umich.instructure.com/files/1823138/download?download_frd=1",
            "beijing_pm25.csv")
} else NULL

if (is.null(pm_raw)) {
  # Synthetic substitute with the same structure: daily cycle, annual cycle,
  # AR(1) persistence, right skew, and a realistic missingness rate.
  set.seed(41); n_h <- 69335
  tt <- seq_len(n_h)
  base <- 90 + 25 * sin(2*pi*tt/24 - 1.2) + 40 * sin(2*pi*tt/8766 + 0.5)
  ar <- as.numeric(arima.sim(list(ar = 0.92), n_h, sd = 18))
  val <- pmax(0, round(base + ar))
  val[sample(n_h, round(0.065 * n_h))] <- -999
  pm_raw <- data.frame(
    Date..LST. = format(seq(as.POSIXct("2008-04-08 15:00:00", tz = "UTC"),
                            by = "hour", length.out = n_h), "%m/%d/%Y %H:%M"),
    Value = val)
  message("Note: Beijing file unavailable; using a synthetic substitute.")
}
c(rows = nrow(pm_raw), columns = ncol(pm_raw))

#>    rows columns 
#>   69335      11

library(lubridate)

# The timestamp format uses %H (24-hour), not %h -- which is an alias for the
# abbreviated MONTH name and would silently drop the time component.
pm <- pm_raw |>
  mutate(datetime = as.POSIXct(Date..LST., format = "%m/%d/%Y %H:%M", tz = "UTC"),
         value = ifelse(Value == -999, NA_real_, Value)) |>
  filter(!is.na(datetime)) |>
  arrange(datetime) |>
  select(datetime, value)

c(rows = nrow(pm),
  missing = sum(is.na(pm$value)),
  missing_pct = round(100 * mean(is.na(pm$value)), 2),
  range_start = format(min(pm$datetime)), range_end = format(max(pm$datetime)),
  negative_values_remaining = sum(pm$value < 0, na.rm = TRUE))

#>                      rows                   missing               missing_pct 
#>                   "69335"                    "4408"                    "6.36" 
#>               range_start                 range_end negative_values_remaining 
#>     "2008-04-08 15:00:00"     "2016-04-30 23:00:00"                       "8"

The \(-999\) codes are converted to NA once; no replacement sentinel is introduced. Every downstream step treats missingness as missingness.

9.1 Imputing temporal gaps

Common misconception: “fill the gaps with the mean, then proceed.” For cross-sectional data mean imputation is merely inefficient. For a time series it is destructive: it inserts flat runs, breaks the local autocorrelation at every imputed point, and biases the sample ACF toward zero, precisely the quantity the whole analysis is about to estimate.

Temporal imputation must respect temporal structure. Linear interpolation is the minimum; Kalman smoothing (imputeTS::na_kalman()) is the principled choice, since it returns the conditional mean of the missing values under a fitted state-space model and therefore preserves both the level and the dependence.

library(imputeTS)

pm_mean <- pm$value; pm_mean[is.na(pm_mean)] <- mean(pm$value, na.rm = TRUE)
pm_interp <- imputeTS::na_interpolation(pm$value)
pm_kalman <- dspa_try(imputeTS::na_kalman(ts(pm$value, frequency = 24),
                                          model = "auto.arima"),
                      fallback = imputeTS::na_seadec(ts(pm$value, frequency = 24)),
                      label = "Kalman imputation")
pm_kalman <- as.numeric(pm_kalman)

acf_at <- function(x, h) as.numeric(acf(x, lag.max = h, plot = FALSE)$acf[h + 1])
data.frame(
  method = c("mean imputation", "linear interpolation", "Kalman smoothing"),
  acf_lag1  = round(c(acf_at(pm_mean, 1),  acf_at(pm_interp, 1),  acf_at(pm_kalman, 1)), 4),
  acf_lag24 = round(c(acf_at(pm_mean, 24), acf_at(pm_interp, 24), acf_at(pm_kalman, 24)), 4),
  sd        = round(c(sd(pm_mean), sd(pm_interp), sd(pm_kalman)), 2))

Mean imputation attenuates the lag-1 autocorrelation and shrinks the variance, because 6.5% of the series has been replaced by a constant. Kalman smoothing preserves both.

win <- 3000:3400
bind_rows(
  data.frame(t = win, v = pm_mean[win],   method = "Mean imputation"),
  data.frame(t = win, v = pm_kalman[win], method = "Kalman smoothing")) |>
  ggplot(aes(t, v, color = method)) +
  geom_line(linewidth = 0.5) +
  facet_wrap(~ method, ncol = 1) +
  scale_color_manual(values = c(`Mean imputation` = "#D8433B",
                                 `Kalman smoothing` = "#3B7DD8"), guide = "none") +
  labs(title = "The same window under two imputation schemes",
       subtitle = "Mean imputation inserts flat segments at the imputed points; Kalman smoothing does not",
       x = "Hour index", y = "PM2.5 AQI") +
  theme_dspa(10)

pm$value_imp <- pm_kalman

# Modeling uses DAILY means; the hourly series is retained for the daily-cycle
# diagnostics. Aggregating is a modeling choice, stated once, not a filter
# applied silently to the series being modeled.
pm_daily <- pm |>
  mutate(date = as.Date(datetime)) |>
  summarise(value = mean(value_imp), .by = date) |>
  arrange(date)

series <- if (PM25_RESOLUTION == "daily") {
  ts(pm_daily$value, frequency = 7)      # weekly period for the daily series
} else {
  ts(pm$value_imp, frequency = 24)       # daily period for the hourly series
}
c(resolution = PM25_RESOLUTION, n = length(series),
  frequency = frequency(series),
  mean = round(mean(series), 2), sd = round(sd(series), 2))

#> resolution          n  frequency       mean         sd 
#>    "daily"     "2890"        "7"     "93.7"    "73.95"

ggplot(pm_daily, aes(date, value)) +
  geom_line(linewidth = 0.3, color = "grey45") +
  geom_smooth(method = "loess", span = 0.15, se = FALSE,
              color = "#3B7DD8", linewidth = 0.9) +
  labs(title = "Beijing daily mean PM2.5, 2008-2016",
       subtitle = "Grey: daily means. Blue: a loess smooth, shown for ORIENTATION only -- the models below are fitted to the grey series",
       x = NULL, y = "PM2.5 AQI") +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly(pm_daily, x = ~date, y = ~value, type = "scatter", mode = "lines",
        name = "Daily mean", line = list(width = 1)) |>
  layout(title = "Beijing daily mean PM2.5",
         xaxis = list(title = "Date"), yaxis = list(title = "PM2.5 AQI"),
         hovermode = "x unified")

The loess curve is labelled as orientation, and the models below are fitted to the unsmoothed series. §12.15 shows what happens when that distinction is not maintained.

9.2 The daily cycle

# Hour extracted from the TIMESTAMP, not from row position -- with 6.5% of
# hours missing, a positional modular index would mis-assign every hour after
# the first gap.
hourly_profile <- pm |>
  mutate(hour = lubridate::hour(datetime)) |>
  summarise(mean = mean(value_imp),
            q20 = quantile(value_imp, 0.20),
            q80 = quantile(value_imp, 0.80),
            n = dplyr::n(), .by = hour) |>
  arrange(hour)
head(hourly_profile, 4)

ggplot(hourly_profile, aes(hour, mean)) +
  geom_ribbon(aes(ymin = q20, ymax = q80), fill = "grey86") +
  geom_line(linewidth = 1, color = "steelblue") +
  geom_point(size = 1.8) +
  scale_x_continuous(breaks = seq(0, 23, 3)) +
  labs(title = "Diurnal PM2.5 profile, averaged over eight years",
       subtitle = "Ribbon: 20th to 80th percentile. Hours taken from the timestamp, so gaps do not shift the alignment",
       x = "Hour of day (local)", y = "PM2.5 AQI") +
  theme_dspa()

10 Differencing and integration

The “I” in ARIMA is integration order \(d\): the number of differences needed to reach stationarity.

\[\nabla X_t=(1-B)X_t=X_t-X_{t-1}, \qquad \nabla^2X_t=X_t-2X_{t-1}+X_{t-2}\]

An ARIMA(\(p,d,q\)) is an ARMA(\(p,q\)) fitted to \(\nabla^dX_t\):

\[\phi(B)\,(1-B)^d X_t=\theta(B)\,\varepsilon_t.\]

Over-differencing has a signature and a cost. Differencing a series that is already stationary introduces a unit root into the MA polynomial: if \(X_t\) is white noise, \(\nabla X_t\) has \(\rho(1)=-0.5\) and is non-invertible. The tell-tale signs are a large negative lag-1 autocorrelation and an inflated residual variance.

Difference the minimum number of times the tests require, usually \(d\le2\), almost never more.

set.seed(51)
wn <- rnorm(500)
data.frame(
  series = c("white noise", "differenced once", "differenced twice"),
  acf_lag1 = round(c(acf_at(wn, 1), acf_at(diff(wn), 1), acf_at(diff(wn, differences = 2), 1)), 4),
  variance = round(c(var(wn), var(diff(wn)), var(diff(wn, differences = 2))), 4),
  ndiffs_says = c(ndiffs(wn), ndiffs(diff(wn)), ndiffs(diff(wn, differences = 2))))

Each unnecessary difference drives \(\rho(1)\) toward \(-0.5\) and roughly doubles the variance. ndiffs() correctly says zero differences are needed in all three cases.

d_needed <- forecast::ndiffs(series)
D_needed <- forecast::nsdiffs(series)
c(regular_differences = d_needed, seasonal_differences = D_needed,
  adf_p = round(suppressWarnings(tseries::adf.test(series)$p.value), 4),
  kpss_p = round(suppressWarnings(tseries::kpss.test(series)$p.value), 4))

#>  regular_differences seasonal_differences                adf_p 
#>               0.0000               0.0000               0.0100 
#>               kpss_p 
#>               0.0624

plot_acf(as.numeric(series), lag.max = 40,
         "Beijing daily PM2.5: ACF and PACF of the unsmoothed series")

11 Seasonal ARIMA

Series with a repeating cycle of period \(m\) need seasonal terms. \(\text{SARIMA}(p,d,q)(P,D,Q)_m\) is

\[\boxed{\;\Phi(B^m)\,\phi(B)\,(1-B^m)^{D}(1-B)^{d}\,X_t=\Theta(B^m)\,\theta(B)\,\varepsilon_t\;}\]

with \(\Phi,\Theta\) polynomials in \(B^m\) acting across cycles while \(\phi,\theta\) act within them.

Common misconception: “just use a large MA order to capture the seasonality.” Fitting \(\text{ARIMA}(1,1,24)\) to hourly data to capture a daily cycle estimates 24 free MA parameters where \(\text{SARIMA}(1,1,1)(0,1,1)_{24}\) captures the same structure with three.

The seasonal parameterization is not merely more parsimonious. It encodes the right structure, that lag 24, 48, 72 are related to each other in a way lags 5, 11, 17 are not, and it generalizes to horizons beyond one cycle, which a long non-seasonal MA cannot.

# Fit on a training span; the final year is held out for honest evaluation.
n_ser <- length(series)
h_out <- if (PM25_RESOLUTION == "daily") 180 else 24 * 30
train <- window(series, end = time(series)[n_ser - h_out])
test  <- window(series, start = time(series)[n_ser - h_out + 1])
c(train_length = length(train), test_length = length(test))

#> train_length  test_length 
#>         2710          180

fit_auto <- forecast::auto.arima(train, seasonal = TRUE, stepwise = TRUE,
                                 approximation = TRUE)
fit_auto

#> Series: train 
#> ARIMA(3,1,2)(0,0,2)[7] 
#> 
#> Coefficients:
#>          ar1    ar2     ar3     ma1     ma2    sma1   sma2
#>       -0.273  0.417  -0.129  -0.048  -0.921  -0.033  0.029
#> s.e.   0.038  0.032   0.023   0.033   0.032   0.020  0.020
#> 
#> sigma^2 = 3366:  log likelihood = -14842.4
#> AIC=29700.8   AICc=29700.8   BIC=29748

fit_ns <- forecast::Arima(train, order = arimaorder(fit_auto)[1:3],
                          seasonal = c(0, 0, 0))
data.frame(
  model = c(sprintf("auto: %s", paste(arimaorder(fit_auto), collapse = ",")),
            "same orders, no seasonal terms"),
  n_params = c(length(coef(fit_auto)), length(coef(fit_ns))),
  AIC = round(c(AIC(fit_auto), AIC(fit_ns)), 2),
  BIC = round(c(BIC(fit_auto), BIC(fit_ns)), 2),
  sigma2 = round(c(fit_auto$sigma2, fit_ns$sigma2), 3))

12 Estimation and search cost

ARIMA parameters are estimated by maximum likelihood, and the likelihood is evaluated by putting the model in state-space form and running the Kalman filter:

\[\ell(\phi,\theta,\sigma^2)=-\frac12\sum_{t=1}^{T}\left[\log(2\pi F_t)+\frac{v_t^2}{F_t}\right],\]

where \(v_t\) is the one-step prediction error and \(F_t\) its variance, both produced by the filter recursion. This costs \(O(T)\) per likelihood evaluation for fixed \((p,q)\), and it handles missing values exactly, the filter simply skips the update step, which is why state-space methods are the right tool for gappy series.

Task	Cost
One likelihood evaluation (Kalman)	\(O(T(p+q+m)^2)\)
One model fit (optimizer, \(I\) iterations)	\(O(I\,T(p+q+m)^2)\)
`auto.arima` stepwise search	\(\approx O\big((p_{\max}+q_{\max})\big)\) fits
`auto.arima` exhaustive (`stepwise = FALSE`)	\(O\big(p_{\max}q_{\max}P_{\max}Q_{\max}\big)\) fits
`stl` decomposition	\(O(T)\) per inner iteration

The stepwise search is the default for a reason. An exhaustive search over \(p,q\le5\) and \(P,Q\le2\) is \(5\times5\times3\times3\times2=450\) fits; the stepwise algorithm typically evaluates a few dozen. It can miss the global optimum, and stepwise = FALSE is worth the cost on a final model.

13 Diagnostics

If the orders are right, the residuals should be indistinguishable from white noise. The omnibus test is Ljung–Box:

\[Q(m)=T(T+2)\sum_{h=1}^{m}\frac{\hat\rho_\varepsilon(h)^2}{T-h}\ \overset{H_0}{\sim}\ \chi^2_{m-p-q}\]

The degrees of freedom must be reduced by the number of fitted ARMA parameters. Testing residuals with \(m\) degrees of freedom instead of \(m-p-q\) makes the test conservative, it under-rejects, and so fails to flag models that are genuinely inadequate. Box.test(..., fitdf = p + q) supplies the correction.

res <- residuals(fit_auto)
n_par <- length(coef(fit_auto))

lb <- do.call(rbind, lapply(c(10, 20, 30), function(m) {
  bt <- Box.test(res, lag = m, type = "Ljung-Box", fitdf = n_par)
  data.frame(lag = m, statistic = round(bt$statistic, 2),
             df = bt$parameter, p_value = round(bt$p.value, 4),
             residuals_look_white = bt$p.value > 0.05)
}))
lb

plot_acf(as.numeric(res), 36,
         sprintf("Residual ACF/PACF: ARIMA(%s)", paste(arimaorder(fit_auto), collapse = ",")))

rdf <- data.frame(t = seq_along(res), r = as.numeric(res))
p_ts <- ggplot(rdf, aes(t, r)) +
  geom_hline(yintercept = 0, color = "grey60") +
  geom_line(linewidth = 0.3, color = "steelblue") +
  labs(title = "Residuals over time", x = NULL, y = NULL) + theme_dspa(9)
p_qq <- ggplot(rdf, aes(sample = r)) +
  stat_qq(size = 0.5, alpha = 0.5) + stat_qq_line(color = "firebrick") +
  labs(title = "Normal Q-Q", x = NULL, y = NULL) + theme_dspa(9)
p_hi <- ggplot(rdf, aes(r)) +
  geom_histogram(aes(y = after_stat(density)), bins = 50,
                 fill = "steelblue", color = "white") +
  stat_function(fun = dnorm, args = list(mean = 0, sd = sd(rdf$r)),
                color = "firebrick", linewidth = 0.8) +
  labs(title = "Residual distribution vs. Normal", x = NULL, y = NULL) + theme_dspa(9)

(p_ts / (p_qq | p_hi)) +
  patchwork::plot_annotation(
    title = "ARIMA residual diagnostics",
    subtitle = "Heavy tails are common in air-quality data; they widen true prediction intervals beyond the Gaussian ones")

# --- Interactive equivalent ------------------------------------------------
plot_ly(rdf, x = ~t, y = ~r, type = "scatter", mode = "lines",
        name = "Residuals") |>
  add_lines(x = ~range(t), y = rep(1.96 * sd(rdf$r), 2), name = "+1.96 SD",
            line = list(dash = "dash")) |>
  add_lines(x = ~range(t),

DSPA3: Data Science and Predictive Analytics (UMich HS650)

Chapter 12: Longitudinal Data, Time Series, and Forecasting

SOCR/MCAIM (Ivo Dinov)

August 28, 2026

How this chapter uses graphics