Assignment 1 · Question 2 · AR(p) order selection

How many lags is enough?

Two autoregressive models generate the data. M1 has five lags with decaying weights and M2 has two strong ones. The task was to fit AR(p) for p = 1…10 by least squares and let three information criteria choose p, first on single samples, then over 1,000 simulated series at n = 100 and n = 15, and finally through closed-form overfitting probabilities.

IC₁ · AIC-like

lo g σ^_{p}^{2} + \frac{2 ( p + 1 )}{T}

IC₂ · AICc-like

lo g σ^_{p}^{2} + \frac{T + p}{T - p - 2}

IC₃ · BIC-like

lo g σ^_{p}^{2} + \frac{p lo g T}{T}

with $σ^_{p}^{2} = RSS_{p} / T$ and $T = n - p$ usable observations; M1: $ϕ = (0.434, 0.217, 0.145, 0.108, 0.087)$ , M2: $ϕ = (0.682, 0.346)$ , $η_{t} \sim N (0, 1)$ .

[Q2.4–2.5] IC(p = 1:10, Y = y1, n = 100)

One sample of each model, three criteria

Simulate 100 observations from M1 and M2, fit AR(p) by least squares for p = 1…10 and evaluate each criterion. The arg-min is the selected order.

sample seed

Simulated series (salary scale, $000s)

y₁ from M1y₂ from M2

submitted code (Q2.4)

y1[6:100] <- sapply(6:100, function(t) y[(t-1):(t-length(phi.m1)) %*% phi.m1 + rnorm(1)])
y2[3:100] <- sapply(3:100, function(t) y[(t-1):(t-length(phi.m2))] %*% phi.m2 + rnorm(1))

M1 · AR(5) sample

M2 · AR(2) sample

IC₁ · AIC-likeIC₂ · AICc-likeIC₃ · BIC-like

Selected order for M1 (true p = 5)

IC₁ →

p^

= 3IC₂ →

p^

= 3IC₃ →

p^

= 2

Selected order for M2 (true p = 2)

IC₁ →

p^

= 7IC₂ →

p^

= 3IC₃ →

p^

= 3

Known bug in the submission: the series were built from baseball salaries

Question 2 reused y, which still held the Hitters salary vector from Question 1. For M1 the noisy linear combination of time indices became an index into the salaries; for M2 the series is a weighted sum of two neighbouring salaries. That is why the 2023 plots needed ylim = c(11.5, 14): $lo g σ^^{2}$ is on a log-salary scale. With seed 10 this view reproduces the printed tables to five decimals.

The Monte Carlo study below was unaffected: its simulation function generated each series correctly.

[Q2.6–2.9] simulate_and_evaluate_IC(n, phi, num_simulations = 1000)

Monte Carlo: how often does each criterion pick the true order?

Generate many series from the true model, select p by each criterion, and count. Runs in a Web Worker on R's Mersenne-Twister stream: seed 10 with M1 and n = 100 reproduces Q2.6. The four 2023 tables were drawn one after another from a single stream, so use Replay the 2023 run to reproduce all of them.

replicates

Monte Carlo seed

M1 · AR(5), n = 100submitted 2023 output · seed 10

IC₁ · AIC-likeIC₂ · AICc-likeIC₃ · BIC-like

Selection rates (true p = 5)
Criterion	under	correct	over
IC₁ · AIC-like	76.9%	7.8%6.3–9.6	15.3%
IC₂ · AICc-like	82.0%	7.8%6.3–9.6	10.2%
IC₃ · BIC-like	97.8%	1.4%0.8–2.3	0.8%

Under each correct rate: its Wilson 95% interval. With R = 1,000 replicates the Monte Carlo standard error of a rate is at most 1.6 points (here IC₁ 0.8, IC₂ 0.8, IC₃ 0.4).

What the submission observed (Q2.9)

With n = 100 the BIC-like IC₃ concentrates near the true order while IC₁ and IC₂ spread into larger p. With n = 15 every criterion overfits to p ≥ 7, almost always to p ≥ 8: with eight or more coefficients only five to seven rows remain, so the residual variance collapses towards zero.

2026 note, with intervals on the same counts. The reading holds for M2: IC₃ picks p = 2 in 79.8% (77.2–82.2%) of replicates. It does not hold for M1. There IC₃ under-fits in 97.8% of replicates and is the criterion least likely to pick p = 5, at 1.4% (0.8–2.3%), against 7.8% (6.3–9.6%) for IC₁. For M1, “near the true order” meant p = 2 or 3, below it.

Convergence: P(selects true p) as replicates accumulate

Running estimate after r replicates with its Wilson 95% band. The band narrows like 1/√r, so the last few hundred replicates move the estimate by well under a percentage point.

IC₁ · AIC-likeIC₂ · AICc-likeIC₃ · BIC-like

P(selects true p) with 95% intervals

Wilson intervals from 1,000 replicates. Each interval describes one criterion on its own. The criteria are scored on the same replicates, so compare them with the paired test below rather than by whether these intervals overlap.

Paired comparison on the same replicates

Only replicates where exactly one of the two criteria finds the true p carry information about the difference. Exact McNemar test on those; difference in P(selects true p) with Newcombe's paired 95% interval.

Pair	only A / B	A − B, points	McNemar
IC₁ vs IC₂	9 / 9	0.0−0.9 to +0.9	p = 1.00
IC₂ vs IC₃	66 / 2	+6.4+4.9 to +8.1	p < 0.001
IC₁ vs IC₃	67 / 3	+6.4+4.9 to +8.1	p < 0.001

Explain this simulationoptional · your own key

Add your own Anthropic or OpenAI key in AI settings to enable. Nothing is sent without one.

[Q2.10–2.13] prob_IC1(n, p0, L), prob_IC2, prob_IC3

Probability of overfitting by L extra lags

The submission derived closed forms for P(IC at p₀ + L beats IC at p₀) and tabulated them for n = 25 and n = 100. Below: those formulas evaluated by the TypeScript port (identical to R's pf to 10 decimals), next to a direct simulation.

As submitted (closed form, Q2.11)p₀ = 8 for n = 25 and p₀ = 3 for n = 100, as coded
n	L	IC₁	IC₂	IC₃
25	1	0.09880	0.26448	0.10639
	2	0.07109	0.87001	0.08197
	3	0.10591	1.00000	0.12734
	4	0.22782	1.00000	0.27063
	5	0.50174	0.00000	0.55469
	6	0.82707	0.00000	0.84675
	7	0.97455	0.00000	0.97272
	8	0.99818	0.00000	0.99682
100	1	0.01451	0.01915	0.01983
	2	0.00135	0.00240	0.00255
	3	0.00024	0.00059	0.00063
	4	0.00006	0.00022	0.00024
	5	0.00002	0.00011	0.00012
	6	0.00001	0.00008	0.00008
	7	0.00001	0.00007	0.00006
	8	0.00000	0.00007	0.00006

Estimate the same probabilities directly: simulate M1 1,000 times for each n and count how often the criterion at p₀ + L beats p₀.

prob_IC*: port of the submitted R functionspf(): regularised incomplete beta

Reading the submitted table with care

Overfitting means RSS_p₀/RSS_p₀+L exceeds a threshold, an upper-tail event, but the code evaluates pf(threshold, …), the lower tail, with the scale factor L/(n − p₀ + L) where the F-statistic needs (n − p₀ − L)/L. The IC₂ penalty is also mis-simplified, as (n + 2L)/(n − 2p₀ − 2L) instead of n/(n − 2p₀ − 2L − 2). That is why IC₂ jumps between exactly 1 and exactly 0 at n = 25.

The simulation shows the probabilities under M1. IC₃'s chance of adding lags shrinks quickly as n grows, while IC₁ and IC₂ keep a sizeable chance of adding one lag even at n = 100. That is the classic contrast between AIC-type and BIC-type penalties, so the submitted limit of 0 (Q2.13) holds only for IC₃.