Despite the overwhelming empirical evidence, in most real-life datasets, of non-Gaussian and non-symmetric (non-elliptical) distributions, multivariate time series analysis remains very deeply marked by explicit or implicit Gaussian assumptions: correlogram-based and spectral methods, pseudo-Gaussian tests, Gaussian quasi-likelihoods, and their possible robustifications, are ubiquitous in methodological developments as well as in daily practice.
The theoretical justification for this prevalence of Gaussian methods looks quite solid at first sight and is seldom questioned: Gaussian quasi-likelihood methods are asymptotically valid (viz., pseudo-Gaussian tests have correct asymptotic nominal size and Gaussian quasimaximum-likelihood estimators (QMLEs) are root-n consistent) under a broad range F of non-Gaussian innovation densities f (typically, under finite fourth-order moments).
A closer look at this asymptotic validity argument, however, reveals a fundamental weakness: while it holds, pointwise, for any given f ∈ F, it fails to hold uniformly over F.
Consider, for instance, the sequence φ (n) of Gaussian level-α tests rejecting a null hypothesis of the form θ = θ 0 whenever some test statistic S (n) exceeds the standard normal quantile z α of order (1 − α). Denoting by P (n) θ,f the distribution of the observation under parameter value θ and innovation density f ∈ F, assume that φ (n) is a valid pseudo-Gaussian test at asymptotic level α, that is, satisfies lim n→∞ P (n) θ 0 ,f S (n) > z α ≤ α (1.1) for any f ∈ F. This, however, is not sufficient if f remains unspecified. The null, then, is the semiparametric hypothesis P (n) θ,f | θ = θ 0 , f ∈ F and, in order for φ (n) to qualify as an asymptotically valid level-α test, the stronger condition
which, in general, does not follow from (1.1), is required.
This lack of uniformity in the asymptotics of pseudo-Gaussian procedures is not just a mathematical detail: in the context of VAR estimation, Hallin et al. (2019) have shown that, depending on the actual innovation density f , the finite-sample performance of a valid QMLE, for given sample size n, can be extremely far from its asymptotic performance. That discrepancy between finite-sample and asymptotic performance is particularly significant (see the motivating example in the introduction of the same reference) in the presence of non-elliptical innovation densities such as mixture or skewed densities-a type of innovation distribution which is likely to occur in the case of omitted variables or covariates-and only gets worse as the dimension of the observation space increases. The same phenomenon is bound to take place in testing problems.
A remedy to this lack of uniformity of pseudo-Gaussian asymptotics is distributionfreeness. If indeed the test statistic S (n) is distribution-free under P (n) θ,f | θ ∈ Θ 0 , f ∈ F , then (1.1) and (1.2) are equivalent. Rank-based inference then naturally enters the picture: it follows indeed from Basu's classical theorems on the relation between minimal sufficiency and distribution-freeness (see Appendix E in Hallin et al. (2020a) or Chapters 5 and 10 in Lehmann and Romano (2005) ) that, in the univariate case, the σ-field generated by residual ranks is essentially 1 maximal ancillary-intuitively, "maximal distribution-free".
Moreover, it follows from Hallin and Werker (2003) that the semiparametric efficiency bounds at given densityf ∈ F-which, for ARMA models, coincide with the parametric onescan be reached by rank-based tests: normal-score (van der Waerden) rank tests in ARMA models, in particular, are achieving, under Gaussian f , the same asymptotic performance as the (pseudo-)Gaussian ones while doing strictly and uniformly better under non-Gaussian f (Hallin, 1994) . Pseudo-Gaussian tests, thus, are not admissible and uniformly dominated by their normal-score rank-based counterparts.
All this, which strongly suggests abandoning pseudo-Gaussian methods in favor of rankbased ones, unfortunately, until recently was limited to univariate ARMA models due to the lack of an adequate multivariate extension of the univariate concept of ranks. Many proposals have been made in the literature-see Hallin et al. (2020a) for a commented bibliography. None of them, however, is enjoying distribution-freeness-let alone maximal ancillarity-except for the so-called Mahalanobis ranks and signs Paindaveine, 2002, 2004a) under the quite restrictive assumption, however, of elliptical symmetry. This lack of a satisfactory concept of multivariate ranks has been dealt with recently with the introduction of the measure-transportation-based Monge-Kantorovich (Chernozhukov et al., 2017) and center-outward (Hallin, 2017) ranks and signs, the maximal ancillarity of which is established in Hallin et al. (2020a) . Building on that concept, Hallin et al. (2020b) are proposing distribution-free center-outward rank-based tests for multiple-output regression and MANOVA models, Ghosal and Sen (2019) for goodness-of-fit, Deb and Sen (2019) and Shi et al. (2020) for independence between vectors, while Hallin et al. (2019) are constructing center-outward R-estimators for VARMA models. The finite-sample performance of the new rank-based methods in all cases appears to be quite remarkable. In this paper, we similarly study a class of center-outward rank-based testing procedures for VAR models and illustrate their role in the problem of sequential VAR order selection (see Pötscher 1983) in the presence of unspecified innovation density.
The paper is organized as follows. Section 2 introduces notation and recalls some methodological facts about local asymptotic normality for VARs (Section 2.2) and the measure-transportation approach to center-outward ranks and signs (Section 2.3). Section 3 contains the main theoretical result of this paper: a Hájek representation theorem for serial statistics based on center-outward ranks and signs. This key result, of independent interest, extends to general scores the spherical-score result of Hallin et al. (2019) , as well as the univariate results of Hallin et al. (1985) and Hallin and Vermandele (1996) . Based on LAN and the Hájek representation of Section 3, Section 4 introduces the notion of rank-based central sequence from which the rank-based test statistics are constructed in Section 5.
Focusing on two particular problems, we develop, in Section 5.1, a rank-based extension of the classical Durbin-Watson procedure to multiple-output and possibly non-linear regression and, in Section 5.2 a sequential rank-based procedure for VAR order identification. Section 6 derives the Gaussian counterparts of the tests developed in Section 5 while Section 7 presents a Monte Carlo investigation of the finite-sample performance of the rank-based methods we are proposing (Sections 7.1 and 7.2) and concludes with a real-data application (Section 7.3).
Our approach is combining two basic methodological ingredients:
(a) the local asymptotic normality of VAR models, as established in Garel and Hallin (1995) and Hallin and Paindaveine (2004a) ; the explicit statement (Section 2.2 and Proposition 2.1) of this property (with closed-form expressions 2 for central sequences and information matrices) requires some algebraic preparation, which is provided Appendix A and can be skipped at first reading;
(b) the recently proposed (Chernozhukov et al. (2017) ; Hallin (2017) ; Hallin et al. (2020a) ) measure-transportation-based concepts of multivariate center-outward ranks and signs: see Section 2.3 for precise definitions .
In this section, we first settle the notation (Section 2.1) before turning (Sections 2.2 and 2.3) to (a) and (b).
innovation process with mean 0 and probability density f , and θ := ((vecA 1 ) , . . . , (vecA p 1 ) ) denotes the (p 1 d 2 )-dimensional vector of autoregressive parameters. For θ of the form
of order p 0 < p 1 . For p 0 = 0, (2.2), with θ 0 = (0, . . . , 0) yields the white noise model X t = t .
Denote by Θ the set of parameter values θ ∈ R p 1 d 2 such that all solutions of the determinantal equation det A θ (z) = 0, z ∈ C lie outside the unit ball (the standard VAR stationarity condition) and by Θ 0 the set of parameter values θ 0 ∈ Θ of the form (2.2) with |A p 0 | = 0. The notation A 0 (L) will be used for a null VAR(p 0 ) operator of the form (2.3) with unspecified θ 0 ∈ Θ 0 , the notation A 1 (L) for an alternative VAR(p) operator A θ 1 (L) of order p, that is, with unspecified θ 1 ∈ Θ and A p = 0 for some p 0 < p ≤ p 1 . Denoting by X (n) := {X (n) 1 , . . . , X (n) n } (superscript (n) omitted whenever possible) an observed triangular array of solutions of (2.1), our objective is to construct rank-based tests of (a) the null
under which X (n) was generated by a VAR operator of the form A θ 0 (L) (specified θ 0 ∈ Θ 0 , unspecified innovation density) and (b) the null hypothesis H (n) 0 under which X (n) was generated by a VAR operator of the form A 0 (L) (unspecified θ 0 ∈ Θ 0 , unspecified innovation density) against the alternative of a VAR operator of the form A 1 (L) (unspecified θ 1 ∈ Θ 1 , unspecified innovation density). The sequential order identification procedure then consists of performing such tests for p 0 = 0, 1, ... and p 1 = p 0 + 1 until the null hypothesis H (n) 0 no longer gets rejected.
In this section, we state the local asymptotic normality (LAN) result on which the construction of our tests heavily relies. LAN, of course, requires some regularity assumptions on f ; we throughout assume the following (essentially borrowed from Garel and Hallin (1995) ).
The innovation density f belongs to the class F d of non-vanishing (with respect to the Lebesgue measure µ d on R d ) densities 3 i.e., for all c ∈ R + , there
(iii) f 1/2 is mean-square differentiable with mean-square gradient Df 1/2 , that is, there exists a square-integrable vector Df 1/2 such that, for all sequence h ∈ R d such
(v) the function x → ϕ f (x) is piecewise Lipschitz, i.e., there exists a finite measurable partition of R d into J non-overlapping subsets I 1 , . . . , I J such that, for all x, y
Denote by F * d the family of innovation densities satisfying Assumption 1.
denote the residuals computed from the initial values X −p+1 , . . . , X 0 , the parameter value θ 0 , and the observations X (n) .
Clearly, X (n) is the finite realization of a solution of (2.3) with parameter value θ 0 iff the residuals Z (n) 1 (θ 0 ), . . . , Z (n) n (θ 0 ) and 1 , . . . , n coincide. Denoting by P (n) θ 0 ;f the distribution of X (n) under parameter value θ 0 and innovation density f , the residuals Z
θ 0 +n −1/2 τ (n) /θ 0 ;f for the log-likelihood ratio log dP (n) 
where M θ 0 , P θ 0 , and Q (n) θ 0 are given in Appendix A ((A.1) and (A.2)) and
with the so-called f -cross-covariance matrices
We then have the following LAN result (Hallin et al., 2019) .
Proposition 2.1. Let Assumption 1 hold. Then, for any θ 0 ∈ Θ 0 and any bounded sequence τ (n) 
Parametrically efficient (in the Hájek-Le Cam asymptotic sense) rank-and sign-based (hence, distribution-free) inference procedures in LAN families (with given f ) are possible when the LAN central sequence (2.4) can be expressed in terms of signs and ranks. More precisely
θ 0 ;f terms), when there exists some ∆ (n) f (θ 0 ) measurable with respects to the ranks and signs of the residuals Z
In the univariate setting, this idea was exploited in Hallin and Puri (1988 , 1991 , 1994 , Hallin and Mélard (1988) , and Garel and Hallin (1999) , leading to a fairly complete toolbox of distribution-free procedures based on traditional ranks or signed ranks for ARMA models, in Hallin (1992, 1996) and La Vecchia (2017, 2020) for nonlinear time series models. In the multivariate context, under the assumption of ellipticity, Hallin and Paindaveine (2002 , 2006 , 2008 for location, regression, and shape, Hallin et al. (2010 Hallin et al. ( , 2013 Hallin et al. ( , 2014 for principal and common principal components, and Paindaveine (2004a,b, 2005) for VAR and VARMA models, proposed test procedures based on the so-called Mahalanobis ranks and signs. However, the assumption of elliptic innovation density (which reduces to symmetry in the univariate case) severely restricts the applicability of these tests, specially in the analysis of economic data, where actual innovations, typically, are skewed and/or leptokurtic. Thus, with the aim of accommodating a more general and realistic setting of non-elliptical distribution, we develop here novel test 
Let P d denote the family of all distributions P with densities in F d . For P in this family, the center-outward distribution functions defined below are continuous: see Hallin et al. (2020a) .
More general cases are studied in del Barrio et al. (2020) , but require more cautious and less intuitive definitions which, for the sake of simplicity, we do not consider here. Denote by U d the spherical uniform distribution over S d , that is, the product of a uniform measure over the hypersphere S d−1 and a uniform over the unit interval of distances to the origin.
The center-outward distribution function F ± of P is defined as the a.e. unique gradient of convex function mapping R d to S d and pushing P forward to U d . 4 For P ∈ P d , such mapping is a homeomorphism between S d \ {0} and R d \ F −1 ± ({0}) (Figalli 2018) and the corresponding center-outward quantile function is defined (letting, with a small abuse of notation, Q ± (0) := F −1 ± ({0})) as Q ± := F −1 ± . For any given distribution P, Q ± induces a collection of continuous, connected, and nested quantile contours and regions; the centeroutward median Q ± (0) is a uniquely defined compact set of Lebesgue measure zero. We refer to Hallin et al. (2020a) for details.
Turning to the sample, the residuals Z
with density f ∈ F d and center-outward distribution function F ± . For the empirical counterpart F
of F ± , let n factorize into n = n R n S + n 0 , for n R , n S , n 0 ∈ N and 0 ≤ n 0 < min{n R , n S }, where n R → ∞ and n S → ∞ as n → ∞, and consider a sequence G (n) of grids, where each grid consists of the intersection between an n S -tuple (u 1 , . . . u n S ) of unit vectors, and the n R hyperspheres with radii 1/(n R + 1), . . . , n R /(n R + 1) centered at the origin, along with n 0 copies of the origin. The only requirement 5 is that the sequence G (n) of grids is such that the discrete distribution with probability masses 1/n at each gridpoint and probability mass n 0 /n at the origin converges weakly to the uniform U d over the ball S d . Then, we define F (n) ± (Z (n) t ), for t = 1, . . . , n as the solution (optimal mapping) of a coupling problem between the residuals and the grid. Specifically, the empirical center-outward distribution function is the (random) discrete mapping
. . , n} coincides with the n points of the grid, · stands for the Euclidean norm, and T denotes the set of all possible bijective mappings between Z (n) 1 , . . . , Z (n) n and the n gridpoints. Based on this empirical center-outward distribution function, the center-outward ranks are defined as
(2.9)
the center-outward signs as
(2.10)
(2.11)
Those ranks and signs are jointly distribution-free (for f ∈ F d ): more precisely, under P (n)
n ) is uniformly distributed over the n!/n 0 ! permutations with repetition of the n underlying gridpoints (the origin having multiplicity n 0 ). Moreover, the center-outward distribution functions, ranks, and signs inherit, from the invariance properties of Euclidean distances, elementary but remarkable invariance and equivariance properties: Hallin et al. (2020b) show that center-outward quantities enjoy invariance/equivariance with respect to shift, global scale, and orthogonal transformations.
center-outward rank statistics
Throughout this section, denote by F ±,t := F ± (Z (n) t (θ 0 )) the value of the center-outward distribution function F ± associated with innovation density f computed at Z t = Z (n) t (θ 0 ) and by F (n)
where the function a :
where E θ 0 ;f denotes expectation under P ±,n ) is uniform over the permutations of the grid it was constructed from (that is, under any innovation density
θ 0 ;f only; none of them depends on i.
An essential step in Hájek's approach to the asymptotics of univariate ranks is the so-called
Hájek asymptotic representation of linear rank statistics which, contrary to earlier approaches based on empirical processes, allows for unbounded square-integrable score functions. Such results first were obtained (Hájek 1961; Hájek and Šidák 1967) for the nonserial rank statistics used in the context of linear models (single-output regression, ANOVA, etc.) ; Hallin et al. (1985) extend them to serial rank statistics-with scores involving several ranks at a time-appearing in the analysis of time-series models. The matrices Γ ∼ (n) i,J 1 ,J 2 are of that serial type, hence require an extension of Hallin et al. (1985) and Hallin and Vermandele (1996) (which only deal with classical "univariate ranks") to the present case of multivariate center-outward ranks and signs. Proposition 3.1 establishes such asymptotic representation (n) ) under the following assumption on the score functions J 1 and J 2 .
Assumption 2. (i) J 1 and J 2 are continuous over S d ; (ii) J 1 and J 2 are square-integrable, that is, S d J (u) 2 dU d < ∞ for = 1, 2, and, for any
n } of n-tuples in S d such that the uniform discrete distribution over s (n) converges weakly to U d as n → ∞,
When J (u) is of the form (spherical score) J ( u ) u u with J : [0, 1) → R continuous, a sufficient condition for (3.4) is the traditional assumption of bounded variation (that is, J is the continuous difference of two nondecreasing functions).
Proposition 3.1. Let Assumption 2 hold. Then, for any positive integer i,
See Section 8 for the proof. 
and (still with E θ 0 ;f denoting expectation under P
We then have, for a a and a e , the following lemma.
and ξ t,s := vec(J 1 (F ±,t )J 2 (F ±,s )).
and, for 1 ≤ i ≤ n − 1,
the following proposition follows as a corollary to Proposition 3.1 and Lemma 3.1.
and, consequently,
, as n → ∞. Assume, without loss of generality, that J 1 and J 2 are such that m = 0 (a sufficient
a and m (n) e both are o(n −1/2 ) (see Lemma 1 in Hallin and La Vecchia (2017) ) and m, m (n) a , and m (n) e can be omitted in Proposition 3.2. That simplification, however, is no longer valid when d > 1.
The proofs of Proposition 3.1 and Lemma 3.1 are given in Appendix B.
The basic idea underlying the construction of our rank-based tests is the definition of rank-
Writing
where J 1 and J 2 : S d → R are score functions satisfying Assumption 2. Call Γ
a (residual) lag-i rank-based cross-covariance matrix. As an approximate-score rank-based
structed from the center-outward distribution function and under P (n)
, hence constitutes a rank-based version of the same central sequence. These scores depend on f : if f stands for the actual innovation density, thus, they are oracle scores. They can be computed, however, at any chosen reference density such as the spherical or skew Gaussian; see Section 5.3. With a slight abuse of terminology, irrespective of J 1 and J 2 , we therefore call ∆ (n)
f (θ 0 ) conditional on the σ-field generated by the centeroutward ranks and signs-that is, the projection of ∆ (n) f (θ 0 ) onto the space of center-outward ranks and signs. We thus can state the following corollary. 
The joint asymptotic normality of the statistics (
and their asymptotic linearity, which entail the asymptotic normality and asymptotic linearity of ∆ (n) J 1 J 2 (θ 0 ), are essential in the construction of our rank-based tests. In view of the Hájek asymptotic representation result of Section 3, it is sufficient to derive these asymptotic results
under Assumptions 1 and 2, K J 1 ,J 2 ,f exists and is finite by the Cauchy-Schwarz inequality.
Decompose the matrix Q
(note that these blocks depend on i, not on n). Lemma 4.1 states the asymptotic normality
Lemma 4.1. Let Assumption 2 hold. Then, for any positive integers i = j, the vectors
Piecing together Lemma 4.1, the definition of ∆ (n) J 1 J 2 (θ 0 ), and Corollary 4.1, we then readily obtain the asymptotic normality of ∆ (n)
where the limits exist due to the exponential decrease of the Green matrices; see Appendix A for details. Note that when
Finally, in order to construct our test statistics, we still need the asymptotic linearity of ∆ (n) J 1 J 2 (θ 0 ); the latter is an immediate consequence of the asymptotic linearity for all i
(θ 0 )), which we now assume; the form of the linear term in the right-hand side of (4.7) follows from the form of the asymptotic shift in Lemma 4.1.
Assumption 3. For any positive integer i, as n → ∞,
The asymptotic linearity of ∆ (n) J 1 J 2 (θ 0 )) readily follows.
Proposition 4.2. Let Assumptions 1, 2, and 3 hold. Then, as n → ∞,
In this section, based on central sequences of the form ∆ (n) J 1 J 2 (θ 0 ), we propose centeroutward rank-based tests for VAR models and, depending on the scores J 1 and J 2 , derive their asymptotic properties. We separately consider two cases: (a) testing the null hypothesis θ = θ 0 , with specified θ 0 (VAR order p 0 ), against the alternative of the form θ = θ 0 with VAR order p (p potentially larger than p 0 ); since θ 0 is specified, rank-based statistics are fully distribution-free under the null (any f ∈ F d );
(b) testing the null hypothesis of a VAR of order p 0 (with unspecified parameter θ 0 ) against the alternative of a VAR of order p > p 0 (with p ≤ p 1 and unspecified parameter θ).
For p 0 > 0, the null value of θ 0 will be estimated; rank-based statistics then are only strongly asymptotically distribution-free, that is, asymptotically equivalent, under the null (any f ∈ F * d ), to a strictly distribution-free statistic. Case (a) is essentially of theoretical interest, and a preparation for case (b). The latter has an obvious application in the sequential identification of the order of VAR models via a sequence of tests of the null hypotheses of VAR(p 0 ) against VAR(p 0 + 1) dependence, p 0 = 0, 1, . . .
We first consider the case of the null hypothesis of θ = θ 0 , with θ 0 being specified, against the alternative of VAR(p), p 0 ≤ p ≤ p 1 with parameter θ = θ 0 . More precisely, the null
The center-outward rank-based test statistics we are proposing are of form 6 S (n)
The following proposition establishes the asymptotic distribution of S (n)
θ 0 and contiguous alternatives of the form P (n)
. Also, in order to discuss local optimality properties, we consider scores J 1 = J 1;f := ϕ f • Q ± and J 2 = J 2;f := Q ± and denote by S
Proposition 5.1. Let Assumption 2 hold. Then,
with d 2 p 1 degrees of freedom and non-centrality parameter
Most null hypotheses of practical interest involve incompletely specified values of the parameter, though. In the problem of testing a VAR of order p 0 against a VAR of order p 6 See Hallin and Paindaveine (2004a) for a similar test statistic based on Mahalanobis ranks and signs.
with 0 < p 0 < p ≤ p 1 , for instance, the value θ 0 of the parameter under the null remains unspecified. A natural idea then consists in replacing the unspecified θ 0 with a root-n consistent (constrained) estimatorθ (n) in (5.1), yielding the aligned rank test statistic S (n)
Such plug-in, however, has an impact on the asymptotic distribution of S (n) J 1 J 2 ; that impact is neutralized via the following classical construction inspired from the traditional theory of Neyman C(α) tests; see Section 11.9 of Le Cam (1986) for details.
Considering the null hypothesis H (n) 0
where θ 0 is an arbitrary parameter value of the form (2.2), and ∆ (n) 
It is easy to see that ∆ Therefore, letθ (n) := ((vecÂ
) denote a constrained-that is, of the form (2.2)-root-n consistent (under P (n) θ 0 ;f ) 8 estimator of θ 0 . For technical reasons, assume thatθ (n) moreover is asymptotically discrete. 9 Also, denote byΥ (n) a consistent (under P (n) θ 0 ;f ) estimator of Υ J 1 ,J 2 ,f (θ 0 ). Such an estimator can be obtained, e.g., from the asymptotic linearity in (4.8): see Section 8 for details. Our test statistic is obtained by plugging these estimators into (5.6), yielding
). We then have the following asymptotic results. 
is asymptotically non-central chi-square with d 2 (p 1 −p 0 ) degrees of freedom and non-centrality parameter
In view of contiguity, this also holds under P (n) θ0+n −1/2 τ ;f with τ of the form (2.2). 9 An estimatorθ (n) of θ 0 is called asymptotically discrete if, for n large enough, it takes at most K(c) distinct values in balls of the form {θ ∈ R p1d 2 : n −1/2 θ − θ 0 ≤ c}, c > 0 fixed, K(c) independent of n. Any root-n consistent estimator is easily discretized into an asymptotically discrete one. Asymptotic discreteness, however, is only a theoretical requirement since, in practice,θ (n) anyway only has a finite number of digits. While discretization is necessary in asymptotic statements, it is not needed in finite-n implementation of tests; see Le Cam and Yang (2000, Chapter 6 ) and van der Vaart (1998, Section 5.7) for details. J 1 J 2 that can be used to construct uniformly valid (in the sense of (1.2)) asymptotic critical values. Propositions 5.1 (ii) and 5.2 (iii) provides asymptotic distributions under contiguous alternatives, allowing for local power evaluations. Propositions 5.1 (iii) and 5.2 (iv) establish the optimality properties of the same tests.
For p 0 = 0, that is, when testing white noise against VAR(p 1 ) dependence, ∆ (n)
and W (n) J 1 J 2 (0) coincide and both take the form
(5.8) If p 1 = 1, Proposition 5.2 implies that the null hypothesis of white noise is rejected in favor of VAR(1) dependence whenever S (n)
J 1 J 2 exceeds the (1 − α) quantile χ 2 d 2 ;1−α of the chi-square distribution with d 2 degrees of freedom.
The rank-based cross-covariance matrices Γ ∼ (n) J 1 J 2 , hence also the test statistics S (n)
and W (n) J 1 J 2 (θ (n) ), depend on the choice of score functions J 1 and J 2 , to be selected by the practitioner. In principle, in view of the maximin and stringency properties in Propositions 5.1 (iii) and 5.2 (iv), an optimal choice is J 1 = ϕ f • Q ± and J 2 = Q ± , where f is the actual innovation density and Q ± the corresponding center-outward quantile function.
Such a choice, unfortunately, is unfeasible since f is unspecified. Alternatives can privilege either simplicity, robustness, or efficiency at some chosen reference density. Here are three examples of sensible choices extending scores that are widely applied in the univariate (see e.g. Hallin and La Vecchia (2017) ) and the elliptical multivariate setting (see Hallin and Paindaveine (2004a) ). Note that for these scores, the centering term m (n) a for Γ
disappears as soon as the grid G (n) is symmetric with respect to the origin while m = 0 since J (−u) = −J (u), = 1, 2.
±,t , = 1, 2 yields the center-outward sign-based cross-covariance matrices
The resulting ∆ (n) sign (θ) entirely relies on the center-outward signs S (n) ±,t (θ), which thus explains the terminology sign test scores. Since the ranks do not enter the score function, there is no need for the grid G (n) to fill the unit ball as n → ∞, and we can set n S = n, producing a grid over the unit hypersphere rather than the unit ball.
Example 2 (Spearman scores). Another simple choice is J 1 (u) = u = J 2 (u). The corresponding rank-based cross-covariance matrices are
reducing, for d = 1, to Spearman autocorrelations, whence the terminology Spearman scores.
, where F χ 2 d denotes the chi-square distribution function with d degrees of freedom. This yields the spherical van der Waerden (vdW) rank scores, with cross-covariance matrices
In view of Propositions 5.1 and 5.2, the resulting tests are asymptotically optimal under spherical Gaussian innovations.
One also may wish to consider more sophisticated reference densities, such as skew-normal or skew-t ones. The problem then is the absence, for non-spherical densities, of a closed-form expression for Q ± . 10 That problem, in principle, can be overcome by means of a numerical evaluation of Q ± : increasingly powerful algorithms indeed exist for the numerical computation of optimal transports. Assume that the (unspecified) innovation covariance matrix Σ exists and is finite and consider, for instance, the scores associated with a Gaussian N (0, Σ) reference density (for which ϕ Σ (x) = Σ −1 x):
where Q N (0,Σ);± denotes the center-outward quantile function associated with the N (0, Σ) density. That quantile function is analytically unknown, but can be evaluated (at the G (n) gridpoints) with arbitrary precision by (a) estimating Σ with some consistent estimator Σ (n) measurable with respect to the order statistic Z ± and the order statistic, does not affect the distribution of the ranks and the signs. As n and N tend to infinity, these tests are locally asymptotically optimal under N (0, Σ) innovation density.
The same idea can be implemented to construct tests based on skew-normal or skew-t scores. Exploring 11 the feasibility and the benefits (relative efficiencies only can be obtained numerically) of such an approach, however, is beyond the scope of this paper, and is left for further research.
In this section, we briefly introduce the routinely-applied Gaussian test procedure (the hypothesis-testing counterpart of quasi-maximum-likelihood estimation (QMLE)), the performance of which will serve as a benchmark for our center-outward rank-based test procedure in the numerical analysis of Section 7. These tests are developed in Hallin and Paindaveine (2004a) . We are focusing on the two cases described in Sections 5.1 and 5.2. 11 A very preliminary investigation of the Gaussian case described above shows that the efficiency gains over spherical-Gaussian-score-based methods is quite small.
We first consider the Gaussian solution to the problem considered in Section 5.1, i.e., the test of the null hypothesis H 
N (θ 0 ) exceeds the (1 − α) quantile of the chi-square distribution with d 2 p 1 degree of freedom. That test is locally asymptotically maximin at asymptotic level α against Gaussian alternatives; see Section 6.3 of Hallin and Paindaveine (2004a) for details and a proof.
Next let us consider the Gaussian tests required for the sequential VAR order identification problem. The null hypothesis of interest is thus H 
We illustrate numerically the performance of our test statistics in two benchmark inferential problems: testing for serial dependency in a multivariate time series and selecting the order of autoregression in a VAR model. As discussed in Hallin and Paindaveine (2004a) , these testing problems are common to many econometric applications.
To elaborate further, consider the d-variate general regression model with VAR error
terms. Under this model, the observation is an n × d array
where the regression function g : R q × R k → R d has known functional form and depends on an unknown parameter β β β ∈ R k and on observed covariates v t ∈ R q . Instead of the traditional assumption that the error term is white noise, we rather assume {X t ; t = 1, . . . , n} to be the finite realization (of length n) of some VAR(p) process generated by
where { t ; t ∈ Z} is d-dimensional white noise with density f , satisfying Assumption 1.
A similar framework is considered in Paindaveine (2004b, 2005) , where the function g is linear in β β β (hence, without any loss of generality, in the covariates). Hallin and Paindaveine moreover assume that the innovation density f is elliptical and consider test procedures based on Mahalanobis ranks and signs. This very restrictive assumption of elliptical symmetry is precisely what we are dropping here as we only assume that f ∈ F * d .
A classical problem in single-output linear regression is the so-called Durbin-Watson problem of testing the null hypothesis under which X t , in the single-output linear version of (7.1) with unspecified β β β, is second-order white noise against the alternative of VAR(1) dependence: see Watson (1950, 1951) . Here, we extend that problem to nonlinear multiple-output regression and propose a rank-based solution. For the sake of simplicity, we limit ourselves to VAR(1) alternatives (writing A instead of A 1 ), but extensions to higher-order VAR(p) are straightforward.
It is easy to see that, under Assumption 1 and appropriate conditions on g and the asymptotic behavior of the covariates, model (7.1) is LAN with block-diagonal information matrix (a nonserial β β β-block and a serial one related with θ = vecA). LAN then implies that the unknown β β β safely can be replaced by any root-n consistent (and, in principle, locally discrete) estimatorβ β β (n) without having any local impact on the θ-part of the central sequence (see Hallin and Paindaveine (2004b) for the case of a linear g; the only difference is that the β β β part of the central sequence here has an additional g factor which does not affect that fact that its covariance with the θ part vanishes). As a consequence, one can construct a test as of β β β were known, then safely replace β β β withβ β β (n) in the resulting test statistic.
The traditional implementation of Durbin-Watson tests involves least-squares estimators of β β β. We rather suggest robust root-n consistent M-estimatorsβ β β (n) such as the boundedinfluence M-estimator for dynamic location models proposed in Mancini et al. (2005) .
The Durbin-Watson problem thus is about testing θ = 0 against θ = 0, on the basis of
Under the null hypothesis, the X t 's are white noise; under the alternative, they are VAR(1).
The Gaussian solution (which is routinely applied) is based on the Gaussian test statistic (6.5),
its rank-based competitors on test statistics of the form (5.8), for which we consider the spherical van der Waerden, Spearman, and sign scores as described in Section 5.3.
To simplify the numerical implementation and avoid specifying any particular g, our simulations directly proceed with the residuals X t in (7.1). We thus generated N = 1000
replications of size n = 800 (n R = 20, n S = 40, and n 0 = 0) from the bivariate models 
with µ 1 = (−5, 0) , µ 2 = (5, 0) , µ 3 = (0, 0) and Σ 1 = 7 5 5 5 , Σ 2 = 7 −6 6 6 , Σ 3 = 4 0 0 3
The d-dimensional skew-t ν distribution has density f (z; ξ, Σ, α, ν)
(see Azzalini and Capitanio 2003) where x = w −1 (z − ξ); w = (diag(Σ)) 1/2 , ξ ∈ R d , and α ∈ R d are scale, location, and shape parameters respectively, T (y; ν) denotes the univariate t ν distribution function, and
Here (d = 2), the values ξ = 0, α = (5, 2) , and Σ = 7 4 4 5 were selected. (1); the sample size is n = 800 (n R = 20, n S = 40, and n 0 = 0); the nominal level is α = 5%; rejection frequencies of the sign tests with n R = 1 are also included. Permutational critical values are based on M = 5000 random permutations. To investigate the robustness of the center-outward rank-based tests, we also considered the case of spherical Gaussian X t 's contaminated by additive outliers (AOs). More precisely, we generated observations of the form {X * t = X t + 1 {t=h} s}, where h and s denote the location and size of the AOs, respectively. We set h in order to obtain 5% of equally spaced AOs and put s = (6, 6) and (9, 9) . All contaminated observations were demeaned prior to the implementation of the testing procedures. Rejection frequencies at 5% nominal level are reported in Table 1 for the Gaussian, vdW, Spearman and sign tests based on the asymptotic critical values provided in Sections 5.1 and 6.1 but also for the vdW and Spearman tests based on bias-corrected critical values as described below. Inspection of Table 1 reveals that, under normal innovations, the rejection frequencies under the null hypothesis of the Gaussian test are very close to the nominal 5% size. But the vdW and Spearman tests are quite below that nominal size, indicating a severe bias. This bias is confirmed by Figure 1 , which shows the QQ plots of the values of the vdW test statistic across the N = 1000 replications with Gaussian f (but f has no impact here, as the ranks are distribution-free under the null). The tails of the vdW statistic, for instance, very clearly do not match those of the chi-square they are converging to. The reason for this is the relatively slow convergence of the distribution of center-outward rank statistics: n R = 20 for n = 800 is still a rather small n R value, which explains the phenomenon. Note that the same QQ plot for the sign test score statistic, which does not depend on the ranks, are more Table 1 under the label "bias-corrected," for the vdW and Spearman tests.
The vdW and Spearman tests have greater power than the sign test and slightly less power than the Gaussian one under Gaussian innovations. Under t 3 innovations, however, the vdW and Spearman tests, despite their bias, both outperform the Gaussian test. The bias-corrected vdW and Spearman tests have correct size under the null and even higher powers; under the mixture and skew-t 3 innovations, they outperform the Gaussian test by a landslide. The tests based on signs have been implemented with the same factorization (viz., n R = 20, n S = 40, and n 0 = 0) as vdW and Spearman; the factorization (n R = 1, n S = 800, and n 0 = 0) also has been considered, but provides little improvement. Turning to robustness issues, the resistance to additive outliers of the rank-based tests over the Gaussian Table 1 against their asymptotically chi-square null distributions: vdW (left), Gaussian, Spearman, and sign test scores (right). one is extremely good, while the Gaussian test literally crashes, with exploding size under the null-the corresponding rejection frequencies under the alternative thus are meaningless.
Despite their simplicity, the sign tests do extremely well (confirming univariate robustness results by Boldin (2012 Boldin ( , 2013 ) but remain less powerful than the bias-corrected vdW tests.
Bias-corrected tests uniformly outperform their uncorrected asymptotic counterparts.
Turning to VAR order identification, let us consider the bivariate VAR(1) model with autoregression matrix coefficient vec(A) = (0.30, −0.06, 0.12, 0.24) , the same innovation densities and the same contamination schemes as in Section 7.1, and sample size n = 800.
The sequential method of Sections 5.2 and 6.2 was applied to N = 1000 replications thereof.
Each step consists in testing (at 5% nominal level) a VAR(p) against a VAR(p + 1), which requires the estimation of the null VAR(p) model; throughout, this estimation is performed via the center-outward R-estimator based on vdW scores as proposed by Hallin et al. (2019) , where it is shown to be root-n consistent. The numbers of under-, correct, and over-identification for the Gaussian, (corrected) vdW, (corrected) Spearman, and (corrected) sign tests are reported in Table 2 .
All tests are rejecting the hypothesis of white noise, irrespective of the underlying innovation density or the presence of outliers: no under-identifications, thus. Also, the rankbased procedures yield uniformly lower over-identification rates than the Gaussian one. Under Table 2 : Numbers of under-, correct, and over-identification for the VAR(1) model under various innovation densities and additive outlier schemes, of the Gaussian, vdW, bias-corrected vdW, Spearman, bias-corrected Spearman, and sign tests out of N = 1000 replications; the sample size is n = 800 (n R = 20, n S = 40, and n 0 = 0); results of the sign tests with n S = n are also included; permutational crirical values are based on M = 1000 random permutations.
Order identification Order identification the t 3 distribution, the Gaussian procedure yields an almost 10% over-identification rate,
indicating that its performance is badly affected by heavy tails. The rank-based procedures, on the other hand, do much better under heavy tails, skewness, and additive outliers.
To illustrate the applicability of our VAR order selection procedure, we consider the joint To start with, in the top panels of Figure 2 , we display the seasonally adjusted time series.
Even a visual inspection reveals some interesting characteristics. A first clearly visible aspect is that both time series display a trend in time. To deal with this, we differentiate the series and plot the resulting outputs in the bottom panels of Figure 2 . We notice that the differentiated series seem to be correlated and display common movements (mainly in opposite directions).
A second noticeable aspect is that, starting from 2005, the trajectories of the differentiated series look increasingly asymmetric and spiky. The period March-September-2020 reveals a cluster of outlying values, most likely due to COVID-19-related policy decisions. we also ran, on the same series, the rank-based order-identification procedures described in Section 7.2. The results are shown in Table 3 . Table 3 reveals a sharp contrast between the conclusions of the Gaussian and rank-nased methods. The Gaussian procedure indeed very clearly selects (at nominal level 5%) a VAR(2) model while all rank-based procedures, whether bias-corrected (with M = 1000 permutations) or not, agree on a larger number of 6-7 lags.
Our educated guess is that the combination of skewness, kurtosis, and outliers are blurring the conclusions of Gaussian tests; the graphical diagnostics in Figure 3 bring some evidence in favor of that guess. Figure 3 displays the residuals of a VAR (7) 
This paper introduces rank-based tests for VAR models with unspecified innovation density.
Based on the residual center-outward ranks and signs recently proposed by Chernozhukov et al. (2017) and (Hallin et al., 2020a , where they are shown to be essentially maximal ancillary), they constitute the testing counterpart of the R-estimators proposed in Hallin et al. (2019) . When testing a VAR with specified parameter θ 0 , they are fully distribution-free;
when θ 0 remains unspecified under the null, they are strongly asymptotically distribution-free in the sense of being asymptotically equal to a fully-distribution-free test. When based on appropriate score functions, they achieve parametric efficiency at selected reference densities.
Monte Carlo experiments and an empirical example demonstrate the excellent performance of the proposed tests, which significantly outperform the traditional pseudo-Gaussian methods under skew, heavy-tailed, and contaminated innovations.
Some linear difference equation algebra is required in order to obtain the explicit form of the central sequences and information matrices in the LAN result of Section 2.2-to be used in the construction of our center-outward rank-based test statistics.
Denote by G u , u ∈ Z the Green's matrices associated with the difference operator A θ 0 (L), with θ 0 ∈ Θ 0 ⊂ Θ 1 ⊂ R p 1 d 2 for some p 0 < p 1 ; those matrices are defined as the solutions of the homogeneous linear recursions
with initial values I d , 0, . . . , 0 at u = 0, −1, . . . , −p 0 + 1.
For θ 0 ∈ Θ, the decrease of { G u , u ∈ N} as u → ∞ is exponential. Specifically, there exists some ε > 0 (depending on the smallest root of the equation det A θ (z) = 0) such that G u (1 + ε) u converges to 0 as u → ∞. Denoting by X −p 0 +1 , . . . , X 0 the initial values (which are typically unobservable) of a solution of the homogenous difference equation A θ 0 (L)X t = 0, this exponential decrease ensures that these initial values have no asymptotic influence on asymptotic results and, therefore, safely can be set to zero in the sequel. This allows us to invert the autoregressive polynomial; the Green matrices G u then are the matrix coefficients 12 of the inverted operator A θ 0 (L) −1 :
, z ∈ C, |z| < 1.
More generally, associated with an arbitrary d-dimensional linear difference operator C(L) := ∞ i=0 C i L i (this of course includes operators of finite order s), define, for any integers v ≤ u, the d 2 u × d 2 v matrices (the form of which is motivated by the vecA u form under which the VAR matrix coefficients A u enter θ 0 and the d 2 v-dimension nature of the 12 In the econometric terminology, the impulse response matrix coefficients.
perturbations τ on the LAN property )
With this notation, G u,u (associated with C(L) = A θ 0 (L) −1 ) is the inverse of A u,u (associated with C(L) = A θ 0 (L)). Denoting by C u,v the matrices associated with the transposed operator C (L) := ∞ i=0 C i L i , we also have that G u,u = (A u,u ) −1 . Then, for any θ 0 ∈ Θ 0 and any p 1 > p 0 , the d 2 p 1 × d 2 p 1 matrix
is of full rank.
Still associated with A θ 0 (L), consider the operator D θ 0 (L) :
(1) t , . . . , ψ (p 0 ) t } be any set of d × d matrices forming a fundamental system of solutions of the homogeneous linear difference equation associated with D θ 0 (L). Such a system can be obtained, for instance, from the Green matrices of D θ 0 (L) (see, e.g., Hallin 1986). Defining
Proof. The proof follows from an enumeration of the n(n − 1) terms in Proof. Lemma B1 implies
The desired result then follows from the fact that C 
±,s )), ξ t,s := vec(J 1 (F ±,t )J 2 (F ±,s )), γ (n) := vec(m (n) ), and γ := vec(m).
We have
Noting that D
In order to show that D (.) = z (n) ∈ R nd , the distribution of α (n) 2,1 thus is that of a score function (the form of which depends on z (n) ) computed at the empirical center-outward ranks and signs. Hence, in view of Lemma B2,
where the second term above is o(n −1 ) since D
It remains to show that D
2,1 and γ = E ξ 2,1 ,
(.) ) − γ) = E E(ξ 2,1 |Z (n) (.) )E(ξ 2,1 |Z (n) (.) ) − E(ξ 2,1 )E(ξ 2,1 ).
Noting that E(ξ 2,1 |Z 1 n 2 (n − 1) 2 1≤j 1 =j 2 ≤n,1≤j 3 =j 4 ≤n E(ξ j 1 ,j 2 ξ j 3 ,j 4 ) − E(ξ 2,1 )E(ξ 2,1 ) = 1 n 2 (n − 1) 2 1≤j 1 =j 2 ≤n,1≤j 3 =j 4 ≤n E [(J 1 (F ±,j 1 )J 1 (F ±,j 3 ))(J 2 (F ±,j 2 )J 2 (F ±,j 4 ))]
− E(J 1 (F ±,2 )) 2 E(J 2 (F ±,1 )) 2 .
Since Z (n) 1 , . . . , Z (n) n are i.i.d., only 2n(n − 1) terms with j 1 = j 2 and j 3 = j 4 contribute to this latter summation. Therefore, in view of the square-integrability of J 1 and J 2 , we have D (n) 31 = O(n −2 ). The result follows.
Proof of Lemma 3.1. In view of LAN, it is sufficient to prove the result under P (n) θ 0 ;f (with f ∈ F d ). It follows from the continuity of J 1 and J 2 and the Glivenko-Cantelli theorem in Hallin et al. (2020a) that ξ (n) 2,1;a converges to ξ 2,1 a.s.. Moreover, square-integrability of J 1 and J 2 and independence of Z 1 and Z 2 entail E ξ 2,1;a 2 = E J 1 (F ±,2 ) 2 E J 2 (F ±,1 ) 2 < ∞.
Hence, E ξ (n) 2,1;a 2 − E ξ 2,1;a 2 = o(1). It follows (see, e.g., part (iv) of Theorem 5.7 in Chapter 3 of Shorack (2000)) that E ξ (n) 2,1;a − ξ 2,1 2 = o(1).
Turning to the second part of (3.5), put ξ 2,1 (ω) = vec(J 1 (F ±,2 (ω))J 2 (F ±,1 (ω))), ξ (n) 2,1;a (ω) = vec(a a (F (n) ±,2 (ω), F (n) ±,1 (ω))) for any point ω in the sample space Ω. Since ξ (n) 2,1;a − ξ 2,1 converges to zero a.s., it follows from the Egorov theorem (see, e.g., part (ii) of Exercise 5.8 in Shorack (2000, Chapter 3)) that, for any ε > 0, there exists a subset A ⊂ Ω such that P(A) > 1 − ε and sup where, due to the square-integrability of J 1 and J 2 and the independence between Z 1 and Z 2 , I (n) 2 is arbitrarily small as ε → 0.
It remains to prove that I (n) 1 = o(1). Denoting by G (n) the regular grid the empirical center-outward distribution function is mapping to, define, for any g 1 , g 2 ∈ G (n) , B (n) g 1 g 2 := {ω : F (n) ±,1 (ω) = g 1 , F
(n) ±,2 (ω) = g 2 }.
We have
where the second equality follows from the fact that 1 B (n)
h 1 h 2 = 1 (g 1 =h 1 ,g 2 =h 2 ) . Note that
±,1 = 1 B (n) g 1 g 2 P(B (n) g 1 g 2 ) η∈B (n) g 1 g 2 ξ 2,1 (η)1 A (η)dP(η).
Therefore,
(ω) η∈B (n) g 1 g 2 ξ 2,1 (ω)1 A (ω) −ξ 2,1 (η)1 A (η) dP(η) P(B (ω) = vec(J 1 (g 2 )J 2 (g 1 )) = ξ (n) 2,1;a (η)1 A (η)1 B (n) g 1 g 2 (η).
In view of (B1), for anyε > 0, there exists n(ε) such that ξ 2,1;a (ω) − ξ (n) 2,1 (ω) 2 <ε for all n > n(ε) and all ω ∈ A. Then, using Jensen's inequality again, I (n) 1 ≤ g 1 ∈G (n) g 2 ∈G (n) Ω 1 B (n) g 1 g 2
(ω) η∈B (n) g 1 g 2 2 ξ 2,1 (ω) − ξ (n) 2,1;a (ω) 2 1 A (ω) +2 ξ (n) 2,1;a (η) − ξ 2,1 (η) 2 1 A (η) dP(η) P(B (n) g 1 g 2 ) dP(ω) ≤ g 1 ∈G (n) g 2 ∈G (n) Ω 1 B (n) g 1 g 2
(ω)4ε 2 P(B (n) g 1 g 2 ) P(B (n) g 1 g 2 ) dP(ω) = 4ε 2 E g 1 ∈G (n) g 2 ∈G (n) 1 B (n) g 1 g 2 = 4ε 2 .
The result follows. with Gaussian f . The same relatively slow asymptotics of center-outward ranks are to be blamed-the value of n R = 21 for n = 1400 indeed is quite small-and the problem similarly is solved by resorting to permutational critical values.
be constructed, though, such as the low-discrepancy sequences (see, e.g., Niederreiter (1992) , Judd (1998) , or Dick and Pillichshammer (2014) ) considered in numerical integration and the design of computer experiments; we suggest the use of the function UnitSphere in R package mvmesh.
(ii) The empirical center-outward distribution function F (n) ± is obtained as the solution of an optimal coupling problem. Many efficient algorithms have been proposed in the measure transportation literature (see, e.g., Peyré and Cuturi (2019)). We followed Hallin et al.
(2020a), using a Hungarian algorithm (see the clue R package). For a detailed account of the Hungarian algorithm and the complexity of different implementations, see, e.g., Chapter 4
in Burkhard et al. (2009) ). Faster algorithms are available, though, as Bertsekas' auction algorithm or its variant, the forward/reverse auction algorithm, (Chapter 4 in Bertsekas (1991)), implemented in the R package transport.
(iii) The computation of the test statistic W (n) J 1 J 2 (θ (n) ) involves two basic ingredients: a root-n consistent estimatorθ (n) of θ 0 and a consistent estimator of the matrix Υ J 1 ,J 2 ,f (θ 0 ). Forθ (n) , a great number of candidates are available, e.g., the center-outward R-estimator of Hallin et al. (2019) , the reweighted multivariate least trimmed squares estimator of Croux (2008) , and the QMLE (provided that fourth-order moments finite). Turning to the estimation of Υ J 1 ,J 2 ,f (θ 0 ), the issue is that this matrix depends on the unknown actual density f . A simple consistent estimator is obtained by letting τ = e i , i = 1, . . . , p 0 d 2 in (4.8) where e i denotes the ith vector of the canonical basis in the parameter space R p 1 d 2 : the difference ∆ 
For the sake of implementation, we summarize, in Algorithm 1 and Algorithm 2, respectively, the algorithms of the rank-based test procedures for the null of θ = θ 0 (with specified θ 0 ) and the VAR order identification (with unspecified θ 0 ).