Resampling and time-series conventions

Part of the numerical-validation suite. The recorded conventions behind permutation tests, VAR / Granger causality and the LOESS smoother.

LOESS smoother on the scatterplot — the v1 conventions

Decided 2026-08-18, before any code, as sub-task 6a of docs/releases/release_0.19.0.md item 6 (kept private). The smoother gives "is this relationship non-linear?" a visual answer the linear-only fit line cannot. Every convention below was verified against R loess (stats, R 4.5.3) by reproducing its predictions from the hand formulas at machine precision (≤ 6e-14) on the cars data (n = 50, tied x values, spans 0.3–1.0), a 6-level × 4-replicate ties set, an n = 8 set, an n = 200 uniform set and n = 4…7 straight lines — not recorded from memory. The pinned values live in tool/validation/generate_loess_references.R and test/validation/loess_reference_values.json (sub-task 6f).

The shape

Two new arguments on the existing scatter tool and analyze kind — smoother (boolean, default false) and span (number in (0, 1], default 0.75, R's own default) — catalog +0. A span outside (0, 1] errors as data, and so does span given without smoother=true (an argument that would be silently ignored is an error, not a shrug). The curve is computed engine-side (the numeric core beside the chart kinds in rust/src/analyze/graphs.rs) and emitted as its own result table — internal chart geometry on the violin-density precedent — which the scatter painter draws as a curve over the points. One engine call produces points and curve together, so the journaled command reproduces both.

The estimator (the verified reference algorithm)

The reference is R loess(y ~ x, span = s, degree = 1, family = "gaussian", surface = "direct") with the default loess.control — an exact local-linear tricube fit, no interpolation surface. The hand algorithm that reproduces it at machine precision:

  1. Neighborhood size q = ⌊span·n + 10⁻⁷⌋, clamped to [2, n]. The 10⁻⁷ is load-bearing and verified: 0.6·50 is 29.999999… in binary, and R rounds it up to 30.
  2. At each evaluation point x₀: distances dᵢ = |xᵢ − x₀|, bandwidth h = the q-th smallest distance, tricube weights over all n points wᵢ = (1 − uᵢ³)³ for uᵢ = dᵢ/h < 1, else 0 — with tied distances more than q points can fall inside the bandwidth (all are weighted), and the q-th-nearest point itself gets weight 0.
  3. Weighted local line: fitted = ȳ_w + β̂ (x₀ − x̄_w) with β̂ = Σw(x−x̄_w)(y−ȳ_w) / Σw(x−x̄_w)².

Degree 1 is deliberate (R loess defaults to degree 2): the row's recorded purpose is a trend-reading aid, and the local-linear fit is the Cleveland lowess convention the exploratory-graphics audience knows. family = "gaussian" means no robustness iterations (a recorded cut).

The curve is evaluated on a fixed grid of 101 equally spaced points over [min x, max x] of the plotted (complete-pair, filtered) points — endpoints included, each split group over its own range — dense enough that the polyline is smooth at chart resolution, and the same grid shape the pins use.

Degenerate windows (recorded divergence)

Where the local fit is singular, R's implementation takes a pseudoinverse path that warns and prints artifacts — verified: on a 10-fold mass of ties at x = 5 (neighborhood mean 25.1) loess predicts 0.0. That is deliberately not copied. ChakataStat's recorded behaviour: h = 0 (the q-th nearest distance is zero) fits the unweighted mean of the points at distance zero; a singular window (weighted x-variance 0, e.g. every in-bandwidth point at one x) fits the tricube-weighted local mean. Both are well-defined limits of the estimator; the pins avoid singular windows, and the parity sweeps above all ran on non-degenerate data.

The other R smoother (lowess) — recorded

stats::lowess with iter = 0, delta = 0 coincides with the pinned algorithm exactly (verified: 1.7e-13 on cars) — the two share Cleveland's estimator. Its defaults differ twice and are the recorded alternative, not what we draw: iter = 3 robustness iterations (up to 4.5 y-units apart on cars) and delta interpolation between anchor points. One more reference quirk, recorded so nobody trips on it when regenerating: loess.control(statistics = "exact") changes R's own predictions (up to ~1.0 on cars) — the pins are against the documented default control only.

Small n and the omission note

The smoother is drawn only when the plotted points have at least 6 complete pairs and at least 2 distinct x values. Below six, R's own reference implementation distrusts its fits (it refuses n = 3 — "span too small" — and takes the warning pseudoinverse path at n = 4–5, even on a straight line); a "curve" through fewer points is decoration. The scatter itself still draws: the points table gains a note naming what was omitted and why ("LOESS smoother not drawn: …"), the omitted-not-failed convention. This is a data condition; a bad span argument stays an error as data.

Presentation

  • The smoother table"LOESS Smoother: <points-table title>", columns the x variable's name and Fitted, 101 rows; its note names the estimator, the span and the pair count ("Local-linear LOESS (tricube), span 0.75, over N complete pairs."). Like the violin density tables it is chart geometry: consumed by the chart layer, not shown as an output table.
  • The curve keeps the theme accent and the single-color colors override recolors the points only — the Bland–Altman rule: the annotation keeps its own accent so it stays readable over recolored data. The curve participates in the y-axis domain (a local line can overshoot the data range slightly; it is never clipped away).
  • The chart editor's scatter surface stays {colors, title}smoother/span are data arguments, not style, and the editor's re-journaled call preserves them untouched (asserted by test).
  • Weights are not applied — scatter's recorded semantics is that weights do not move points, and the smoother summarizes the plotted points; a declared survey design is not covered (the standard ignores-design note; the item-1 list is unchanged).

Out of scope for v1 (recorded cuts)

  • No robustness iterations (lowess's iter) and no degree option (degree 1 only) — both would change every number; revisit only if asked.
  • No confidence/error band around the curve.
  • Scatter only: no smoother on the scatterplot matrix or the bubble chart in v1.
  • The gam tool (stage b of the bucket entry) stays in the bucket — the release plan's own cut, recorded up front.

Permutation tests — the v1 conventions

Decided 2026-08-18, before any code, as sub-task 7a of docs/releases/release_0.19.0.md item 7 (kept private). The permutation (randomization) test answers "how unusual is this statistic under random relabeling?" without distributional assumptions — the methods-course and small-sample companion to the t test and the correlation. Every convention below was verified against R coin 1.4-5 and perm 1.0.0.4 (installed into the item-1 R env: coin from conda-forge, perm from CRAN) by reproducing their p-values from hand enumeration in R — not recorded from memory. The pinned values live in tool/validation/generate_permutation_references.R and test/validation/permutation_reference_values.json (sub-task 7f).

The shape

A permutation boolean (plus replicates and seed) on the two existing tools where the convention is settled — independent_ttest (the two-sample mean difference) and correlate (the bivariate coefficient) — catalog +0. Engine-side the options ride the existing ttest_independent and correlate kinds, which append a "Permutation Test" table after their standard output; the asymptotic tables are unchanged. For correlate, permutation requires exactly two variables and method pearson or spearman — a wider matrix or Kendall errors as data (the matrix stays asymptotic; recorded cuts below).

The test

  1. The reported statistic is the mean difference (group 1 − group 2, the t-table's direction) or the correlation coefficient. Under permutation with fixed group sizes these are monotone-equivalent to the t statistic and to Σxy respectively, which is why the counts agree with coin's standardized linear statistic — verified, not assumed.
  2. Two-sided p counts |T_b| ≥ |T_obs|; with the tools' tails = one the count is one-sided in the observed direction (the existing t_significance convention). Verified exactly on a 4 + 4 example: two-sided 6/70 ≡ coin oneway_test(distribution = "exact")perm::permTS(method = "exact.ce"); one-sided 3/70 ≡ coin alternative = "greater".
  3. Exact enumeration is automatic when the number of distinct arrangements is at most 1,000,000C(n, n₁) for the two-sample test (balanced designs up to n = 22), n! for correlation (n ≤ 9). The exact p is b/C with the identity arrangement included (so p ≥ 1/C, never 0), no +1 correction — it is the whole distribution — and no CI (the p is exact). Verified: Pearson |r| enumeration 12/5040 on an n = 7 pair (an independent hand implementation — coin's exact algorithms refuse the correlation problem, "not a two-sample problem"); the Spearman |ρ| enumeration coincides exactly with cor.test(method = "spearman")'s exact p (62/5040) on tie-free data.
  4. Monte Carlo above the threshold: default m = 10,000 resamples (coin's nresample default), clamped at 1,000,000 (the bootstrap cap, with the same clamp note), and p̂ = (b + 1)/(m + 1) — the Davison & Hinkley / Phipson & Smyth estimator, and exactly perm's exact.mc (verified: b = 109, m = 999 → 0.11). Recorded divergence: coin reports the raw b/m (verified: p·m integer, p·(m+1) not); the +1 form is used because a Monte-Carlo p of exactly zero misstates the evidence.
  5. The MC confidence interval is 99% Clopper–Pearson on b of m — exactly coin's conf.int (verified equal to binom.test(b, m, conf.level = 0.99) to machine precision). It describes simulation uncertainty about the permutation p, not sampling uncertainty about the effect — the note says so.
  6. Reproducibility: replicates ride the shared SplitMix64 generator (resample.rs's contract; seed option, default 1). A two-sample replicate draws a random n₁-subset by partial Fisher–Yates (n₁ below draws); a correlation replicate is a full Fisher–Yates shuffle of y (n − 1 draws) — fixed draw counts in fixed replicate order, so a pinned seed reproduces the table byte-for-byte on every platform. Cancellation is polled at replicate boundaries (the resampling loops are what the §4 channel is for), exactly like the bootstrap.
  7. Tie counting under floating point: an arrangement counts as extreme when |T_b| ≥ |T_obs| − 10⁻¹² · max(1, |T_obs|) — the regenerated observed arrangement must count itself even when summation order perturbs the last bits.
  8. Spearman permutes the midranks: ranks are computed once on the observed y and permuted (a permutation of y permutes its ranks). Ties are allowed — the exact p is then over midrank arrangements, where cor.test would fall back to asymptotics (recorded; the pins use tie-free data so the cross-check stays exact).

Errors as data

  • Weight Cases + permutation errors as data: rearranging weighted cases compounds two case multiplicities (the propensity / survey-design precedent); the message says to remove the weight.
  • permutation with one-sample or paired t tests, with Kendall, or with a correlation matrix wider than two variables: errors as data, each naming the unsupported combination.
  • A degenerate correlation pair (fewer than two complete pairs, or no variation on either variable) errors as data — there is no coefficient to permute. A t-test variable with an empty group keeps its row with null cells instead (the asymptotic table's own NaN-row convention), since the other variables in the same run may be fine.

The table

"Permutation Test" — one row per test variable (t test) or the single pair (correlate): the statistic, a Method cell naming what actually ran ("Exact, C arrangements" or "Monte Carlo, m resamples"), the permutation Sig. under the tails label, and Lower/Upper of the 99% CI of Sig. (Monte Carlo rows only; empty for exact). The note names the estimator ((b+1)/(m+1)), the seed when any row is Monte Carlo, and what the CI means. Exact-vs-MC is decided per row (missing data can put one variable under the threshold and another over it).

Out of scope for v1 (recorded cuts)

  • No sign-flip (paired/one-sample) permutation test — a different convention (2ⁿ sign assignments), recorded for a later item if asked.
  • No method-force option: exact-vs-Monte-Carlo is automatic at the recorded threshold and named in the output.
  • No Kendall permutation and no permutation p's for a correlation matrix — the matrix stays asymptotic.

VAR / Granger causality — the v1 conventions

Decided 2026-08-19, before any code, as sub-task 8a of docs/releases/release_0.19.0.md item 8 (kept private). The vector autoregression is the standard multi-series companion to the shipped single-series tools: K series regressed on p lags of all of them, the lag-order criteria that choose p, and the pairwise Granger tests that ask "do the lags of one series improve the prediction of another?". Every convention below was verified against R vars 1.6-1 and lmtest 0.9-40 (vars installed into the item-1 R env from conda-forge — the CRAN source build needs a C toolchain the env does not carry) — the VARselect scalings read from its source and every estimation identity reproduced numerically to machine precision, not recorded from memory. The pinned values live in tool/validation/generate_var_references.R and test/validation/var_reference_values.json (sub-task 8g).

The shape

A new var_model tool (catalog +1): K ≥ 2 numeric series taken in case order, fitted as a VAR(p) by per-equation OLS, with the lag order either fixed (lags, default 1) or selected by an information criterion (auto, over 1..max_lags, default max 10 — VARselect's lag.max default) under the ic option (aic default, or hq, bic, fpe). Deterministic terms ride a trend option with the stationarity tool's vocabulary extended to vars' four: constant (default), trend, both, none.

  • Alignment is listwise over the K series: the analyzed sample is the included rows where all K series are non-missing, in case order (the CCF's aligned_pair convention widened to K), with the aligned count named in the output.
  • The trend regressor is the 1-based time index of the aligned series, running p+1..n over the fitted observations — exactly vars::VAR / VARselect (seq(lag.max + 1, length = sample); verified by reproducing the trend-model coefficients).
  • Deliberately unweighted — the standing time-series convention (a frequency weight has no defined meaning for an ordered series): the weight only selects cases, surfaced by the usual unweighted note. Not design-covered (the generic design-ignored note applies).

Estimation (verified identities)

  1. Per-equation OLS on the common regressor matrix [y₁,ₜ₋₁ … y_K,ₜ₋₁, y₁,ₜ₋₂ …, deterministic] — lag-major, variable-minor ordering (x.l1, y.l1, x.l2, y.l2, …, R's embed layout), deterministic terms last, over the T = n − p effective observations. Coefficients and standard errors are the equationwise OLS ones — s²ⱼ (X'X)⁻¹ with s²ⱼ = RSSⱼ/(T − k) and k = pK + d regressors per equation (d deterministic columns). Verified equal to summary(VAR(...)) to machine precision, and t tests use the Student t on T − k df.
  2. The residual covariance matrix is the df-adjusted, column-centered one: Σ̂ᵢⱼ = Σₜ(eᵢₜ − ēᵢ)(eⱼₜ − ēⱼ)/(T − k) — read from summary.varest's source (cov(resids)·(obs−1)/(obs−k)) and verified on the type = "none" case, where the centering genuinely differs (no constant absorbs the residual means; with a constant the two forms coincide). The correlation matrix is its cov2cor. The ML covariance Σ̃ = E'E/T (raw cross-products, denominator T) is used only inside the criteria and the log-likelihood.
  3. Log-likelihood is the Gaussian one at the ML estimate: −(TK/2)·ln 2π − (T/2)·ln det Σ̃ − TK/2 — verified equal to logLik(varest).

Lag-order criteria (VARselect's scalings, from source)

With T the estimation sample and d = detint deterministic columns:

  • AIC(p) = ln det Σ̃ₚ + (2/T)·(pK² + Kd)
  • HQ(p) = ln det Σ̃ₚ + (2 ln ln T/T)·(pK² + Kd)
  • BIC(p) = ln det Σ̃ₚ + (ln T/T)·(pK² + Kd) — Lütkepohl's SC; labeled BIC in the output (the app's name for it), with the note naming the equivalence.
  • FPE(p) = ((T + k)/(T − k))^K · det Σ̃ₚ, k = pK + d.

Two sample conventions, both recorded:

  • The selection table holds the sample fixed at T = n − max_lags for every candidate order (each candidate fits on the same observations — VARselect's loop, read from source), so the criteria are comparable; the selected order per criterion is the first minimum (which.min).
  • The fitted model's own criteria (model-summary table) use the same formulas on the model's own sample T = n − p — equal to VARselect(lag.max = p)'s column p, and deliberately not equal to the selection table's row p when max_lags > p (different sample; the note says so).

Granger tests (the pairwise convention)

For every ordered pair cause → effect (K(K−1) rows): the F test that all p lag coefficients of the cause are zero in the effect's equation of the full fitted system — the per-equation Wald W = (Rβ)'(R V̂ R')⁻¹(Rβ) with the equation's OLS covariance, reported as F = W/p on (p, T − k) df. Verified: in a bivariate system this is exactly lmtest::grangertest's restricted-vs-unrestricted F (equal to machine precision, same df) — which is why the pins use grangertest for the pairwise case. Recorded divergence: vars::causality tests the joint hypothesis cause → all-others with a system-wide covariance — a different (blockwise) test, out of scope.

Errors as data

  • Fewer than two series, a non-numeric series, lags < 1, or max_lags < 1 — each named plainly.
  • Too few observations: the fit needs T − k ≥ 1 (T > pK + d), and under auto the fixed selection sample needs n − max_lags > max_lags·K + d; the message names the counts.
  • Collinear regressors (a constant series, duplicated series, n too small for the lag matrix) surface as the shared singular-matrix error.

The tables

  1. "VAR Model Summary" — Statistic/Value rows (the ARIMA Model Fit shape): aligned observations, effective observations T, lag order p, parameters per equation k, log-likelihood, AIC, HQ, BIC, FPE (own sample). The note names the series, the deterministic terms, and how p was chosen (fixed, or selected by which criterion over what range).
  2. "Lag Order Selection" (auto only) — one row per candidate order: Lag, AIC, HQ, BIC, FPE, on the fixed n − max_lags sample; the note names the per-criterion selections and the order actually fitted.
  3. "VAR Coefficients: 〈series〉" per equation — Term, Estimate, Std. Error, t, Sig.; terms in the verified x.l1, y.l1, …, const, trend order; the note carries the equation's residual SE and R².
  4. "Residual Covariance Matrix" and "Residual Correlation Matrix" — K×K, df-adjusted denominator named in the note.
  5. "Granger Causality Tests" — Cause, Effect, F, df1, df2, Sig.; the note states the hypothesis form and the bivariate grangertest equivalence.

Out of scope for v1 (recorded cuts)

  • No IRFs and no FEVD — deferred until asked (the bucket pointer recorded when the queue item moved into the release).
  • No forecasting from the fitted VAR — the single-series tools forecast; a VAR predict is a later ask.
  • No seasonal dummies and no exogenous regressors (VARselect's season/exogen arguments).
  • No joint or instantaneous causality (vars::causality's two tests, including its bootstrap option) — the pairwise matrix is the v1 answer.
  • No stability (companion-root) table — the companion matrix is not symmetric, and the engine carries no general eigensolver; recorded rather than approximated.
  • No SVAR, no VECM/Johansen, no coefficient restrictions.