In this work we provide a review of basic ideas and novel developments about Conformal Prediction — an innovative distribution-free, non-parametric forecasting method, based on minimal assumptions — that is able to yield in a very straightforward way prediction sets that are valid in a statistical sense also in the finite sample case. The discussion provided in the paper covers the theoretical underpinnings of Conformal Prediction, and then proceeds to list the more advanced developments and adaptations of the original idea.

In this paper, we explore the two star Exponential Random Graph Model, which is a two parameter exponential family on the space of simple labeled graphs. We introduce auxiliary variables to express the two star model as a mixture of the β model on networks. Using this representation, we study asymptotic distribution of the number of edges, and the sampling variance of the degrees. In particular, the limiting distribution for the number of edges has similar phase transition behavior to that of the magnetization in the Curie-Weiss Ising model of Statistical Physics. Using this, we show existence of consistent estimates for both parameters. Finally, we prove that the centered partial sum of degrees converges as a process to a Brownian bridge in all parameter domains, irrespective of the phase transition.

We consider the problem of matching a template to a noisy signal. Motivated by some recent proposals in the signal processing literature, we suggest a rank-based method and study its asymptotic properties using some well-established techniques in empirical process theory combined with Hájek’s projection method. The resulting estimator of the shift is shown to achieve a parametric rate of convergence and to be asymptotically normal. Some numerical simulations corroborate these findings.

We provide quantitative upper bounds on the total variation mixing time of the Markov chain corresponding to the unadjusted Hamiltonian Monte Carlo (uHMC) algorithm. For two general classes of models and fixed time discretization step size h, the mixing time is shown to depend only logarithmically on the dimension. Moreover, we provide quantitative upper bounds on the total variation distance between the invariant measure of the uHMC chain and the true target measure. As a consequence, we show that an ε-accurate approximation of the target distribution μ in total variation distance can be achieved by uHMC: (i) for a broad class of models with $O\left(d^{3∕4}{\mathit{\epsilon }}^{-1∕2}log(d∕\mathit{\epsilon })\right)$ gradient evaluations; and (ii) for mean field models with weak interactions with $O\left(d^{1∕2}{\mathit{\epsilon }}^{-1∕2}log(d∕\mathit{\epsilon })\right)$ gradient evaluations. The proofs are based on the construction of successful couplings for uHMC that realize the upper bounds.

High dimensional generalized linear models are widely applicable in many scientific fields with data having heavy tails. However, little is known about statistical guarantees on the estimates of such models in a time series setting. In this article, we establish statistical error bounds and support recovery guarantees of the classical ${\ell }_1$ regularized procedure for generalized linear model with temporal dependent data. We also propose a new robust M-estimator for high dimensional time series. Properties of the proposed robust procedure are investigated both theoretically and numerically. As an extension, we introduce a robust estimator for linear regression and show that the proposed robust estimator achieves nearly the optimal rate as that for i.i.d sub-Gaussian data. Simulation results show that the proposed method performs well numerically in the presence of heavy-tailed and serially dependent covariates and/or errors, and it significantly outperforms the classical Lasso method. For applications, we demonstrate, in the supplementary material, the regularized robust procedure via analyzing high-frequency trading data in finance.

We define and prove limit results for a class of dominant Pólya sequences, which are randomly reinforced urn processes with color-specific random weights and unbounded number of possible colors. Under fairly mild assumptions on the expected reinforcement, we show that the predictive and the empirical distributions converge almost surely (a.s.) in total variation to the same random probability measure $\tilde{P}$; moreover, $\tilde{P}(\mathcal{D})=1$ a.s., where $\mathcal{D}$ denotes the set of dominant colors for which the expected reinforcement is maximum. In the general case, the predictive probabilities and the empirical frequencies of any δ-neighborhood of $\mathcal{D}$ converge a.s. to one. That is, although non-dominant colors continue to be regularly observed, their distance to $\mathcal{D}$ converges in probability to zero. We refine the above results with rates of convergence. We further hint potential applications of dominant Pólya sequences in randomized clinical trials and species sampling, and use our central limit results for Bayesian inference.

We consider admissibility of generalized Bayes estimators of the mean of a p-variate normal distribution when the scale is unknown, and the loss is quadratic. The priors considered put the improper invariant prior on the scale while the prior on the mean has a hierarchical normal structure conditional on the scale. This conditional hierarchical prior is indexed by a hyperparameter, a. In earlier studies, the authors established admissibility/inadmissibility of the generalized Bayes estimator under the proper/improper conditional prior ($a>-1$ / $a<-2$), respectively. In this paper we complete the admissibility/inadmissibility characterization for this class of priors by establishing admissibility for the improper conditional prior ($-2\le a\le -1$). This boundary, $a=-2$, with admissibility for $a\ge -2$ and inadmissibility for $a<-2$ corresponds exactly to that in the known scale case for this class of conditional priors, and which follows from Brown’s 1971 paper. As a notable benefit of this enlargement of the class of admissible generalized Bayes estimators, we give admissible and minimax estimators for $p\ge 3$ as opposed to an earlier study which required $p\ge 5$. In one particularly interesting special case, we establish that the joint Stein prior for the unknown scale case leads to a minimax admissible estimator for $p\ge 3$.

We introduce a projection-based class of uniformity tests on the hypersphere. The class employs statistics that integrate, along all possible directions, a weighted quadratic discrepancy between the empirical cumulative distribution of the projections and the projected uniform distribution. Simple expressions for several test statistics are obtained for the circle and the sphere, as well as relatively tractable forms for higher dimensions. Despite their different origins, variants of the proposed class are shown to contain and be contained in variants of the Sobolev class of uniformity tests. Our new class proves itself advantageous by allowing the derivation of new tests that neatly extend the circular-only tests by Watson, Ajne, and Rothman, and by introducing the first instance of an Anderson–Darling-like test for directional data. We obtain usable asymptotic distributions and the local asymptotic optimality against certain alternatives of the new tests. A simulation study evaluates the theoretical findings and provides evidence that the new testing proposals are competitive. An application to the study of the crater distribution on Rhea illustrates the usage of the new tests.

Recently, Radulovic and Wegkamp introduced a new technique to show convergence in distribution of the empirical process indexed by functions of bounded variation. This method of proof allows to directly extend convergence results known for the canonical empirical process to convergence in distribution of the empirical process indexed by functions of bounded variation. The purpose of this article is twofold. First, we extend the mentioned technique to index functions of locally bounded variation. Second, and more importantly, we demonstrate that this technique provides a new approach to show convergence in distribution of the smoothed empirical process based on kernel density estimators. Using this approach we can prove to the best of our knowledge the first results on convergence in distribution of the smoothed empirical process of dependent data. Our results cover both weak and strong dependence as well as index sets of functions of locally bounded variation. Moreover our results cover an MISE optimal choice of the bandwidth for the kernel density estimator which to some extent is the plug-in property in the Bickel–Ritov sense. In the case of i.i.d. data our results extend a seminal result of Giné and Nickl.

We develop a class of tests for semiparametric vector autoregressive (VAR) models with unspecified innovation densities based on the recent measure-transportation-based concepts of multivariate center-outward ranks and signs. We show that these concepts, combined with Le Cam’s asymptotic theory of statistical experiments, yield novel testing procedures, which (a) are valid under a broad class of innovation densities (possibly non-elliptical, skewed, and/or with infinite moments), (b) are optimal (locally asymptotically maximin or most stringent) at selected ones, and (c) are robust against additive outliers. In order to show this, we establish, for a general class of center-outward rank-based serial statistics, a Hájek asymptotic representation result, of independent interest, which allows for a rank-based reconstruction of central sequences. As an illustration, we consider the problems of testing the absence of serial correlation in multiple-output and possibly non-linear regression (an extension of the classical Durbin-Watson problem) and the sequential identification of the order p of a VAR(p) model. A Monte Carlo comparative study of our tests and their routinely-applied Gaussian competitors demonstrates the benefits (in terms of size, power, and robustness) of our methodology; these benefits are particularly significant in the presence of asymmetric and leptokurtic innovation densities. A real-data application concludes the paper.

In Stein’s method, the exchangeable pair approach is commonly used to estimate the approximation errors in normal approximation. In this paper, we establish a Cramér-type moderate deviation theorem of normal approximation for unbounded exchangeable pairs. As applications, Cramér-type moderate deviation theorems for the sums of local statistics and general Curie–Weiss model are obtained.

The asymptotic variance is an important criterion to evaluate the performance of Markov chains, especially for the central limit theorems. We give the variational formulas for the asymptotic variance of discrete-time (non-reversible) Markov chains on general state space. The variational formulas provide many applications, extending the classical Peskun’s comparison theorem to non-reversible Markov chains, and obtaining several comparison theorems between Markov chains with various perturbations.

Principled nonparametric tests for regression curvature in ${\mathbb{R}}^d$ are often statistically and computationally challenging. This paper introduces the stratified incomplete local simplex (SILS) tests for joint concavity of nonparametric multiple regression. The SILS tests with suitable bootstrap calibration are shown to achieve simultaneous guarantees on dimension-free computational complexity, polynomial decay of the uniform error-in-size, and power consistency for general (global and local) alternatives. To establish these results, we develop a general theory for incomplete U-processes with stratified random sparse weights. Novel technical ingredients include maximal inequalities for the supremum of multiple incomplete U-processes.

The first essential ingredient to build up Stein’s method for a continuous target distribution is to identify a so-called Stein operator, namely a linear differential operator with polynomial coefficients. In this paper, we introduce the notion of algebraic Stein operators (see Definition 3.4), and provide a novel algebraic method to find all the algebraic Stein operators up to a given order and polynomial degree for a target random variable of the form $Y=h(X)$, where $X=(X_1,\dots ,X_d)$ has i.i.d. standard Gaussian components and $h\in \mathbb{K}[X]$ is a polynomial with coefficients in the ring $\mathbb{K}$. Our approach links the existence of an algebraic Stein operator with null controllability of a certain linear discrete system. A $\mathtt{MATLAB}$ code checks the null controllability up to a given finite time T (the order of the differential operator), and provides all null control sequences (polynomial coefficients of the differential operator) up to a given maximum degree m. This is the first paper that connects Stein’s method with computational algebra to find Stein operators for highly complex probability distributions, such as $H_{20}(X_1)$, where $H_p$ is the p-th Hermite polynomial. Some examples of Stein operators for $H_p(X_1)$, $p=3,4,5,6$, are gathered in the Appendix and many other examples are given in the Supplementary Information.

In many applications, we have access to the complete dataset but are only interested in the prediction of a particular region of predictor variables. A standard approach is to find the globally best modeling method from a set of candidate methods. However, it is perhaps rare in reality that one candidate method is uniformly better than the others. A natural approach for this scenario is to apply a weighted $L_2$ loss in performance assessment to reflect the region-specific interest. We propose a targeted cross-validation (TCV) to select models or procedures based on a general weighted $L_2$ loss. We show that the TCV is consistent in selecting the best performing candidate under the weighted $L_2$ loss. Experimental studies are used to demonstrate the use of TCV and its potential advantage over the global CV or the approach of using only local data for modeling a local region.

Previous investigations on CV have relied on the condition that when the sample size is large enough, the ranking of two candidates stays the same. However, in many applications with the setup of changing data-generating processes or highly adaptive modeling methods, the relative performance of the methods is not static as the sample size varies. Even with a fixed data-generating process, it is possible that the ranking of two methods switches infinitely many times. In this work, we broaden the concept of the selection consistency by allowing the best candidate to switch as the sample size varies, and then establish the consistency of the TCV. This flexible framework can be applied to high-dimensional and complex machine learning scenarios where the relative performances of modeling procedures are dynamic.

We study a Bayesian approach to the inference of an event time distribution in the current status model where observation times are supported on a grid of potentially unknown sparsity and multiple subjects share the same observation time. The model leads to a very simple likelihood, but statistical inferences are non-trivial due to the unknown sparsity of the grid. In particular, for an inference based on the maximum likelihood estimator, one needs to estimate the density of the event time distribution which is challenging because the event time is not directly observed. We consider Bayes procedures with a Dirichlet prior on the event time distribution. With this prior, the Bayes estimator and credible sets can be easily computed via a Gibbs sampler algorithm. Our main contribution is to provide thorough investigation of frequentist’s properties of the posterior distribution. Specifically, it is shown that the posterior convergence rate is adaptive to the unknown sparsity of the grid. If the grid is sufficiently sparse, we further prove the Bernstein–von Mises theorem which guarantees frequentist’s validity of Bayesian credible sets. A numerical study is also conducted for illustration.

Although the operator (spectral) norm is one of the most widely used metrics for covariance estimation, comparatively little is known about the fluctuations of error in this norm. To be specific, let $\hat{\mathrm{\Sigma }}$ denote the sample covariance matrix of n i.i.d. observations in ${\mathbb{R}}^p$ that arise from a population matrix Σ, and let $T_n=\sqrt{n}\Vert \hat{\mathrm{\Sigma }}-\mathrm{\Sigma }{\Vert }_{\text{op}}$. In the setting where the eigenvalues of Σ have a decay profile of the form ${\mathit{\lambda }}_j(\mathrm{\Sigma })\asymp j^{-2\mathit{\beta }}$, we analyze how well the bootstrap can approximate the distribution of $T_n$. Our main result shows that up to factors of $log(n)$, the bootstrap can approximate the distribution of $T_n$ with respect to the Kolmogorov metric at the rate of $n^{-\frac{\mathit{\beta }-1∕2}{6\mathit{\beta }+4}}$, which does not depend on the ambient dimension p. In addition, we offer a supporting result of independent interest that establishes a high-probability upper bound for $T_n$ based on flexible moment assumptions. More generally, we discuss the consequences of our work beyond covariance matrices, and show how the bootstrap can be used to estimate the errors of sketching algorithms in randomized numerical linear algebra (RandNLA). An illustration of these ideas is also provided with a climate data example.

Modelling a large collection of functional time series arises in a broad spectral of real applications. Under such a scenario, not only the number of functional variables can be diverging with, or even larger than the number of temporally dependent functional observations, but each function itself is an infinite-dimensional object, posing a challenging task. In this paper, we propose a three-step procedure to estimate high-dimensional functional time series models. To provide theoretical guarantees for the three-step procedure, we focus on multivariate stationary processes and propose a novel functional stability measure based on their spectral properties. Such stability measure facilitates the development of some useful concentration bounds on sample (auto)covariance functions, which serve as a fundamental tool for further convergence analysis in high-dimensional settings. As functional principal component analysis (FPCA) is one of the key dimension reduction techniques in the first step, we also investigate the non-asymptotic properties of the relevant estimated terms under a FPCA framework. To illustrate with an important application, we consider vector functional autoregressive models and develop a regularization approach to estimate autoregressive coefficient functions under the sparsity constraint. Using our derived non-asymptotic results, we investigate convergence properties of the regularized estimate under high-dimensional scaling. Finally, the finite-sample performance of the proposed method is examined through both simulations and a public financial dataset.

We study the problem of predicting as well as the best linear predictor in a bounded Euclidean ball with respect to the squared loss. When only boundedness of the data generating distribution is assumed, we establish that the least squares estimator constrained to a bounded Euclidean ball does not attain the classical $O(d∕n)$ excess risk rate, where d is the dimension of the covariates and n is the number of samples. In particular, we construct a bounded distribution such that the constrained least squares estimator incurs an excess risk of order $\mathrm{\Omega }(d^{3∕2}∕n)$ hence refuting a recent conjecture of Ohad Shamir [JMLR 2015]. In contrast, we observe that non-linear predictors can achieve the optimal rate $O(d∕n)$ with no assumptions on the distribution of the covariates. We discuss additional distributional assumptions sufficient to guarantee an $O(d∕n)$ excess risk rate for the least squares estimator. Among them are certain moment equivalence assumptions often used in the robust statistics literature. While such assumptions are central in the analysis of unbounded and heavy-tailed settings, our work indicates that in some cases, they also rule out unfavorable bounded distributions.

We study the concept of quadratic variation of a continuous path along a sequence of partitions and its dependence with respect to the choice of the partition sequence. We introduce the concept of quadratic roughness of a path along a partition sequence and show that for Hölder-continuous paths satisfying this roughness condition, the quadratic variation along balanced partitions is invariant with respect to the choice of the partition sequence. Typical paths of Brownian motion are shown to satisfy this quadratic roughness property almost-surely along any partition with a required step size condition. Using these results we derive a formulation of the pathwise Föllmer-Itô calculus which is invariant with respect to the partition sequence. We also derive an invarience of local time under quadratic roughness.

We establish a Chung-type law of the iterated logarithm and the exact local and uniform moduli of continuity for a large class of anisotropic Gaussian random fields with a harmonizable-type integral representation and the property of strong local nondeterminism. Compared with the existing results in the literature, our results do not require the assumption of stationary increments and provide more precise upper and lower bounds for the limiting constants. The results are applicable to the solutions of a class of linear stochastic partial differential equations driven by a fractional-colored Gaussian noise, including the stochastic heat equation.

We consider situations where the applicability of sequential Monte Carlo particle filters is compromised due to the expensive evaluation of the particle weights. To alleviate this problem, we propose a new particle filter algorithm based on the multilevel approach. We show that the resulting multilevel bootstrap particle filter (MLBPF) retains the strong law of large numbers as well as the central limit theorem of classical particle filters under mild conditions. Our numerical experiments demonstrate up to 85% reduction in computation time compared to the classical bootstrap particle filter, in certain settings. While it should be acknowledged that this reduction is highly application dependent, and a similar gain should not be expected for all applications across the board, we believe that this substantial improvement in certain settings makes MLBPF an important addition to the family of sequential Monte Carlo methods.

We obtain explicit criteria for both exponential ergodicity and strong ergodicity for one-dimensional time-changed symmetric stable processes with $\mathit{\alpha }\in (1,2)$. Explicit lower bounds for ergodic convergence rates are given.

A continuous-time particle system on the real line satisfying the branching property and an exponential integrability condition is called a branching Lévy process, and its law is characterized by a triplet $({\mathit{\sigma }}^2,a,\mathrm{\Lambda })$. We obtain a necessary and sufficient condition for the convergence of the derivative martingale of such a process to a non-trivial limit in terms of $({\mathit{\sigma }}^2,a,\mathrm{\Lambda })$. This extends previously known results on branching Brownian motions and branching random walks. To obtain this result, we rely on the spinal decomposition and establish a novel zero-one law on the perpetual integrals of centred Lévy processes conditioned to stay positive.

We consider the Riemannian random wave model of Gaussian linear combinations of Laplace eigenfunctions on a general compact Riemannian manifold. With probability one with respect to the Gaussian coefficients, we establish that, both for large and small band models, the process properly rescaled and evaluated at an independently and uniformly chosen point X on the manifold, converges in distribution under the sole randomness of X towards an universal Gaussian field as the frequency tends to infinity. This result extends the celebrated central limit Theorem of Salem–Zygmund for trigonometric polynomials series to the more general framework of compact Riemannian manifolds. We then deduce from the above convergence the almost-sure asymptotics of the nodal volume associated with the random wave. To the best of our knowledge, in the real Riemannian case, these asymptotics were only known in expectation and not in the almost sure sense due to the lack of sufficiently accurate variance estimates. This in particular addresses a question of S. Zelditch regarding the almost sure equidistribution of nodal volume.

We study randomized designs that minimize the asymptotic variance of a debiased lasso estimator when a large pool of unlabeled data is available but measuring the corresponding responses is costly. The optimal sampling distribution arises as the solution of a semidefinite program. The improvements in efficiency that result from these optimal designs are demonstrated via simulation experiments.

We introduce a general approach which unifies some previous attempts for modeling the dynamic of multivariate time series or for regression analysis when the data are of mixed type (binary/count/continuous). Our approach is quite flexible since conditionally on past values, each coordinate at time t can have a distribution compatible with a standard univariate time series model such as GARCH, ARMA, INGARCH or logistic models whereas past values of the other coordinates play the role of exogenous covariates in the dynamic. The simultaneous dependence in the multivariate time series can be modeled with a copula. Additional exogenous covariates are also allowed in the dynamic. We first study some usual stability properties of these models and then show that autoregressive parameters can be consistently estimated equation-by-equation using a pseudo-maximum likelihood method, leading to a fast implementation even when the number of time series is large. Moreover, we prove consistency results when a parametric copula model is fitted to the time series and in the case of Gaussian copulas, we show that the likelihood estimator of the correlation matrix is strongly consistent. We carefully check all our assumptions for two prototypical examples: a GARCH/INGARCH model and logistic/log-linear INGARCH model. Our results are illustrated with numerical experiments as well as two real data sets.

We analyze the bootstrap percolation process on the stochastic block model (SBM), a natural extension of the Erdős–Rényi random graph that incorporates the community structure observed in many real systems. In the SBM, nodes are partitioned into two subsets, which represent different communities, and pairs of nodes are independently connected with a probability that depends on the communities they belong to. Under mild assumptions on the system parameters, we prove the existence of a sharp phase transition for the final number of active nodes and characterize the sub-critical and the super-critical regimes in terms of the number of initially active nodes, which are selected uniformly at random in each community.

We formulate a trajectorial version of the relative entropy dissipation identity for McKean–Vlasov diffusions, extending recent results which apply to non-interacting diffusions. Our stochastic analysis approach is based on time-reversal of diffusions and Lions’ differential calculus over Wasserstein space. It allows us to compute explicitly the rate of relative entropy dissipation along every trajectory of the underlying diffusion via the semimartingale decomposition of the corresponding relative entropy process. As a first application, we obtain a new interpretation of the gradient flow structure for the granular media equation, generalizing a formulation developed recently for the linear Fokker–Planck equation. Secondly, we show how the trajectorial approach leads to a new derivation of the HWBI inequality, which relates relative entropy (H), Wasserstein distance (W), barycenter (B) and Fisher information (I).

For $d\ge 3$, we prove that time-inhomogeneous stochastic differential equations driven by additive noises with drifts in critical Lebesgue space $L^q([0,T];L^p({\mathbb{R}}^d))$, where $(p,q)\in (d,\mathrm{\infty }]\times [2,\mathrm{\infty })$ and $d∕p+2∕q=1$, or $(p,q)=(d,\mathrm{\infty })$ and $divb\in L^{\mathrm{\infty }}([0,T];L^{d∕2+\mathit{\epsilon }}({\mathbb{R}}^d))$, are well-posed. The weak uniqueness is obtained by solving corresponding Kolmogorov backward equations in some second-order Sobolev spaces, which is analytically interesting in itself.

We introduce the notion of linear multifractional stable sheets in the broad sense (LMSS) with $\mathit{\alpha }\in (0,2]$, to include both linear multifractional Brownian sheets ($\mathit{\alpha }=2$) and linear multifractional stable sheets ($\mathit{\alpha }<2$). The purpose of the present paper is to study the existence and joint continuity of the local times of LMSS, and also the local Hölder condition of the local times in the set variable. Among the main results of this paper, Theorem 2.4 provides a sufficient and necessary condition for the existence of local times of LMSS; Theorem 3.1 shows a sufficient condition for the joint continuity of local times; and Theorem 4.1 proves a sharp local Hölder condition for the local times in the set variable. All these theorems improve significantly the existing results for the local times of multifractional Brownian sheets and linear multifractional stable sheets in the literature.

Non-homogeneous Poisson processes are used in a wide range of scientific disciplines, ranging from the environmental sciences to the health sciences. Often, the central object of interest in a point process is the underlying intensity function. Here, we present a general model for the intensity function of a non-homogeneous Poisson process using measure transport. The model is built from a flexible bijective mapping that maps from the underlying intensity function of interest to a simpler reference intensity function. We enforce bijectivity by modeling the map as a composition of multiple bijective maps that have increasing triangular structure, and show that the model exhibits an important approximation property. Estimation of the flexible mapping is accomplished within an optimization framework, wherein computations are efficiently done using tools originally designed to facilitate deep learning, and a graphics processing unit. Point process simulation and uncertainty quantification are straightforward to do with the proposed model. We demonstrate the potential benefits of our proposed method over conventional approaches to intensity modeling through various simulation studies. We also illustrate the use of our model on a real data set containing the locations of seismic events near Fiji since 1964.

The main scope of this paper is to give some explicit classes of examples of $L^1$-optimal couplings. Optimal transportation w.r.t. the Kantorovich metric ${\ell }_1$ (resp. the Wasserstein metric $W_1$) between two absolutely continuous measures is known since the basic papers of Kantorovich and Rubinstein (Dokl. Akad. Nauk SSSR 115 (1957) 1058–1061) and Sudakov (Proc. Steklov Inst. Math. 141 (1979) 1–178) to occur on rays induced by a decomposition of the basic space (and more generally to higher dimensional decompositions in the case of general measures) induced by the corresponding dual potentials. Several papers have given this kind of structural result and established existence and uniqueness of solutions in varying generality. Since the dual problems pose typically too strong challenges to be solved in explicit form, these structural results have so far been applied for the solution of few particular instances.

First, we give a self-contained review of some basic optimal coupling results and we propose and investigate in particular some basic principles for the construction of $L^1$-optimal couplings given by a reduction principle and some usable forms of the decomposition method. This reduction principle, together with symmetry properties of the reduced measures, gives a hint to the decomposition of the space into sectors and via the non crossing property of optimal transport leads to the choice of transportation rays. The optimality of the induced transports is then a consequence of the characterization results of optimal couplings.

Then, we apply these principles to determine in explicit form $L^1$-optimal couplings for several classes of examples of elliptical distributions. In particular, we give for the first time a general construction of $L^1$-optimal couplings between two bivariate Gaussian distributions. We also discuss optimality of special constructions like shifts and scalings, and provide an extended class of dual functionals allowing for the closed-form computation of the ${\ell }_1$-metric or of accurate lower bounds of it in a variety of examples.