<figure><img src="image_1.png" /></figure>

HAL
open science

# Estimation for the convolution of several multidimensional densities

Fabienne Comte, Bianca Neubert

<figure><img src="image_2.png" /></figure>

## To cite this version:

Fabienne Comte, Bianca Neubert. Estimation for the convolution of several multidimensional densities. 2025. hal-05109668

HAL Id: hal-05109668

https://hal.science/hal-05109668v1

Preprint submitted on 12 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Estimation for the convolution of several multidimensional densities

FABIENNE COMTE*
Université Paris Cité

BIANCA NEUBERT **
Ruprecht-Karls-Universität Heidelberg

## Abstract

This work is concerned with the problem of estimating the *m*-fold convolution of the densities of *m* independent random vectors. Two nonparametric estimators are proposed, a kernel and a projection estimator, and their integrated quadratic risk is studied. We use Fourier analysis to bound the variance and consider standard Sobolev classes to discuss the convergence rates for the kernel estimator. In addition, we propose a bandwidth selection method for the kernel estimator and a study model selection for the projection estimator. Finally, we illustrate the results in simulation experiments.

**Keywords:** Nonparametric density estimation, Convolutions of densities, Multivariate observations, Fourier transform, Data-driven estimation

AMS subject classifications: Primary 62G07, Secondary 62H12.

## 1 Introduction

In this work, we assume that we observe *p* independent samples with possibly different sample sizes, of *d*-dimensional vectors i.e. for $j \in [p]$

$$X_k^{(j)} = (X_{k,1}^{(j)}, \dots, X_{k,d}^{(j)}) \in \mathbb{R}^d,$$

where for $k \in [n_j]$, the $X_k^{(j)}$ are independent and identically distributed (i.i.d.) with common density $f_j$. We consider the problem of estimating the *p*-fold convolution of the densities $f_j$; in other words, we study the estimation of $g^{(p)}: \mathbb{R}^d \to \mathbb{R}$ defined by, for $x \in \mathbb{R}^d$,

$$ \begin{aligned} g^{(p)}(\boldsymbol{x}) &:= f_1 * \dots * f_p(\boldsymbol{x}) \\ &= \int_{\mathbb{R}^d} \dots \int_{\mathbb{R}^d} f_1(\boldsymbol{x} - \boldsymbol{x}_2 - \dots - \boldsymbol{x}_p) f_2(\boldsymbol{x}_2) \dots f_p(\boldsymbol{x}_p) d\boldsymbol{x}_2 \dots d\boldsymbol{x}_p. \end{aligned} \quad (1.1) $$

Note that the convolution describes the density of the sum of components, i.e. $\sum_{j \in [p]} X_1^{(j)} \sim g^{(p)}$. The sum of independent random variables play a role in many fields. For example, in actuarial sciences the distribution of aggregate claims or ruin probabilities are of interest, see Panjer and Willmot (1992). For an extensive overview of further fields of application see Midhu et al. (2023). Furthermore, in the estimation of compound Poisson processes, it is possible to relate its density to the convolution of the

<footnote>*Université Paris Cité, MAP5 UMR 8145, F-75006 Paris, France, e-mail: fabienne.comte@u-paris.fr</footnote>

<footnote>**Institut für Mathematik, MATHEMATIKON, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany, e-mail: neubert@math.uni-heidelberg.de</footnote>
density of increments. In this context, estimators for convolutions of densities are of interest; see Duval (2013). Comte et al. (2014) precisely plug-in estimates of convolutions of any order of a density in order to build a specific Poisson-density estimator.

**Previous works** For some classes of distributions, it is possible to analytically derive the distribution of their convolution, e.g. Gamma distributions (Hu et al. (2019)), in particular exponentially distributed random variables, which are of interest in insurance risk theory (Willmot and Woo (2007)). In the case when the convolution of densities can not be derived analytically, a nonparametric estimator is useful, as a method of representation of an untractable analytical computation.

First, nonparametric estimators of kernel type were proposed; see notably Frees (1994) and Saavedra and Cao (2000). More precisely, Frees (1994) considers the problem of estimating $q(X_1, \dots, X_m)$ for some symmetric measurable function $q$ given i.i.d. observations $(X_i)_{i \in [n]}$ of some random variable $X$. The work was motivated by the measure of the degree of spatial randomness in a point pattern with $q(x_1, x_2) = |x_1 - x_2|$ but the results hold for a broader class of functions $q$, including the sum of observations. Consequently, this corresponds to our setting with $d = 1$ and $f_j = f$ for all $j \in [m]$. The proposed estimator is based on a U-statistic and is proven to be point-wise $\sqrt{n}$-consistent. This estimator has been extensively studied, see e.g. Ahmad and Mugdadi (2003), Bandyopadhyay (2012). Giné and Mason (2007a,b) prove functional central limit theorems in various function spaces and provide a law of iterated logarithm for the estimator. Ahmad and Fan (2001) deal with practical aspects of the choice of optimal bandwidths. Prakasa Rao (2004) derives moment inequalities for the deviation of the estimator from the true density. Further works deriving convergence in law for this estimator include a series of papers by Schick and Wefelmeyer (2004b, 2008, 2009a,b, 2012).

Motivated by estimating the density of a moving average process of order one, Saavedra and Cao (2000) consider the case of $q(X_1, X_2) = X_1 + aX_2$. They propose a plug-in kernel estimator for the density, which can be written as the convolution of two densities; they obtain $\sqrt{n}$-consistency and an expression for mean squared and integrated mean squared error. Further, they propose an extension to higher order moving average processes which corresponds to estimating the $m$-fold convolution of a density. The estimator proposed in Saavedra and Cao (2000) has also been extensively studied, e.g. Schick and Wefelmeyer (2004a) prove asymptotic normality and efficiency. In Schick and Wefelmeyer (2007) it is shown that the estimator proposed by Saavedra and Cao (2000) is asymptotically equivalent to the estimator proposed by Frees (1994), consequently, a series of results hold for both estimators. We mention that Midhu et al. (2023) propose an approximation for the distribution function of the sum of two independent nonnegative random variables using quantile-based representation. In the case of dependent variables, there is no formula equivalent to Equation (1.1) for the density of their sum, and very different tools are required; e.g. Cherubini et al. (2011) use the concept of Copula-based convolution. Here, we focus on the case of independence between variables.

Chesneau et al. (2013) consider the problem of estimating the $m$-fold convolution of a density on $\mathbb{R}$. Relying on Fourier methods, they propose an adaptive estimator based on kernel methods and prove that bandwidth selection by Lepski method offers a general squared bias-variance compromise. Nickl (2007) obtains as a by-product of general results on nonparametric MLE estimators the $\sqrt{n}$-convergence in law of the convolution of two such estimators; see further Nickl (2009). Although all former works focus on a density on $\mathbb{R}$, the latter results can also be applied to the estimation of the convolution of two densities on $\mathbb{R}^d$.

An alternative to kernel estimation is to rely on projection methods. This was considered by Guo et al. (2018) for the $m$-fold convolution of a density; an estimate is defined as the $m$-fold convolution of a wavelet density estimator, but only $L^p$-consistency is shown.
**Contribution** In this work, we study two approaches for the estimation of the density of the sum of $p$ different $d$-dimensional independent random variables.

First, we show the link between the convolution of kernel estimators in the spirit of the plug-in approach of Saavedra and Cao (2000) and the Fourier formulation used in Chesneau et al. (2013). Then, we extend this estimator to the case of possibly different densities in $d$-dimensions. We show upper bounds of the mean integrated squared error and discuss the resulting convergence rates.

Second, we propose a general projection estimator and prove an upper bound on the integrated squared error. Based on the Hermite basis as described in Belomestny et al. (2019), we also build an adaptive estimator that enjoys a parsimonious development in the basis. This is the first time general results for projection estimators in terms of risk, rates, and model selection are proved, as far as we know.

The paper is structured as follows. In Section 2 we consider the kernel estimator, starting with the case of the convolution of two independent univariate variables in Section 2.1 and extending the estimator and results to the convolution of the densities of $p$ independent multivariate variables in Section 2.2. In addition, we discuss the convergence rates in Section 2.3 and propose an adaptive procedure in Section 2.4. Then, we turn to the projection estimator in Section 3. We first consider the case of two independent univariate variables in Section 3.1 and then give the generalization in Section 3.3. In Section 3.4, we propose an adaptive procedure and show an upper bound for its integrated squared error. Finally, Section 4 illustrates the results through a simulation study. Proofs can be found in Appendix A and Appendix B.

**Notation** We use multi-index notation relying on bold letters for variables in $\mathbb{R}^d$, that is, $\boldsymbol{x} = (x_1, \dots, x_d)$. For $\boldsymbol{x} = (x_1, \dots, x_d)$, $\boldsymbol{y} = (y_1, \dots, y_d) \in \mathbb{R}^d$ when using standard computation op- erators for $\mathbb{R}$, they are to be interpreted component-wise. For example, $\boldsymbol{x}\boldsymbol{y} = (x_1y_1, \dots, x_dy_d)$; $\boldsymbol{x} \le \boldsymbol{y} \Leftrightarrow (x_k \le y_k \forall k \in \llbracket d \rrbracket)$ and taking the sum over $\boldsymbol{j} \in \llbracket \boldsymbol{m} \rrbracket$ for $\boldsymbol{m} \in \mathbb{N}^d$ is to be understood as the sum over $j_k = 1$ to $m_k$ for all $k \in \llbracket d \rrbracket$.

Let $u, v: \mathbb{R}^d \to \mathbb{R}$ be two functions, we define the convolution $u * v$, as soon as it is defined (e.g. for $u$ and $v$ square-integrable on $\mathbb{R}^d$) by $u * v(\boldsymbol{x}) = \int_{\mathbb{R}^d} u(\boldsymbol{x} - \boldsymbol{y})v(\boldsymbol{y})d\boldsymbol{y}$. The Fourier transform of an integrable function $u$ is, for $t \in \mathbb{R}^d$, given by

$$u^*(\boldsymbol{t}) := \int_{\mathbb{R}^d} e^{i\langle \boldsymbol{t}, \boldsymbol{x} \rangle} u(\boldsymbol{x})d\boldsymbol{x}.$$

The $L^2(\mathbb{R}^d)$-norm of a square-integrable function $u$ is denoted by $\|u\|^2 = \int_{\mathbb{R}^d} |u(\boldsymbol{x})|^2d\boldsymbol{x}$, the $L^\infty(\mathbb{R}^d)$-norm of a bounded function $u$ by $\|u\|_\infty = \sup_{\boldsymbol{x} \in \mathbb{R}^d} |u(\boldsymbol{x})|$.

# 2 Kernel estimator

In this section, we propose a kernel estimator for the estimation of the convolution of two one-dimensional densities in a first step, and several multidimensional densities in a second step. For the proofs of this section see Appendix A.1.

## 2.1 Convolution of the densities of two independent univariate variables

To explain the estimation problem of interest and the proposed estimation method, we first consider the simple case of two independent one-dimensional variables and two samples of the same size $n$.

Given i.i.d. observations $(X_k^{(1)})_{k \in \llbracket n \rrbracket}$ of $X^{(1)} \sim f_1$ and $(X_k^{(2)})_{k \in \llbracket n \rrbracket}$ of $X^{(2)} \sim f_2$, let us look for an estimate of their convolution $g: \mathbb{R} \to \mathbb{R}$, $g(x) := f_1 * f_2(x) := \int f_1(y)f_2(x - y)dy$. Given a kernel
function $K: \mathbb{R} \to \mathbb{R}$ (i.e. an integrable function with $\int K = 1$), we set for $j \in \{1, 2\}$

$$ \forall x \in \mathbb{R}, \quad \hat{f}_{j,h}(x) := \frac{1}{nh} \sum_{k \in [n]} K \left( \frac{X_k^{(j)} - x}{h} \right). $$

A natural proposal to estimate $f_1 * f_2$ is $\hat{f}_{1,h} * \hat{f}_{2,h}$, which we can rewrite:

$$ \begin{aligned} \hat{f}_{1,h} * \hat{f}_{2,h}(x) &= \frac{1}{n^2 h^2} \sum_{k_1 \in [n], k_2 \in [n]} \int K \left( \frac{X_{k_1}^{(1)} - (x-y)}{h} \right) K \left( \frac{X_{k_2}^{(2)} - y}{h} \right) dy \\ &= \frac{1}{n^2 h} \sum_{k_1 \in [n], k_2 \in [n]} \int K \left( \frac{X_{k_1}^{(1)} + X_{k_2}^{(2)} - x}{h} - u \right) K(u) du \\ &= \frac{1}{n^2 h} \sum_{k_1 \in [n], k_2 \in [n]} L \left( \frac{X_{k_1}^{(1)} + X_{k_2}^{(2)} - x}{h} \right), \end{aligned} $$

where $L(x) := K * K(x)$ for $x \in \mathbb{R}$. For $K$ integrable and $\int K = 1$, we also have $L$ integrable and $\int L = 1$, therefore $L$ is also a kernel. Thus, given a kernel $L$, we obtain, for $h \in \mathbb{R}_+^*$ and $x \in \mathbb{R}$,

$$ \hat{g}_h(x) := \frac{1}{n^2} \sum_{k_1 \in [n]} \sum_{k_2 \in [n]} L_h(X_{k_1}^{(1)} + X_{k_2}^{(2)} - x), \quad L_h(x) = \frac{1}{h} L\left(\frac{x}{h}\right). \quad (2.1) $$

The estimator in Equation (2.1) can be seen as a kernel estimator with kernel $L = K * K$ from observations $(X_{k_1}^{(1)} + X_{k_2}^{(2)})_{k_1 \in [n], k_2 \in [n]}$ with the important feature that the observations are not independent. Relying only on the representation eq. (2.1), we can prove the following result.

**Proposition 2.1 (Bound on the integrated quadratic risk).** If $f_1, f_2 \in L^2(\mathbb{R})$ and $L \in L^1(\mathbb{R}) \cap L^2(\mathbb{R})$ is an even kernel, then the estimator $\hat{g}_h$ of $g$ defined by (2.1) satisfies

$$ \mathbb{E} \left[ \|\hat{g}_h - g\|^2 \right] \le \|g_h - g\|^2 + \frac{\|L\|^2}{n^2 h} + \|L\|_1^2 \frac{\|f_1\|^2 + \|f_2\|^2}{n}, $$

where $g_h := L_h * g$.

**Proposition 2.1** provides a squared-bias variance decomposition, plus a residual term of order $1/n$; note that the order of the variance corresponds to what would be obtained with $n^2$ independent observations in density estimation. If we had considered the $n$ sample $(X_k^{(1)} + X_k^{(2)})_{k \in [n]}$, we would have preserved the independence between observations and would have directly applied the standard kernel estimation result, see Comte (2017), chap.3; however, we get a variance bound of order $1/(nh)$ instead of $1/(n^2h)$, as above.

We also provide a point-wise risk bound:

**Proposition 2.2 (Bound on the point-wise quadratic risk).** If $f_1, f_2 \in L^\infty(\mathbb{R}) \cap L^2(\mathbb{R})$ and $L \in L^1(\mathbb{R}) \cap L^2(\mathbb{R})$ is an even kernel, then the estimator $\hat{g}_h$ of $g$ defined by (2.1) satisfies for $x \in \mathbb{R}$ that

$$ \mathbb{E} \left[ (\hat{g}_h(x) - g(x))^2 \right] \le (g_h(x) - g(x))^2 + \|L\|_2^2 \frac{\|g\|_\infty}{n^2 h} + \frac{C_f \|L\|_1^2}{n}, $$

where $C_f := \|f_1\|_\infty \|f_2\|^2 + \|f_2\|_\infty \|f_1\|^2$.

**Remark 2.3.** Proposition 2.2 improves the bound of Proposition 1 for the case $m = 2$ in Chesneau et al. (2013), who obtain a variance term of order $1/(n^2h^2)$ instead of $1/(n^2h)$ here.
## 2.2 Convolution of the densities of p independent multivariate variables

In this section we generalize the results from Section 2.1 to the convolution of multiple multivariate densities. The goal is now the estimation of $g^{(p)}: \mathbb{R}^d \to \mathbb{R}$ defined by (1.1) from independent vector observations which are all d-dimensional but with possibly different sample sizes, i.e. for $j \in \llbracket p \rrbracket$, we observe $X_k^{(j)} = (X_{k,1}^{(j)}, \dots, X_{k,d}^{(j)}) \in \mathbb{R}^d$, with $X_k^{(j)} \sim f_j$ i.i.d., for $k \in \llbracket n_j \rrbracket$ and independent for $j \in \llbracket p \rrbracket$.

Given a kernel $L: \mathbb{R}^d \to \mathbb{R}$, we define for given $h = (h_1, \dots, h_d) \in \mathbb{R}_+^d$ and $n = (n_1, \dots, n_p) \in \mathbb{N}^p$ the estimator of $g^{(p)}$ at $x \in \mathbb{R}^d$ as

$$ \widehat{g}_h^{(p)}(\boldsymbol{x}) = \frac{1}{h_1 \cdots h_d n_1 \cdots n_p} \sum_{k_1 \in \llbracket n_1 \rrbracket} \cdots \sum_{k_p \in \llbracket n_p \rrbracket} L \left( \frac{\boldsymbol{X}_{k_1}^{(1)} + \cdots + \boldsymbol{X}_{k_p}^{(p)} - \boldsymbol{x}}{\boldsymbol{h}} \right). \quad (2.2) $$

This extends the estimator from Equation (2.1) on the one hand to the convolution of p different densities and on the other hand to the d-dimensional case using independent samples of different sizes. Our direct computations to prove Proposition 2.1 would be tedious to extend, more specifically, the variance bound. Thus, we use Fourier methods and, to that aim, we rewrite the estimator in the Fourier domain. Recall that $f_j^*$ denotes the Fourier transform of density $f_j$, and define its empirical counterpart by

$$ \forall \boldsymbol{t} \in \mathbb{R}^d, \quad \widehat{f}_j^*(\boldsymbol{t}) := \frac{1}{n_j} \sum_{k \in \llbracket n_j \rrbracket} e^{i \langle \boldsymbol{t}, \boldsymbol{X}_k^{(j)} \rangle}. $$

Then, $\widehat{g}_h^{(p)}$ defined by Equation (2.2) can also be written as

$$ \widehat{g}_h^{(p)}(\boldsymbol{x}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{-i\langle \boldsymbol{t}, \boldsymbol{x} \rangle} L^*(\boldsymbol{h}\boldsymbol{t}) \prod_{j \in \llbracket p \rrbracket} \widehat{f}_j^*(\boldsymbol{t}) d\boldsymbol{t}. \quad (2.3) $$

By the Fourier inversion formula, Equation (2.2) and Equation (2.3) are equivalent if $L$ and $L^*$ are integrable, see Appendix A.1.

If we start from formula (2.3) and take $L$ as the sinus cardinal kernel in dimension $d = 1$, that is, for $x \in \mathbb{R}$, $L(x) := \sin(\pi x)/(\pi x)$, we find that $L$ (although not integrable) admits a Fourier transform given by $L^*(hu) = \mathbb{1}_{[-\pi, \pi]}(hu) = \mathbb{1}_{[-\pi/h, \pi/h]}(u)$. With this kernel choice, the estimator writes $\widehat{g}_h^{(p)}(x) = \frac{1}{2\pi} \int_{[-\pi/h, \pi/h]} e^{-iux} \prod_{j \in \llbracket p \rrbracket} \widehat{f}_j^*(u) du$. This shows that Equation (2.3) extends the estimator (7) of Chesneau et al. (2013) in three directions: first, to the case of differing densities; second, for a general kernel; third, to multidimensional densities. Note that we can use both kernel methods to analyze the estimator using Equation (2.2) or Fourier analysis using Equation (2.3), the assumptions required on $L$ may differ, but with no main drawback, as the kernel is user-chosen.

We consider estimator (2.3) with a general kernel $L$ and determine the integrated quadratic risk. We consider the following assumptions.

**Assumption A1.** (i) The kernel $L$ is an even function which belongs to $\mathbb{L}^2(\mathbb{R}^d)$.

(ii) Its Fourier transform $L^*$ belongs to $\mathbb{L}^2(\mathbb{R}^d)$ and $L^*$ is bounded.

Under Assumption A1, we can define

$$ g_h^{(p)}(\boldsymbol{x}) := L_h * g^{(p)}(\boldsymbol{x}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{-i\langle \boldsymbol{t}, \boldsymbol{x} \rangle} L^*(\boldsymbol{h}\boldsymbol{t}) \prod_{j \in \llbracket p \rrbracket} f_j^*(\boldsymbol{t}) d\boldsymbol{t}. $$
where

$$L_h(\boldsymbol{x}) = \frac{1}{h_1 \dots h_d} L\left(\frac{x_1}{h_1}, \dots, \frac{x_d}{h_d}\right).$$

**Proposition 2.4 (Bound on the integrated quadratic risk - Kernel estimator).** Assume that $L$ satisfies Assumption A1, and that $f_j^* \in L^1(\mathbb{R}^d) \cap L^2(\mathbb{R}^d)$ for all $j \in [p]$. Then the estimator $\tilde{g}_h^{(p)}(x)$ of $g^{(p)}$ defined in Equation (2.3) satisfies

$$\mathbb{E} \left[ \| \widehat{g_h^{(p)}} - g^{(p)} \|^2 \right] \le \| g_h^{(p)} - g^{(p)} \|^2 + \frac{2 \| L^* \|^2}{h_1 \cdots h_d n_1 \cdots n_p} + C \sum_{j \in [p]} \frac{1}{n_j},$$

with $C := 2p3^{2(p-2)} \|L^*\|_{\infty}^2 \max_{j \in [p]} \|f_j^*\|^2$.

The upper bound on the integrated quadratic risk consists of the squared bias $\|g_h^{(p)} - g^{(p)}\|^2$, an extended variance of order $1/(h_1 \cdots h_d n_1 \cdots n_p)$ and a residual term of order $\max_{j \in [p]} 1/n_j$ which does not depend on the bandwidth.

Note that if one sample size is small the residual term causes the order of the upper bound to be bad. If all sample sizes are of the same order $n_j = c_j n$, then the variance term is simplified to be of the order $1/(h_1 \cdots h_d n^p)$ and the residual term of the order $1/n$.

Note also that, if $L(\boldsymbol{x}) = \prod_{k \in [d]} L_k(x_k)$, then $L^*(\boldsymbol{t}) = \prod_{k \in [d]} L_k^*(t_k)$, and all $L_k$ can be taken as sinus cardinal. In that case $\|g_h^{(p)} - g^{(p)}\|^2 = \int_{\mathbb{R}^d \setminus \prod_{k \in [d]} [-\pi/h_k, \pi h_k]} |(g^{(p)})^*(\boldsymbol{t})|^2 dt$.

The point-wise risk is omitted for the sake of brevity but may also be handled similarly.

## 2.3 Convergence rates

Let us consider convergence rates for the kernel estimator defined in Equation (2.3). More precisely, we discuss the order of the upper bound given in Proposition 2.4. We focus on the case of anisotropic Sobolev classes with an additional boundedness assumption generalizing the one proposed in Chesneau et al. (2013) as follows:

**Definition (D1) (Anisotropic Sobolev space).** A function $f: \mathbb{R}^d \to \mathbb{R}_+$ belongs to Sobolev space $S(\alpha, R)$ for $\alpha \in \mathbb{R}_+^d$ and $R \in \mathbb{R}_+$ if $f \in L^2(\mathbb{R}^d)$ and for its Fourier transform $f^*$ it holds that

$$\sum_{k \in [d]} \int_{\mathbb{R}^d} (1 + t_k^2)^{\alpha_k} |f^*(\boldsymbol{t})|^2 dt \le R^2 \quad \text{and} \quad \sup_{k \in [d]} \sup_{\boldsymbol{t} \in \mathbb{R}^d} (1 + t_k^2)^{\alpha_k} |f^*(\boldsymbol{t})|^2 \le R^2.$$

To be able to control the bias, we need to ensure that the Fourier transform of kernel $L$ is a real valued function. Further, we consider the following assumptions as introduced by Taupin (2001) for dimension $d = 1$.

**Assumption A2.** The Fourier transform of $L$ satisfies $L^*(\boldsymbol{t}) = 1$ for any $\boldsymbol{t} \in [-1, 1]^d$.

For dimension $d = 1$, the sinus cardinal kernel satisfies Assumption A1 and Assumption A2. Taupin (2001) points out that the de La Vallée-Poussin kernel $V$ defined by

$$V(x) = \frac{\cos(x) - \cos(2x)}{\pi x^2}$$

also satisfies Assumption A1 and Assumption A2.
Now we can discuss the convergence rates for the quadratic integrated risk. We denote by $\bar{\alpha}$ the harmonic mean of the regularity components $\alpha_j$, i.e.

$$ \frac{1}{\bar{\alpha}} = \frac{1}{d} \left( \frac{1}{\alpha_1} + \dots + \frac{1}{\alpha_d} \right). $$

**Proposition 2.5.** Let Assumption A1 and Assumption A2 be satisfied. Assume that all $n_j$ are of the same order $n \in \mathbb{N}$, that is, $n_j = c_j n$ for $c_j$ a constant independent of $n$, $j \in [\lfloor p \rfloor]$ and that $g^{(p)} \in S(\alpha, R)$, for $\alpha \in \mathbb{R}_+$ and $R \in \mathbb{R}_+$. Then, with $h_{k,\text{opt}} = n^{-\frac{p\bar{\alpha}}{\alpha_k(2\bar{\alpha}+d)}}$ for $k \in [\lfloor d \rfloor]$ and some $C > 0$ the resulting rate is

$$ \mathrm{E} \left[ \| \widehat{g}_{h_{\text{opt}}}^{(p)} - g^{(p)} \|^2 \right] \le C \begin{cases} n^{-1}, & \text{if } \bar{\alpha} \ge \frac{d}{2(p-1)}, \\ n^{-\frac{2p\bar{\alpha}}{2p\bar{\alpha}+d}}, & \text{otherwise.} \end{cases} $$

*Proof of Proposition 2.5.* First, due to Assumption A1 and $g^{(p)} \in S(\alpha, R)$, we have that

$$ \|g_h^{(p)} - g^{(p)}\|^2 \le C(h_1^{2\alpha_1} + \dots + h_d^{2\alpha_d}), \qquad (2.4) $$

where $C = \frac{R^2(1+\|L^*\|_{\infty})^2}{(2\pi)^d}$. With Proposition 2.4, we obtain for some constant $C > 0$

$$ \mathrm{E} \left[ \| \widehat{g}_h^{(p)} - g^{(p)} \|^2 \right] \le C \left( h_1^{2\alpha_1} + \dots + h_d^{2\alpha_d} + \frac{1}{h_1 \dots h_d n^p} + \frac{1}{n} \right) $$

and

$$ \mathrm{E} \left[ \| \widehat{g}_{h_{\text{opt}}}^{(p)} - g^{(p)} \|^2 \right] \le C \max \left( n^{-\frac{2\bar{\alpha}p}{2\bar{\alpha}p+d}}, n^{-1} \right) \\ = C \begin{cases} n^{-1}, & \text{if } \bar{\alpha} \ge \frac{d}{2(p-1)}, \\ n^{-\frac{2\bar{\alpha}p}{2\bar{\alpha}p+d}}, & \text{otherwise.} \end{cases} $$

This concludes the proof. $\square$

Note that for increasing components $p$, keeping all other parameters constant, the rate is $n^{-1}$. However, since the factor $p^2 3^2 (p-2)$ is included in the constant $C$, $p$ has to grow sufficiently slow in comparison to $n$. On the other hand, if $d$ is increasing, assuming that all components have the same regularity, then we have rate $n^{-\frac{2\bar{\alpha}p}{2\bar{\alpha}p+d}}$, which gets slower. Again $d$ is included in the constant $C$ such that it has to grow sufficiently slower than the sample size $n$.

**Remark 2.6.** We can deduce the regularity of $g^{(p)}$ from the regularities of the densities $f_j$. Indeed, if $f_j \in S((\alpha_{j,1}, \dots, \alpha_{j,d}), R_j)$ for $j \in [\lfloor p \rfloor]$, that is each component density belongs to a Sobolev space, then the convoluted density $g^{(p)}$ also belongs to a Sobolev space. More precisely, $g^{(p)} \in S((\sum_{j \in [\lfloor p \rfloor]} \alpha_{j,1}, \dots, \sum_{j \in [\lfloor p \rfloor]} \alpha_{j,d}), \prod_{j \in [\lfloor p \rfloor]} R_j)$. Consequently, the regularity in the rates is the harmonic mean of the directional regularities of the function $g^{(p)}$, i.e.

$$ \frac{1}{\bar{\alpha}} = \frac{1}{d} \sum_{k \in [\lfloor d \rfloor]} \frac{1}{\sum_{j \in [\lfloor p \rfloor]} \alpha_{j,k}}. $$

Furthermore, if all directions and all components are of the same regularity $a \in \mathbb{R}_+$, that is $\alpha_{j,k} = a$ for all $j \in [\lfloor p \rfloor]$, $k \in [\lfloor d \rfloor]$, then the regularity appearing in the convergence rate simplifies to $\bar{\alpha} = pa$.

Alternatively, the squared bias also has the form (2.4), if the Kernel $L$ is of order $\sum_{k \in [\lfloor d \rfloor]} \lfloor \alpha_k \rfloor$ and $g^{(p)}$ belongs to an anisotropic Nikolski space, see Comte and Lacour (2013) for a definition. Therefore in such cases, the rates would be as given by Proposition 2.5. It may also be possible to get a link between the regularity of single densities $f_j$ and their convolution $g^{(p)}$ similar to Remark 2.6.
## 2.4 Bandwidth selection with the PCO-method

Since the optimal choice of the bandwidth parameter $h$ for the kernel estimator (2.2) depends on the regularity of the density $g^{(p)}$ and is asymptotic, which is unknown in practice, a natural question is to find a finite-sample data-driven approach to select these parameters.

We propose to use the approach presented in Lacour et al. (2017), called Penalized Comparison to Overfitting (PCO) method. However, theoretical results are lengthy and since our small simulation study (see Section 4 for more details) suggests that the performance of the model selection method for projection estimation proposed in Section 3.4 is slightly better, we show upper bounds for the integrated squared risk in the latter case only. Let us explain the PCO-method in case $d = 1$, starting from the bias-variance decomposition $\mathbb{E}[\|g_h^{(p)} - g^{(p)}\|^2] = \mathbb{E}[\|g_h^{(p)} - g_h^{(p)}\|^2] + \|g_h^{(p)} - g^{(p)}\|^2 =: v_h + b_h$. It is natural to propose a criterion of the form $\text{Crit}(h) = \hat{b}_h + \hat{v}_h$. From the variance bound proved in Proposition 2.4, we set

$$ \hat{v}_h = \frac{\lambda \|L\|^2}{hn_1 \cdots n_p} = \frac{\lambda \|L_h\|^2}{n_1 \cdots n_p} $$

with tuning parameter $\lambda$. For the estimation of the bias, assume that $h_{\min}$, the minimum of the bandwidths in the grid $\mathcal{H}$, is very small. Plugging $\hat{g}_{h_{\min}}^{(p)}$ and $\hat{g}_h^{(p)}$ in the bias term introduces a bias, consequently, a corrected term is proposed, and we get

$$ \hat{b}_h \approx \|\hat{g}_{h_{\min}}^{(p)} - \hat{g}_h^{(p)}\|^2 - \frac{\|L_{h_{\min}} - L_h\|^2}{n_1 \cdots n_p}. $$

Choosing $\lambda = 1$, as suggested in the remark following the oracle inequality in Lacour et al. (2017), this yields the criterion

$$ \text{Crit}(h) = \|\hat{g}_{h_{\min}}^{(p)} - \hat{g}_h^{(p)}\|^2 + \frac{2\langle L_h, L_{h_{\min}} \rangle}{n_1 \cdots n_p}, \qquad (2.5) $$

and the estimator $\hat{h} = \arg \min_{h \in \mathcal{H}} \text{Crit}(h)$ for the bandwidth.

# 3 Projection estimator

We propose now a projection strategy for the estimation of the convolution, first of two one-dimensional densities and next, of several multidimensional densities. Proofs of this section can be found in Appendix A.2.

## 3.1 Convolution of the densities of two independent univariate variables

Let us explain the estimation method in the simple case of two independent one-dimensional random variables, given i.i.d. observations $(X_k^{(1)})_{k \in \llbracket n \rrbracket}$ of $X^{(1)} \sim f_1$ and $(X_k^{(2)})_{k \in \llbracket n \rrbracket}$ of $X^{(2)} \sim f_2$. Let $A \subseteq \mathbb{R}$ and assume that $g \in \mathbb{L}^2(A)$. Given an orthonormal basis $(\varphi_\ell)_{\ell \ge 0}$ of $\mathbb{L}^2(A)$ we can write $g = \sum_{\ell \ge 0} a_\ell \varphi_\ell$ with $a_\ell = \langle g, \varphi_\ell \rangle$. We define

$$ B(m) := \sup_{x \in A} \sum_{\ell=0}^{m-1} \varphi_{\ell}^{2}(x). $$
Let $g_m$ denote the orthogonal projection of $g$ on $S_m := \text{vect}\{\varphi_0, \dots, \varphi_{m-1}\}$, i.e. $g_m = \sum_{\ell=0}^{m-1} a_\ell \varphi_\ell$. To build an estimator of $g$, we estimate the first coefficients of its development in the basis and get, for $m \ge 1$ the proposal

$$ \hat{g}_m := \sum_{\ell=0}^{m-1} \hat{a}_\ell \varphi_\ell, \quad \hat{a}_\ell := \frac{1}{n^2} \sum_{k_1 \in [n]} \sum_{k_2 \in [n]} \varphi_\ell (X_{k_1}^{(1)} + X_{k_2}^{(2)}). \qquad (3.1) $$

As $\mathbb{E}[\hat{a}_\ell] = \langle \varphi_\ell, g \rangle = a_\ell$, this provides an unbiased estimator of $g_m$. Note that this estimator is different from the convolution of two projection density estimators. We prove the following bound.

**Proposition 3.1.** If $f_1, f_2 \in \mathbb{L}^\infty(\mathbb{R}) \cap \mathbb{L}^2(\mathbb{R})$, it holds for the estimator defined in Equation (3.1) that

$$ \mathbb{E} [\|\hat{g}_m - g\|^2] \le \|g_m - g\|^2 + \frac{B(m) \wedge m \|g\|_\infty}{n^2} + \frac{\|f_1\|^2 + \|f_2\|^2}{n}. $$

Similarly to the kernel estimator, we cannot apply standard methods for the variance since the $n^2$ variables involved in the definition of $\hat{a}_\ell$ are not independent but direct computation is still possible. Proposition 3.1 holds for any orthonormal basis. For example, the trigonometric basis satisfies $B(m) = m$ (Comte (2017), Remark 1.4, Chap. 2), for the Legendre polynomial basis $B(m) \sim m^2$ (Cohen et al. (2013), p. 831) and for the Hermite basis, it holds $B(m) \sim \sqrt{m}$ (Comte and Lacour (2023), Lemma 1 and Section 3.2 hereafter). Consequently, for the Hermite basis $B(m)$ dominates in the second summand and the order of the variance term improves to $\sqrt{m}/n^2$ which is coherent with the influence of the dimension of the Hermite space on the order of the bias term (see Belomestny et al. (2019)).

## 3.2 Hermite basis

To extend the definition of the projection estimator to the convolution of more than two variables, it is convenient to rely on a basis admitting a regular Fourier transform. This is why we develop extensions relying on the Hermite basis, and thus start by giving its definition and useful properties.

We start from the orthonormal Hermite basis for $A = \mathbb{R}$, and then build a basis for $\mathbb{R}^d$ by tensorization. The Hermite basis $(\varphi_\ell)_{\ell \ge 0}$ is a basis on $\mathbb{L}^2(\mathbb{R})$ defined from the Hermite polynomials $(H_\ell)_{\ell \ge 0}$: $H_\ell(x) = (-1)^\ell \exp(x^2)(d^\ell/dx^\ell)[\exp(-x^2)]$ and

$$ \varphi_\ell(x) := c_\ell H_\ell(x) \exp(-x^2/2), \quad c_\ell := (2^\ell \ell! \sqrt{\pi})^{-1/2}. \qquad (3.2) $$

The Hermite basis has the following two properties which will be used in the proofs. The Fourier transform verifies

$$ \varphi_{\ell}^{*} = \sqrt{2\pi i^{\ell}} \varphi_{\ell}, \qquad (3.3) $$

and we have according to Askey and Wainger (1965) that

$$ |\varphi_\ell(x)| \le \kappa_\varphi \exp(-\xi x^2), \quad |x| \ge \sqrt{2\ell + 1}, \quad \kappa_\varphi > 0. \qquad (3.4) $$

See Sacko (2020) for more details on projection estimation with the Hermite basis. Using this basis, we build a basis on $\mathbb{L}^2(\mathbb{R}^d)$. That is, for $\ell = (\ell_1, \dots, \ell_d) \in \mathbb{N}_0^d$, define

$$ \phi_\ell(\mathbf{x}) = \phi_{\ell_1, \dots, \ell_d}(\mathbf{x}) = (\varphi_{\ell_1} \otimes \dots \otimes \varphi_{\ell_d})(\mathbf{x}) = \varphi_{\ell_1}(x_1) \cdots \varphi_{\ell_d}(x_d). \qquad (3.5) $$
### 3.3 Convolution of the densities of $p$ independent multivariate variables

We now extend the projection estimator proposed in Section 3.1 to the multivariate multiple density case. Assume that $g^{(p)} \in L^2(\mathbb{R}^d)$. Given the Hermite orthonormal basis $(\phi_\ell)_{\ell \ge \mathbb{N}_0^d}$ of $L^2(\mathbb{R}^d)$, we can write $g^{(p)} = \sum_{\ell_1 \ge 0} \cdots \sum_{\ell_d \ge 0} a_{\ell_1, \dots, \ell_d} \phi_{\ell_1, \dots, \ell_d} = \sum_{\ell \ge 0} a_\ell \phi_\ell$ with $a_\ell = a_{\ell_1, \dots, \ell_d} := \langle g^{(p)}, \phi_{\ell_1, \dots, \ell_d} \rangle$. Let $g_m^{(p)} = \sum_{\ell=0}^{m-1} a_\ell \phi_\ell$ denote the orthogonal projection of $g^{(p)}$ on

$$S_m := \text{vec}(\{\phi_\ell | \forall k \in \llbracket d \rrbracket : \ell_k \le m - 1\}).$$

For $m = (m_1, \dots, m_d) \in \mathbb{N}_0^d$, this motivates analogously to Equation (3.1) the estimator

$$\widehat{g_m^{(p)}} := \sum_{\ell=0}^{m-1} \hat{a}_\ell \phi_\ell \qquad (3.6)$$

$$\hat{a}_{\ell} := \frac{1}{n_1 \cdots n_p} \sum_{k_1 \in \llbracket n_1 \rrbracket} \cdots \sum_{k_p \in \llbracket n_p \rrbracket} \phi_{\ell} (X_{k_1}^{(1)} + \cdots + X_{k_p}^{(p)}) = \frac{1}{(2\pi)^d} \langle \widehat{f}_1^* \cdots \widehat{f}_p^*, \phi_{\ell}^* \rangle. \qquad (3.7)$$

This defines an unbiased estimator of $g_m^{(p)}$. The general requirement is that the basis is such that the Fourier transform $\phi_\ell^*$ exists and is integrable. To consider the quadratic risk, the functions $\phi_\ell^*$ should also be square-integrable. These properties hold for the chosen Hermite basis. We show the following upper bound.

**Proposition 3.2 (Bound on the integrated quadratic risk - Projection estimator).** For some $C_0 > 0$ let $m \in \mathbb{N}^d$ satisfy

$$\frac{1}{\xi} \log \left( \prod_{j \in \llbracket p \rrbracket} n_j \right) \le m_k, \forall k \in \llbracket d \rrbracket \quad \text{and} \quad \frac{\prod_{k \in \llbracket d \rrbracket} \sqrt{2m_k + 1}}{\prod_{j \in \llbracket p \rrbracket} n_j} \le C_0. \qquad (3.8)$$

If $f_j^* \in L^1(\mathbb{R}^d) \cap L^2(\mathbb{R}^d)$ for all $j \in \llbracket p \rrbracket$, then, the estimator $\widehat{g_m}^{(p)}$ of $g^{(p)}$ defined in (3.6) with $\phi_\ell$ as in (3.2)-(3.5), satisfies

$$\mathrm{E} \left[ \| \widehat{g_m}^{(p)} - g^{(p)} \|^2 \right] \le \| g_m^{(p)} - g^{(p)} \|^2 + C_1 \frac{\prod_{k \in \llbracket d \rrbracket} \sqrt{2m_k + 1}}{\prod_{j \in \llbracket p \rrbracket} n_j} + C_2 \sum_{j \in \llbracket p \rrbracket} \frac{1}{n_j},$$

with $C_1 = \frac{8}{(2\pi)^d}$ and $C_2 = 4p3^{2(p-2)} \max_{j \in \llbracket p \rrbracket} \|f_j^*\|^2 + \frac{C_0 2^d 6^{2(d-1)} \kappa_\varphi^2}{\xi^2}$, with $\kappa_\varphi$ and $\xi$ defined in Equation (3.4).

Condition (3.8) is composed of two bounds on $m$ which are rather weak. The upper bound part ensures the boundedness of the variance and of cross terms. The lower bound part set on the $m_k$'s imposes minimal values more that log of the sample sizes and is not a strong restriction; in particular, it would allow for the optimal choice of $m$, see Remark 3.3.

**Remark 3.3.** Analogously to the discussion of convergence rates in Section 2.3, it is also possible to consider $d$-dimensional Hermite-Sobolev spaces as an extension of such one-dimensional spaces defined and studied in Bongioanni and Torrea (2006), introducing a condition on the coefficient leading to

$$\| g_m^{(p)} - g^{(p)} \|^2 \le C(m_1^{-\alpha_1} + \dots + m_d^{-\alpha_d}).$$

As for Proposition 3.2, choosing $\sqrt{m_{k,\text{opt}}}$, we get the same rates as in Section 2.3. However, the link between the behavior of the coefficients and functional regularity is not proven in the multivariate setting.
## 3.4 Model selection

In this section, we consider the case where all sample sizes are the same, i.e. $n_1 = \dots = n_p = n \in \mathbb{N}$. The results are analogous to assuming that the sample sizes are of the same order up to changes in constants. We continue to use the Hermite basis $(\phi_\ell)_{\ell \ge 0}$ of $L^2(\mathbb{R}^d)$ and the estimator $\hat{g}_m^{(p)}$ given by (3.6). For $t \in S_m$, we set

$$ \gamma_n(t) = \|t\|^2 - \frac{2}{n^p} \sum_{k_1, \dots, k_p \in \llbracket n \rrbracket} t(\mathbf{X}_{k_1}^{(1)} + \dots + \mathbf{X}_{k_p}^{(p)}), $$

and define

$$ \widehat{m} = \arg \min_{\mathbf{m} \in \mathcal{M}_n} \{\gamma_n(\widehat{g}_{\mathbf{m}}^{(p)}) + \text{pen}(\mathbf{m})\} \quad (3.9) $$

where $\mathcal{M}_n \subset \mathbb{N}^d$, and for $\mathbf{m} \in \mathbb{N}^d$ the penalty sequence $(\text{pen}(\mathbf{m}))_{\mathbf{m} \in \mathbb{N}^d}$ is given by

$$ \text{pen}(\mathbf{m}) = \kappa \frac{\prod_{k \in \llbracket d \rrbracket} \sqrt{2m_k + 1} \log^d(n)}{n^d}. $$

We take $\mathcal{M}_n$ as the collection of models such that, for $c_0 := (d^2 - p + 1)/(2\xi)$,

$$ \mathcal{M}_n = \left\{ \mathbf{m} \in \mathbb{N}^d : \forall k \in \llbracket d \rrbracket, m_k \ge c_0 \log(n) \text{ and } \prod_{k \in \llbracket d \rrbracket} \sqrt{2m_k + 1} \le n^d \right\}. \quad (3.10) $$

Note the similarity of the collection to the assumption made on the dimension parameter in Equation (3.8) and see the comments herewith. It is slightly reinforced taking in account that the sample sizes are assumed to be equal in this section.

**Theorem 3.4.** Assume that $f_j^* \in L^1(\mathbb{R}^d) \cap L^2(\mathbb{R}^d)$ for all $j \in \llbracket p \rrbracket$, and consider the estimator $\widehat{g}_{\widehat{m}}^{(p)}$ of $g^{(p)}$ defined by (3.6)-(3.9) with $\phi_\ell$ as in (3.2)-(3.5). Then, there exists $\kappa_0$ such that for any $\kappa \ge \kappa_0$,

$$ \mathbb{E} \left[ \| \widehat{g}_{\widehat{m}}^{(p)} - g^{(p)} \|^2 \right] \le \inf_{\mathbf{m} \in \mathcal{M}_n} \left( 9 \mathbb{E} \left[ \| \widehat{g}_m^{(p)} - g^{(p)} \|^2 \right] + 2 \text{pen}(\mathbf{m}) \right) + \frac{C}{n}, $$

where $\kappa_0 = 11 \frac{2^d (4pd - 1)^p}{\pi^d}$ suits, and C is a constant depending on p, d, $\xi$, $\kappa_\varphi$ and $\|f_j\|$.

Proposition 3.2 gives a decomposition of $\mathbb{E}[\|\widehat{g}_m^{(p)} - g^{(p)}\|^2]$. Therefore, Theorem 3.4 shows that, up to a multiplicative constant, the risk of the final estimator minimizes over all $\mathbf{m}$'s in the collection, the squared-bias $\|\widehat{g}_m^{(p)} - g^{(p)}\|^2$ plus the penalty which has the order of the variance up to a $\log^d(n)$ term. However, this log-loss does not change the rates.

Note that the constant $\kappa$ has to be calibrated and to be fixed, see Baudry et al. (2012), and the theoretical value $\kappa_0$ is likely to be too large in practice.

## 4 Simulation study

Our simulation experiments deal with the cases where $d = 1$ and $p = 2$ or $p = 3$. First, we consider the following examples for which we can analytically derive the convolution to be able to evaluate the procedures:
<figure><img src="image_4.png" /><figcaption>Figure 1: Kernel estimator with $n = 200$ for examples (i)-(iii) and $n = 60$ for (iv); Left column: Plot of $\hat{g}_h$ for $h \in \mathcal{H}$ (green lines) and true function $g$ (black line); Right column: Estimated MISE (black line) using $M = 50$ Monte Carlo simulations and estimated values $\hat{h}$ using the PCO-method (green histograms).</figcaption></figure>
(i) Two independent Gaussian distributed random variables.

(ii) Two independent Claw distributed random variables.

(iii) Two independent Unif[-1, 1] random variables.

(iv) Three independent Unif[-1, 1] random variables.

The first three are the examples given in Chesneau et al. (2013). For all cases we consider the kernel estimator defined in Equation (2.2) and the projection estimator defined in Equation (3.6). We choose $L(x) = 3n_1(x) - 3n_2(x) + n_3(x)$, a kernel of order 5. Here, $n_j$ denotes the density of a centered Gaussian random variable with variance $j \in \mathbb{N}$. Note that one can build kernels of higher order; see Comte and Marie (2020) for details. For the projection method, we use the Hermite basis throughout the simulations.

Figure 1 shows the behavior of kernel estimators. In the left column of Figure 1, kernel estimators corresponding to 20 values of the bandwidth $h \in \mathcal{H}$ are depicted, where $\mathcal{H}$ denotes the set of 20 equidistant points starting from 0.005 and ending with 1.25. It visualizes the necessity of a data-driven choice of the parameter $h$ to obtain a trade-off between over-fitting (dark green lines) and a too rough estimator (light green lines). In the right column of Figure 1, the curve represents the MISE computed over 50 simulations as a function of $h$. One can see that in the given examples, there seem to be intervals of appropriate choices for $h$ of different sizes depending on the example; furthermore, the theory suggests choosing the smallest bandwidth $h = 1/n$, but this choice is revealed to be too small in these examples (first point of each black line). The histograms depict the choices for $h$ made with the PCO-method described in Section 2.4. They show that the data-driven method performs quite well on these examples. Values of bandwidth outside the optimal interval occur only for the Claw-distribution example, for two out of 50 Monte Carlo iterations.

Figure 2 depicts similar results for the projection estimator and the model selection method proposed in Section 3.4. Here, apart from finding a sufficiently good estimator, the goal is to find a parsimonious decomposition of $g^{(p)}$, i.e. a small number $m$ of coefficients for estimating the function $g^{(p)}$ appropriately well. For examples (i), (ii) and (iv), this goal seems to be achieved.

Figure 3 and Figure 4 compare the performance of the projection and the kernel method. The slightly better results for the projection method over all examples and the significantly shorter computation time motivated the focus of the theoretical consideration in Section 3.4 on the selection method.

<table><thead><tr><th rowspan="2">n</th><th colspan="3">(i) Gauss + Gauss</th><th colspan="3">(ii) Claw + Claw</th><th colspan="3">(iii) Unif + Unif</th></tr><tr><th>K</th><th>P</th><th>A</th><th>K</th><th>P</th><th>A</th><th>K</th><th>P</th><th>A</th></tr></thead><tbody><tr><td>50</td><td>4.04</td><td>0.18</td><td>0.05</td><td>4.15</td><td>0.36</td><td>0.27</td><td>4.56</td><td>0.24</td><td>2.01</td></tr><tr><td>100</td><td>0.45</td><td>0.07</td><td>0.01</td><td>0.48</td><td>0.23</td><td>0.08</td><td>0.62</td><td>0.15</td><td>1.98</td></tr><tr><td>200</td><td>0.08</td><td>0.04</td><td>0.00</td><td>0.15</td><td>0.08</td><td>0.02</td><td>0.13</td><td>0.05</td><td>1.98</td></tr></tbody></table>

<table><caption>Table 1: Errors for Examples (i)–(iii), a Gamma(2, 1) and a Gamma(1, 1) distributed random variable (Gam + Gam) $\mathcal{T} = [0, 20]$ and two Exp(1) distributed random variables (Exp + Exp) $\mathcal{T} = [0, 10]$; using kernel (K), projection (P), and numerical approximation (A) methods.</caption><thead><tr><th rowspan="2">n</th><th colspan="3">Gam + Gam</th><th colspan="3">Exp + Exp</th></tr><tr><th>K</th><th>P</th><th>A</th><th>K</th><th>P</th><th>A</th></tr></thead><tbody><tr><td>50</td><td>10.37</td><td>2.35</td><td>4.92</td><td>4.77</td><td>6.53</td><td>0.24</td></tr><tr><td>100</td><td>6.21</td><td>1.82</td><td>4.78</td><td>0.83</td><td>5.85</td><td>0.10</td></tr><tr><td>200</td><td>5.53</td><td>1.54</td><td>4.72</td><td>0.21</td><td>5.68</td><td>0.06</td></tr></tbody></table>

In Table 1 and Table 2 the performance of the kernel (K) and projection method (P) is compared to the
<figure><img src="image_5.png" /><figcaption>Figure 2: Projection estimator with $n = 200$ for examples (i)-(iii) and $n = 60$ for (iv); Left column: Plot of $\hat{g}_m$ for $m \in [20], [50]$ and $[150]$, respectively (green lines) and true function $g$ (black line); Right column: Estimated MISE plot with $M = 50$ Monte Carlo simulations and estimated values $\hat{m}$ using the selection method (green histograms).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 3: Plots of $\hat{g}_{\hat{h}}$ using the PCO-method (green lines, first row) and $\hat{g}_{\hat{m}}$ using selection method (second line) for $M = 50$ Monte Carlo simulations and $n = 200$ for Examples (i)-(iii) and $n = 60$ for (iv).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Estimated integrated squared error choosing $h$ with PCO-Method and projection method with Hermite basis for $M = 50$ Monte Carlo simulations; $n = 200$ for examples (i)-(iii) and $n = 60$ for example (iv).</figcaption></figure>
<table><caption>Table 2: Errors for three Unif[-1, 1]-distributed random variables estimated on $\mathcal{T} = [-5, 5]$ and three Exp(1) distributed random variables on $\mathcal{T} = [0, 10]$; using kernel (K), projection (P), and numerical approximation (A) methods with sample sizes $n = 10, 30, 60, 90$.</caption><thead><tr><th rowspan="2">n</th><th colspan="3">(iv) Unif[-1, 1]</th><th colspan="3">Exp(1)</th></tr><tr><th>K</th><th>P</th><th>A</th><th>K</th><th>P</th><th>A</th></tr></thead><tbody><tr><td>10</td><td>2.32</td><td>0.48</td><td>23.13</td><td>1.76</td><td>3.95</td><td>4.11</td></tr><tr><td>30</td><td>0.78</td><td>0.11</td><td>1.83</td><td>0.58</td><td>2.56</td><td>0.84</td></tr><tr><td>60</td><td>0.34</td><td>0.04</td><td>0.43</td><td>0.33</td><td>2.23</td><td>0.50</td></tr><tr><td>90</td><td>0.26</td><td>0.02</td><td>0.19</td><td>0.22</td><td>2.15</td><td>0.45</td></tr></tbody></table>

numerical approximation (A) for different choices of the sample size $n$. Here, for function $\hat{f}$ evaluated at points $t_1, \dots, t_N$ with interval steps $\Delta_t$ approximating function $f$ we consider the error

$$100 * \sum_{i=1}^{N} \Delta_t (\hat{f}(t_i) - f(t_i))^2,$$

where we scale by 100 for better readability. For all examples, we chose $N = 100$ and $\{t_1, \dots, t_N\}$ as the set of equidistant points on $\mathcal{T}$ depending on the example. For examples (i), (ii) and (iv) we set $\mathcal{T} = [-5, 5]$ and for (iii) we set $\mathcal{T} = [-2.5, 2.5]$. When considering estimation, we additionally take the mean over 50 Monte Carlo simulations. For the numerical approximation (A), the Fast Fourier Transform was used (in R implemented with the function `fft`) with $n$ evaluation points.

As expected, increasing the sample size $n$ reduces the error. However, there is not one method out-performing the others over all examples. The numerical approximation has smaller errors, e.g., for two exponentially distributed random variables. On the other hand, for the case of three exponentially distributed random variables, the kernel method performs best. Conversely, for two uniform distributed random variables the projection method appears to be most fitting. In addition, the relationships might change with increasing $n$. For instance, for three uniform distributed random variables, the kernel methods seems better than the numerical approximation for small $n$ but for large values the roles switch.

Lastly, in Table 3 the results are given for two examples where the analytical derivation of the convolution is not easy to do. Again, they indicate a better performance of the projection method.

<table><caption>Table 3: Errors for estimating the distribution of the sum of a Unif[-1, 1] distributed variable and either a standard Gaussian or a Claw-distributed random variable. For both examples, the function was estimated on $\mathcal{T} = [-5, 5]$ and for the error FFT was used as the ground truth.</caption><thead><tr><th rowspan="2">n</th><th colspan="2">Unif + Gaussian</th><th colspan="2">Unif + Claw</th></tr><tr><th>K</th><th>P</th><th>K</th><th>P</th></tr></thead><tbody><tr><td>50</td><td>4.38</td><td>0.26</td><td>4.86</td><td>0.25</td></tr><tr><td>100</td><td>0.49</td><td>0.13</td><td>0.77</td><td>0.18</td></tr><tr><td>200</td><td>0.17</td><td>0.11</td><td>0.19</td><td>0.13</td></tr></tbody></table>

**Concluding Remarks** In this work, we have studied a kernel and a projection estimator for the estimation of the density of the sum of $p$ different $d$-dimensional independent random vectors, proving upper bounds for the integrated squared error, discussing the convergence rates and proposing data-driven bandwidths or model selection methods.
In the kernel setting, the convolution of kernel density estimators provides an estimator of the convolution of the densities; this is not true for projection estimators. This is due to the fact that, when $L^2$-orthonormal bases are considered, the convolution of two basis functions does not have a simple development in the basis, whereas the convolution of two kernels is a kernel. However, this would occur for the Laguerre basis (see, e.g. Comte and Lacour (2023)), which is an orthonormal basis on $L^2(\mathbb{R}^+)$: the convolution of two Laguerre functions can be written as a linear combination of two other Laguerre functions. Thus, it would be interesting to propose a specific study for the estimation of the density of the sum of nonnegative random variables, as one could consider the estimator resulting from the convolution of two projection estimators.

Further work of interest may include the consideration of measurement errors on the observations, which would include several parameters in the rates of convergence; in particular, the parameters associated with the rate of decay of the characteristic function of the noise would play a role.

# A Proofs

## A.1 Proofs of Section 2

For the proofs of Section 2.1 and Section 3.1 we first see that for a suitable function $F$ we have that

$$ 
\begin{aligned}
& \operatorname{var}\left(\frac{1}{n^2} \sum_{i \in \llbracket n \rrbracket} \sum_{k \in \llbracket n \rrbracket} F(X_i^{(1)} + X_k^{(2)})\right) \\
& = \frac{1}{n^4} \sum_{i_1, i_2 \in \llbracket n \rrbracket} \sum_{k_1, k_2 \in \llbracket n \rrbracket} \operatorname{cov}\left(F(X_{i_1}^{(1)} + X_{k_1}^{(2)}), F(X_{i_2}^{(1)} + X_{k_2}^{(2)})\right) \\
& = \frac{1}{n^4} \left( \sum_{i \in \llbracket n \rrbracket} \sum_{\substack{k_1, k_2 \in \llbracket n \rrbracket \\ k_1 \neq k_2}} \operatorname{cov}\left(F(X_i^{(1)} + X_{k_1}^{(2)}), F(X_i^{(1)} + X_{k_2}^{(2)})\right) \right. \\
& \qquad \left. + \sum_{k \in \llbracket n \rrbracket} \sum_{\substack{i_1, i_2 \in \llbracket n \rrbracket \\ i_1 \neq i_2}} \operatorname{cov}\left(F(X_{i_1}^{(1)} + X_{k}^{(2)}), F(X_{i_2}^{(1)} + X_{k}^{(2)})\right) \sum_{i \in \llbracket n \rrbracket} \sum_{k \in \llbracket n \rrbracket} \operatorname{var}\left(F(X_i^{(1)} + X_k^{(2)})\right) \right).
\end{aligned}
$$

Due to the independence assumption on the observations, we further get

$$ \operatorname{var}\left(\frac{1}{n^2} \sum_{i \in \llbracket n \rrbracket} \sum_{k \in \llbracket n \rrbracket} F(X_i^{(1)} + X_k^{(2)})\right) = T_1 + T_2 + T_3, \qquad (A.1) $$

with

$$ T_1 := \frac{1}{n} \frac{n-1}{n} \operatorname{cov}\left(F(X_1^{(1)} + X_1^{(2)}), F(X_1^{(1)} + X_2^{(2)})\right), $$

$$ T_2 := \frac{1}{n} \frac{n-1}{n} \operatorname{cov}\left(F(X_1^{(1)} + X_1^{(2)}), F(X_2^{(1)} + X_1^{(2)})\right), $$

$$ T_3 := \frac{1}{n^2} \operatorname{var}(F(X_1^{(1)} + X_1^{(2)})). $$

Choosing $F$ correspondingly and calculating $T_1, T_2, T_3$, we show Proposition 2.1, Proposition 2.2 and Proposition 3.1.
*Proof of Proposition 2.1.* First we see with the elementary triangle inequality that

$$ \mathbb{E} [\|\hat{g}_h - g\|^2] \le \|g_h - g\|^2 + \mathbb{E} [\|\hat{g}_h - g_h\|^2]. $$

We have that

$$ \begin{aligned} \mathbb{E} [\|\hat{g}_h - g_h\|^2] &= \mathbb{E} \left[ \left| \int \left| \frac{1}{n^2} \sum_{k_1 \in \llbracket n \rrbracket} \sum_{k_2 \in \llbracket n \rrbracket} L_h (X_{k_1}^{(1)} + X_{k_2}^{(2)} - x) - L_h * g(x) \right|^2 dx \right] \\ &= \int \operatorname{var} \left( \frac{1}{n^2} \sum_{k_1 \in \llbracket n \rrbracket} \sum_{k_2 \in \llbracket n \rrbracket} L_h (X_{k_1}^{(1)} + X_{k_2}^{(2)} - x) \right) dx. \end{aligned} $$

With $F(y) = L_h(y - x)$ we use Equation (A.1) and bound $T_1, T_2, T_3$ integrated over $x$. For $T_3$, we see that

$$ \begin{aligned} \int \frac{1}{n^2} \operatorname{var}(L_h(X_1^{(1)} + X_1^{(2)} - x)) dx &\le \frac{1}{n^2 h^2} \int \int g(y) \left| L\left(\frac{y-x}{h}\right) \right|^2 dy dx \\ &= \frac{1}{n^2 h} \int g(y) dy \int |L(x)|^2 dx = \frac{\|L\|^2}{n^2 h}. \end{aligned} $$

For $T_1$, we have using Young's inequality that

$$ \begin{aligned} & \frac{1}{n} \int \operatorname{cov} (L_h(X_1^{(1)} + X_1^{(2)} - x), L_h(X_1^{(1)} + X_2^{(2)} - x)) dx \\ &= \frac{1}{n} \iint \iint L_h(u+v-x)L_h(u+z-x)f_1(u)f_2(v)f_2(z)dudvdzdx \\ & \quad - \frac{1}{n} \int \left( \mathbb{E}[L_h(X_1^{(1)} + X_1^{(2)} - x)] \right)^2 dx \\ & \le \frac{1}{n} \iint \left( \int L_h(u+v-x)f_2(v)dv \right)^2 f_1(u)dudx \\ & \le \frac{1}{n} \int f_1(u)du \int (L_h * f_2(x-u))^2 dx \le \frac{1}{n} \|L_h\|_1^2 \|f_2\|^2 = \frac{1}{n} \|L\|_1^2 \|f_2\|^2 \end{aligned} $$

Due to the symmetry of $T_1$ and $T_2$, we get for $T_2$ the bound $\frac{1}{n}\|L\|_1^2\|f_1\|^2$. It follows that

$$ \mathbb{E} [\|\hat{g}_h - g_h\|^2] \le \frac{\|L\|^2}{n^2 h} + \|L\|_1^2 \frac{\|f_1\|^2 + \|f_2\|^2}{n} $$

which concludes the proof.

*Proof of Proposition 2.2.* We have that

$$ \mathbb{E} [((\hat{g}_h(x) - g(x))^2)] = \mathbb{E} [((\hat{g}_h(x) - g_h(x))^2) + (g_h(x) - g(x))^2]. $$

Analogously to the proof of Proposition 2.1 we have that

$$ \mathbb{E} [((\hat{g}_h(x) - g_h(x))^2)] = \operatorname{var} \left( \frac{1}{n^2} \sum_{k_1 \in \llbracket n \rrbracket} \sum_{k_2 \in \llbracket n \rrbracket} L_h (X_{k_1}^{(1)} + X_{k_2}^{(2)} - x) \right). $$
So we can use again Equation (A.1) with $F(y) = L_h(y-x)$ but we do not integrate the terms $T_1, T_2, T_3$ over $x$. We get for $T_3$ that

$$ \frac{1}{n^2} \mathrm{var}(L_h(X_1^{(1)} + X_1^{(2)} - x)) \leq \int g(y) \left| L\left(\frac{y-x}{h}\right) \right|^2 dy \leq \frac{1}{n^2 h} \|g\|_{\infty} \|L\|^2. $$

For $T_1$, it holds using again Young's inequality

$$ \begin{aligned} & \frac{1}{n} \operatorname{cov} (L_h(X_1^{(1)} + X_1^{(2)} - x), L_h(X_1^{(1)} + X_2^{(2)} - x)) \\ & = \frac{1}{n} \int \left( \int L_h(u+v-x) f_2(v) dv \right)^2 f_1(u) du \\ & \leq \frac{1}{n} \|f_1\|_{\infty} \|L_h * f_2(x-\cdot)\|^2 \leq \frac{1}{n} \|f_1\|_{\infty} \|L\|_1^2 \|f_2\|^2. \end{aligned} $$

Analogously, $T_2$ can be upper bounded by $\frac{1}{n}\|f_2\|_{\infty}\|L\|_1^2\|f_1\|^2$. Consequently, it holds that

$$ \mathbb{E} \left[ \left( \hat{g}_h(x) - g_h(x) \right)^2 \right] \leq \frac{\|g\|_{\infty} \|L\|^2}{n^2 h} + \|L\|_1^2 \frac{\|f_1\|_{\infty} \|f_2\|^2 + \|f_2\|_{\infty} \|f_1\|^2}{n} $$

which concludes the proof. $\square$

*Proof of equivalence of Equation (2.2) and Equation (2.3).* If $\hat{g}_h^{(p)} \in \mathbb{L}^1(\mathbb{R}^d)$, we can rewrite Equation (2.2) taking the Fourier inverse of its Fourier transform, i.e.

$$ \begin{aligned} \hat{g}_h^{(p)}(\boldsymbol{x}) &= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{-i\langle \boldsymbol{u}, \boldsymbol{x} \rangle} (\hat{g}_h^{(p)})^* (\boldsymbol{u}) d\boldsymbol{u} \\ &= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{-i\langle \boldsymbol{u}, \boldsymbol{x} \rangle} \frac{1}{h_1 \cdots h_d n_1 \cdots n_p} \\ & \qquad \cdot \sum_{j_1 \in [n_1]} \cdots \sum_{j_p \in [n_p]} \int_{\mathbb{R}^d} e^{i\langle \boldsymbol{u}, \boldsymbol{y} \rangle} L \left( \frac{\boldsymbol{X}_{j_1}^{(1)} + \cdots + \boldsymbol{X}_{j_p}^{(p)} - \boldsymbol{y}}{h} \right) dy d\boldsymbol{u} \\ &= \frac{1}{(2\pi)^d} \int e^{-i\langle \boldsymbol{u}, \boldsymbol{x} \rangle} \frac{1}{n_1 \cdots n_p} \sum_{j_1 \in [n_1]} \cdots \sum_{j_p \in [n_p]} \int L(\boldsymbol{w}) e^{i\langle \boldsymbol{u}, \boldsymbol{X}_{j_1}^{(1)} + \cdots + \boldsymbol{X}_{j_p}^{(p)} - h\boldsymbol{w} \rangle} dw d\boldsymbol{u} \\ &= \frac{1}{(2\pi)^d} \int e^{-i\langle \boldsymbol{u}, \boldsymbol{x} \rangle} L^*(h\boldsymbol{u}) \prod_{j \in [p]} \hat{f}_j^*(\boldsymbol{u}) d\boldsymbol{u}. \end{aligned} $$

*Proof of Proposition 2.4.* First we see that with the elementary triangle inequality that

$$ \mathbb{E} \left[ \| \hat{g}_h^{(p)} - g^{(p)} \|^2 \right] \leq \| g_h^{(p)} - g^{(p)} \|^2 + \mathbb{E} \left[ \| \hat{g}_h^{(p)} - g_h^{(p)} \|^2 \right]. $$

Having representation Equation (2.3) in mind, Parseval's theorem gives for the variance term

$$ \mathbb{E} \left[ \| \hat{g}_h^{(p)} - g_h^{(p)} \|^2 \right] = \frac{1}{(2\pi)^d} \mathbb{E} \left[ \int_{\mathbb{R}^d} \left| L^*(h\boldsymbol{t}) \prod_{j \in [p]} \hat{f}_j^*(\boldsymbol{t}) - \prod_{j \in [p]} f_j^*(\boldsymbol{t}) \right|^2 dt \right] $$
Using Lemma B.1 and independence between the random variables $\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(p)}$ we get that

$$
\begin{aligned}
\mathrm{E} \left[ \left| \prod_{j \in [p]} \widehat{f}_j^*(\mathbf{t}) - \prod_{j \in [p]} f_j^*(\mathbf{t}) \right|^2 \right] &\le 2 \prod_{j \in [p]} \mathrm{E} \left[ \left| \frac{1}{n_j} \sum_{k \in [n_j]} e^{i \langle \mathbf{t}, \mathbf{x}_k^{(j)} \rangle} - f_j^*(\mathbf{t}) \right|^2 \right] \\
&\quad + 2p \cdot 3^{2(p-2)} \sum_{j \in [p]} |f_{j+1}^*(\mathbf{t})|^2 \mathrm{E} \left[ \left| \frac{1}{n_j} \sum_{k \in [n_j]} e^{i \langle \mathbf{t}, \mathbf{x}_k^{(j)} \rangle} - f_j^*(\mathbf{t}) \right|^2 \right].
\end{aligned}
$$

Further, we see that

$$
\mathrm{E} \left[ \left| \frac{1}{n_j} \sum_{k \in [n_j]} e^{i \langle \mathbf{t}, \mathbf{X}_k^{(j)} \rangle} - f_j^*(\mathbf{t}) \right|^2 \right] = \frac{1}{n_j} \mathrm{var} \left( e^{i \langle \mathbf{t}, \mathbf{X}_k^{(j)} \rangle} \right) \le \frac{1}{n_j} \quad (A.2)
$$

with which follows

$$
\mathrm{E} \left[ \left| \prod_{j \in [p]} \left( \frac{1}{n_j} \sum_{k \in [n_j]} e^{i \langle \mathbf{t}, \mathbf{X}_k^{(j)} \rangle} \right) - \prod_{j \in [p]} f_j^*(\mathbf{t}) \right|^2 \right] \le \frac{2}{n_1 \cdots n_p} + 2p \cdot 3^{2(p-2)} \sum_{j \in [p]} \frac{|f_{j+1}^*(\mathbf{t})|^2}{n_j}. \quad (A.3)
$$

In addition, we have for any $j \in [p]$ using again Parseval's theorem that

$$
\int_{\mathbb{R}^d} |L^*(\mathbf{h}\mathbf{t}) f_j^*(\mathbf{t})|^2 d\mathbf{t} \le \|L^*\|_{\infty}^2 \|f_j^*\|^2 \le \|L^*\|_{\infty}^2 \max_{j \in [p]} \|f_j^*\|^2.
$$

Consequently,

$$
\mathrm{E} \left[ \| \widehat{g}_{\mathbf{h}}^{(p)} - g_{\mathbf{h}}^{(p)} \|^2 \right] \le \frac{2}{h_1 \cdots h_d n_1 \cdots n_p} \int_{\mathbb{R}^d} |L^*(\mathbf{t})|^2 d\mathbf{t} + C \sum_{j \in [p]} \frac{1}{n_j},
$$

$$
\text{with } C := 2p3^{2(p-2)} \|L^*\|_{\infty}^2 \max_{j \in [p]} \|f_j^*\|^2.
$$

□

**A.2 Proofs of Section 3**

*Proof of Proposition 3.1.* Using the squared bias variance decomposition it remains to show an upper bound for the variance. First, we have that

$$
\mathrm{E} \left[ \| \widehat{g}_m - g_m \|^2 \right] = \sum_{\ell=0}^{m-1} \mathrm{E} \left[ (\hat{a}_\ell - a_\ell)^2 \right] = \sum_{\ell=0}^{m-1} \mathrm{var} \left( \frac{1}{n^2} \sum_{i \in [n]} \sum_{k \in [n]} \varphi_\ell (X_i^{(1)} + X_k^{(2)}) \right)
$$

Consequently, using Equation (A.1) with $F = \varphi_\ell$, we have that

$$
\begin{aligned}
\mathrm{E} \left[ \| \widehat{g}_m - g_m \|^2 \right] \\
&= \frac{1}{n} \sum_{\ell=0}^{m-1} \frac{n-1}{n} \left( \begin{aligned}[t]
& \operatorname{cov} \left( \varphi_\ell (X_1^{(1)} + X_1^{(2)}), \varphi_\ell (X_1^{(1)} + X_2^{(2)}) \right) \\
& + \operatorname{cov} \left( \varphi_\ell (X_1^{(1)} + X_1^{(2)}), \varphi_\ell (X_2^{(1)} + X_1^{(2)}) \right)
\end{aligned} \right) \\
&\quad + \frac{1}{n^2} \sum_{\ell=0}^{m-1} \operatorname{var}(\varphi_\ell (X_1^{(1)} + X_1^{(2)})) = T_1 + T_2 + T_3.
\end{aligned}
$$
For the third summand it holds using that the basis is orthonormal and that $g$ is a density

$$T_3 \le \frac{1}{n^2} \sum_{\ell=0}^{m-1} \mathbb{E}[\varphi_{\ell}^2(X_1^{(1)} + X_1^{(2)})] = \frac{1}{n^2} \sum_{\ell=0}^{m-1} \int \varphi_{\ell}^2(x)g(x)dx \le \frac{1}{n^2} (L(m) \wedge \|g\|_{\infty}m).$$

For the first summand, we see that

$$\begin{align*}
\operatorname{cov}(\varphi_{\ell}(X_1^{(1)} + X_1^{(2)}), \varphi_{\ell}(X_1^{(1)} + X_2^{(2)}))
&= \mathbb{E}[\varphi_{\ell}(X_1^{(1)} + X_1^{(2)})\varphi_{\ell}(X_1^{(1)} + X_2^{(2)})] - (\mathbb{E}[\varphi_{\ell}(X_1^{(1)} + X_2^{(2)})])^2 \\
&= \int \left(\int \varphi_{\ell}(x+y)f_2(y)dy\right) \left(\int \varphi_{\ell}(x+z)f_2(z)dz\right) f_1(x)dx - \left(\int \varphi_{\ell}(x)g(x)dx\right)^2 \\
&= \int \left(\int \varphi_{\ell}(x)f_2(u-x)dy\right)^2 f_1(x)dx - \left(\int \varphi_{\ell}(x)g(x)dx\right)^2.
\end{align*}$$

Denoting by $\Pi_{S_m}$ the orthogonal projection on $S_m$ we get that

$$\begin{align*}
\sum_{\ell=0}^{m-1} \operatorname{cov}(\varphi_{\ell}(X_1^{(1)} + X_1^{(2)}), \varphi_{\ell}(X_1^{(1)} + X_2^{(2)}))
&\le \sum_{\ell=0}^{m-1} \int \left(\int \varphi_{\ell}(x)f_2(u-x)dy\right)^2 f_1(x)dx \\
&= \int \sum_{\ell=0}^{m-1} \langle \varphi_{\ell}, f_2(\cdot-x) \rangle^2 f_1(x)dx = \int \|\Pi_{S_m} f_2(\cdot-x)\|_2^2 f_1(x)dx \\
&\le \int \|f_2(\cdot-x)\|^2 f_1(x)dx = \int f_2^2(v)dv \int f_1(x)dx = \int f_2^2(v)dv.
\end{align*}$$

Analogously, one can bound $\sum_{\ell=0}^{m-1} \operatorname{cov}(\varphi_{\ell}(X_1^{(1)} + X_1^{(2)}), \varphi_{\ell}(X_2^{(1)} + X_1^{(2)})) \le \int f_1^2(v)dv$. With this, it follows that

$$T_1 + T_2 \le \frac{1}{n} \left( \int f_2^2(v)dv + \int f_1^2(v)dv \right)$$

which concludes the proof.

<figure><img src="image_8.png" /></figure>

*Proof of Proposition 3.2.* First we see that with the elementary triangle inequality that

$$\mathbb{E} \left[ \| \widehat{g_m^{(p)}} - g^{(p)} \|^2 \right] \le \| g_m^{(p)} - g^{(p)} \|^2 + \mathbb{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right].$$

Further, it holds that

$$\begin{align*}
\mathbb{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right]
&= \sum_{\ell=0}^{m-1} \mathbb{E} \left[ |\hat{a}_{\ell} - a_{\ell}|^2 \right] \\
&= \frac{1}{(2\pi)^{2d}} \sum_{\ell=0}^{m-1} \mathbb{E} \left[ \left| \langle \widehat{f}_1^* \cdots \widehat{f}_p^* - f_1^* \cdots f_p^*, \phi_{\ell}^* \rangle \right|^2 \right]
\end{align*}$$
Note that with Equation (3.3) we have

$$
\begin{align*}
\phi_{\ell}^{*} = \phi_{\ell_1, \dots, \ell_d}^{*}(\boldsymbol{x}) &= \phi_{\ell_1}^{*}(x_1) \cdots \phi_{\ell_d}^{*}(x_d) = (\sqrt{2\pi})^d i^{\ell_1 + \dots + \ell_d} \varphi_{\ell_1}(x_1) \cdots \varphi_{\ell_d}(x_d) \\
&= (\sqrt{2\pi})^d i^{\ell_1 + \dots + \ell_d} \phi_{\ell_1, \dots, \ell_d}(\boldsymbol{x}) = (\sqrt{2\pi})^d i^{\ell_1 + \dots + \ell_d} \phi_{\ell}(\boldsymbol{x})
\end{align*}
$$

It follows that

$$
\begin{align*}
\mathrm{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right] &\le \frac{1}{(2\pi)^d} \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \langle \widehat{f}_1^* \cdots \widehat{f}_p^* - f_1^* \cdots f_p^*, \phi_\ell \rangle \right|^2 \right] \\
&= \frac{1}{(2\pi)^d} \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \int_{\mathbb{R}^d} \left( \prod_{j \in [p]} \widehat{f}_j^*(t) - \prod_{j \in [p]} f_j^*(t) \right) \phi_\ell(t) dt \right|^2 \right].
\end{align*}
$$

Denote $\widehat{\Delta f}^*(t) := \prod_{j \in [p]} \widehat{f}_j^*(t) - \prod_{j \in [p]} f_j^*(t)$. We split the integral on the set $|t| \le \sqrt{2m+1}$ (understood component wise) and its complement and bound each summand separately. More precisely, we have that

$$
\mathrm{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right] \le \frac{2}{(2\pi)^d} (T_1 + T_2)
$$

with

$$
T_1 := \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \int_{|t| \le \sqrt{2m+1}} (\widehat{\Delta f}^*(t)) \phi_\ell(t) dt \right|^2 \right],
$$

$$
T_2 := \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \int_{\mathbb{R}^d} \mathbb{1}_{\{\exists k \in [d] : |t_k| > \sqrt{2m_k+1}\}} (\widehat{\Delta f}^*(t)) \phi_\ell(t) dt \right|^2 \right].
$$

For $T_1$ we see using Equation (A.3) that

$$
\begin{align*}
T_1 &= \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \langle \widehat{\Delta f}^* \mathbb{1}_{\{|\cdot| \le \sqrt{2m+1}\}}, \phi_\ell \rangle \right|^2 \right] \le \mathrm{E} \left[ \int_{\mathbb{R}^d} \left| \widehat{\Delta f}^*(t) \right|^2 \mathbb{1}_{\{|t| \le \sqrt{2m+1}\}} dt \right] \\
&\le \frac{4 \prod_{k \in [d]} \sqrt{2m_k + 1}}{\prod_{j \in [p]} n_j} + 2p \cdot 3^{2(p-2)} \max_{j \in [p]} \|f_j^*\|^2 \sum_{j \in [p]} \frac{1}{n_j}.
\end{align*}
$$

For $T_2$, denote

$$
(\widehat{\Delta f}^*)_1(t) := \prod_{j \in [p]} (\widehat{f}_j^*(t) - f_j^*(t)) \qquad (A.4)
$$

$$
(\widehat{\Delta f}^*)_2(t) := (\widehat{\Delta f}^*)(t) - (\widehat{\Delta f}^*)_1(t) \qquad (A.5)
$$

and use this to further split $T_2$. That is,

$$
T_2 \le 2 \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \int_{\mathbb{R}^d} \mathbb{1}_{\{\exists k \in [d] : |t_k| > \sqrt{2m_k+1}\}} (\widehat{\Delta f}^*(t))_1 \phi_\ell(t) dt \right|^2 \right] \\
+ 2 \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \int_{\mathbb{R}^d} \mathbb{1}_{\{\exists k \in [d] : |t_k| > \sqrt{2m_k+1}\}} (\widehat{\Delta f}^*(t))_2 \phi_\ell(t) dt \right|^2 \right]
$$
$$
\begin{align*}
T_2 &\le 2d \sum_{k \in \llbracket d \rrbracket} \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \int_{\mathbb{R}^d} \mathbf{1}_{\{|t_k| > \sqrt{2m_k+1}\}} \left| (\widehat{\Delta f^*}(\mathbf{t}))_1 \right| |\phi_\ell(\mathbf{t})| dt \right|^2 \right] + 2\mathrm{E} \left[ \int_{\mathbb{R}^d} \left| (\widehat{\Delta f^*}(\mathbf{t}))_2 \right|^2 dt \right] \\
&=: T_{2,1} + T_{2,2}
\end{align*}
$$

With Remark B.2 we get that

$$
T_{2,2} \le 2p3^{2(p-2)} \max_{j \in \llbracket p \rrbracket} \|f_j^*\|^2 \sum_{k \in \llbracket p \rrbracket} \frac{1}{n_k}.
$$

Consider $T_{2,1}$. Applying Cauchy-Schwarz, it holds

$$
\begin{align*}
T_{2,1} &\le 2d \sum_{k \in \llbracket d \rrbracket} \sum_{\ell=0}^{m-1} \mathrm{E} \left[ \left| \int_{\mathbb{R}^d} \mathbf{1}_{\{|t_k| > \sqrt{2m_k+1}\}} \left| (\widehat{\Delta f^*}(\mathbf{t}))_1 \right|^2 |\phi_\ell(\mathbf{t})| dt \int_{\mathbb{R}^d} \mathbf{1}_{\{|t_k| > \sqrt{2m_k+1}\}} |\phi_\ell(\mathbf{t})| dt \right] \\
&\le 2d \left( \prod_{j \in \llbracket p \rrbracket} \frac{1}{n_k} \right) \sum_{k \in \llbracket d \rrbracket} \sum_{\ell=0}^{m-1} \left( \prod_{\tilde{k} \in \llbracket d \rrbracket, \tilde{k} \neq k} \int |\varphi_{\tilde{k}}(t_{\tilde{k}})| dt_{\tilde{k}} \int \mathbf{1}_{\{|t_k| > \sqrt{2m_k+1}\}} |\varphi_{\ell_k}(t_k)| dt_k \right)^2
\end{align*}
$$

using also Equation (A.2). To consider the remaining integrals consider first the following remark.

**Remark** For $A \gg 0$ it holds that

$$
\int_A^\infty \exp(-cx^2) dx = \int_A^\infty x \exp(-cx^2) \frac{1}{x} dx \le \frac{1}{A} \int_A^\infty x \exp(-cx^2) dx = \frac{\exp(-cA)}{2cA}.
$$

Together with Equation (3.4), it follows for $k \in \llbracket d \rrbracket$ and $\ell_k \in \llbracket 0, m_k - 1 \rrbracket$ that

$$
\int_{|t_k| > \sqrt{2m_k+1}} |\varphi_{\ell_k}(t_k)| dt_k \le \int_{|t_k| > \sqrt{2m_k+1}} \kappa_{\varphi} \exp(-\xi t_k^2) dt_k \le \frac{\kappa_{\varphi} \exp(-\xi(2m_k + 1))}{\xi \sqrt{2m_k + 1}}.
$$

In addition, we have that

$$
\begin{align*}
\left( \int |\varphi_{\ell_k}(t_k)| dt_k \right)^2 &= \left( \int_{|t_k| < \sqrt{2m_k+1}} |\varphi_{\ell_k}(t_k)| dt_k + \int_{|t_k| > \sqrt{2m_k+1}} |\varphi_{\ell_k}(t_k)| dt_k \right)^2 \\
&\le 2 \left( \int \mathbf{1}_{\{|t_k| < \sqrt{2m_k+1}\}} dt_k \int |\varphi_{\ell_k}(t_k)|^2 dt_k \right) + \frac{2 \exp(-2\xi(2m_k+1))}{2m_k+1} \\
&= 2 \left( 2\sqrt{2m_k+1} + \frac{\exp(-2\xi(2m_k+1))}{2m_k+1} \right) \le 6\sqrt{2m_k+1}.
\end{align*}
$$

Consequently, using that *m* satisfying Equation (3.8) it follows

$$
\begin{align*}
T_{2,1} &\le \frac{2d6^{2(d-1)}\kappa_{\varphi}^2}{\xi^2} \left( \prod_{j \in \llbracket p \rrbracket} \frac{1}{n_j} \right) \sum_{k \in \llbracket d \rrbracket} \sum_{\ell=0}^{m-1} \prod_{\tilde{k} \in \llbracket d \rrbracket, \tilde{k} \neq k} \sqrt{2m_{\tilde{k}} + 1} \frac{\exp(-2\xi(2m_k + 1))}{2m_k + 1} \\
&\le \frac{2d6^{2(d-1)}\kappa_{\varphi}^2}{\xi^2} \left( \prod_{j \in \llbracket p \rrbracket} \frac{1}{n_j} \right) \sum_{k \in \llbracket d \rrbracket} \prod_{\tilde{k} \in \llbracket d \rrbracket} \sqrt{2m_{\tilde{k}} + 1} \exp(-4\xi m_k) \\
&\le \frac{2d6^{2(d-1)}\kappa_{\varphi}^2}{\xi^2} \left( \prod_{j \in \llbracket p \rrbracket} \frac{1}{n_j} \right) \sum_{k \in \llbracket d \rrbracket} \prod_{\tilde{k} \in \llbracket d \rrbracket} \sqrt{2m_{\tilde{k}} + 1} \exp(-2 \log \left( \prod_{j \in \llbracket p \rrbracket} n_j \right)) \\
&\le \frac{C_0 2d^2 6^{2(d-1)} \kappa_{\varphi}^2}{\xi^2} \left( \prod_{j \in \llbracket p \rrbracket} \frac{1}{n_j^2} \right),
\end{align*}
$$
where $C_0$ denotes the constant resulting from the assumption in Equation (3.8). Putting all bounds together it follows that

$$
\begin{align*}
\mathbb{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right] &\le \frac{2}{(2\pi)^d} (T_1 + T_{2,2} + T_{2,1}) \\
&\le \frac{2}{(2\pi)^d} \left( \frac{4 \prod_{k \in \llbracket d \rrbracket} \sqrt{2m_k + 1}}{\prod_{j \in \llbracket p \rrbracket} n_j} + 2p \cdot 3^{2(p-2)} \max_{j \in \llbracket p \rrbracket} \|f_j^*\|^2 \sum_{j \in \llbracket p \rrbracket} \frac{1}{n_j} \right. \\
&\qquad \left. + 2p 3^{2(p-2)} \max_{j \in \llbracket p \rrbracket} \|f_j^*\|^2 \sum_{j \in \llbracket p \rrbracket} \frac{1}{n_j} + \frac{C_0 2d^2 6^{2(d-1)} \kappa_\varphi^2}{\xi^2} \prod_{j \in \llbracket p \rrbracket} \frac{1}{n_j^2} \right) \\
&\le C_1 \frac{\prod_{k \in \llbracket d \rrbracket} \sqrt{2m_k + 1}}{\prod_{j \in \llbracket p \rrbracket} n_j} + C_2 \sum_{j \in \llbracket p \rrbracket} \frac{1}{n_j},
\end{align*}
$$

with $C_1 = \frac{8}{(2\pi)^d}$ and $C_2 = 4p3^{2(p-2)} \max_{j \in \llbracket p \rrbracket} \|f_j^*\|^2 + \frac{C_0 2d^2 6^{2(d-1)} \kappa_\varphi^2}{\xi^2}$, which proves the result. $\square$

**Proof of Theorem 3.4.** We consider a nesting space, i.e. a space $S_M$ such that, $\forall \mathbf{m} \in \mathcal{M}_n$, $S_m \subset S_M$, $\mathbf{M} = (M_1, \dots, M_d)$. As our worst case assumption for the spaces $S_m$ is $\sqrt{2m_j + 1} \le n^d$ in case $\mathbf{m} = (m_1, \dots, m_d) \in \mathcal{M}_n$ is such that all $m_k$ are equal to 1 except $m_j$, we set $M_j = \lfloor (n^{2d} - 1)/2 \rfloor$ for all $j \in \llbracket p \rrbracket$, where $\lfloor x \rfloor$ denotes the integer part of $x \in \mathbb{R}$.

**Step 1.** Using the nesting space $S_M$, we prove that, $\forall \mathbf{m} \in \mathcal{M}_n$,

$$
\begin{align*}
\mathbb{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right] &\le 3\mathbb{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right] + 6\mathbb{E} \left[ \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 \right] + 2\text{pen}(\mathbf{m}) \\
&\quad + 11\mathbb{E} \left[ \left( \| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 - \frac{2}{11} \text{pen}(\widehat{\mathbf{m}}) \right)_+ \right].
\end{align*}
$$

See the proof of Theorem 2 in Brenner Miguel (2022). Let us set

$$
\mathbb{E}_0 := \mathbb{E} \left[ \left( \| \widehat{g_m}^{(p)} - g_m^{(p)} \|^2 - \frac{2}{11} \mathrm{pen}(\widehat{m}) \right)_+ \right],
$$

which is the term we aim now to bound. We will show $\mathbb{E}_0 \lesssim n^{-1}$ which then yields the result.

**Step 2.** Decomposition of $\mathbb{E}_0$. It holds that

$$
\| \widehat{g_m^{(p)}} - g_m^{(p)} \|^2 = \sum_{\ell=0}^{m-1} (\widehat{a}_\ell - a_\ell)^2 = \sum_{\ell=0}^{m-1} \nu_n^2(\varphi_\ell) = \sup_{t \in S_m, \|t\|=1} \nu_n^2(t)
$$

where

$$
\nu_n(t) = \frac{1}{n^p} \sum_{i_k \in \llbracket n_k \rrbracket, k \in \llbracket d \rrbracket} \left( t(\mathbf{X}_{i_1}^{(1)} + \dots + \mathbf{X}_{i_p}^{(p)}) - \langle t, g^{(p)} \rangle \right).
$$

Therefore we need to bound

$$
\mathbb{E}_0 = \mathbb{E} \left[ \left( \sup_{t \in S_{\widehat{m}}, \|t\|=1} \nu_n^2(t) - \frac{2}{11} \mathrm{pen}(\widehat{m}) \right)_+ \right].
$$

We use the other writing of $\nu_n$, namely,

$$
\nu_n(t) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \left[ \prod_{j \in \llbracket p \rrbracket} \widehat{f}_j^*(\boldsymbol{x}) - \prod_{j \in \llbracket p \rrbracket} f_j^*(\boldsymbol{x}) \right] t^*(-\boldsymbol{x}) d\boldsymbol{x}.
$$
We define for some $\tau \in \mathbb{R}_+$, which will be chosen later:

$$ \Omega(\boldsymbol{x}) = \bigcap_{j=1}^{p} \Omega_j(\boldsymbol{x}), \quad \Omega_j(\boldsymbol{x}) = \left\{ |\widehat{f}_j^*(\boldsymbol{x}) - f_j^*(\boldsymbol{x})| \le \tau \sqrt{\frac{\log(n)}{n}} \right\}. \qquad (A.6) $$

Then we consider the decomposition $\nu_n(t) = \nu_{n,1}(t) + \nu_{n,2}(t) + \nu_{n,3}(t) + \nu_{n,4}(t)$ with

$$ \nu_{n,1}(t) = \frac{1}{(2\pi)^d} \int_{|\boldsymbol{x}|\le\sqrt{2\widehat{m}+1}} \prod_{j\in\llbracket p\rrbracket} [\widehat{f}_j^*(\boldsymbol{x}) - f_j^*(\boldsymbol{x})] t^*(-\boldsymbol{x}) \mathbf{1}_{\Omega(\boldsymbol{x})} d\boldsymbol{x}, $$

$$ \nu_{n,2}(t) = \frac{1}{(2\pi)^d} \int_{|\boldsymbol{x}|\le\sqrt{2\widehat{m}+1}} \prod_{j\in\llbracket p\rrbracket} [\widehat{f}_j^*(\boldsymbol{x}) - f_j^*(\boldsymbol{x})] t^*(-\boldsymbol{x}) \mathbf{1}_{\Omega(\boldsymbol{x})^c} d\boldsymbol{x}, $$

$$ \nu_{n,3}(t) = \frac{1}{(2\pi)^d} \int_{\exists j, |x_j| > \sqrt{2\widehat{m}+1}} \prod_{j\in\llbracket p\rrbracket} [\widehat{f}_j^*(\boldsymbol{x}) - f_j^*(\boldsymbol{x})] t^*(-\boldsymbol{x}) d\boldsymbol{x} $$

and

$$ \nu_{n,4}(t) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} (\widehat{\Delta f^*})_2(\boldsymbol{x}) t^*(-\boldsymbol{x}) d\boldsymbol{x}, \quad (\widehat{\Delta f^*})_2 = \left[ \prod_{j\in\llbracket p\rrbracket} \widehat{f}_j^* - \prod_{j\in\llbracket p\rrbracket} f_j^* \right] - \prod_{j\in\llbracket p\rrbracket} [\widehat{f}_j^* - f_j^*], $$

where we recall that $u \le v$ is to be understood component-wise if $u$ and $v$ are vectors. We get

$$ \mathbb{E}_0 \le 2\mathbb{E} \left[ \left( \sup_{t \in S_{\widehat{m}}, \|t\|=1} \nu_{n,1}^2(t) - \frac{1}{11} \mathrm{pen}(\widehat{m}) \right)_+ \right] + 6 \sum_{k=2}^4 \mathbb{E} \left( \sup_{t \in S_{\widehat{m}}, \|t\|=1} \nu_{n,k}^2(t) \right) \\ := 2\mathbb{E}_1 + 6(\mathbb{E}_2 + \mathbb{E}_3 + \mathbb{E}_4), $$

and we successively study the 4 terms.

**Step 3.** The bound on the main term $\mathbb{E}_1$ is straightforward. Indeed, using first the Cauchy-Schwarz inequality, and then taking the definition of $\Omega(\boldsymbol{x})$ into account, we get

$$ \begin{aligned} \sup_{t \in S_{\widehat{m}}, \|t\|=1} \nu_{n,1}^2(t) &\le \sup_{t \in S_{\widehat{m}}, \|t\|=1} \frac{1}{(2\pi)^{2d}} \int_{\mathbb{R}^d} |t^*(\boldsymbol{x})|^2 d\boldsymbol{x} \int_{|\boldsymbol{x}|\le\sqrt{2\widehat{m}+1}} \left| \prod_{j\in\llbracket p\rrbracket} (\widehat{f}_j^*(\boldsymbol{x}) - f_j^*(\boldsymbol{x})) \right|^2 \mathbf{1}_{\Omega(\boldsymbol{x})} d\boldsymbol{x} \\ &\le \frac{1}{(2\pi)^d} \int_{|\boldsymbol{x}|\le\sqrt{2\widehat{m}+1}} \left| \prod_{j\in\llbracket p\rrbracket} \tau \sqrt{\frac{\log(n)}{n}} \right|^2 d\boldsymbol{x} = \frac{\tau^{2p}}{(2\pi)^d} \prod_{k \in \llbracket d \rrbracket} 2\sqrt{2\widehat{m}_k + 1} \frac{\log^p(n)}{n^p} \\ &\le \frac{2^d \tau^{2p}}{(2\pi)^d} \prod_{k \in \llbracket d \rrbracket} \sqrt{2\widehat{m}_k + 1} \frac{\log^d(n)}{n^d} \end{aligned} $$

Now, if $\kappa$ is such that the above quantity is less than $(1/11)\mathrm{pen}(\widehat{m})$, that is

$$ \frac{2^d \tau^{2p}}{(2\pi)^d} \le \kappa/11, \quad \text{i.e.} \quad \kappa \ge \kappa_0 := 11 \frac{2^d \tau^{2p}}{(2\pi)^d} $$

we obtain that $\mathbb{E}_1 \le 0$ and therefore $\mathbb{E}_1 = 0$.
**Step 4.** The bound on $\mathbb{E}_2$ relies on Lemma 1 p.12 in Ammous et al. (2024), which states that for any $\tau > 0$, $\boldsymbol{x} \in \mathbb{R}^d$, $d \ge 1$,

$$ \mathbb{P}(\Omega_j(\boldsymbol{x})^c) \le 4n^{-\tau^2/4}. $$

Therefore, $\mathbb{P}(\Omega(\boldsymbol{x})^c) \le 4pn^{-\tau^2/4}$. With standard rough bounds we get

$$ \nu_{n,2}^2(t) \le \frac{1}{(2\pi)^d} \int_{|\boldsymbol{x}|\le\sqrt{2\widehat{m}+1}} 2^{2p} \mathbf{1}_{\Omega(\boldsymbol{x})^c} d\boldsymbol{x} $$

and

$$ \mathbb{E}_2 \le \frac{2^{2p}}{(2\pi)^d} \int_{|\boldsymbol{x}|\le\sqrt{2\widehat{M}+1}} \mathbb{P}(\Omega(\boldsymbol{x})^c) d\boldsymbol{x} \le 4 \frac{p2^{2p+2d}}{(2\pi)^d} \prod_{k\in\llbracket d\rrbracket} \sqrt{2M_k+1} n^{-\tau^2/4}. $$

Using that per definition $\sqrt{2\widehat{M}_k+1} \le n^d$ for all $k \in \llbracket d\rrbracket$, we have $\mathbb{E}_2 \lesssim n^{d^2-\tau^2/4}$. Thus $\mathbb{E}_2 \lesssim n^{-1}$ for

$$ \tau^2 \ge 4(d^2 - 1). $$

**Step 5.** The study of $\mathbb{E}_3$ follows the line of the study of $T_{2,1}$ in the proof of Proposition 3.2 (in Appendix A). For $t \in S_{\widehat{m}}$, $\|t\| = 1$ writing $t^* = \sum_{\ell=0}^{\widehat{m}-1} a_\ell \phi_\ell^*$, using that $\sum_{\ell=0}^{\widehat{m}-1} |a_\ell|^2 = 1$ and applying Cauchy-Schwarz twice, we get that

$$ \begin{aligned} & \mathbb{E} \left[ \sup_{t \in S_{\widehat{m}}, \|t\|=1} \nu_{n,3}^2(t) \right] \\ &= \frac{1}{(2\pi)^d} \mathbb{E} \left[ \sup_{t \in S_{\widehat{m}}, \|t\|=1} \left( \sum_{\ell=0}^{\widehat{m}-1} |a_\ell|^2 \right) \left( \sum_{\ell=0}^{\widehat{m}-1} \left| \int_{\exists i, |x_i| > \sqrt{2\widehat{m}_i+1}} \prod_{j \in \llbracket p \rrbracket} [\widehat{f}_j^*(\boldsymbol{x}) - f_j^*(\boldsymbol{x})] \phi_\ell(\boldsymbol{x}) d\boldsymbol{x} \right|^2 \right) \right] \\ & \le 2d \sum_{k \in \llbracket d \rrbracket} \mathbb{E} \left[ \sum_{\ell=0}^{\widehat{m}-1} \int_{\mathbb{R}^d} \mathbf{1}_{\{|x_k| > \sqrt{2\widehat{m}_k+1}\}} \left| (\Delta \widehat{f}^*)_1(\boldsymbol{x}) \right|^2 |\phi_\ell(\boldsymbol{x})| d\boldsymbol{x} \int_{\mathbb{R}^d} \mathbf{1}_{\{|x_k| > \sqrt{2\widehat{m}_k+1}\}} |\phi_\ell(\boldsymbol{x})| d\boldsymbol{x} \right]. \end{aligned} $$

Now we use that $\widehat{m}_k \ge c_0 \log(n)$ and $\prod_{k \in \llbracket d \rrbracket} \sqrt{2\widehat{m}_k+1} \le n^d$ as $\widehat{m} \in \mathcal{M}_n^*$ given by (3.10) and we get with $\kappa_\varphi$ defined in (3.4)

$$ \mathbb{E}_3 \le \frac{2d^2 6^{2(d-1)} \kappa_\varphi^2 \prod_{k \in \llbracket d \rrbracket} \sqrt{2M_k+1}}{\xi^2 n^p} \sum_{k \in \llbracket d \rrbracket} \exp(-2\xi c_0 \log(n)) = \frac{2d^2 6^{2(d-1)} d \kappa_\varphi^2}{\xi^2} \frac{1}{n}, $$

due to the definition of $c_0$. This ensures that $\mathbb{E}_3 \lesssim 1/n$.

**Step 6.** The study of $\mathbb{E}_4$ is also rather straightforward. Indeed, by Cauchy-Schwarz Inequality

$$ \mathbb{E}_4 \le \mathbb{E} \left( \sup_{t \in S_{\widehat{m}}, \|t\|=1} \frac{1}{(2\pi)^{2d}} \int_{\mathbb{R}^d} |t^*(\boldsymbol{x})|^2 d\boldsymbol{x} \int_{\mathbb{R}^d} (\widehat{\Delta f}^*)_2^2(\boldsymbol{x}) d\boldsymbol{x} \right) \le \frac{1}{(2\pi)^d} \mathbb{E} \left( \int_{\mathbb{R}^d} (\widehat{\Delta f}^*)_2^2(\boldsymbol{x}) d\boldsymbol{x} \right). $$

Now using Remark B.2 (revisited), we get

$$ \mathbb{E}_4 \le \frac{p3^{2(p-2)}}{(2\pi)^d} \sum_{j \in \llbracket p \rrbracket} \frac{1}{n} \int_{\mathbb{R}^d} |f_{j+1}^*(\boldsymbol{x})|^2 d\boldsymbol{x} \le p^2 3^{2(p-2)} \max_{j \in \llbracket p \rrbracket} (\|f_j\|^2) \frac{1}{n}, $$

as we assumed that all $n_j$ are equal to $n$. This concludes the proof. $\square$
B Auxiliary results

The following result is an extension of Lemma 1 in Chesneau et al. (2013) and is used in several proofs.

**Lemma B.1.** Let $(u, v) \in \mathbb{C}^p \times \mathbb{C}^p$ such that $|u_j|, |v_j| \le 1$ for all $j \in \llbracket p \rrbracket$. Then, for any integer $p \ge 1$, we have

$$ \left| \prod_{j \in \llbracket p \rrbracket} u_j - \prod_{j \in \llbracket p \rrbracket} v_j \right| \le \prod_{j \in \llbracket p \rrbracket} |u_j - v_j| + \sum_{j \in \llbracket p \rrbracket} C_j^p |u_j - v_j| |v_{j+1}| $$

setting $v_{p+1} := v_1$, where $C_1^1 = 0$ and for $p \ge 2$ we define $C_j^{p+1} := 3C_j^p$ for $j \in \llbracket p-1 \rrbracket$, $C_p^{p+1} := C_p^p + 2^{p-1}$ and $C_{p+1}^{p+1} := 1 + 2C_p^p$. In addition, for $p \ge 2$ and $j \in \llbracket p \rrbracket$ it holds that $C_j^p \le 3^{p-2}$.

**Remark B.2.** Similarly one shows that

$$ \left| \prod_{j \in \llbracket p \rrbracket} u_j - \prod_{j \in \llbracket p \rrbracket} v_j - \prod_{j \in \llbracket p \rrbracket} (u_j - v_j) \right| \le \sum_{j \in \llbracket p \rrbracket} C_j^p |u_j - v_j| |v_{j+1}|. $$

*Proof of Lemma B.1.* We proof both assertions simultaneously by induction.

**Base case** For $p=1$ this is obviously satisfied with $C_1^1 = 0$. For $p=2$ we see that for $(u, v) \in \mathbb{C}^2 \times \mathbb{C}^2$

$$ \begin{aligned} u_1 u_2 - v_1 v_2 &= (u_1 - v_1)(u_2 - v_2) + v_2 u_2 + v_2 u_1 - 2v_1 v_2 \\ &= (u_1 - v_1)(u_2 - v_2) + v_2(u_1 - v_1) + v_1(u_2 - v_2). \end{aligned} $$

Consequently, by the triangle inequality, the inequality holds for $C_1^2 = C_2^2 = 1$. Further $C_1^2 = C_1^2 \le 3^0 = 3^{2-2}$.

**Induction step** $p \mapsto p+1$ Assume that the inequality holds for $p \ge 2$. Let $(u, v) \in \mathbb{C}^{p+1} \times \mathbb{C}^{p+1}$ with $|u_j|, |v_j| \le 1$ for all $j \in \llbracket p+1 \rrbracket$. First, note that

$$ \begin{aligned} \prod_{j \in \llbracket p+1 \rrbracket} u_j - \prod_{j \in \llbracket p+1 \rrbracket} v_j &= u_{p+1} \left( \prod_{j \in \llbracket p \rrbracket} u_j - \prod_{j \in \llbracket p \rrbracket} v_j \right) + (u_{p+1} - v_{p+1}) \prod_{j \in \llbracket p \rrbracket} v_j \\ &= (u_{p+1} - v_{p+1}) \left( \prod_{j \in \llbracket p \rrbracket} u_j - \prod_{j \in \llbracket p \rrbracket} v_j \right) + v_{p+1} \left( \prod_{j \in \llbracket p \rrbracket} u_j - \prod_{j \in \llbracket p \rrbracket} v_j \right) + (u_{p+1} - v_{p+1}) \prod_{j \in \llbracket p \rrbracket} v_j. \end{aligned} $$

It follows using the inequality for $p$ on $(u_1, \dots, u_p, v_1, \dots, v_p)$ and that $|v_j|, |u_j| \le 1$ for all $j \in \llbracket p+1 \rrbracket$

$$ \begin{aligned} & \left| \prod_{j \in \llbracket p+1 \rrbracket} u_j - \prod_{j \in \llbracket p+1 \rrbracket} v_j \right| \\ & \le |u_{p+1} - v_{p+1}| \left| \prod_{j \in \llbracket p \rrbracket} u_j - \prod_{j \in \llbracket p \rrbracket} v_j \right| + |v_{p+1}| \left| \prod_{j \in \llbracket p \rrbracket} u_j - \prod_{j \in \llbracket p \rrbracket} v_j \right| + |u_{p+1} - v_{p+1}| \prod_{j \in \llbracket p \rrbracket} |v_j| \\ & \le |u_{p+1} - v_{p+1}| \left( \prod_{j \in \llbracket p \rrbracket} |u_j - v_j| + \sum_{j \in \llbracket p \rrbracket} C_j^p |u_j - v_j| |v_{j+1}| \right) \\ & \qquad + |v_{p+1}| \left( \prod_{j \in \llbracket p \rrbracket} |u_j - v_j| + \sum_{j \in \llbracket p \rrbracket} C_j^p |u_j - v_j| |v_{j+1}| \right) + |u_{p+1} - v_{p+1}| \prod_{j \in \llbracket p \rrbracket} |v_j|. \end{aligned} $$
Further,

$$
\begin{align*}
& \left| \prod_{j \in [p+1]} u_j - \prod_{j \in [p+1]} v_j \right| \\
& \le \prod_{j \in [p+1]} |u_j - v_j| + 2 \sum_{j \in [p-1]} C_j^p |u_j - v_j| |v_{j+1}| + 2C_p^p |u_{p+1} - v_{p+1}| |v_1| \\
& \quad + 2^{p-1} |v_{p+1}| |u_p - v_p| + \sum_{j \in [p-1]} C_j^p |u_j - v_j| |v_{j+1}| + C_p^p |v_{p+1}| |u_p - v_p| \\
& \quad + |u_{p+1} - v_{p+1}| |v_1| \\
& = \prod_{j \in [p+1]} |u_j - v_j| + \sum_{j \in [p-1]} 3C_j^p |u_j - v_j| |v_{j+1}| + (C_p^p + 2^{p-1}) |u_p - v_p| |v_{p+1}| \\
& \quad + (1 + 2C_p^p) |u_{p+1} - v_{p+1}| |v_1| = \prod_{j \in [p+1]} |u_j - v_j| + \sum_{j \in [p+1]} C_j^{p+1} |u_j - v_j| |v_{j+1}|,
\end{align*}
$$

where $C_j^{p+1} = 3C_j^p$ for $j \in [p-1]$, $C_p^{p+1} = C_p^p + 2^{p-1}$ and $C_{p+1}^{p+1} = 1 + 2C_p^p$.

For the second assertion, assume that $C_j^p \le 3^{p-2}$ for $p$. Then it follows that

$$
C_j^{p+1} = 3C_j^p \le 3 \cdot 3^{p-2} = 3^{(p+1)-2}, \quad \text{for } j \in [(p+1)-1],
$$

$$
C_p^{p+1} = C_p^p + 2^{p-1} \le 3^{p-2} + 2^{p-1} = 3^{(p+1)-2} \left( \frac{1}{3} + \left(\frac{2}{3}\right)^{p-1} \right) \le 3^{(p+1)-2}
$$

$$
C_{p+1}^{p+1} = 1 + 2C_p^p \le 1 + 2 \cdot 3^{p-2} \le 3^{(p+1)-2} \left( 3^{-p+1} + \frac{2}{3} \right) \le 3^{(p+1)-2}.
$$

This concludes the proof.

<figure><img src="image_9.png" /></figure>

References

Ahmad, I. and Fan, Y. (2001). Optimal bandwidths for kernel density estimators of functions of obser-
vations. *Statistics & Probability Letters*, 51(3):245–251.

Ahmad, I. A. and Mugdadi, A. R. (2003). Analysis of kernel density estimation of functions of random variables. *Journal of Nonparametric Statistics*, 15(4–5):579–605.

Ammous, S., Dedecker, J., and Duval, C. (2024). Adaptive directional estimator of the density in rd for independent and mixing sequences. *Journal of Multivariate Analysis*, 203:105332.

Askey, R. and Wainger, S. (1965). Mean convergence of expansions in laguerre and hermite series. *American Journal of Mathematics*, 87(3):695–708.

Bandyopadhyay, S. (2012). A note on efficient density estimators of convolutions. *Journal of Statistical Planning and Inference*, 142(11):3056–3060.

Baudry, J.-P., Maugis, C., and Michel, B. (2012). Slope heuristics: overview and implementation. *Stat. Comput.*, 22(2):455–470.
Belomestny, D., Comte, F., and Genon-Catalot, V. (2019). Sobolev-Hermite versus Sobolev nonparametric density estimation on $\mathbb{R}$. *Ann. Inst. Statist. Math.*, 71(1):29-62.

Bongioanni, B. and Torrea, J. L. (2006). Sobolev spaces associated to the harmonic oscillator. *Proc. Indian Acad. Sci. Math. Sci.*, 116(3):337-360.

Brenner Miguel, S. (2022). Anisotropic spectral cut-off estimation under multiplicative measurement errors. *Journal of Multivariate Analysis*, 190:104990.

Cherubini, U., Mulinacci, S., and Romagnoli, S. (2011). On the distribution of the (un)bounded sum of random variables. *Insurance: Mathematics and Economics*, 48(1):56-63.

Chesneau, C., Comte, F., and Navarro, F. (2013). Fast nonparametric estimation for convolutions of densities. *Canadian Journal of Statistics*, 41(4):617-636.

Cohen, A., Davenport, M. A., and Leviatan, D. (2013). On the stability and accuracy of least squares approximations. *Found. Comput. Math.*, 13(5):819-834.

Comte, F. (2017). *Nonparametric Estimation*. Spartacus Supérieur. Spartacus-Idh.

Comte, F., Duval, C., and Genon-Catalot, V. (2014). Nonparametric density estimation in compound Poisson processes using convolution power estimators. *Metrika*, 77(1):163-183.

Comte, F. and Lacour, C. (2013). Anisotropic adaptive kernel deconvolution. *Annales de l'IHP Probabilités et statistiques*, 49:569-609.

Comte, F. and Lacour, C. (2023). Non compact estimation of the conditional density from direct or noisy data. *Ann. Inst. Henri Poincaré, Probab. Stat.*, 59(3):1463-1507.

Comte, F. and Marie, N. (2020). Bandwidth selection for the wolverton-wagner estimator. *Journal of Statistical Planning and Inference*, 207:198-214.

Duval, C. (2013). Density estimation for compound poisson processes from discrete data. *Stochastic Processes and their Applications*, 123(11):3963-3986.

Frees, E. W. (1994). Estimating densities of functions of observations. *Journal of the American Statistical Association*, 89(426):517-525.

Giné, E. and Mason, D. M. (2007a). Laws of the iterated logarithm for the local u-statistic process. *Journal of Theoretical Probability*, 20(3):457-485.

Giné, E. and Mason, D. M. (2007b). On local u-statistic processes and the estimation of densities of functions of several sample variables. *The Annals of Statistics*, 35(3).

Guo, H., Wang, J., and Tian, X. (2018). The mean consistency of wavelet estimators for convolutions of the density functions. *Journal of Computational and Applied Mathematics*, 343:1-11.

Hu, C., Pozdnyakov, V., and Yan, J. (2019). Density and distribution evaluation for convolution of independent gamma variables. *Computational Statistics*, 35(1):327-342.

Lacour, C., Massart, P., and Rivoirard, V. (2017). Estimator selection: a new method with applications to kernel density estimation. *Sankhya A*, 79:298-335.
Midhu, N. N., Dewan, I., Sudheesh, K. K., and Sreedevi, E. P. (2023). On approximation and estimation of distribution function of sum of independent random variables. *Statistical Papers*, 65(3):1437–1467.

Nickl, R. (2007). Donsker-type theorems for nonparametric maximum likelihood estimators. *Probability Theory and Related Fields*, 138(3–4).

Nickl, R. (2009). On convergence and convolutions of random signed measures. *Journal of Theoretical Probability*, 22(1):38–56.

Panjer, H. H. and Willmot, G. E. (1992). *Insurance risk models*. Society of Acturaries.

Prakasa Rao, B. L. S. (2004). Moment inequalities for supremum of empirical processes of U-statistic structure and application to density estimation. *J. Iran. Stat. Soc. JIRSS*, 3(1):59–68.

Saavedra, A. and Cao, R. (2000). On the estimation of the marginal density of a moving average process. *Canadian Journal of Statistics*, 28(4):799–815.

Sacko, O. B. (2020). Hermite density deconvolution. *ALEA : Latin American Journal of Probability and Mathematical Statistics*, 17:419–443.

Schick, A. and Wefelmeyer, W. (2004a). Root n consistent and optimal density estimators for moving average processes. *Scandinavian Journal of Statistics*, 31(1):63–78.

Schick, A. and Wefelmeyer, W. (2004b). Rootnconsistent density estimators for sums of independent random variables. *Journal of Nonparametric Statistics*, 16(6):925–935.

Schick, A. and Wefelmeyer, W. (2007). Root-n consistent density estimators of convolutions in weighted 11-norms. *Journal of Statistical Planning and Inference*, 137(6):1765–1774.

Schick, A. and Wefelmeyer, W. (2008). Root-n consistency in weighted 1 1-spaces for density estimators of invertible linear processes. *Statistical Inference for Stochastic Processes*, 11(3):281–310.

Schick, A. and Wefelmeyer, W. (2009a). Convergence rates of density estimators for sums of powers of observations. *Metrika*, 69(2–3):249–264.

Schick, A. and Wefelmeyer, W. (2009b). Non-standard behavior of density estimators for sums of squared observations. *Statistics & Decisions*, 27(1):55–73.

Schick, A. and Wefelmeyer, W. (2012). On efficient estimation of densities for sums of squared observations. *Statistics & Probability Letters*, 82(9):1637–1640.

Taupin, M.-L. (2001). Semi-parametric estimation in the nonlinear structural errors-in-variables model. *The Annals of Statistics*, 29(1):66–93.

Willmot, G. E. and Woo, J.-K. (2007). On the class of erlang mixtures with risk theoretic applications. *North American Actuarial Journal*, 11(2):99–115.