Lecture notes. Written from the material of Berkeley Stat 153, licensed CC BY 4.0. These are notes, not a transcript: the material has been reorganised and rewritten. This adaptation carries the same licence, and the original is linked above.

69. Autoregression and Yule's AR(2)¶

What this covers¶

Every regression model built so far in the course — trend regression, broken-stick regression, harmonic (sinusoidal) regression — predicts $y_t$ from functions of the time index $t$: the design matrix $X$ is made of columns like $t$, $(t-c)_+$, $\cos 2\pi ft$. This chapter introduces autoregression (AR), where $X$ is instead made of the series' own past values: $y_{t-1}$, $y_{t-2}$, and so on. It defines $AR(1)$, $AR(2)$ and the general $AR(p)$ model, shows that fitting one is still ordinary least squares once the design matrix is relabelled, and works out how to forecast several steps ahead by feeding earlier forecasts back in as if they were data. It then asks where this idea came from: Udny Yule's 1927 attempt to model the 11-year sunspot cycle, which is worked through in enough detail to see why a recursive model succeeds where a fixed sinusoid does not. Assumes: linear regression by least squares and by maximum likelihood, and harmonic (sinusoidal) regression $y_t = \beta_0+\beta_1\cos2\pi ft+\beta_2\sin2\pi ft+\varepsilon_t$.

From regression on time to regression on the past¶

The family $ARIMA$ — AutoRegressive Integrated Moving Average — names the three ingredients this part of the course assembles: $AR$, $MA$ (moving average) and $I$ (integration, i.e. differencing). This chapter is about the first one.

Every model considered up to now has the generic linear-regression shape $$y = X\beta + \varepsilon,$$ for example a broken-stick trend $y_t = \beta_0+\beta_1 t+\beta_2(t-c)_++\varepsilon_t$, with $$y = \begin{pmatrix} y_1\\ \vdots\\ y_n\end{pmatrix}, \qquad X = \begin{bmatrix} 1 & t_1 & (t_1-c)_+ \\ \vdots & \vdots & \vdots \\ 1 & t_n & (t_n-c)_+\end{bmatrix}.$$ In every one of these, $X$ is built out of the time index $t$ (or functions of it, such as a frequency). Autoregression keeps exactly the same equation $y=X\beta+\varepsilon$ and only changes what goes into $X$: its columns are the series' own past values, $y_{t-1}, y_{t-2},\dots$, rather than functions of $t$. That single change is the whole idea, and it is also why fitting an AR model needs no new machinery: it is ordinary least squares on a differently-built $X$.

The $AR(p)$ model¶

$AR(1)$. The simplest case regresses $y_t$ on the single preceding value: $$y_t = \phi_0 + \phi_1 y_{t-1} + \varepsilon_t, \qquad \varepsilon_t \overset{\text{iid}}{\sim} N(0,\sigma^2), \qquad t=2,3,\dots,n.$$ The "$1$" is the order: the number of lags used. Writing it as $y=X\beta+\varepsilon$, $$y = \begin{pmatrix} y_2\\ \vdots\\ y_n\end{pmatrix}, \qquad X = \begin{bmatrix} 1 & y_1 \\ \vdots & \vdots \\ 1 & y_{n-1}\end{bmatrix}, \qquad \beta = \begin{pmatrix} \phi_0\\ \phi_1\end{pmatrix}.$$ Notice the fitting sample starts at $t=2$: there is no $y_0$ to serve as a regressor for $y_1$, so the first observation is used only as a predictor, never as a response.

$AR(2)$ adds a second lag, $$y_t = \phi_0+\phi_1y_{t-1}+\phi_2y_{t-2}+\varepsilon_t, \qquad t=3,4,\dots,n,$$ $$y=\begin{pmatrix} y_3\\ \vdots\\ y_n\end{pmatrix}, \qquad X=\begin{bmatrix} 1 & y_2 & y_1\\ \vdots & \vdots & \vdots\\ 1 & y_{n-1} & y_{n-2}\end{bmatrix},$$ and now the sample loses its first two observations to the role of predictors.

General $AR(p)$: $$y_t = \phi_0+\phi_1y_{t-1}+\dots+\phi_py_{t-p}+\varepsilon_t, \qquad t=p+1,\dots,n,$$ $$y=\begin{pmatrix} y_{p+1}\\ \vdots\\ y_n\end{pmatrix}, \qquad X=\begin{bmatrix} 1 & y_p & y_{p-1} & \dots & y_1\\ \vdots & \vdots & \vdots & & \vdots\\ 1 & y_{n-1} & y_{n-2} & \dots & y_{n-p}\end{bmatrix}.$$ Each row of $X$ holds the $p$ values immediately preceding the response in that row, oldest last; the price of using $p$ lags is always the first $p$ observations, which appear only inside $X$ and never as a $y$.

Fitting: still least squares¶

The unknown parameters are $\phi_0,\phi_1,\dots,\phi_p$ and $\sigma$. Because $y=X\beta+\varepsilon$ is again an ordinary linear model with iid Gaussian errors, the same equivalence used for trend and harmonic regression applies unchanged: maximizing the Gaussian likelihood over $\phi_0,\dots,\phi_p$ is the same problem as minimizing the residual sum of squares, $$\min_{\phi_0,\dots,\phi_p} \sum_{t=p+1}^n \big(y_t-\phi_0-\phi_1y_{t-1}-\dots-\phi_py_{t-p}\big)^2,$$ so sm.OLS(y, X).fit() on the $y,X$ built above gives $\hat\phi_0,\dots,\hat\phi_p$ directly — no new estimation method is needed, only a new design matrix. The MLE for the noise scale is the usual plug-in formula, $$\hat\sigma_{\text{MLE}} = \sqrt{\frac{1}{n-p}\sum_{t=p+1}^n\big(y_t-\hat\phi_0-\hat\phi_1y_{t-1}-\dots-\hat\phi_py_{t-p}\big)^2}.$$

Forecasting: feed the predictions back in¶

Forecasting $y_{n+1},\dots,y_{n+k}$ from a fitted $AR(p)$ runs into a problem one step ahead: the first forecast, $$\hat y_{n+1} = \hat\phi_0+\hat\phi_1y_n+\hat\phi_2y_{n-1}+\dots+\hat\phi_py_{n+1-p},$$ only needs observed data, since every lag it calls on is at or before time $n$. The second, $$\hat y_{n+2} = \hat\phi_0+\hat\phi_1\hat y_{n+1}+\hat\phi_2y_n+\dots+\hat\phi_py_{n+2-p},$$ already needs $\hat y_{n+1}$, which does not exist yet as data — so the model's own first forecast is substituted in its place. Continuing this way, forecasting $i$ steps ahead uses $$\boxed{\hat y_{n+i} = \hat\phi_0+\hat\phi_1\hat y_{n+i-1}+\dots+\hat\phi_p\hat y_{n+i-p}}, \qquad \boxed{\hat y_j = y_j \text{ if } j\le n}, \qquad i=1,2,\dots$$ i.e. wherever a lag falls on or before time $n$, use the real observation; wherever it falls after $n$, use the forecast already computed for it. This "plug the forecast back in as if it were data" recursion is the only way to reach beyond one step ahead, since an $AR(p)$ model only ever knows how to predict the very next value.

Where $AR(2)$ came from: Yule and the sunspot cycle (1927)¶

The autoregressive idea was not introduced as an abstract generalization of regression — it was invented to solve a specific problem. Udny Yule was trying to model the number of sunspots $y_t$ observed per year, which rises and falls in a cycle roughly $11$ years long but not on a fixed schedule: the peaks are not evenly spaced and not of even height.

The first attempt was the harmonic regression already familiar from earlier in the course: $$\text{Model ①:} \qquad y_t = \underbrace{\beta_0+\beta_1\cos2\pi ft+\beta_2\sin2\pi ft}_{s(t)} + \varepsilon_t, \qquad \varepsilon_t\overset{\text{iid}}{\sim}N(0,\sigma^2),$$ with unknowns $f,\beta_0,\beta_1,\beta_2,\sigma$ — the frequency $f$ would come out near $1/11$ per year. This treats the sunspot count as a perfectly regular signal $s(t)$ plus independent noise on top of it.

Turning the signal into a recursion. The signal $s(t)$ satisfies a second-order differential equation: differentiating twice, $$s''(t) = -(2\pi f)^2\big[\beta_1\cos2\pi ft+\beta_2\sin2\pi ft\big] = -(2\pi f)^2\big(s(t)-\beta_0\big),$$ using $s(t)-\beta_0=\beta_1\cos2\pi ft+\beta_2\sin2\pi ft$. This is the equation of a simple harmonic oscillator: the "acceleration" of $s$ is always pulling it back toward $\beta_0$, proportional to its current displacement.

The discrete analogue of this — replacing the second derivative by a second difference — turns out to hold exactly, not just approximately, for any sinusoid at frequency $f$: $$\big(s_t-s_{t-1}\big)-\big(s_{t-1}-s_{t-2}\big) = 2\big[\cos2\pi f-1\big]\big(s_{t-1}-\beta_0\big) \iff s_t = \beta_0+\beta_1\cos2\pi ft+\beta_2\sin2\pi ft.$$ Adding the observation noise back in (a sinusoid observed with error induces noise $\eta_t$ in this recursion too) gives $$\big(y_t-y_{t-1}\big)-\big(y_{t-1}-y_{t-2}\big) = 2\big(\cos2\pi f-1\big)\big(y_{t-1}-\beta_0\big)+\eta_t,$$ which rearranges into $$y_t = \phi_0+\phi_1y_{t-1}-y_{t-2}+\eta_t.$$ This is a special case of $AR(2)$ — special because the coefficient on $y_{t-2}$ is forced to be exactly $-1$ by the sinusoid identity it came from. Letting that coefficient be a free parameter $\phi_2$ instead of fixing it at $-1$ gives the general model $$\text{Model ②:} \qquad y_t = \phi_0+\phi_1y_{t-1}+\phi_2y_{t-2}+\varepsilon_t,$$ which is the $AR(2)$ model proper.

Why the recursion, and not just the sinusoid. Models ① and ② look like alternative routes to the same curve, but they describe very different kinds of irregularity, and a physical analogy makes the difference concrete. Picture a mass on a spring, obeying Hooke's law $$s''(t) = -\tfrac{k}{m}\big(s(t)-\beta_0\big),$$ a clean oscillation around the rest position $\beta_0$ with a fixed period set by the spring constant $k$ and the mass $m$. Now introduce randomness by throwing stones at it, in one of two ways:

Scenario 1 (matching Model ①): the spring swings exactly as the equation says, and the randomness is measurement noise added afterward, $y_t = s(t)+\varepsilon_t$. Removing the noise would reveal a perfectly regular, unchanging oscillation.
Scenario 2 (matching Model ②, in continuous time): the stones hit the object itself, so the noise enters the equation of motion, $y''(t) = -k\big(y(t)-\beta_0\big)+\eta_t$. Now $y(t)$ is smooth — it still solves a differential equation, so it cannot jump — but every impact nudges its future trajectory, so its amplitude and effective period keep drifting rather than repeating exactly.

Sunspot activity looks like Scenario 2, not Scenario 1: the 11-year cycle wanders in height and timing rather than repeating a fixed waveform. That is the reason Yule reached for the recursive model rather than stopping at the harmonic regression he started with.

Left to right in each panel is time. Model ① measures a perfectly regular sinusoid with independent noise on each observation, so the underlying cycle never changes. Model ② is an AR(2) recursion in which fresh noise enters the recursion itself every step, so it stays smooth but its cycle length and amplitude drift — the behaviour Yule needed to match the sunspot record.

The Yule model, as originally written, keeps the link to the frequency $f$ that motivated it: $$y_t = \phi_0+\phi_1y_{t-1}-y_{t-2}+\varepsilon_t, \qquad \phi_1 = 2\cos2\pi f,$$ or, moving the lag-2 term to the left, $$y_t+y_{t-2} = \phi_0+\phi_1y_{t-1}+\varepsilon_t.$$ The fitted coefficient $\phi_1$ can be read back as an implied cycle frequency $f=\tfrac{1}{2\pi}\arccos(\phi_1/2)$ — the trace of the sinusoid the recursion grew out of, even though the general $AR(2)$ no longer requires the coefficient on $y_{t-2}$ to be exactly $-1$.

Sources¶

Handwritten lecture notes (fall 2025), reconstructed by a model from a PDF with no text layer — every equation there is flagged unverified in the source and has been cross-checked here where possible (the sinusoid difference-equation identity in the history section was independently verified against the trigonometric sum formula): 01-lagged-regression.md (the "From regression on time to regression on the past" and "$AR(1)$/$AR(2)$" sections), 02-ar.md (general $AR(p)$, estimation, and the forecasting recursion), and 03-history-how-ar-models-were-invented.md (Yule 1927 and the spring analogy) — all from berkeley-stat153/fall-2025, HandwrittenNotesLectureSixteen153248Fall2025.pdf, CC BY 4.0.
No slide deck, transcript, or problem set was supplied for this lecture, so none is drawn on here.
The notes' opening agenda lists "Neural Networks (RNN, LSTM)" as the lecture's second topic; the supplied material contains no content on it, so it is not covered in this chapter.
Not contained in the supplied material: the original Yule (1927) paper itself, and any numerical fit of the sunspot series (no data or fitted coefficients appear in the notes).

← 68. Regression with an Unknown Changepoint · Contents · 70. Bayesian Shrinkage and Variance Models →