Checking the Model, Seeing the Season

// Payload cell for makeSSI: owned by the including document because Quarto resolves
// FileAttachment relative to THIS file, not the partial (widgets/README.md).
ssiData = FileAttachment("cache/ssi_unratensa.json").json()

// Unique-id generator in its own cell: a factory may not reference its own name (that
// is a circular definition in OJS). Same idiom as Module 3's nextUid and Module 4's m4Uid.
misspecUid = (function () { let n = 0; return () => ++n; })()

// The numeric core. Everything between the BEGIN and END markers is byte-identical to
// slides/module_05/widgets/misspec_core.mjs (misspec_check.mjs asserts it). Edit one
// copy, then paste it over the other.
misspecCore = (function () {
// ===================== misspec core BEGIN =====================
// Numeric core for the Module 5 misspecification widget (makeMisspec; decisions 50, 56).
// Plain ES2017, no imports, no globals, no DOM. Course sign convention throughout
// (overview/notation_dictionary.md): eps_t is the DGP innovation, e_t the residual.
//
// Contents: RNG kit (mulberry32, hashSeed, normF); simARMA; sampleACF; solveLinear;
// olsLags; hannanRissanen; cssResid; cssJacobian; lmCSS; fitCSS (conditional sum of
// squares for ARMA(p, q) with a mean); gammln / gammaPQ / chisqUpper; ljungBox.

// --- RNG kit, ported verbatim from slides/module_03/module_03_slides.qmd (fpk) ---
const mulberry32 = a => () => {
  a |= 0; a = a + 0x6D2B79F5 | 0;
  let t = Math.imul(a ^ a >>> 15, 1 | a);
  t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
  return ((t ^ t >>> 14) >>> 0) / 4294967296;
};

const hashSeed = a => {
  a = a >>> 0;
  a = Math.imul(a ^ (a >>> 16), 2246822507);
  a = Math.imul(a ^ (a >>> 13), 3266489909);
  return (a ^ (a >>> 16)) >>> 0;
};

const normF = rng => () => {
  let u = 0, v = 0;
  while (u === 0) u = rng();
  while (v === 0) v = rng();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
};

// --- Simulator: Module 3's simARMA with optional single seasonal AR / MA coefficients at
// lag s (additive). The widget no longer uses the seasonal arguments (decision 56); with
// Phi = Theta = 0 this is exactly fpk.simARMA:
//   y_t = sum_j phi_j y_{t-j} + sum_l theta_l eps_{t-l} + eps_t,   eps_t ~ N(0, sigma^2)
// Zero-mean DGP; burn-in default 300.
function simARMA(phi, theta, T, seed, sigma, burn, Phi, Theta, s) {
  sigma = (sigma == null) ? 1 : sigma;
  burn  = (burn  == null) ? 300 : burn;
  Phi   = (Phi   == null) ? 0 : Phi;
  Theta = (Theta == null) ? 0 : Theta;
  s     = (s     == null) ? 12 : s;
  const nrm = normF(mulberry32(hashSeed(seed)));
  const p = phi.length, q = theta.length, n = T + burn;
  const eps = new Float64Array(n), y = new Float64Array(n);
  for (let t = 0; t < n; t++) eps[t] = sigma * nrm();
  for (let t = 0; t < n; t++) {
    let v = eps[t];
    for (let i = 1; i <= p; i++) if (t - i >= 0) v += phi[i - 1]   * y[t - i];
    for (let j = 1; j <= q; j++) if (t - j >= 0) v += theta[j - 1] * eps[t - j];
    if (Phi   !== 0 && t - s >= 0) v += Phi   * y[t - s];
    if (Theta !== 0 && t - s >= 0) v += Theta * eps[t - s];
    y[t] = v;
  }
  return Array.from(y.slice(burn));
}

// --- Sample ACF, divisor T (biased), demeaned by the sample mean — matches
// stats::acf(), which is what Box.test() calls.
function sampleACF(y, K) {
  const T = y.length;
  let m = 0; for (let i = 0; i < T; i++) m += y[i]; m /= T;
  const g = new Float64Array(K + 1);
  for (let k = 0; k <= K; k++) {
    let s = 0;
    for (let t = k; t < T; t++) s += (y[t] - m) * (y[t - k] - m);
    g[k] = s / T;
  }
  const r = new Array(K + 1);
  for (let k = 0; k <= K; k++) r[k] = g[0] === 0 ? 0 : g[k] / g[0];
  return r;
}

// --- Gaussian elimination with partial pivoting.
// A: array of n rows (each an array of n numbers), b: length-n array. Both are copied.
function solveLinear(A, b) {
  const n = b.length;
  const M = A.map((row, i) => row.slice().concat([b[i]]));
  for (let c = 0; c < n; c++) {
    let piv = c;
    for (let r = c + 1; r < n; r++) if (Math.abs(M[r][c]) > Math.abs(M[piv][c])) piv = r;
    if (Math.abs(M[piv][c]) < 1e-300) throw new Error("solveLinear: singular system");
    if (piv !== c) { const tmp = M[c]; M[c] = M[piv]; M[piv] = tmp; }
    for (let r = c + 1; r < n; r++) {
      const f = M[r][c] / M[c][c];
      if (f === 0) continue;
      for (let k = c; k <= n; k++) M[r][k] -= f * M[c][k];
    }
  }
  const x = new Array(n).fill(0);
  for (let r = n - 1; r >= 0; r--) {
    let s = M[r][n];
    for (let k = r + 1; k < n; k++) s -= M[r][k] * x[k];
    x[r] = s / M[r][r];
  }
  return x;
}

// --- OLS of y_t on an intercept and the lags listed in `lags` (e.g. [1, 2] for AR(2)).
// The first max(lags) observations are dropped, so the residual vector has n = T - max(lags)
// entries. Returns { coef, resid, n, lags } with coef[0] the intercept and coef[i] the
// coefficient on y_{t - lags[i-1]}. lags = [] is the intercept-only fit.
function olsLags(y, lags) {
  const T = y.length, m = lags.length, maxlag = m ? Math.max.apply(null, lags) : 0;
  const n = T - maxlag, k = m + 1;
  const XtX = Array.from({ length: k }, () => new Array(k).fill(0));
  const Xty = new Array(k).fill(0);
  const row = new Array(k);
  for (let t = maxlag; t < T; t++) {
    row[0] = 1;
    for (let i = 0; i < m; i++) row[i + 1] = y[t - lags[i]];
    for (let a = 0; a < k; a++) {
      Xty[a] += row[a] * y[t];
      for (let b = a; b < k; b++) XtX[a][b] += row[a] * row[b];
    }
  }
  for (let a = 0; a < k; a++) for (let b = 0; b < a; b++) XtX[a][b] = XtX[b][a];
  const coef = solveLinear(XtX, Xty);
  const resid = new Array(n);
  for (let t = maxlag; t < T; t++) {
    let f = coef[0];
    for (let i = 0; i < m; i++) f += coef[i + 1] * y[t - lags[i]];
    resid[t - maxlag] = y[t] - f;
  }
  return { coef, resid, n, lags: lags.slice() };
}

// --- Hannan-Rissanen start values for ARMA(p, q) with a mean. Step 1: a long AR(m) by OLS
// gives proxy innovations (zero before observation m). Step 2: OLS of y_t on a constant,
// p lags of y and q lags of the proxies over t >= m + q. Deterministic; m fixed (>= p).
function hannanRissanen(y, p, q, m) {
  const T = y.length;
  m = Math.max(m || 10, p);
  const longLags = []; for (let j = 1; j <= m; j++) longLags.push(j);
  const longAR = olsLags(y, longLags);
  const eh = new Array(T).fill(0);
  for (let t = m; t < T; t++) eh[t] = longAR.resid[t - m];
  const k = 1 + p + q, t0 = m + q;
  const XtX = Array.from({ length: k }, () => new Array(k).fill(0));
  const Xty = new Array(k).fill(0);
  const row = new Array(k);
  for (let t = t0; t < T; t++) {
    row[0] = 1;
    for (let j = 1; j <= p; j++) row[j] = y[t - j];
    for (let l = 1; l <= q; l++) row[p + l] = eh[t - l];
    for (let a = 0; a < k; a++) {
      Xty[a] += row[a] * y[t];
      for (let b = a; b < k; b++) XtX[a][b] += row[a] * row[b];
    }
  }
  for (let a = 0; a < k; a++) for (let b = 0; b < a; b++) XtX[a][b] = XtX[b][a];
  let coef;
  try { coef = solveLinear(XtX, Xty); } catch (err) { return null; }
  const phi = coef.slice(1, 1 + p), theta = coef.slice(1 + p);
  let sphi = 0; for (let j = 0; j < p; j++) sphi += phi[j];
  let ybar = 0; for (let t = 0; t < T; t++) ybar += y[t]; ybar /= T;
  const mu = Math.abs(1 - sphi) > 1e-8 ? coef[0] / (1 - sphi) : ybar;
  return { mu, phi, theta };
}

// --- Conditional-sum-of-squares residuals of ARMA(p, q) with mean mu:
//   e_t = (y_t - mu) - sum_j phi_j (y_{t-j} - mu) - sum_l theta_l e_{t-l},   t = p ... T-1,
// with pre-sample residuals zero (e_{t-l} = 0 whenever t - l < p). This is the recursion
// R's arima(method = "CSS") minimises (stats:::C_ARIMA_CSS with ncond = p), so the residual
// vector and the optimum coincide with R's. Returns { e, sse, n } with e of length T - p.
function cssResid(y, mu, phi, theta) {
  const T = y.length, p = phi.length, q = theta.length, n = T - p;
  const e = new Array(n);
  let sse = 0;
  for (let t = p; t < T; t++) {
    let v = y[t] - mu;
    for (let j = 1; j <= p; j++) v -= phi[j - 1] * (y[t - j] - mu);
    for (let l = 1; l <= q; l++) { const s = t - p - l; if (s >= 0) v -= theta[l - 1] * e[s]; }
    e[t - p] = v; sse += v * v;
  }
  return { e, sse, n };
}

// Analytic Jacobian of the residual vector (n x k, k = 1 + p + q), the same recursions
// differentiated; derivative of a pre-sample residual is zero:
//   de_t/dmu      = -(1 - sum_j phi_j) - sum_l theta_l de_{t-l}/dmu
//   de_t/dphi_j   = -(y_{t-j} - mu)    - sum_l theta_l de_{t-l}/dphi_j
//   de_t/dtheta_l = -e_{t-l}           - sum_m theta_m de_{t-m}/dtheta_l
function cssJacobian(y, mu, phi, theta, e) {
  const T = y.length, p = phi.length, q = theta.length, n = T - p, k = 1 + p + q;
  const J = Array.from({ length: n }, () => new Array(k).fill(0));
  let sphi = 0; for (let j = 0; j < p; j++) sphi += phi[j];
  for (let t = p; t < T; t++) {
    const r = t - p, row = J[r];
    row[0] = -(1 - sphi);
    for (let j = 1; j <= p; j++) row[j] = -(y[t - j] - mu);
    for (let l = 1; l <= q; l++) { const s = r - l; row[p + l] = s >= 0 ? -e[s] : 0; }
    for (let l = 1; l <= q; l++) {
      const s = r - l; if (s < 0) continue;
      const th = theta[l - 1], prev = J[s];
      for (let c = 0; c < k; c++) row[c] -= th * prev[c];
    }
  }
  return J;
}

// --- Levenberg-Marquardt on the CSS sum of squares from one start. Deterministic: fixed
// damping schedule (lambda / 10 on success, x 10 on failure), fixed iteration budget, no
// randomness. Stops when a step reduces the sum of squares by less than 1e-12 relative or
// moves no coefficient by more than 1e-10, or when no damping up to 1e12 yields a descent
// step (a stationary point to machine precision). Returns null if the start is unusable.
function lmCSS(y, start, maxIter) {
  maxIter = maxIter || 300;
  const p = start.phi.length, q = start.theta.length, k = 1 + p + q;
  const unpack = v => ({ mu: v[0], phi: v.slice(1, 1 + p), theta: v.slice(1 + p) });
  let x = [start.mu].concat(start.phi, start.theta);
  let cur = unpack(x), R = cssResid(y, cur.mu, cur.phi, cur.theta);
  if (!isFinite(R.sse)) return null;
  let lambda = 1e-3, iters = 0, converged = false, stationary = false;
  for (; iters < maxIter && !converged; iters++) {
    const J = cssJacobian(y, cur.mu, cur.phi, cur.theta, R.e);
    const A = Array.from({ length: k }, () => new Array(k).fill(0)), g = new Array(k).fill(0);
    for (let r = 0; r < R.n; r++) {
      const row = J[r], er = R.e[r];
      for (let a = 0; a < k; a++) {
        g[a] += row[a] * er;
        for (let b = a; b < k; b++) A[a][b] += row[a] * row[b];
      }
    }
    for (let a = 0; a < k; a++) for (let b = 0; b < a; b++) A[a][b] = A[b][a];
    let accepted = false;
    for (let tries = 0; tries < 40; tries++) {
      const M = A.map((row, i) => row.map((v, j) => i === j ? v * (1 + lambda) + 1e-12 : v));
      let d;
      try { d = solveLinear(M, g.map(v => -v)); } catch (err) { lambda *= 10; continue; }
      const xn = x.map((v, i) => v + d[i]);
      const nxt = unpack(xn), Rn = cssResid(y, nxt.mu, nxt.phi, nxt.theta);
      if (isFinite(Rn.sse) && Rn.sse < R.sse) {
        const drop = R.sse - Rn.sse;
        let step = 0; for (let i = 0; i < k; i++) step = Math.max(step, Math.abs(d[i]));
        x = xn; cur = nxt; R = Rn; lambda = Math.max(lambda / 10, 1e-15); accepted = true;
        if (drop <= 1e-12 * R.sse || step <= 1e-10) converged = true;
        break;
      }
      lambda *= 10;
      if (lambda > 1e12) break;
    }
    if (!accepted) { stationary = true; converged = true; }
  }
  return { mu: cur.mu, phi: cur.phi, theta: cur.theta, e: R.e, sse: R.sse, n: R.n, iters, converged, stationary };
}

// --- CSS fit of ARMA(p, q) with a mean.
//   q = 0: the problem is linear in (alpha, phi) with alpha = mu (1 - sum phi), so the fit is
//          plain OLS with intercept (olsLags) and mu = alpha / (1 - sum phi); residuals are
//          identical to the CSS recursion at that optimum.
//   q > 0: Levenberg-Marquardt from two deterministic starts, Hannan-Rissanen and R's own
//          start (mu = ybar, phi = theta = 0); the lower sum of squares wins.
// Returns { p, q, mu, phi, theta, resid, n, sse, iters, converged, start, method, starts }
// with resid of length T - p (t = p ... T-1) and n = T - p.
function fitCSS(y, p, q) {
  const T = y.length;
  let ybar = 0; for (let t = 0; t < T; t++) ybar += y[t]; ybar /= T;
  if (q === 0) {
    const lags = []; for (let j = 1; j <= p; j++) lags.push(j);
    const f = olsLags(y, lags);
    const phi = f.coef.slice(1);
    let sphi = 0; for (let j = 0; j < p; j++) sphi += phi[j];
    const mu = Math.abs(1 - sphi) > 1e-12 ? f.coef[0] / (1 - sphi) : ybar;
    let sse = 0; for (let i = 0; i < f.n; i++) sse += f.resid[i] * f.resid[i];
    return { p, q, mu, phi, theta: [], resid: f.resid, n: f.n, sse, iters: 0, converged: true,
             start: "ols", method: "ols", starts: [{ name: "ols", sse }] };
  }
  const zeros = m => { const v = new Array(m); for (let i = 0; i < m; i++) v[i] = 0; return v; };
  const starts = [];
  const hr = hannanRissanen(y, p, q);
  if (hr) starts.push({ name: "hr", mu: hr.mu, phi: hr.phi, theta: hr.theta });
  starts.push({ name: "zero", mu: ybar, phi: zeros(p), theta: zeros(q) });
  let best = null; const tried = [];
  for (let i = 0; i < starts.length; i++) {
    const s = starts[i], r = lmCSS(y, s);
    if (!r) { tried.push({ name: s.name, sse: NaN }); continue; }
    tried.push({ name: s.name, sse: r.sse, iters: r.iters, converged: r.converged });
    if (!best || r.sse < best.sse) { best = r; best.start = s.name; }
  }
  return { p, q, mu: best.mu, phi: best.phi, theta: best.theta, resid: best.e, n: best.n, sse: best.sse,
           iters: best.iters, converged: best.converged, start: best.start, method: "css-lm", starts: tried };
}

// --- Chi-square upper tail via the regularized incomplete gamma function.
// Standard numerical recipe: Lanczos ln-gamma, series for x < a + 1, continued
// fraction (modified Lentz) otherwise. Double-precision tolerances.
function gammln(xx) {
  const cof = [57.1562356658629235, -59.5979603554754912, 14.1360979747417471,
    -0.491913816097620199, 0.339946499848118887e-4, 0.465236289270485756e-4,
    -0.983744753048795646e-4, 0.158088703224912494e-3, -0.210264441724104883e-3,
    0.217439618115212643e-3, -0.164318106536763890e-3, 0.844182239838527433e-4,
    -0.261908384015814087e-4, 0.368991826595316234e-5];
  let x = xx, y = xx;
  let tmp = x + 5.24218750000000000;
  tmp = (x + 0.5) * Math.log(tmp) - tmp;
  let ser = 0.999999999999997092;
  for (let j = 0; j < 14; j++) ser += cof[j] / ++y;
  return tmp + Math.log(2.5066282746310005 * ser / x);
}

function gammaSeries(a, x) {          // P(a, x), converges fast for x < a + 1
  const gln = gammln(a);
  let ap = a, sum = 1 / a, del = sum;
  for (let n = 0; n < 5000; n++) {
    ap += 1; del *= x / ap; sum += del;
    if (Math.abs(del) < Math.abs(sum) * 1e-16) break;
  }
  return sum * Math.exp(-x + a * Math.log(x) - gln);
}

function gammaContFrac(a, x) {        // Q(a, x), converges fast for x >= a + 1
  const gln = gammln(a), FPMIN = 1e-300, EPS = 1e-16;
  let b = x + 1 - a, c = 1 / FPMIN, d = 1 / b, h = d;
  for (let i = 1; i < 5000; i++) {
    const an = -i * (i - a);
    b += 2;
    d = an * d + b; if (Math.abs(d) < FPMIN) d = FPMIN;
    c = b + an / c; if (Math.abs(c) < FPMIN) c = FPMIN;
    d = 1 / d;
    const del = d * c; h *= del;
    if (Math.abs(del - 1) < EPS) break;
  }
  return Math.exp(-x + a * Math.log(x) - gln) * h;
}

// Regularized lower (P) and upper (Q) incomplete gamma, P + Q = 1.
function gammaPQ(a, x) {
  if (!(x > 0)) return { P: 0, Q: 1 };
  if (x < a + 1) { const P = gammaSeries(a, x); return { P, Q: 1 - P }; }
  const Q = gammaContFrac(a, x); return { P: 1 - Q, Q };
}

// Upper-tail probability of a chi-square with df degrees of freedom: pchisq(x, df, lower.tail = FALSE).
function chisqUpper(x, df) { return gammaPQ(df / 2, x / 2).Q; }

// --- Ljung-Box on a residual vector: Q(h) = n(n+2) sum_{k=1}^{h} rho_k^2 / (n - k),
// compared with chi-square on h - fitdf degrees of freedom. Matches
// Box.test(e, lag = h, type = "Ljung-Box", fitdf = fitdf) term for term.
function ljungBox(e, h, fitdf) {
  fitdf = fitdf || 0;
  const n = e.length, r = sampleACF(e, h);
  let Q = 0;
  for (let k = 1; k <= h; k++) Q += r[k] * r[k] / (n - k);
  Q *= n * (n + 2);
  const df = h - fitdf;
  return { Q, df, p: df > 0 ? chisqUpper(Q, df) : NaN, h, fitdf, n };
}

// --- The widget's DGP menu (decision 56). Shared coefficients: phi = (0.5, 0.3) has real
// inverse AR roots 0.85 and -0.35; theta = (0.5, 0.4) has a complex pair of inverse MA roots
// of modulus 0.63, so the ARMA(2,2) has no near-cancelling AR/MA pair. Zero mean, sigma = 1.
// Fingerprints at T = 500: AR(2) PACF cuts off at 2 (rho_1 = 0.71, rho_2 = 0.66); MA(2) ACF
// cuts off at 2 (rho_1 = 0.50, rho_2 = 0.28); ARMA(2,2) tails off on both sides.
const MISSPEC_DGPS = {
  ar2:    { key: "ar2",    label: "AR(2)",     p: 2, q: 0, phi: [0.5, 0.3], theta: [],
            eq: "yₜ = 0.5·yₜ₋₁ + 0.3·yₜ₋₂ + εₜ" },
  ma2:    { key: "ma2",    label: "MA(2)",     p: 0, q: 2, phi: [],         theta: [0.5, 0.4],
            eq: "yₜ = εₜ + 0.5·εₜ₋₁ + 0.4·εₜ₋₂" },
  arma22: { key: "arma22", label: "ARMA(2,2)", p: 2, q: 2, phi: [0.5, 0.3], theta: [0.5, 0.4],
            eq: "yₜ = 0.5·yₜ₋₁ + 0.3·yₜ₋₂ + εₜ + 0.5·εₜ₋₁ + 0.4·εₜ₋₂" }
};
// Default base seed: the first base >= 6376 whose draw 1 tells every beat cleanly (see README).
const MISSPEC_DEFAULT_SEED = 6383;

const MISSPEC_CORE = { mulberry32, hashSeed, normF, simARMA, sampleACF, solveLinear, olsLags,
                       hannanRissanen, cssResid, cssJacobian, lmCSS, fitCSS,
                       gammln, gammaPQ, chisqUpper, ljungBox,
                       DGPS: MISSPEC_DGPS, DEFAULT_SEED: MISSPEC_DEFAULT_SEED };
// ===================== misspec core END =====================
return MISSPEC_CORE;
})()

makeMisspec = function (opts) {
  opts = opts || {};
  const C = misspecCore;
  const uid = misspecUid();
  const T = opts.T == null ? 500 : Number(opts.T);
  const K = opts.K == null ? 36 : Number(opts.K);
  const H = opts.h == null ? 24 : Number(opts.h);
  const PMAX = 3, QMAX = 3, BURN = 300, ALPHA = 0.05;
  const baseSeed = opts.seed0 == null ? C.DEFAULT_SEED : Number(opts.seed0) >>> 0;
  const DGPS = C.DGPS, KEYS = Object.keys(DGPS);
  const clampOrd = (v, dflt, mx) => v == null ? dflt : Math.max(0, Math.min(mx, Number(v) | 0));
  const defaults = {
    dgp: KEYS.indexOf(opts.dgp0) >= 0 ? opts.dgp0 : "arma22",
    p: clampOrd(opts.p0, 2, PMAX),
    q: clampOrd(opts.q0, 0, QMAX),
    draw: opts.draw0 == null ? 1 : Math.max(1, Number(opts.draw0) | 0)
  };
  let dgpKey = defaults.dgp, p = defaults.p, q = defaults.q, draw = defaults.draw;

  // Reseed steps a counter; the series seed is a hash of it, so every student sees the
  // same draw 1, draw 2, ... in the same order.
  const seedOf = d => C.hashSeed((baseSeed + d) >>> 0);

  // ---- tiny DOM kit (Module 4's m4k, inlined so the partial is self-contained) ----
  const NS = "http://www.w3.org/2000/svg";
  const el = (tag, cls, text) => { const n = document.createElement(tag); if (cls) n.className = cls; if (text != null) n.textContent = text; return n; };
  const sEl = (tag, attrs, text) => { const n = document.createElementNS(NS, tag); Object.keys(attrs || {}).forEach(k => n.setAttribute(k, String(attrs[k]))); if (text != null) n.textContent = text; return n; };
  const button = (label, action, cls) => { const b = el("button", "btn" + (cls ? " " + cls : ""), label); b.type = "button"; b.dataset.action = action; return b; };
  const clear = n => { while (n.firstChild) n.removeChild(n.firstChild); };
  const fmt = (x, d) => { const s = Number(x).toFixed(d == null ? 2 : d); return s === "-0." + "0".repeat(d == null ? 2 : d) ? s.slice(1) : s; };
  const SUB = ["₀", "₁", "₂", "₃", "₄", "₅", "₆", "₇", "₈", "₉"];
  const sub = n => String(n).split("").map(ch => SUB[+ch]).join("");
  const HAT = "̂";

  // ---- DOM ----
  const root = el("div", "widget card misspec");
  root.dataset.widget = "misspec"; root.dataset.uid = String(uid);

  const bar1 = el("div", "wbar");
  const grpD = el("div", "grp"); grpD.appendChild(el("span", "lab", "the data (DGP)"));
  const dgpBtn = {};
  KEYS.forEach(k => { const d = DGPS[k]; const b = button(d.label, "dgp-" + k); b.title = d.eq; dgpBtn[k] = b; grpD.appendChild(b); });
  const reseedB = button("Reseed", "reseed");
  const drawTag = el("span", "tag");
  bar1.append(grpD, el("div", "sep"), reseedB, drawTag);

  const bar2 = el("div", "wbar");
  const grpF = el("div", "grp"); grpF.appendChild(el("span", "lab", "your model"));
  const mkSlider = (label, mx, val, action) => {
    const ctl = el("label", "ctl"); ctl.append(label + " ");
    const s = document.createElement("input"); s.type = "range"; s.min = 0; s.max = mx; s.step = 1; s.value = val; s.dataset.action = action;
    const out = el("span", "cval"); ctl.append(s, out); return { ctl, s, out };
  };
  const pC = mkSlider("AR order p =", PMAX, p, "p");
  const qC = mkSlider("MA order q =", QMAX, q, "q");
  grpF.append(pC.ctl, qC.ctl);
  const resetB = button("Reset", "reset", "ghost");
  bar2.append(grpF, el("div", "sep"), resetB);

  const chart = sEl("svg", { width: 1120, height: 320, viewBox: "0 0 1120 320", role: "img",
    "aria-label": "Residual autocorrelation function to lag 36 with plus or minus two over root n bands" });
  const readout = el("div", "readout");
  const caveat = el("div", "caveat");
  caveat.innerHTML = "Fits here are conditional sum of squares; R’s <code>Arima()</code> uses maximum likelihood — the residual story is the same.";
  root.append(bar1, bar2, chart, readout, caveat);

  // ---- compute ----
  function state() {
    const d = DGPS[dgpKey];
    const seed = seedOf(draw);
    const y = C.simARMA(d.phi, d.theta, T, seed, 1, BURN);
    const fit = C.fitCSS(y, p, q);
    const r = C.sampleACF(fit.resid, K);
    const band = 2 / Math.sqrt(fit.n);
    const lb = C.ljungBox(fit.resid, H, p + q);
    const outs = []; for (let k = 1; k <= K; k++) if (Math.abs(r[k]) > band) outs.push(k);
    let kmax = 0; for (let i = 0; i < outs.length; i++) if (kmax === 0 || Math.abs(r[outs[i]]) > Math.abs(r[kmax])) kmax = outs[i];
    return { d, seed, fit, r, band, lb, outs, kmax };
  }

  // ---- chart ----
  function drawACF(st) {
    clear(chart);
    // Type inside the viewBox: every text element is >= 22 units. The deck caps the chart at
    // 1500 CSS px; reveal.js scales the 1920 canvas by (1 - margin) at 1080p, 0.98 for the
    // Module 5 deck (margin 0.02), so one viewBox unit is 1500 * 0.98 / 1120 = 1.3125 screen
    // px and 22 units render at 28.9 px. Even at Quarto's default margin 0.1 it is 26.5 px.
    const W = 1120, Hh = 320, L = 72, R = 26, Tm = 44, B = 46;
    const iw = W - L - R, ih = Hh - Tm - B, YMAX = 0.8;
    const clampY = v => Math.max(-YMAX, Math.min(YMAX, v));
    const sy = v => Tm + (YMAX - clampY(v)) / (2 * YMAX) * ih;
    const slot = iw / K, sx = k => L + slot * (k - 0.5);
    const bw = Math.min(16, slot * 0.5);
    const c = { ink: "#16202b", muted: "#5d6b78", line: "#d7dee5", soft: "#e7ecf1", shade: "#eef2f6",
                grey: "#8a94a0", blue: "#2c6fbb", orange: "#e4572e", navy: "#0b1f3a" };
    // band
    chart.appendChild(sEl("rect", { x: L, y: sy(st.band), width: iw, height: sy(-st.band) - sy(st.band), fill: c.shade }));
    [-0.5, 0.5].forEach(v => chart.appendChild(sEl("line", { x1: L, x2: L + iw, y1: sy(v), y2: sy(v), stroke: c.soft, "stroke-width": 1 })));
    [st.band, -st.band].forEach(v => chart.appendChild(sEl("line", { x1: L, x2: L + iw, y1: sy(v), y2: sy(v), stroke: c.grey, "stroke-width": 1.5, "stroke-dasharray": "7 6" })));
    // zero line
    chart.appendChild(sEl("line", { x1: L, x2: L + iw, y1: sy(0), y2: sy(0), stroke: c.line, "stroke-width": 1.5 }));
    // y labels
    [-0.5, 0, 0.5].forEach(v => chart.appendChild(sEl("text", { x: L - 12, y: sy(v) + 8, "text-anchor": "end", "font-size": 22, fill: c.muted }, fmt(v, 1))));
    // bars
    const outSet = {}; st.outs.forEach(k => { outSet[k] = true; });
    for (let k = 1; k <= K; k++) {
      const v = st.r[k], out = !!outSet[k];
      chart.appendChild(sEl("line", { x1: sx(k), x2: sx(k), y1: sy(0), y2: sy(v), stroke: out ? c.orange : c.blue,
        "stroke-width": bw, "stroke-linecap": "butt", "data-lag": k, "data-rho": v.toFixed(4) }));
      if (Math.abs(v) > YMAX) chart.appendChild(sEl("text", { x: sx(k), y: v > 0 ? Tm - 6 : Tm + ih + 20, "text-anchor": "middle", "font-size": 22, fill: c.orange }, "▲"));
    }
    // annotate the largest excursions (at most three); tall bars get the label beside the
    // tip so it never runs into the chart title or the x labels
    st.outs.slice().sort((a, b) => Math.abs(st.r[b]) - Math.abs(st.r[a])).slice(0, 3).forEach(k => {
      const v = st.r[k], tall = Math.abs(v) > 0.55;
      const attrs = tall
        ? { x: sx(k) + bw / 2 + 6, y: sy(v) + (v >= 0 ? 18 : -6), "text-anchor": "start" }
        : { x: sx(k), y: v >= 0 ? sy(v) - 8 : sy(v) + 24, "text-anchor": "middle" };
      attrs["font-size"] = 22; attrs.fill = c.orange; attrs["font-weight"] = 700;
      chart.appendChild(sEl("text", attrs, fmt(v, 2)));
    });
    // x labels; 1, 2, 3 are where AR/MA truncation shows, 12/24/36 stay marked for the hinge
    [1, 2, 3, 6, 12, 18, 24, 30, 36].filter(k => k <= K).forEach(k => {
      const strongLag = k <= 3 || k % 12 === 0;
      chart.appendChild(sEl("text", { x: sx(k), y: Hh - 12, "text-anchor": "middle", "font-size": strongLag ? 26 : 22,
        fill: strongLag ? c.navy : c.muted, "font-weight": strongLag ? 700 : 400 }, String(k)));
    });
    // titles
    chart.appendChild(sEl("text", { x: L, y: 28, "font-size": 26, fill: c.ink }, "residual ACF  ρ" + HAT + "ₑ(k) by lag k"));
    chart.appendChild(sEl("text", { x: L + iw, y: 28, "text-anchor": "end", "font-size": 22, fill: c.muted },
      "±2/√n = ±" + fmt(st.band, 3) + " (n = " + st.fit.n + ")"));
  }

  // ---- words ----
  function modelName(pp, qq) {
    if (pp === 0 && qq === 0) return "ARMA(0,0), mean only";
    if (qq === 0) return "AR(" + pp + ")";
    if (pp === 0) return "MA(" + qq + ")";
    return "ARMA(" + pp + "," + qq + ")";
  }
  // Every readout string is kept short enough to stay on ONE line at the deck's width
  // (about 150 characters at the deck's 0.82rem readout), so the widget height is fixed.
  function coefLine(st) {
    const parts = ["μ" + HAT + " = " + fmt(st.fit.mu, 2)];
    st.fit.phi.forEach((v, i) => parts.push("φ" + HAT + sub(i + 1) + " = " + fmt(v, 2)));
    st.fit.theta.forEach((v, i) => parts.push("θ" + HAT + sub(i + 1) + " = " + fmt(v, 2)));
    return parts.join(", ");
  }
  function fitNote(st) {
    if (st.fit.method === "ols") return "exact OLS";
    return "LM " + st.fit.iters + " steps, " + (st.fit.start === "hr" ? "HR" : "zero") + " start" +
      (st.fit.converged ? "" : ", budget reached (flat objective)");
  }
  // Verdict: Ljung-Box decides adequacy; the largest excursion is named. It says only that
  // structure is missing: which block is missing is read off the fitted orders against the
  // DGP's, so it is named in the why() line under the "DGP known here" label (in practice
  // the residual ACF/PACF shape only proposes a candidate).
  function missingBlock(st) {
    const underAR = p < st.d.p, underMA = q < st.d.q;
    return underAR && underMA ? "AR and MA structure" : underAR ? "AR structure" : underMA ? "MA structure" : null;
  }
  function verdict(st) {
    const d = st.d, rej = st.lb.p < ALPHA;
    const exact = p === d.p && q === d.q;
    const over = p >= d.p && q >= d.q && p + q > d.p + d.q;
    const largest = st.kmax ? ", largest at lag " + st.kmax + " (ρ" + HAT + "ₑ = " + fmt(st.r[st.kmax], 2) + ")" : "";
    const noise = st.outs.length ? " The excursion at lag " + st.outs.slice(0, 3).join(", ") + " is sampling noise (1 lag in 20)." : "";
    if (rej) {
      if (missingBlock(st)) return { warn: true, text: "Residuals show serial correlation" + largest + " → this fit is missing structure." };
      return { warn: true, text: "Serial correlation" + largest + " although the orders cover the DGP: a 5% false alarm or an off-optimum fit — reseed." };
    }
    if (over) return { warn: false, text: "Residuals look like white noise — adequate, but you are paying for parameters the data do not need (Module 4)." };
    if (exact) return { warn: false, text: "Residuals look like white noise — this fit is adequate." + noise };
    return { warn: false, text: "Residuals look like white noise — adequate though not the DGP’s orders: a finite AR or MA approximates the other side (duality)." };
  }
  // The why() line carries the "DGP known here" label, so it is where the missing block is named.
  function why(st) {
    const miss = missingBlock(st);
    return (miss ? "missing " + miss + " — " : "") + whyBody(st);
  }
  function whyBody(st) {
    const k = st.d.key;
    const flat = "; a cancelling root pair flattens the objective: estimates wander, residuals do not.";
    if (k === "ar2") {
      if (p < 2 && q === 0) return "an AR(2) needs two lags of y; AR(" + p + ") leaves the second lag’s memory in eₜ, so the residual ACF still decays.";
      if (p < 2) return "an AR(2) is an MA(∞) with slowly fading weights; MA(" + q + ") truncates that series, so structure remains (Module 3 duality).";
      if (q === 0) return p === 2 ? "AR(2) is the DGP; eₜ is εₜ up to estimation error." : "AR(2) nests inside AR(" + p + "); the extra lag estimates ≈ 0 and costs a parameter (Module 4).";
      return "the MA terms explain nothing" + flat;
    }
    if (k === "ma2") {
      if (p === 0 && q < 2) return "an MA(2) has exactly two non-zero autocorrelations; MA(" + q + ") leaves lag " + (q + 1) + " in eₜ.";
      if (q < 2) return "an MA(2) is an AR(∞) with weights fading like 0.63ᵏ; AR(" + p + ") truncates it, so the leftover shrinks as p grows (duality).";
      if (p === 0) return q === 2 ? "MA(2) is the DGP; eₜ is εₜ up to estimation error." : "MA(2) nests inside MA(" + q + "); the extra term estimates ≈ 0 and costs a parameter (Module 4).";
      return "the AR terms explain nothing" + flat;
    }
    if (q === 0) return "AR(" + p + ") truncates the AR(∞) form, leaving MA structure in eₜ at lags 1–2.";
    if (p === 0) return "MA(" + q + ") truncates the MA(∞) form, leaving slowly fading autocorrelation in eₜ.";
    if (p < 2 || q < 2) return "the leftover is the truncated remainder of the missing term (Module 3 duality).";
    if (p === 2 && q === 2) return "ARMA(2,2) is the DGP; eₜ is εₜ up to estimation error.";
    return "the extra terms explain nothing (≈ 0, a parameter each)" + flat;
  }

  // ---- render ----
  function render() {
    Object.keys(dgpBtn).forEach(k => dgpBtn[k].classList.toggle("on", k === dgpKey));
    pC.s.value = p; pC.out.textContent = String(p);
    qC.s.value = q; qC.out.textContent = String(q);
    drawTag.textContent = "draw " + draw;
    const st = state();
    drawACF(st);
    const v = verdict(st);
    readout.innerHTML =
      "<div><b>Fitted:</b> " + modelName(p, q) + " by CSS · " + coefLine(st) +
      " · n = " + st.fit.n + " · " + fitNote(st) + " · DGP " + st.d.label + ", draw " + draw + "</div>" +
      "<div><b>Ljung-Box</b> Q(" + H + ") = <b>" + fmt(st.lb.Q, 2) + "</b> · h = " + H + ", fitdf = " + (p + q) +
      " · χ²" + sub(st.lb.df) + " · p = <b>" + (st.lb.p < 0.0005 ? "&lt; 0.001" : fmt(st.lb.p, 3)) + "</b> " +
      "<span class=\"tag " + (st.lb.p < ALPHA ? "warn" : "good") + "\">" + (st.lb.p < ALPHA ? "reject at 5%: serial correlation" : "fail to reject at 5%") + "</span></div>" +
      "<div class=\"verdict\" data-verdict=\"" + (v.warn ? "warn" : "ok") + "\">" + v.text + "</div>" +
      "<div class=\"why\">Why (DGP known here, not in practice): " + why(st) + "</div>";
  }

  // ---- wire ----
  Object.keys(dgpBtn).forEach(k => dgpBtn[k].addEventListener("click", () => { dgpKey = k; render(); }));
  pC.s.addEventListener("input", () => { p = Number(pC.s.value) | 0; render(); });
  qC.s.addEventListener("input", () => { q = Number(qC.s.value) | 0; render(); });
  reseedB.addEventListener("click", () => { draw += 1; render(); });
  resetB.addEventListener("click", () => { dgpKey = defaults.dgp; p = defaults.p; q = defaults.q; draw = defaults.draw; render(); });
  render();
  return root;
}