In recent years, rough volatility models have gained considerable attention in quantitative finance. In this paradigm, the stochastic volatility of the price of an asset has quantitative properties similar to that of a fractional Brownian motion with small Hurst index $\mathit{H}<1/2$. In this work, we provide the first rigorous statistical analysis of the problem of estimating H from historical observations of the underlying asset. We establish minimax lower bounds and design optimal procedures based on adaptive estimation of quadratic functionals based on wavelets. We prove in particular that the optimal rate of convergence for estimating H based on price observations at n time points is of order ${\mathit{n}}^{-1/(4\mathit{H}\mathbf{+}2)}$ as n grows to infinity, extending results that were known only for $\mathit{H}>1/2$. Our study positively answers the question whether H can be inferred, although it is the regularity of a latent process (the volatility); in rough models, when H is close to 0, we even obtain an accuracy comparable to usual $\sqrt{\mathit{n}}$-consistent regular statistical models.

The coefficients in a second order parabolic linear stochastic partial differential equation (SPDE) are estimated from multiple spatially localised measurements. Assuming that the spatial resolution tends to zero and the number of measurements is nondecreasing, the rate of convergence for each coefficient depends on its differential order and is faster for higher order coefficients. Based on an explicit analysis of the reproducing kernel Hilbert space of a general stochastic evolution equation, a Gaussian lower bound scheme is introduced. As a result, minimax optimality of the rates as well as sufficient and necessary conditions for consistent estimation are established.

We consider the twin problems of estimating the effective rank and the Schatten norms $\Vert \mathbf{A}{\Vert }_{\mathit{s}}$ of a rectangular $\mathit{p}\times \mathit{q}$ matrix A from noisy observations. When s is an even integer, we introduce a polynomial-time estimator of $\Vert \mathbf{A}{\Vert }_{\mathit{s}}$ that achieves the minimax rate ${(\mathit{p}\mathit{q})}^{1/4}$. Interestingly, this optimal rate does not depend on the underlying rank of the matrix A. When s is not an even integer, the optimal rate is much slower. A simple thresholding estimator of the singular values achieves the rate $(\mathit{q}\wedge \mathit{p}){(\mathit{p}\mathit{q})}^{1/4}$, which turns out to be optimal up to a logarithmic multiplicative term. The tight minimax rate is achieved by a more involved polynomial approximation method. This allows us to build estimators for a class of effective rank indices. As a byproduct, we also characterize the minimax rate for estimating the sequence of singular values of a matrix.

While batching methods have been widely used in simulation and statistics, their higher-order coverage behaviors and relative advantages in this regard remain open. We develop techniques to obtain higher-order coverage errors for batching methods by building Edgeworth-type expansions on t-statistics. The coefficients in these expansions are intricate analytically, but we provide algorithms to estimate the coefficients of the ${\mathit{n}}^{-1}$ error terms via Monte Carlo simulation. We provide insights on the effect of the number of batches on the coverage error, where we demonstrate generally nonmonotonic relations. We also compare different batching methods both theoretically and numerically, and argue that none of the methods is uniformly better than others in terms of coverage. However, when the number of batches is large, sectioned jackknife has the best coverage among all.

The main computational challenge in Bayesian inference is to compute integrals against a high-dimensional posterior distribution. In the past decades, variational inference (VI) has emerged as a tractable approximation to these integrals, and a viable alternative to the more established paradigm of Markov chain Monte Carlo. However, little is known about the approximation accuracy of VI. In this work, we bound the TV error and the mean and covariance approximation error of Gaussian VI in terms of dimension and sample size. Our error analysis relies on a Hermite series expansion of the log posterior whose first terms are precisely cancelled out by the first order optimality conditions associated to the Gaussian VI optimization problem.

We study worst-case-growth-rate-optimal (GROW) $\mathit{\text{e}}$-statistics for hypothesis testing between two group models. It is known that under a mild condition on the action of the underlying group G on the data, there exists a maximally invariant statistic. We show that among all $\mathit{\text{e}}$-statistics, invariant or not, the likelihood ratio of the maximally invariant statistic is GROW, both in the absolute and in the relative sense, and that an anytime-valid test can be based on it. The GROW $\mathit{\text{e}}$-statistic is equal to a Bayes factor with a right Haar prior on G. Our treatment avoids nonuniqueness issues that sometimes arise for such priors in Bayesian contexts. A crucial assumption on the group G is its amenability, a well-known group-theoretical condition, which holds, for instance, in scale-location families. Our results also apply to finite-dimensional linear regression.

We propose a new Bayesian strategy for adaptation to smoothness in nonparametric models based on heavy-tailed series priors. We illustrate it in a variety of settings, showing in particular that the corresponding Bayesian posterior distributions achieve adaptive rates of contraction in the minimax sense (up to logarithmic factors) without the need to sample hyperparameters. Unlike many existing procedures, where a form of direct model (or estimator) selection is performed, the method can be seen as performing a soft selection through the prior tail. In Gaussian regression, such heavy-tailed priors are shown to lead to (near-)optimal simultaneous adaptation both in the ${\mathit{L}}^2$- and ${\mathit{L}}^{\infty }$-sense. Results are also derived for linear inverse problems, for anisotropic Besov classes, and for certain losses in more general models through the use of tempered posterior distributions. We present numerical simulations corroborating the theory.

We consider the problem of estimating the factors of a low-rank $\mathit{n}\times \mathit{d}$ matrix, when this is corrupted by additive Gaussian noise. A special example of our setting corresponds to clustering mixtures of Gaussians with equal (known) covariances. Simple spectral methods do not take into account the distribution of the entries of these factors and are therefore often suboptimal. Here, we characterize the asymptotics of the minimum estimation error under the assumption that the distribution of the entries is known to the statistician.

Our results apply to the high-dimensional regime $\mathit{n},\mathit{d}\to \infty$ and $\mathit{d}/\mathit{n}\to \infty$ (or $\mathit{d}/\mathit{n}\to 0$) and generalize earlier work that focused on the proportional asymptotics $\mathit{n},\mathit{d}\to \infty$, $\mathit{d}/\mathit{n}\to \mathit{\delta }\in (0,\infty )$. We outline an interesting signal strength regime in which $\mathit{d}/\mathit{n}\to \infty$ and partial recovery is possible for the left singular vectors while impossible for the right singular vectors.

We illustrate the general theory by deriving consequences for Gaussian mixture clustering and carrying out a numerical study on genomics data.

In their seminal work, Polyak and Juditsky showed that stochastic approximation algorithms for solving smooth equations enjoy a central limit theorem. Moreover, it has since been argued that the asymptotic covariance of the method is best possible among any estimation procedure in a local minimax sense of Hájek and Le Cam. A long-standing open question in this line of work is whether similar guarantees hold for important nonsmooth problems, such as stochastic nonlinear programming or stochastic variational inequalities. In this work, we show that this is indeed the case.

Westling and Carone (Ann. Statist. 48 (2020) 1001–1024) proposed a framework for studying the large sample distributional properties of generalized Grenander-type estimators, a versatile class of nonparametric estimators of monotone functions. The limiting distribution of those estimators is representable as the left derivative of the greatest convex minorant of a Gaussian process whose monomial mean can be of unknown order (when the degree of flatness of the function of interest is unknown). The standard nonparametric bootstrap is unable to consistently approximate the large sample distribution of the generalized Grenander-type estimators even if the monomial order of the mean is known, making statistical inference a challenging endeavour in applications. To address this inferential problem, we present a bootstrap-assisted inference procedure for generalized Grenander-type estimators. The procedure relies on a carefully crafted, yet automatic, transformation of the estimator. Moreover, our proposed method can be made “flatness robust” in the sense that it can be made adaptive to the (possibly unknown) degree of flatness of the function of interest. The method requires only the consistent estimation of a single scalar quantity, for which we propose an automatic procedure based on numerical derivative estimation and the generalized jackknife. Under random sampling, our inference method can be implemented using a computationally attractive exchangeable bootstrap procedure. We illustrate our methods with examples and we also provide a small simulation study. The development of formal results is made possible by some technical results that may be of independent interest.

We present estimators for smooth Hilbert-valued parameters, where smoothness is characterized by a pathwise differentiability condition. When the parameter space is a reproducing kernel Hilbert space, we provide a means to obtain efficient, root-n rate estimators and corresponding confidence sets. These estimators correspond to generalizations of cross-fitted one-step estimators based on Hilbert-valued efficient influence functions. We give theoretical guarantees even when arbitrary estimators of nuisance functions are used, including those based on machine learning techniques. We show that these results naturally extend to Hilbert spaces that lack a reproducing kernel, as long as the parameter has an efficient influence function. However, we also uncover the unfortunate fact that, when there is no reproducing kernel, many interesting parameters fail to have an efficient influence function, even though they are pathwise differentiable. To handle these cases, we propose a regularized one-step estimator and associated confidence sets. We also show that pathwise differentiability, which is a central requirement of our approach, holds in many cases. Specifically, we provide multiple examples of pathwise differentiable parameters and develop corresponding estimators and confidence sets. Among these examples, four are particularly relevant to ongoing research by the causal inference community: the counterfactual density function, dose-response function, conditional average treatment effect function, and counterfactual kernel mean embedding.

This work investigates multiple testing by considering minimax separation rates in the sparse sequence model, when the testing risk is measured as the sum FDR+FNR (False Discovery Rate plus False Negative Rate). First, using the popular beta-min separation condition, with all nonzero signals separated from 0 by at least some amount, we determine the sharp minimax testing risk asymptotically and thereby explicitly describe the transition from “achievable multiple testing with vanishing risk” to “impossible multiple testing.” Adaptive multiple testing procedures achieving the corresponding optimal boundary are provided: the Benjamini–Hochberg procedure with a properly tuned level, and an empirical Bayes ℓ-value (‘local FDR’) procedure. We prove that the FDR and FNR make nonsymmetric contributions to the testing risk for most optimal procedures, the FNR part being dominant at the boundary. The multiple testing hardness is then investigated for classes of arbitrary sparse signals. A number of extensions, including results for classification losses and convergence rates in the case of large signals, are also investigated.

The vast majority of literature on evaluating the significance of a treatment effect based on observational data has been confined to discrete treatments. These methods are not applicable to drawing inference for a continuous treatment, which arises in many important applications. To adjust for confounders when evaluating a continuous treatment, existing inference methods often rely on discretizing the treatment or using (possibly misspecified) parametric models for the effect curve. Recently, Kennedy et al. (J. R. Stat. Soc. Ser. B. Stat. Methodol. 79 (2017) 1229–1245) proposed nonparametric doubly robust estimation for a continuous treatment effect in observational studies. However, inference for the continuous treatment effect is a harder problem. To the best of our knowledge, a completely nonparametric doubly robust approach for inference in this setting is not yet available. We develop such a nonparametric doubly robust procedure in this paper for making inference on the continuous treatment effect curve. Using empirical process techniques for local U- and V-processes, we establish the test statistic’s asymptotic distribution. Furthermore, we propose a wild bootstrap procedure for implementing the test in practice. In addition, we define a version of the test procedure based on sample splitting. We illustrate the new method(s) via simulations and a study of a constructed dataset relating the effect of nurse staffing hours on hospital performance. We implement our doubly robust dose response test in the R package DRDRtest on CRAN.

The Gromov–Wasserstein (GW) distance, rooted in optimal transport (OT) theory, quantifies dissimilarity between metric measure spaces and provides a framework for aligning heterogeneous datasets. While computational aspects of the GW problem have been widely studied, a duality theory and fundamental statistical questions concerning empirical convergence rates remained obscure. This work closes these gaps for the quadratic GW distance over Euclidean spaces of different dimensions ${\mathit{d}}_{\mathit{x}}$ and ${\mathit{d}}_{\mathit{y}}$. We treat both the standard and the entropically regularized GW distance, and derive dual forms that represent them in terms of the well-understood OT and entropic OT (EOT) problems, respectively. This enables employing proof techniques from statistical OT based on regularity analysis of dual potentials and empirical process theory, using which we establish the first GW empirical convergence rates. The derived two-sample rates are ${\mathit{n}}^{-2/max\{min\{{\mathit{d}}_{\mathit{x}},{\mathit{d}}_{\mathit{y}}\},4\}}$ (up to a log factor when $min\{{\mathit{d}}_{\mathit{x}},{\mathit{d}}_{\mathit{y}}\}=4$) for standard GW and ${\mathit{n}}^{-1/2}$ for entropic GW (EGW), which matches the corresponding rates for standard and entropic OT. The parametric rate for EGW is evidently optimal, while for standard GW we provide matching lower bounds, which establish sharpness of the derived rates. We also study stability of EGW in the entropic regularization parameter and prove approximation and continuity results for the cost and optimal couplings. Lastly, the duality is leveraged to shed new light on the open problem of the one-dimensional GW distance between uniform distributions on n points, illuminating why the identity and anti-identity permutations may not be optimal. Our results serve as a first step towards a comprehensive statistical theory as well as computational advancements for GW distances, based on the discovered dual formulations.

This paper proposes a novel methodology for the online detection of changepoints in the factor structure of large matrix time series. Our approach is based on the well-known fact that, in the presence of a changepoint, the number of spiked eigenvalues in the second moment matrix of the data increases (e.g., in the presence of a change in the loadings, or if a new factor emerges). Based on this, we propose two families of procedures—one based on the fluctuations of partial sums, and one based on extreme value theory—to monitor whether the first nonspiked eigenvalue diverges after a point in time in the monitoring horizon, thereby indicating the presence of a changepoint. Our procedure is based only on rates; at each point in time, we randomise the estimated eigenvalue, thus obtaining a normally distributed sequence which is i.i.d. with mean zero under the null of no break, whereas it diverges to positive infinity in the presence of a changepoint. We base our monitoring procedures on such sequence. Extensive simulation studies and empirical analysis justify the theory. An R package implementing the procedure is available on CRAN. ( https://cran.r-project.org/web/packages/OLCPM/index.html.)

Privacy-preserving data analysis has become more prevalent in recent years. In this study, we propose a distributed group differentially private Majority Vote mechanism, for the sign selection problem in a distributed setup. To achieve this, we apply the iterative peeling to the stability function and use the exponential mechanism to recover the signs. For enhanced applicability, we study the private sign selection for mean estimation and linear regression problems, in distributed systems. Our method recovers the support and signs with the optimal signal-to-noise ratio as in the nonprivate scenario, which is better than contemporary works of private variable selections. Moreover, the sign selection consistency is justified by theoretical guarantees. Simulation studies are conducted to demonstrate the effectiveness of the proposed method.

We study the multivariate deconvolution problem of recovering the distribution of a signal from independent and identically distributed observations additively contaminated with random errors having known distribution. For errors with ordinary smooth distribution, we recast the multidimensional problem as a one-dimensional problem leveraging the equivalence between the ${\mathit{L}}^1$-Wasserstein and the max-sliced ${\mathit{L}}^1$-Wasserstein metrics and derive an inversion inequality relating the ${\mathit{L}}^1$-Wasserstein distance between two distributions of the signal to the ${\mathit{L}}^1$-distance between the corresponding mixture densities of the observations. This smoothing inequality outperforms existing inversion inequalities. We apply it to derive ${\mathit{L}}^1$-Wasserstein rates of convergence for the distribution of the signal. As an application to the Bayesian framework, we consider ${\mathit{L}}^1$-Wasserstein deconvolution with the Laplace noise in dimension one using a Dirichlet process mixture of normal densities as a prior measure for the mixing distribution. We construct an adaptive approximation of the sampling density by convolving the Laplace density with a well-chosen mixture of Gaussian densities and show that the posterior measure contracts at a nearly minimax-optimal rate, up to a log-factor, in the ${\mathit{L}}^1$-distance. The rate automatically adapts to the unknown Sobolev regularity of the mixing density, thus leading to a new Bayesian adaptive estimation procedure over the full scale of regularity levels. We illustrate the utility of the inversion inequality also in a frequentist setting by showing that a minimum distance estimator attains the minimax convergence rates for ${\mathit{L}}^1$-Wasserstein deconvolution in any dimension $\mathit{d}\ge 1$, lower bounds being derived here.

We consider detecting and estimating breaks in heterogenous mean functions of high-dimensional functional time series which are allowed to be cross-sectionally correlated. A new test statistic combining the functional CUSUM statistic and power enhancement component is proposed with asymptotic null distribution comparable to the conventional CUSUM theory derived for a single functional time series. In particular, the extra power enhancement component enlarges the region where the proposed test has power, and results in stable power performance when breaks are sparse in the alternative hypothesis. Furthermore, we impose a latent group structure on the subjects with heterogenous break points and introduce an easy-to-implement clustering algorithm with an information criterion to consistently estimate the unknown group number and membership. The estimated group structure improves the convergence property of the break point estimate. Monte Carlo simulation studies and empirical applications show that the proposed estimation and testing techniques have satisfactory performance in finite samples.

Smooth backfitting has been proposed and proved as a powerful nonparametric estimation technique for additive regression models in various settings. Existing studies are restricted to cases with a moderate number of covariates and are not directly applicable to high dimensional settings. In this paper, we develop new kernel estimators based on the idea of smooth backfitting for high dimensional additive models. We introduce a novel penalization scheme, combining the idea of functional Lasso with the smooth backfitting technique. We investigate the theoretical properties of the functional Lasso smooth backfitting estimation. For the implementation of the proposed method, we devise a simple iterative algorithm where the iteration is defined by a truncated projection operator. The algorithm has only an additional thresholding operator over the projection-based iteration of the smooth backfitting algorithm. We further present a debiased version of the proposed estimator with implementation details, and investigate its theoretical properties for statistical inference. We demonstrate the finite sample performance of the methods via simulation and real data analysis.

Gaussian mixture models form a flexible and expressive parametric family of distributions that has found a variety of applications. Unfortunately, fitting these models to data is a notoriously hard problem from a computational perspective. Currently, only moment-based methods enjoy theoretical guarantees while likelihood-based methods are dominated by heuristics such as Expectation-Maximization that are known to fail in simple examples. In this work, we propose a new algorithm to compute the nonparametric maximum likelihood estimator (NPMLE) in a Gaussian mixture model. Our method is based on gradient descent over the space of probability measures equipped with the Wasserstein–Fisher–Rao geometry for which we establish convergence guarantees. In practice, it can be approximated using an interacting particle system where the weight and location of particles are updated alternately. We conduct extensive numerical experiments to confirm the effectiveness of the proposed algorithm compared not only to classical benchmarks but also to similar gradient descent algorithms with respect to simpler geometries. In particular, these simulations illustrate the benefit of updating both weight and location of the interacting particles.

Statistical machine learning methods often face the challenge of limited data available from the population of interest. One remedy is to leverage data from auxiliary source populations, which share some conditional distributions or are linked in other ways with the target domain. Techniques leveraging such dataset shift conditions are known as domain adaptation or transfer learning. Despite extensive literature on dataset shift, limited works address how to efficiently use the auxiliary populations to improve the accuracy of risk evaluation for a given machine learning task in the target population.

In this paper, we study the general problem of efficiently estimating target population risk under various dataset shift conditions, leveraging semiparametric efficiency theory. We consider a general class of dataset shift conditions, which includes three popular conditions—covariate, label and concept shift—as special cases. We allow for partially nonoverlapping support between the source and target populations. We develop efficient and multiply robust estimators along with a straightforward specification test of these dataset shift conditions. We also derive efficiency bounds for two other dataset shift conditions, posterior drift and location-scale shift. Simulation studies support the efficiency gains due to leveraging plausible dataset shift conditions.

We consider a nonlinear Bayesian data assimilation model for the periodic two-dimensional Navier–Stokes equations with initial condition modelled by a Gaussian process prior. We show that if the system is updated with sufficiently many discrete noisy measurements of the velocity field, then the posterior distribution eventually concentrates near the ground truth solution of the time evolution equation, and in particular that the initial condition is recovered consistently by the posterior mean vector field. We further show that the convergence rate can in general not be faster than inverse logarithmic in sample size, but describe specific conditions on the initial conditions when faster rates are possible. In the proofs, we provide an explicit quantitative estimate for backward uniqueness of solutions of the two-dimensional Navier–Stokes equations.

This paper investigates the stability of deep ReLU neural networks for nonparametric regression under the assumption that the noise has only a finite pth moment. We unveil how the optimal rate of convergence depends on p, the degree of smoothness and the intrinsic dimension in a class of nonparametric regression functions with hierarchical composition structure when both the adaptive Huber loss and deep ReLU neural networks are used. This optimal rate of convergence cannot be obtained by the ordinary least squares but can be achieved by the Huber loss with a properly chosen parameter that adapts to the sample size, smoothness, and moment parameters. A concentration inequality for the adaptive Huber ReLU neural network estimators with allowable optimization errors is also derived. To establish a matching lower bound within the class of neural network estimators using the Huber loss, we employ a different strategy from the traditional route: constructing a deep ReLU network estimator that has a better empirical loss than the true function and the difference between these two functions furnishes a low bound. This step is related to the Huberization bias, yet more critically to the approximability of deep ReLU networks. As a result, we also contribute some new results on the approximation theory of deep ReLU neural networks.