The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
111,188 characters
Calibrated Order-Randomized Rosenblatt Tests
\begin{frontmatter}
\title{Calibrated Order-Randomized Rosenblatt Tests}
\author{Mehrdad Pournaderi}
\ead{[email removed]}
\affiliation{organization={Mofid Securities},
country={Iran}}
\begin{abstract}
We test whether a multivariate vector $\mathbf{X}$ conforms to a specified distribution $F$, a problem in
copula modelling and density forecasting. The Rosenblatt transform reduces it to a test of
uniformity, but depends on an arbitrary coordinate ordering that strongly affects power under
dependence. We study \emph{order randomization}: applying the transform under many random orderings and
merging the evidence with dependence-robust rules. Reordering conserves the total Mahalanobis signal
energy and merely redistributes it, so one ordering is a lucky or unlucky draw. In simulations we observe
significant gains in calibrated power over both the expected single random ordering and
order-invariant references. Two ingredients are essential: a two-sided
base statistic, and a re-estimating parametric bootstrap that restores level under an estimated null
and unlocks the gain. The calibrated pooled tests are robust to the departure's \emph{shape}; no
order-invariant reference we compare is: the symmetric-root test collapses on diffuse
departures, while the shape-flat $\chi^2$ test trails on concentrated ones. We apply it to a
Gaussian foreign-exchange risk model over a decade
of daily data on nine currencies, where it detects episodes such as Brexit and COVID. Though
we focus on Gaussian nulls, the procedure extends to any null whose conditional distributions
can be computed and simulated from.
\end{abstract}
\begin{keyword}
Rosenblatt transformation \sep goodness-of-fit \sep e-values \sep Simes test \sep
parametric bootstrap \sep false discovery rate \sep risk-model validation
\end{keyword}
\end{frontmatter}
\section{Introduction}\label{sec:intro}
Deciding whether a single multivariate vector $\mathbf{X}=(X_1,\dots,X_n)$ conforms to a specified joint
distribution $F$---the hypothesis $H_0:\mathbf{X}\sim F$---is a basic task in copula modeling,
density-forecast evaluation, simulation-based inference and model checking, and anomaly detection.
It is also the per-hypothesis building block of many multiple-testing problems, where one screens a
large collection of such multivariate hypotheses; our focus is the single test, and the multiplicity
layer (Section~\ref{sec:fdr}) is a standard add-on. A classical and very general device for the
single test is the \emph{Rosenblatt} \cite{rosenblatt1952} transformation, which writes $F$ as a
chain of conditionals and applies the probability integral transform (PIT) to each. Under $H_0$ the
resulting values are i.i.d.\ $\mathrm{Unif}[0,1]$, reducing the question to a test of uniformity addressable
by a univariate combination such as the Simes statistic \cite{simes1986}. This places the problem in
the goodness-of-fit / specification-testing tradition \cite{genest2009,dgt1998}; we stress, however,
that each hypothesis concerns a \emph{single} multivariate observation (the uniforms are its $n$
transformed coordinates, not a sample of repeated draws).
The catch is that to use the transform the analyst must first \emph{order} the coordinates, and that
choice is far from innocuous. Under dependence, different orderings whiten the data differently and
yield tests of materially different power: the same observation, ordered two ways, can comfortably
accept under one and decisively reject under the other (Example~\ref{ex:bivar} makes this concrete). The arbitrariness is well documented, and has
been treated as a nuisance to be removed---by fixing an order by convention, or by symmetrizing a
single statistic over orderings \cite{dgt1998,dovern2020}. We read it instead as an
\emph{opportunity}: if no ordering is canonical, committing to one is a gamble, and it is a gamble
one need not take, because a single observation can be transformed under many orderings at once and
the resulting evidence pooled.
Two questions then organize the paper, and the first is prior to the method. \emph{When should
the Rosenblatt reduction be used at all}, rather than a test that never touches an
ordering---the $\chi^2$ energy test, or a symmetrized whitening (\S\ref{sec:baselines})? Our
answer is a robustness argument made quantitative. Every order-invariant reference we examine
has a shape of departure it cannot handle: the symmetric-root test
is strong against concentrated departures and collapses against diffuse ones of the same
energy, while the $\chi^2$ energy test---whose power depends on the departure only through its
energy (Proposition~\ref{prop:chi2}) and is therefore \emph{flat} across shapes---trails the
pooled tests against concentrated departures and is competitive only against diffuse ones
(\S\ref{sec:headline}); the ordering ensemble is the
one contender that need not know the shape in advance. The geometry behind this is the energy
invariance of Lemma~\ref{lem:energy}: reordering conserves the total Mahalanobis signal energy
and only redistributes it across the whitened coordinates, so a fixed ordering's fortunes
depend entirely on how that redistribution happens to fall---committing to one ordering is a
gamble. The second question is then the practical one this paper answers in full:
\begin{quote}\itshape
Does applying the Rosenblatt transform under many random orderings, and combining the resulting
tests, beat using a single (arbitrarily chosen) ordering?
\end{quote}
Whether pooling actually helps is far from obvious, and answering it is most of this paper. Three
things must line up. \emph{First}, one must understand what reordering does to a departure from $F$:
we show (Section~\ref{sec:theory}) that it conserves the total Mahalanobis signal energy and merely
\emph{redistributes} it across the whitened coordinates, so every ordering is a different view of one
fixed signal, a single ordering is a lucky or unlucky draw, and pooling averages over the luck---by
an amount that, in every design we ran, grows with dependence. \emph{Second}, the per-coordinate test must look in both tails, because the whitening routinely
sends a departure into either tail depending on the ordering (Section~\ref{sec:misspec}). \emph{Third}, and this is the crux, the pooled test must be \emph{calibrated}: naive
pooling pays multiplicity penalties that almost exactly cancel its advantage, and under an estimated
null it is not even valid, so a re-estimating parametric bootstrap is needed---and in restoring level
it \emph{unveils} the very power gain those penalties had hidden (Section~\ref{sec:exp}).
When the three align, the payoff is real. Its most useful form is \emph{shape robustness}: no
order-invariant reference is safe in both regimes---the symmetric-root test is strong against
concentrated departures and collapses against diffuse ones of the same energy, where its power
falls by more than $0.30$, while the $\chi^2$ test, whose power is the same for every shape at
fixed energy, is well below the pooled tests against concentrated departures and competitive only
against diffuse ones---while pooling sits near the top of both. Against a single random ordering---the expected
power of an arbitrarily chosen one---the improvement is smaller but systematic: $0.03$--$0.06$ for
sparse departures, growing with dependence in every design we ran, in the direction the
Jensen-gap mechanism of Section~\ref{sec:theory} suggests (from $+0.02$ at
$\rho=0.2$ to $+0.06$ at $\rho=0.8$, known covariance), up to $+0.17$ for dense ones, and, for
sparse departures, unchanged when the covariance must be estimated. The geometry behind this is provable, and Section~\ref{sec:theory} does so: the log-evidence
advantage of order-averaging is exactly the spread of views across orderings. We identify the
combiners to use and the number of orderings, and mark where pooling is merely neutral
(near-independence). Screening many such hypotheses at once needs no new machinery, only
BH or e-BH on top (Section~\ref{sec:fdr}); it is also where shape robustness pays most, since the
hypotheses in a screen need not deviate from the null distribution in the same way. That is also
where a modern alternative meets us head-on, and \S\ref{sec:conformal} compares the procedure with
conformal novelty detection \cite{bates2023,marandon2024}: the two test different nulls---conformity
to a specified $F$ against exchangeability with a reference sample---and when the model is
correct, the one setting in which their power can be compared, the calibrated pooled test leads in
the sparse regime that screening actually presents, across reference samples up to ten
times the size of the screen. Section~\ref{sec:fx} takes the whole pipeline
to real data: validating a
yearly-refitted Gaussian risk model on nine currencies over 2016--2025, where the calibrated test
rejects at close to the nominal rate in calm years and the pooled e-values single out the Brexit, COVID, 2022, and 2025
episodes. Section~\ref{sec:prelim} first assembles the tools; readers fluent in the
Rosenblatt transform, Simes, and e-values may skim it.
\paragraph{Notation}
$\Phi$ and $\varphi$ are the standard normal CDF and density; $z_a=\Phi^{-1}(1-a)$ the upper-$a$
quantile. $S_n$ is the symmetric group (the set of all $n!$ orderings) on $\{1,\dots,n\}$. For
$\sigma\in S_n$, $P_\sigma$ is the permutation matrix with $(P_\sigma\mathbf{x})_k=x_{\sigma(k)}$. For a
vector $\mathbf p\in[0,1]^n$, $p_{(1)}\le\dots\le p_{(n)}$ are its order statistics (the entries
sorted increasingly). ``i.i.d.'' is independent and identically distributed. All distributions are
absolutely continuous unless stated.
\section{Preliminaries}\label{sec:prelim}
Pooling orderings needs only three ingredients: the transform that turns a multivariate hypothesis
into a uniformity question, its explicit Gaussian form (a Cholesky whitening), and the univariate
rules we pool with. We collect them here together with the validity facts the later arguments rely
on; the reader who knows them can move to Section~\ref{sec:method}.
\subsection{The Rosenblatt transformation}
Let $\mathbf{X}\sim F$ on $\mathbb{R}^n$ with absolutely continuous law, and fix an ordering $\sigma\in S_n$. Write
$F_k^\sigma(\cdot\mid\cdot)$ for the conditional CDF of $X_{\sigma(k)}$ given
$X_{\sigma(1)},\dots,X_{\sigma(k-1)}$ under $F$ (with $F_1^\sigma$ the marginal CDF of
$X_{\sigma(1)}$). The \emph{Rosenblatt map} $R_\sigma:\mathbb{R}^n\to[0,1]^n$ is
\begin{equation}\label{eq:rosen}
\big(R_\sigma(\mathbf{x})\big)_k \;=\; U_k \;=\; F_k^\sigma\!\big(x_{\sigma(k)}\mid x_{\sigma(1)},\dots,x_{\sigma(k-1)}\big),
\qquad k=1,\dots,n.
\end{equation}
\begin{theorem}[Rosenblatt, 1952]\label{thm:rosen}
If $\mathbf{X}\sim F$ then $\mathbf{U}=R_\sigma(\mathbf{X})$ has i.i.d.\ $\mathrm{Unif}[0,1]$ coordinates, for every fixed
$\sigma\in S_n$.
\end{theorem}
\begin{proof}[Proof sketch]
By the probability integral transform, conditionally on the prefix
$X_{\sigma(1)}$, \dots, $X_{\sigma(k-1)}$ the variable
$U_k=F_k^\sigma(X_{\sigma(k)}\mid\cdot)$ is $\mathrm{Unif}[0,1]$ and hence independent of the prefix.
Therefore $U_k$ is independent of $U_1,\dots,U_{k-1}$ (which are functions of the prefix) and
uniform; induction over $k$ gives joint independence and uniformity.
\end{proof}
In words: take the coordinates one at a time in the chosen order; for each, ask ``given everything
seen so far, at which quantile of its predicted conditional distribution did this coordinate land?''.
If the model is right, those quantile positions are pure noise---independent uniforms; any structure
in them is evidence against the model. Theorem~\ref{thm:rosen} is the engine: it converts
``$\mathbf{X}\sim F$?'' into ``are the $U_k$ i.i.d.\ uniform?''. The map $R_\sigma$ is a measurable bijection
(a triangular, Knothe--Rosenblatt map) whenever the conditionals are strictly increasing, a fact we
use in Section~\ref{sec:misspec}.
\subsection{The Gaussian case: Cholesky whitening}\label{sec:gaussrosen}
Let $F=\mathcal N(\mathbf 0,\Sigma)$. Let $\Sigma_\sigma=P_\sigma\Sigma P_\sigma^\top$ be the permuted
covariance and $L_\sigma$ its lower-triangular Cholesky factor, $\Sigma_\sigma=L_\sigma L_\sigma^\top$.
Define the whitened scores
\begin{equation}\label{eq:whiten}
\mathbf{z}^\sigma \;=\; L_\sigma^{-1}\,P_\sigma\mathbf{X} \;=\; L_\sigma^{-1}\mathbf{X}_\sigma .
\end{equation}
Because $L_\sigma$ is lower triangular, the $k$-th entry $z^\sigma_k$ is an affine function of
$x_{\sigma(1)},\dots,x_{\sigma(k)}$, and a direct computation identifies it with the standardized
Gaussian conditional:
\[
z^\sigma_k=\frac{x_{\sigma(k)}-\mu_{k\mid<k}(\mathbf{x})}{s_{k\mid<k}},\qquad
F_k^\sigma\!\big(x_{\sigma(k)}\mid x_{\sigma(1)},\dots,x_{\sigma(k-1)}\big)=\Phi(z^\sigma_k),
\]
where $\mu_{k\mid<k}$ and $s_{k\mid<k}$ are the conditional mean and standard deviation of
$X_{\sigma(k)}$ given the prefix. Hence the Gaussian Rosenblatt map is
$U^\sigma_k=\Phi(z^\sigma_k)$, and under $H_0$, $\mathbf{z}^\sigma\sim\mathcal N(\mathbf 0,I)$. We work
directly with the scores $\mathbf{z}^\sigma$ (equivalently the normal scores $\Phi^{-1}(U^\sigma_k)=z^\sigma_k$).
In plain terms, \emph{whitening} means transforming correlated coordinates into standardized,
uncorrelated ones: each score $z^\sigma_k$ answers ``by how many (conditional) standard deviations
did coordinate $\sigma(k)$ miss the value its predecessors predicted for it?''---so under the model
every score should look like a fresh standard normal draw.
\begin{example}[Why the ordering matters]\label{ex:bivar}
Take $n=2$, $\Sigma=\left(\begin{smallmatrix}1&\rho\\\rho&1\end{smallmatrix}\right)$ with $\rho=0.9$, and
observe the single point $\mathbf{x}=(2,0)$. The two orderings give:
\[
\begin{aligned}
\text{order }(1,2):\quad
&\mathbf{z}=\Big(x_1,\ \tfrac{x_2-\rho x_1}{\sqrt{1-\rho^2}}\Big)=(2.00,\,-4.13),\\
&\mathbf{U}=\Phi(\mathbf{z})=(0.977,\ 1.8\times10^{-5});\\
\text{order }(2,1):\quad
&\mathbf{z}=\Big(x_2,\ \tfrac{x_1-\rho x_2}{\sqrt{1-\rho^2}}\Big)=(0.00,\,+4.59),\\
&\mathbf{U}=\Phi(\mathbf{z})=(0.500,\ 0.999998).
\end{aligned}
\]
The same observation produces \emph{entirely different} transforms---and not merely a relabeling: the
multisets $\{0.977,\,1.8\times10^{-5}\}$ and $\{0.5,\,0.999998\}$ differ, so even a
permutation-symmetric combination such as Simes depends on the ordering.
\emph{Whitening view.} Both transforms apply the same triangular factor
\[
L^{-1}=\begin{pmatrix}1&0\\[2pt]-\rho/\sqrt{1-\rho^2}&1/\sqrt{1-\rho^2}\end{pmatrix}
\]
to the \emph{ordered} data $\mathbf{X}_\sigma$---here $\Sigma$ is exchangeable, so $L_\sigma$ is unchanged and only
$\mathbf{X}_\sigma$ is permuted. Triangularity makes the slot decisive (slot one passes through; slot two is the
residual after regressing it out), so the swap exchanges the roles $X_2\mid X_1$ and $X_1\mid X_2$. (For
non-exchangeable $\Sigma$, $\Sigma_\sigma$ and hence $L_\sigma$ change with the ordering too.)
The effect is interpretable.
Under strong positive correlation, $\mathbf{x}=(2,0)$ is surprising because $X_2$ ``should'' track $X_1$.
Order $(1,2)$ charges that surprise to the conditional $X_2\mid X_1$ (observed $0$ versus expected
$\rho x_1=1.8$), producing a tiny $U_2$; order $(2,1)$ charges it to $X_1\mid X_2$ (observed $2$
versus expected $0$), producing a $U_2$ near one. The surprise has the \emph{same magnitude}---indeed
$\|\mathbf{z}\|^2=\mathbf{x}^\top\Sigma^{-1}\mathbf{x}=21.05$ under both orderings (Lemma~\ref{lem:energy})---but it surfaces in
\emph{opposite tails}. A one-sided test that flags small uniforms detects it under order $(1,2)$ yet
is blind to it under order $(2,1)$; a two-sided test catches both. This single point already exhibits
the paper's two levers: reordering reshuffles a fixed amount of signal across coordinates and tails
(the ``orbit'' of Section~\ref{sec:theory}), and the base test must be two-sided
(Section~\ref{sec:misspec}) to see it whichever way it lands.
\end{example}
\begin{remark}[Why the \emph{triangular} factor]\label{rem:tri}
The whitening matrix is not unique: any $W$ with $W\Sigma_\sigma W^\top=I$ yields i.i.d.\ scores, and two
such factors differ by an orthogonal rotation. Rosenblatt singles out the lower-triangular Cholesky
factor (unique with positive diagonal) precisely because triangularity \emph{is} the chain of
conditionals: $z^\sigma_k$ standardizes $X_{\sigma(k)}$ given its predecessors. This is what makes
ordering matter---Cholesky does not commute with permutation, $\mathrm{chol}(P_\sigma\Sigma P_\sigma^\top)\neq
P_\sigma\,\mathrm{chol}(\Sigma)$. A symmetric root gives instead the whitening matrix
$W_\sigma=\Sigma^{-1/2}P_\sigma^\top$ (its inverse $P_\sigma\Sigma^{1/2}$ is the corresponding
\emph{non-triangular} covariance factor of $\Sigma_\sigma$), whose scores
$W_\sigma\mathbf{X}_\sigma=\Sigma^{-1/2}\mathbf{X}$ are \emph{order-invariant}---a
legitimate test, but a global mix rather than a chain of conditionals; for general
non-Gaussian $F$ a linear covariance-root transformation does not produce independent uniform
coordinates, and the nonlinear alternatives (optimal-transport maps to a reference law) require
additional construction, whereas the sequential conditional CDFs are always available. Order-invariance is no free lunch here: by
Lemma~\ref{lem:energy} the symmetric whitening is a single fixed profile on the energy sphere---one
point, not the orbit that pooling exploits.
\end{remark}
\subsection{Combining uniforms: Simes, e-values, BH and e-BH}
\begin{definition}[$p$-value, e-value]
A statistic $T\ge 0$ is a $p$-\emph{value} if $\mathbb{P}_{H_0}(T\le t)\le t$ for all $t\in[0,1]$, and an
$e$-\emph{value} if $\mathbb{E}_{H_0}[T]\le 1$. By Markov's inequality, if $E$ is an e-value then
$E\ge 1/\alpha$ is a level-$\alpha$ test, and $1/E$ (capped at $1$) is a $p$-value.
\end{definition}
\begin{definition}[Simes statistic]
For per-coordinate $p$-values $\mathbf p=(p_1,\dots,p_n)$, the Simes combination is
$S(\mathbf p)=\min_{1\le i\le n} n\,p_{(i)}/i$. If the $p_i$ are i.i.d.\ $\mathrm{Unif}[0,1]$ then
$S(\mathbf p)\sim\mathrm{Unif}[0,1]$ \cite{simes1986}; thus $S$ is itself a $p$-value under independence,
exactly the situation produced by Theorem~\ref{thm:rosen}.
\end{definition}
Simes scans all sparsity levels at once: the $i=1$ term $n\,p_{(1)}$ asks ``is the single most
extreme coordinate surprising after a multiplicity correction?'', the $i=n$ term asks ``are all
coordinates a bit off?'', and the minimum takes the strongest of these views---so one statistic is
sensitive both to one large departure and to many small ones.
\begin{lemma}[Merging rules]\label{lem:merge}
Let $E_1,\dots,E_M$ be e-values and $p_1,\dots,p_M$ be $p$-values, with \emph{arbitrary} dependence.
Then: (i) the average $\frac1M\sum_m E_m$ is an e-value \cite{vovkwang2021}; (ii) twice the average
$\min(1,\frac2M\sum_m p_m)$ is a $p$-value (arithmetic-mean merging, \cite{vovkwang2020}); and
(iii) $\min(1, M\min_m p_m)$ is a $p$-value (Bonferroni / union bound).
\end{lemma}
\begin{proof}
(i) $\mathbb{E}_{H_0}[\frac1M\sum_m E_m]=\frac1M\sum_m\mathbb{E}_{H_0}[E_m]\le 1$ by linearity, with no independence
needed. (iii) $\mathbb{P}_{H_0}(\min_m p_m\le \alpha/M)\le\sum_m\mathbb{P}(p_m\le\alpha/M)\le M\cdot\alpha/M=\alpha$.
(ii) is the arithmetic-mean member ($r=1$) of the generalized-mean merging family of
\cite{vovkwang2020}, where it is shown that twice the average of arbitrarily dependent $p$-values is
again a $p$-value, and that the factor $2$ cannot be improved in general.
\end{proof}
Finally, when \emph{many} hypotheses are tested at once (Sections~\ref{sec:fdr} and~\ref{sec:fx}
screen thousands of days), controlling the per-test false-alarm rate is not enough: with $N=1{,}000$
true nulls at $\alpha=0.05$ one expects $50$ false alarms. The standard remedy controls the
\emph{false discovery rate} (FDR): the expected fraction of the rejected hypotheses that are false
alarms.
\begin{definition}[BH and e-BH]
Given $p$-values $p_1,\dots,p_N$ for $N$ hypotheses and target FDR $q$, the Benjamini--Hochberg
procedure \cite{bh1995} rejects the $k^\ast$ hypotheses with smallest $p$, where
$k^\ast=\max\{k: p_{(k)}\le qk/N\}$; it controls FDR under independence or positive dependence
(PRDS). Given e-values $E_1,\dots,E_N$, the \emph{e-BH} procedure \cite{wangramdas2022} rejects the
$k^\ast$ with largest $E$, where $k^\ast=\max\{k: E_{(k)}\ge N/(qk)\}$ ($E_{(1)}\ge\dots\ge E_{(N)}$);
it controls FDR at $q$ under \emph{arbitrary} dependence among the $E_i$, with no correction.
\end{definition}
\subsection{Order-invariant baselines}\label{sec:baselines}
The ordering problem has an obvious response: do not order. Two standard tests take it, and they
are the references against which everything in this paper is measured, so we fix them here. Both
are stated for the Gaussian null $\mathcal N(\mathbf 0,\Sigma)$ of \S\ref{sec:gaussrosen}; both are
invariant to relabelling the coordinates; and both are carried through every experiment below,
bootstrap-calibrated by the same Procedure~1 as the Rosenblatt-based tests, so that all comparisons
are made at matched size.
\begin{definition}[$\chi^2$ energy test]\label{def:chi2}
Reject for large values of the quadratic form
\begin{equation}\label{eq:chi2}
T_{\chi^2}(\mathbf{X}) \;=\; \mathbf{X}^\top\Sigma^{-1}\mathbf{X} \;\sim\; \chi^2_n \quad\text{under } H_0 .
\end{equation}
This is the Mahalanobis norm of $\mathbf{X}$, the classical omnibus test of a Gaussian null. It is order
invariant because the quadratic form does not see the labels at all.
\end{definition}
\begin{definition}[Symmetric-root test]\label{def:sym}
Whiten with the symmetric square root of Remark~\ref{rem:tri}, $\mathbf z=\Sigma^{-1/2}\mathbf{X}$, and apply
the two-sided Simes statistic of \S\ref{sec:base} to the per-coordinate $p$-values
$p_j=2\Phi(-|z_j|)$. Order invariance holds because the symmetric root commutes with permutation,
$(P_\sigma\Sigma P_\sigma^\top)^{-1/2}=P_\sigma\Sigma^{-1/2}P_\sigma^\top$, which is exactly what the
triangular Cholesky factor fails to do.
\end{definition}
The two differ in what they are sensitive to, and this will matter throughout. The $\chi^2$ test
aggregates the whitened coordinates by summing their squares, so it spreads its $n$ degrees of
freedom evenly: its power is the same for every shape of departure at fixed energy, which makes
it competitive only where the shape-sensitive tests do badly, on diffuse departures; the symmetric-root test aggregates by
a min-type rule, so it is at its best when a departure is concentrated in few coordinates. Neither
is uniformly preferable, and \S\ref{sec:headline} shows that the ranking between them inverts
between the two regimes.
There is a second, sharper reason to keep $T_{\chi^2}$ in view. By Lemma~\ref{lem:energy} the total
Mahalanobis signal energy is conserved under reordering, and \eqref{eq:chi2} is a function of that
energy alone. Proposition~\ref{prop:chi2} makes the consequence precise: the $\chi^2$ test is
\emph{orbit-flat}, its power identical under every ordering. That is simultaneously its guarantee
and its ceiling, and it is the cleanest way to see what order randomization can and cannot buy---a
test that already ignores the ordering has nothing to gain from averaging over it.
Symmetrizing a Rosenblatt-based statistic to remove the ordering is the route taken by Dovern and
Manner \cite{dovern2020} for multivariate density-forecast calibration. The present paper reaches
the same symmetry by a different road: rather than removing the ordering analytically, it draws
many orderings and pools the resulting evidence---whose $M\to\infty$ average is itself a
symmetrization (Section~\ref{sec:disc})---so the difference is \emph{what} is aggregated and
with what finite-$M$ guarantee, not aggregation versus symmetrization.
\section{Method: an order-randomized test for one multivariate hypothesis}\label{sec:method}
The construction is short: whiten $\mathbf{X}$ under each of $M$ random orderings, score each whitening, and
pool the scores. Two design choices in it are not arbitrary but are the two things that, by the
introduction's account, have to be right for pooling to pay off. We state them here and foreshadow
their justification: the base statistic is \emph{two-sided} (Section~\ref{sec:misspec} shows a
one-sided one is blind to half the signal), and the pooled test is \emph{bootstrap-calibrated}
(Section~\ref{sec:exp} shows this is both necessary for validity and the step that uncovers the power
gain). With those flagged, everything below is routine.
\subsection{Two-sided base statistics}\label{sec:base}
Given the scores $\mathbf{z}^\sigma\sim\mathcal N(\mathbf 0,I)$ under $H_0$, we form a base statistic that is
sensitive to departures in \emph{either} tail (the reason is the directional blind spot of
Section~\ref{sec:misspec}). The toolkit has three base statistics, differing in the signal
\emph{shape} they are tuned to: Simes and the mixture e-value for concentrated (sparse) departures,
Fisher for diffuse (dense) ones.
\paragraph{Two-sided Simes}
Per-coordinate two-sided $p$-values $p^\sigma_i=2\,\Phi(-|z^\sigma_i|)$, combined by
$S(\mathbf p^\sigma)=\min_i n\,p^\sigma_{(i)}/i$. Under $H_0$ this is a $p$-value.
\paragraph{Two-sided mixture e-value}
With a small grid $\mathcal T\subset(0,\infty)$ of ``bet sizes,''
\begin{equation}\label{eq:evalue}
E^\sigma=\frac1n\sum_{k=1}^n\frac1{|\mathcal T|}\sum_{\tau\in\mathcal T} e_\tau(z^\sigma_k),
\qquad e_\tau(z):=e^{-\tau^2/2}\cosh(\tau z).
\end{equation}
\begin{lemma}[Validity of the mixture e-value]\label{lem:evalue}
$\mathbb{E}_{\mathcal N(0,1)}[e_\tau(Z)]=1$ for every $\tau$, and hence $\mathbb{E}_{H_0}[E^\sigma]=1$.
\end{lemma}
\begin{proof}
$\mathbb{E}[e_\tau(Z)]=e^{-\tau^2/2}\,\mathbb{E}[\cosh(\tau Z)]=e^{-\tau^2/2}\cdot\tfrac12(\mathbb{E}[e^{\tau Z}]+\mathbb{E}[e^{-\tau Z}])
=e^{-\tau^2/2}\cdot\tfrac12(e^{\tau^2/2}+e^{\tau^2/2})=1$, using the Gaussian MGF
$\mathbb{E}[e^{tZ}]=e^{t^2/2}$. Averaging unit-mean terms over $k$ and $\tau$ preserves unit mean.
\end{proof}
One caveat must be stated plainly, because everything in Sections~\ref{sec:calib}
and~\ref{sec:fdr} turns on it. Lemma~\ref{lem:evalue} requires the scores to be
\emph{exactly} standard normal, which they are when the null---the distribution being
tested---is fully specified, as when one tests an issued forecast (Section~\ref{sec:fx}).
It is \emph{not} preserved by plugging in an estimated covariance: if the null is an unknown
Gaussian population whose $\Sigma$ is estimated from a reference sample, the plug-in scores are
heavier-tailed than normal and $E^\sigma$ need not be an e-value at all. The one-dimensional
case makes the failure vivid: with an independent unregularized variance estimate,
$Z=X/\hat\sigma$ is Student-$t$, and $\mathbb{E}[e^{-\tau^2/2}\cosh(\tau Z)]=\infty$ for every
$\tau>0$, since the $t$ distribution has no moment generating function. Under an estimated
null the mixture statistic is therefore used only as a \emph{test statistic}, calibrated into
a $p$-value by Procedure~1; the e-value reading, and everything that relies on it (Markov
thresholds, e-BH), is reserved for fully specified nulls.
The statistic $e_\tau$ is the likelihood ratio of $\mathcal N(0,1)$ against the symmetric mixture
$\tfrac12\mathcal N(\tau,1)+\tfrac12\mathcal N(-\tau,1)$; \eqref{eq:evalue} bets on a sparse,
sign-agnostic mean departure. The ``bet size'' $\tau$ sets the stakes: $e_\tau(z)$ is
increasing in $|z|$ for every $\tau$, growing like $e^{\tau|z|}$, but it starts from
$e_\tau(0)=e^{-\tau^2/2}$, so a large $\tau$ pays off faster on a large departure and loses more
on a null coordinate ($\mathbb{E}_{\mathcal N(0,1)}\log e_\tau(Z)\approx-0.13$, $-0.94$, $-2.7$ at
$\tau=1,2,3$). Averaging over the small grid $\mathcal T=\{1,2,3\}$ hedges between the two without
tuning to an unknown truth. The coordinate \emph{mixture} (a sum over $k$, not a product) is
the other design choice: a product $\prod_k e_\tau(z^\sigma_k)$ pays that ``betting tax'' on every
noise coordinate, and under a sparse alternative the few signal coordinates must repay $n-1$ of
them; for moderate departures and the larger bet sizes they do not (at $n=10$ and
$\mathrm{ncp}=12$, in the ordering that whitens the shifted coordinate last so that the whitened
profile is $(0,\dots,0,\sqrt{12})$, the $\tau=3$ product has expected log-value
$5.2-9\times2.7<0$, and it is negative under every ordering of that design), whereas in the
mixture a noise coordinate merely dilutes the evidence.
\paragraph{Two-sided Fisher}
Fisher's combination on the two-sided $p$-values,
\[
p^\sigma_F \;=\; \bar F_{\chi^2_{2n}}\!\Big(-2\sum_{k=1}^n \log p^\sigma_k\Big),
\]
where $\bar F_{\chi^2_{2n}}$ is the $\chi^2_{2n}$ survival function. Under $H_0$ the $p^\sigma_k$
are i.i.d.\ uniform, so $p^\sigma_F$ is an exact $p$-value.
\begin{remark}[Where Fisher fits---and where it cannot go]\label{rem:fisher}
Fisher aggregates \emph{every} coordinate, which makes it the dense-friendly member of the toolkit:
strong when the departure is spread across many coordinates, weak when it is concentrated
(empirically the weakest method in the study against sparse alternatives, behind even the $\chi^2$
energy test, while for dense alternatives it is competitive; Figure~\ref{fig:power}).
On \emph{one-sided} $p$-values it would inherit the directional blind spot of
Section~\ref{sec:misspec} in full, so the two-sided form above is essential. One boundary is
non-negotiable: Fisher's $\chi^2$ null law requires \emph{independent} $p$-values, so it may only be
used \emph{within} an ordering---never to merge across orderings, where the per-ordering statistics
are strongly dependent recomputations of one observation. The Fisher $p$-values
$p^{\sigma_m}_F$ can be merged across orderings with the dependence-robust rules of
Lemma~\ref{lem:merge}, and we ran both aggregation bases through every pooled combiner of
\S\ref{sec:combine} in every configuration of \S\ref{sec:exp}. The outcome is one-sided enough to
settle the design. Pooling adds almost nothing to Fisher
itself---the calibrated arithmetic-mean merge over $M=12$ orderings matches single-ordering Fisher
to within $0.005$, and the Bonferroni merge gains at most $+0.02$---so the Fisher-based combiners
inherit Fisher's shape profile: against every sparse departure they trail the corresponding
Simes-based combiner by $0.04$--$0.26$, and on the dense departure, Fisher's home ground, the
Simes-based Bonferroni still leads at known $F$ ($0.691$ vs.\ $0.642$) while the Fisher-based
$p$-merge edges its Simes counterpart ($0.623$ vs.\ $0.577$), and under heavy estimation
($N_{\mathrm{tr}}=4n$) both Fisher-based combiners lead there ($0.502$ vs.\ $0.415$ and $0.498$ vs.\
$0.421$). We therefore
aggregate only the Simes base and report Fisher as a single-ordering test. The mechanism is the
orbit theory. A sum over all coordinates is only weakly profile-sensitive---much closer to the
orbit-flat $\chi^2$ of Proposition~\ref{prop:chi2} than to the min-type Simes---so its per-ordering
values barely vary, and there is
nothing for pooling to average over---the same mechanism that makes the Jensen gap of
Proposition~\ref{prop:main} small for a profile-insensitive e-value, offered here as an
interpretation: the proposition concerns e-value averaging, whereas the Fisher $p$-values are
merged through their mean or minimum, so it does not itself deliver this conclusion. The pooling gain is a property of \emph{profile-sensitive} base
statistics, and Fisher is the demonstration of the converse.
\end{remark}
\subsection{Combining across orderings}\label{sec:combine}
Draw $\sigma_1,\dots,\sigma_M$ i.i.d.\ uniform on $S_n$, independent of the data. The per-ordering
statistics are recomputations on the \emph{same} $\mathbf{X}$ and are therefore strongly dependent;
Lemma~\ref{lem:merge} lets us combine them regardless. Write $P^\sigma$ for a generic per-ordering
base $p$-value---either the two-sided Simes $S(\mathbf p^\sigma)$ or the two-sided Fisher
$p^\sigma_F$ of \S\ref{sec:base}. We consider:
\begin{itemize}
\item \textbf{e-value averaging:} $\bar E_M=\frac1M\sum_{m=1}^M E^{\sigma_m}$; reject if
$\bar E_M\ge 1/\alpha$.
\item \textbf{$p$-merging:} $P^{\mathrm{merge}}_M=\min(1,\frac2M\sum_{m}P^{\sigma_m})$.
\item \textbf{Bonferroni-over-orders:} $P^{\mathrm{Bonf}}_M=\min(1,M\min_m P^{\sigma_m})$.
\end{itemize}
The capped forms are the nominal merged $p$-values. Whenever a combiner is
\emph{calibrated} (Procedure~1 below), the statistic
handed to the calibration is the uncapped summary---$\frac1M\sum_m P^{\sigma_m}$,
$\min_m P^{\sigma_m}$, or $\bar E_M$---of which the capped $p$-value is a monotone function below
the cap. The distinction matters under a nonrandomized rank rule: the cap places an atom at $1$, an
observation on the atom can never reject, and the calibrated test's rejection probability is
bounded by the null mass below the cap once that mass falls under $\alpha$ (with independent
coordinates and known $\Sigma$ every ordering returns the same Simes value $P$, so at $M=128$
the capped Bonferroni statistic equals $1$ with probability $127/128$ and the test could reject
at most $0.8\%$ of the time). Capping is therefore reserved for reporting a nominal merged
$p$-value. (The single-test experiments and the application calibrate the uncapped summaries; the
screening experiments of Section~\ref{sec:fdr} calibrated the capped forms, which there gives
identical decisions, because at $M=12$ the null mass below the Bonferroni cap is at least $1/12$,
far above every threshold used.)
By Lemma~\ref{lem:merge}, $\bar E_M$ is an e-value and $P^{\mathrm{merge}}_M,P^{\mathrm{Bonf}}_M$ are
$p$-values under $H_0$ for \emph{every} $M$, any dependence among orderings, and either choice of
base $p$-value; randomizing the orderings independently of the data preserves these properties
(condition on $\sigma_{1:M}$ and average). The $M\to\infty$ limit of $\bar E_M$ is the
order-invariant symmetrized e-value $\mathbb{E}_\sigma[E^\sigma(\mathbf{X})\mid\mathbf{X}]$; finite $M$ is a Monte-Carlo
approximation. We evaluated the pooled combiners on both bases; the Simes base won in every
sparse configuration at both training sizes, while on the dense departure the Fisher base leads
under heavy estimation and, at known $F$, for the $p$-merge (Remark~\ref{rem:fisher}), so
the experiments report Simes-based pooling and Fisher as a single-ordering test.
\subsection{Calibration under an estimated null}\label{sec:calib}
In practice $\Sigma$ is unknown and estimated as $\hat\Sigma$ from a reference sample of size $N_{\mathrm{tr}}$.
The plug-in scores $\hat\mathbf{z}^\sigma=\hat L_\sigma^{-1}\mathbf{X}_\sigma$ are then only approximately
$\mathcal N(\mathbf 0,I)$ under $H_0$, so the nominal thresholds of Section~\ref{sec:combine} are no
longer exact and the test is miscalibrated (Section~\ref{sec:exp}). We restore level by a
\emph{re-estimating parametric bootstrap}. The idea is simple: since we cannot compute the exact
null distribution of our statistic under estimation noise, we \emph{manufacture} it---simulate many
worlds in which the fitted model is true, run the \emph{entire} pipeline (including re-estimating the
covariance from a fresh training sample, so the simulated worlds carry the same estimation noise as
the real one) in each world, and use the simulated statistics' quantile as the rejection threshold.
This is the standard remedy for goodness-of-fit testing with estimated parameters
\cite{stute1993,genestremillard2008}, an instance of bootstrap calibration (prepivoting) in the sense
of \cite{beran1988}; see \cite{davisonhinkley1997} for parametric and Monte-Carlo bootstrap tests
generally.
\begin{procedure}{1: Bootstrap calibration of a combiner $T$}
Given the fitted model $\hat\Sigma$, target level $\alpha$, orderings $\sigma_{1:M}$, replicate count
$B$:\\
1. For $b=1,\dots,B$: draw a fresh ``training'' sample of size $N_{\mathrm{tr}}$ from $\mathcal N(\mathbf 0,\hat\Sigma)$
and re-estimate $\hat\Sigma^{\ast}_b$; draw one calibration vector $\mathbf{X}^\ast_b\sim\mathcal N(\mathbf 0,\hat\Sigma)$;
compute $T^\ast_b=T(\mathbf{X}^\ast_b;\hat\Sigma^{\ast}_b)$ using the \emph{re-estimated} whitening.\\
2. Rank-based Monte-Carlo $p$-value: for a $p$-value combiner (small values are evidence),
$p_{\mathrm{MC}}=\bigl(1+\#\{b:T^\ast_b\le T_{\mathrm{obs}}\}\bigr)/(B+1)$; for an e-value
combiner (large values are evidence),
$p_{\mathrm{MC}}=\bigl(1+\#\{b:T^\ast_b\ge T_{\mathrm{obs}}\}\bigr)/(B+1)$. Reject if
$p_{\mathrm{MC}}\le\alpha$.
\end{procedure}
Two remarks on Procedure~1. First, on what is claimed. The observed statistic
$T_{\mathrm{obs}}=T(\mathbf{X};\hat\Sigma)$ whitens by the fixed fitted covariance while each replicate
$T^\ast_b=T(\mathbf{X}^\ast_b;\hat\Sigma^\ast_b)$ whitens by a re-estimated one, so even under the fitted
null $\mathbf{X}\sim\mathcal N(\mathbf 0,\hat\Sigma)$ the two are \emph{not} identically distributed and
$p_{\mathrm{MC}}$ is not an exact $p$-value for that null (in one dimension, a fixed-variance
score is $\mathcal N(0,1)$ while a re-estimated one is Student-$t$; calibrating the former against
the latter, two-sided at ten degrees of freedom, rejects with probability $0.026$, not $0.05$,
however large $B$). Nor is the fitted null the target. The re-estimating bootstrap aims at the
\emph{estimated-population} null: the fitted model stands in for the unknown population, and each
replicate re-enacts the whole pipeline---estimation from $N_{\mathrm{tr}}$ draws, then whitening of an
independent vector by that estimate---so that $\{T^\ast_b\}$ mimics the sampling distribution of
$T_{\mathrm{obs}}$ when both the training sample and $\mathbf{X}$ come from the population. That
approximation is asymptotic in $N_{\mathrm{tr}}$, under the regularity conditions of the parametric-bootstrap
literature \cite{stute1993,genestremillard2008,beran1988}; we claim no finite-sample exactness,
and Section~\ref{sec:calib-exp} measures how close the realized size comes at $N_{\mathrm{tr}}=4n$ and $8n$.
(The implementation adds a fixed ridge $10^{-3}I$ to every covariance estimate for numerical
stability; the asymptotic argument presumes a ridge that vanishes with $N_{\mathrm{tr}}$, and at the sizes
used its effect is negligible.)
The rank form of step~2 is a separate, finite-$B$ matter. When the simulated null is exact---a
fully specified $F$, where step~1 reduces to drawing $\mathbf{X}^\ast_b$ from $F$ and whitening it by the
same fixed factor as the data, the situation of the issued-forecast null in
Section~\ref{sec:fx}---$T_{\mathrm{obs}},T^\ast_1,\dots,T^\ast_B$ are exchangeable and
$p_{\mathrm{MC}}$ is an exact finite-sample $p$-value \cite{davisonhinkley1997}; an interpolated
empirical $\alpha$-quantile of $\{T^\ast_b\}$ is instead not exact in general at finite $B$,
and can be anti-conservative. For a
statistic that is uniform under the null, such as a $p$-value, the expected rejection probability
of the interpolated rule is $((B-1)\alpha+1)/(B+1)$, e.g.\ $0.0545$ at $B=199$ and $\alpha=0.05$
(for other continuous statistics the bias is of the same order but not distribution-free); it was
visible in earlier drafts of our own experiments and is removed by the rank rule, whose only cost
is the granularity $1/(B+1)$. Second, re-estimating
$\hat\Sigma^\ast_b$ \emph{inside} each replicate is what makes the bootstrap null carry
the same estimation noise as the real pipeline; this is essential for the calibration to remove the
estimation-induced bias (Section~\ref{sec:exp}).
\subsection{The full procedure}
\begin{procedure}{2: Order-randomized Rosenblatt test (single statistic)}
Input: data vector $\mathbf{X}$; fitted $\hat\Sigma$; level $\alpha$; $M$, $B$, base statistic, combiner.\\
1. Draw $\sigma_1,\dots,\sigma_M\stackrel{\text{i.i.d.}}{\sim}\mathrm{Unif}(S_n)$ and form
$\hat L_{\sigma_m}=\mathrm{chol}(P_{\sigma_m}\hat\Sigma P_{\sigma_m}^\top)$.\\
2. Whiten: $\hat\mathbf{z}^{\sigma_m}=\hat L_{\sigma_m}^{-1}P_{\sigma_m}\mathbf{X}$; compute the base statistic
(\S\ref{sec:base}) for each $m$.\\
3. Combine across orderings (\S\ref{sec:combine}) to get $T_{\mathrm{obs}}$.\\
4. Calibrate the threshold by Procedure~1 and reject accordingly.
\end{procedure}
The procedure is for a single multivariate hypothesis. When a collection of such hypotheses is
screened simultaneously, the multiplicity layer depends on which null is in force: under a fully
specified null the per-hypothesis e-values $\bar E_M^{(j)}$ are valid and e-BH applies directly
(Section~\ref{sec:fx}); under an estimated null the per-hypothesis outputs are the calibrated
$p$-values of Procedure~1 and the layer is BH (Section~\ref{sec:fdr}). This layer is a standard add-on and is treated as a secondary
result in Section~\ref{sec:fdr}.
\section{Theory: power as a profile functional on the orbit of orderings}\label{sec:theory}
Why should pooling orderings help at all? Reordering shuffles the data through different whitenings,
but it cannot create or destroy information, so any gain must come from \emph{how} a fixed departure
is presented to the test, not from how much departure there is. This section makes that intuition
exact. We show that reordering moves the alternative around a sphere of fixed radius---the
\emph{orbit}---changing only the shape of the signal across coordinates (Lemma~\ref{lem:energy}, with
Proposition~\ref{prop:kl} as its general, non-Gaussian form: the transform's total information is
conserved under reordering); that
a pure energy test is therefore blind to ordering and cannot be helped by pooling
(Proposition~\ref{prop:chi2}); but that a shape-sensitive test sees a lucky or unlucky view under each
ordering, so averaging the views strictly improves the expected evidence by exactly the orbit's
variability---which is zero under independence (Proposition~\ref{prop:main}) and is large,
empirically, under strong dependence.
A short calculation then gives the optimal number of orderings for the max-type combiner
(Proposition~\ref{prop:optM}).
Throughout this section $\mathbf{X}\sim\mathcal N(\boldsymbol{\mu},\Sigma)$ (the alternative), and we assume the
whitening uses the true $\Sigma$ (the role of estimation is studied empirically in
Section~\ref{sec:exp}). Under ordering $\sigma$, $\mathbf{z}^\sigma\sim\mathcal N(\mathbf{m}^\sigma,I)$ with whitened
mean $\mathbf{m}^\sigma=L_\sigma^{-1}P_\sigma\boldsymbol{\mu}$.
\begin{lemma}[Energy invariance]\label{lem:energy}
For every $\sigma\in S_n$, $\ \|\mathbf{m}^\sigma\|^2=\boldsymbol{\mu}^\top\Sigma^{-1}\boldsymbol{\mu}$, and $\|\mathbf{z}^\sigma\|^2=\mathbf{X}^\top\Sigma^{-1}\mathbf{X}$
is independent of $\sigma$.
\end{lemma}
\begin{proof}
$\Sigma_\sigma=P_\sigma\Sigma P_\sigma^\top=L_\sigma L_\sigma^\top$ gives
$L_\sigma^{-1}\Sigma_\sigma L_\sigma^{-\top}=I$, hence $\Sigma_\sigma^{-1}=L_\sigma^{-\top}L_\sigma^{-1}$.
Therefore
\[
\begin{aligned}
\|\mathbf{m}^\sigma\|^2
&=(P_\sigma\boldsymbol{\mu})^\top L_\sigma^{-\top}L_\sigma^{-1}(P_\sigma\boldsymbol{\mu})
=(P_\sigma\boldsymbol{\mu})^\top\Sigma_\sigma^{-1}(P_\sigma\boldsymbol{\mu})\\
&=\boldsymbol{\mu}^\top P_\sigma^\top(P_\sigma\Sigma P_\sigma^\top)^{-1}P_\sigma\boldsymbol{\mu}
=\boldsymbol{\mu}^\top\Sigma^{-1}\boldsymbol{\mu},
\end{aligned}
\]
using $P_\sigma^{-1}=P_\sigma^\top$. The same computation with $\mathbf{X}$ replacing $\boldsymbol{\mu}$ gives the second
claim.
\end{proof}
Thus reordering is a \emph{norm-preserving} reparametrization: the family
$\{\mathbf{m}^\sigma:\sigma\in S_n\}$ lies on a sphere of fixed radius $\|\mathbf{m}\|=(\boldsymbol{\mu}^\top\Sigma^{-1}\boldsymbol{\mu})^{1/2}$,
which we call the \emph{orbit} of the alternative. What changes across $\sigma$ is the \emph{profile}
of $\mathbf{m}^\sigma$ (how the fixed energy is distributed across coordinates) and its sign pattern. The
quantity $\boldsymbol{\mu}^\top\Sigma^{-1}\boldsymbol{\mu}$ is the squared \emph{Mahalanobis distance} of the alternative from
the null---the natural ``how far is the truth from the model'' scale that accounts for correlations;
we call it the signal's \emph{energy}. Lemma~\ref{lem:energy} says reordering can neither create nor
destroy energy; it only re-apportions it among the whitened coordinates. One ordering may pile the
energy onto a single coordinate (easy for a sparse-friendly test to see) while another smears the
same energy thinly over all coordinates (nearly invisible to it): same distance from the model,
very different visibility. The effect is large. With $n=50$ independent standard normal coordinates
and a mean departure of fixed total $\mathrm{KL}=4$, the two-sided Simes test of \S\ref{sec:base}
detects the departure $36\%$ of the time at $\alpha=0.05$ when the whole budget sits on one
coordinate and $10\%$ when it is spread evenly over all fifty (one-sided: $43\%$ and $17\%$)---equally informative alternatives, at the same distance from the null, differing only in
profile.
The energy lemma is the Gaussian instance of a fully general conservation law, which does not even
require the null to be Gaussian: however the coordinates are ordered, the transform's departure from
uniformity carries exactly the model's error, no more and no less.
\begin{proposition}[KL invariance]\label{prop:kl}
Suppose the truth is $G$ but the transform uses $F$ (both absolutely continuous, with the
conditional distribution functions of $F$ strictly increasing on their supports, so that
$R_\sigma^F$ is almost surely bijective), and let
$\mathbf{U}=R_\sigma^F(\mathbf{X})$ for $\mathbf{X}\sim G$. Then $\mathrm{KL}(\mathcal L(\mathbf{U})\,\|\,\mathrm{Unif}[0,1]^n)=\mathrm{KL}(G\,\|\,F)$, and in
particular this is invariant to $\sigma$.
\end{proposition}
\begin{proof}
$R_\sigma^F$ is a measurable bijection and $R_\sigma^F(\mathbf{X})\sim\mathrm{Unif}[0,1]^n$ when $\mathbf{X}\sim F$
(Theorem~\ref{thm:rosen}). KL divergence is invariant under a common bijective transformation of its
two arguments; applying $R_\sigma^F$ to $G$ gives $\mathcal L(\mathbf{U})$ and to $F$ gives $\mathrm{Unif}[0,1]^n$,
so $\mathrm{KL}(\mathcal L(\mathbf{U})\,\|\,\mathrm{Unif})=\mathrm{KL}(G\,\|\,F)$. The right side does not involve $\sigma$.
\end{proof}
(The Kullback--Leibler divergence $\mathrm{KL}(G\,\|\,F)$ measures, in information units, how distinguishable
the truth $G$ is from the model $F$; it is the best possible average log-likelihood ratio, so it
bounds what \emph{any} test can achieve.) Together, Lemma~\ref{lem:energy} and
Proposition~\ref{prop:kl} say that ordering changes \emph{how} a fixed departure is presented to the
test---never how much departure there is. Any gain from pooling orderings must therefore come from
presentation, and any loss from a single ordering is a presentation failure. The next two results
make both directions precise.
\begin{proposition}[Orbit-flatness of energy tests]\label{prop:chi2}
Any test whose statistic depends on $\mathbf{z}^\sigma$ only through $\|\mathbf{z}^\sigma\|^2$ has power independent
of $\sigma$: the statistic is literally the same number under every ordering, so ordering, and any
combination of orderings applied to the statistic itself, changes nothing. In particular the
$\chi^2$ test based on $\mathbf{X}^\top\Sigma^{-1}\mathbf{X}$ is unaffected by ordering, and its power under a
mean-shift alternative depends on the departure only through its energy
$\boldsymbol{\mu}^\top\Sigma^{-1}\boldsymbol{\mu}$---it is \emph{flat} across departure shapes at fixed energy.
\end{proposition}
\begin{proof}
By Lemma~\ref{lem:energy}, $\|\mathbf{z}^\sigma\|^2$ is the same for all $\sigma$; hence the statistic and
its distribution do not depend on $\sigma$, and under $\mathbf{X}\sim N(\boldsymbol{\mu},\Sigma)$ the statistic is
noncentral $\chi^2_n$ with noncentrality $\boldsymbol{\mu}^\top\Sigma^{-1}\boldsymbol{\mu}$, whatever the profile of $\boldsymbol{\mu}$.
\end{proof}
(One clarification: this concerns combining the identical statistic. Pushing $M$ duplicated
copies through a \emph{penalized} nominal merger is a different operation and can only lose---
Bonferroni's $M\min_m p$ applied to $M$ equal $p$-values multiplies the $p$-value by $M$;
calibrated combination of the duplicates, on the uncapped summaries of \S\ref{sec:combine},
is exactly the original test.)
Profile-sensitive tests behave differently. The headline positive result, stated for the e-value
combiner, uses expected log-evidence (growth rate)---the canonical optimality criterion for e-values
\cite{grunwald2024}. In the betting reading, $\mathbb{E}_Q[\log E]$ is the long-run exponential growth rate
of repeatedly staking one's evidence capital on the bet $E$ when the truth is $Q$; maximizing it is
the evidential analogue of maximizing power, and it is the natural criterion here because log-evidence,
unlike fixed-level power, behaves cleanly under the averaging that pooling performs.
\begin{proposition}[Order-averaging dominates a single ordering]\label{prop:main}
Let $\sigma\sim\mathrm{Unif}(S_n)$ be drawn independently of $\mathbf{X}$, and for each $\sigma$ let $E^\sigma(\mathbf{X})\ge0$
be an e-value, $\mathbb{E}_{H_0}[E^\sigma]=1$. Let $\bar E(\mathbf{X})=\mathbb{E}_\sigma[E^\sigma(\mathbf{X})\mid\mathbf{X}]$, and for
(ii) assume $E^\sigma(\mathbf{X})>0$ $Q$-a.s.\ with $\mathbb{E}_Q|\log E^\sigma(\mathbf{X})|<\infty$ for each $\sigma$
and $\mathbb{E}_Q|\log\bar E(\mathbf{X})|<\infty$ (satisfied by the mixture e-value \eqref{eq:evalue}, which is
bounded below by $e^{-\max\mathcal T^2/2}/(n|\mathcal T|)>0$ and has Gaussian log-moments). Then:
\begin{enumerate}
\item[(i)] $\bar E$ is an e-value: $\mathbb{E}_{H_0}[\bar E]=1$.
\item[(ii)] for any alternative $Q$,
$\ \mathbb{E}_Q[\log\bar E(\mathbf{X})]\ \ge\ \mathbb{E}_\sigma\,\mathbb{E}_Q[\log E^\sigma(\mathbf{X})]$, with equality iff $E^\sigma(\mathbf{X})$ is
$Q$-a.s.\ constant in $\sigma$;
\item[(iii)] if $\Sigma$ is diagonal and the base e-value is exchangeable (permutation-symmetric in the
coordinatewise scores), then $E^\sigma(\mathbf{X})$ is constant in $\sigma$ and (ii) holds with equality:
order combination yields no growth-rate gain.
\end{enumerate}
\end{proposition}
\begin{proof}
(i) Fubini and the e-value property: $\mathbb{E}_{H_0}[\bar E]=\mathbb{E}_\sigma\mathbb{E}_{H_0}[E^\sigma]=1$.
(ii) Condition on $\mathbf{X}$; $\log$ is concave, so Jensen gives
$\log\mathbb{E}_\sigma[E^\sigma(\mathbf{X})\mid\mathbf{X}]\ge\mathbb{E}_\sigma[\log E^\sigma(\mathbf{X})\mid\mathbf{X}]$, with equality iff
$E^\sigma(\mathbf{X})$ is a.s.\ constant over $\sigma$ given $\mathbf{X}$; take $\mathbb{E}_Q$ and apply Fubini.
(iii) With diagonal $\Sigma$, $\mathbf{z}^\sigma$ is a permutation of the fixed coordinatewise standardized
scores, and an exchangeable base e-value---e.g.\ the coordinate-mixture \eqref{eq:evalue}---is
invariant to the order, so $E^\sigma(\mathbf{X})$ is constant in $\sigma$.
\end{proof}
The content of (ii) is elementary but worth unpacking. Conditional on the data, the per-ordering
e-values $E^\sigma(\mathbf{X})$ are a collection of numbers---one view of the same observation per
ordering---and the theorem compares two policies: average the numbers first and take the log
(pooling), or take the log first and average (a single random ordering, in expectation). Because the
logarithm is concave, averaging first always wins, and the margin---the \emph{Jensen gap}---is
governed by how much the views disagree. Under independence the views are all identical
(part~(iii)) and pooling buys nothing; under strong dependence different orderings whiten the same
observation very differently, the views disagree a lot, and pooling's margin is large. This
is the mechanism we read behind the empirical pattern that the boost grows with $\rho$ in the
designs of Section~\ref{sec:exp}, but it should not be overread. The proposition is a statement
about expected log-evidence: it does not imply that the \emph{calibrated power} gain is positive
or that it grows with dependence, and the Jensen gap itself need not be monotone in $\rho$---for
$n=2$ and $\rho\to1$ the two Cholesky whitenings agree up to the signs of the coordinates, so a
sign-invariant base statistic such as \eqref{eq:evalue} has a Jensen gap that vanishes at fixed
energy. The empirical sections verify that the power boost indeed materializes, quantify it, and
map where it is large or negligible.
\begin{remark}[When expected log-evidence controls power]\label{rem:bridge}
Write $\mu=\mathbb{E}_Q[\log E]$ and $s^2=\operatorname{Var}_Q(\log E)$, and reject when $E\ge1/\alpha$,
equivalently $\log E\ge t:=\log(1/\alpha)$---the threshold Markov's inequality licenses for any
e-value. Two bridges lead from $\mu$ to power.
\emph{(i) Repeated observations.} Suppose $T$ i.i.d.\ observations from $Q$ are tested and
their e-values multiplied---the natural use of evidence, valid with no correction---and that
$\mu>0$ with $\mathbb{E}_Q|\log E|<\infty$. Then $\log\prod_{i\le T}E_i$ has mean $T\mu$, and by the law of
large numbers the number of observations needed to reach any fixed target power $1-\varepsilon$,
$T_{\varepsilon}(\alpha):=\min\{T:\mathbb{P}_Q(\sum_{i\le T}\log E_i\ge t)\ge1-\varepsilon\}$,
satisfies $T_{\varepsilon}(\alpha)\sim t/\mu$ as $\alpha\downarrow0$ (i.e.\ $t\to\infty$).
Asymptotically, then, the required sample size is \emph{inversely proportional} to expected
log-evidence: a procedure with the larger $\mu$ needs proportionally less data for the same power.
The approximation is informative well before the limit. In the design of
Figure~\ref{fig:power} with the dependence swept over $\rho\in\{0.2,0.5,0.8\}$, pooling raises $\mu$
by factors $1.027$, $1.108$ and $1.114$, and the measured reduction in the number of observations
needed for power $0.80$ (sample sizes interpolated linearly between consecutive integers $T$)
is by factors $1.035$, $1.117$ and $1.125$---agreement to within one percent
at $T$ between two and three.
\emph{(ii) A single observation, strong signal.} At $T=1$, Cantelli's one-sided Chebyshev
inequality \cite{cantelli1928,blm2013} gives
\[
\mathbb{P}_Q(\log E\ge t)\ \ge\ \frac{(\mu-t)_+^2}{s^2+(\mu-t)^2},
\]
increasing in $\mu$ and decreasing in $s^2$: raising expected log-evidence without inflating its
dispersion raises a lower bound on power. The bound is informative only when $\mu>t$---when the
expected evidence already exceeds the level one is testing at, $2.996$ nats at $\alpha=0.05$. That
condition fails through most of our designs, and this is the honest reason our fixed-level gains are
a fraction of the evidence gains rather than proportional to them. Where it holds it orders the two
procedures correctly, and pooling is then measured to reduce $s^2$ as well as raise $\mu$, so both
terms move the same way: at $\rho=0.8$ and $\mathrm{ncp}=20$ the bound rises from $0.12$ to $0.30$
under pooling while true power rises from $0.59$ to $0.72$, and at $\mathrm{ncp}=30$ from $0.64$ to
$0.72$ against true power $0.90$ to $0.95$.
\end{remark}
\begin{proposition}[Optimal number of orderings for Bonferroni]\label{prop:optM}
Fix the data at an alternative configuration and suppose that, over the random draw of $\sigma$, the
per-ordering Simes $p$-value equals a small value $p_g$, $0<p_g\le\alpha$, with probability
$q\in(0,1)$ (a ``concentrating'' ordering) and is $>\alpha$ otherwise. Draw $M$ i.i.d.\ orderings. Then the
Bonferroni-over-orders test (reject iff $M\min_m p^{\sigma_m}\le\alpha$) has power
\[
\beta(M)=\begin{cases}1-(1-q)^M, & M\le \alpha/p_g,\\[2pt] 0, & M>\alpha/p_g,\end{cases}
\]
so $\beta$ is increasing on the feasible region and the optimal count is
$M^\ast=\lfloor\alpha/p_g\rfloor$; the optimal power $1-(1-q)^{M^\ast}$ reaches
$1-\varepsilon$ once $M^\ast\ge\lceil\log\varepsilon/\log(1-q)\rceil$, i.e.\ roughly
$p_g\le\alpha\log(1-q)/\log\varepsilon$
(at $M^\ast\approx1/q$ it is only about $1-e^{-1}\approx0.63$ for small $q$).
\end{proposition}
\begin{proof}
A ``bad'' ordering has $p>\alpha\ge\alpha/M$ for $M\ge1$, so it never triggers; only ``good''
orderings (value $p_g$) can, and they do iff $p_g\le\alpha/M$, i.e.\ $M\le\alpha/p_g$. Conditional on
that, rejection occurs iff at least one of the $M$ i.i.d.\ orderings is good, with probability
$1-(1-q)^M$, which is increasing in $M$. Hence within the feasible region the largest admissible $M$
is optimal, $M^\ast=\lfloor\alpha/p_g\rfloor$.
\end{proof}
\begin{remark}\label{rem:dome}
Proposition~\ref{prop:optM} is a statement about the \emph{nominal} $\alpha/M$ threshold, and its
hard cliff is an artifact of the stylized two-point model. With real alternatives the per-ordering
$p$-value has a continuous distribution with substantial mass near zero, so as $M$ grows the best of
$M$ draws shrinks almost as fast as the threshold tightens, and the cliff softens into a shallow
dome: power rises steeply, peaks at a moderate $M^\ast$, then declines slowly. The decline is
eventual for every fixed signal: the smallest per-ordering $p$-value over the finite set of all
orderings, $p_{\min}(\mathbf{X})=\min_{\sigma\in S_n}p^\sigma(\mathbf{X})$, is positive almost surely, so
$M\min_{m\le M}p^{\sigma_m}\ge Mp_{\min}\to\infty$ and the nominal Bonferroni test eventually
loses all power as $M\to\infty$. Within the range of $M$ we tested, the decline was visible only for
weak signals (in a hub design with a weak one-coordinate shift, power peaks at $M=2$ with $0.933$ and
decays to $0.877$ by $M=512$; at stronger signals the curve saturates at $1$ and no decline is
visible up to $M=512$). The averaging combiner shows no dome in the same sweep---its nominal-threshold power
was nondecreasing in $M$ in every configuration we ran; we note that this is an empirical
statement about these designs, not a monotonicity theorem for fixed-level power, which
Proposition~\ref{prop:main} does not supply. None of this survives calibration, which is what
we recommend in practice: with
bootstrap-calibrated thresholds the $\times M$ penalty is absorbed into the critical value and
calibrated power is essentially \emph{flat} in $M$ for every combiner. The
practical reading is therefore simply that a modest $M$ suffices, and that in every calibrated
configuration we ran, ``too many'' orderings did not hurt.
\end{remark}
\begin{remark}[Two orthogonal levers; calibration unlocks the orbit lever]
Improving a profile-sensitive test has two independent axes: (a) choosing a profile- and sign-robust
base statistic (the two-sided bet)---a property of the \emph{test}; and (b) combining across
orderings to escape an unfavorable profile---a property of the \emph{aggregation}. With nominal
thresholds, the combiners carry multiplicity/conservativeness penalties (Bonferroni's $\times M$, the
merge's $\times2$, the e-value's Markov slack) that can erase the orbit advantage and make a single
ordering look as good or better. Bootstrap calibration (\S\ref{sec:calib}) removes these
penalties---and is independently required for validity under estimation---so the richer combined
statistic realizes its advantage. This is the mechanism behind the headline result.
\end{remark}
\section{The base statistic must be two-sided}\label{sec:misspec}
The Rosenblatt scores are standardized \emph{conditional} quantities, and the sign with which a
departure surfaces in them depends on the ordering: in Example~\ref{ex:bivar} the same observation
produced $U=1.8\times10^{-5}$ under one ordering and $U=0.999998$ under the other---the identical
surprise, delivered once to the lower and once to the upper tail. So under reordering, wrong-tail
signal is the norm, not an edge case, and the classical one-sided Simes combination---which fixes
one tail in advance; in Figure~\ref{fig:blindspot} it is the upper tail, $p_k=1-\Phi(z_k)$, so
only large \emph{positive} scores register as evidence---is unusable when the sign of the
departure is not known in advance. Its failure is worse than
lost power: a
wrong-sign departure sends the affected $p$-value toward $1$, where the statistic quietly
discards it while the critical constants $\alpha i/n$ remain scaled for all $n$ coordinates, so the
rejection probability falls strictly \emph{below} the level and stays there however extreme the data
become (as the shift grows it converges to $\alpha(n-1)/n$, the level carried by the $n-1$ null
coordinates---an exact consequence of the Simes null identity, and visible as the plateau at
$0.047\approx0.05\times19/20$ in Figure~\ref{fig:blindspot}A). Nor can the failure be repaired by
simply testing the other tail: whitened departures generically carry \emph{both} signs at once (the
hub example of \S\ref{sec:exp} whitens to one positive and nineteen negative components), and
Figure~\ref{fig:blindspot}B makes the point cleanly by shifting half the coordinates by $+c/\sqrt n$
and half by $-c/\sqrt n$---either one-sided orientation then sees at most half the signal (by
symmetry both give the same power curve), and the two-sided test dominates uniformly in $c$. We
therefore fold every coordinate two-sided, $p_k=2\,\Phi(-|z_k|)=2\min(U_k,1-U_k)$: the fold acts
coordinatewise, so the folded values remain exactly i.i.d.\ uniform under $H_0$, and evidence
registers whichever tails the whitening chooses.
\begin{figure}[t]\centering
\includegraphics[width=\linewidth]{onesided_blindspot.png}
\caption{One- vs.\ two-sided Simes ($n=20$; independent null, where all orderings coincide). (A)
Signed single-coordinate shift: the one-sided test is blind to $\delta<0$ and sits strictly below
its level (dashed: the limit $\alpha(n-1)/n$); the two-sided fold restores symmetric detection at
negligible cost for $\delta>0$. (B) Mixed-sign shift, $\pm c/\sqrt n$ on half the coordinates each:
no one-sided orientation can see more than half the signal, and the two-sided test dominates
uniformly---the deficiency is structural, not a matter of choosing the right tail.}
\label{fig:blindspot}
\end{figure}
\section{Calibration and the payoff: Gaussian experiments}\label{sec:exp}
With the orbit established (Section~\ref{sec:theory}) and the base test fixed
(Section~\ref{sec:misspec}), only the third ingredient remains. We use the experiments to do two
things at once: show that under an estimated null the naive pooled test is
miscalibrated and that the bootstrap repairs it (\S\ref{sec:calib-exp})---and then, on the level
playing field calibration creates, measure the payoff the theory motivates: a power gain over the
expected single random ordering, one that grows with dependence in these designs, together with a shape-robustness
advantage over order-invariant tests (\S\ref{sec:headline}); a gain that persists under
estimation. Practical guidance on the choice of combiner is deferred to \S\ref{sec:which}, since
it turns on evidence from both this section and the screen of \S\ref{sec:fdr-main}.
\paragraph{Protocol}
Unless stated otherwise: the null is $\mathcal N(\mathbf 0,\Sigma)$ with $\Sigma$ equicorrelated,
$\Sigma_{ij}=\rho$ for $i\neq j$ and $1$ on the diagonal; the alternative shifts a single coordinate
(``sparse'') or all coordinates equally (``dense''), with the squared Mahalanobis distance of the
alternative from the null---equivalently the \emph{noncentrality parameter}
$\mathrm{ncp}:=\boldsymbol{\mu}^\top\Sigma^{-1}\boldsymbol{\mu}$, the ``signal energy'' of Lemma~\ref{lem:energy} and the
noncentrality of the $\chi^2$ oracle---fixed at
$\mathrm{ncp}=12$ in \emph{every} experiment so that detectability is
comparable across configurations (holding energy fixed is what makes ``sparse vs.\ dense'' or
``$\rho=0.2$ vs.\ $\rho=0.8$'' comparisons meaningful: every alternative is equally far from the
null, only its shape changes); level $\alpha=0.05$; bet grid $\mathcal T=\{1,2,3\}$; $M=12$ random
orderings; calibration by Procedure~1 with $B$ replicates. We estimate $\Sigma$ by the sample
covariance from $N_{\mathrm{tr}}$ null vectors (``heavy estimation'' uses $N_{\mathrm{tr}}=4n$; ``known $F$'' uses the true
$\Sigma$). Each reported number averages over $R$ Monte-Carlo realizations of the fitted model, each
with a batch of fresh test vectors; all seeds are fixed integers. Power is compared \emph{at matched
calibrated level}, so differences are genuine power differences rather than calibration artifacts.
For the headline runs, $B\in\{199,299\}$, $R\in\{150,300\}$; full settings are in the scripts
(Appendix~\ref{app:repro}).
Two design points deserve emphasis. \emph{First}, the single-ordering baseline is drawn \emph{afresh
in every realization}, so its average estimates the expected power of an arbitrarily chosen
ordering---the quantity the theory speaks to. This matters: the power of one fixed ordering depends
strongly on which coordinates precede the departure, so freezing a single unlucky permutation
across realizations can inflate the apparent gain from pooling several-fold. Because the
single-ordering baseline is a Simes $p$-value while the e-value average pools mixture e-values,
we also carry the calibrated mixture e-value of one random ordering (``single-ordering e-value''),
so that the e-value path can be compared like with like. \emph{Second}, the two order-invariant references of
\S\ref{sec:baselines}---the $\chi^2$ energy test and the symmetric-root test---are carried through
every table, both bootstrap-calibrated by the same Procedure~1 with $\hat\Sigma$ in place of $\Sigma$. The
alternatives we study are \emph{location} (mean-shift) departures; a systematic evaluation against
scale, dependence, and shape alternatives is left to future work.
\subsection{Calibration restores level---and is necessary}\label{sec:calib-exp}
Under an estimated covariance the na\"ive plug-in test is badly miscalibrated, in \emph{opposite}
directions for different combiners (Figure~\ref{fig:boot}A), and each direction has a clear cause.
The min-type statistics are \emph{anti}-conservative---at $N_{\mathrm{tr}}=4n$ the single-ordering Simes
realizes size $0.20$, $\chi^2$ and Fisher $0.28$, the symmetric-root test $0.19$---because whitening
with a noisy $\hat\Sigma$ leaves residual structure in the scores and a min-type statistic latches onto
whichever coordinate the estimation error happens to make extreme. The e-value average errs the
other way, spending $0.029$ of its level at $N_{\mathrm{tr}}=4n$ and $0.001$ at known $F$: its Markov threshold
$1/\alpha$ is loose, so most of the level goes unused. Bonferroni-over-orders and the $p$-merge are
conservative once the estimation noise that inflates every min-type statistic subsides
(from $N_{\mathrm{tr}}=8n$ on) for a related reason, their union bound and factor $2$ ignoring that the $M$ orderings
are strongly dependent recomputations of one observation. The re-estimating bootstrap (Procedure~1)
brings \emph{every} method to size $\approx\alpha$ at every training size (Figure~\ref{fig:boot}B,
drawn on the same scale as panel A). Calibration is therefore both a validity requirement and the
enabler of a fair power comparison: without it, the power ranking of the methods would mostly
reflect how much level each one wastes. Size control after calibration is not confined to this
configuration: across all $64$ configuration$\times$method cells of the wider robustness sweep the
calibrated size averages $0.050$ (range $0.042$--$0.059$).
Figure~\ref{fig:boot} shows the whole pipeline---miscalibration, repair, and payoff---in a single
setting: $n=20$, $\rho=0.5$, two same-sign shifted coordinates (one step off the one-coordinate
knife edge where the symmetric-root test is at its best, \S\ref{sec:headline}). At matched size the
pooled $p$-merge outperforms both order-invariant references of \S\ref{sec:baselines}, the single
random ordering and the single-ordering Fisher test, at every training size and at known $F$.
\begin{figure}[t]\centering
\includegraphics[width=\linewidth]{bootstrap_recalibrated.png}
\caption{One setting, whole pipeline ($n=20$, $\rho=0.5$, two same-sign shifted coordinates,
$\mathrm{ncp}=12$, $M=12$). (A) Na\"ive thresholds are miscalibrated in both directions under an
estimated null ($0.28$ for $\chi^2$ and Fisher down to $0.001$ for the e-value average). (B) After
bootstrap calibration every method sits at $\alpha$ at every training size (same $y$-scale as A).
(C) On that level playing field the pooled p-merge outperforms every order-invariant reference at
every training size and at known $F$ ($0.437$, $0.467$, $0.495$ against the best reference, the
symmetric-root test, at $0.410$, $0.438$, $0.475$, which $\chi^2$ ties at known $F$). The dashed line is the calibrated
\emph{single-ordering} mixture e-value; its distance to the e-value average ($+0.024$, $+0.038$,
$+0.032$) is the gain from aggregation alone, the base statistic held fixed.}
\label{fig:boot}
\end{figure}
\subsection{Shape robustness, and a gain that grows with dependence}\label{sec:headline}
The practical case for pooling is visible in Figure~\ref{fig:power}, and it is a
statement about \emph{robustness to the shape of the departure} rather than a headline power number.
Against a \emph{sparse} departure the dense-friendly tests trail badly (Fisher $0.578$, $\chi^2$
$0.629$) while the symmetric-root test is the strongest single competitor ($0.723$). Against a
\emph{dense} departure of the same energy the ranking inverts violently: the symmetric-root test
collapses to $0.376$---its whitening spreads a dense signal thinly over all coordinates, exactly
where a Simes-type statistic is weakest, whereas the Cholesky whitenings concentrate it---and
$\chi^2$ moves to the front of the invariant references, with the single-ordering Fisher test
close behind. No reference is safe in both regimes, whereas the pooled combiners hold
the top two places in each: $p$-merge $0.730$ and the e-value average $0.728$ against the sparse
departure, Bonferroni-over-orders $0.691$ and the e-value average $0.662$ against the dense one.
Which \emph{pair} leads changes with the shape---$p$-merge drops to $0.577$ on the dense departure
and Bonferroni to $0.713$ on the sparse---and the e-value average is the only test in the top two of
both, which is the finding \S\ref{sec:which} takes up. An analyst who must fix a test before seeing the departure is therefore
better served by pooling than by any single whitening.
The gain over the \emph{expected} single random ordering---the baseline
Proposition~\ref{prop:main} speaks to, in expected-log-evidence terms---grows with dependence
in this sweep, an empirical regularity the proposition motivates but does not assert: at known
$F$ it rises from $+0.017$ at $\rho=0.2$ to $+0.042$ and $+0.063$ at
$\rho=0.5$ and $0.8$ for the e-value average, and $+0.018$, $+0.041$, $+0.059$ for the p-merge.
The gains are modest---strong dependence enlarges the orbit, but most orderings of an
equicorrelated $\Sigma$ are reasonably favorable---and statistically unambiguous, with paired
$t\ge5$ in every known-$F$ cell. A cleaner decomposition compares the e-value average with
the calibrated mixture e-value of a \emph{single} ordering (Figures~\ref{fig:boot} and
\ref{fig:power}), which holds the base statistic fixed and isolates aggregation: at known $F$
that gap is $+0.012$, $+0.033$ and $+0.049$ at $\rho=0.2$, $0.5$ and $0.8$, $+0.025$ on the sparse
and $+0.095$ on the dense departure of Figure~\ref{fig:power}, and $+0.021$ to $+0.040$ across
$n=5$, $10$, $20$ (paired $t\ge4.7$ in every known-$F$ cell, $t\ge3.3$ under estimation).
\begin{figure}[t]\centering
\includegraphics[width=0.78\linewidth]{power_summary.png}
\caption{Calibrated power against departure shape at equal energy ($n=10$, $\rho=0.5$,
$\mathrm{ncp}=12$, $M=12$, known $F$); all tests are bootstrap-calibrated to size $\alpha$, so the
comparison is at matched level. The order-invariant references swap places between the two regimes:
the symmetric-root test is strong on the sparse departure and collapses on the dense one, Fisher
does the reverse, and $\chi^2$---flat across shapes by Proposition~\ref{prop:chi2}---moves up the
ranking on the dense departure only because the others move. The pooled combiners are at or near the top of both. The light-blue
bar is the calibrated single-ordering mixture e-value, the base statistic of the e-value average
without aggregation.}
\label{fig:power}
\end{figure}
\section{Screening many hypotheses}\label{sec:fdr}
\subsection{Calibrated pooling in the multiplicity layer}\label{sec:fdr-main}
When a collection of multivariate hypotheses $H_0^{(j)}:\mathbf{X}^{(j)}\sim F$, $j=1,\dots,N$, is screened
simultaneously, the single-hypothesis test of Section~\ref{sec:method} plugs into a standard
multiplicity layer---BH on per-hypothesis $p$-values, or e-BH on per-hypothesis e-values
$\bar E_M^{(j)}$. It introduces no new methodology, but it is the setting in which the shape
robustness of \S\ref{sec:headline} pays most clearly, because a screen rarely faces one kind of
departure: we let the non-null hypotheses fail heterogeneously, half on a single coordinate and half
on all of them, at common energy $\mathrm{ncp}=20$ ($n=20$, $\rho=0.5$, $200$ hypotheses,
$\pi_0=0.95$, $q=0.10$, $R=200$ replications). One design point matters: the bootstrap replicate
count should satisfy $B>N/q$, since a bootstrap $p$-value cannot fall below $1/(B+1)$ while BH's
\emph{smallest} threshold is $q/N$---the condition guarantees the whole BH ladder is reachable
(it is not necessary for making rejections: the $k$-th threshold is $qk/N$, so with $B=999$
ten non-null $p$-values of $0.001$ would all be rejected at the tenth rung, $0.005$; what a smaller
$B$ forfeits is the bottom of the ladder);
we use $B=2500$.
Two findings, both in Figure~\ref{fig:estF}. First, validity again hinges on calibration, and
again fails in both directions: with the na\"ive plug-in, BH realizes FDR $0.80$ for $\chi^2$,
$0.66$ for a single ordering and $0.61$ for the symmetric-root test against a target of $0.10$ at
$N_{\mathrm{tr}}=4n$, while e-BH on the conservative e-value average spends nothing at all (realized FDR
$0.000$ at known $F$). Calibrating each per-hypothesis statistic by Procedure~1 and applying BH to
the resulting bootstrap $p$-values brings every method to target---as an empirical matter:
the $N$ test vectors are independent, but their calibrated $p$-values share the estimated
$\hat\Sigma$ and the bootstrap pool, so BH's positive-dependence condition is not established here
and the realized-FDR figures are evidence of control under that shared calibration rather than a
theorem---with one exception worth naming:
$p$-merge realizes $0.137$ at $N_{\mathrm{tr}}=4n$, three standard errors above the $q\pi_0=0.095$ that BH
guarantees. The re-estimating bootstrap is an asymptotic device, and at four observations per
dimension it is not yet exact; a sparse screen, where few rejections make the realized proportion
sensitive to the tail of the calibrating distribution, is where that shows. It is gone by
$N_{\mathrm{tr}}=8n$. Second, on that level playing field the pooled combiners dominate: at known $F$,
Bonferroni-over-orders detects $0.642$ of the non-nulls and the e-value average $0.619$, against
$0.474$ for a single ordering, $0.353$ for the symmetric-root test and $0.349$ for $\chi^2$; the
same ordering holds under estimation ($0.490$ and $0.482$ against $0.323$, $0.320$ and $0.240$ at
$N_{\mathrm{tr}}=8n$). Bonferroni-over-orders leads here although it trails in the single-hypothesis
comparison of Figure~\ref{fig:boot}; \S\ref{sec:which} takes that up next, the short version being
that BH interrogates each statistic at $qk/N$ rather than at $\alpha$. The reason the pooled combiners
lead as a group is the shape robustness of
\S\ref{sec:headline}: the symmetric-root test recovers the sparse half of the non-nulls and misses
the dense half, $\chi^2$ is mediocre on both, and only the pooled tests are competitive on each. A
final advantage of e-BH, under a fully specified null where the $\bar E_M^{(j)}$ are e-values,
is that it controls FDR under \emph{arbitrary} dependence among the
hypotheses---relevant when the $\mathbf{X}^{(j)}$ share data or structure---whereas BH on $p$-values
requires independence or positive dependence.
\begin{figure}[t]\centering
\includegraphics[width=\linewidth]{estimated_F.png}
\caption{Multiplicity layer with \emph{heterogeneous} departures ($n=20$, $\rho=0.5$, $M=12$,
$\mathrm{ncp}=20$, $200$ hypotheses, $\pi_0=0.95$, non-nulls half sparse and half dense, target
$q=0.10$, $B=2500$, $R=200$). (A) With na\"ive per-hypothesis statistics the realized FDR is
uncontrolled for the min-type and order-invariant tests and near zero for the e-value average.
(B) After bootstrap calibration every method sits at target at every training size, apart from
$p$-merge at $N_{\mathrm{tr}}=4n$ (same scale as A). (C) Only then is a power comparison meaningful: the
pooled combiners lead at the two larger training sizes, because no shape-tuned competitor handles
both kinds of departure; at $N_{\mathrm{tr}}=4n$ every method is weak and the differences are within two
standard errors. Na\"ive power is deliberately not shown, being bought with excess FDR.}
\label{fig:estF}
\end{figure}
\subsection{Which combiner, and how many orderings}\label{sec:which}
\textbf{e-value averaging} is the combiner to reach for by default. It is never the weakest of the
three pooled rules in any setting we ran, and it is never far from the best. Figures~\ref{fig:boot}
and \ref{fig:estF} report three training sizes each, $N_{\mathrm{tr}}=4n$, $8n$ and known $F$; across those six
settings the e-value average's largest shortfall against the leading combiner is
$0.040$, against $0.076$ for Bonferroni-over-orders and $0.212$ for $p$-merge. That is the whole
argument for it. The alternatives win in one regime and lose in the other---$p$-merge leads the
single-hypothesis comparison at $\alpha=0.05$ (Figure~\ref{fig:boot}) and trails
Bonferroni-over-orders by $0.212$ in the screen (Figure~\ref{fig:estF}); Bonferroni leads the screen
and is last of the three at $\alpha=0.05$---whereas the e-value average is second in five of the six
settings and first in the remaining one, and in the threshold sweep reported below it is first or
second at every level from $0.05$ down to $0.0005$. A combiner chosen before one knows whether the departure will
be sparse or dense, or whether the threshold will be conventional or a screen's, should be the one
without a bad case.
Choosing it costs nothing in the form of the output. Under bootstrap calibration the e-value average
is treated exactly like any other statistic---Procedure~1 returns a calibrated $p$-value for it
directly, its accuracy being that of the bootstrap itself (Section~\ref{sec:calib-exp}) plus
the $1/(B+1)$ granularity that every rank-based $p$-value shares. Under a \emph{fully specified} null, where $\bar E_M$ is a genuine e-value
(Lemma~\ref{lem:evalue} and the caveat after it), no calibration is needed at all: Markov's
inequality gives the valid $p$-value $1/\bar E_M$ at once. The converse does not hold: a $p$-value cannot be turned into
an e-value without a calibrator and a loss. So under a fully specified null the e-value
average yields both currencies, feeds e-BH directly for FDR control under arbitrary dependence
among the (valid) e-values (Section~\ref{sec:fdr}), and composes across independent datasets by
multiplication (or sequentially, when each factor is an e-value conditionally on the past); under
an estimated-population null it is a calibrated test statistic like the others, and the e-value
reading is unavailable.
Why the two specialists trade places is worth setting out, since it is what makes the default
safe. A ranking of combiners stated at $\alpha=0.05$ does not survive a move to the far smaller
thresholds a screen imposes. Benjamini--Hochberg over $N$ hypotheses reports a rejection
only when some $p$-value reaches $qk/N$, which for the screen of \S\ref{sec:fdr-main} means
$0.0005$ to $0.005$ when the rejection set consists of the non-nulls (false positives raise
$k$ and loosen the rung): one to two orders of magnitude below the level of a single test. The
combiners are not equally placed to supply such values. Bonferroni-over-orders, $M\min_m P^\sigma_m$,
keys on the single most favourable ordering, so its calibrated $p$-value has a heavier left tail
than the averaging rules and a worse typical value; $p$-merge and the e-value average use all $M$
views, which pays at a moderate threshold and costs at an extreme one. Holding the design of
Figure~\ref{fig:boot} completely fixed---same covariance, dimension, alternative (two shifted
coordinates, $\mathrm{ncp}=12$) and $M$, known $F$---and sweeping only $\alpha$ makes the mechanism
plain (the sweep is a separate, larger Monte-Carlo run with one fixed set of orderings and
$2\times10^5$ draws, so its $\alpha=0.05$ values differ a little from Figure~\ref{fig:boot}'s):
$p$-merge leads Bonferroni-over-orders at $\alpha=0.05$ ($0.510$ against $0.487$) and trails
it at $\alpha=0.0005$ ($0.045$ against $0.066$), having lost over nine-tenths of its power to
Bonferroni's six-sevenths, with nothing but the threshold changing. Against the dense half of the
Figure~\ref{fig:estF} screen (all coordinates shifted, $\mathrm{ncp}=20$) the same sweep is
starker: $p$-merge falls from $0.757$ to $0.094$ while Bonferroni falls only from $0.892$ to
$0.393$, their gap widening from $0.135$ to $0.299$. This is why Bonferroni-over-orders is the weakest of the three pooled
combiners in Figure~\ref{fig:boot} and first in Figure~\ref{fig:estF}. It is also why the e-value average is the
safe default: alone among the three it is competitive at both ends of the threshold range, so it
does not require the user to know in advance which end they are at. If that is known, the
specialists are slightly better---$p$-merge for a single hypothesis at a conventional level,
Bonferroni-over-orders in a sparse screen.
\subsection{Comparison with conformal novelty detection}\label{sec:conformal}
A reader who screens many multivariate observations for anomalies will ask how this compares with
conformal novelty detection, which has developed rapidly in exactly this space
\cite{bates2023,marandon2024}. The honest answer begins with a distinction, because the two
literatures test different nulls.
Our null is \emph{conformity to a specified $F$}: $H_0^{(j)}:\mathbf{X}^{(j)}\sim F$, with $F$ supplied by
the modeller and, when unknown, estimated from a training sample and accounted for by
Procedure~1. The conformal null is \emph{exchangeability with a reference sample}: given inliers
$Y_1,\dots,Y_{N_{\mathrm{tr}}}$ and test points $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(N)}$, one computes a nonconformity score
$S$, holds out $\ell$ of the reference points for calibration, and forms
\begin{equation}\label{eq:confp}
p_j \;=\; \frac{1+\#\{i \in \text{cal}: S_i \ge S(\mathbf{X}^{(j)})\}}{\ell+1},
\end{equation}
which is super-uniform under exchangeability alone. Under the independence and scoring
assumptions of Bates et al.\ \cite{bates2023}---test inliers and reference data jointly
independent, with conditions on the score distribution, which our simulation design
satisfies---these $p$-values are PRDS, so BH controls FDR; Marandon et al.\ \cite{marandon2024} observe that a
prespecified score wastes power and propose \emph{AdaDetect}, which splits the reference sample,
trains a probabilistic classifier to separate the fit part from the calibration part pooled with the
test sample, and uses the predicted probability of belonging to the latter as $S$. Because that
score is invariant to permutations within calibration-plus-test, \eqref{eq:confp} remains valid, so
the score may be \emph{learned} without any distributional assumption.
The trade is therefore explicit. By conformal methods we mean, as in the works cited,
procedures calibrated against observed inliers alone (a rank $p$-value of the same form computed
against draws \emph{simulated} from a specified model is the Monte-Carlo test of Procedure~1, and
belongs to the parametric route). They buy validity without distributional assumptions, at the
price of needing a reference sample at test time, and they answer a different question from ours:
whether a new observation resembles the reference data, not whether a given model is correct. If a
model misdescribes the reference data and the new data in the same way, the new data remain
exchangeable with the reference sample and a conformal test has nothing to detect, whereas the
parametric route asks directly the question a risk manager asks---is \emph{this} model right?
Because the two routes test different nulls, their power can be compared only where the nulls
coincide, that is, when the specified model is correct, and there one should expect access to the
model to help. The comparison below therefore asks not \emph{whether} the model helps but
\emph{how much}: how large the advantage is once the model must be estimated from the same
reference sample the conformal methods receive, in the sparse regime in which screens operate, and
whether a learned conformal score closes the gap. We do not study misspecification, where the two
routes answer different questions and power is not the relevant comparison.
\paragraph{Design} We reuse the screening problem above ($n=20$, $N=200$ hypotheses,
$\mathrm{ncp}=20$, half sparse and half dense non-nulls, $\rho=0.5$, $q=0.10$, $M=12$, $B=2500$),
supply both routes with the same reference sample, and compare five conformal scores: the Mahalanobis
form $\mathbf{y}^\top\hat\Sigma_{\mathrm{fit}}^{-1}\mathbf{y}$, the nearest-neighbour distance, the
negative log of a Gaussian kernel density estimate, and AdaDetect with a random forest and with
penalized logistic regression on quadratic features.
\paragraph{Held-out calibration size} The resolution consideration behind $B>N/q$ on the
bootstrap binds far harder on the conformal calibration set once novelties are rare. A conformal
$p$-value is a multiple of $1/(\ell+1)$, so reaching the $m_1$-th BH threshold $qm_1/N$
requires $\ell+1\ge N/(qm_1)$, with $m_1$ a lower bound on the number of rejections
(Remark~2.2 of \cite{marandon2024} states the condition in the strict form $\ell>N/(qm_1)$); with
the ten non-nulls of \S\ref{sec:fdr-main} this is $\ell\ge199$, strictly $\ell>200$. We start the
reference-sample grid at $N_{\mathrm{tr}}=400$ and use a balanced split $\ell=N_{\mathrm{tr}}/2$, which clears the
condition at $N_{\mathrm{tr}}=1000$ and $2000$; at $N_{\mathrm{tr}}=400$ the split $\ell=200$ sits exactly on the
resolution boundary, where any nonempty BH rejection set contains at least ten hypotheses (if
no true null is rejected, all ten non-nulls must be rejected together), so the $N_{\mathrm{tr}}=400$
conformal column of Table~\ref{tab:conformal} is resolution-limited by construction and should be
read as such. Where a method's authors recommend a value we use theirs: for AdaDetect
the same remark recommends $\ell=N$, which we adopt wherever its strict condition admits it
and replace by the balanced split where it does not (the balanced split is a design choice
fixed in advance; the resolution condition only bounds $\ell$ from below). All figures are averages over $R=200$
replications, with Monte-Carlo standard errors given in the captions.
\begin{table}[t]\centering\small\setlength{\tabcolsep}{4pt}
\caption{Comparison with conformal novelty detection in the sparse screening regime
($n=20$, $N=200$ hypotheses, $\pi_0=0.95$, so ten novelties, $\mathrm{ncp}=20$, $q=0.10$,
$R=200$). Entries are calibrated power; every method controls FDR, the largest realized value being
$0.107$ against a target of $0.10$, within one Monte-Carlo standard error. Standard errors are at
most $0.012$ on the realized FDR and $0.018$ on the powers. Every tuning constant is fixed in
advance and none is selected on the results: the distance to the nearest neighbour as the
nonparametric score, and a balanced calibration split $\ell=N_{\mathrm{tr}}/2$, a design choice fixed
in advance (the resolution condition only bounds $\ell$ from below). The data are generated from
the specified Gaussian model, the setting in which the two routes test the same null.}
\label{tab:conformal}
\begin{tabular}{lccc}
\toprule
& $N_{\mathrm{tr}}=400$ & $N_{\mathrm{tr}}=1000$ & $N_{\mathrm{tr}}=2000$\\
\midrule
Pooled, Bonferroni-over-orders (this paper) & \textbf{0.629} & \textbf{0.644} & \textbf{0.663}\\
Pooled, e-value average (this paper) & 0.605 & 0.613 & 0.630\\
Pooled, $p$-merge (this paper) & 0.459 & 0.445 & 0.465\\
$\chi^2$ energy test & 0.336 & 0.338 & 0.369\\
\addlinespace
Conformal, nearest-neighbour distance & 0.073 & 0.405 & 0.520\\
Conformal, Mahalanobis & 0.038 & 0.183 & 0.336\\
Conformal, kernel density & 0.037 & 0.141 & 0.265\\
AdaDetect, random forest & 0.005 & 0.116 & 0.235\\
AdaDetect, penalized logistic & 0.003 & 0.026 & 0.063\\
\bottomrule
\end{tabular}
\end{table}
\paragraph{What the specified model is worth when it is right}
Table~\ref{tab:conformal} answers the quantitative question. With the model estimated from the same
reference sample, the pooled tests lead across the whole range of reference-sample sizes: at
$N_{\mathrm{tr}}=1000$---five times the size of the screen itself---Bonferroni-over-orders detects $0.644$ of
the ten novelties against $0.405$ for the best conformal method, and at $N_{\mathrm{tr}}=2000$, ten times the
screen, $0.663$ against $0.520$, gaps of $10.6$ and $8.4$ Monte-Carlo standard errors. The rows
compare complete procedures, each a score together with a calibration, and the specified model
enters the pooled tests through both: through the Rosenblatt score, and through a calibration that
simulates from the fitted model rather than holding out inliers. One pair of rows uses the same
quadratic-score form. The conformal Mahalanobis score and the $\chi^2$ energy test are the same quadratic form (with
the covariance estimated from the fitting part of the reference sample in the conformal case),
calibrated respectively against held-out inliers and by Procedure~1, so the gap between
them---$0.038$ against $0.336$ at $N_{\mathrm{tr}}=400$, $0.183$ against $0.338$ at $1000$, $0.336$ against
$0.369$ at $2000$---illustrates the combined benefit of model-based calibration and of retaining
the whole reference sample for fitting (the conformal version fits on half of it and calibrates
on the other half, with the resolution of a held-out set): decisive when reference data are
scarce and small once they are plentiful. The gaps narrow over the reference-sample sizes
examined; performance beyond $N_{\mathrm{tr}}=2000$ was not assessed, and different scores need not converge
to equal power (the conformal Mahalanobis score approaches an energy test, which at known $F$
remains well below the pooled combiners, Figure~\ref{fig:estF}). The reference samples of the
larger sizes are in any case larger than the approximately $250$ observations used in our FX
application, which calibrates on a single year of returns and tests the next.
It is worth also comparing the conformal methods with one another. AdaDetect exists because a
prespecified nonconformity score wastes power, and it replaces one with a score learned from the
data; here that reverses. At $N_{\mathrm{tr}}=2000$ its random forest reaches $0.235$ and its penalized
logistic version $0.063$, while the three prespecified scores reach $0.520$, $0.336$ and $0.265$:
the two learned scores are the weakest of the five. This is not a criticism of AdaDetect, which is
designed for a regime it is not being given here. A classifier trained to separate inliers from a
test sample can only learn what the novelties in that test sample reveal, and ten novelties diluted
among a hundred and ninety inliers reveal very little.
\section{Application: validating a Gaussian foreign-exchange risk model}\label{sec:fx}
We close the empirical work with a real risk-management application in which the coordinate ordering
is genuinely arbitrary: daily foreign-exchange (FX) returns. The FX market is where currencies are
traded; an \emph{exchange rate} is the price of one currency in another (e.g.\ dollars per euro),
and a currency's daily \emph{return} is the relative change in that price from one trading day to
the next. A bank or fund holding positions in several currencies summarizes a day of market movement
as one vector $\mathbf{X}_t\in\mathbb{R}^n$---the day's returns of the $n$ currencies it is exposed to---and its
\emph{risk model} is a probability distribution for that vector, used to set capital buffers and
loss limits. The simplest and still most widely used such model is multivariate Gaussian with mean
and covariance estimated from recent history. \emph{Validating} the model means asking, day after
day, exactly this paper's question: is today's observed vector consistent with the model's
distribution? In the language of the forecasting literature, this is multivariate density-forecast
evaluation \cite{dgt1998,dovern2020}: the fitted $\mathcal N(\hat\boldsymbol{\mu},\hat\Sigma)$ is a one-day-ahead
density forecast, held fixed within each year, and each day's realization is scored against it. And the coordinates of that vector are currencies: no convention orders them, so any
fixed choice (alphabetical tickers, portfolio weights, legacy code) is exactly the hidden analyst
degree of freedom this paper is about.
\paragraph{Setup}
Data are daily spot exchange rates for nine major currencies against the U.S.\ dollar (euro,
Japanese yen, British pound, Canadian dollar, Swiss franc, Australian and New Zealand dollars,
Swedish and Norwegian kronor; all expressed as USD per foreign unit), published by the Federal
Reserve (H.10 release, via FRED), 2015--2025: $2{,}747$ days of log-returns with $n=9$. The model to
validate is refitted in the standard \emph{walk-forward} way: at the start of each calendar year,
$\hat\boldsymbol{\mu}$ and $\hat\Sigma$ are estimated from the preceding year's $\approx250$ trading days (so
$N_{\mathrm{tr}}\approx28n$---moderate, realistic estimation), and every day of the new year is tested against
$\mathcal N(\hat\boldsymbol{\mu},\hat\Sigma)$ out of sample. The null tested each day deserves a precise
statement, because two readings are possible and they license different guarantees. We test the
\emph{issued forecast}, stated conditionally: $H_0^{(t)}:\mathbf{X}_t\mid\mathcal F_{t-1}\sim\mathcal
N(\hat\boldsymbol{\mu}_t,\hat\Sigma_t)$, with $(\hat\boldsymbol{\mu}_t,\hat\Sigma_t)$ the numbers the model actually produces
from information available before day $t$---the density-forecast-evaluation
null of \cite{dgt1998}. Under this null the whitened scores are \emph{exactly} standard normal,
Lemma~\ref{lem:evalue} applies as stated, and the raw pooled e-values are valid without
adjustment; rejections mean ``the forecast, as issued, is wrong about this day.'' This differs
from the estimated-population null of \S\ref{sec:calib}, under which the plug-in scores are not
e-values (see the caveat after Lemma~\ref{lem:evalue}) and only calibrated $p$-values are
available. The two nulls call for different calibrations, and we run both. Under the
issued-forecast null the whitened scores are exactly standard normal, so each statistic's null
distribution is obtained by simulating from the issued forecast itself and whitening by the same
fixed factors as the data (Procedure~1 with step~1 replaced by exact draws from
$\mathcal N(\hat\boldsymbol{\mu}_t,\hat\Sigma_t)$; no re-estimation). Under the estimated-population null we use
the re-estimating bootstrap exactly as in Section~\ref{sec:exp}, with mean and covariance
re-estimated inside every replicate. The per-year rejection rates below are reported under both;
the e-value ranking and the e-BH screen are statements about the issued forecasts only. For each
test year $2016$--$2025$ we run, at calibrated level $\alpha=0.05$ ($B=999$, $M=12$ fresh
orderings per day): the single-random-ordering Simes test, the three pooled combiners, and the two
order-invariant references. Every day is one multivariate hypothesis---the paper's basic setting---and the ten years
form a $2{,}497$-day screening problem for e-BH. If the Gaussian model were exactly right, every
method would flag $\approx5\%$ of days each year; systematic exceedance of $5\%$ is evidence of model
failure, and \emph{which} days carry the evidence is exactly what a risk manager wants to know.
\paragraph{The test behaves as designed}
In calm years the calibrated rejection rates sit near the nominal $5\%$ (2017: $0.06$; 2021: $0.03$;
2023: $0.06$; 2024: $0.04$ for the e-value average under the estimated-population calibration;
$0.06$, $0.03$, $0.08$, $0.05$ under the issued-forecast calibration, which absorbs no
estimation noise and so runs at or a little above the other in every year). One caution in reading this: with the truth
unknown, a near-nominal rejection rate is \emph{consistent with} adequacy rather than proof of
it---a misspecified family can produce the same rate for this statistic. What the calm years do
establish is that the calibrated procedure behaves sensibly out of sample: last year's
covariance is not visibly contradicted by this year's ordinary days. In stress years the model is flagged
massively: $30\%$ of days in 2020 (the COVID pandemic), $49\%$ in 2022 (the U.S.\ Federal Reserve's
aggressive rate hikes and the resulting dollar surge), $16\%$ in 2025 (the April tariff shock),
$14\%$ in 2016 (the United Kingdom's Brexit referendum); $32\%$, $55\%$, $17\%$ and $14\%$
under the issued-forecast calibration. The pooled e-values also \emph{rank} the
history interpretably: an e-value is the factor by which the day's evidence multiplies a bet
against the model, so $\bar E_M=20$ is borderline ($=1/\alpha$) while
$\bar E_M=10^{13}$ is overwhelming. The largest values in ten years, in order, are the Brexit
referendum result day 2016-06-24 ($\bar E_M\approx10^{20}$: sterling fell $8\%$ against a model
fitted on placid 2015), the Bank of Japan's surprise loosening of its yield-curve control on
2022-12-20 ($\approx5\times10^{17}$: a yen-specific jump), the COVID ``dash-for-dollars'' days
2020-03-18/20/23 ($5\times10^{12}$--$10^{14}$), the day the Swiss National Bank unexpectedly raised rates
alongside a large Fed hike, 2022-06-16, and the 2025-04-04 tariff announcement
(Figure~\ref{fig:fx}A). None of these dates was supplied to the procedure; they emerge from the
statistics alone. Finally, applying na\"ive e-BH at $q=0.10$ across all $2{,}497$
days---consecutive days are serially dependent, which e-BH tolerates by construction, no
adjustment needed---flags $78$ days, essentially the union of these episodes: a defensible list
of ``days on which the issued forecasts demonstrably failed,'' with FDR control that is valid
under the forecast nulls (under the estimated-population reading no such guarantee is claimed;
see the caveat after Lemma~\ref{lem:evalue}).
\paragraph{What pooling buys here---and honestly, what it does not}
The day-level paired gain of the e-value average over a single random ordering is $+0.004\pm0.004$
under both calibrations: statistically indistinguishable from zero (the days are serially dependent, so the standard
error is computed with a Newey--West adjustment; the conclusion of a null gain is unchanged
either way). This is the operating map of \S\ref{sec:headline} speaking,
not a failure: real FX stress days are \emph{far} from the borderline-power regime---when the
departure is a $10$-sigma event every ordering rejects, and when the day is ordinary none does. The
practical value of pooling in this application is different and threefold: it removes the arbitrary
ordering from the \emph{definition} of the procedure (the reported statistic no longer depends
on a ticker convention: at finite $M$ it is a Monte-Carlo approximation to the fully
order-symmetric statistic, whose distribution is the same whichever labelling the data arrive in;
the realized number is reproducible for a given input labelling and seed, and exact invariance
holds in the $M\to\infty$ average), it supplies valid forecast-null
e-values for dependence-robust day screening via e-BH, and its observed rejection rates are
indistinguishable from a single ordering's---which, without ground-truth labels for the days, is a
statement about how often the methods flag, not a demonstration of equal power. The
episode-by-episode pattern of \S\ref{sec:headline} also appears here, with the same caveat that
these are observed flag rates rather than measured power: the
best-performing reference \emph{changes across episodes} (the symmetric-root test leads in 2020,
$\chi^2$ leads in 2022; Figure~\ref{fig:fx}B), while the pooled test tracks within
$0.01$--$0.04$ of whichever is winning without knowing the episode's shape in advance.
\begin{figure}[t]\centering
\includegraphics[width=\linewidth]{fx_application.png}
\caption{FX application (9 currencies vs.\ USD, 2016--2025, yearly walk-forward refit). (A) Pooled
e-value per day; red dots are the $78$ days flagged by e-BH at $q=0.10$; annotated spikes are the
Brexit referendum, the COVID dollar scramble, the SNB/Fed and BoJ surprises of 2022, and the 2025
tariff shock. (B) Calibrated per-day rejection rates by year (bars: estimated-population null,
re-estimating bootstrap; diamonds: e-value average under the issued-forecast null, exact
simulation from the issued forecast): near-nominal in calm years, large in
stress years; the leading order-invariant reference changes across episodes while the pooled test
stays near the top.}
\label{fig:fx}
\end{figure}
\section{Discussion}\label{sec:disc}
The paper followed a single question---whether one must gamble on a coordinate ordering---to a
three-part answer: reordering only reshuffles a fixed signal, so pooling \emph{can} help; the base
test must look in both tails, or the gamble is the least of one's problems; and calibration is what
converts the latent advantage into a usable one. We close by stating what this buys, where it stops,
and how it connects.
\paragraph{What the evidence supports}
The reason to pool is not that any one ordering is bad but that one cannot know in advance which
ordering, or which fixed whitening, suits the departure at hand. That is why the gain over a single
random ordering is real yet modest, while the gain over a \emph{committed} alternative---an
order-invariant test tuned, implicitly, to one shape of departure---is large and, in a screen whose
hypotheses deviate from the null in different ways, decisive. The theory says the same thing in
another language: reordering conserves the signal energy and only redistributes it
(Lemma~\ref{lem:energy}), so pooling cannot manufacture evidence, it can only stop one from
betting everything on a single view of it (Proposition~\ref{prop:main}). The foreign-exchange
application shows the machinery surviving contact with data it was not designed around: rejection
rates stay near nominal out of sample in calm years (consistent with adequacy, not proof of it),
the largest e-values are the episodes a risk manager would name
unprompted, and e-BH screens a decade of serially dependent days without adjustment.
\paragraph{An honest operating map}
The boost over a single ordering is small at near-independence (small orbit), and it is masked
entirely by na\"ive thresholds---calibration is not optional. The symmetric-root test is a
competitive order-invariant alternative for sparse departures, and it collapses for dense ones
(Figure~\ref{fig:power}); the case for pooling is that it
does not need to know the signal's shape in advance. One methodological point deserves emphasis:
pooling must be evaluated against the \emph{expected} power of a random ordering, not against one
fixed ordering, whose value depends on which coordinates happen to precede the departure and can
overstate the gain several-fold.
Guidance on the choice of combiner is collected in \S\ref{sec:which} and not repeated here; the
short version is that the e-value average is the default, being the only one of the three pooled
rules without a regime in which it does badly, and that it costs nothing to adopt since bootstrap
calibration returns a $p$-value for it directly. The one prohibition worth restating is that the
arithmetic mean of $p$-values must not be used without calibration. To that map
\S\ref{sec:conformal} adds a boundary against conformal novelty detection. A learned conformal score
must learn the alternative from the novelties present in the test sample, and where those are rare
it has little to learn from; since screens are usually sparse, that leaves a good deal of ground to
a procedure that brings its own null, when that null is right.
\paragraph{Limitations and scope}
We deliberately restrict to the Gaussian null. Behavior under heavier-tailed and non-parametric nulls
is more intricate---the orbit geometry, calibration quality, and even the sign of the gain can change
with the tail behavior of $F$---and a full treatment is left to future work. It bears saying plainly
that the error control established here is contingent on that family being the right one: like every
test of a specified null, and unlike the conformal procedures of \S\ref{sec:conformal}, our
guarantee is a statement about departures from $F$ and not about departures from whatever generated
the data. Our calibration
machinery is general (Procedure~1 needs only the ability to simulate from $\hat F$ and recompute the
transform), so extending to copula and non-parametric nulls, where the orbit may be even richer, is a
natural next step. Our power study uses location (mean-shift) alternatives at fixed energy under
equicorrelated covariances (with a hub design used only for the nominal-threshold sweep of
Remark~\ref{rem:dome}); broader covariance designs, scale/dependence/shape alternatives,
estimated nuisance parameters, and skewed nulls remain open. The location-style Simes base statistic
is also not omnibus: for instance, an \emph{under}-dispersed coordinate ($\sigma=0.5$ against an
assumed $1$) leaves the power of the one- and two-sided versions alike at or below the nominal
level (about $9\times10^{-5}$ for the two-sided test at $\alpha=0.05$ in the one-dimensional
case), so scale-aware
base statistics are a natural addition to the toolkit.
\paragraph{Relation to prior work}
The closest predecessor is \cite{dovern2020}, who construct order-invariant tests from
conditional PITs across orderings and address estimation uncertainty; our full-order average
$\mathbb{E}_\sigma[E^\sigma\mid\mathbf{X}]$ is of course itself a symmetrization, so the distinction is not
aggregation versus symmetrization but \emph{what} is aggregated and with what guarantee. What
this paper adds over that construction is: aggregation at the level of \emph{evidence}
(e-values and $p$-values) with finite-$M$ validity under the arbitrary dependence that
recomputation creates (Lemma~\ref{lem:merge})---a randomized $M$-ordering test is valid as a
test, not only in a symmetrized limit; the two-sided base statistic, which removes the
directional blind spot of one-sided scoring when the sign of the departure is unknown
(Section~\ref{sec:misspec}); the re-estimating calibration that makes
the aggregate usable and powerful under estimation; and the multiplicity layer. Specification testing via the
Rosenblatt transform is classical for copulas \cite{genest2009} and for density-forecast evaluation
\cite{dgt1998}. The e-value averaging used to combine orderings, and the e-BH multiplicity layer of
Section~\ref{sec:fdr}, connect the test to recent e-value theory
\cite{vovkwang2021,vovkwang2020,wangramdas2022,grunwald2024}. Contemporary
PIT-of-order-statistics methods combine the transform with Cauchy combination \cite{covington2025}; our
contribution is orthogonal---the ordering-aggregation geometry, its calibration, and the power
characterization. The closest recent alternative for the screening problem of
Section~\ref{sec:fdr} is conformal novelty detection \cite{bates2023} and its adaptive extension
\cite{marandon2024}, which replace a specified null by exchangeability with a reference sample;
\S\ref{sec:conformal} sets out the different hypotheses the two answer and measures what the
specified model is worth when it is correct.
\paragraph{What this paper is not}
It is worth forestalling a natural misreading: the recipe here is not ``whenever several tests
exist, apply all of them and merge the evidence.'' Validity is not what separates the two:
Jensen's inequality and the merging lemmas apply to any collection of valid e-values. What the
generic recipe lacks is the structure that makes the average worth taking here, and this paper
itself furnishes the counterexamples. The $M$ views pooled here are not
different methods but the \emph{same} statistic under a symmetry of the null, and that structure is
what the results use: each view carries exactly the same total signal (Lemma~\ref{lem:energy}), so
averaging is over equally informed presentations, never over procedures of unequal merit; and the
choice being averaged over is \emph{forced}---the transform cannot be computed without an ordering---so
the alternative to pooling is not ``pick the best method'' but ``commit to an arbitrary one.''
Where those conditions fail, so does the slogan. Pooling adds nothing to the Fisher base statistic
and pooled Fisher can fall below a single ordering (Remark~\ref{rem:fisher}); against the
exchangeable dense alternative every ordering carries the identical \emph{mean} profile
($P_\sigma\boldsymbol{\mu}=\boldsymbol{\mu}$, with an exchangeable $\Sigma$), so there is no ordering lottery for
randomization to fix and any pooling gain there operates through the noise rather than the
signal; and without calibration the multiplicity cost of merging cancels the gain
almost exactly (Remark~\ref{rem:dome}). What the paper actually argues is narrower: when a
\emph{nuisance choice} indexes equally informed views of one observation, do not commit---average,
and calibrate so that the averaging is free.
\paragraph{Acknowledgments}
Claude Opus 4.7 and Claude Fable 5 (Anthropic) were used in the preparation of this manuscript; the
author thanks Javad Alipanah for providing access to them. The author retains full responsibility
for the mathematical content and conclusions of the paper.
\paragraph{Declaration of generative AI in scientific writing}
During the preparation of this work the author used Claude (Anthropic) in the design and
verification of the experiments, in drafting, and in checking the mathematical arguments. After
using this tool, the author reviewed and edited the content as needed and takes full
responsibility for the content of the published article.