---
title: "CSSS/POLS 512 - Time Series and Panel Data Methods"
subtitle: "Lab 2: Time Series Diagnostics — Box-Jenkins Identification"
author: "Ramses Llobet"
date: "Spring 2026"
output:
  html_document:
    theme: flatly
    toc: true
    toc_depth: 4
    toc_float: true
    code_folding: show
    fig_width: 8
    fig_height: 3.7
  pdf_document:
    toc: true
    toc_depth: 4
    number_sections: true
    fig_width: 7
    fig_height: 3.4
    latex_engine: xelatex
linkcolor: blue
urlcolor: blue
fontsize: 11pt
header-includes:
  - \usepackage{booktabs}
  - \usepackage{amsmath}
  - \usepackage{amssymb}
  - \usepackage{xcolor}
  - \usepackage{float}
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      error = FALSE,
                      warning = FALSE,
                      message = FALSE,
                      fig.align = 'center',
                      fig.width = 8,
                      fig.height = 8 / 1.618)

# clean environment
rm(list = ls())

# load packages
library(forecast)    # Acf, Pacf, auto.arima, Arima, checkresiduals
library(tseries)     # adf.test, kpss.test, PP.test, jarque.bera.test
library(patchwork)   # composing multi-panel figures
library(tidyverse)   # load last to avoid conflicts

set.seed(512)

# Helper: choose htmlreg or texreg based on output format
# This ensures regression tables render correctly in both HTML and PDF
show_reg <- function(...) {
  if (knitr::is_html_output()) {
    texreg::htmlreg(...)
  } else {
    texreg::texreg(...)
  }
}
```


# Disclaimer

These materials have been revised and updated by successive teaching assistants over recent years. 

**Revision history:**

- Ramses Llobet (Spring 2026); [rllobet@uw.edu](mailto:rllobet@uw.edu)



# Packages and Functions Used in This Lab

This lab uses four packages. Load them at the start of your session (the setup chunk above does this automatically).

| Package | Purpose | Key functions |
|---------|---------|---------------|
| `forecast` | Time series modeling and diagnostics | `Arima()`, `auto.arima()`, `checkresiduals()`, `Acf()`, `Pacf()`, `autoplot()` |
| `tseries` | Stationarity and specification tests | `adf.test()`, `kpss.test()`, `PP.test()`, `jarque.bera.test()` |
| `patchwork` | Composing multi-panel ggplot figures | `+`, `/` operators to arrange plots |
| `tidyverse` | Data manipulation and plotting | `ggplot()`, `tibble()`, `mutate()`, pipes (`|>`) |

Several tasks can be done with either base R or the `forecast` package. We use the `forecast` versions throughout for cleaner output, but both are valid:

| Task | Base R | `forecast` |
|------|--------|------------|
| Plot a series | `plot(ts)` | `autoplot(ts)` |
| ACF / PACF | `acf(x)` / `pacf(x)` | `Acf(x)` / `Pacf(x)` (drops lag 0) |
| Estimate ARIMA | `arima(x, order)` | `Arima(x, order)` (better constant handling) |
| Auto model selection | — | `auto.arima(x)` |
| Residual diagnostics | `Box.test()` + manual plots | `checkresiduals(fit)` (all-in-one) |
| Decomposition | `stl(ts, s.window)` | `mstl(ts)` (multiple seasonal periods) |


# Part 1: Time Series Building Blocks

In cross-sectional regression, we typically assume that observations are independent of each other. Time series data violate this assumption fundamentally: because observations are ordered in time, the value of a variable today is often related to its own past values and to past shocks. This creates **serial correlation** — a pattern of dependence across time that, if ignored, leads to biased standard errors, inefficient estimates, and misleading inference.

The **Box-Jenkins method** provides a systematic framework for identifying, estimating, and validating the dynamic structure of a time series. The core idea is that most time series can be decomposed into a combination of recognizable components:

$$
y_t =  \underbrace{\beta_0}_\text{level} + \underbrace{\beta_1 t}_\text{trend} + \underbrace{S_t}_\text{seasonal} + \underbrace{\phi_1 y_{t-1}}_\text{AR(1)} + \underbrace{\theta_1 \varepsilon_{t-1}}_\text{MA(1)} + \underbrace{\varepsilon_t}_\text{white noise}
$$

Not every series contains all of these components — the goal of diagnostics is to determine which ones are present and how strong they are. Once we correctly specify the dynamic form, the model's residuals $\hat{e}_t$ should resemble **white noise**: uncorrelated, zero-mean, constant-variance shocks with no remaining temporal structure. If they don't, we revise and re-estimate — this iterative loop is the essence of Box-Jenkins.

This part of the lab introduces each component in isolation so that you learn to recognize its signature. Part 2 addresses complications that arise in practice (trends, unit roots, seasonality). Part 3 covers estimation, diagnostics, and model selection. Part 4 puts it all together with practice exercises.


## 1.1 What is a time series?

A **time series** is a sequence of observations indexed by time. The ordering matters: unlike cross-sectional data, where each row is exchangeable, the position of each observation in a time series carries information. The value of GDP in Q3 depends on GDP in Q2; today's approval rating is related to yesterday's. This temporal ordering is what distinguishes time series from cross-sectional analysis and is the source of both the challenges and the opportunities in this course.

In R, the `ts()` function creates a time series object that stores metadata about the frequency and time index alongside the data.

```{r ts-basics}
# Create a ts object from a vector
# frequency: number of observations per unit of time (1 = annual, 4 = quarterly, 12 = monthly)
y_annual   <- ts(rnorm(30), start = 1990, frequency = 1)
y_quarterly <- ts(rnorm(80), start = c(2000, 1), frequency = 4)
y_monthly  <- ts(rnorm(120), start = c(2010, 1), frequency = 12)

# Inspect the time series attributes
# Annual series:
tsp(y_annual)    # start, end, frequency

# Quarterly series:
tsp(y_quarterly)
```

Two ways to plot a time series:

```{r ts-plot, fig.height=4}
# Base R: plot() recognizes ts objects and formats the x-axis automatically
plot(y_quarterly, main = "Base R: plot()", ylab = "Value", col = "steelblue")

# forecast package: autoplot() produces a ggplot2 object — easier to customize
autoplot(y_quarterly) +
  labs(title = "forecast::autoplot()", y = "Value") +
  theme_minimal(base_size = 13)
```



## 1.2 What to look for in a time series

Because observations are ordered in time, the first question with any time series is: **what kind of temporal structure does this series have?** The answer determines how we model it, what estimators are valid, and what inferences we can draw. Getting this wrong — for example, ignoring a trend or missing an autoregressive component — leads to biased coefficients, spurious significance, or poor predictions.

The Box-Jenkins method provides a practical ordering for this diagnostic work:

1. **Trend** — Is the series drifting upward or downward over time? If so, is the trend deterministic (a straight line we can remove) or stochastic (a unit root that requires differencing)?
2. **Seasonality** — Are there regular, repeating patterns at fixed intervals (monthly, quarterly)?
3. **Stationarity** — After removing trend and seasonality, does the series have a constant mean and variance over time? If not, we need to transform it before proceeding.
4. **Dynamic dependence** — Does the (stationary) series exhibit autoregressive (AR) or moving average (MA) behavior? This is where the ACF and PACF become central.
5. **Residual diagnostics** — After fitting a model, do the residuals look like white noise? If not, return to step 4.

This progression — from visual inspection to formal testing to model specification to residual checks — is not a one-shot procedure. It is an iterative loop that requires judgment and practice. The rest of this lab builds the tools you need at each step.


## 1.3 Decomposing trends and seasonality (preview)

Before we can diagnose the dynamic structure of a series (i.e., its AR/MA components), we often need to separate out two "large-scale" features — **trend** and **seasonality** — that can mask or distort the patterns we care about. The idea is that any time series can be written as an additive combination of three components:

$$
y_t = \underbrace{T_t}_{\text{trend}} + \underbrace{S_t}_{\text{seasonal}} + \underbrace{R_t}_{\text{remainder}}
$$

where the remainder $R_t$ is what we ultimately analyze for dynamic dependence. **STL** (*Seasonal-Trend decomposition using Loess*) is a standard tool for performing this decomposition. Let's see what it looks like on a series that actually has trend and seasonality built in:

```{r decompose-demo}
# Simulate a series with all three components:
# trend (upward drift) + seasonality (quarterly cycle) + AR(1) noise
set.seed(512)
n_demo <- 80  # 20 years of quarterly data

trend_demo    <- 0.15 * (1:n_demo)                           # linear trend
seasonal_demo <- rep(c(2, -1, -1.5, 0.5), length.out = n_demo) # quarterly pattern
ar_demo       <- arima.sim(model = list(ar = 0.5), n = n_demo)  # AR(1) noise

y_demo <- ts(as.numeric(trend_demo + seasonal_demo + ar_demo),
             start = c(2000, 1), frequency = 4)
```

```{r decompose-demo-plot, fig.height=3.5}
autoplot(y_demo) +
  labs(title = "Quarterly Series: Trend + Seasonality + AR(1)", y = expression(y[t])) +
  theme_minimal(base_size = 13)
```

The upward drift (trend) and the repeating quarterly oscillation (seasonality) are both visible. Now decompose:

```{r decompose-stl, fig.height=6}
# stl() requires a ts object with frequency > 1
# s.window = "periodic" assumes the seasonal pattern is the same every year
stl_demo <- stl(y_demo, s.window = "periodic")
plot(stl_demo, main = "STL Decomposition")
```

STL has separated the three components. The **remainder** panel at the bottom is what we would analyze for AR/MA structure — it should reflect the AR(1) process we built in. In Section 2.3 we will work through this workflow in detail, including the important question of *how* to remove seasonality properly for Box-Jenkins diagnostics.

```{r cleanup-ts-basics, include=FALSE}
rm(y_annual, y_quarterly, y_monthly, y_demo, trend_demo, seasonal_demo,
   ar_demo, n_demo, stl_demo)
```


## 1.4 White noise

White noise is the simplest time series process: a sequence of uncorrelated, identically distributed random variables with mean zero and constant variance:

$$
\varepsilon_t \sim WN(0, \sigma^2), \quad \text{Cov}(\varepsilon_t, \varepsilon_s) = 0 \text{ for } t \neq s
$$

Three properties define white noise: (1) the mean is zero at every time point, (2) the variance $\sigma^2$ is the same at every time point, and (3) there is no correlation between any two observations at different times. If your data satisfy these three properties, there is no temporal structure left to model.

This is our **null model** — what "no structure" looks like. Every diagnostic we learn in this lab is essentially asking: *does this series look different from white noise?* When we fit an ARMA model and check its residuals, we are asking whether the model has successfully captured all the temporal structure, leaving behind only white noise. If the residuals still show patterns, the model is incomplete.

```{r wn-simulate}
n <- 200
wn <- ts(rnorm(n, mean = 0, sd = 1))
```

```{r wn-plot, fig.height=3.5}
autoplot(wn) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  labs(title = "White Noise Process", y = expression(epsilon[t])) +
  theme_minimal(base_size = 13)

ggsave("output/white_noise_series.pdf", width = 8, height = 3.5)
```

### ACF and PACF

The **autocorrelation function (ACF)** measures how correlated a series is with its own past values at each lag $k$. Think of it as the Pearson correlation, but instead of measuring the association between two different variables, we measure the association between the *same* variable at two different points in time. The ACF at lag $k$ is the correlation between $y_t$ and $y_{t-k}$.

The **partial autocorrelation function (PACF)** measures the *direct* correlation between $y_t$ and $y_{t-k}$ after removing the indirect correlation transmitted through intermediate lags $y_{t-1}, \dots, y_{t-k+1}$. If you are familiar with partial regression coefficients from multiple regression, the logic is the same: the PACF at lag $k$ is the additional information in $y_{t-k}$ beyond what is already explained by lags $1$ through $k-1$.

This distinction is what makes Box-Jenkins identification possible: AR and MA processes leave different *signatures* in the ACF vs. PACF, allowing us to tell them apart. The formal definitions and sample estimation formulas are in Appendix B.

### ACF and PACF of white noise

For white noise, both the ACF and PACF should be zero at all lags $k \geq 1$ (since there is no correlation between observations at different times).

**The significance bands.** The blue dashed lines in ACF/PACF plots are approximate 95% confidence bands under the null hypothesis of white noise. For the ACF, these are computed as $\pm 1.96 / \sqrt{n}$, where $n$ is the length of the series. The logic: under $H_0$ (no autocorrelation), each sample autocorrelation $r_k$ is approximately normally distributed with mean 0 and standard error $1/\sqrt{n}$. The $\pm 1.96$ gives us the 95% critical value.

This means the **width of the bands depends entirely on the sample size**:

| $n$ | Band width ($\pm 1.96/\sqrt{n}$) | Minimum detectable $|r_k|$ |
|-----|----------------------------------|---------------------------|
| 50 | $\pm 0.277$ | Correlations below 0.28 are "insignificant" |
| 100 | $\pm 0.196$ | Threshold drops to 0.20 |
| 200 | $\pm 0.139$ | Threshold drops to 0.14 |
| 500 | $\pm 0.088$ | Even small correlations become detectable |
| 1000 | $\pm 0.062$ | Very narrow bands |

The practical implication: with short series ($n < 100$), only strong autocorrelation will be detected — you may fail to identify a weak AR or MA component. With long series ($n > 500$), even trivially small correlations become "significant," and you may see spikes that are statistically significant but substantively negligible. In either case, use judgment: a single borderline spike in an otherwise clean correlogram is likely noise (we expect about 5% of lags to exceed the bands by chance), while a systematic pattern of decaying spikes is meaningful regardless of whether each individual bar crosses the line.

R provides two interfaces for ACF/PACF plots: the base R functions `acf()`/`pacf()` (include lag 0, base graphics) and the `forecast` package functions `ggAcf()`/`ggPacf()` (drop lag 0, ggplot2 output). We use the `forecast` versions throughout this lab for cleaner output.

```{r wn-acf-forecast, fig.height=4}
p_acf  <- ggAcf(wn) + labs(title = "ACF — White Noise") + theme_minimal(base_size = 12)
p_pacf <- ggPacf(wn) + labs(title = "PACF — White Noise") + theme_minimal(base_size = 12)

p_acf + p_pacf

ggsave("output/white_noise_acf_pacf.pdf", width = 8, height = 3.5)
```

All bars sit within the confidence bands — as expected for a series with no autocorrelation.

### Ljung-Box test

The ACF and PACF plots let us inspect individual lags visually, but sometimes we want a single, formal test of whether a series has *any* autocorrelation across multiple lags simultaneously. The **Ljung-Box test** does exactly this.

- $H_0$: The first $h$ autocorrelations are all zero (the data are white noise)
- $H_1$: At least one of the first $h$ autocorrelations is nonzero (serial correlation exists)

The test works by summing the squared sample autocorrelations across the first $h$ lags. If the series is truly white noise, all autocorrelations should be near zero, so the test statistic will be small. If any autocorrelation is present, the statistic grows and we reject. Always use `type = "Ljung-Box"` (preferred over the older Box-Pierce variant for its better small-sample performance). The test statistic formula and lag selection guidelines are in Appendix B.

```{r wn-ljungbox}
Box.test(wn, lag = 10, type = "Ljung-Box")
```

The $p$-value is large ($> 0.05$), so we fail to reject the null of white noise. This is correct — we simulated white noise.

**Why this test matters for regression users:** If you run an OLS regression and the Ljung-Box test on the residuals fails to reject, your error term is consistent with white noise — meaning OLS standard errors, $t$-tests, and confidence intervals are valid. If the test *rejects*, there is autocorrelation in the residuals, which means your standard errors are **too small** (you think you have more independent information than you actually do) and you are over-rejecting null hypotheses. In that case, you need to either model the dynamic structure of the errors (ARMA errors, which we cover in Lab 3) or use autocorrelation-robust standard errors (e.g., Newey-West). The Ljung-Box test is the bridge between time series diagnostics and the regression inference you already know.

```{r cleanup-wn, include=FALSE}
rm(wn, p_acf, p_pacf)
```


## 1.5 Autoregressive processes

An **autoregressive** (AR) process generates the current value as a linear function of its own past values plus a white noise shock. The logic is straightforward: what happened yesterday influences what happens today. This is perhaps the most natural form of temporal dependence in the social sciences — a country's GDP this quarter depends on last quarter's, today's presidential approval is related to yesterday's.

### AR(1)

The simplest autoregressive process is AR(1):

$$
y_t = \phi \, y_{t-1} + \varepsilon_t, \quad \varepsilon_t \sim WN(0, \sigma^2)
$$

The parameter $\phi$ controls the **persistence**: how strongly the current value depends on the previous one. To build intuition, we can substitute the AR(1) equation into itself repeatedly (recursive substitution):

$$
y_t = \phi \, y_{t-1} + \varepsilon_t = \phi(\phi \, y_{t-2} + \varepsilon_{t-1}) + \varepsilon_t = \phi^2 y_{t-2} + \phi \varepsilon_{t-1} + \varepsilon_t = \cdots
$$

Continuing this indefinitely, we get the **moving average representation** of the AR(1):

$$
y_t = \sum_{j=0}^{\infty} \phi^j \, \varepsilon_{t-j} = \varepsilon_t + \phi \, \varepsilon_{t-1} + \phi^2 \, \varepsilon_{t-2} + \phi^3 \, \varepsilon_{t-3} + \cdots
$$

This reveals something important: an AR(1) process is the sum of *all* past shocks, each weighted by $\phi^j$. The current value of $y$ "remembers" every random shock back to the beginning of time. However, when $|\phi| < 1$, the weights $\phi^j$ shrink geometrically — each past shock matters less than the one before it. A shock that occurred 10 periods ago has weight $\phi^{10}$; if $\phi = 0.8$, that weight is $0.8^{10} \approx 0.11$, so only about 11% of the shock remains. If $\phi = 0.3$, the weight after 10 periods is $0.3^{10} \approx 0.000006$ — essentially gone.

This is the **stationarity condition**: the process is stationary when $|\phi| < 1$, because past shocks decay and the series fluctuates around a fixed mean. When $|\phi| = 1$ (unit root / random walk), the weights never decay — every shock permanently shifts the level of the series. When $|\phi| > 1$, the process is explosive. In this course, we focus on the stationary case ($|\phi| < 1$) and on detecting when the unit root case ($|\phi| = 1$) applies.

### What $\phi$ tells you: the half-life of shocks

A useful way to interpret $\phi$ substantively is through the **half-life** of a shock: how many periods it takes for a perturbation to lose half its impact. For an AR(1), the half-life is $t_{1/2} = -\ln(2) / \ln(|\phi|)$:

| $\phi$ | Half-life | Interpretation |
|--------|-----------|----------------|
| 0.3 | 0.6 periods | Shock virtually gone in 1 period |
| 0.7 | 1.9 periods | Shock halves roughly every 2 periods |
| 0.9 | 6.6 periods | Shock persists for many periods |
| 0.95 | 13.5 periods | Very persistent — hard to distinguish from unit root |
| 1.0 | $\infty$ | Shock **never** dissipates |

### Why $\phi \geq 1$ breaks estimation: the variance explosion

The variance of a stationary AR(1) process is (see Appendix B for derivation):

$$
\text{Var}(y_t) = \frac{\sigma^2_\varepsilon}{1 - \phi^2}
$$

This formula reveals the fundamental problem with unit roots. As $\phi \to 1$, the denominator $(1 - \phi^2) \to 0$ and the variance $\to \infty$. The variance of the series is not a fixed number — it grows without bound because shocks accumulate without dissipating. A process whose variance is infinite cannot have a stable mean or a well-defined autocorrelation function, so the usual estimation and inference tools break down. This is why we test for unit roots (Section 2.2) before attempting to estimate ARMA models.

```{r variance-explosion, fig.height=3.5}
# Visualize how the variance of an AR(1) grows as φ → 1
phi_seq <- seq(0, 0.99, by = 0.01)
var_ratio <- 1 / (1 - phi_seq^2)  # Var(y) / sigma^2

plot(phi_seq, var_ratio, type = "l", lwd = 2, col = "steelblue",
     xlab = expression(phi), ylab = expression(Var(y[t]) / sigma[epsilon]^2),
     main = expression("Variance of AR(1) as " * phi %->% 1),
     ylim = c(0, 30), las = 1)
abline(v = 1, lty = 2, col = "firebrick")
text(0.92, 25, expression(phi == 1: ~ Var %->% infinity),
     col = "firebrick", cex = 1.1, adj = 1)

ggsave("output/ar1_variance_explosion.pdf", width = 8, height = 3.5)
```

The plot makes the point visually: the variance is modest for small $\phi$ but accelerates sharply as $\phi$ approaches 1. At $\phi = 0.9$, the variance is already about 5 times the innovation variance $\sigma^2_\varepsilon$; at $\phi = 0.99$, it is about 50 times; at $\phi = 1$, it is undefined.

```{r ar1-simulate}
ar1 <- arima.sim(model = list(ar = 0.85), n = 200)
```

```{r ar1-plot, fig.height=3.5}
autoplot(ar1) +
  labs(title = "AR(1) with \u03c6 = 0.85", y = expression(y[t])) +
  theme_minimal(base_size = 13)
```

Notice the smooth, persistent patterns — the series "remembers" where it has been. Compared to white noise, the series wanders more slowly and stays above or below zero for extended stretches. This is the visual hallmark of strong positive autocorrelation.

```{r ar1-acf, fig.height=4}
p1 <- ggAcf(ar1) + labs(title = "ACF — AR(1)") + theme_minimal(base_size = 12)
p2 <- ggPacf(ar1) + labs(title = "PACF — AR(1)") + theme_minimal(base_size = 12)

p1 + p2

ggsave("output/ar1_acf_pacf.pdf", width = 8, height = 3.5)
```

**Key pattern:** The ACF **tails off** — it decays geometrically but never drops abruptly to zero. This is because the AR(1) process generates correlation at *all* lags (recall the infinite sum representation: $y_t$ contains traces of every past shock). The theoretical ACF of an AR(1) is $\rho_k = \phi^k$, so at lag 1 the ACF equals $\phi$, at lag 2 it equals $\phi^2$, and so on — a smooth geometric decay.

The PACF, in contrast, **cuts off sharply after lag 1**. This is because once we control for $y_{t-1}$, there is no *additional* direct information in $y_{t-2}$, $y_{t-3}$, etc. — all their influence on $y_t$ is already captured through $y_{t-1}$. The PACF at lag 1 gives a direct estimate of $\phi$ — here it should be close to 0.85.

This "tailing off ACF + sharp PACF cutoff" is **the** signature of an AR process. The lag at which the PACF cuts off tells you the order $p$.

### Understanding the PACF: a regression perspective

The PACF can feel abstract — what does "partialing out intermediate lags" actually mean? If you understand multiple regression, you already have the intuition. The PACF at lag $k$ is simply the **coefficient on $y_{t-k}$** when you regress $y_t$ on all lags from 1 through $k$.

Let's see this directly. Using the AR(1) series we just simulated, we fit a sequence of regressions, each adding one more lag:

```{r pacf-regression-demo}
# Build a data frame with the series and its lags
ar1_df <- tibble(
  y    = as.numeric(ar1),
  lag1 = dplyr::lag(y, 1),
  lag2 = dplyr::lag(y, 2),
  lag3 = dplyr::lag(y, 3),
  lag4 = dplyr::lag(y, 4),
  lag5 = dplyr::lag(y, 5)
) |> na.omit()

# Fit sequential regressions — the coefficient on the LAST lag is the PACF
reg1 <- lm(y ~ lag1, data = ar1_df)
reg2 <- lm(y ~ lag1 + lag2, data = ar1_df)
reg3 <- lm(y ~ lag1 + lag2 + lag3, data = ar1_df)
reg4 <- lm(y ~ lag1 + lag2 + lag3 + lag4, data = ar1_df)
reg5 <- lm(y ~ lag1 + lag2 + lag3 + lag4 + lag5, data = ar1_df)

# Compare: coefficient on the last lag added ≈ PACF at that lag
pacf_from_regression <- c(
  coef(reg1)["lag1"],
  coef(reg2)["lag2"],
  coef(reg3)["lag3"],
  coef(reg4)["lag4"],
  coef(reg5)["lag5"]
)

# Compare to R's PACF function
pacf_from_r <- pacf(ar1, lag.max = 5, plot = FALSE)$acf[1:5]

tibble(
  Lag = 1:5,
  `PACF (from lm)` = round(pacf_from_regression, 3),
  `PACF (from pacf())` = round(pacf_from_r, 3)
)
```

The values match. The PACF at lag 1 ($\approx 0.85$) is the coefficient on $y_{t-1}$ in a simple regression — this is just the AR(1) parameter $\phi$. But when we add $y_{t-2}$ to the regression, its coefficient is near zero — it has nothing to add once $y_{t-1}$ is already included. The same is true for lags 3, 4, 5. This is why the PACF *cuts off after lag 1* for an AR(1): only the first lag has a direct relationship with $y_t$; all other lags' apparent correlation (visible in the ACF) is entirely **inherited** through the chain $y_t \leftarrow y_{t-1} \leftarrow y_{t-2} \leftarrow \cdots$.

### Direct vs. inherited correlation

This gives us a powerful way to understand the difference between ACF and PACF. For an AR(1) with $\phi = 0.85$:

- The **ACF** at lag 2 is $\phi^2 = 0.72$ — a large value. But this correlation is entirely *inherited* through lag 1. There is no direct link between $y_t$ and $y_{t-2}$.
- The **PACF** at lag 2 is $\approx 0$ — after accounting for $y_{t-1}$, lag 2 contributes nothing.

We can visualize this decomposition:

```{r direct-vs-inherited, fig.height=4.5}
# Decompose the ACF of an AR(1) into direct (PACF) and inherited (ACF - PACF)
n_lags <- 10
acf_vals  <- acf(ar1, lag.max = n_lags, plot = FALSE)$acf[-1]  # drop lag 0
pacf_vals <- pacf(ar1, lag.max = n_lags, plot = FALSE)$acf

decomp_df <- tibble(
  Lag = 1:n_lags,
  ACF = acf_vals,
  PACF = pacf_vals,
  Direct = ifelse(abs(PACF) > 1.96 / sqrt(length(ar1)), PACF, 0),
  Inherited = ACF - Direct
)

# Stacked bar chart
decomp_df |>
  pivot_longer(cols = c(Direct, Inherited),
               names_to = "Component", values_to = "Correlation") |>
  mutate(Component = factor(Component, levels = c("Inherited", "Direct"))) |>
  ggplot(aes(x = factor(Lag), y = Correlation, fill = Component)) +
  geom_col(position = "stack", width = 0.6) +
  geom_hline(yintercept = c(-1, 1) * 1.96 / sqrt(length(ar1)),
             linetype = "dashed", color = "steelblue", linewidth = 0.5) +
  scale_fill_manual(values = c("Inherited" = "grey70", "Direct" = "firebrick")) +
  labs(title = "AR(1): Decomposing ACF into Direct and Inherited Correlation",
       subtitle = "Direct (red) = PACF contribution; Inherited (grey) = passed through intermediate lags",
       x = "Lag", y = "Correlation", fill = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

ggsave("output/ar1_direct_vs_inherited.pdf", width = 8, height = 4.5)
```

Only lag 1 has a red (direct) component — this is the $\phi$ parameter. Every other lag's correlation is entirely grey (inherited through the chain of intermediate lags). This is *why* the PACF cuts off after lag 1: there is only one direct connection, and the rest is inherited.

**A note on spurious spikes:** You may notice that one or two bars at higher lags (e.g., lag 5 or lag 9) barely cross the significance band. This is expected — recall that the confidence bands are set at 95%, so roughly **1 in 20 lags will cross the band by pure chance** even in a true white noise series. A single isolated spike at a high lag, surrounded by insignificant lags, is almost certainly noise. What matters is the *overall pattern*: a systematic decay (tailing off) or a clean drop to zero (cutoff), not individual borderline bars.

**For an AR(2)**, you would see red bars at lags 1 *and* 2 (two direct connections), with everything else inherited. For an **MA(1)**, the picture is different: the ACF itself cuts off at lag 1 (no inheritance chain exists, since shocks don't propagate through past values), while the PACF tails off.

This decomposition makes the identification table feel inevitable rather than arbitrary: it's just asking "how many direct connections are there?" (read from the PACF cutoff for AR) or "how long do shocks persist?" (read from the ACF cutoff for MA).


### Varying $\phi$: how persistence shapes behavior

The AR(1) parameter $\phi$ controls both the persistence and the visual character of the series. Understanding how different values of $\phi$ affect the series is essential for visual diagnosis:

```{r ar1-grid, fig.height=6}
set.seed(512)
phis <- c(0.3, 0.85, -0.7, 0.95)
labels <- paste0("\u03c6 = ", phis)

# Simulate four AR(1) series with different phi values
ar_sims <- map(phis, ~ arima.sim(model = list(ar = .x), n = 200))

# Compute a shared y-axis range across all four series
y_range <- range(unlist(map(ar_sims, range)))

ar_plots <- map2(ar_sims, labels, function(y, lab) {
  autoplot(y) +
    labs(title = lab, y = expression(y[t])) +
    ylim(y_range) +
    theme_minimal(base_size = 11)
})

# Arrange in a 2x2 grid — shared y-axis makes amplitude differences visible
(ar_plots[[1]] + ar_plots[[2]]) / (ar_plots[[3]] + ar_plots[[4]])

ggsave("output/ar1_varying_phi.pdf", width = 8, height = 6)
```

- **Small positive $\phi$** (0.3): Weak persistence. Shocks decay quickly ($0.3^5 \approx 0.002$), so the series looks nearly like white noise — a slight "smoothing" effect is the only visible difference.
- **Large positive $\phi$** (0.85, 0.95): Strong persistence. Shocks decay slowly, producing smooth, wandering patterns. As $\phi \to 1$, the series becomes increasingly difficult to distinguish visually from a random walk — this is why formal stationarity tests (Section 2.2) are essential.
- **Negative $\phi$** ($-0.7$): Oscillating, choppy patterns — the visual opposite of the smooth positive-$\phi$ series. See the discussion below.

### A note on negative AR(1): why it's rare and what it means

The $\phi = -0.7$ panel in the grid above looks strikingly different from the others — choppy, erratic, crossing the mean constantly. This is because a negative $\phi$ means each observation is pulled in the **opposite direction** from the previous one: a positive deviation is systematically followed by a negative one. The ACF reflects this, alternating in sign: $\rho_1 = -0.7$, $\rho_2 = 0.49$, $\rho_3 = -0.34$, etc.

**Negative AR is rare in the social sciences.** Most social, economic, and political processes exhibit positive persistence — today's state tends to resemble yesterday's. Positive AR reflects the inertia of institutions, habits, and structural conditions. A negative AR requires a specific **overcorrection mechanism** that pushes the system *past* its mean, not just *back toward* it. This does arise in some settings — the "cobweb model" in agricultural economics (high prices → overplanting → price crash → underplanting → price spike), budget cycles where institutional rules force overcorrection, and some ecological predator-prey dynamics — but these are exceptions. In most applied social science work, expect positive $\phi$ between 0.3 and 0.95. If you encounter a series that crosses the mean far more frequently than white noise, consider a negative AR component.

### AR(2)

The AR(2) process extends the dependence to two lags:

$$
y_t = \phi_1 \, y_{t-1} + \phi_2 \, y_{t-2} + \varepsilon_t
$$

Now $y_t$ depends directly on both $y_{t-1}$ and $y_{t-2}$. The interplay between two coefficients can produce richer dynamics, including quasi-cyclical behavior when the characteristic roots of the process are complex.

```{r ar2-simulate}
ar2 <- arima.sim(model = list(ar = c(0.5, 0.3)), n = 200)
```

```{r ar2-acf, fig.height=4}
p1 <- ggAcf(ar2) + labs(title = "ACF — AR(2)") + theme_minimal(base_size = 12)
p2 <- ggPacf(ar2) + labs(title = "PACF — AR(2)") + theme_minimal(base_size = 12)

p1 + p2

ggsave("output/ar2_acf_pacf.pdf", width = 8, height = 3.5)
```

**Key pattern:** The ACF tails off (decays, possibly with oscillation depending on the roots). The PACF **cuts off after lag 2** — lags 1 and 2 are significant, but lag 3 and beyond are not. The PACF values at lags 1 and 2 approximate $\phi_1$ and $\phi_2$. The general rule holds: for an AR($p$) process, the PACF cuts off after lag $p$.

```{r cleanup-ar, include=FALSE}
rm(ar1, ar2, ar_plots, phis, labels, p1, p2)
```


## 1.6 Moving average processes

While AR processes model dependence on past *values* of the series (the system's own inertia), **moving average** (MA) processes model the system's **response to external shocks**. The current observation is a linear combination of the current and recent white noise innovations. The key difference is one of **memory**: an AR process has infinite memory (all past shocks contribute, though with declining weights), while an MA($q$) process has **finite memory** — only the last $q$ shocks matter, and their influence disappears completely after $q$ periods.

Think of the distinction substantively:

- **AR = persistence / inertia.** GDP this quarter depends on GDP last quarter. A country at war tends to stay at war. Presidential approval today is related to yesterday's — the system carries its own state forward.
- **MA = shock response.** A policy announcement moves financial markets for a few days, then the effect fades. A political scandal depresses approval ratings for $q$ months, then disappears entirely. The system responds to external perturbations, but doesn't carry its own state forward — it "forgets" the shock after $q$ periods.

This distinction has an important implication for regression that we will develop in Lab 3: **AR errors in a regression can bias your coefficient estimates** (because the lagged dependent variable is correlated with the error term), while **MA errors do not bias coefficients** — they only affect the *efficiency* of standard errors. In other words, ignoring AR dynamics can give you the wrong answer; ignoring MA dynamics gives you the right answer but with the wrong confidence interval. This is one reason why correctly identifying the dynamic structure matters for inference, not just for forecasting.

### MA(1)

The simplest moving average process is:

$$
y_t = \varepsilon_t + \theta \, \varepsilon_{t-1}
$$

Here, $y_t$ depends on the current shock $\varepsilon_t$ and the immediately preceding shock $\varepsilon_{t-1}$. The parameter $\theta$ governs how strongly the previous shock carries over. After one period, the effect of any given shock is gone — the MA(1) process has a memory of exactly one period. This makes MA processes inherently stationary for any finite order $q$: since only a fixed number of past shocks matter, the statistical properties of the series do not change over time.

### What "MA(1) equals AR(∞)" means in plain language

In Section 1.5 we showed that an AR(1) can be rewritten as an infinite sum of past shocks (the moving average representation). The reverse is also true: an MA(1) can be rewritten as a regression on infinitely many past values of $y$. Formally, an MA(1) with $|\theta| < 1$ is equivalent to:

$$y_t = (-\theta) \, y_{t-1} + (-\theta)^2 \, y_{t-2} + (-\theta)^3 \, y_{t-3} + \cdots + \varepsilon_t$$

What does this mean in practice? Suppose you have a series where shocks fade out after one period (MA(1) behavior). If you tried to capture this pattern using an AR model instead, you would need to include *infinitely many lags* of $y$, each with a geometrically shrinking coefficient. With 3 lags you'd get a rough approximation; with 10 lags, a better one; but you'd never capture it exactly. The MA(1) does it with a single parameter $\theta$.

This is why MA representations exist: they provide a **parsimonious** way to capture shock-response dynamics that would require an unwieldy number of AR terms. Conversely, AR representations are parsimonious for persistence dynamics that would need many MA terms. Choosing between AR and MA is not about which is "correct" — both can represent the same process — but about which gives a **simpler description** of the dynamics you observe.

```{r ma1-simulate}
ma1 <- arima.sim(model = list(ma = 0.9), n = 200)
```

```{r ma1-plot, fig.height=3.5}
autoplot(ma1) +
  labs(title = "MA(1) with \u03b8 = 0.9", y = expression(y[t])) +
  theme_minimal(base_size = 13)
```

```{r ma1-acf, fig.height=4}
p1 <- ggAcf(ma1) + labs(title = "ACF — MA(1)") + theme_minimal(base_size = 12)
p2 <- ggPacf(ma1) + labs(title = "PACF — MA(1)") + theme_minimal(base_size = 12)

p1 + p2

ggsave("output/ma1_acf_pacf.pdf", width = 8, height = 3.5)
```

**Key pattern:** The ACF **cuts off sharply after lag 1** — only the first lag is significant, and all subsequent lags are essentially zero. This is because an MA(1) process only "remembers" one past shock; the correlation between $y_t$ and $y_{t-2}$ is zero since they share no common shock. The PACF, in contrast, **tails off** (geometric or oscillating decay).

Notice that this is the **mirror image** of the AR(1) pattern, where the ACF tails off and the PACF cuts off. This duality is a core insight of Box-Jenkins identification: the ACF behavior of an AR process mirrors the PACF behavior of an MA process, and vice versa. Recognizing this symmetry is the key to telling AR and MA processes apart.

### MA(2)

The MA(2) extends the memory to two past shocks:

$$
y_t = \varepsilon_t + \theta_1 \, \varepsilon_{t-1} + \theta_2 \, \varepsilon_{t-2}
$$

Now $y_t$ is influenced by both the previous shock and the one before it. After two periods, any given shock has no further effect — the process has a memory of exactly two periods.

```{r ma2-simulate}
ma2 <- arima.sim(model = list(ma = c(0.8, -0.5)), n = 200)
```

```{r ma2-acf, fig.height=4}
p1 <- ggAcf(ma2) + labs(title = "ACF — MA(2)") + theme_minimal(base_size = 12)
p2 <- ggPacf(ma2) + labs(title = "PACF — MA(2)") + theme_minimal(base_size = 12)

p1 + p2
```

The ACF cuts off after lag 2 — lags 1 and 2 are significant, lag 3 and beyond are not. The PACF tails off. This is exactly the MA($q$) pattern from the identification table, with $q = 2$. Note the symmetry with AR(2): for AR(2), it was the PACF that cut off after lag 2 while the ACF tailed off. For MA(2), the roles are reversed.

```{r cleanup-ma, include=FALSE}
rm(ma1, ma2, p1, p2)
```


## 1.7 ARMA processes

When a series has both autoregressive and moving average components, we get an **ARMA** process. This is common in practice: many real-world series exhibit both the persistence of AR dynamics (current values depend on past values) and the shock-absorption of MA dynamics (current values respond to recent innovations). For ARMA(1,1):

$$
y_t = \phi \, y_{t-1} + \varepsilon_t + \theta \, \varepsilon_{t-1}
$$

The $\phi$ parameter captures the autoregressive persistence, while $\theta$ captures the moving average response to past shocks. Both operate simultaneously.

```{r arma11-simulate}
arma11 <- arima.sim(model = list(ar = 0.7, ma = 0.4), n = 200)
```

```{r arma11-plot, fig.height=3.5}
autoplot(arma11) +
  labs(title = "ARMA(1,1) with \u03c6 = 0.7, \u03b8 = 0.4", y = expression(y[t])) +
  theme_minimal(base_size = 13)
```

```{r arma11-acf, fig.height=4}
p1 <- ggAcf(arma11) + labs(title = "ACF — ARMA(1,1)") + theme_minimal(base_size = 12)
p2 <- ggPacf(arma11) + labs(title = "PACF — ARMA(1,1)") + theme_minimal(base_size = 12)

p1 + p2

ggsave("output/arma11_acf_pacf.pdf", width = 8, height = 3.5)
```

**Key pattern:** Both the ACF and the PACF **tail off** — neither cuts off cleanly. This is the signature of a mixed ARMA process. Intuitively, the AR component makes the PACF decay (rather than cut off, as it would in a pure MA), and the MA component makes the ACF decay (rather than cut off, as it would in a pure AR). When you see both tailing off, you cannot determine the exact orders $p$ and $q$ from the correlograms alone. Instead, you need to try several candidate $(p, q)$ combinations and compare them using information criteria (AIC/BIC), which we cover in Part 3.

### Why mixed ARMA processes are hard to identify

There is a deep mathematical reason why identification is challenging: **AR and MA processes are equivalent representations of each other.** Any stationary AR(1) can be rewritten as an infinite-order MA process ($y_t = \sum_{j=0}^{\infty} \phi^j \varepsilon_{t-j}$, which we saw in Section 1.5), and conversely, any invertible MA(1) can be rewritten as an infinite-order AR process. This is called **invertibility**, and it means that in finite samples — especially short ones — an AR process can mimic an MA process and vice versa (Box-Steffensmeier et al., 2014, pp. 33–34).

This equivalence has two practical implications. First, it explains why pure ARMA($p,q$) processes with both $p > 0$ and $q > 0$ are relatively uncommon in practice: many mixed specifications can be reduced to a more parsimonious pure AR or pure MA model. Second, it means that when the ACF and PACF are ambiguous, you should always **over-fit and under-fit** your candidate model — try adding and removing terms to see whether the extra parameters are justified — rather than relying solely on the correlograms.

```{r cleanup-arma, include=FALSE}
rm(arma11, p1, p2)
```


## 1.8 The identification table

The following table summarizes the diagnostic signatures we have just covered. This is the core reference for Box-Jenkins identification — you will return to it repeatedly throughout this course.

| Process | ACF pattern | PACF pattern | Identification logic |
|---------|------------|--------------|---------------------|
| White noise | No significant lags | No significant lags | No temporal structure to model |
| AR($p$) | Tails off (geometric/oscillating decay) | **Cuts off after lag $p$** | PACF tells you the order $p$ |
| MA($q$) | **Cuts off after lag $q$** | Tails off | ACF tells you the order $q$ |
| ARMA($p,q$) | Tails off | Tails off | Cannot determine orders from correlograms alone; use AIC |

**How to use it:** Plot the ACF and PACF of your (stationary) series. Ask two questions: (1) Does the ACF cut off or tail off? (2) Does the PACF cut off or tail off? The combination of answers points you to the right row of the table. If one cuts off and the other tails off, you have a pure AR or pure MA process, and the cutoff lag gives you the order. If both tail off, you have a mixed ARMA and need to estimate several candidate models and compare them with information criteria.

**Important caveat:** The table assumes the series is already stationary. If trends or seasonality are present, the ACF will be dominated by those features (typically very slow decay), making it impossible to read the AR/MA structure underneath. This is why the Box-Jenkins procedure starts with detrending and seasonal adjustment *before* examining the ACF and PACF for dynamic dependence.


---

\newpage


# Part 2: Common Complications

Part 1 introduced the building blocks of time series in their "clean" forms — pure AR, pure MA, or mixed ARMA processes with no trend or seasonality. In practice, real-world series almost always come with additional complications that must be addressed before we can apply the identification table. This part covers the three most common: deterministic trends, unit roots (stochastic trends), and seasonality. Learning to recognize and handle these features is essential because they distort the ACF and PACF, making it impossible to read the underlying dynamic structure without first removing them.


## 2.1 Deterministic trends

Many real-world series have trends — GDP grows over time, temperatures rise, populations increase. A **deterministic trend** is the simplest form: a systematic, predictable increase (or decrease) over time that can be captured by a linear function of $t$. In the context of a dynamic model:

$$
y_t = \alpha + \delta \, t + \phi \, y_{t-1} + \varepsilon_t
$$

Here, $\delta$ is the slope of the trend. For every one-period increase in $t$, the expected value of $y_t$ increases by $\delta$ (holding the autoregressive dynamics constant). The key feature of a deterministic trend is that it is a *fixed, predictable* function of time — once you know the slope $\delta$, you can account for the trend exactly by regressing $y$ on $t$ and working with the residuals.

```{r trend-simulate}
set.seed(512)
n <- 200
trend_component <- 2 + 0.05 * (1:n)
ar_component <- arima.sim(model = list(ar = 0.6), n = n)
y_trend <- ts(as.numeric(trend_component + ar_component))
```

```{r trend-plot, fig.height=3.5}
autoplot(y_trend) +
  labs(title = "AR(1) with Deterministic Trend", y = expression(y[t])) +
  theme_minimal(base_size = 13)
```

The upward trend is clear. Now look at what this does to the ACF:

```{r trend-acf, fig.height=4}
p_raw <- ggAcf(y_trend) +
  labs(title = "ACF — Original (with trend)") +
  theme_minimal(base_size = 12)
```

The ACF decays very slowly — it looks like a unit root or a very persistent AR process. But the series is actually stationary around a trend. The slow ACF decay is an **artifact of the trend**, not genuine persistence. This is one of the most common pitfalls in time series diagnostics: a trending series produces high autocorrelation at all lags simply because adjacent observations share a common upward (or downward) drift. If you try to read AR/MA structure from this ACF, you will be misled.

### Detrending

The fix is simple: regress the series on time and work with the residuals.

```{r detrend}
# Fit a linear trend
time_index <- 1:length(y_trend)
trend_fit <- lm(y_trend ~ time_index)

# Extract residuals — these are the detrended series
y_detrended <- ts(residuals(trend_fit))
```

```{r detrend-compare, fig.height=4}
p_detrended <- ggAcf(y_detrended) +
  labs(title = "ACF — After Detrending") +
  theme_minimal(base_size = 12)

p_raw + p_detrended

ggsave("output/trend_detrending_acf.pdf", width = 8, height = 3.5)
```

After removing the trend, the ACF shows the geometric decay characteristic of the underlying AR(1) process. The PACF confirms:

```{r detrend-pacf, fig.height=4}
ggPacf(y_detrended) +
  labs(title = "PACF — Detrended Series") +
  theme_minimal(base_size = 12)
```

The PACF cuts off after lag 1 with a value close to 0.6 — exactly the AR(1) structure we built in.

```{r cleanup-trend, include=FALSE}
rm(y_trend, y_detrended, trend_component, ar_component, time_index, trend_fit,
   p_raw, p_detrended, n)
```


## 2.2 Unit roots and stationarity tests

The distinction between a deterministic trend and a **stochastic trend** (unit root) is one of the most important in time series analysis, and one of the most consequential for inference.

A **deterministic trend** is a straight line added on top of a stationary process. The randomness fluctuates *around* a predictable path, and after removing the trend, the series is well-behaved. A **stochastic trend** (unit root) is fundamentally different — the process wanders without any tendency to return to a fixed level:

$$
y_t = y_{t-1} + \varepsilon_t \quad \text{(random walk)}
$$

This is an AR(1) with $\phi = 1$. Recall from Section 1.5 that the moving average representation of an AR(1) is $y_t = \sum_{j=0}^{\infty} \phi^j \varepsilon_{t-j}$. When $\phi = 1$, all the weights equal 1 — every past shock contributes with full weight to the current value. Each shock **permanently shifts the level** of the series. Unlike a stationary process, which reverts to its mean after a shock, a random walk has no fixed mean to revert to.

This distinction matters for inference because:

- **Deterministic trend:** The "fix" is simple — regress on time, work with residuals. Standard regression tools apply to the detrended series.
- **Stochastic trend:** Regressing a unit root process on time (or on another unit root process) can produce **spurious regression** — high $R^2$ and significant coefficients for variables that are actually unrelated. The fix is **differencing**, not detrending.

```{r unitroot-simulate}
set.seed(512)
n <- 200

# Stationary AR(1) for comparison
ar1_stat <- ts(arima.sim(model = list(ar = 0.8), n = n))

# Random walk (unit root)
rw <- ts(cumsum(rnorm(n)))
```

```{r unitroot-plot, fig.height=4.5}
p1 <- autoplot(ar1_stat) +
  labs(title = "Stationary AR(1), \u03c6 = 0.8", y = expression(y[t])) +
  theme_minimal(base_size = 12)

p2 <- autoplot(rw) +
  labs(title = "Random Walk (unit root)", y = expression(y[t])) +
  theme_minimal(base_size = 12)

p1 / p2

ggsave("output/stationary_vs_unitroot.pdf", width = 8, height = 5)
```

Both can look trending over short windows — this is the fundamental challenge. A stationary AR(1) with $\phi = 0.95$ and a random walk with $\phi = 1$ are visually very similar, yet they require completely different treatments. This is why visual inspection alone is insufficient and formal tests are essential.

### Augmented Dickey-Fuller (ADF) test

The ADF test is the most widely used unit root test.

- $H_0$: The series has a **unit root** (non-stationary)
- $H_1$: The series is **stationary**

The intuition: if a series is stationary, it exhibits **mean reversion** — a high value of $y_{t-1}$ should predict a *negative* change in the next period (the series corrects back toward its mean). The ADF test checks whether this mean-reversion pattern is significantly present. **Rejection** of $H_0$ (small $p$-value) is evidence for stationarity. Failure to reject means we cannot rule out a unit root — but this is not the same as confirming one, especially in small samples where ADF has low power. The regression behind the test and the Phillips-Perron alternative are detailed in Appendix B.

### KPSS test

The **KPSS test** reverses the burden of proof:

- $H_0$: The series is **stationary**
- $H_1$: The series has a **unit root**

This is the *opposite* of ADF. Rejection of $H_0$ (small $p$-value) is evidence for non-stationarity. Why do we need a test with the opposite null? Because **hypothesis tests can only provide evidence against the null, not for it**. If ADF fails to reject the unit root null, that could mean either (a) there really is a unit root, or (b) the test simply lacks power (common when $\phi$ is close to 1, or when the sample is small). KPSS addresses this by letting stationarity be the null — if KPSS also fails to reject, we have positive evidence *in favor of* stationarity, not just absence of evidence against it.

### Using ADF, PP, and KPSS together

**Best practice:** Use at least ADF and KPSS together. Their opposite null hypotheses make them complementary — convergent results strengthen your conclusion substantially.

| Scenario | ADF result | KPSS result | Interpretation |
|----------|-----------|-------------|----------------|
| Stationary | Reject $H_0$ | Fail to reject $H_0$ | Both agree: stationary |
| Unit root | Fail to reject $H_0$ | Reject $H_0$ | Both agree: non-stationary |
| Ambiguous | Fail to reject | Fail to reject | Insufficient evidence either way; consider larger sample or alternative tests |
| Contradictory | Reject | Reject | Rare; may indicate a structural break or fractional integration |

```{r stationarity-tests}
# Test the stationary AR(1)
cat("=== Stationary AR(1) ===\n")
adf.test(ar1_stat)
kpss.test(ar1_stat)

# Test the random walk
cat("\n=== Random Walk ===\n")
adf.test(rw)
kpss.test(rw)
```

For the AR(1): ADF rejects (evidence against unit root), KPSS fails to reject (no evidence against stationarity). Both tests agree — the series is stationary.

For the random walk: ADF fails to reject (no evidence against unit root), KPSS rejects (evidence against stationarity). Both tests agree — the series has a unit root.

A summary table of all hypothesis tests introduced in this lab (Ljung-Box, ADF, PP, KPSS) and their null hypotheses is collected in Appendix B for quick reference.

### Differencing

If a series has a unit root, **differencing** — taking $\Delta y_t = y_t - y_{t-1}$ — removes it:

```{r differencing}
rw_diff <- diff(rw)

cat("ADF test on differenced random walk:\n")
adf.test(rw_diff)
```

The differenced random walk is white noise (since $\Delta y_t = y_t - y_{t-1} = \varepsilon_t$). This is the key insight behind **ARIMA**($p, d, q$) models: the "I" stands for *integrated*, meaning the series needs to be differenced to become stationary. If a series needs to be differenced $d$ times to become stationary, then $d$ is the integration order. The ARIMA model applies ARMA($p, q$) to the differenced series. We will cover ARIMA estimation in depth in later labs, but the diagnostic logic is already here: if ADF/KPSS tests indicate a unit root, difference the series, then apply the identification table to the differenced data.

```{r cleanup-unitroot, include=FALSE}
rm(ar1_stat, rw, rw_diff, n, p1, p2)
```


## 2.3 Seasonality

Many economic and social time series exhibit **seasonal patterns** — regular fluctuations that repeat at fixed intervals (monthly, quarterly, etc.). Unemployment claims spike in winter, ice cream sales peak in summer, retail spending surges in December. Unlike trends, which are long-run movements, seasonality operates at a fixed periodicity: every 12 months for monthly data, every 4 quarters for quarterly data, etc.

Seasonality matters for diagnostics because, like trends, it inflates the ACF at specific lags (multiples of the seasonal period), masking the AR/MA structure we are trying to identify. The solution is the same principle as detrending: remove the seasonal component first, then examine the remainder for dynamic dependence. But *how* we remove seasonality matters — different approaches can distort the dynamics in the remainder. We will compare two approaches below.

### Simulating a seasonal series

Let's build a series with known components so we can check whether our deseasonalization procedure recovers them. The DGP has three parts: a gentle upward trend, a realistic (non-sinusoidal) monthly seasonal pattern, and an AR(1) noise component with $\phi = 0.4$.

```{r seasonal-simulate}
set.seed(512)
n_season <- 240  # 20 years of monthly data

# Trend component
trend_s <- 0.02 * (1:n_season)

# Seasonal component: irregular monthly constants (not a smooth sinusoid)
# These represent, e.g., monthly fluctuations in unemployment claims
monthly_pattern <- c(1.8, 1.3, 0.5, -0.4, -1.0, -1.5, -1.2, -0.5, 0.1, 0.6, 1.0, 1.5)
seasonal <- rep(monthly_pattern, length.out = n_season)

# AR(1) component
ar_base <- arima.sim(model = list(ar = 0.4), n = n_season)

# Combined series
y_seasonal <- ts(as.numeric(trend_s + seasonal + ar_base),
                 frequency = 12, start = c(2005, 1))
```

```{r seasonal-plot, fig.height=3.5}
autoplot(y_seasonal) +
  labs(title = "Monthly Series: Trend + Seasonality + AR(1)",
       y = expression(y[t]), x = "Year") +
  theme_minimal(base_size = 13)
```

The repeating annual pattern is clearly visible — the series oscillates up and down with a regular 12-month cycle, layered on top of a slight upward trend.

### How seasonality distorts the ACF and PACF

Before removing anything, let's see what the raw ACF and PACF look like:

```{r seasonal-acf-raw, fig.height=4}
p_acf_raw <- ggAcf(y_seasonal, lag.max = 48) +
  labs(title = "ACF — Raw Series (with seasonality)") +
  theme_minimal(base_size = 11)

p_pacf_raw <- ggPacf(y_seasonal, lag.max = 48) +
  labs(title = "PACF — Raw Series (with seasonality)") +
  theme_minimal(base_size = 11)

p_acf_raw + p_pacf_raw

ggsave("output/seasonal_acf_raw.pdf", width = 8, height = 3.5)
```

The ACF shows the characteristic seasonal signature: strong spikes at lags 12, 24, 36 (January is correlated with next January, and the January after that). The slow overall decay in the ACF reflects the trend. The PACF shows a sharp spike at lag 1 (the AR component) and another spike at lag 12 (the seasonal component) — notice that the PACF is somewhat cleaner than the ACF for identifying individual components, because it partials out intermediate lags.

In this state, the ACF is dominated by the seasonality and trend — we cannot read the AR(1) structure clearly. We need to remove these features first.

### Approach 1: STL decomposition (quick diagnostic)

The `stl()` function uses Loess (a nonparametric smoother) to separate trend, seasonal, and remainder components. This is useful for a quick visual diagnostic:

```{r stl-decompose, fig.height=6}
decomp_stl <- stl(y_seasonal, s.window = "periodic")
plot(decomp_stl, main = "STL Decomposition")
```

We can extract the remainder and check its ACF/PACF:

```{r stl-remainder-acf, fig.height=4}
remainder_stl <- ts(decomp_stl$time.series[, "remainder"], frequency = 12)

p1 <- ggAcf(remainder_stl) +
  labs(title = "ACF — STL Remainder") +
  theme_minimal(base_size = 11)

p2 <- ggPacf(remainder_stl) +
  labs(title = "PACF — STL Remainder") +
  theme_minimal(base_size = 11)

p1 + p2
```

The PACF cuts off after lag 1 — it recovers the AR(1) signature. However, there is a subtlety: **STL uses Loess to estimate the trend**, which is a nonparametric smoother. Because Loess can absorb some of the autocorrelation into its trend estimate, the remainder may have slightly attenuated dynamics compared to the true residual process. For quick visual diagnostics, STL is fine. But for the Box-Jenkins workflow — where we want to specify and estimate a parametric model — a more controlled approach is preferable.

### Approach 2: Linear detrending + seasonal means (recommended for Box-Jenkins)

This approach uses tools you already know (linear regression) and gives you more control over what gets removed and what stays in the residuals.

**Step 1: Detrend with `lm()`.** This removes a *linear* trend (a deterministic, parametric function of time) rather than a flexible Loess curve. The advantage is that `lm()` only removes a straight line, so any curvature or persistence in the data stays in the residuals — exactly where we want it for AR/MA diagnosis.

```{r manual-detrend}
# Step 1: Remove linear trend with lm()
trend_fit_s <- lm(y_seasonal ~ time(y_seasonal))
y_detrended_s <- ts(residuals(trend_fit_s), frequency = 12, start = c(2005, 1))
```

**Step 2: Estimate and remove the seasonal pattern.** We compute the average value for each month across all years — this gives us the typical monthly pattern. Then we subtract it. This is essentially what `decompose()` does under the hood.

```{r manual-deseasonalize}
# Step 2: Compute monthly means across years and subtract them
# Arrange the detrended series into a matrix: rows = years, columns = months
seasonal_means <- y_detrended_s |>
  matrix(ncol = 12, byrow = TRUE) |>
  colMeans()

# Repeat the monthly pattern for all years and subtract
seasonal_component <- ts(rep(seasonal_means, length.out = length(y_detrended_s)),
                         frequency = 12, start = c(2005, 1))
y_clean <- y_detrended_s - seasonal_component
```

The estimated seasonal means should be close to the true monthly constants we built in (not exact, because the AR noise and finite sample size introduce estimation error).

**Step 3: Examine the ACF and PACF of the cleaned series.**

```{r manual-clean-acf, fig.height=4}
p1 <- ggAcf(y_clean) +
  labs(title = "ACF — After lm() Detrend + Seasonal Removal") +
  theme_minimal(base_size = 11)

p2 <- ggPacf(y_clean) +
  labs(title = "PACF — After lm() Detrend + Seasonal Removal") +
  theme_minimal(base_size = 11)

p1 + p2

ggsave("output/seasonal_manual_clean_acf.pdf", width = 8, height = 3.5)
```

The PACF should cut off after lag 1 with a value close to 0.4, and the ACF should show the geometric decay characteristic of an AR(1). This matches the DGP we simulated.

### Comparing the two approaches

Both approaches recover the AR(1) structure, but they differ in how they treat the trend:

- **STL (Approach 1):** Quick and visual. However, because Loess is a nonparametric smoother, it can absorb some of the autocorrelation into the trend estimate, slightly attenuating the AR dynamics in the remainder.
- **lm() + seasonal means (Approach 2):** More controlled. The linear detrending removes only a straight line, so the full AR persistence is preserved in the residuals. This is the recommended approach for Box-Jenkins diagnostics when the trend is approximately linear.

Use STL for quick visual exploration or when the trend is clearly nonlinear. Use the `lm()` approach when your goal is to specify and estimate a parametric dynamic model.

```{r cleanup-seasonal, include=FALSE}
rm(y_seasonal, ar_base, seasonal, trend_s, monthly_pattern, n_season,
   decomp_stl, remainder_stl, trend_fit_s, y_detrended_s,
   seasonal_means, seasonal_component, y_clean,
   p_acf_raw, p_pacf_raw, p1, p2)
```


---

\newpage


# Part 3: Estimation and Model Selection

Parts 1 and 2 gave us the tools to *identify* candidate models — to look at a series and say "this looks like it might be AR(2)" or "both ACF and PACF tail off, so probably ARMA." Now we turn to *estimation*: fitting the proposed model to data, checking whether it adequately captures the temporal structure, and choosing among competing specifications when the diagnostics are ambiguous.


## 3.1 Estimating ARIMA models

Let's simulate an ARMA(1,1) and practice fitting models to it.

```{r estimation-data}
set.seed(2026)
arma_data <- arima.sim(model = list(ar = 0.7, ma = 0.4), n = 300)
```

R provides multiple interfaces for ARIMA estimation:

```{r estimation-compare}
# 1. Base R: arima()
fit_base <- arima(arma_data, order = c(1, 0, 1))
fit_base

# 2. forecast package: Arima() — same engine, more features
fit_forecast <- Arima(arma_data, order = c(1, 0, 1))
fit_forecast

# 3. Automatic selection: auto.arima() — searches over (p,d,q) and picks by AIC
fit_auto <- auto.arima(arma_data)
fit_auto
```

All three should recover coefficients close to the true values ($\phi = 0.7$, $\theta = 0.4$). The `auto.arima()` function searches over a grid of $(p, d, q)$ values and selects the model with the lowest AIC (or AICc by default).

**A word of caution about `auto.arima()`:** It is a useful tool for checking your manual identification, but it is not a substitute for understanding the Box-Jenkins procedure. Be aware of its limitations:

- It searches a **restricted grid** (by default, $p \leq 5$, $q \leq 5$, $d \leq 2$) using a stepwise algorithm that can get stuck in local optima. Setting `stepwise = FALSE` and `approximation = FALSE` gives a more thorough search but is slower.
- It may propose models that are **statistically adequate but not substantively interpretable**. For example, it might fit an ARIMA(0,1,1) when the data are actually stationary AR(1) with high persistence — both fit well, but they imply very different DGPs.
- When seasonal data are involved, `auto.arima()` can produce complex seasonal ARIMA specifications (e.g., ARIMA(1,0,1)(1,1,1)$_{12}$) that are difficult to interpret. In many applied settings, handling seasonality with dummy variables or STL decomposition and then fitting a simpler non-seasonal ARIMA is more transparent.

Use `auto.arima()` as a **sanity check** — if it picks a very different model from your manual identification, that is a signal to investigate further. But always verify its choice with residual diagnostics and substantive reasoning.

**Reading the output:** The key elements are:

- **Coefficients:** The estimated $\hat{\phi}$ and $\hat{\theta}$ (labeled `ar1` and `ma1`). The `intercept` is the estimated mean $\hat{\mu}$ of the process (not the same as a regression intercept — it is $\mu(1 - \phi)$ in the notation of some textbooks, but R reports $\hat{\mu}$ directly).
- **Standard errors:** In parentheses below the coefficients. A coefficient is "significant" if it is roughly twice its SE.
- **sigma^2:** The estimated innovation variance $\hat{\sigma}^2_\varepsilon$.
- **AIC / AICc / BIC:** Information criteria for model comparison. Lower is better.

```{r estimation-coefs}
# Extract coefficients and standard errors explicitly
coef(fit_forecast)
sqrt(diag(vcov(fit_forecast)))  # standard errors
```


## 3.2 Model diagnostics: residual checks

This is where the Box-Jenkins loop closes. Recall from Part 1 that the goal of model specification is to capture *all* the temporal structure in the data, leaving behind only white noise. A well-specified model should therefore produce residuals $\hat{e}_t$ that are indistinguishable from white noise: no autocorrelation, no patterns, no remaining structure. If the residuals still show temporal dependence, the model is incomplete — it has failed to account for some component of the DGP, and we need to go back and revise.

### Good fit: ARMA(1,1) fitted to ARMA(1,1) data

```{r diagnostics-good, fig.height=5}
# checkresiduals() does everything in one call:
# 1. Plots residuals over time
# 2. Shows ACF of residuals
# 3. Histogram of residuals
# 4. Runs Ljung-Box test
checkresiduals(fit_forecast)
```

The residuals look like white noise: no visible patterns in the time plot, ACF bars within the confidence bands, and the Ljung-Box test fails to reject the null of no autocorrelation (large $p$-value). This is a good fit — the model has captured the temporal structure, leaving only noise behind.

### Checking normality: Jarque-Bera test

The `checkresiduals()` output includes a histogram of the residuals, which gives a visual impression of whether they are approximately normal. For a formal test, we can use the **Jarque-Bera test**, which tests whether the residuals have the skewness and kurtosis consistent with a normal distribution:

- $H_0$: The residuals are normally distributed
- $H_1$: The residuals are not normally distributed (excess skewness or kurtosis)

```{r jarque-bera}
jarque.bera.test(residuals(fit_forecast))
```

If the $p$-value is large ($> 0.05$), the residuals are consistent with normality. Non-normality in residuals does not necessarily invalidate the model — ARIMA estimation by maximum likelihood is reasonably robust to moderate departures from normality, especially in larger samples. However, severe non-normality (heavy tails, strong skewness) can affect confidence intervals and hypothesis tests on the coefficients.

Together, the **Ljung-Box test** (no autocorrelation) and the **Jarque-Bera test** (normality) form the core residual diagnostic battery for Box-Jenkins modeling.

### Bad fit: AR(1) fitted to ARMA(1,1) data

What happens if we fit the wrong model? Let's fit a pure AR(1) to data that actually has an MA component:

```{r diagnostics-bad, fig.height=5}
fit_wrong <- Arima(arma_data, order = c(1, 0, 0))
checkresiduals(fit_wrong)
```

The Ljung-Box test should reject — the residuals from the AR(1) model still contain MA structure that the model failed to capture. When you see this, go back and revise your model specification.

You can also run the diagnostic steps individually:

```{r diagnostics-manual, fig.height=4}
# Residual ACF
resid_acf <- ggAcf(residuals(fit_wrong)) +
  labs(title = "ACF of Residuals — AR(1) fit to ARMA(1,1) data") +
  theme_minimal(base_size = 12)

resid_acf

# Ljung-Box test at specific lag
Box.test(residuals(fit_wrong), lag = 10, type = "Ljung-Box")
```

```{r cleanup-diagnostics, include=FALSE}
rm(fit_base, fit_wrong, resid_acf)
```


## 3.3 Model comparison with AIC

When the identification table is ambiguous — both ACF and PACF tail off, suggesting a mixed ARMA process — we cannot determine the exact orders $(p, q)$ from visual inspection alone. Instead, we fit several candidate models and compare their **AIC** (Akaike Information Criterion) values. AIC balances goodness of fit against model complexity: it rewards models that explain the data well but penalizes each additional parameter. Lower AIC = better model. BIC (Bayesian Information Criterion) works similarly but applies a heavier penalty for extra parameters, tending to favor simpler models.

```{r aic-comparison}
# Fit several candidate models
candidates <- list(
  "AR(1)"     = Arima(arma_data, order = c(1, 0, 0)),
  "AR(2)"     = Arima(arma_data, order = c(2, 0, 0)),
  "MA(1)"     = Arima(arma_data, order = c(0, 0, 1)),
  "MA(2)"     = Arima(arma_data, order = c(0, 0, 2)),
  "ARMA(1,1)" = Arima(arma_data, order = c(1, 0, 1)),
  "ARMA(2,1)" = Arima(arma_data, order = c(2, 0, 1)),
  "ARMA(1,2)" = Arima(arma_data, order = c(1, 0, 2))
)

# Extract AIC values
aic_table <- tibble(
  Model = names(candidates),
  AIC   = map_dbl(candidates, AIC),
  BIC   = map_dbl(candidates, BIC)
) |>
  arrange(AIC) |>
  mutate(delta_AIC = AIC - min(AIC))

aic_table
```

The ARMA(1,1) should have the lowest (or near-lowest) AIC, since it matches the true data-generating process. Models that are too simple (AR(1), MA(1)) will have higher AIC because they leave structure in the residuals. Models that are too complex (ARMA(2,1)) will have slightly higher AIC due to the complexity penalty.

**Over-fitting and under-fitting as a diagnostic strategy.** Beyond comparing AIC across a grid of candidates, a useful heuristic is to take your best model and deliberately **over-fit** it (add one more AR or MA term) and **under-fit** it (remove a term). If the extra parameter in the over-fit model is not significant and its AIC is higher, you have confirmation that your model is not too simple. If the under-fit model's residuals fail the Ljung-Box test, you have confirmation that the removed term was necessary. This "bracket" approach — checking one step above and one step below — is a practical safeguard against misspecification (Box-Steffensmeier et al., 2014, pp. 44–47).

**Rule of thumb:** Models within $\Delta$AIC $< 2$ of the best model are considered equally supported by the data. When two models are close in AIC, prefer the more **parsimonious** specification — the one with fewer parameters.

```{r aic-visual, fig.height=3.5}
# Visualize AIC comparison
aic_table |>
  mutate(Model = fct_reorder(Model, AIC)) |>
  ggplot(aes(x = Model, y = delta_AIC)) +
  geom_col(fill = "steelblue", alpha = 0.8) +
  geom_hline(yintercept = 2, linetype = "dashed", color = "firebrick") +
  annotate("text", x = 6, y = 3.5, label = "\u0394AIC = 2 threshold",
           color = "firebrick", size = 3.5) +
  labs(title = "Model Comparison: \u0394AIC from Best Model",
       y = "\u0394AIC", x = NULL) +
  coord_flip() +
  theme_minimal(base_size = 12)

ggsave("output/aic_comparison.pdf", width = 8, height = 3.5)
```

```{r cleanup-estimation, include=FALSE}
rm(arma_data, fit_forecast, fit_auto, candidates, aic_table)
```


---

\newpage


# Part 4: Practice — Diagnose Unknown Series

## Instructions

Now it is your turn. Load `data/practice_series.csv`. The file contains 6 time series (`series_A` through `series_F`), each with 120 observations. Some resemble real economic data; others are purely simulated. Your task: apply the full Box-Jenkins procedure to identify each one.

**For each series, follow the diagnostic workflow:**

1. **Plot the series** — look for trends, seasonality, level shifts, changing variance
2. **Check stationarity** — run ADF and KPSS tests; if a trend is present, detrend; if a unit root is present, difference
3. **Examine ACF and PACF** — match what you see to the identification table (Section 1.8)
4. **Propose an ARMA($p,q$) model** — based on whether ACF/PACF cut off or tail off
5. **Estimate the model** with `Arima()` or `auto.arima()`
6. **Check residuals** — do they pass the Ljung-Box test? Does the ACF of residuals show any remaining structure?
7. **If not, revise** — return to step 3, reconsider the orders, and re-estimate

```{r load-practice}
practice <- read_csv("data/practice_series.csv")
glimpse(practice)
```


## Exercise 1: Series A (guided example)

Series A has 120 observations. Work through the full diagnostic workflow step by step. This is a walk-through of the Box-Jenkins procedure — the remaining exercises follow the same logic with less scaffolding.

**Step 1: Plot the series.** Always start with a visual. Look for trends, level shifts, changing variance, and obvious patterns.

```{r ex1-plot}
# Your code here:
# Hint: convert to ts object first
#   y_a <- ts(practice$series_A)
#   autoplot(y_a) + labs(title = "Series A") + theme_minimal()
# Does it look stationary? Any visible trend?
# What is the approximate mean — is it centered at zero?
```

**Step 2: Check stationarity.** Run both ADF and KPSS. Remember their null hypotheses are opposite.

```{r ex1-stationarity}
# Your code here:
# Hint:
#   adf.test(y_a)    # H0: unit root. Want to REJECT for stationarity.
#   kpss.test(y_a)   # H0: stationary. Want to FAIL TO REJECT.
# If both agree the series is stationary, proceed to ACF/PACF.
# If there is a trend, detrend first:
#   trend_fit <- lm(y_a ~ seq_along(y_a))
#   y_a_detrended <- ts(residuals(trend_fit))
```

**Step 3: ACF and PACF.** This is the identification step. Match what you see to the table in Section 1.8.

```{r ex1-acf}
# Your code here:
# Hint:
#   p1 <- ggAcf(y_a) + labs(title = "ACF — Series A") + theme_minimal()
#   p2 <- ggPacf(y_a) + labs(title = "PACF — Series A") + theme_minimal()
#   p1 + p2
#
# Questions to answer:
#   - Does the ACF tail off (gradual decay) or cut off (drops to zero after lag q)?
#   - Does the PACF tail off or cut off (drops to zero after lag p)?
#   - What (p, q) does this suggest?
```

**Step 4: Propose and estimate a model.** Based on the ACF/PACF pattern, choose candidate values of $(p, q)$.

```{r ex1-estimate}
# Your code here:
# Based on what you saw in the ACF/PACF, choose (p, q):
#   fit_a <- Arima(y_a, order = c(?, 0, ?))
#   fit_a
# Also try auto.arima() and see if it agrees:
#   auto.arima(y_a)
# Compare the AIC of your model vs. auto.arima's choice.
```

**Step 5: Check residuals.** A good model leaves white noise residuals.

```{r ex1-residuals}
# Your code here:
# Hint: checkresiduals(fit_a)
# Three things to verify:
#   1. Residual time plot — no visible patterns
#   2. Residual ACF — all bars within the bands
#   3. Ljung-Box p-value > 0.05
# If residuals still show structure, go back to Step 4 and try a different (p,q).
```

**Step 6: Summarize your findings.**

Write a brief sentence: *"Series A appears to be a [process type] with estimated parameters [values]. The Ljung-Box test [passes/fails] with p = [value]."*


## Exercise 2: Series B

Series B has monthly structure. **This is important** — you must tell R about the frequency when creating the `ts` object, or the seasonal diagnostics will not work.

```{r ex2-plot}
# Your code here:
# IMPORTANT: set frequency = 12 for monthly data
#   y_b <- ts(practice$series_B, frequency = 12)
# Step 1: autoplot(y_b) — do you see a repeating annual pattern?
```

```{r ex2-acf}
# Your code here:
# Step 2: ACF with lag.max = 36 — look for spikes at lags 12, 24, 36
#   ggAcf(y_b, lag.max = 36) + theme_minimal()
# These seasonal spikes are the diagnostic signature.
```

```{r ex2-model}
# Your code here:
# Step 3: auto.arima() can handle seasonal ARIMA models
#   fit_b <- auto.arima(y_b)
#   fit_b
# Step 4: checkresiduals(fit_b)
```


## Exercise 3: Series C

```{r ex3-diagnose}
# Your code here:
# Hint: this one is a pure MA process
#   y_c <- ts(practice$series_C)
# Look at ACF — does it cut off? At which lag?
# Look at PACF — does it tail off?
```

```{r ex3-model}
# Your code here:
# Based on your ACF/PACF diagnosis:
#   fit_c <- Arima(y_c, order = c(0, 0, ?))
#   checkresiduals(fit_c)
```


## Exercise 4: Series D

```{r ex4-diagnose}
# Your code here:
# Hint: look carefully at the raw plot — is there a trend?
#   y_d <- ts(practice$series_D)
#   autoplot(y_d) + theme_minimal()
# If trend: detrend first with lm(), then examine ACF/PACF of residuals
#   trend_fit <- lm(y_d ~ seq_along(y_d))
#   y_d_detrended <- ts(residuals(trend_fit))
```

```{r ex4-model}
# Your code here:
# After detrending, both ACF and PACF may tail off — what does that suggest?
# (Hint: check the identification table for ARMA)
# Try a few candidate models and compare AIC:
#   fit_d1 <- Arima(y_d_detrended, order = c(1, 0, 0))
#   fit_d2 <- Arima(y_d_detrended, order = c(1, 0, 1))
#   AIC(fit_d1, fit_d2)
#   checkresiduals() on the best model
```


## Exercise 5: Series E

This series is tricky. Pay close attention to the stationarity tests.

```{r ex5-diagnose}
# Your code here:
# Hint: this series is very persistent
#   y_e <- ts(practice$series_E)
#   autoplot(y_e) + theme_minimal()
# Run ADF and KPSS tests — what do they say?
#   adf.test(y_e)
#   kpss.test(y_e)
# The ACF will decay very slowly. Is this a unit root or just high persistence?
```

```{r ex5-model}
# Your code here:
# Try differencing: diff(y_e). Re-examine ACF/PACF of the differenced series.
#   ggAcf(diff(y_e)) + theme_minimal()
#   ggPacf(diff(y_e)) + theme_minimal()
# Compare models:
#   fit_e1 <- Arima(y_e, order = c(1, 0, 0))  # stationary AR(1)
#   fit_e2 <- Arima(y_e, order = c(0, 1, 0))  # random walk (differenced)
#   AIC(fit_e1, fit_e2)
# Which specification do you prefer? Why?
```


## Exercise 6: Series F

```{r ex6}
# Your code here:
# Hint: what does "no structure" look like?
#   y_f <- ts(practice$series_F)
#   autoplot(y_f) + theme_minimal()
# Plot, ACF, PACF, Ljung-Box
# What is your diagnosis?
# If a series is white noise, any ARMA model is unnecessary — the best model
# is the mean: Arima(y_f, order = c(0, 0, 0))
```


**Summary checklist for the exercises:**

- [ ] Series A: identified process, estimated model, residuals pass Ljung-Box
- [ ] Series B: frequency set to 12, seasonal pattern detected, model fitted
- [ ] Series C: pure MA process identified via ACF cutoff
- [ ] Series D: trend detected and removed, ARMA structure identified
- [ ] Series E: near-unit-root behavior diagnosed, stationarity tests applied
- [ ] Series F: white noise correctly identified


---

\newpage


# Appendix A: The Box-Jenkins Flowchart (Self-Study)

The **Box-Jenkins method** is an iterative procedure for identifying, estimating, and diagnosing time series models. It was developed by George Box and Gwilym Jenkins in the 1970s and remains the standard framework for univariate time series modeling. The key insight of the method is that model specification is not a one-shot decision — it is a cycle of proposing, testing, and revising until the model adequately captures the temporal structure of the data.

Here is the full workflow:

**Step 1: Plot the series.**
Visual inspection comes first. Look for trends (upward/downward drift), seasonality (repeating cycles), level shifts (abrupt changes in the mean), outliers, and non-constant variance (heteroskedasticity). This step gives you a qualitative sense of what features the model will need to account for.

**Step 2: Check for stationarity.**
Use both ADF ($H_0$: unit root) and KPSS ($H_0$: stationary). Their opposite null hypotheses make them complementary — when both agree, you can be more confident. If the series is non-stationary, proceed to Step 3 before attempting to read the ACF/PACF.

**Step 3: Transform if needed.**

- If there is a **deterministic trend**: detrend by regressing on time and taking residuals. This is appropriate when the trend is a fixed, predictable function of $t$.
- If there is a **stochastic trend** (unit root): difference the series with `diff()`. Apply once for a single unit root, twice for a double unit root (rare in practice). The number of differences is $d$ in ARIMA($p, d, q$).
- If there is **seasonality**: consider seasonal differencing (`diff(y, lag = s)` where $s$ is the seasonal period) or STL decomposition. In some cases, both regular and seasonal differencing are needed.
- If there is **non-constant variance**: consider a log transformation or Box-Cox transformation to stabilize the variance before proceeding.

**Step 4: Examine ACF and PACF of the (transformed) series.**
Apply the identification table:

- ACF cuts off at lag $q$, PACF tails off $\Rightarrow$ MA($q$)
- ACF tails off, PACF cuts off at lag $p$ $\Rightarrow$ AR($p$)
- Both tail off $\Rightarrow$ ARMA($p,q$) — try small values of $p$ and $q$ (typically $p, q \leq 3$)
- Neither shows significant lags $\Rightarrow$ white noise — no model needed

**Step 5: Estimate candidate models.**
Fit the proposed model(s) using `arima()` or `Arima()`. Check that coefficients are statistically significant (roughly, coefficient / SE $> 2$). If a coefficient is not significant, consider dropping it — a more parsimonious model is preferred.

**Step 6: Check residuals.**

- Plot residuals over time — they should look like white noise (no visible patterns, no trends, no changing variance).
- ACF of residuals — no significant lags. All bars should fall within the confidence bands.
- Ljung-Box test — fail to reject $H_0$ of no autocorrelation ($p > 0.05$).
- `checkresiduals()` from the `forecast` package performs all three checks in a single call.

**Step 7: Iterate if needed.**
If residuals fail the diagnostic checks, return to Step 4 using the *residuals* as your guide. Look at the ACF/PACF of the residuals to identify what structure the model missed. Common scenarios: (a) the residual ACF shows a significant spike at lag 1, suggesting an MA(1) component is missing; (b) the residual PACF shows spikes at higher lags, suggesting higher-order AR terms are needed. Augment your model accordingly and re-estimate. Repeat until residuals pass.

**Step 8: Compare and select.**
If multiple candidate models pass the residual checks, compare them using AIC (or BIC). The model with the lowest information criterion is preferred. Remember the $\Delta$AIC $< 2$ rule: models within 2 AIC units of the best are considered roughly equally supported. When two models are close in AIC, prefer the more parsimonious one (fewer parameters).

This iterative loop — identify, estimate, diagnose, revise — is the essence of the Box-Jenkins approach. Mastering it requires practice: the more series you diagnose, the better you become at reading the signatures in the ACF and PACF.


---

\newpage


# Appendix B: Mathematical Details (Self-Study)

This appendix collects the formal definitions and derivations referenced in the main lab text. These are for your reference and deeper understanding — you do not need to memorize these formulas, but understanding them will help you interpret R output and troubleshoot when things go wrong.


## B.1 Autocorrelation function (ACF)

The **population autocorrelation** at lag $k$ is defined as:

$$
\rho_k = \frac{\text{Cov}(y_t, \, y_{t-k})}{\text{Var}(y_t)} = \frac{\gamma_k}{\gamma_0}
$$

where $\gamma_k = \text{Cov}(y_t, y_{t-k})$ is the **autocovariance** at lag $k$, and $\gamma_0 = \text{Var}(y_t)$ is the variance of the series (the autocovariance at lag 0). Key properties:

- $\rho_k$ is bounded between $-1$ and $1$, like any correlation.
- $\rho_0 = 1$ always (every observation is perfectly correlated with itself).
- $\rho_k = \rho_{-k}$ (the correlation is symmetric in leads and lags).

In practice, we estimate $\rho_k$ with the **sample autocorrelation**:

$$
r_k = \frac{\sum_{t=k+1}^{n}(y_t - \bar{y})(y_{t-k} - \bar{y})}{\sum_{t=1}^{n}(y_t - \bar{y})^2}
$$

This is the sample autocovariance at lag $k$ divided by the sample variance. As the sample size $n$ grows, $r_k$ converges to $\rho_k$. This is exactly what R's `acf()` function computes.


## B.2 Partial autocorrelation function (PACF)

The ACF at lag $k$ captures the *total* correlation between $y_t$ and $y_{t-k}$, including any indirect correlation transmitted through the intermediate lags $y_{t-1}, y_{t-2}, \dots, y_{t-k+1}$. For example, if $y_t$ depends on $y_{t-1}$ and $y_{t-1}$ depends on $y_{t-2}$, then $y_t$ will be correlated with $y_{t-2}$ even if there is no direct link between them.

The **partial autocorrelation function (PACF)** removes this indirect correlation. It measures the correlation between $y_t$ and $y_{t-k}$ *after controlling for* the linear effects of the intermediate lags. Formally, the PACF at lag $k$ is the coefficient $\alpha_k$ in the regression:

$$
y_t = \alpha_1 y_{t-1} + \alpha_2 y_{t-2} + \cdots + \alpha_k y_{t-k} + \varepsilon_t
$$

Equivalently, it is the correlation between the residuals of regressing $y_t$ on $(y_{t-1}, \dots, y_{t-k+1})$ and the residuals of regressing $y_{t-k}$ on the same set. This is the time series analog of partial regression coefficients in multiple regression.


## B.3 Ljung-Box test statistic

The **Ljung-Box** test statistic for testing whether the first $h$ autocorrelations are jointly zero is:

$$
Q(h) = n(n+2) \sum_{k=1}^{h} \frac{r_k^2}{n - k}
$$

Under $H_0$ (no autocorrelation), $Q(h) \sim \chi^2(h)$. The intuition: if the series is truly white noise, all the $r_k$ should be close to zero, so $Q$ will be small. If any autocorrelation is present, the squared terms inflate $Q$ and we reject.

The older **Box-Pierce** variant uses the simpler statistic $Q^* = n \sum r_k^2$, but the Ljung-Box adds the $(n+2)/(n-k)$ correction that improves the $\chi^2$ approximation in small samples.

**Choosing the number of lags $h$:** Too few lags and you may miss higher-order autocorrelation; too many and the test loses power. Common rules of thumb:

- For non-seasonal data: $h = \min(10, \, n/5)$
- For seasonal data with period $s$: $h = 2s$ (e.g., $h = 24$ for monthly data)
- When testing *residuals* from an ARIMA($p,d,q$) model, the degrees of freedom should be adjusted to $h - p - q$ (since $p + q$ autocorrelation parameters have already been estimated). The `checkresiduals()` function in `forecast` handles this automatically.


## B.4 Augmented Dickey-Fuller (ADF) test: the regression

The ADF test is based on a rearrangement of the AR(1) model:

$$
\begin{aligned}
y_t &= \rho \, y_{t-1} + \varepsilon_t \\
y_t - y_{t-1} &= (\rho - 1) \, y_{t-1} + \varepsilon_t \\
\Delta y_t &= \gamma \, y_{t-1} + \varepsilon_t, \quad \text{where } \gamma = \rho - 1
\end{aligned}
$$

Under $H_0$: $\gamma = 0$ (equivalently, $\rho = 1$, a unit root). Under $H_1$: $\gamma < 0$ (stationary). The "augmented" version includes lagged differences to account for higher-order serial correlation:

$$
\Delta y_t = \gamma \, y_{t-1} + \sum_{j=1}^{p} \psi_j \, \Delta y_{t-j} + \varepsilon_t
$$

The test statistic is the $t$-ratio on $\gamma$, but under the null of a unit root, this statistic does *not* follow a standard $t$-distribution — it follows the **Dickey-Fuller distribution**, which has more negative critical values than the usual $-1.96$. This is why a specialized test is needed rather than just running `lm()` and checking the $p$-value.


## B.5 Variance of an AR(1) process

The variance of a stationary AR(1) process can be derived from the moving average representation $y_t = \sum_{j=0}^{\infty} \phi^j \varepsilon_{t-j}$. Since the shocks $\varepsilon_{t-j}$ are independent with variance $\sigma^2_\varepsilon$:

$$
\text{Var}(y_t) = \sum_{j=0}^{\infty} \phi^{2j} \, \sigma^2_\varepsilon = \sigma^2_\varepsilon \sum_{j=0}^{\infty} (\phi^2)^j = \frac{\sigma^2_\varepsilon}{1 - \phi^2}
$$

The last step uses the geometric series formula $\sum_{j=0}^{\infty} r^j = 1/(1-r)$ for $|r| < 1$, which requires $|\phi| < 1$ (stationarity). When $\phi = 1$, the geometric series diverges and the variance is infinite — this is the mathematical reason why unit root processes are non-stationary.

The **half-life** of a shock in an AR(1) — the number of periods for the impact to halve — is $t_{1/2} = -\ln(2) / \ln(|\phi|)$.


## B.6 Phillips-Perron (PP) test

The **Phillips-Perron test** has the same null hypothesis as ADF (unit root) but handles serial correlation differently. Instead of adding lagged differences to the regression (the ADF approach), the PP test estimates the simpler regression $\Delta y_t = \gamma y_{t-1} + u_t$ and then applies a **nonparametric correction** to the test statistic to account for any serial correlation and heteroskedasticity in $u_t$.

In R: `PP.test(x)` from the `tseries` package.

In practice, ADF and PP usually agree. When they disagree, it may indicate the presence of a **structural break** — a large, one-time shock that creates the *appearance* of a unit root even if the underlying process is stationary around the break. ADF is particularly susceptible to this because structural breaks inflate the estimated persistence of the series.

```{r pp-test-appendix, eval=FALSE}
# Phillips-Perron test: same null as ADF, nonparametric correction
PP.test(y)    # H0: unit root. Small p-value → reject → stationary.
```


## B.7 Diagnostic toolkit: summary of hypothesis tests

| Test | R function | $H_0$ | $H_1$ | Use when... |
|------|-----------|--------|--------|-------------|
| **Ljung-Box** | `Box.test(x, type="Ljung-Box")` | No autocorrelation (white noise) | Serial correlation exists | Testing raw series or model residuals for remaining structure |
| **Jarque-Bera** | `jarque.bera.test(x)` | Normal distribution | Non-normal | Checking residual normality |
| **ADF** | `adf.test(x)` | Unit root (non-stationary) | Stationary | Checking whether a series needs differencing |
| **Phillips-Perron** | `PP.test(x)` | Unit root (non-stationary) | Stationary | Alternative to ADF; more robust to serial correlation |
| **KPSS** | `kpss.test(x)` | Stationary | Unit root | Complementing ADF — provides positive evidence *for* stationarity |


---

\newpage


# Appendix C: Further Reading

The following references provide deeper coverage of the topics introduced in this lab.

**Primary textbook for this course:**

- Box-Steffensmeier, J. M., Freeman, J. R., Hitt, M. P., & Pevehouse, J. C. W. (2014). *Time Series Analysis for the Social Sciences*. Cambridge University Press.
  - **Chapter 2** (pp. 22–67) covers the same ground as this lab: white noise, AR and MA processes, ACF/PACF derivations, the Box-Jenkins cycle, model estimation and residual diagnosis, AIC/BIC comparison, and seasonality. The applied examples (U.S. homicide rates, U.S.–Middle East relations) illustrate the full diagnostic workflow on real data.
  - **Section 2.1.1** (pp. 23–26): Why deterministic trend terms in OLS are problematic — sensitivity to outliers, first/last observation leverage, and the case for stochastic modeling.
  - **Section 2.2.3** (pp. 33–34): The mathematical equivalence (invertibility) of AR and MA processes — why identification is hard and mixed ARMA models are rare.
  - **Section 2.3** (pp. 34–43): Full derivations of the ACF for AR(1) ($\rho_k = \phi^k$) and MA(1), and the PACF via Yule-Walker equations.
  - **Section 2.4** (pp. 44–48): Residual analysis, the Q statistic (Ljung-Box), over-fitting/under-fitting, and metadiagnosis with AIC/BIC.

**Additional references:**

- Hyndman, R. J., & Athanasopoulos, G. (2021). *Forecasting: Principles and Practice* (3rd ed.). OTexts. Available free at [otexts.com/fpp3](https://otexts.com/fpp3). Chapters 8–9 cover ARIMA modeling with R code examples using the `fable` package.
- Shumway, R. H., & Stoffer, D. S. (2017). *Time Series Analysis and Its Applications* (4th ed.). Springer. Chapter 3 covers ARIMA models in greater mathematical depth.
