In fitting linear regression models for functional data, a complicating factor with regressors as random curves is that regression estimators have complex distributions, due to issues in bias and scaling. Bias arises because the target slope function is infinite-dimensional, while finite-sample estimators necessarily involve truncations. To approximate sampling distributions, we develop a residual bootstrap method. Despite the parametric regression problem, the bootstrap for functional data requires a development that resembles resampling for nonparametric regression with multivariate regressors. Essentially, original- and bootstrap-data estimators require coordination in the truncation levels to remove bias (akin to tuning parameter choices). The resulting bootstrap has wide applicability for constructing both confidence and prediction regions at target regressor points, and with coverage properties even holding conditionally on data regressors; the method also extends to simultaneous regions. Establishment of the bootstrap further involves generalizing, refining, and correcting a foundational central limit theorem for functional linear regression. Numerical studies verify our theory, showing that the bootstrap performs better than normal approximations, and also suggest a rule of thumb for setting the truncation levels. The bootstrap method is illustrated with an application to wheat spectrum data.

In this work, we generalise the stochastic local time space integration introduced in (Potential Anal. 13 (2000) 303–328) to the case of Brownian sheet. This allows us to prove a generalised two-parameter Itô formula and derive Davie type inequalities for the Brownian sheet. Such estimates are useful to obtain regularity bounds for some averaging type operators along Brownian sheet curves.

This paper studies the problem of model selection in a large class of causal time series models that includes ARMA or AR (∞) processes as well as GARCH or ARCH (∞), APARCH, ARMA-GARCH - and many other processes. First, we study the asymptotic behavior of the ideal penalty that minimizes the risk defined from a quasi-likelihood estimation among a finite family of models containing the true model. We then establish general conditions on the penalty term to obtain properties of consistency and efficiency. In particular, we prove that consistent model selection criteria outperform the classical AIC criterion in terms of efficiency. Finally, we derive the usual BIC criterion from a Bayesian approach and, retaining all second-order terms of the Laplace approximation, a data-driven criterion, which we call KC’. Monte Carlo experiments illustrate the asymptotic results obtained and show that the KC’ criterion performs better than the AIC and BIC criteria in terms of consistency and efficiency.

The main result of this article is that we obtain an elementwise error bound for the Fused Lasso estimator for any general convex loss function ρ. We then focus on the special cases when either ρ is the square loss function (for mean regression) or is the quantile loss function (for quantile regression) for which we derive new pointwise error bounds. Even though error bounds for the usual Fused Lasso estimator and its quantile version have been studied before; our bound appears to be new. This is because all previous works bound a global loss function like the sum of squared error, or a sum of huber losses in the case of quantile regression in Padilla and Chatterjee (Biometrika 109 (2022) 751–768). Clearly, element wise bounds are stronger than global loss error bounds as it reveals how the loss behaves locally at each point. Our element wise error bound also has a clean and explicit dependence on the tuning parameter λ which informs the user of a good choice of λ. In addition, our bound is nonasymptotic with explicit constants and is able to recover almost all the known results for Fused Lasso (both mean and quantile regression) with additional improvements in some cases.

The computation of the mean perimeter density via the notion of mean covariogram for non-stationary Boolean models has been proposed as further work in Galerne (Image Anal. Stereol. 30 (2011) 39–51). We address this issue by considering here more general germ-grain models. Furthermore, we discuss similarities and differences with respect to the computation of the mean boundary density by means of the outer Minkowski content notion.

Aldous, Evans and Pitman (1998) studied the behavior of the fragmentation process derived from deleting the edges of a uniform random tree on n labelled vertices. In particular, they showed that, after proper rescaling, the above fragmentation process converges as $n\to \mathrm{\infty }$ to the fragmentation process of the Brownian CRT obtained by cutting-down the Brownian CRT along its skeleton in a Poisson manner.

In this work, we continue the above investigation and study the fragmentation process obtained by deleting randomly chosen edges from a critical Galton–Watson tree ${\mathbf{t}}_n$ conditioned on having n vertices, whose offspring distribution belongs to the domain of attraction of a stable law of index $\mathit{\alpha }\in (1,2]$. Our main results establish that, after rescaling, the fragmentation process of ${\mathbf{t}}_n$ converges as $n\to \mathrm{\infty }$ to the fragmentation process obtained by cutting-down at a rate proportional to the length measure on the skeleton of an α-stable Lévy tree. We further show that the latter can be constructed by considering the partitions of the unit interval induced by the normalized α-stable Lévy excursion with a deterministic drift studied by Miermont (2001). This extends the result of Bertoin (2000) on the fragmentation process of the Brownian CRT.

We consider a stationary random field indexed by an increasing sequence of subsets of ${\mathbb{Z}}^d$ obeying a very broad geometrical assumption on how the sequence expands. Under certain mixing and local conditions, we show how the tail distribution of the individual variables relates to the tail behavior of the maximum of the field over the index sets in the limit as the index sets expand.

In a framework where we let the increasing index sets be scalar multiplications of a fixed set C, potentially with different scalars in different directions, we use two cluster definitions to define associated cluster counting point processes on the rescaled index set C; one cluster definition divides the index set into more and more boxes and counts a box as a cluster if it contains an extremal observation. The other cluster definition that is more intuitive considers extremal points to be in the same cluster, if they are close in distance. We show that both cluster point processes converge to a Poisson point process on C. Additionally, we find a limit of the mean cluster size. Finally, we pay special attention to the case without clusters.

The angular measure on the unit sphere characterizes the first-order dependence structure of the components of a random vector in extreme regions and is defined in terms of standardized margins. Its statistical recovery is an important step in learning problems involving observations far away from the center. In the common situation that the components of the vector have different distributions, the rank transformation offers a convenient and robust way of standardizing data in order to build an empirical version of the angular measure based on the most extreme observations. However, the study of the sampling distribution of the resulting empirical angular measure is challenging. It is the purpose of the paper to establish finite-sample bounds for the maximal deviations between the empirical and true angular measures, uniformly over classes of Borel sets of controlled combinatorial complexity. The bounds are valid with high probability and, up to logarithmic factors, scale as the square root of the effective sample size. The bounds are applied to provide performance guarantees for two statistical learning procedures tailored to extreme regions of the input space and built upon the empirical angular measure: binary classification in extreme regions through empirical risk minimization and unsupervised anomaly detection through minimum-volume sets of the sphere.

We study a semiparametric model for robust analysis of panel count data with an informative terminal event. To explore the explicit effect of the terminal event on recurrent events of interest, we propose a conditional mean model for a reversed counting process anchoring at the terminal event. Treating the distribution function of the terminal event as a nuisance functional parameter, we develop a predicted least squares-based two-stage estimation procedure with the spline-based sieve estimation technique, and derive the convergence rate of the proposed estimator. Furthermore, overcoming the difficulties caused by the convergence rate slower than $1∕\sqrt{n}$, we establish the asymptotic normality for the estimator of the finite-dimensional parameter and a functional of the estimator of the infinite-dimensional parameter. The proposed method is evaluated through extensive simulation studies and illustrated with an application to the Longitudinal Healthy Longevity Survey study on elder people in China.

The Grenander estimator (nonparametric maximum likelihood estimator) of a decreasing density function is not consistent at boundaries of support. This phenomenon has great influences on the global measures of deviation between density function and its Grenander estimator. In this article, by strong approximation technique and comparison method, we establish Cramér type moderate deviations for the Grenander estimator near the boundaries of the support. The obtained results provide a nice uniform comparison between tail probability of the estimator and limiting fluctuation distribution.

This paper is concerned with estimation and inference for the location of a change point in the mean of independent high-dimensional data. Our change point location estimator maximizes a new U-statistic based objective function, and its convergence rate and asymptotic distribution after suitable centering and normalization are obtained under mild assumptions. Our estimator turns out to have better efficiency as compared to the least squares based counterpart in the literature. Based on the asymptotic theory, we construct a confidence interval by plugging in consistent estimates of several quantities in the normalization. We also provide a bootstrap-based confidence interval and state its asymptotic validity under suitable conditions. Through simulation studies, we demonstrate favorable finite sample performance of the new change point location estimator as compared to its least squares based counterpart, and our bootstrap-based confidence intervals, as compared to several existing competitors. The asymptotic theory based on high-dimensional U-statistic is substantially different from those developed in the literature and is of independent interest. The usefulness of our bootstrap-based confidence interval is illustrated in a genomics data set.

We prove that the random simple connected cubic planar graph ${\mathsf{C}}_n$ with an even number n of vertices admits a novel uniform infinite cubic planar graph (UICPG) as quenched local limit. We describe how the limit may be constructed by a series of random blow-up operations applied to the dual map of the type III Uniform Infinite Planar Triangulation established by Angel and Schramm (Comm. Math. Phys. 241 (2003) 191–213). Our main technical lemma is a contiguity relation between ${\mathsf{C}}_n$ and a model where the networks inserted at the links of the largest 3-connected component of ${\mathsf{C}}_n$ are replaced by independent copies of a specific Boltzmann network. We prove that the number of vertices of the largest 3-connected component concentrates at $\mathrm{\kappa }n$ for $\mathrm{\kappa }\approx 0.85085$, with Airy-type fluctuations of order $n^{2∕3}$. The second-largest component is shown to have significantly smaller size $O_p(n^{2∕3})$.

In Euclidean spaces, the empirical mean vector as an estimator of the population mean is known to have polynomial concentration unless a strong tail assumption is imposed on the underlying probability measure. The idea of median-of-means tournament has been considered as a way of overcoming the sub-optimality of the empirical mean vector. In this paper, to address the sub-optimal performance of the empirical mean in a more general setting, we consider general Polish spaces with a general metric, which are allowed to be non-compact and of infinite-dimension. We discuss the estimation of the associated population Fréchet mean, and for this we extend the existing notion of median-of-means to this general setting. We devise several new notions and inequalities associated with the geometry of the underlying metric, and using them we study the concentration properties of the extended notions of median-of-means as the estimators of the population Fréchet mean. We show that the new estimators achieve exponential concentration under only a second moment condition on the underlying distribution, while the empirical Fréchet mean has polynomial concentration. We focus our study on spaces with non-positive Alexandrov curvature since they afford slower rates of convergence than spaces with positive curvature. We note that this is the first work that derives non-asymptotic concentration inequalities for extended notions of the median-of-means in non-vector spaces with a general metric.

We construct a generalization of the Ornstein–Uhlenbeck processes on the cone of covariance matrices endowed with the Log-Euclidean and the Affine-Invariant metrics. Our development exploits the Riemannian geometric structure of symmetric positive definite matrices viewed as a differential manifold. We then provide Bayesian inference for discretely observed diffusion processes of covariance matrices based on an MCMC algorithm built with the help of a novel diffusion bridge sampler accounting for the geometric structure. Our proposed algorithm is illustrated with a real data financial application.

In this article, we consider a stochastic partial differential equation (SPDE) driven by Gaussian colored noise with Dirichlet, Neumann or mixed nonhomogeneous random boundary conditions when the drift and diffusion coefficients are non-Lipschitz. We prove the existence of a unique strong solution to this SPDE and obtain a comparison theorem between such SPDEs. We also study the Hölder continuity of the solution in both time and space variables, and find the dependence of the Hölder exponent on that of the Dirichlet boundary.

We consider the i.i.d. Bernoulli field ${\mathrm{\mu }}_p$ on ${\mathbb{Z}}^d$ with occupation density $p\in [0,1]$. To each realization of the set of occupied sites we apply a thinning map that removes all occupied sites that are isolated in graph distance. We show that, while this map seems non-invasive for large p, as it changes only a small fraction $p{(1-p)}^{2d}$ of sites, there is $p(d)<1$ such that for all $p\in (p(d),1)$ the resulting measure is a non-Gibbsian measure, i.e., it does not possess a continuous version of its finite-volume conditional probabilities. On the other hand, for small p, the Gibbs property is preserved.

We consider the estimation of a non-linear reaction term in the stochastic heat or more generally in a semi-linear stochastic partial differential equation (SPDE). Consistent inference is achieved by studying a small diffusivity level, which is realistic in applications. Our main result is a central limit theorem for the estimation error of a parametric estimator, from which confidence intervals can be constructed. Statistical efficiency is demonstrated by establishing local asymptotic normality. The estimation method is extended to local observations in time and space, which allows for non-parametric estimation of a reaction intensity varying in time and space. Furthermore, discrete observations in time and space can be handled. The statistical analysis requires advanced tools from stochastic analysis like Malliavin calculus for SPDEs, the infinite-dimensional Gaussian Poincaré inequality and regularity results for SPDEs in $L^p$-interpolation spaces.

Additive models have been studied as a way of overcoming theoretical and practical difficulties in estimating a multivariate nonparametric regression function. Several methods have been proposed that ensure the optimal univariate rate one can achieve in estimating univariate nonparametric functions. In this paper a new method is proposed which reduces the constant factor in the first-order approximation of the average squared error of the most successful existing method. The new estimator is based on an orthogonal decomposition of the underlying regression function, with an arbitrarily chosen parametric family, under a special inner product structure arising from the bias formula of the estimator. It is shown that the proposed method entails reduction in the constant factor of the leading bias of the existing method while it retains the same first-order variance. These theoretical findings are confirmed in Monte Carlo experiments.

We propose a Langevin diffusion-based algorithm for non-convex optimization and sampling on a product manifold of spheres. Under a logarithmic Sobolev inequality, we establish a guarantee for finite iteration convergence to the Gibbs distribution in terms of Kullback–Leibler divergence. We show that with an appropriate temperature choice, the suboptimality gap to the global minimum is guaranteed to be arbitrarily small with high probability.

As an application, we consider the Burer–Monteiro approach for solving a semidefinite program (SDP) with diagonal constraints, and analyze the proposed Langevin algorithm for optimizing the non-convex objective. In particular, we establish a logarithmic Sobolev inequality for the Burer–Monteiro problem when there are no spurious local minima, but under the presence saddle points. Combining the results, we then provide a global optimality guarantee for the SDP and the Max-Cut problem. More precisely, we show that the Langevin algorithm achieves ϵ accuracy with high probability in $\tilde{\mathrm{\Omega }}({\mathit{ϵ}}^{-4.5})$ iterations.

Gaussian couplings of partial sum processes are derived for the high-dimensional regime $d=o(n^{1∕3})$. The coupling is derived for sums of independent random vectors and subsequently extended to nonstationary time series. Our inequalities depend explicitly on the dimension and on a measure of nonstationarity, and are thus also applicable to arrays of random vectors. To enable high-dimensional statistical inference, a feasible Gaussian approximation scheme is proposed. Applications to sequential testing and change-point detection are described.

We introduce a location statistic for distributions on non-linear geometric spaces, the diffusion mean, serving as an extension and an alternative to the Fréchet mean. The diffusion mean arises as the generalization of Gaussian maximum likelihood analysis to non-linear spaces by maximizing the likelihood of a Brownian motion. The diffusion mean depends on a time parameter t, which admits the interpretation of the allowed variance of the diffusion. The diffusion t-mean of a distribution X is the most likely origin of a Brownian motion at time t, given the end-point distribution X. We give a detailed description of the asymptotic behavior of the diffusion estimator and provide sufficient conditions for the diffusion estimator to be strongly consistent. Particularly, we present a smeary central limit theorem for diffusion means and we show that joint estimation of the mean and diffusion variance rules out smeariness in all directions simultaneously in general situations. Furthermore, we investigate properties of the diffusion mean for distributions on the sphere ${\mathcal{S}}^m$. Experimentally, we consider simulated data and data from magnetic pole reversals, all indicating similar or improved convergence rate compared to the Fréchet mean. Here, we additionally estimate t and consider its effects on smeariness and uniqueness of the diffusion mean for distributions on the sphere.

In this paper, we investigate the asymptotic behavior of spiked eigenvalues of the noncentral Fisher matrix defined by ${\mathbf{F}}_p={\mathbf{C}}_n{({\mathbf{S}}_N)}^{-1}$, where ${\mathbf{C}}_n$ is a noncentral sample covariance matrix defined by $({\mathbf{\Xi }}_n+{\mathbf{X}}_n){({\mathbf{\Xi }}_n+{\mathbf{X}}_n)}^{\ast }∕n$ and ${\mathbf{S}}_N={\mathbf{Y}}_N{\mathbf{Y}}_N^{\ast }∕N$. The matrices ${\mathbf{X}}_n$ and ${\mathbf{Y}}_N$ are two independent random matrices with standard Gaussian entries, and ${\mathbf{\Xi }}_n$ is a nonrandom matrix. When the dimensions of ${\mathbf{X}}_n$ and ${\mathbf{Y}}_N$ grow to infinity proportionally, we establish a phase transition of the spiked eigenvalues of ${\mathbf{F}}_p$. Furthermore, we derive the central limit theorem (CLT) for these spiked eigenvalues. As a byproduct of the proof of the above results, we provide the fluctuations of the spiked eigenvalues of ${\mathbf{C}}_n$, which are also of interest. Furthermore, we develop the limits and CLT for the sample canonical correlation coefficients using the results of the spiked noncentral Fisher matrix and present consistent estimators of the population spiked eigenvalues and the population canonical correlation coefficients.

It is well known that most of the existing theoretical results in statistics are based on the assumption that the sample is generated with replacement from an infinite population. However, in practice, available samples are almost always collected without replacement. If the population is a finite set of real numbers, whether we can still safely use the results from samples drawn without replacement becomes an important problem. In this paper, we focus on the eigenvalues of high-dimensional sample covariance matrices generated without replacement from finite populations. Specifically, we derive the Tracy-Widom laws for their largest eigenvalues and apply these results to parallel analysis. We provide new insight into the permutation methods proposed by Buja and Eyuboglu in [Multivar Behav Res. 27(4) (1992) 509–540]. Simulation and real data studies are conducted to demonstrate our results.

We consider the hypothesis testing problem that two vertices i and j of a generalized random dot product graph have the same latent positions, possibly up to scaling. Special cases of this hypothesis testing problem include testing whether two vertices in a stochastic block model or degree-corrected stochastic block model graph have the same block membership vectors, or testing whether two vertices in a popularity adjusted block model have the same community assignment. We propose several test statistics based on the empirical Mahalanobis distances between the ith and jth rows of either the adjacency or the normalized Laplacian spectral embedding of the graph. We show that, under mild conditions, these test statistics have limiting chi-square distributions under both the null and local alternative hypothesis, and we derive explicit expressions for the non-centrality parameters under the local alternative. Using these limiting results, we address the model selection problems including choosing between the standard stochastic block model and its degree-corrected variant, and choosing between the Erdős–Rényi model and stochastic block model. The effectiveness of our proposed tests is illustrated via both simulation studies and real data applications.

We investigate a family of discrete-time stationary processes defined by multiple stable integrals and renewal processes with infinite means. The model may exhibit behaviors of short-range or long-range dependence, respectively, depending on the parameters. The main contribution is to establish a phase transition in terms of the tail processes that characterize local clustering of extremes. Moreover, in the short-range dependence regime, the model provides an example where the extremal index is different from the candidate extremal index.

We consider the stochastic behavior of a class of local U-statistics of Poisson processes—which include subgraph and simplex counts as special cases, and amounts to quantifying clustering behavior—for point clouds lying in diverging halfspaces. We provide limit theorems for distributions with light and heavy tails. In particular, we prove finite-dimensional central limit theorems. In the light tail case we investigate tails that decay at least as slow as exponential and at least as fast as Gaussian. These results also furnish as a corollary that U-statistics for halfspaces diverging at different angles are asymptotically independent, and that there is no asymptotic independence for heavy-tailed densities. Using state-of-the-art bounds derived from recent breakthroughs combining Stein’s method and Malliavin calculus, we quantify the rate of this convergence in terms of Kolmogorov distance. We also investigate the behavior of local U-statistics of a Poisson process conditioned to lie in a diverging halfspace and find that the upper bound on the Kolmogorov distance to a standard normal distribution is smaller the lighter the tail of the density is.

In this paper we derive non-classical Tauberian asymptotics at infinity for the tail, the density and its derivatives of a large class of exponential functionals of subordinators. More precisely, we consider the case for which the Lévy measure of the subordinator satisfies the well-known and mild condition of positive increase. This is achieved via a convoluted application of the saddle point method to the Mellin transform of these exponential functionals which is given in terms of Bernstein–gamma functions. To apply the saddle point method, we improved the Stirling type asymptotics for Bernstein–gamma functions in the complex plane. As an application we derive the asymptotics of the density and its derivatives for all exponential functionals of a broad class of non-decreasing, potentially killed compound Poisson processes, which turns out to be precisely as that of an exponentially distributed random variable. We show further that a large class of densities are even analytic in a cone of the complex plane.

We derive non-asymptotic minimax bounds for the Hausdorff estimation of d-dimensional submanifolds $M\subset {\mathbb{R}}^D$ with (possibly) non-empty boundary $\partial M$. The model reunites and extends the most prevalent ${\mathcal{C}}^2$-type set estimation models: manifolds without boundary, and full-dimensional domains. We consider both the estimation of the manifold M itself and that of its boundary $\partial M$ if non-empty. Given n samples, the minimax rates are of order $O\left({(logn∕n)}^{2∕d}\right)$ if $\partial M=\varnothing$ and $O\left({(logn∕n)}^{2∕(d+1)}\right)$ if $\partial M\ne \varnothing$, up to logarithmic factors. In the process, we develop a Voronoi-based procedure that allows to identify enough points $O\left({(logn∕n)}^{2∕(d+1)}\right)$-close to $\partial M$ for reconstructing it. Explicit constant derivations are given, showing that these rates do not depend on the ambient dimension $D\gg d$.

This paper establishes necessary and sufficient conditions for the asymptotic normality of higher order Turing estimators. It further gives several easy to verify sufficient conditions. These conditions are then used to show that the assumptions hold for large classes of distributions with regularly varying tails. This includes classes for which asymptotic normality had not been previously verified.

Let L be a linear differential operator acting on functions defined over an open set $\mathcal{D}\subset {\mathbb{R}}^d$. In this article, we characterize the measurable second order random fields $U={(U(x))}_{x\in \mathcal{D}}$ whose sample paths all verify the partial differential equation (PDE) $L(u)=0$, solely in terms of their first two moments. When compared to previous similar results, the novelty lies in that the equality $L(u)=0$ is understood in the sense of distributions, which is a powerful functional analysis framework mostly designed to study linear PDEs. This framework enables to reduce to the minimum the required differentiability assumptions over the first two moments of ${(U(x))}_{x\in \mathcal{D}}$ as well as over its sample paths in order to make sense of the PDE $L(U_{\mathit{\omega }})=0$. In view of Gaussian process regression (GPR) applications, we show that when ${(U(x))}_{x\in \mathcal{D}}$ is a Gaussian process (GP), the sample paths of ${(U(x))}_{x\in \mathcal{D}}$ conditioned on pointwise observations still verify the constraint $L(u)=0$ in the distributional sense. We finish by deriving a simple but instructive example, a GP model for the 3D linear wave equation, for which our theorem is applicable and where the previous results from the literature do not apply in general.

Given n samples from a population of individuals belonging to different species, what is the number U of hitherto unseen species that would be observed if $\mathit{\lambda }n$ new samples were collected? This is the celebrated unseen-species problem, which has been the subject of recent breakthrough studies introducing non-parametric estimators of U that are minimax near-optimal and consistent all the way up to $\mathit{\lambda }\asymp logn$. These works do not rely on assumptions on the underlying unknown distribution p of the population, and therefore, while providing a theory in its greatest generality, worst-case distributions may hamper the estimation of U in concrete settings. In this paper, we strengthen the non-parametric framework for estimating U, making use of suitable assumptions on p. Inspired by the estimation of rare probabilities in extreme value theory, and motivated by the ubiquitous power-law type distributions in many natural and social phenomena, we make use of a semi-parametric assumption of regular variation of index $\mathit{\alpha }\in (0,1)$ for the tail behaviour of p. Under this assumption, we introduce an estimator of U that is simple, linear in the sampling information, computationally efficient, and scalable to massive datasets. Then, uniformly over our class of regularly varying tail distributions, we show that the proposed estimator has provable guarantees: i) it is minimax near-optimal, up to a power of $logn$ factor; ii) it is consistent all of the way up to $log\mathit{\lambda }\asymp n^{\mathit{\alpha }∕2}∕\sqrt{logn}$, and this range is the best possible. This is the first study on the estimation of the unseen under regularly varying tail distributions p. Our results rely on a novel approach, of independent interest, which combines the renowned method of the two fuzzy hypotheses for minimax estimation of discrete functionals, with Bayesian arguments under Poisson-Kingman priors for p. An illustration of our method is presented for synthetic and real data.

This article uses a combination of three ideas from simulation to establish a nearly optimal polynomial upper bound for the joint density of the stable process and its associated supremum at a fixed time on the entire support of the joint law. The representation of the concave majorant of the stable process and the Chambers-Mallows-Stuck representation for stable laws are used to define an approximation of the random vector of interest. An interpolation technique using multilevel Monte Carlo is applied to accelerate the approximation, allowing us to establish the infinite differentiability of the joint density as well as nearly optimal polynomial upper bounds for the joint mixed derivatives of any order.