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\)? It does adjust for degrees of freedom, but it was designed for linear-regression summaries, not for comparing likelihood-based ARMA specifications. It does not carry the KL-risk motivation of AIC or the model-selection motivation of BIC, and it is not the score reported by the ARMA workflow we use. Friends do not let friends rank ARMA candidates by adjusted \(R^2\) when AIC, AICc, BIC, diagnostics, and out-of-sample evaluation are available.
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.
Read those four steps as a procedure rather than as a description. Each one hands something specific to the next, and the reason for running all four in order is that the model you end up with is then the output of a stated process instead of a judgment call about which specification looked right.
Skipping a step breaks something specific. The common version is going from an ACF plot straight to auto.arima() and reporting whatever comes back: the candidate set was never stated, so nothing is known about what the winner beat, and the residuals were never inspected, so nothing is known about whether it fits. A ranking gets reported as a finished model.
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. Now ask, “which parameter values would have made this the most likely thing to happen?” Those are the maximum likelihood estimates.
Likelihood freezes the data and moves the parameters. The data you observed stay exactly where they are — nothing about them is being redrawn or reimagined. What moves is the parameter vector. You walk through candidate parameter values and, at each one, ask how much density the model places at the data you already have. The parameter values that place the most density there are the maximum likelihood estimates.
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:
intercept) — these are the MLE values of the parameters under R’s convention.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\).
R convention to memorize. For an undifferenced stationary ARMA fit, stats::arima() and forecast::Arima() label the estimated process mean \(\hat{\mu}\) as intercept. It is not the intercept \(\hat{\alpha}\) in the course’s recursive estimating equation. R parameterizes the fitted model as
\[y_t - \hat{\mu} = \sum_{j=1}^{p}\hat{\phi}_j(y_{t-j}-\hat{\mu}) + \sum_{l=1}^{q}\hat{\theta}_l e_{t-l} + e_t.\]
Expanding it gives the translation
\[\hat{\alpha} = \hat{\mu}\left(1-\sum_{j=1}^{p}\hat{\phi}_j\right).\]
So report the printed coefficient as \(\hat{\mu}\), and convert it before placing it into the master-equation recursion. For a zero-mean DGP, use include.mean = FALSE explicitly rather than relying on a default.
The generic structure of the classical likelihood-based criteria used in Core 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, R’s process-mean term) 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:
These are conditional theorems, not software warranties. AIC’s predictive-efficiency results and BIC’s selection-consistency results require standard likelihood regularity and appropriate candidate-set conditions; BIC additionally needs the true finite-dimensional model to be represented in the candidate set. Neither theorem guarantees that the selected model wins on the particular future sample you care about. The applied rule remains AIC/AICc for a prediction-first shortlist and BIC for a parsimony-first shortlist; Block 4 checks actual predictive performance with rolling-origin out-of-sample losses. Hamilton Chapter 5 supplies the likelihood regularity behind the fitting step, and Enders Chapter 2, pp. 79–87, develops forecast evaluation.
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, you have more confidence that the penalty choice is not driving the shortlist. When they disagree, AIC usually picked a bigger model whose extra parameters improved the likelihood enough to clear its prediction-oriented penalty but not BIC’s parsimony-oriented penalty. Whether those terms improve actual forecasts is checked out of sample in Block 4. The disagreement 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 | Posterior mean deviance + effective complexity | Bayesian hierarchical-model criterion; not the plug-in classical formula above |
| WAIC | Watanabe-Akaike IC | Pointwise predictive fit + effective complexity | Fully Bayesian predictive criterion; not a fixed-\(k\) penalty |
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. AIC, AICc, BIC, and HQIC share the simple maximized-log-likelihood-plus-penalty structure shown here. DIC and WAIC use posterior summaries and an effective complexity rather than merely inserting another fixed penalty into that classical formula. The names sound alike; the objects are not interchangeable.
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:
# The simulated DGP has mean zero, so every comparison uses that policy.
# Manual with stats::arima — minimal output
fit_m <- arima(y_arma, order = c(1, 0, 1), include.mean = FALSE)
# Manual with forecast::Arima — IC fields populated
fit_f <- Arima(y_arma, order = c(1, 0, 1), include.mean = FALSE)
# Automatic — let the algorithm choose p and q under the same zero-mean policy
fit_auto <- auto.arima(y_arma, ic = "aicc", stepwise = TRUE,
approximation = FALSE, d = 0, seasonal = FALSE,
allowmean = FALSE, allowdrift = FALSE)
summary(fit_f)Series: y_arma
ARIMA(1,0,1) with zero mean
Coefficients:
ar1 ma1
0.5398 0.4265
s.e. 0.0501 0.0574
sigma^2 = 0.9855: log likelihood = -705.3
AIC=1416.61 AICc=1416.65 BIC=1429.25
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.04983849 0.9907508 0.7826894 71.16302 195.8851 0.8997499
ACF1
Training set -0.01510177summary(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.002863087The mean policy is now held fixed: the manual fits use include.mean = FALSE, and the automatic fit uses allowmean = FALSE. On this seed, auto.arima() selects zero-mean ARMA(2,2), with AICc 1415.84 versus 1416.65 for the true ARMA(1,1). Rerun with a different seed and the answer can change; sometimes it finds (1,1), sometimes it does not. Automation is a labor-saver, not a conscience. It applies one IC along a particular search path. If you override its answer, 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),
include.mean = FALSE),
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 5 -702.8606 1415.721 1415.843 1436.794
6 2 1 4 -704.2464 1416.493 1416.574 1433.351
5 1 1 3 -705.3030 1416.606 1416.654 1429.250
8 1 2 4 -704.5013 1417.003 1417.083 1433.861
3 2 0 3 -710.6757 1427.351 1427.400 1439.995
7 0 2 3 -719.4077 1444.815 1444.864 1457.459
2 1 0 2 -726.1914 1456.383 1456.407 1464.812
4 0 1 2 -746.9366 1497.873 1497.897 1506.302
1 0 0 1 -914.6322 1831.264 1831.272 1835.479Look at the table this single 500-observation realization produced under the common zero-mean policy. AIC and AICc put ARMA(2,2) first; BIC puts the true ARMA(1,1) first by more than four points. The criteria disagree even though the truth is in the grid. One realization cannot tell us whether AIC’s miss is a fluke or a systematic tendency. 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. Step 7 waits for Module 5; Step 8 waits for Block 4 (Modules 6–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 <- rep(NA_character_, n_reps)
pick_bic <- rep(NA_character_, n_reps)
rejected_candidate_fits <- 0L
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 under the known zero-mean DGP policy.
# Errors return NULL; warnings are muffled here because validity is
# checked explicitly below using the optimizer code and finite IC.
fits <- lapply(seq_len(nrow(grid)), function(i) {
suppressWarnings(tryCatch(
Arima(y, order = c(grid$p[i], 0, grid$q[i]),
include.mean = FALSE),
error = function(e) NULL
))
})
# 3. A returned object is not enough. Keep a candidate only if the
# fit exists, both IC values are finite, and optim converged (code 0).
ok <- vapply(fits, function(f) {
!is.null(f) &&
isTRUE(f$code == 0) &&
all(is.finite(c(f$aic, f$bic)))
}, logical(1))
rejected_candidate_fits <- rejected_candidate_fits + sum(!ok)
# Guard the replication in the unlikely event that every candidate fails.
if (!any(ok)) next
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. Report fit screening, then tabulate successful replications.
cat("Rejected candidate fits:", rejected_candidate_fits, "of",
n_reps * nrow(grid), "\n")Rejected candidate fits: 72 of 4500 cat("Replications with no valid candidate:", sum(is.na(pick_aic)), "\n")Replications with no valid candidate: 0 table(pick_aic, useNA = "ifany")pick_aic
(1,1) (1,2) (2,0) (2,1) (2,2)
372 36 29 27 36 table(pick_bic, useNA = "ifany")pick_bic
(1,1) (1,2) (2,0) (2,1)
452 3 42 3 # 6. Express as percentages
round(100 * prop.table(table(pick_aic, useNA = "no")), 1)pick_aic
(1,1) (1,2) (2,0) (2,1) (2,2)
74.4 7.2 5.8 5.4 7.2 round(100 * prop.table(table(pick_bic, useNA = "no")), 1)pick_bic
(1,1) (1,2) (2,0) (2,1)
90.4 0.6 8.4 0.6 The exact numbers depend on the seed — the table below collects the percentages from the screened run above (seed 1999). The validity checks reject 72 of the 4,500 candidate fits, but every replication retains at least one converged candidate. Change the seed and the percentages move a little; the pattern is the point of the exercise:
| Model | AIC selects (%) | BIC selects (%) |
|---|---|---|
| (1,1) — the truth | 74.4 | 90.4 |
| (2,2) | 7.2 | 0.0 |
| (1,2) | 7.2 | 0.6 |
| (2,0) | 5.8 | 8.4 |
| (2,1) | 5.4 | 0.6 |
1. Neither scoring rule is a sure thing, even under laboratory conditions. AIC misses the true model 25.6% of the time and BIC 9.6% — 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 as the sample shrinks (Practice Problems 3 and 4 let you watch that happen). On real data, where the candidate models are approximations rather than known truths, IC is even less useful as a “truth finder.”
2. AIC systematically leans larger. Look at where the “mistakes” go. AIC selects (2,2), (1,2), or (2,1) — models with more parameters than the truth — in 19.8% of all replications. BIC makes those picks in only 1.2%. This is the finite-sample echo of the asymptotic discussion from Section 4.4.3: you can see the “AIC overfits” half of the slogan in the numbers. One honest wrinkle: at this sample size, BIC’s main miss is (2,0), an AR(2) with the same number of AR/MA coefficients as the truth that approximates its dynamics. The familiar shorthand describes a tendency, not where every miss must land; Practice Problem 4 asks you to examine a shorter sample.
3. The disagreement between AIC and BIC is informative. They pick different models in a nontrivial share of replications. 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 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:
Commit to the mean policy before scoring. For this exercise we impose a zero process mean on monthly changes: \(\mathbb{E}[\Delta y_t]=0\). Economically, that says the unemployment rate does not have a permanent deterministic monthly drift over this long window. It also makes the comparison fair: every manual candidate and the non-seasonal automatic search face the same restriction. In R that means include.mean = FALSE for Arima() and allowmean = FALSE for auto.arima(). This is a modeling decision, not a theorem; if your application implies a nonzero mean change, rerun every candidate with that alternative policy.
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, include.mean = FALSE)
})
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) -549.3327 -549.2767 -531.0212
(2,0,2) (2,0,2) -548.8514 -548.7672 -525.9621
(2,0,1) (2,0,1) -545.2027 -545.1467 -526.8913
(1,0,1) (1,0,1) -517.8257 -517.7921 -504.0921
(1,0,0) (1,0,0) -453.0474 -453.0306 -443.8917
(0,0,1) (0,0,1) -450.6092 -450.5924 -441.4534Under the stated zero-mean policy, all three criteria select ARMA(1,2). ARMA(2,2) sits only 0.48 AIC points behind, but 5.06 BIC points behind: the extra AR lag buys almost enough likelihood to clear AIC’s light hurdle, but not BIC’s heavier one. This is exactly why we committed to the mean restriction before looking at the scoreboard; otherwise an apparent order-selection result could actually be a hidden mean-term comparison.
Compare with auto.arima():
# Non-seasonal (our toolkit today), under the same zero-mean policy
auto.arima(d_unrate, d = 0, seasonal = FALSE,
allowmean = FALSE, allowdrift = FALSE,
approximation = FALSE)Series: d_unrate
ARIMA(1,0,2) with zero mean
Coefficients:
ar1 ma1 ma2
0.8638 -0.8963 0.2319
s.e. 0.0327 0.0473 0.0379
sigma^2 = 0.02707: log likelihood = 278.67
AIC=-549.33 AICc=-549.28 BIC=-531.02# Seasonal preview on the levels, also excluding deterministic drift
auto.arima(unrate, seasonal = TRUE,
allowmean = FALSE, allowdrift = FALSE,
approximation = FALSE)Series: unrate
ARIMA(1,1,2)(2,0,1)[12]
Coefficients:
ar1 ma1 ma2 sar1 sar2 sma1
0.9057 -0.9369 0.2187 0.4784 -0.0756 -0.7194
s.e. 0.0267 0.0441 0.0396 0.0867 0.0495 0.0794
sigma^2 = 0.02476: log likelihood = 310.83
AIC=-607.67 AICc=-607.51 BIC=-575.62The non-seasonal call is now comparable to the manual grid on the mean restriction: both model \(\Delta y_t\) with zero process mean, and both select ARMA(1,2). The seasonal call on levels selects ARIMA\((1,1,2)(2,0,1)_{12}\) without drift. That is a preview rather than an IC comparison with the differenced table: it is evidence to add seasonal AR/MA candidates to the candidate set, not proof that a seasonal difference is required.
Choose — publicly, with reasoning — a working candidate. For the course’s explanation/parsimony goal in this section, use the BIC winner under the zero-mean policy, provided its estimates and roots are sensible. Here that is ARMA(1,2), also selected by AIC, AICc, and the comparable automatic search. Hold it up as a working model, not a verdict: “This is our best guess today. It is not final. Next module 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 working ARMA has absorbed the short-run dynamics. But see that negative spike at lag 12, with an echo at lag 24? That is a hint of seasonal structure, not a verdict. Our non-seasonal ARMA may not explain it. That is the cliffhanger for Module 5: residual diagnostics will tell us whether the pattern is material, and seasonal modeling will separate two distinct possibilities — seasonal AR/MA dynamics versus a seasonal unit root. A lag-12 spike alone does not justify seasonal differencing; \(D=1\) needs unit-root evidence.
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 sign 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 and separates seasonal AR/MA dynamics from seasonal differencing. The latter is used only when evidence supports a seasonal unit root.
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 in the candidate set.
“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 must agree once the sample is reasonably large.” There is no universal sample-size or decimal-place rule. Their gap depends on persistence, proximity to nonstationary or noninvertible boundaries, initial conditions, and the sample itself. The default method = "CSS-ML" uses conditional sum of squares for starting values and then performs the ML optimization; use it as the course default, but inspect convergence and sensitivity when the fit is difficult.
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), include.mean = FALSE, method = "CSS") and with a no-intercept lm() (after constructing the lagged variable by hand), and verify that the coefficients agree. 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 ML fit under the same zero-mean policy and discuss the observed gap without invoking a fixed decimal-place rule.
Single-realization IC table. Simulate a zero-mean 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\}\) with include.mean = FALSE. 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 under the stated zero-mean policy and reproduce the IC table. Which candidate would you pick, and why? (Two to three sentences, including whether AIC and BIC agree.) Then refit the grid with a process mean and explain whether the policy changes the ranking.
Build your own grid search. Without using auto.arima(), write a function that takes a time series, maximum \(p\) and \(q\), and a mean-policy argument; fits the full \((p, q)\) grid with Arima(); rejects null fits, non-finite IC values, and optimizer codes other than zero; and returns a ranked IC table with AIC, AICc, and BIC columns. Test it on the simulated ARMA(1,1) from Problem 2 and compare its top-ranked model with auto.arima() under the same mean restriction. If they disagree, check the searched order set and stepwise path before deciding which result you believe.
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.
Classical likelihood 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 among models fitted to comparable data and likelihoods.
AIC tends to select larger models; BIC tends to select smaller ones. Neither is automatically “right.” Report both and use the disagreement as information. Under their respective regularity conditions, AIC has a prediction-efficiency result and BIC has a selection-consistency result; actual forecasts still face rolling-origin evaluation.
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 selected under an explicit mean policy. Its residual ACF is already whispering that something may be missing at lag 12. Module 5 tests that whisper rather than treating it as a verdict:
After that, the Box-Jenkins loop is complete — identify, estimate, score, check — and Block 4 puts the surviving model to work forecasting.