Econ 6376 — Applied Time Series Econometrics
This chapter has one main reading path and three optional layers.
If you missed lecture, read the main text, run the Core code, and complete the Core Practice problems. Then return to the optional sections that address your questions or interests.
By the end of this module, you should be able to:
arima() is usually numerical.stats::arima() and forecast::Arima(), compare them on AIC/BIC/AICc, and contrast manual selection with forecast::auto.arima().Module 3 ended with a powerful tool: the ACF/PACF identification table.
| Process | ACF | PACF |
|---|---|---|
| AR(p) | Tails off | Cuts off at lag \(p\) |
| MA(q) | Cuts off at lag \(q\) | Tails off |
| ARMA(p,q) | Tails off | Tails off |
For pure AR and pure MA processes, the fingerprint is decisive. A sharp PACF cutoff at lag 2 says AR(2). A sharp ACF cutoff at lag 1 says MA(1). These rules are clean, reliable, and testable in the simulation gallery. But look at the third row. For ARMA(p,q), both correlograms tail off. The fingerprint tells you the class — “this is probably an ARMA” — but it cannot tell you the orders. Six students can look at the same ARMA(1,1) correlogram and propose six different models, and none of them would be obviously wrong from the correlogram alone.
This is the correlogram trap. It is not a failure of the identification rules; it is the boundary of what those rules were designed to do. And it is the problem that motivates everything in this module.
One thing has not changed since last module: the mixing board. Every non-seasonal slider on the master equation was already turned on by the end of Module 3. Today is about how we choose the settings, not about which sliders exist — we restate the board in §4.2.3.
Recall the ARMA(1,1) from the Module 3 simulation gallery (\(\phi = 0.6, \theta = 0.4\)). We can regenerate it using the arma_simulator() function built in Module 3 — per the course convention, from Module 4 onward we source() it from its canonical home in helpers/simulators.R instead of re-defining it:
library(forecast); library(ggplot2); library(patchwork)
source("../helpers/simulators.R") # arma_simulator(), built live in Module 3
set.seed(1985)
y_arma <- arma_simulator(n = 500, phi = 0.6, theta = 0.4)
autoplot(y_arma) + ggtitle("ARMA(1,1): phi = 0.6, theta = 0.4") + theme_bw()
(ggAcf(y_arma, lag.max = 20) | ggPacf(y_arma, lag.max = 20))
Look at the two correlograms. The ACF decays — it tails off, with no sharp cutoff. The PACF also tails off, with no clean cutoff at any particular lag. If you were forced to read these as a fingerprint, you might call it AR(1), or MA(1), or ARMA(1,1), or maybe AR(2). All four guesses are defensible from the picture alone, and the data came from a DGP we wrote ourselves. On real data, where the DGP is unknown and the sample correlograms are noisier, the ambiguity is worse.
The rest of the module fills the gap with three tools, introduced in order:
We will use (1) to generate candidates, (2) to score them, and (3) next module to confirm the winner. That three-step loop is the core of the Box-Jenkins methodology, which we meet in Section 4.2.
A natural question from anyone trained in cross-sectional regression: why not pick the model with the highest \(R^2\)? Because adding parameters always weakly increases \(R^2\), and for a flexible family like ARMA you can drive \(R^2\) arbitrarily high by piling on lags. \(R^2\) has no complexity penalty and is not a model-selection tool — it is a descriptive statistic.
What about adjusted \(R^2\)? Its penalty is ad hoc and weak — it was designed as a quick degrees-of-freedom correction for linear regression, not as a principled model-selection criterion for time series. In simulation studies comparing IC performance across candidate sets, adjusted \(R^2\) is consistently among the worst-performing criteria. It underpunishes complexity relative to BIC and lacks the theoretical grounding (KL divergence for AIC, Bayesian posterior approximation for BIC) that gives the proper information criteria their asymptotic properties. Friends do not let friends use adjusted \(R^2\) for model selection. If you encounter it in an ARMA context, replace it with AIC, AICc, or BIC — any of them. This is not a close call.
We need a scoring rule that builds a principled fit-versus-complexity trade-off into its structure. That is what AIC and BIC do.
The most influential framework for time series modeling is the Box-Jenkins methodology, introduced in Box and Jenkins (1970) and still the organizing template for applied work more than fifty years later. It has four steps:
Today’s module covers Step 2 in full, adds information criteria as the intermediate scoring step between Steps 2 and 3, and previews Step 3 just enough to set up Module 5.
This is the organizing analogy for the module. A detective does not solve a case by pointing at one suspect. They:
No step skipped. A detective who goes straight from “I have fingerprints” to “I have a verdict” is the detective who convicts the wrong person. A time series modeler who goes straight from an ACF plot to auto.arima() without looking at the diagnostics is the same detective.
It is worth seeing explicitly where each lecture sits in the workflow:
| Box-Jenkins Step | Course Location |
|---|---|
| 1. Identification | Module 2 (\(d\)), Module 3 (\(p, q\) from ACF/PACF) |
| 2. Estimation | Module 4 (MLE, this module) |
| IC scoring (between 2 and 3) | Module 4 (AIC, BIC, AICc) |
| 3. Diagnostic checking | Module 5 (residual diagnostics, Ljung-Box, seasonality) |
| 4. Forecasting | Block 4 |
Reproduce the master equation from Modules 1 and 3 one more time:
\[y_t = \alpha + \delta t + \sum_{j=1}^{p} \phi_j y_{t-j} + \sum_{l=1}^{q} \theta_l \epsilon_{t-l} + \epsilon_t\]
All of the non-seasonal terms are now on and fittable. The seasonal block is still off — that is Module 5’s job. What changes in this module is not which knobs are available on the mixing board, but how we pick the settings. Module 3 wrote down the DGP. Module 4 fits it.
Plain-language definition. The likelihood is the probability of seeing the data you saw, treated as a function of the parameters. It asks: “Given these particular parameter values, how surprising is it that I got this exact data set?” Bigger likelihood means less surprising means the parameters fit the data well.
The hindsight framing. Likelihood is the frequentist’s version of hindsight. You already saw the data. You freeze it in place and then ask, “which parameter values would have made this the most likely thing to happen?” Those are the maximum likelihood estimates.
The eyewitness lineup analogy. You are the witness to a crime. You already saw the event happen — the data is fixed. The lineup is the set of possible parameter values. MLE is the answer to “which suspect in the lineup would have made what I saw most likely?” You are not changing the crime; you are walking down the lineup and asking, for each suspect, “how plausible is it that this person produced what I witnessed?” The suspect who scores highest is the MLE.
For an ARMA(p,q) with Gaussian innovations \(\epsilon_t \sim N(0, \sigma^2)\), the (exact) likelihood of observing \(y_1, y_2, \ldots, y_T\) as a function of the parameter vector \(\boldsymbol{\vartheta} = (\alpha, \phi_1, \ldots, \phi_p, \theta_1, \ldots, \theta_q, \sigma^2)\) is:
\[L(\boldsymbol{\vartheta} \mid y_1, \ldots, y_T) = f(y_1, \ldots, y_T \mid \boldsymbol{\vartheta})\]
where \(f\) is the joint density of the observations under the model. We take logs for two reasons: products become sums, and the numerical range becomes manageable. The log-likelihood is:
\[\ell(\boldsymbol{\vartheta}) = \log L(\boldsymbol{\vartheta} \mid y_1, \ldots, y_T)\]
The maximum likelihood estimate is:
\[\hat{\boldsymbol{\vartheta}}_{\text{MLE}} = \arg\max_{\boldsymbol{\vartheta}} \; \ell(\boldsymbol{\vartheta})\]
That is the entire statement. We do not derive the closed form — it does not have one for general ARMA, and that is why we hand the problem to a numerical optimizer.
Two flavors, both names only:
arima() calls method = "ML".method = "CSS-ML" in arima() initializes with CSS and polishes with ML.When you see method = "CSS-ML" in an arima() call, that is what it means. For the homework and for this course, always use the default unless you have a reason not to.
This is the bridge between the regression world you already know and the likelihood world we are entering.
Claim. For a pure AR(p) with Gaussian innovations, the conditional MLE / CSS fit (conditioning on \(y_1, \ldots, y_p\)) is OLS-like: it gives the same coefficient estimates as OLS under the usual setup, so the AR(p) fit can be read like a regression.
\[y_t = \alpha + \phi_1 y_{t-1} + \ldots + \phi_p y_{t-p} + \epsilon_t, \qquad t = p+1, \ldots, T.\]
Intuition, not proof. When the innovations are Gaussian, maximizing the likelihood is equivalent to minimizing the sum of squared errors of the one-step-ahead predictions. For a pure AR, those one-step-ahead predictions are linear in past \(y\)’s with no unobserved quantities, so the conditional likelihood and CSS objective collapse to the same least-squares problem you solve with OLS. For MA or ARMA models, the one-step-ahead prediction depends on past innovations, which are unobserved and must be recursively reconstructed — that is what breaks the OLS equivalence and forces the numerical optimizer.
Live verification. Simulate an AR(1) and fit it both ways:
set.seed(1985)
y <- arma_simulator(n = 500, phi = 0.6)
fit_css <- arima(y, order = c(1, 0, 0), include.mean = FALSE, method = "CSS")
fit_ols <- lm(y[-1] ~ y[-length(y)] - 1)
coef(fit_css)["ar1"] # approx 0.6 ar1
0.5573313 coef(fit_ols) # approx 0.6 — same numbery[-length(y)]
0.5573313 The coefficients agree to every printed digit. That is the conditional AR fit behaving like OLS. The default arima() fit adds the exact-ML polish, which is why the right classroom phrase here is “OLS-like” rather than “literally OLS.” For MA or ARMA, MLE is genuinely new — because the past errors \(\epsilon_{t-1}, \epsilon_{t-2}, \ldots\) are not observable and the optimizer has to reconstruct them as it goes.
The briefest possible numerical moment, to demystify what log likelihood = -234.1 means in an arima() printout.
Take an AR(1) with \(\phi = 0.6\), \(\sigma^2 = 1\), and two observations \(y_1 = 1.0\), \(y_2 = 0.8\). Under the model, \(y_2 \mid y_1 \sim N(\phi y_1, \sigma^2) = N(0.6, 1)\), so:
\[\log f(y_2 \mid y_1; \phi = 0.6) = -\tfrac{1}{2}\log(2\pi) - \tfrac{1}{2}(0.8 - 0.6)^2 = -0.9189 - 0.02 = -0.939.\]
Now try \(\phi = 0.9\): the conditional mean becomes \(0.9\), and the log density becomes \(-0.9189 - \tfrac{1}{2}(0.8 - 0.9)^2 = -0.924\). Slightly higher. That is MLE on two points — walk the parameter, watch the log density move, pick the top. Scale that up to 500 points and a multi-dimensional \(\boldsymbol{\vartheta}\) and you have the optimization problem arima() is solving.
The log likelihood in an arima() printout is just a sum of \(T\) log densities like the ones above. It is not mysterious; it is arithmetic.
arima() PrintoutFit a model on the simulated ARMA(1,1) from Section 4.1 and walk through every field in the output:
fit <- arima(y_arma, order = c(1, 0, 1))
fit
Call:
arima(x = y_arma, order = c(1, 0, 1))
Coefficients:
ar1 ma1 intercept
0.5349 0.4287 -0.1538
s.e. 0.0503 0.0573 0.1356
sigma^2 estimated as 0.9791: log likelihood = -704.67, aic = 1417.33The printout shows:
sigma^2 — the estimated innovation variance \(\hat{\sigma}^2\).log likelihood — the value \(\ell(\hat{\boldsymbol{\vartheta}})\) from Section 4.3.2. This is the sum of all the log densities evaluated at the MLE.aic — the Akaike Information Criterion, which is the next topic.Remember our notation convention: \(\epsilon_t\) is the true DGP innovation; \(e_t = y_t - \hat{y}_t\) is the residual from an estimated model. The sigma^2 in the printout estimates \(\text{Var}(\epsilon_t)\), but the residuals you extract with residuals(fit) are \(e_t\), not \(\epsilon_t\).
The generic structure of every information criterion we will use in this course is:
\[\text{IC} = -2 \log L(\hat{\boldsymbol{\vartheta}}) + \text{penalty}(k, T)\]
where \(k\) is the number of estimated parameters (including \(\sigma^2\) and, if present, the intercept) and \(T\) is the sample size.
Read it piece by piece:
The Olympic all-around analogy. Think of IC scoring like gymnastics. You get points for your routine (log-likelihood, the fit component) and then a difficulty tariff is applied (penalty, the complexity component). A simple routine done perfectly can beat a difficult routine done sloppily. The whole scoring system exists because otherwise gymnasts would all attempt triple twisting double backflips and fall on their faces — and every ARMA modeler would fit ARMA(10,10) on 200 data points.
Three penalties, three stories.
AIC (Akaike, 1974):
\[\text{AIC} = -2 \log L + 2k\]
The penalty is \(2\) per parameter, regardless of the sample size. Akaike motivated AIC (intuitively) by minimizing Kullback-Leibler divergence between the fitted model and the (unknown) true DGP — “how wrong, on average, is the model as a description of nature?” AIC does not ask whether the true model is in the candidate set; it asks which candidate is the least wrong approximation. This is the philosophical underpinning for why AIC tends to prefer slightly larger models: a slightly-too-large model is a better approximation than a slightly-too-small one, even if the extra parameters are not “truly” present.
BIC (or SBC) (Schwarz, 1978):
\[\text{BIC} = -2 \log L + k \log T\]
The penalty is \(\log T\) per parameter. Schwarz motivated BIC (intuitively) by a Bayesian approximation to the marginal likelihood — “how well does this model explain the data under a prior over parameters?” The penalty is heavier than AIC for any \(T \geq 8\) (since \(\log 8 \approx 2.08 > 2\)), and it gets heavier as \(T\) grows. This makes BIC more aggressive at eliminating extra parameters.
AICc (Hurvich & Tsai, 1989):
\[\text{AICc} = -2 \log L + 2k + \frac{2k(k+1)}{T - k - 1}\]
The extra term \(2k(k+1)/(T - k - 1)\) vanishes as \(T \to \infty\), so AICc and AIC agree in large samples. In small samples, AICc penalizes extra parameters more heavily than AIC, which helps prevent AIC’s well-known overfitting tendency in short series. Rule of thumb: use AICc whenever \(T / k < 40\).
A numerical moment. Suppose we have \(T = 200\) and we are comparing a model with \(k = 3\) parameters to a model with \(k = 5\). The BIC penalty increases by \(2 \log 200 \approx 10.6\) when we add two parameters. The AIC penalty increases by \(2 \times 2 = 4\). So AIC is happy to add parameters that improve \(-2 \log L\) by more than 4; BIC needs the improvement to exceed 10.6. BIC is the tougher grader.
This is the question students ask about most, and it deserves an honest answer rather than a slogan.
Two different questions, two different answers:
| Property | AIC | BIC |
|---|---|---|
| Motivating criterion | KL divergence from truth | Bayesian posterior / marginal likelihood |
| Penalty weight | \(2\) per parameter | \(\log T\) per parameter |
| Asymptotic property | Efficient for forecasting | Consistent for selection |
| Finite-sample tendency | Overfits — picks too large a model | Underfits — picks too small a model |
| When to prefer | Forecasting, prediction, dense DGPs | Parsimony, sparse DGPs, identifying the “true” model |
What “efficient” and “consistent” mean here:
This is the key practical guidance. You should choose your primary IC before fitting, based on your modeling goal — not after, based on which one gives you the answer you prefer.
If your goal is to explain — you believe the true DGP is in (or near) your candidate set and you want the most parsimonious description of the data-generating mechanism — use BIC. Its consistency property means it will find the true model in the limit, and its heavier penalty keeps you from overparameterizing.
If your goal is to predict — you want the model that will forecast best out of sample, and you do not particularly care whether it is “true” — use AIC (or AICc). Its efficiency property means it targets minimum forecast error, and its connection to cross-validation gives it a direct link to the out-of-sample loss you actually care about.
Or, as the in-class mnemonic goes: AIC likes big models and it cannot lie. BIC thinks flat models are the thing.
The choice is not about which criterion is “better” — it is about what question you are asking the data. Explanation and prediction are different tasks with different optimal trade-offs between fit and parsimony. Picking one a priori based on your task is the disciplined move; reporting both lets you see how sensitive the choice is.
The philosopher and the accountant. AIC is the philosophy professor — “all models are wrong, but some are useful; pick the one that will be least wrong on average when forecasting.” BIC is the accountant — “tell me exactly which model generated this data, assuming you really are in my set.” Both grade the exam; they grade on different scales because they are asking different questions.
The honest message. In practice, reporting both is the professional default. When they agree, the choice does not matter for this dataset — you have more confidence. When they disagree, AIC picked a bigger model whose extra parameters help prediction a little but are not needed for a parsimonious explanation. That disagreement tells you you are in the region where the explain/predict trade-off is live, which is useful information, not a failure.
The textbook and software conventions vary. Enders reports AIC and SBC (= BIC). stats::arima() reports AIC by default. forecast::auto.arima() uses AICc by default. forecast::Arima() reports all three. Students should not panic if different software shows different fields — know what each field is and read it. (The Core field guide in §4.5.1 repeats the parts of this you need for the problem sets.)
AIC, BIC, and AICc are the three you will use in this course — but they are far from the only letters in the soup. Frame this the way you would the government’s three-letter agencies: there are more of them than you think, and each has its own jurisdiction.
| Abbreviation | Full Name | Penalty | One-Liner |
|---|---|---|---|
| AIC | Akaike Information Criterion | \(2k\) | KL divergence, the baseline |
| BIC / SBC | Bayesian / Schwarz Bayesian IC | \(k \log T\) | Bayesian motivation, harsher penalty |
| AICc | Corrected AIC | \(2k + 2k(k+1)/(T-k-1)\) | Small-sample fix for AIC |
| HQIC | Hannan-Quinn IC | \(2k \log \log T\) | A hybrid between AIC and BIC |
| DIC | Deviance IC | — | Bayesian, for hierarchical models |
| WAIC | Watanabe-Akaike IC | — | Fully Bayesian, \(\approx\) leave-one-out CV |
| CIC | Conditional IC | — | Conditional on random effects |
The HQIC form explicitly (Hannan & Quinn, 1979):
\[\text{HQIC} = -2 \log L + 2k \log \log T\]
The penalty \(2 \log \log T\) sits between AIC’s flat \(2\) and BIC’s \(\log T\) for moderate sample sizes. For \(T = 200\), AIC penalizes \(2\) per parameter, HQIC penalizes \(2 \log \log 200 = 2 \times 1.667 \approx 3.33\), and BIC penalizes \(\log 200 \approx 5.30\). HQIC is a compromise grader — not as generous as AIC, not as strict as BIC.
You will encounter HQIC in Block 7 when vars::VARselect() reports it alongside AIC, BIC, and a fourth criterion (FPE, the final prediction error) for VAR lag selection. When you see it in that printout, you will know what it is: the same \(-2 \log L + \text{penalty}\) structure with yet another penalty.
The message. We will use AIC, AICc, and BIC for the rest of this course. When you encounter the others in the wild — in a paper, in a software printout, in a job interview — now you know they exist and that they all share the same \(-2 \log L + \text{penalty}\) structure with different penalties. The differences are in what question the penalty is answering.
R gives us three functions for fitting ARMA/ARIMA models, and students often confuse them. Here is the field guide.
stats::arima() — base R. Fits ARIMA models by MLE (default method = "CSS-ML"). Returns parameter estimates, standard errors, log-likelihood, and AIC. Does not return AICc or BIC by default. Quiet, no package loads, always available. Use when you want a clean minimal fit.
forecast::Arima() — Rob Hyndman’s wrapper around stats::arima(). Same fitting engine, but (a) returns AIC, AICc, and BIC in one object; (b) handles drift terms more gracefully; (c) plays nicely with the rest of the forecast package (e.g., forecast(), checkresiduals()). This is the course default for the rest of Block 3 and Block 4.
forecast::auto.arima() — the Hyndman-Khandakar algorithm (2008). Searches over \((p, d, q)\) (and seasonal orders if seasonal data is detected) and returns the model that minimizes an IC — AICc by default. Options to force BIC (ic = "bic") or AIC (ic = "aic"). Also performs unit-root tests (KPSS by default) to pick \(d\) automatically. It is fast because it uses a stepwise search rather than a full grid, and it is mostly reliable but not perfect.
Return to the simulated ARMA(1,1) from Section 4.1 and fit it three ways:
# Manual with stats::arima — minimal output
fit_m <- arima(y_arma, order = c(1, 0, 1))
# Manual with forecast::Arima — IC fields populated
fit_f <- Arima(y_arma, order = c(1, 0, 1))
# Automatic — let the algorithm choose
fit_auto <- auto.arima(y_arma, ic = "aicc", stepwise = TRUE,
approximation = FALSE)
summary(fit_f)Series: y_arma
ARIMA(1,0,1) with non-zero mean
Coefficients:
ar1 ma1 mean
0.5349 0.4287 -0.1538
s.e. 0.0503 0.0573 0.1356
sigma^2 = 0.985: log likelihood = -704.67
AIC=1417.33 AICc=1417.41 BIC=1434.19
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -8.318111e-05 0.9894949 0.7833288 70.86365 197.2721 0.9004849
ACF1
Training set -0.01206216summary(fit_auto)Series: y_arma
ARIMA(2,0,2) with zero mean
Coefficients:
ar1 ar2 ma1 ma2
-0.2208 0.4601 1.1830 0.2643
s.e. 0.1194 0.0571 0.1271 0.0931
sigma^2 = 0.9797: log likelihood = -702.86
AIC=1415.72 AICc=1415.84 BIC=1436.79
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.04827924 0.9858417 0.7787737 77.70084 199.1255 0.8952485
ACF1
Training set -0.002863087Does auto.arima() land on (1,1)? On this particular seed, no — it reports an ARMA(2,2) with zero mean, whose AICc comes out about a point and a half better than the true specification’s. Rerun with a different seed and the answer changes; sometimes it finds (1,1), sometimes it does not. Automation is a labor-saver, not a conscience. It runs the same logic we will run by hand in the next section, with one specific IC and a specific search path. If you disagree with its answer, you are allowed to override it — but you should be able to say why.
Now do it manually. Fit all nine candidates in the \(\{0, 1, 2\} \times \{0, 1, 2\}\) grid and compile the IC table:
set.seed(1985)
y_arma <- arma_simulator(n = 500, phi = 0.6, theta = 0.4)
candidates <- expand.grid(p = 0:2, q = 0:2)
ic_table <- do.call(rbind, lapply(seq_len(nrow(candidates)), function(i) {
p <- candidates$p[i]; q <- candidates$q[i]
fit <- tryCatch(Arima(y_arma, order = c(p, 0, q)),
error = function(e) NULL)
if (is.null(fit)) return(NULL)
data.frame(
p = p, q = q,
k = length(fit$coef) + 1, # +1 for sigma^2
loglik = fit$loglik,
AIC = fit$aic,
AICc = fit$aicc,
BIC = fit$bic
)
}))
ic_table[order(ic_table$AIC), ] p q k loglik AIC AICc BIC
9 2 2 6 -702.2585 1416.517 1416.687 1441.805
5 1 1 4 -704.6655 1417.331 1417.412 1434.189
6 2 1 5 -703.6746 1417.349 1417.471 1438.422
8 1 2 5 -703.9274 1417.855 1417.976 1438.928
3 2 0 4 -709.9999 1428.000 1428.081 1444.858
7 0 2 4 -718.2895 1444.579 1444.660 1461.437
2 1 0 3 -725.7688 1457.538 1457.586 1470.181
4 0 1 3 -745.3764 1496.753 1496.801 1509.397
1 0 0 2 -912.2906 1828.581 1828.605 1837.010Look at the table this single 500-observation realization produced. AIC (and AICc) put ARMA(2,2) first, with the true ARMA(1,1) less than a point behind; BIC puts ARMA(1,1) first by a comfortable margin. One realization, one disagreement — AIC’s lighter penalty lets the two extra parameters pay for themselves, while BIC’s \(\log T\) penalty cuts them. Is that a fluke of this seed or a systematic pattern? You cannot tell from one draw. That is the hook for the Monte Carlo in Section 4.6.
auto.arima() and When to OverrideTrust it when:
Override it when:
auto.arima() does have root safeguards, but they can still slip through.)The honest framing: auto.arima() is a very good intern. You still read its work.
This is the workflow students will use in every problem set for the rest of the course:
Arima(). Record AIC, AICc, BIC.We will do Steps 1–6 in this module on real data in Section 4.7. Steps 7 and 8 wait for Modules 5 and 7. This checklist is your compass for the rest of the course.
The clean story from Section 4.4 is that IC rewards the true model: fit a grid, sort by AIC or BIC, and the true ARMA(1,1) rises to the top. Section 4.5.3 already dented that story — on our single realization, AIC crowned ARMA(2,2) and only BIC found the truth.
The honest story is that IC is a scoring rule with measurable error rates, and both AIC and BIC get it wrong a nontrivial fraction of the time — even when the true model is in the candidate set, even with \(T = 500\), even with a perfectly behaved ARMA(1,1) DGP. This is the assumption-break moment of the module: the place in the four-step ladder (concept \(\to\) math \(\to\) simulate \(\to\) real data) where we stress the machinery until we can see where it bends.
The design: 500 Monte Carlo replications of an ARMA(1,1) with \(\phi = 0.6\), \(\theta = 0.4\), \(\sigma = 1\), and \(T = 500\). For each replication, fit ARMA(p,q) for \((p, q) \in \{0, 1, 2\} \times \{0, 1, 2\}\) — nine candidates — and record which model AIC selects and which model BIC selects. Tabulate the selection frequencies.
The full code, with line-by-line comments so that students can replicate the experiment (fair warning: 500 replications \(\times\) 9 fits is 4,500 maximum-likelihood estimations, so this chunk takes a couple of minutes to run):
# --- Monte Carlo: IC selection frequencies for ARMA(1,1) ---
# True DGP: y_t = 0.6 * y_{t-1} + 0.4 * eps_{t-1} + eps_t
# We fit all 9 models in {0,1,2} x {0,1,2} on each replication
# and record which (p,q) is chosen by AIC and by BIC.
set.seed(1999) # reproducible
n_reps <- 500 # number of Monte Carlo replications
grid <- expand.grid(p = 0:2, q = 0:2) # 9 candidate (p,q) pairs
# Pre-allocate character vectors to store the winner for each rep
pick_aic <- character(n_reps)
pick_bic <- character(n_reps)
for (r in seq_len(n_reps)) {
# 1. Simulate one ARMA(1,1) realization, T = 500
y <- arma_simulator(n = 500, phi = 0.6, theta = 0.4)
# 2. Fit all 9 candidates; tryCatch guards against occasional
# convergence failures on extreme draws
fits <- lapply(seq_len(nrow(grid)), function(i) {
tryCatch(
Arima(y, order = c(grid$p[i], 0, grid$q[i])),
error = function(e) NULL
)
})
# 3. Drop any candidates that failed to converge
ok <- !sapply(fits, is.null)
aics <- sapply(fits[ok], function(f) f$aic)
bics <- sapply(fits[ok], function(f) f$bic)
# 4. Record the (p,q) pair with the smallest AIC and BIC
pick_aic[r] <- paste0("(", grid$p[ok][which.min(aics)], ",",
grid$q[ok][which.min(aics)], ")")
pick_bic[r] <- paste0("(", grid$p[ok][which.min(bics)], ",",
grid$q[ok][which.min(bics)], ")")
}
# 5. Tabulate: how often did each criterion pick each model?
table(pick_aic)pick_aic
(1,1) (1,2) (2,0) (2,1) (2,2)
363 36 29 28 44 table(pick_bic)pick_bic
(1,1) (1,2) (2,0) (2,1) (2,2)
448 4 44 3 1 # 6. Express as percentages
round(100 * prop.table(table(pick_aic)), 1)pick_aic
(1,1) (1,2) (2,0) (2,1) (2,2)
72.6 7.2 5.8 5.6 8.8 round(100 * prop.table(table(pick_bic)), 1)pick_bic
(1,1) (1,2) (2,0) (2,1) (2,2)
89.6 0.8 8.8 0.6 0.2 The exact numbers depend on the seed — the table below collects the percentages from the run above (seed 1999); change the seed and they move a little, but the pattern is stable across seeds, and the pattern is the whole point of the exercise:
| Model | AIC selects (%) | BIC selects (%) |
|---|---|---|
| (1,1) — the truth | 72.6 | 89.6 |
| (2,2) | 8.8 | 0.2 |
| (1,2) | 7.2 | 0.8 |
| (2,0) | 5.8 | 8.8 |
| (2,1) | 5.6 | 0.6 |
1. Neither scoring rule is a sure thing, even under laboratory conditions. AIC misses the true model about 27% of the time and BIC about 10% — and the true model is literally in the candidate set, the sample is a comfortable \(T = 500\), and the DGP is a perfectly behaved Gaussian ARMA(1,1). These error rates grow quickly as the sample shrinks (Practice Problems 3 and 4 let you watch that happen). On real data, where no model in the set is the true DGP (because there is no true ARMA DGP for macroeconomic data), IC is even less reliable as a “truth finder.”
2. AIC systematically leans larger. Look at where the “mistakes” go. When AIC misses, its picks go overwhelmingly to (2,2), (1,2), and (2,1) — models with more parameters than the truth — about 22% of all replications. BIC almost never makes those picks (under 2%). This is the finite-sample echo of the asymptotic theorems from Section 4.4.3: you can see the “AIC overfits” half of the slogan emerge from the numbers in front of you. One honest wrinkle: at this sample size, BIC’s own misses do not go to smaller models — they go to (2,0), an AR(2) that carries the same number of parameters as the truth and approximates its dynamics. With \(T = 500\), dropping a genuinely needed term costs so much likelihood that neither criterion underfits; the classic “BIC underfits” tendency shows up in shorter samples, which is exactly what Practice Problem 4 asks you to demonstrate at \(T = 100\).
3. The disagreement between AIC and BIC is informative. In roughly a fifth of the replications, AIC and BIC pick different models. If both rules pick the same model, you have some confidence. If they disagree, you have a live decision. That decision should be informed by what you plan to do with the model — forecasting (lean AIC) vs. description and parsimony (lean BIC) — and confirmed by the Module 5 residual diagnostics.
The takeaway for real data. If AIC picks the “wrong” model a quarter of the time on data it should nail, how should that change how much you trust it on the differenced UNRATE series in ten minutes? Good answer: IC is one line of evidence, not a verdict; pair it with residual diagnostics (Module 5) and with judgment about what the model is for.
We return to the FRED US civilian unemployment rate series that has been our running example since Module 2. In Module 2 we established that the level series is \(I(1)\) and should be differenced. In Module 3 we looked at the ACF and PACF of the differenced series and guessed “something low-order.” Now we have the tools to pin it down.
A note on the sample window. We estimate on January 1960 through December 2019 — deliberately stopping before the pandemic. In April 2020 the unemployment rate jumped from 4.4 to 14.8 in a single month, a one-observation shock dozens of times the size of a typical monthly change. An outlier that extreme dominates every sample autocovariance: include it and the ACF and PACF of the differenced series flatten toward zero, and every fitted candidate collapses toward a near-cancelling ARMA. Handling outliers and structural breaks properly is beyond this module’s toolkit, so the honest move is to model the pre-COVID sample and say explicitly that we did so. (Instructive exercise: rerun this section on the full cached sample and watch the correlogram flatten.)
To pull the series live, store your FRED key in the environment (never hardcode it in course files) and run:
library(fredr)
fredr_set_key(Sys.getenv("FRED_KEY"))
unrate_raw <- fredr(series_id = "UNRATE",
observation_start = as.Date("1960-01-01"),
observation_end = as.Date("2019-12-31"))
unrate <- ts(unrate_raw$value, start = c(1960, 1), frequency = 12)
d_unrate <- diff(unrate)So that these notes render identically for everyone — with or without an API key or an internet connection — the code below loads the same series from the course’s cached copy in data/UNRATE.csv (the same cache Modules 2 and 3 used) and applies the same sample window. The workflow from here on is identical either way.
unrate_csv <- read.csv("../data/UNRATE.csv") # cached FRED pull; see data/README.md
# Same estimation window as the live pull above: 1960-01 through 2019-12
# (see the sample-window note in the text).
unrate_csv <- subset(unrate_csv,
observation_date >= "1960-01-01" &
observation_date <= "2019-12-01")
unrate <- ts(unrate_csv$UNRATE, start = c(1960, 1), frequency = 12)
d_unrate <- diff(unrate)autoplot(d_unrate) +
ggtitle("First-differenced UNRATE, 1960-2019") + theme_bw()
(ggAcf(d_unrate, lag.max = 36) | ggPacf(d_unrate, lag.max = 36))
Before we proceed, a caveat that is both honest and pedagogically deliberate. The unemployment rate is a monthly series, and the correlograms show bumps at lags 12 and 24 — negative, modest, but visible — which is a clue worth checking for annual seasonality. A seasonal ARIMA is a candidate to compare against, but the ACF by itself is not proof. We are going to proceed with a non-seasonal ARMA today knowingly, as a pedagogical exercise, and treat the seasonal ACF bumps as evidence that our model may be incomplete. This is a deliberate teaching choice: we want to run the Box-Jenkins workflow on a real series inside the Module 4 toolkit, and then discover that the toolkit is insufficient, which motivates Module 5.
Reading the ACF/PACF of the differenced series: both tail off. There is positive autocorrelation spread over roughly the first six lags — the largest spikes sit at lags 2 through 4 rather than in a single clean spike at lag 1 — and the PACF shows the same early structure. Neither correlogram has a clean cutoff anywhere. This is the ARMA correlogram trap from Section 4.1, live on real data: the class is “low-order ARMA,” and the orders are anyone’s guess. A defensible candidate list:
candidates <- list(
"(1,0,0)" = c(1, 0, 0),
"(0,0,1)" = c(0, 0, 1),
"(1,0,1)" = c(1, 0, 1),
"(2,0,1)" = c(2, 0, 1),
"(1,0,2)" = c(1, 0, 2),
"(2,0,2)" = c(2, 0, 2)
)
fits <- lapply(candidates, function(ord) Arima(d_unrate, order = ord))
ic_table <- data.frame(
model = names(fits),
AIC = sapply(fits, AIC),
AICc = sapply(fits, function(f) f$aicc),
BIC = sapply(fits, BIC)
)
ic_table[order(ic_table$AIC), ] model AIC AICc BIC
(1,0,2) (1,0,2) -547.3492 -547.2651 -524.4599
(2,0,2) (2,0,2) -546.8717 -546.7538 -519.4046
(2,0,1) (2,0,1) -543.2168 -543.1327 -520.3275
(1,0,1) (1,0,1) -515.8414 -515.7854 -497.5299
(1,0,0) (1,0,0) -451.1455 -451.1120 -437.4120
(0,0,1) (0,0,1) -448.7133 -448.6798 -434.9798On the 1960–2019 window, the three criteria agree: ARMA(1,2) ranks first on AIC, AICc, and BIC. When the criteria agree, the choice is easy and you have more confidence in it — this is the “agreement” branch of the Section 4.4.4 guidance. But look one row down: ARMA(2,2) sits within about half a point of the winner on AIC, essentially a tie, while BIC separates the two decisively. That widening gap as you read across the columns is the \(k \log T\) penalty doing exactly the work Section 4.4.2 said it would: the extra AR lag buys almost no likelihood, so the tougher grader punishes it. Had we chosen AIC alone, the (1,2)-vs-(2,2) call would have been a coin flip; parsimony breaks the tie.
Compare with auto.arima():
# Non-seasonal (our toolkit today)
auto.arima(d_unrate, seasonal = FALSE)Series: d_unrate
ARIMA(5,0,0) with zero mean
Coefficients:
ar1 ar2 ar3 ar4 ar5
-0.0216 0.1939 0.1669 0.1704 0.0728
s.e. 0.0373 0.0370 0.0373 0.0371 0.0377
sigma^2 = 0.0271: log likelihood = 279.29
AIC=-546.59 AICc=-546.47 BIC=-519.12# Seasonal (what a full toolkit would do — preview of Module 5)
auto.arima(unrate)Series: unrate
ARIMA(2,1,2)(1,0,2)[12] with drift
Coefficients:
ar1 ar2 ma1 ma2 sar1 sma1 sma2 drift
1.1717 -0.2581 -1.1927 0.4346 0.3809 -0.6243 -0.0919 -0.0014
s.e. 0.1758 0.1674 0.1643 0.1346 0.1322 0.1316 0.0613 0.0077
sigma^2 = 0.02477: log likelihood = 311.61
AIC=-605.23 AICc=-604.97 BIC=-564.03The non-seasonal auto.arima() lands on an AR(5) with zero mean — a model that was not even on our candidate list, and (check against the table above) one whose AIC is slightly worse than our ARMA(1,2). That is the stepwise search path showing through: auto.arima() is fast precisely because it does not exhaust the grid, and here the shortcut walked past the model our manual six-candidate search found. A very good intern; read its work. The seasonal auto.arima() applied to the undifferenced levels moves to a full SARIMA — seasonal AR and MA terms at the 12-month frequency — confirming that once seasonal structure is allowed, the algorithm wants it. That is still a model-comparison result, not proof. We are using it as a reminder to compare seasonal candidates in Module 5, then verify the residuals.
Choose — publicly, with reasoning — a working candidate. A defensible choice is whichever model wins BIC (parsimony wins ties) and does not have weird parameter estimates. Here that is ARMA(1,2), which happens to win all three columns. Fit it, save it, and hold it up: “This is our best guess today. It is not final. Next week we will check whether the residuals look like white noise, and if they do not, seasonality is one of the first things we will inspect.”
Show the residual plot and the residual ACF — but do not diagnose them formally. That is Module 5.
e_hat <- residuals(fits[["(1,0,2)"]])
(autoplot(e_hat) +
ggtitle("Residuals from ARMA(1,2) on diff(UNRATE)") + theme_bw()) /
(ggAcf(e_hat, lag.max = 36) + ggtitle("Residual ACF") + theme_bw())
Look at the residual ACF. The low lags are quiet — the ARMA(1,2) has absorbed the short-run dynamics. But see that negative spike at lag 12, with an echo at lag 24? That is a seasonal clue, not a verdict. Our non-seasonal ARMA may not explain it. That is the cliffhanger for Module 5: the model we just picked is probably incomplete, and the evidence is staring at us from the residual ACF. Next module, we will learn how to read that evidence (Ljung-Box test, residual ACF/PACF) and how to fix it (seasonal differencing and the full SARIMA specification).
Three important topics remain open after this module. All of them will be addressed before the midpoint of the course.
Residual diagnostics. We previewed the residual ACF in Section 4.7.6 but did not interpret it formally. The tools for doing so — the Ljung-Box \(Q\) statistic, residual ACF/PACF reading, and the logic of “residuals should be white noise if the model is correct” — are the subject of Module 5.
Seasonality. The lag-12 spike in the UNRATE residual ACF is a clue that seasonal structure may still be present in our non-seasonal ARMA fit. Module 5 introduces the SARIMA\((p,d,q)(P,D,Q)_s\) specification, seasonal differencing, and the “airline model” as the canonical seasonal example.
Forecasting. The Box-Jenkins workflow ends with “forecast using the fitted model.” We have not done that yet. Forecasting methods, out-of-sample evaluation, and forecast comparison are Block 4.
After Module 5, students will be fitting the full (S)ARIMA\((p,d,q)(P,D,Q)_s\) specification end-to-end, with diagnostics to confirm the fit. After Block 4, they will be forecasting with it. The workflow we built in this module — identify, estimate, score, check — is the template for all of that. The only thing that changes is the range of models on the suspect list.
“The model with the lowest AIC is the right model.” No. IC is one line of evidence among several. Residual diagnostics (Module 5) can overturn an IC winner. On real data, “the right model” usually does not exist.
“AIC and BIC always agree.” They do not. When they disagree, AIC usually picks the larger model and BIC usually picks the smaller one. The disagreement is informative, not a problem.
“BIC is better because its penalty is stronger.” No. BIC is better for selection consistency under a “true model exists” assumption. AIC is better for forecasting when the true model is not in the set. Different goals, different tools.
“MLE is something different from OLS.” For pure AR(p), the conditional fit is OLS-like under Gaussian innovations. The exact/default arima() fit is still numerical, and MLE becomes genuinely new only once the model has moving-average or ARMA structure.
“auto.arima() picks the correct model.” It picks the model that minimizes its chosen IC over its chosen search path. That is often, but not always, the model a thoughtful analyst would pick. Read its work.
“If I add more parameters, my log-likelihood always goes up, so I should add parameters.” True that the log-likelihood weakly increases with parameters. That is exactly why the penalty exists. The log-likelihood is the race car; the penalty is the speed limit.
“Log-likelihood is a probability.” No. Likelihood is a density evaluated at the data, treated as a function of the parameters. It can be greater than 1. What matters is its comparison across parameter values, not its absolute magnitude.
“Conditional and exact MLE give different answers in any meaningful way.” For \(T \geq 100\), they usually agree to the third or fourth decimal. The default method = "CSS-ML" in arima() initializes with conditional sum of squares and polishes with exact ML — that is the right default for this course and for most practical work.
A convention warning when you cross-reference. Enders reports the BIC under the name SBC (Schwarz Bayesian Criterion) — same formula, different label — and, as in Module 3, writes AR coefficients as \(a_i\) where we write \(\phi_j\). Hamilton Chapter 5 is the rigorous treatment of exact and conditional maximum likelihood for ARMA models, if you want the derivations we deliberately skipped. Hyndman & Athanasopoulos cover the same selection workflow from the forecasting angle (their Chapter 9), including why auto.arima() defaults to AICc and how its stepwise search works.
Conditional AR fit \(\approx\) OLS verification. Fit an AR(1) with arima(y, order = c(1,0,0), method = "CSS") and with lm() (after constructing the lagged variable by hand), and verify that the coefficients agree to several decimal places. Write one paragraph explaining why the conditional fit is OLS-like for AR models but not for MA models. Then compare that result with the default arima() fit and note how close it is.
Single-realization IC table. Simulate a single ARMA(1,1) with \(\phi = 0.6\), \(\theta = 0.4\), and \(T = 500\). Fit all nine candidates in \(\{0, 1, 2\} \times \{0, 1, 2\}\). Report the full IC table (AIC, AICc, BIC). Which model is ranked first by each criterion? Do they agree?
Scaled-down Monte Carlo. Repeat Problem 2 as a Monte Carlo with 200 replications (a scaled-down version of the in-class assumption-break from Section 4.6). Report AIC and BIC selection frequencies. Write one paragraph on the overfit-vs-underfit pattern: when AIC misses, does it tend to pick larger or smaller models? Same question for BIC.
AICc vs. AIC on a short series. Take an AR(2) DGP with \(\phi_1 = 0.7\), \(\phi_2 = 0.2\), and \(T = 100\) (a small sample). Run a Monte Carlo with 200 replications, fitting a \(\{0, 1, 2, 3\} \times \{0, 1, 2, 3\}\) grid. Compare AIC and AICc selection frequencies. Does AICc do better than AIC at avoiding overfit? How does BIC compare to both?
UNRATE IC table. For the differenced FRED UNRATE series (1960–2019, as in Section 4.7), fit the six candidates used in class (Section 4.7.4) and reproduce the IC table. Which candidate would you pick, and why? (Two to three sentences, including whether AIC and BIC agree.)
Build your own grid search. Without using auto.arima(), write a function that takes a time series and a maximum \(p\) and \(q\), fits the full \((p, q)\) grid with Arima(), and returns a ranked IC table with AIC, AICc, and BIC columns. Test it on the simulated ARMA(1,1) from Problem 2 and verify that your function’s top-ranked model matches what auto.arima() picks — and if the two disagree, explain which one you believe and why. (This is essentially auto.arima() without the stepwise shortcut — a full grid search.)
The correlogram trap. For ARMA(p,q), both the ACF and the PACF tail off, so the identification rules from Module 3 cannot determine the orders \((p, q)\) by themselves. Information criteria fill the gap.
Box-Jenkins is the workflow. Identification \(\rightarrow\) Estimation \(\rightarrow\) Diagnostics \(\rightarrow\) Forecasting. The rest of the course runs on this template (Box & Jenkins, 1970).
Maximum likelihood answers “which parameters make the data most likely?” For pure AR(p), the conditional fit is OLS-like on the regression with lagged values. The exact/default arima() fit is usually numerical. For MA/ARMA, the optimizer is genuinely doing new work because past innovations are unobserved.
Information criteria: \(\text{IC} = -2 \log L + \text{penalty}\). Three main flavors — AIC (\(2k\), Akaike 1974), BIC (\(k \log T\), Schwarz 1978), AICc (small-sample correction, Hurvich & Tsai 1989). Smaller is better.
AIC overfits, BIC underfits. Neither is “right.” Report both. Use the disagreement as information. AIC is efficient for forecasting; BIC is consistent for selection.
auto.arima() is a very good intern — read its work, override it when you have reasons, and always confirm with diagnostics.
IC is not an oracle. The Monte Carlo showed that even with the true model in the candidate set and \(T = 500\), AIC missed it about a quarter of the time (and BIC about one in ten) — and both do worse in shorter samples. IC is one line of evidence among several.
Sign conventions, one more time. \(\Phi(L) = 1 - \sum \phi_j L^j\) (minus), \(\Theta(L) = 1 + \sum \theta_l L^l\) (plus). \(\epsilon_t\) for DGP innovations, \(e_t\) for residuals of an estimated model. When in doubt, consult notation_dictionary.md.
Module 4 hands Module 5 a fitted but unchecked model: the ARMA(1,2) on differenced UNRATE won every IC column, and its residual ACF is already whispering that something is missing at lag 12. Module 5 turns that whisper into a formal verdict:
After that, the Box-Jenkins loop is complete — identify, estimate, score, check — and Block 4 puts the surviving model to work forecasting.