ARIMA Is Not Three Models Bolted Together, It Is One Autocovariance Read Six Different Ways

Everybody meets ARIMA as a shopping list of three hyperparameters. Pick p, pick d, pick q, call .fit(), ship the chart with the shaded band on it.

That framing is backwards. The three letters are the last thing you choose, and every one of them is read off a single object: the sample autocovariance γ_k. Identification reads it, estimation solves a system built out of it, the innovation variance falls out of the same numbers, and the forecast band is a recursion on the coefficients it produced. One bad index in γ_k poisons all three at once, silently.

I built the whole thing from scratch in one browser page, no libraries. Here is what matters in it.

Every model before this one could have had its rows shuffled

This is the first estimator in the series where the order of the file is the data. Everything from linear regression on Day 1 to the one-class SVM on Day 58 returns identical numbers with the rows shuffled, because they all assume independent draws. That assumption does enormous quiet work: it is what licences you to split randomly. A time series violates it deliberately, and a random split trains on Friday to predict Wednesday.

# the mistake that makes every time-series project look brilliant
X_tr, X_te = train_test_split(X, y, shuffle=True)   # NO. leaks the future.

# the only honest split, and the baseline you must beat
X_tr, y_tr = X[:t], y[:t]
X_te, y_te = X[t:t+h], y[t:t+h]      # refit as t advances (rolling origin)
yhat_naive = y[t-1]                  # "tomorrow equals today"

The dependence is not noise to be scrubbed out. It is the signal.

Stationarity is a licence, and the test points the other way

ARIMA estimates fixed numbers from history and carries them into the future, which is only defensible if the rules did not change while you were not looking. The augmented Dickey–Fuller test asks whether the level has any restoring force, by regressing the change on the level plus a few lagged changes.

// Delta y_t = alpha + rho*y_{t-1} + sum gamma_i * Delta y_{t-i} + e_t
const dy = diff1(y), rows = [], targ = [];
for (let t = lags; t < dy.length; t++) {
  const row = [1, y[t]];          // dy[t] = y[t+1]-y[t], so y[t] IS the lagged level
  for (let i = 1; i <= lags; i++) row.push(dy[t - i]);
  rows.push(row); targ.push(dy[t]);
}
const f = olsFit(rows, targ);     // se from s^2 * (X'X)^-1, computed here
return { stat: f.beta[1] / f.se[1], rho: f.beta[1], se: f.se[1], lags };
// MacKinnon critical values (constant, no trend): 1% -3.46  5% -2.88  10% -2.57

The standard error is a real one from s²(XᵀX)⁻¹ and the statistic is ρ̂/se(ρ̂). Read the direction carefully, because this is where results get misreported: ρ = 0 means a unit root, the null hypothesis is non-stationarity, and a statistic above about −2.88 fails to reject it. "Not rejected" is the default state of the world, not a finding. Those critical values are the only numbers on the page that are tabulated rather than computed.

The repair is the I, and it is one subtraction: ∇y_t = y_t − y_{t−1}, applied d times. Keep the first value of each round and the operation is exactly invertible, which makes the round trip the best correctness test in the pipeline. Two details bite. Reconstructing history cumulates from the first level while pushing a forecast back up cumulates from the last. And over-differencing leaves a fingerprint: it injects a spurious lag-1 autocorrelation of exactly −0.5 into an already stationary series.

One number, divided by n on purpose

Here is the object everything else reads from.

function acovf(y, K) {
  const n = y.length, m = mean(y), g = new Array(K + 1).fill(0);
  for (let k = 0; k <= K; k++) {
    let s = 0;
    for (let t = 0; t + k < n; t++) s += (y[t] - m) * (y[t + k] - m);
    g[k] = s / n;                   // divide by n, NOT by (n-k): keeps the matrix PSD
  }
  return g;
}
const acf     = (y, K) => { const g = acovf(y, K); return g.map(v => v / g[0]); };
const acfBand = n => 1.959963984540054 / Math.sqrt(n);

Only n−k products are summed and you divide by n anyway. That looks like a bug. It is a deliberate bias toward zero at long lags, and it is what keeps the covariance matrix positive semi-definite — which is what makes the Yule–Walker route below structurally incapable of returning a non-stationary model. The ±1.96/√n band comes from the fact that under white noise each ρ̂_k is roughly N(0, 1/n). It is a per-lag test, not a family-wise one, so with twenty lags plotted about one pokes outside by luck. A stray spike at lag 13 is not a seasonal discovery.

The PACF, or the same answer at O(K²) instead of O(K⁴)

In an AR(1) with φ = 0.8 the lag-2 autocorrelation is 0.64, and none of it is direct — it is the lag-1 effect applied twice. The partial autocorrelation strips the middlemen out, and it equals the last coefficient φ_kk of a fitted AR(k). Brute force means solving a fresh k×k Yule–Walker system for every k, at O(K⁴). Durbin–Levinson gets the identical answer by updating the previous solution in place.

function pacfDL(y, K) {
  const r = acf(y, K), pac = new Array(K + 1).fill(0); pac[0] = 1;
  let prev = [];                                   // phi_{k-1, 1..k-1}
  for (let k = 1; k <= K; k++) {
    let num = r[k], den = 1;
    for (let j = 1; j < k; j++) { num -= prev[j-1] * r[k-j]; den -= prev[j-1] * r[j]; }
    const pkk = num / den;                         // <-- the partial autocorrelation
    const cur = new Array(k); cur[k-1] = pkk;
    for (let j = 1; j < k; j++) cur[j-1] = prev[j-1] - pkk * prev[k-j-1];
    pac[k] = pkk; prev = cur;
  }
  return pac;
}

Read the two correlograms together and you get the Box–Jenkins table: two rules and three warnings. AR(p) tails off in the ACF and cuts off after lag p in the PACF, because once you have conditioned on p lags nothing direct is left. MA(q) is the mirror image, and its ACF is mathematically zero beyond lag q, because points further apart than q share no shocks. Everything inside the band means white noise, so stop. A near-linear decay still high at lag 20 means you forgot to difference. Neither plot cutting off means eyeballing has failed and you should run a grid.

Fit the AR part twice, and learn from the disagreement

Multiply the AR equation by y_{t−k}, take expectations, and the noise drops out. What is left is ρ_k = Σφ_iρ_{k−i} for k = 1…p, a system whose matrix is Toeplitz because entry (i,j) is ρ_{|i−j|}.

function yuleWalker(y, p) {                     // route 1: the Toeplitz system
  const g = acovf(y, p), r = g.map(v => v / g[0]), A = [], b = [];
  for (let i = 1; i <= p; i++) {
    const row = []; for (let j = 1; j <= p; j++) row.push(r[Math.abs(i - j)]);
    A.push(row); b.push(r[i]);
  }
  const phi = solveLin(A, b);
  let s = 0; for (let i = 1; i <= p; i++) s += phi[i-1] * r[i];
  return { phi, sigma2: g[0] * (1 - s), mu: mean(y) };
}
// route 2 stacks rows [1, y_{t-1}, ..., y_{t-p}] for t = p..n-1 and runs OLS.
// p = 1 has a closed form both routes must match: phi_hat = gamma_1 / gamma_0.

Fitting the same model two genuinely different ways is one of the most useful habits in numerical work: agreement is evidence, disagreement is a bug report. They differ in finite samples for a structural reason, and the page's table prints each route's row count so you can see why. Yule–Walker uses all n observations through γ_k and always returns a stationary model; OLS conditions on the first p points, throws those rows away, and can hand back an explosive root. Push n up and the gap shrinks roughly like 1/n. A gap that stays large means a near unit root, or a p too big for the sample.

The MA regressors do not exist

An MA(q) says y_t = c + e_t + θ₁e_{t−1} + … + θ_qe_{t−q}, and the comfortable regression story ends immediately, because the regressors are the innovations and innovations are not in your data. Full maximum likelihood answers that with a Kalman filter and a nonlinear optimiser. Hannan–Rissanen answers it with two ordinary regressions, which is why it runs live on a slider drag.

// STAGE 1 — a long AR whose residuals stand in for the unobservable shocks
const st1 = arOLS(z, mLong), eh = new Array(n).fill(0);
for (let t = mLong; t < n; t++) eh[t] = st1.fit.resid[t - mLong];

// STAGE 2 — now the MA columns EXIST, so it is just another regression
for (let t = mLong + Math.max(p, q); t < n; t++) {
  const row = [1];
  for (let i = 1; i <= p; i++) row.push(z[t - i]);       // AR columns
  for (let j = 1; j <= q; j++) row.push(eh[t - j]);      // MA columns (estimated!)
  rows.push(row); targ.push(z[t]);
}

Manufacture the missing column, then regress on it. It is consistent, needs no starting values and no optimiser, and in practice it seeds a proper likelihood fit. The cost is real: stage one burns mLong observations before anything is fitted at all, and that burn is exactly what wrecks the order-selection grid below.

The forecast is a recursion fed its own output

One step ahead is the fitted equation evaluated at the end of the data. Two steps ahead you no longer have y_{n+1}, so you substitute the forecast you just made, and future shocks are replaced by their expectation, which is zero.

for (let s = 0; s < h; s++) {
  let v = m.c;
  for (let i = 1; i <= p; i++) v += m.phi[i-1] * zz[n + s - i];   // may be a FORECAST
  for (let j = 1; j <= q; j++) { const idx = n + s - j;
    v += m.theta[j-1] * (idx < n ? e[idx] : 0); }                 // E[future shock] = 0
  zz.push(v); e.push(0); out.push(v);
}

That one substitution explains everything people find surprising here. An AR(1) gives ŷ_{n+h} = μ + φ^h(y_n − μ), so the path decays geometrically and a long-horizon forecast is just the mean. The MA part goes silent after exactly q steps, because beyond that every shock in the expression is a future one. And a random walk has φ = 1, nothing decays, and the optimal forecast is a flat line at the last value — not a failure of the method but the correct answer, and the reason the naive baseline is so hard to beat on financial data.

The band must widen, and for d > 0 you have to expand the operator first

Every stationary ARMA can be rewritten as an infinite moving average, y_t = μ + Σψ_je_{t−j}, where the ψ-weights are the impulse response: how much of one shock is still visible j steps later.

function psiWeights(phi, theta, H) {
  const psi = new Array(H + 1).fill(0); psi[0] = 1;
  for (let j = 1; j <= H; j++) {
    let v = j <= theta.length ? theta[j-1] : 0;
    for (let i = 1; i <= Math.min(j, phi.length); i++) v += phi[i-1] * psi[j - i];
    psi[j] = v;
  }
  return psi;
}
// se(h) = sqrt( sigma2 * sum_{j<h} psi_j^2 )      95% band = forecast +/- 1.96*se(h)
// for d > 0 feed it polyMulARd(phi, d): the operator phi(B) * (1-B)^d

The h-step forecast error is the sum of the shocks that have not happened yet, so Var(e_h) = σ²Σ_{j<h}ψ_j² — strictly increasing, which is why an honest band widens. For a stationary AR(1) the weights are φ^j and the sum converges to σ²/(1−φ²), so the band flattens at the unconditional variance. For a random walk every ψ_j = 1, the variance is σ²h, and the band fans out like √h forever.

Folding (1−B)^d into the AR operator is the step people skip. Compute ψ-weights from the fitted coefficients alone on a differenced model and they decay when they should not, so the intervals come out confidently, badly too narrow.

Choose the order honestly, then let the backtest overrule you

More parameters always fit the training data better, so residual variance cannot pick an order. AIC = −2ℓ + 2k and BIC = −2ℓ + k·log n add a penalty, BIC charging more per parameter as the sample grows. But the trap is not the formula.

let burn = 0;                              // the SAME startup for every candidate
for (let p = 0; p <= pmax; p++) for (let q = 0; q <= qmax; q++)
  burn = Math.max(burn, startupFor(z, p, q));
for (let p = 0; p <= pmax; p++) for (let q = 0; q <= qmax; q++) {
  const cs = css(z, hannanRissanen(z, p, q), burn);   // SSE over the COMMON window
  rows.push({ p, d, q, ...infoCriteria(cs.sse, cs.n, p + q + 2) });
}
// TRAP 1: AIC across different d is meaningless -- different data, different likelihood.
// TRAP 2: unequal effective samples make the whole table a lie. Force a common burn.
// TRAP 3: the grid minimum is an IN-SAMPLE choice. It is a shortlist, not a verdict.

That first loop is the point. A conditional fit burns startup observations, a bigger model burns more, and a grid that lets each candidate use whatever window it likes will crown whichever one got the easier subset. The live grid runs p from 0 to 3 and q from 0 to 2 and prints the effective sample in its own column.

Then trap three, the one that decides everything.

for (let t0 = minTrain; t0 + h <= y.length; t0++) {
  const train = y.slice(0, t0);                       // NOTHING after t0 is visible
  const fit = fitARIMA(train, p, d, q);               // REFIT from scratch every origin
  const fm = forecastARIMA(train, fit, h).mean;
  const fn = naiveForecast(train, h);                 // flat at the last value
  const fd = driftForecast(train, h);                 // straight line through the ends
  const fs = seasonalNaive(train, h, season);         // last season, repeated
}
// skill = 1 - rmse_model / rmse_naive.   <= 0 means the model added nothing.

Refit at every origin — not "fit once and slide the window", which quietly lets tomorrow's data into today's coefficients — and score against three baselines that cost nothing. Drift catches a trend the naive forecast misses; seasonal naive is savage on anything with a real period. The honest outcome is built into the demo: load the random-walk preset and the model loses every time, because there the naive forecast is provably optimal.

What the page actually checks

The maths sits between two markers in the file, ARIMA-ENGINE-START and ARIMA-ENGINE-END, with no DOM and no library calls between them, and a Node suite extracts that exact block from the shipped HTML and runs it against independently written baselines:

• differencing then integrating recovers the original series to 1e-9, for d = 1, 2 and 3
• pacfDL agrees with a fresh k×k Yule–Walker solve at every lag to 1e-12, and pacf[1] === acf[1] exactly
• acovf against an independent O(n²) loop, with ρ_0 = 1 and |ρ_k| ≤ 1 as free assertions
• the Gaussian elimination satisfies max|A·x − b| < 1e-10, and at p = 1 Yule–Walker returns exactly γ_1/γ_0
• the ψ-weights match an impulse-response simulation with e_0 = 1, and satisfy Σφ_iψ_{j−i} + θ_j − ψ_j = 0 term by term
• the closed forms: AR(1) forecasts against μ + φ^h(y_n − μ), AR(1) weights against φ^j, ARIMA(0,1,0) against ψ_j = 1 with Var = σ²h

That ψ-weight identity is my favourite. It is not a spot check against a magic number; it is φ(B)ψ(B) = θ(B), the defining polynomial equation of the model, asserted coefficient by coefficient.

The takeaway

One object explains the model. γ_k identifies the order through two correlograms, estimates the coefficients through a Toeplitz solve, sets the innovation variance, and feeds the recursion that produces both the forecast and the band around it. Durbin–Levinson, Yule–Walker, Hannan–Rissanen, the ψ-weights — each is another way of reading the same numbers.

ARIMA is not state of the art any more. It is still worth learning first, because stationarity, differencing, the correlograms, the widening band and rolling-origin evaluation are the vocabulary every newer method assumes you already have.

Drag the sliders, watch Yule–Walker and OLS converge as n grows, then try to beat the naive baseline on a random walk: https://dev48v.infy.uk/ml/day61-arima-time-series.html