3.1 Forecast representation

Every score in corrscore takes a forecast as a plain dict tagged by a "kind" key, rather than a class hierarchy, so that a caller’s own forecasting code needs no dependency on the package beyond producing this shape:

{"kind": "point", "Q": Q}
{"kind": "mixture", "components": [(p_1, Q_1), ..., (p_J, Q_J)]}
{"kind": "isotropic_gaussian_mixture", "components": [(p_1, Q_1, sigma_1), ...]}
{"kind": "ensemble", "draws": [Q_1, ..., Q_M]}

These four kinds span a genuine closed-form tractability spectrum, not an arbitrary API choice. A "point" forecast is a single deterministic matrix ($M = 1$). A "mixture" is any number of discrete atoms with associated probabilities – exactly determined, so every expectation in Equation (1) reduces to a finite weighted sum, at $O(J^2)$ cost, where $J$ is the number of atoms, with no simulation. An "isotropic_gaussian_mixture" adds an isotropic Gaussian scatter around each atom in the $K(K-1)/2$-dimensional space of free (upper-triangle) matrix entries; for the energy score this remains exactly computable, via the mean-norm formula for a noncentral chi distribution (Johnson et al., 1994), evaluated through the confluent hypergeometric function ${}_1F_1$ (scipy.special.hyp1f1) rather than simulated. The convention is as follows. A draw is a symmetric, zero-diagonal perturbation of an atom: a vector $u \sim \mathcal{N}(\mu, \sigma^2 I_n)$ with $n = K(K-1)/2$, placed into both triangles. The Frobenius norm of such a matrix counts each free entry twice, so $\lVert A \rVert_F = \sqrt{2}\,\lVert u \rVert_2$. The noncentral-chi mean is therefore taken with $n$ degrees of freedom and scale $\sigma$, and multiplied by $\sqrt{2}$. Because the scatter lives in the free-entry space, a draw need not be positive semidefinite or lie in $[-1, 1]$. The forecast is thus a distribution over symmetric matrices that may leave the set of valid correlation matrices, and the score remains well defined because the realized outcome $y$ is valid. An "ensemble" is a general Monte Carlo draw set with no assumed structure – the only kind with no closed form, falling back to the direct $O(M^2)$ pairwise sum in Equation (1).

3.2 The energy score

matrix_energy_score(forecast, y) implements Equation (1) with $\psi$ the Frobenius distance between two $K \times K$ matrices, dispatched across all four forecast kinds of Section 3.1. For a "point" forecast this reduces to plain Frobenius distance to the realized outcome $y$; for a "mixture" forecast with components $(p_j, Q_j)$, $j = 1, \dots, J$, $$\begin{equation} S_{\text{ES}}(F, y) \;=\; \sum_j p_j \lVert Q_j - y \rVert_F \;-\; \tfrac{1}{2} \sum_{j,l} p_j p_l \lVert Q_j - Q_l \rVert_F. \label{eq:esmixture} \end{equation}$$

3.3 The variogram score

matrix_variogram_score(forecast, y, p=0.5, weights=None,
                        n_samples=500, random_state=None)

implements the variogram score (Scheuerer & Hamill, 2015), $$\begin{equation} S_{\text{VS}}(F, y) \;=\; \sum_{i < j} w_{ij} \Big( |y_i - y_j|^p - \mathsf{E}_F|X_i - X_j|^p \Big)^2, \label{eq:vs} \end{equation} \qquad(2)$$ adapted to index $i, j$ over the $K(K-1)/2$ free upper-triangle entries of the correlation matrix, rather than over the original $K$-vector the published formula was defined for: the diagonal (always exactly $1$ for a correlation matrix) and the mirrored lower triangle both carry no information and are excluded. p defaults to $0.5$, matching Scheuerer & Hamill (2015)’s own recommended default; weights defaults to uniform. The score is exact (no sampling) for "point" and "mixture" forecasts, since both are fully determined by a finite set of deterministic atoms. For "ensemble" and "isotropic_gaussian_mixture" forecasts, $\mathsf{E}_F|X_i - X_j|^p$ has no known closed form for general $p$ and is instead estimated by Monte Carlo with n_samples draws (default $500$); cost and memory scale as $O(n_{\text{samples}} \times K^2(K-1)^2/4)$, which the default keeps modest through $K = 16$.

The Monte Carlo path is slightly biased upward. Equation (2) squares a residual that contains the estimated $\mathsf{E}_F|X_i - X_j|^p$, and squaring a noisy estimate adds its variance, so the expected Monte Carlo score exceeds the exact one by $\sum_{i<j} w_{ij}\,\mathrm{Var}_F(|X_i - X_j|^p)/M$, where $M$ is the number of draws (n_samples, or the ensemble size for an "ensemble" forecast). This is always nonnegative and shrinks as $1/M$. For an outcome drawn from the forecast itself, the exact expected score equals $\sum_{i<j} w_{ij}\,\mathrm{Var}_F(|X_i - X_j|^p)$, so the relative inflation is $1/M$, or $0.2\%$ at the default $M = 500$. A simulation with $K = 4$ and an isotropic Gaussian mixture matched this to within Monte Carlo error at $M = 50$ to $2000$. The bias is negligible against typical score gaps, but it tilts a backtest that compares a Monte Carlo–scored forecast against a closed-form one ("point" or "mixture") against the former. When the forecasts being compared differ by less than about $1/M$ in relative score, give every candidate the same kind or raise $M$.

3.4 A geometric variogram score

A correlation matrix is constrained to a curved region (the elliptope), not an ordinary Euclidean cube, and that constraint becomes increasingly binding as any entry approaches $\pm 1$. Pinson & Tastu (2013) show the energy score is close to blind to this kind of boundary-proximate danger for Gaussian forecasts; Nguyen (2026) shows the flat variogram score of Section 3.3 inherits a related blind spot – its concave power transform saturates once two entries are already far apart, which is exactly what happens near the boundary. A second, unrelated blind spot follows directly from Equation (2)’s reindexing over free entries: because the score is built entirely from differences between entries, it is exactly invariant to a uniform shift in the overall level of correlation – a forecast that is systematically biased high or low by the same amount across every pair scores identically to one that is not, since only the relative structure across pairs, never the shared level, is compared. Nguyen (2026) derives a fix for the saturation blind spot: matrix_geodesic_variogram_score runs the identical machinery of Equation (2), but with every free entry $\rho$ first passed through $$\begin{equation} \phi(\rho) \;=\; \operatorname{sign}(\rho) \sqrt{\tfrac{1}{2} \Big[ \log(1+|\rho|)^2 + \log(1-|\rho|)^2 \Big]}, \label{eq:phi} \end{equation}$$ the signed affine-invariant (Fisher–Rao) geodesic distance from independence, treating each entry as its own isolated $2 \times 2$ correlation matrix. $\phi$ is smooth, odd, and strictly increasing on $(-1, 1) \to \mathbb{R}$, and stretches the scale specifically near $\pm 1$ where the untransformed score has already saturated – a close cousin of the century-old Fisher $z$-transform $\operatorname{arctanh}(\rho)$. Both diverge logarithmically as $\rho \to \pm 1$, at the same order: $\phi(\rho)/\operatorname{arctanh}(\rho) \to \sqrt{2}$. Because $\phi$ is nonlinear, matrix_geodesic_variogram_score does not leave the level-shift blind spot above fully intact either, but the repair is asymmetric: it is strongest for a shift that pushes correlations further from independence, where $\phi$’s slope is steepest, and weakest for one that pulls them toward independence, where $\phi$ is nearly linear and the flat score’s exact invariance is closely approximated. Because Equation (2)’s propriety argument holds for any fixed measurable per-entry transform, not only ones with a metric interpretation (Nguyen, 2026), matrix_geodesic_variogram_score needs no separate propriety proof: it is matrix_variogram_score itself, unmodified, called on $\phi$-transformed entries, and shares that function’s exact/Monte-Carlo split by forecast kind and its p/weights/ n_samples semantics. Nguyen (2026) validates this construction on real, walk-forward regime-switching correlation forecasts. The geometric variant’s advantage over the flat variogram score is real but conditional: it is concentrated in market episodes where the realized correlation genuinely approaches the boundary, and is not present uniformly across an unconditional average. We report this precisely rather than claim an unconditional improvement. For the full propriety proof, the open questions it leaves, and the complete empirical picture, we refer the reader to that paper. One such open question is that a companion geometric energy score does not, to date, demonstrate a comparable real-data advantage.

3.5 Choosing between the scores

All three scores are proper, so none rewards a dishonest forecast, but they penalize different errors, and the choice should follow what the evaluation is meant to catch. The energy score compares whole matrices through the Frobenius distance, so it is sensitive to the overall level of correlation: a forecast biased uniformly high or low is penalized. It is the natural choice when the level itself matters, as in a risk forecast that should not be systematically too optimistic, and it is the only one of the three with an exact form for an isotropic Gaussian mixture (Section 3.1). The flat variogram score compares differences between entries, so it is built to judge the relative structure across pairs and is exactly blind to a uniform shift in level (Section 3.4). It suits a question about which pairs move together more than others, and it should not be used alone where a common bias would matter. The geometric variogram score keeps that focus on relative structure and corrects the flat score’s saturation near the boundary, so it is the better candidate when forecasts are compared during episodes where correlations approach $\pm 1$. Its advantage there is conditional (Nguyen, 2026), and it only partly repairs the level-shift blind spot. Because the scores answer different questions, a backtest that reports the energy score alongside one variogram score is more informative than either alone, and a disagreement between them is itself a finding about where a forecast is wrong. Where a forecast has to be scored by sampling, the variogram scores’ small upward bias (Section 3.3) is a further reason to compare candidates of the same forecast kind.

4.1 The zero-overlap backtest harness

backtest_zero_overlap runs a walk-forward evaluation across a sequence of forecast origins and scores one or more named forecasting methods against a shared, model-independent ground truth:

backtest_zero_overlap(forecast_fns, ground_truth_fn, origins, horizon,
                       purge_gap=0, score_fn=matrix_energy_score,
                       severity_fn=None)

For each origin in origins, the harness computes ground truth from a window starting strictly after the origin:

y = ground_truth_fn(origin + 1 + purge_gap,
                     origin + 1 + purge_gap + horizon)

then scores every forecast_fns[name](origin) against that same y. The specific design choice worth calling out is what the harness cannot and does not attempt to police: how much history a caller’s own forecast_fn consults internally is opaque to it (a full-history discounted filter and a short trailing window are both permitted), exactly the same responsibility boundary scikit-learn’s TimeSeriesSplit leaves to its own caller. What the harness does enforce, unconditionally and by construction, is that ground_truth_fn is only ever called with a start point strictly after the origin. This specific guarantee exists because of a real bug. An earlier walk-forward evaluation, in the companion regime-switching correlation study this package grew out of, computed “realized” ground truth from a window ending near the forecast origin, not one starting strictly after it. That silently leaked information back into the evaluation. It made a naive persistence forecast look decisively better than a genuine forecasting model, the opposite of the study’s eventual, correctly-measured finding. backtest_zero_overlap’s API makes that specific mistake structurally impossible to reproduce: the caller never controls, and never sees, the start point passed to ground_truth_fn.

One cost of this design deserves stating. The realized outcome that ground_truth_fn returns is itself a windowed sample estimate of the correlation, not the true correlation, so it carries its own sampling error, which is larger the shorter the horizon. Scoring-rule theory treats the outcome as a draw from the true law, whereas here that estimation error is added to the forecast’s own error and enters every reported score. It also adds noise to the per-origin score differences, so a short horizon calls for more origins before two forecasts can be separated. The horizon=10 window of $K = 4$ series in Section 5 is short in exactly this sense.

4.2 Significance testing

Three complementary significance tests operate on the per-origin score arrays a BacktestResult returns.

circular_block_bootstrap(scores_a, scores_b, block_lengths,
                          n_boot=2000, seed=None)

runs a percentile circular block bootstrap (Politis & Romano, 1992) on the paired differential $d_i = a_i - b_i$, where $a_i, b_i$ are the $i$th entries of scores_a, scores_b, via arch’s CircularBlockBootstrap (Sheppard, 2024), swept across a list of candidate block lengths as a sensitivity check rather than relying on one automatically “optimal” length. The range of $p$-values and intervals across the swept lengths is the result to report. If a single headline entry is wanted, the largest-$p$ one is a conservative choice, but it is a choice made after looking at the sweep, and a block length picked by a data-driven rule (for example arch’s optimal_block_length) is the cleaner alternative to state in advance.

diebold_mariano(loss_a, loss_b, h=1, varestimator="acf")

implements the Diebold–Mariano test (Diebold & Mariano, 1995) on the same paired differential, including the Harvey–Leybourne–Newbold small-sample correction (Harvey et al., 1997) and a Student-$t$ (not asymptotic normal) reference distribution with $n-1$ degrees of freedom, matching R’s forecast::dm.test (Hyndman & Khandakar, 2008) formula for formula rather than following the original 1995 paper’s asymptotic-normal reference directly. The h argument sets how many autocovariance lags enter the long-run variance ($h - 1$), and it also drives the correction factor. It should reflect the serial dependence in the per-origin score differences, not be set mechanically to the forecast horizon. When origins are spaced so that evaluation windows are adjacent and do not overlap, the differences are close to serially independent and h near $1$ is appropriate, while a larger h widens the variance estimate and costs power. This vendored implementation is deliberately cross-checked against forecast::dm.test’s own output on fixed synthetic data as a development-time oracle: doing so caught a genuine normalization bug during development, where the package’s long-run-variance estimator initially normalized each sample autocovariance by $n - \text{lag}$ (the textbook “unbiased” convention) rather than by the full sample size $n$ at every lag, the convention R’s own acf() function – and hence forecast::dm.test – actually uses. The two conventions coincide at lag $0$, so this discrepancy was invisible at horizon $h=1$ and only surfaced once the test suite exercised $h > 1$ against the oracle. The oracle is stored with the package as fixed data (tests/_reference/): six cases at $h \in \{1, 5, 10\}$ with both variance estimators, giving R’s own random draws together with the statistic and $p$-value forecast::dm.test returned, which the test suite matches to $10^{-6}$. The comparison therefore runs from Python alone, with no R installation. With R 4.6.0 and forecast 9.0.2, regenerating the fixture reproduces it byte for byte (make check-oracles).

model_confidence_set(scores, alpha=0.10, block_len=5,
                      n_boot=1000, seed=None)

implements the Model Confidence Set (Hansen et al., 2011) via its range-statistic elimination algorithm: repeatedly test whether the current candidate set of models is statistically distinguishable from its own best member, and if so, drop the single worst-performing model and repeat, until the surviving set cannot be rejected at level alpha. The null distribution and each round’s studentizing standard errors both come from the same joint circular block bootstrap of the full loss matrix, reusing the identical arch machinery Section 4.2’s bootstrap function depends on. This is the package’s own reasonable operationalization of Hansen et al.’s procedure, not a line-by-line port of any specific existing implementation’s internal choices – R’s MCS package’s automatic block-length selection, for instance, is not reproduced, and block_len is instead a required, caller-chosen parameter, swept manually if robustness to it matters. The test suite accordingly cross-checks this implementation’s verdict (which models survive) against R’s MCS::MCSprocedure on fixed synthetic data where the correct verdict is unambiguous by construction, rather than requiring byte-exact statistic or $p$-value agreement. The fixed loss matrix, drawn in R, and R’s full output for it (MCS 0.2.0) are stored alongside the Diebold–Mariano fixture. On that matrix both implementations keep the same two models and drop the same two, and the final Model Confidence Set $p$-value is about $0.89$ here against $0.90$ in R, although the two eliminate the dropped models in a different order.

5.1 Two point forecasts, end to end

The following example simulates a $K = 4$ asset panel whose true pairwise correlation drifts slowly between roughly $0.15$ and $0.50$ over 900 trading days, and compares two deliberately simple forecasting rules: persistence, which forecasts tomorrow’s correlation matrix as the sample correlation of the trailing 60-day window, and shrinkage, which additionally shrinks that estimate 30% of the way toward the identity matrix. Because the underlying process has real, substantial, slowly-varying correlation, shrinking toward independence is the wrong thing to do here by construction – the example is deliberately built so that one forecast is genuinely worse, to show what a decisive verdict from the full toolchain looks like end to end.

The output below was produced with corrscore 0.1.2 on Python 3.12.7 with NumPy 1.26.4, SciPy 1.13.1 and arch 7.0.0. Every random draw is seeded (the data with default_rng(0), the bootstrap with seed=1, the Model Confidence Set with seed=2) and every forecast is a "point" forecast, so the output is deterministic for fixed library versions on a fixed platform. It does depend on both. NumPy 2.x draws different returns from the same seed, and the same script then prints a mean variogram score of $3.0639$ for persistence, not $3.0274$. The numbers below were produced on macOS (Apple silicon, arm64). On Linux (aarch64) with the same pinned versions, 13 of the 29 printed lines differ: the mean variogram score of persistence is $3.0974$, and in the near-tie of Section 5.3 the Model Confidence Set drops win35 at all three block lengths, while the verdict of this first example is unchanged. The cause is the draw: Generator.multivariate_normal factorises the correlation matrix by a singular value decomposition, the matrix here has an eigenvalue of multiplicity three, and the macOS and Linux builds of NumPy 1.26.4 return different bases for it on 16 of the 900 simulated days. Passing method="cholesky" makes both platforms print the same lines, but the numbers reported here are those of the default method. The script is examples/section5_example.py in the repository, and running make paper-example from the repository root builds the pinned environment of requirements-paper.txt and prints the numbers below, followed by those of Sections 5.2 and 5.3. Release 0.1.2 is archived at https://doi.org/10.5281/zenodo.23078623, so these results attach to a fixed snapshot.

import numpy as np
from corrscore import (backtest_zero_overlap, circular_block_bootstrap,
                        diebold_mariano, matrix_geodesic_variogram_score,
                        matrix_variogram_score, model_confidence_set)

rng = np.random.default_rng(0)
K, n_days = 4, 900
true_rho = 0.15 + 0.35 * (0.5 + 0.5 * np.sin(np.arange(n_days) / 60.0))
returns = np.empty((n_days, K))
for t in range(n_days):
    corr = np.full((K, K), true_rho[t]); np.fill_diagonal(corr, 1.0)
    returns[t] = rng.multivariate_normal(np.zeros(K), corr)

def sample_corr(start, end):
    return np.corrcoef(returns[start:end], rowvar=False)

def persistence_forecast(origin, window=60):
    return {"kind": "point", "Q": sample_corr(origin - window + 1, origin + 1)}

def shrinkage_forecast(origin, window=60, alpha=0.3):
    q = sample_corr(origin - window + 1, origin + 1)
    return {"kind": "point", "Q": (1 - alpha) * q + alpha * np.eye(K)}

origins, horizon = list(range(120, n_days - 10, 10)), 10
result = backtest_zero_overlap(
    forecast_fns={"persistence": persistence_forecast, "shrinkage": shrinkage_forecast},
    ground_truth_fn=sample_corr, origins=origins, horizon=horizon,
    score_fn=matrix_variogram_score)

print("mean VS, persistence:", result.scores["persistence"].mean())
print("mean VS, shrinkage:  ", result.scores["shrinkage"].mean())

boot = circular_block_bootstrap(result.scores["shrinkage"], result.scores["persistence"],
                                 block_lengths=[5, 10], seed=1)
def fmt_boot_p(r):
    # A percentile bootstrap with n_boot resamples cannot resolve a p-value
    # below 2/n_boot (two tails of one resample each); report an upper bound.
    floor = 2.0 / r.n_boot
    return f"p<{floor:g}" if r.p_value < floor else f"p={r.p_value:.4f}"

for block_len, r in boot.items():
    print(f"bootstrap block_len={block_len}: mean diff={r.obs:.4f}, "
          f"95% CI=({r.ci_lo:.4f}, {r.ci_hi:.4f}), {fmt_boot_p(r)}")

dm = diebold_mariano(result.scores["shrinkage"], result.scores["persistence"], h=horizon)
print(f"Diebold-Mariano: stat={dm.statistic:.3f}, p={dm.p_value:.1e}")

mcs = model_confidence_set(result.scores, alpha=0.10, block_len=10, seed=2)
print("MCS survivors:", mcs.survivors, " eliminated:", mcs.eliminated)

mean VS, persistence: 3.0274186774073626
mean VS, shrinkage:   3.52280921068856
bootstrap block_len=5: mean diff=0.4954, 95% CI=(0.4302, 0.5600), p<0.001
bootstrap block_len=10: mean diff=0.4954, 95% CI=(0.4267, 0.5592), p<0.001
Diebold-Mariano: stat=13.026, p=4.7e-21
MCS survivors: ['persistence']  eliminated: [('shrinkage', 0.0)]

All three significance tools agree that shrinkage scores worse than persistence on this panel, and the Model Confidence Set eliminates it at $\alpha = 0.10$. The bootstrap $p$-values are upper bounds. With n_boot=2000 the smallest nonzero two-sided value is $2/2000 = 0.001$, and a smaller reported value means no resample crossed zero. The Model Confidence Set’s $p$-value of $0.0$ is limited by its own resample count in the same way (its default of $1000$ resamples resolves no finer than $1/1000$). The Diebold–Mariano $p$-value comes from a Student-$t$ reference distribution and is not limited by a resample count.

Three caveats apply to how this example should be read. First, it demonstrates the toolchain end to end and is not evidence about the geometric score. Every off-diagonal entry shares one true level and shrinkage rescales that level, which is close to the level-shift configuration of Section 3.4, and the panel does not exercise the boundary-proximate case where that section says the geometric variant helps. Second, the bootstrap and Diebold–Mariano tests assume a covariance-stationary loss differential, and a correlation that drifts between $0.15$ and $0.50$ makes the per-origin differences plausibly nonstationary in mean. The printed $p$-values therefore demonstrate the toolchain and are not a validated inference. Third, h=10 equals the horizon and origins are spaced $10$ apart, so adjacent evaluation windows just abut. The Harvey–Leybourne–Newbold correction and the $n-1$ degrees of freedom are applied as in forecast::dm.test, but here the $h-1 = 9$ lags estimate near-zero autocovariances (the lag-1 autocorrelation of the differential is $-0.09$), which makes $h=10$ conservative: $h=1$ gives a statistic of $14.206$ in place of $13.026$, with the same verdict.

Passing matrix_geodesic_variogram_score as the score_fn argument to backtest_zero_overlap reuses every other line unchanged, since Section 3.4’s geometric score shares the flat variogram score’s exact call signature – a deliberate design consequence of implementing it as a transform composed with existing machinery rather than a parallel code path.

5.2 Closed form against the Monte Carlo fallback

The example above uses only "point" forecasts. This one exercises the closed-form kinds of Section 3.1. It scores a two-atom forecast (a calm and a stress correlation matrix with probabilities $0.7$ and $0.3$, $K = 4$) against a fixed outcome $y$, once through its exact kind and once as an "ensemble" of $M = 2000$ draws from it. For the "isotropic_gaussian_mixture" kind each draw adds symmetric Gaussian scatter with $\sigma = 0.05$ to the free entries of its atom, following the convention of Section 3.1. The ensemble score is averaged over 50 independent seeds, and the difference from the closed form is expressed in standard errors of that average.

mixture = {"kind": "mixture", "components": [(0.7, q_calm), (0.3, q_stress)]}
iso = {"kind": "isotropic_gaussian_mixture",
       "components": [(0.7, q_calm, 0.05), (0.3, q_stress, 0.05)]}

matrix_energy_score(mixture, y)                 # exact
matrix_energy_score(ensemble_from(mixture), y)  # M = 2000 draws
matrix_energy_score(iso, y)                     # exact (noncentral chi)
matrix_energy_score(ensemble_from(iso), y)
matrix_variogram_score(mixture, y, p=0.5)       # exact
matrix_variogram_score(ensemble_from(mixture), y, p=0.5)
matrix_variogram_score(iso, y, p=0.5, n_samples=2000)  # Monte Carlo
matrix_variogram_score(ensemble_from(iso), y, p=0.5)

K=4, ensemble size M=2000, mean over 50 seeds
energy, mixture                exact=0.26990  ens=0.27003+-0.00027  (+0.5 se)
energy, isotropic Gauss. mix.  exact=0.22600  ens=0.22622+-0.00035  (+0.6 se)
variogram, mixture             exact=0.03327  ens=0.03277+-0.00036  (-1.4 se)
variogram, isotropic (both MC) sampler=0.04213+-0.00030  ens=0.04202+-0.00034  (+0.2 se)

(ensemble_from stands for the draw-and-average loop in the script examples/closed_form_vs_mc.py, which also contains the exact matrices and seeds.) Both energy scores agree with their closed forms to within one standard error, about $0.1\%$ of the score, which checks the isotropic-Gaussian-mixture formula, including its $\sqrt{2}$ Frobenius bookkeeping, against direct simulation. The variogram score of the discrete mixture agrees within $1.4$ standard errors. That row is noisy because the score itself is small, and a run with 600 seeds narrows the gap to $0.8$ standard errors. The Monte Carlo bias of Section 3.3 is a $4 \times 10^{-5}$ effect in this setting, far below the seed noise. The isotropic-Gaussian-mixture variogram score has no closed form, so the last row compares the kind’s own sampler with an ensemble of independent draws, and the two agree within $0.2$ standard errors.

5.3 A near-tie among six models

The verdict in Section 5’s first example was decided by construction. This one is not. It uses the same simulated panel and six candidates that differ only in the length of the trailing window of their persistence forecast ($20$, $25$, $30$, $35$, $40$ and $60$ days). Windows near the best length give score differences that change sign from origin to origin, so the outcome is not known in advance. The win60 model is the persistence model of Section 5, and its mean score of $3.027$ matches. Each model is compared with win30, with the bootstrap swept over block lengths $2$, $5$ and $10$ and the Diebold–Mariano test run at h=1, since the evaluation windows are adjacent and do not overlap (Section 4.2).

scores = backtest_zero_overlap(forecast_fns={f"win{w}": persistence(w) for w in windows},
                               ground_truth_fn=sample_corr, origins=origins,
                               horizon=10, score_fn=matrix_variogram_score).scores
dm = diebold_mariano(scores[m], scores["win30"], h=1)
boot = circular_block_bootstrap(scores[m], scores["win30"],
                                block_lengths=[2, 5, 10], seed=1)
mcs = model_confidence_set(scores, alpha=0.10, block_len=b, seed=2)   # b in 2, 5, 10

77 origins, horizon 10; mean variogram score:
  win20=2.496  win25=2.545  win30=2.550  win35=2.682  win40=2.815  win60=3.027

Each model against win30 (mean diff > 0 means worse than win30):
  model   diff    DM p   bootstrap p at block 2 / 5 / 10
  win20  -0.054   0.604   0.576 / 0.520 / 0.566
  win25  -0.005   0.949   0.906 / 0.906 / 0.943
  win35  +0.132   0.053   0.035 / 0.044 / 0.013
  win40  +0.265   0.003   0.001 / 0.002 / <0.001
  win60  +0.477  <0.001   <0.001 / <0.001 / <0.001

Model Confidence Set, alpha=0.10 (elimination path, with p at each step):
  block  2: survivors ['win20', 'win25', 'win30', 'win35']
            eliminated: win60 (p=0.001) -> win40 (p=0.015)
  block  5: survivors ['win20', 'win25', 'win30', 'win35']
            eliminated: win60 (p<0.001) -> win40 (p=0.008)
  block 10: survivors ['win20', 'win25', 'win30']
            eliminated: win60 (p<0.001) -> win40 (p=0.004) -> win35 (p=0.078)

Three readings follow. First, win20 has the lowest mean score ($2.496$ against $2.550$ for win30), but the gap is far from distinguishable from zero (Diebold–Mariano $p = 0.604$, bootstrap $p$ between $0.52$ and $0.58$), and win25 is a tie as well. The ordering of windows $20$, $25$ and $30$ is not supported by the data. Second, win35 is borderline: the Diebold–Mariano test gives $p = 0.053$ while the bootstrap gives $0.013$ to $0.044$ depending on block length, so whether it crosses the $5\%$ line depends on a choice the user makes. This is the situation in which Section 4.2’s advice to report the whole sweep, and not one entry, matters. Third, win40 and win60 are worse under every test. The Model Confidence Set drops win60 and then win40 at every block length. It keeps win35 at block lengths $2$ and $5$ (final $p$ of $0.17$ and $0.15$) and drops it at block length $10$ ($p = 0.078$). The set at $\alpha = 0.10$ is therefore either four or three models, always containing windows $20$, $25$ and $30$ and never $40$ or $60$, and membership of win35 is not settled by this panel.

The candidates were chosen to produce a near-tie, so this panel shows how to read the output and is not a finding about window lengths. The stationarity caveat of Section 5’s first example applies here as well.