library(ggplot2)
library(forecast)
source("R/arma_simulator.R")
source("R/model_selection.R")
full_mc <- read.csv("data/ic_mc_summary.csv", stringsAsFactors = FALSE)
unrate_raw <- read.csv("data/UNRATE.csv", stringsAsFactors = FALSE)
theme_set(theme_bw(base_size = 13))Module 4 Lab: Choosing and Fitting an ARMA Model
Econ 6376 — Applied Time Series Econometrics · AI Learning Pack
How to use this lab
Every Core chunk runs as shipped, offline. There are no blanks to fill. Your work is to predict, run, modify, and explain—especially when a score table does not crown the model you expected.
Write each prediction in LEARNING_LOG.md, a code comment, or on paper before you run. If an AI tutor is helping, it will pause once per work unit and offer two or three concrete predictions. You can answer or say “just run it”; either way, the code runs next.
The master equation is
\[ y_t=\alpha+\delta t+\sum_{j=1}^{p}\phi_jy_{t-j} +\sum_{l=1}^{q}\theta_l\epsilon_{t-l}+\epsilon_t. \]
Module 4 adds no new non-seasonal slider. Module 3 supplied AR and MA terms; today we choose \((p,q)\) and estimate their settings. The sequence follows the production deck—compression, selection, criteria, likelihood, Box-Jenkins, Monte Carlo, real data—while the canonical notes define every claim and the Core boundary.
1. Setup
arma_simulator() is the Module 3 DGP. safe_fit() and fit_is_valid() are the course’s shared screening rule: an optimizer error returns NULL, and a candidate survives only if it converged and all coefficients and scores are finite. We will not rebuild that machinery in each loop.
So what: the lab starts with one known DGP, one declared fit screen, and two documented caches—nothing depends on a network call or an assistant’s memory.
2. Isomorphism and the compression payoff
Module 3’s fingerprint table ends with an unresolved row:
| Process | ACF | PACF |
|---|---|---|
| AR(\(p\)) | tails off | cuts off at \(p\) |
| MA(\(q\)) | cuts off at \(q\) | tails off |
| ARMA(\(p,q\)) | tails off | tails off |
A stationary AR has an MA(\(\infty\)) representation; an invertible MA has an AR(\(\infty\)) representation. A finite ARMA is therefore compact code for dynamics that may need many terms in a pure representation. Compression is the payoff, and loss of a visual order cutoff is the price.
Predict first
The DGP is a zero-mean ARMA(1,1) with \(\phi=0.6\), \(\theta=0.4\), and \(T=500\). Will (a) the ACF cut off at 1, (b) the PACF cut off at 1, or (c) both tail off?
set.seed(1985)
y_arma <- arma_simulator(n = 500, phi = 0.6, theta = 0.4)
print(autoplot(y_arma) +
labs(title = "Known ARMA(1,1): neither print names the orders",
subtitle = "phi = 0.6, theta = 0.4, T = 500", x = "t", y = "y[t]"))
print(ggAcf(y_arma, lag.max = 20) + ggtitle("Sample ACF: tails off"))
print(ggPacf(y_arma, lag.max = 20) + ggtitle("Sample PACF: tails off"))
Now price the compression. Fit pure AR models from order 1 through 10, then compare their maximized log-likelihood with the two-coefficient ARMA(1,1).
Predict first
How many pure AR coefficients will it take to beat the two-coefficient ARMA’s in-sample log-likelihood: (a) 2, (b) about 5, or (c) more than 10?
ar_orders <- 1:10
ar_fits <- lapply(ar_orders, function(p) safe_fit(y_arma, p, 0))
stopifnot(all(vapply(ar_fits, fit_is_valid, logical(1))))
fit_arma11 <- safe_fit(y_arma, 1, 1)
compression <- data.frame(
model = c(paste0("AR(", ar_orders, ")"), "ARMA(1,1)"),
ar_ma_coefficients = c(ar_orders, 2),
logLik = c(vapply(ar_fits, function(f) f$loglik, numeric(1)),
fit_arma11$loglik)
)
compression$gap_from_ARMA11 <- compression$logLik - fit_arma11$loglik
compression$logLik <- round(compression$logLik, 3)
compression$gap_from_ARMA11 <- round(compression$gap_from_ARMA11, 3)
compression model ar_ma_coefficients logLik gap_from_ARMA11
1 AR(1) 1 -726.191 -20.888
2 AR(2) 2 -710.676 -5.373
3 AR(3) 3 -706.966 -1.663
4 AR(4) 4 -705.461 -0.158
5 AR(5) 5 -704.326 0.977
6 AR(6) 6 -704.158 1.145
7 AR(7) 7 -704.079 1.224
8 AR(8) 8 -703.423 1.880
9 AR(9) 9 -703.414 1.889
10 AR(10) 10 -702.722 2.581
11 ARMA(1,1) 2 -705.303 0.000AR(2) trails the ARMA by 5.37 log-likelihood points; AR(4) is still 0.16 behind; AR(5) finally clears it. AR(10), with five times as many AR/MA coefficients, buys only 2.58 more log-likelihood points. Because AR(\(p\)) nests inside AR(\(p+1\)), likelihood alone never tells us to stop adding lags.
So what: a small ARMA can compress long dynamics, but the representation is not visually unique—so model choice needs an explicit price for complexity.
3. Model selection is a general problem
Raw \(R^2\) cannot supply that price. Adding a regressor can only weakly reduce the residual sum of squares, so raw \(R^2\) can only weakly rise. Adjusted \(R^2\) adds a regression degrees-of-freedom correction, but it is not the likelihood-based ARMA workflow and carries neither AIC’s predictive-risk motivation nor BIC’s selection motivation.
Predict first
Four useless noise regressors are added to a correctly specified cross-section. Will raw \(R^2\) (a) fall, (b) stay the same or rise, or (c) become negative? And must adjusted \(R^2\) move the same way?
set.seed(44)
n <- 120
x1 <- rnorm(n)
noise <- replicate(4, rnorm(n))
y <- 1 + 2 * x1 + rnorm(n)
fit_small <- lm(y ~ x1)
fit_large <- lm(y ~ x1 + noise)
data.frame(
model = c("signal only", "signal + four noise regressors"),
R2 = c(summary(fit_small)$r.squared, summary(fit_large)$r.squared),
adjusted_R2 = c(summary(fit_small)$adj.r.squared,
summary(fit_large)$adj.r.squared)
) model R2 adjusted_R2
1 signal only 0.8434240 0.8420971
2 signal + four noise regressors 0.8507129 0.8441652The raw score rises even though the added variables are generated as noise. The adjusted score can decline, which is useful as a regression summary, but it still is not the score reported and motivated for our ARMA likelihood comparison. Core stops there; the deck’s deeper historical comparison of adjusted-\(R^2\) variants is not part of this lab.
So what: better in-sample fit is automatic when flexibility grows; a useful selection rule must say how much improvement an extra parameter has to buy.
4. Information criteria: declare the goal, then reveal the winner
All three Core criteria start from the same maximized log-likelihood:
\[\mathrm{IC}=-2\ell(\hat{\boldsymbol\vartheta})+\text{penalty}(k,T),\]
with smaller better. The parameter count \(k\) includes every estimated parameter—AR, MA, \(\sigma^2\), and any process mean.
\[ \mathrm{AIC}=-2\ell+2k, \qquad \mathrm{BIC}=-2\ell+k\log T, \]
\[ \mathrm{AICc}=\mathrm{AIC}+\frac{2k(k+1)}{T-k-1}. \]
Suppose \(T=200\). Candidate A has \(\ell=-100\), \(k=3\); candidate B fits better with \(\ell=-96\) but uses \(k=5\).
Commit to a goal before scoring
For forecasting, choose (a) AIC/AICc or (b) BIC as your primary criterion. For parsimonious explanation, choose again. Then predict whether the two criteria will crown the same candidate below.
hypothetical <- data.frame(
model = c("A: smaller", "B: better fit, larger"),
logLik = c(-100, -96),
k = c(3, 5)
)
T_ic <- 200
hypothetical$AIC <- -2 * hypothetical$logLik + 2 * hypothetical$k
hypothetical$AICc <- hypothetical$AIC +
2 * hypothetical$k * (hypothetical$k + 1) /
(T_ic - hypothetical$k - 1)
hypothetical$BIC <- -2 * hypothetical$logLik + hypothetical$k * log(T_ic)
hypothetical[c("logLik", "AIC", "AICc", "BIC")] <-
round(hypothetical[c("logLik", "AIC", "AICc", "BIC")], 3)
hypothetical model logLik k AIC AICc BIC
1 A: smaller -100 3 206 206.122 215.895
2 B: better fit, larger -96 5 202 202.309 218.492AIC and AICc accept B’s four-point deviance improvement after their lighter price; BIC keeps A. Neither criterion is “better” without a task. Use AIC or AICc for prediction; use BIC for explanation/parsimony when a sparse true model is plausibly in the set. Choose before seeing winners and report all three. The course rule of thumb prefers AICc to AIC when \(T/k<40\).
So what: an IC disagreement is not a malfunction—it identifies exactly where the conclusion depends on the chosen price of complexity.
5. Likelihood and the pure-AR bridge
Likelihood freezes the observed data and moves the parameters:
\[ L(\boldsymbol\vartheta\mid\mathbf y)=f_{\boldsymbol\vartheta}(\mathbf y), \qquad \hat{\boldsymbol\vartheta}_{\mathrm{MLE}} =\arg\max_{\boldsymbol\vartheta}\ell(\boldsymbol\vartheta). \]
It is a density evaluated at the fixed data, not the probability of an exact continuous sample and not a probability distribution over the parameter.
For a zero-mean pure AR, conditional Gaussian likelihood/CSS minimizes the same one-step squared residuals as OLS because every lagged \(y\) is observed.
Predict first
Fit the same zero-mean AR(1) by conditional CSS and no-intercept OLS. Will the \(\hat\phi\) values (a) agree to printed precision, (b) differ by about 0.1, or (c) have opposite signs?
set.seed(1985)
y_ar <- arma_simulator(n = 500, phi = 0.6)
fit_css <- stats::arima(y_ar, order = c(1, 0, 0),
include.mean = FALSE, method = "CSS")
fit_ols <- lm(as.numeric(y_ar[-1]) ~ as.numeric(y_ar[-length(y_ar)]) - 1)
fit_ml <- stats::arima(y_ar, order = c(1, 0, 0),
include.mean = FALSE, method = "CSS-ML")
data.frame(
method = c("conditional CSS", "no-intercept OLS", "default CSS-ML"),
phi_hat = c(coef(fit_css)["ar1"], coef(fit_ols)[1], coef(fit_ml)["ar1"])
) method phi_hat
1 conditional CSS 0.5573313
2 no-intercept OLS 0.5573313
3 default CSS-ML 0.5569450CSS and OLS agree: about 0.55733133. The numerical ML polish moves slightly to about 0.55694500. Say conditional pure-AR fitting is OLS-like, not “MLE is always OLS.” With MA terms, lagged innovations are unobserved and must be reconstructed recursively for each proposed parameter vector.
phi_grid <- seq(0.35, 0.76, length.out = 200)
y_now <- as.numeric(y_ar[-1])
y_lag <- as.numeric(y_ar[-length(y_ar)])
sse <- vapply(phi_grid, function(phi) sum((y_now - phi * y_lag)^2), numeric(1))
relative_ll <- -0.5 * length(y_now) * log(sse / min(sse))
ggplot(data.frame(phi = phi_grid, relative_ll), aes(phi, relative_ll)) +
geom_line(colour = "steelblue4", linewidth = 0.9) +
geom_vline(xintercept = coef(fit_ols)[1], linetype = "dashed") +
labs(title = "Minimizing conditional AR squared error maximizes likelihood",
subtitle = "Dashed line: the no-intercept OLS / CSS peak",
y = "relative conditional log-likelihood")
Fit and screen one known ARMA realization
Return to y_arma, the seed-1985 ARMA(1,1), and fit the full nine-candidate grid under the known zero-mean policy.
Predict first
The truth is in the grid. Will (a) every criterion choose (1,1), (b) AIC/AICc choose a larger model while BIC chooses (1,1), or (c) no criterion choose truth?
grid <- expand.grid(p = 0:2, q = 0:2)
one_fits <- lapply(seq_len(nrow(grid)), function(i) {
safe_fit(y_arma, grid$p[i], grid$q[i])
})
ok <- vapply(one_fits, fit_is_valid, logical(1))
stopifnot(all(ok))
one_table <- data.frame(
p = grid$p,
q = grid$q,
k = vapply(one_fits, function(f) length(f$coef) + 1L, integer(1)),
logLik = vapply(one_fits, function(f) f$loglik, numeric(1)),
AIC = vapply(one_fits, function(f) f$aic, numeric(1)),
AICc = vapply(one_fits, function(f) f$aicc, numeric(1)),
BIC = vapply(one_fits, function(f) f$bic, numeric(1))
)
one_table[order(one_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.479AIC and AICc choose ARMA(2,2); BIC chooses the true ARMA(1,1). A returned object still needs screening before ranking—the helper checks convergence and finite values rather than mistaking “R returned something” for success.
Compare the fitting interfaces under the same zero-mean policy:
manual_stats <- stats::arima(y_arma, order = c(1, 0, 1), include.mean = FALSE)
manual_forecast <- forecast::Arima(y_arma, order = c(1, 0, 1),
include.mean = FALSE)
automatic <- forecast::auto.arima(
y_arma, d = 0, seasonal = FALSE,
allowmean = FALSE, allowdrift = FALSE,
ic = "aicc", stepwise = TRUE, approximation = FALSE
)
data.frame(
fit = c("stats::arima ARMA(1,1)", "forecast::Arima ARMA(1,1)",
"auto.arima AICc search"),
model = c("(1,1)", "(1,1)",
paste0("(", arimaorder(automatic)["p"], ",",
arimaorder(automatic)["q"], ")")),
AICc = c(NA, manual_forecast$aicc, automatic$aicc)
) fit model AICc
1 stats::arima ARMA(1,1) (1,1) NA
2 forecast::Arima ARMA(1,1) (1,1) 1416.654
3 auto.arima AICc search (2,2) 1415.843On this realization, the automatic AICc search picks (2,2), AICc 1415.84, versus 1416.65 for the true (1,1). Automation is a labor-saver, not a conscience. For an undifferenced stationary fit, remember that R’s printed intercept is the process mean \(\hat\mu\), not recursive \(\hat\alpha\):
\[\hat\alpha=\hat\mu\left(1-\sum_j\hat\phi_j\right).\]
So what: likelihood supplies comparable fit, screening decides which fits are admissible, and the IC supplies a ranking—not proof that the top row fits.
6. Put the pieces in Box-Jenkins order
| Step | Output | Course location |
|---|---|---|
| Identification | \(d\) and a short \((p,q)\) candidate set | Modules 2–3 |
| Estimation | \(\hat{\boldsymbol\vartheta}\), uncertainty, \(\ell\) | Module 4 |
| IC scoring | ranking under a declared goal | Module 4 |
| Diagnostic checking | residual and fitted-root evidence | Module 5 |
| Forecasting | predictions and out-of-sample loss | Block 4 |
Identification produces a set, not a favorite. Write the set before fitting; hold the data, likelihood method, and mean policy fixed; screen every fit; rank under the chosen criterion. The winner then faces diagnostics. An order never fitted has no score, and a score winner has not yet passed a residual check.
So what: Box-Jenkins makes a model the output of a stated procedure rather than a sequence of undocumented tries until one answer looks comfortable.
7. Break the clean story: two Monte Carlo experiments
7.1 Verified full experiment—cached evidence
The shipped CSV records the complete experiment used in the canonical notes and production deck: seed 1999, 500 Gaussian ARMA(1,1) replications, \(T=500\), nine candidates, zero mean, and screened CSS-ML fits.
full_mc criterion model count percent
1 AIC (0,0) 0 0.0
2 AIC (1,0) 0 0.0
3 AIC (2,0) 29 5.8
4 AIC (0,1) 0 0.0
5 AIC (1,1) 372 74.4
6 AIC (2,1) 27 5.4
7 AIC (0,2) 0 0.0
8 AIC (1,2) 36 7.2
9 AIC (2,2) 36 7.2
10 BIC (0,0) 0 0.0
11 BIC (1,0) 0 0.0
12 BIC (2,0) 42 8.4
13 BIC (0,1) 0 0.0
14 BIC (1,1) 452 90.4
15 BIC (2,1) 3 0.6
16 BIC (0,2) 0 0.0
17 BIC (1,2) 3 0.6
18 BIC (2,2) 0 0.0subset(full_mc, model == "(1,1)") criterion model count percent
5 AIC (1,1) 372 74.4
14 BIC (1,1) 452 90.4AIC selects truth 372/500 = 74.4%; BIC selects truth 452/500 = 90.4%. The generator rejected 72 of 4,500 individual candidate fits, but every replication retained at least one candidate. This is verified cached evidence, not code that silently reruns during render.
7.2 Smaller live experiment—evidence you generate
The loop below uses the same DGP, grid, and screen, but only 50 replications and a different seed. A short sequential build probe found about 0.05 seconds per replication on the build machine; 50 keeps the default interactive on slower laptops. No parallel machinery is needed for this small run.
run_ic_experiment <- function(n_reps, seed) {
set.seed(seed)
mc_grid <- expand.grid(p = 0:2, q = 0:2)
pick_aic <- pick_bic <- rep(NA_character_, n_reps)
rejected <- 0L
for (rep_id in seq_len(n_reps)) {
y <- arma_simulator(n = 500, phi = 0.6, theta = 0.4)
fits <- lapply(seq_len(nrow(mc_grid)), function(i) {
safe_fit(y, mc_grid$p[i], mc_grid$q[i])
})
ok <- vapply(fits, fit_is_valid, logical(1))
rejected <- rejected + sum(!ok)
if (!any(ok)) next
aic <- vapply(fits[ok], function(f) f$aic, numeric(1))
bic <- vapply(fits[ok], function(f) f$bic, numeric(1))
pick_aic[rep_id] <- sprintf("(%d,%d)",
mc_grid$p[ok][which.min(aic)],
mc_grid$q[ok][which.min(aic)])
pick_bic[rep_id] <- sprintf("(%d,%d)",
mc_grid$p[ok][which.min(bic)],
mc_grid$q[ok][which.min(bic)])
}
levels <- sprintf("(%d,%d)", mc_grid$p, mc_grid$q)
summarize <- function(x, criterion) {
counts <- table(factor(x, levels = levels))
data.frame(criterion, model = names(counts),
count = as.integer(counts),
percent = 100 * as.integer(counts) / sum(counts))
}
list(
summary = rbind(summarize(pick_aic, "AIC"),
summarize(pick_bic, "BIC")),
rejected_candidate_fits = rejected,
all_failed_replications = sum(is.na(pick_aic))
)
}
Predict first
In 50 new replications, will (a) both truth rates exactly equal the cached 74.4% and 90.4%, (b) rates move but BIC still usually finds truth more often, or (c) neither criterion ever finds truth?
live_mc <- run_ic_experiment(n_reps = 50, seed = 4040)
live_mc$summary criterion model count percent
1 AIC (0,0) 0 0
2 AIC (1,0) 0 0
3 AIC (2,0) 3 6
4 AIC (0,1) 0 0
5 AIC (1,1) 37 74
6 AIC (2,1) 4 8
7 AIC (0,2) 0 0
8 AIC (1,2) 3 6
9 AIC (2,2) 3 6
10 BIC (0,0) 0 0
11 BIC (1,0) 0 0
12 BIC (2,0) 3 6
13 BIC (0,1) 0 0
14 BIC (1,1) 46 92
15 BIC (2,1) 1 2
16 BIC (0,2) 0 0
17 BIC (1,2) 0 0
18 BIC (2,2) 0 0c(rejected_candidate_fits = live_mc$rejected_candidate_fits,
all_failed_replications = live_mc$all_failed_replications)rejected_candidate_fits all_failed_replications
4 0 This run produces 74% truth selection by AIC and 92% by BIC, rejects four individual fits, and retains at least one fit in all 50 replications. The closeness to the cached rates is incidental. The durable pattern is that neither criterion is certain, AIC’s misses lean toward larger models, and BIC finds the sparse truth more often in this designed experiment.
So what: even with truth in the grid and a comfortable sample, IC has an error rate—so real-data scores are evidence to carry into diagnostics, not an oracle’s verdict.
8. Real data: the pre-COVID UNRATE candidate workflow
The primary design uses January 1960 through December 2019. April 2020 jumps from 4.4 to 14.8 and can flatten the full-sample correlograms by dominating the covariance scale. Outlier and break treatment is outside Module 4, so we state the pre-COVID window rather than silently discarding observations.
unrate_raw$observation_date <- as.Date(unrate_raw$observation_date)
unrate_df <- subset(
unrate_raw,
observation_date >= as.Date("1960-01-01") &
observation_date <= as.Date("2019-12-01")
)
stopifnot(nrow(unrate_df) == 720L, !anyNA(unrate_df$UNRATE))
unrate <- ts(unrate_df$UNRATE, start = c(1960, 1), frequency = 12)
d_unrate <- diff(unrate)
data.frame(levels = length(unrate), differences = length(d_unrate),
first = min(unrate_df$observation_date),
last = max(unrate_df$observation_date)) levels differences first last
1 720 719 1960-01-01 2019-12-01
Predict first
For first-differenced pre-COVID UNRATE, will (a) the ACF cut off cleanly, (b) the PACF cut off cleanly, or (c) both tail off and leave a low-order ARMA candidate set?
print(autoplot(d_unrate) +
labs(title = "First-differenced UNRATE, 1960–2019",
x = "year", y = "monthly change, percentage points"))
print(ggAcf(d_unrate, lag.max = 36) + ggtitle("Pre-COVID difference ACF"))
print(ggPacf(d_unrate, lag.max = 36) + ggtitle("Pre-COVID difference PACF"))
Both tail off. Modest annual-lag structure is visible, but UNRATE is already seasonally adjusted: call this residual annual-lag dependence, not raw calendar seasonality. It suggests later seasonal AR/MA candidates; it does not prove a seasonal unit root or justify \(D=1\).
Write the six candidates and mean policy before scoring:
\[ (p,q)\in\{(1,0),(0,1),(1,1),(2,1),(1,2),(2,2)\}, \qquad \mathbb E[\Delta\mathrm{UNRATE}_t]=0. \]
For this section’s explanation/parsimony goal, BIC is primary. Every candidate uses the same 719 differences, include.mean = FALSE, and CSS-ML method.
Commit before scoring
Which pattern do you expect: (a) all three criteria agree, (b) AIC/AICc choose a larger candidate than BIC, or (c) every fit fails? Record which candidate you expect to win under your primary criterion too.
unrate_orders <- data.frame(
p = c(1, 0, 1, 2, 1, 2),
q = c(0, 1, 1, 1, 2, 2)
)
unrate_fits <- lapply(seq_len(nrow(unrate_orders)), function(i) {
safe_fit(d_unrate, unrate_orders$p[i], unrate_orders$q[i])
})
stopifnot(all(vapply(unrate_fits, fit_is_valid, logical(1))))
unrate_ic <- data.frame(
model = sprintf("ARMA(%d,%d)", unrate_orders$p, unrate_orders$q),
AIC = vapply(unrate_fits, function(f) f$aic, numeric(1)),
AICc = vapply(unrate_fits, function(f) f$aicc, numeric(1)),
BIC = vapply(unrate_fits, function(f) f$bic, numeric(1))
)
unrate_ic[order(unrate_ic$BIC), ] model AIC AICc BIC
5 ARMA(1,2) -549.3327 -549.2767 -531.0212
4 ARMA(2,1) -545.2027 -545.1467 -526.8913
6 ARMA(2,2) -548.8514 -548.7672 -525.9621
3 ARMA(1,1) -517.8257 -517.7921 -504.0921
1 ARMA(1,0) -453.0474 -453.0306 -443.8917
2 ARMA(0,1) -450.6092 -450.5924 -441.4534All three criteria select ARMA(1,2): AIC \(=-549.333\), AICc \(=-549.277\), BIC \(=-531.021\). ARMA(2,2) is only 0.48 AIC points behind but 5.06 BIC points behind. Same fits; only the price changed.
Match the automatic search to the manual policy before comparing it:
unrate_auto <- auto.arima(
d_unrate, d = 0, seasonal = FALSE,
allowmean = FALSE, allowdrift = FALSE,
approximation = FALSE
)
c(model = paste0("ARMA(", arimaorder(unrate_auto)["p"], ",",
arimaorder(unrate_auto)["q"], ")"),
AICc = unrate_auto$aicc) model AICc
"ARMA(1,2)" "-549.276659094274" The comparable automatic search also chooses (1,2). Agreement names a working model; it does not finish Box-Jenkins.
working_fit <- unrate_fits[[which(unrate_orders$p == 1 & unrate_orders$q == 2)]]
ar_coef <- working_fit$coef[grep("^ar", names(working_fit$coef))]
ma_coef <- working_fit$coef[grep("^ma", names(working_fit$coef))]
ar_root_min <- min(Mod(polyroot(c(1, -ar_coef))))
ma_root_min <- min(Mod(polyroot(c(1, ma_coef))))
rho_e <- as.numeric(acf(residuals(working_fit), lag.max = 36,
plot = FALSE, na.action = na.pass)$acf)
data.frame(
phi1 = unname(ar_coef[1]),
theta1 = unname(ma_coef[1]),
theta2 = unname(ma_coef[2]),
convergence_code = working_fit$code,
smallest_AR_root_modulus = ar_root_min,
smallest_MA_root_modulus = ma_root_min,
residual_rho12 = rho_e[13],
residual_rho24 = rho_e[25],
residual_rho36 = rho_e[37]
) phi1 theta1 theta2 convergence_code smallest_AR_root_modulus
1 0.8638129 -0.8962667 0.2318646 0 1.157658
smallest_MA_root_modulus residual_rho12 residual_rho24 residual_rho36
1 2.076743 -0.1449465 -0.1251656 -0.1298691The fitted recursion is approximately
\[ \Delta y_t=0.8638\Delta y_{t-1}-0.8963e_{t-1}+0.2319e_{t-2}+e_t. \]
The optimizer converged and the smallest AR and MA root moduli are about 1.158 and 2.077. Those are sanity checks. Residual autocorrelations near \(-0.145\), \(-0.125\), and \(-0.130\) at lags 12, 24, and 36 remain beyond the approximate \(\pm0.073\) band. Formal residual diagnosis belongs to Module 5.
So what: the declared workflow gives us a plausible zero-mean ARMA(1,2) for pre-COVID changes in unemployment—and hands Module 5 a fitted model to test, not a conclusion to protect.
9. Explorations
Optional and outside the assessable surface
Nothing below is required, and none of it appears on problem sets or exams. Core ends with the pre-COVID working model. Use these as starting points for a question you actually care about.
1. Change the goal, not the answer after the fact. (guided)
Take a candidate table where AIC and BIC disagree. Write two one-paragraph recommendations: one for a forecasting client, one for a researcher seeking a compact mechanism. Keep the numerical evidence fixed and change only the goal.
2. Make AICc matter. (guided)
Repeat the known-DGP grid at \(T=60\) instead of 500. Compare AIC and AICc and explain the correction using \(T/k\), not the slogan “small sample.”
3. Regenerate the full 500×9 experiment. (guided)
This is executable but deliberately excluded from the default render. It runs 4,500 numerical fits and can take several minutes on a student laptop.
full_regenerated <- run_ic_experiment(n_reps = 500, seed = 1999)
full_regenerated$summary
full_regenerated$rejected_candidate_fits
full_regenerated$all_failed_replications
# Write to a new file; do not overwrite the verified shipped cache casually.
write.csv(full_regenerated$summary, "full_mc_regenerated.csv", row.names = FALSE)Compare to data/ic_mc_summary.csv. If results differ, first compare R version, optimizer behavior, helper stamp, seed, grid order, and screening—not just the headline percentages.
4. Put COVID back as a labeled sensitivity analysis.
Extend the sample beyond 2019, difference it, and compare ACF/PACF and the same six candidates. Do not describe the result as proof that either window is “right”; explain what one extreme episode changed and which missing method (outlier or break treatment) would be needed for a primary full-sample model.
5. Change the mean policy coherently.
Refit all six UNRATE candidates with an estimated process mean. Do not mix those scores with the zero-mean table. Did the winner change because of ARMA orders, or because every model gained a mean parameter?
6. Remove truth from the candidate set. (open)
Simulate the same ARMA(1,1) but forbid (1,1). Which approximation does each criterion prefer? This makes AIC’s “least wrong candidate” motivation literal.
10. Checks
R/checks.R is the referee. It independently checks the CSS/OLS bridge, IC arithmetic, the deterministic grid, fit screening, the full cache, the live Monte Carlo pattern, and the UNRATE design.
source("R/checks.R")Module 4 checks
---------------
PASS CSS AR(1) coefficient matches the no-intercept OLS slope [CSS = 0.557331331, OLS = 0.557331336, absolute gap = 4.7e-09]
PASS AIC, AICc, and BIC use the documented k = 3 arithmetic [k = 3, largest stored/formula gap = 0]
PASS Seed 1985 makes AIC/AICc choose (2,2) and BIC choose (1,1) [AIC (2,2), AICc (2,2), BIC (1,1)]
PASS safe_fit() enforces convergence, finite scores, and zero process mean [valid = TRUE; NULL/code/Inf rejected = TRUE/TRUE/TRUE; intercept absent = TRUE]
PASS The verified 500-replication Monte Carlo summary is complete and intact [AIC truth = 372/500; BIC truth = 452/500; rows = 18]
PASS The 50-replication live run has the documented broad IC behavior [AIC truth = 37/50; BIC truth = 46/50; rejected = 4; all failed = 0]
PASS The 1960-2019 zero-mean UNRATE grid is intact and selects ARMA(1,2) [levels = 720, changes = 719, winner rows = 5/5/5, AIC/AICc/BIC = -549.333/-549.277/-531.021]
---------------
7/7 checks passed.Where you are now
- Isomorphism explains why ARMA orders do not fall out of a correlogram and why compact models are valuable.
- Likelihood rewards fit; AIC, AICc, and BIC attach different prices to complexity.
- The modeling goal chooses the primary criterion before the result is visible.
- Conditional pure-AR CSS is OLS-like; MA terms require reconstructed innovations and numerical optimization.
- A fit must survive screening before it enters a ranking.
- The full and live Monte Carlos show error rates, not universal constants.
- Pre-COVID UNRATE leaves a provisional ARMA(1,2) with residual annual-lag dependence for Module 5 to test.
Before closing, record one place where your prediction and the result differed. The revised-reasoning and remaining-uncertainty fields belong in your own words.