Nonparametric estimation in a regression model with additive and multiplicative noise

TL;DR

This paper introduces two wavelet estimators for nonparametric regression models with both additive and multiplicative noise, achieving fast convergence rates under weak conditions, supported by numerical experiments.

Contribution

The paper proposes novel wavelet estimators for complex noise models in nonparametric regression, with theoretical convergence guarantees and practical validation.

Findings

Estimators achieve fast convergence rates under weak conditions.

Numerical results support theoretical claims.

Applicable to models similar to stochastic frontier estimation.

Abstract

In this paper, we consider an unknown functional estimation problem in a general nonparametric regression model with the feature of having both multiplicative and additive noise.We propose two new wavelet estimators in this general context. We prove that they achieve fast convergence rates under the mean integrated square error over Besov spaces. The obtained rates have the particularity of being established under weak conditions on the model. A numerical study in a context comparable to stochastic frontier estimation (with the difference that the boundary is not necessarily a production function) supports the theory.

Equations

$$Y_i=f({\mathbf{X}}_i)U_i+V_i,i\in \{1,\dots ,n\},$$

$$Y_i=g(X_i)+f(X_i)U_i$$

$$h(\mathbf{x})=\sum _{\mathbf{k}\in {\Lambda }_{j_{\ast }}}{\alpha }_{j_{\ast },\mathbf{k}}{\Phi }_{j_{\ast },\mathbf{k}}(\mathbf{x})+\sum _{j=j_{\ast }}^{\infty }\sum _{u=1}^{2^d-1}\sum _{\mathbf{k}\in {\Lambda }_j}{\beta }_{j,\mathbf{k},u}{\Psi }_{j,\mathbf{k},u}(\mathbf{x}),\mathbf{x}\in [0,1{]}^d,$$

$$P_{j}h(\mathbf{x})=\sum _{\mathbf{k}\in {\Lambda }_j}{\alpha }_{j,\mathbf{k}}{\Phi }_{j,\mathbf{k}}(\mathbf{x}),\mathbf{x}\in [0,1{]}^d.$$

$$\Vert f{\Vert }_{B_{p,q}^s}:=\Vert ({\alpha }_{\tau ,.}){\Vert }_p+\Vert (2^{j(s-\frac{d}{p}+\frac{d}{2})}\Vert {\beta }_{j,.,.}{\Vert }_p{)}_{j\ge \tau }{\Vert }_q,\text{where}\Vert {\beta }_{j,.,.}{\Vert }_p^p=\sum _{u=1}^{2^d-1}\sum _{\mathbf{k}\in {\Lambda }_j}|{\beta }_{j,\mathbf{k},u}{|}^p.$$

$${\hat{\alpha }}_{j,\mathbf{k}}:=\frac{1}{n}\sum _{i=1}^n{Y}_i^2{\Phi }_{j,\mathbf{k}}({\mathbf{X}}_i)-v_{j,\mathbf{k}},$$

$$v_{j,\mathbf{k}}:=\{\begin{array}{rlrl}& \mathbf{E}\left[V_1^2\right]2^{-jd/2}& & \text{under refhyp:A5},\\ & {\int }_{[0,1{]}^d}{g}^2(\mathbf{x}){\Phi }_{j,\mathbf{k}}(\mathbf{x})d\mathbf{x}& & \text{under refhyp:A6}.\end{array}$$

$${\hat{\beta }}_{j,\mathbf{k},u}:=\frac{1}{n}\sum _{i=1}^n\left(Y_i^2{\Psi }_{j,\mathbf{k},u}({\mathbf{X}}_i)-w_{j,\mathbf{k},u}\right)\text{mathbbm}{1}_{\left\{|Y_i^2{\Psi }_{j,\mathbf{k},u}({\mathbf{X}}_i)-w_{j,\mathbf{k},u}|\le {\rho }_n\right\}},$$

$$w_{j,\mathbf{k},u}:=\{\begin{array}{rlrl}& 0& & \text{under refhyp:A5},\\ & {\int }_{[0,1{]}^d}{g}^2(\mathbf{x}){\Psi }_{j,\mathbf{k},u}(\mathbf{x})d\mathbf{x}& & \text{under refhyp:A6}.\end{array}$$

$${\hat{r}}_n^{\mathrm{lin}}(\mathbf{x}):=\sum _{\mathbf{k}\in {\Lambda }_{j_{\ast }}}{\hat{\alpha }}_{j_{\ast },\mathbf{k}}{\Phi }_{j_{\ast },\mathbf{k}}(\mathbf{x}),\mathbf{x}\in [0,1{]}^d.$$

$${\hat{r}}_n^{\mathrm{non}}(\mathbf{x}):=\sum _{\mathbf{k}\in {\Lambda }_{j_{\ast }}}{\hat{\alpha }}_{j_{\ast },\mathbf{k}}{\Phi }_{j_{\ast },\mathbf{k}}(\mathbf{x})+\sum _{j=j_{\ast }}^{j_1}\sum _{u=1}^{2^d-1}\sum _{\mathbf{k}\in {\Lambda }_j}{\hat{\beta }}_{j,\mathbf{k},u}\text{mathbbm}{1}_{\{|{\hat{\beta }}_{j,\mathbf{k},u}|\ge \kappa t_n\}}{\Psi }_{j,\mathbf{k},u}(\mathbf{x}),\mathbf{x}\in [0,1{]}^d,$$

$$\mathbf{E}\left[{\int }_{[0,1{]}^d}{\left({\hat{r}}_n^{\mathrm{lin}}(\mathbf{x})-r(\mathbf{x})\right)}^{2}d\mathbf{x}\right]\lesssim n^{-\frac{2s^{\prime }}{2s^{\prime }+d}};$$

$$\mathbf{E}\left[{\int }_{[0,1{]}^d}{\left({\hat{r}}_n^{\mathrm{non}}(\mathbf{x})-r(\mathbf{x})\right)}^{2}d\mathbf{x}\right]\lesssim (\ln n)n^{-\frac{2s}{2s+d}}.$$

$${\hat{\alpha }}_{j,\mathbf{k}}=\frac{1}{n}\sum _{i=1}^n\frac{Y_i^2}{h({\mathbf{X}}_i)}{\Phi }_{j,\mathbf{k}}({\mathbf{X}}_i)-v_{j,\mathbf{k}},$$

$${\hat{\beta }}_{j,\mathbf{k},u}=\frac{1}{n}\sum _{i=1}^n\left(\frac{Y_i^2}{h({\mathbf{X}}_i)}{\Psi }_{j,\mathbf{k},u}({\mathbf{X}}_i)-w_{j,\mathbf{k},u}\right)\text{mathbbm}{1}_{\left\{\left|\frac{Y_i^2}{h({\mathbf{X}}_i)}{\Psi }_{j,\mathbf{k},u}({\mathbf{X}}_i)-w_{j,\mathbf{k},u}\right|\le {\rho }_n\right\}}.$$

$$\mathbf{E}[{\hat{\alpha }}_{j,\mathbf{k}}]={\alpha }_{j,\mathbf{k}},\mathbf{E}\left[\frac{1}{n}\sum _{i=1}^n{Y}_i^2{\Psi }_{j,\mathbf{k},u}({\mathbf{X}}_i)-w_{j,\mathbf{k},u}\right]={\beta }_{j,\mathbf{k},u}.$$

$$\mathbf{E}\left[U_1{V}_{1}f({\mathbf{X}}_1){\Phi }_{j,\mathbf{k}}({\mathbf{X}}_1)\right]$$

$$=0,$$

$$v_{j,\mathbf{k}}$$

$$=\mathbf{E}\left[V_1^2{\Phi }_{j,\mathbf{k}}({\mathbf{X}}_1)\right].$$

$$\mathbf{E}[{\hat{\alpha }}_{j,\mathbf{k}}]$$

$$=\mathbf{E}\left[U_1^{2}r({\mathbf{X}}_1){\Phi }_{j,\mathbf{k}}({\mathbf{X}}_1)\right]+2\mathbf{E}\left[U_1{V}_{1}f({\mathbf{X}}_1){\Phi }_{j,\mathbf{k}}({\mathbf{X}}_1)\right]+\mathbf{E}\left[V_1^2{\Phi }_{j,\mathbf{k}}({\mathbf{X}}_1)\right]-v_{j,\mathbf{k}}$$

$$=\mathbf{E}\left[U_1^2\right]\mathbf{E}\left[r({\mathbf{X}}_1){\Phi }_{j,\mathbf{k}}({\mathbf{X}}_1)\right]={\int }_{[0,1{]}^d}r(\mathbf{x}){\Phi }_{j,\mathbf{k}}(\mathbf{x})d\mathbf{x}={\alpha }_{j,\mathbf{k}}.$$

$$w_{j,\mathbf{k},u}$$

$$=\mathbf{E}\left[V_1^2{\Psi }_{j,\mathbf{k},u}({\mathbf{X}}_1)\right].$$

$$\mathbf{E}\left[({\hat{\alpha }}_{j,\mathbf{k}}-{\alpha }_{j,\mathbf{k}}{)}^2\right]\lesssim \frac{1}{n},\mathbf{E}\left[({\hat{\beta }}_{j,\mathbf{k},u}-{\beta }_{j,\mathbf{k},u}{)}^2\right]\lesssim \frac{\ln n}{n}.$$

$$\mathbf{E}\left[({\hat{\alpha }}_{j,\mathbf{k}}-{\alpha }_{j,\mathbf{k}}{)}^2\right]$$

$$=\mathbf{V}\left[\frac{1}{n}\sum _{i=1}^n{Y}_i^2{\Phi }_{j,\mathbf{k}}({\mathbf{X}}_i)\right]=\frac{1}{n}\mathbf{V}\left[Y_1^2{\Phi }_{j,\mathbf{k}}({\mathbf{X}}_1)\right]\le \frac{1}{n}\mathbf{E}\left[Y_1^4{\Phi }_{j,\mathbf{k}}^2({\mathbf{X}}_1)\right]$$

$$\lesssim \frac{1}{n}\left(\mathbf{E}\left[U_1^4{f}^4({\mathbf{X}}_1){\Phi }_{j,\mathbf{k}}^2({\mathbf{X}}_1)\right]+\mathbf{E}\left[V_1^4{\Phi }_{j,\mathbf{k}}^2({\mathbf{X}}_1)\right]\right)$$

$$=\frac{1}{n}\left(\mathbf{E}\left[U_1^4\right]\mathbf{E}\left[f^4({\mathbf{X}}_1){\Phi }_{j,\mathbf{k}}^2({\mathbf{X}}_1)\right]+\mathbf{E}\left[V_1^4{\Phi }_{j,\mathbf{k}}^2({\mathbf{X}}_1)\right]\right).$$

$$\mathbf{E}\left[V_1^4{\Phi }_{j,\mathbf{k}}^2({\mathbf{X}}_1)\right]=\{\begin{array}{ll}\mathbf{E}\left[V_1^4\right]\mathbf{E}\left[{\Phi }_{j,\mathbf{k}}^2({\mathbf{X}}_1)\right]\lesssim 1& \text{under refhyp:A5},\\ {\int }_{[0,1{]}^d}{g}^4(\mathbf{x}){\Phi }_{j,\mathbf{k}}^2(\mathbf{x})d\mathbf{x}\lesssim {\int }_{[0,1{]}^d}{\Phi }_{j,\mathbf{k}}^2(\mathbf{x})d\mathbf{x}=1& \text{under refhyp:A6}.\end{array}$$

$${\hat{\beta }}_{j,\mathbf{k},u}-{\beta }_{j,\mathbf{k},u}=\frac{1}{n}\sum _{i=1}^n{D}_i-\mathbf{E}\left[K_1\text{mathbbm}{1}_{\{|K_1|>{\rho }_n\}}\right].$$

$$\mathbf{E}\left[({\hat{\beta }}_{j,\mathbf{k},u}-{\beta }_{j,\mathbf{k},u}{)}^2\right]$$

$$\lesssim \frac{1}{n}\mathbf{E}\left[D_1^2\right]+{\left(\mathbf{E}\left[|K_1|\text{mathbbm}{1}_{\{|K_1|>{\rho }_n\}}\right]\right)}^2.$$

$${\left(\mathbf{E}\left[K_1\text{mathbbm}{1}_{\{|K_1|>{\rho }_n\}}\right]\right)}^2\lesssim {\left(\mathbf{E}\left[K_1^2\right]/{\rho }_n\right)}^2\lesssim \frac{\ln n}{n}.$$

Full text

Keywords: Nonparametric regression, multiplicative regression models, nonparametric frontier, rates of convergence, wavelets.

1 Introduction

We consider a nonparametric regression model with both multiplicative and additive noise. It is defined by $n$ random variables $Y_{1},\ldots,Y_{n}$, where

[TABLE]

$f$ is an unknown regression function defined on a subset $\Delta$ of $\mathbb{R}^{d}$, with $d\geq 1$, $\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{n}$ are $n$ identically distributed random vectors with support on $\Delta$, $U_{1},\ldots,U_{n}$ are $n$ identically distributed random variables and $V_{1},\ldots,V_{n}$ are $n$ identically distributed random variables. Moreover, it is supposed that $\boldsymbol{X}_{i}$ and $U_{i}$ are independent, and $U_{i}$ and $V_{i}$ are independent for any $i\in\{1,\ldots,n\}$. We are interested in the estimation of the unknown function $r:=f^{2}$ from $(\boldsymbol{X}_{1},Y_{1}),\ldots,(\boldsymbol{X}_{n},Y_{n})$; the random vectors $(U_{1},V_{1}),\ldots,(U_{n},V_{n})$ form the multiplicative-additive noise. We consider the general formulation of model given by (1) since besides the theoretical interest it embodies several potential applications. For example, for $U_{i}=1$, (1) becomes the standard nonparametric regression model with additive noise. It has been studied in many papers via various nonparametric methods, including kernel, splines, projection and wavelets methods. See, for instance, the books of Härdle et al. (2012), Tsybakov (2009) and Comte (2015), and the references therein. For $V_{i}=0$, (1) becomes the standard nonparametric regression model with multiplicative noise. Recent studies can be found in Chichignoud (2012); Comte (2015) and the references therein. For $V_{i}\neq 0$ with the same variance across $i$ a first study, based on a linear wavelet estimator, was proposed by Chesneau et al. (2019). In the case where $V_{i}$ is a function of $\boldsymbol{X}_{i}$, (1) becomes the popular nonparametric regression model with multiplicative heteroscedastic noises:

[TABLE]

In particular, this model is widely used in financial applications, where the aim is to estimate the variance function $r:=f^{2}$ from the returns of an asset, for instance, to establish confidence intervals/bands for the mean function $g$. Variance estimation is a fundamental statistical problem with wide applications (see Muller et al. (1987); Hall and Marron (1990); Härdle and Tsybakov (1997); Wang et al. (2008); Brown et al. (2007); Cai et al. (2008) for fixed design and recently Kulik et al. (2011); Verzelen et al. (2018); Shen et al. (2019) for random design).

This multiplicative regression model is also popular in various application areas. For example, in econometrics, within deterministic ($V_{i}=0$) and stochastic ($V_{i}\neq 0$) non-parametric frontier models. These models can be interpreted as a special case of the model (1), where the random variable $U_{i}$ represents the technical inefficiency of the company and $V_{i}$ represents noise that disrupts its performance, the nature of which comes from unanticipated events such that machine failure, strikes, staff strikes, etc. Under monotonicity and concavity assumptions, the regression function $r$ can be viewed in this case as a function of the production set of a firm and its estimation is therefore of paramount importance in production econometrics. Specific estimation methods have been developed, see for instance Farrell (1957); De Prins et al. (1984); Gijbels et al. (1999); Daouia and Simar (2005) for deterministic frontiers models and Fan et al. (1996); Kumbhakar et al. (2007); Simar and Zelenyuk (2011) for stochastic frontier models. For general regression function and general nonparametric setting, we refer to Girard and Jacob (2008); Girard et al. (2013) and Jirak et al. (2014) for the definitions and properties of robust estimators.

Applications also exist in signal and image processing (e.g., for Global Positioning System signal detection Huang et al. (2013) as well as in speckle noise reduction encounter in particular in synthetic-aperture radar images Kuan et al. (1985) or in medical ultrasound images Rabbani et al. (2008); Mateo and Fernández-Caballero (2009)), where noise sources can be both additive and multiplicative. In this context, one can also cite Korostelev and Tsybakov (2012) where the author deals with the estimation of the function’s support.

The aim of this paper is to develop wavelet methods for the general model (1), with a special focus on mild assumptions on the distributions of $U_{1}$ and $V_{1}$ (moments of order $4$ will be required, including Gaussian noise). Wavelet methods are of interest in nonparametric statistics thanks to their ability to estimate efficiently a wide variety of unknown functions, including those with spikes and bumps. We refer to Abramovich et al. (2000) and Härdle et al. (2012), and the references therein. To the best of our knowledge, their development for (1) taking in full generality is new in the literature. First of all, we construct a linear wavelet estimator using projections of wavelet coefficients estimators. We evaluate its rate of convergence under the mean integrated square error (MISE) under mild assumptions on the smoothness of $r$; it is assumed that $r$ belongs to Besov spaces. The linear wavelet estimator has the advantage to be simple, but the knowledge of the smoothness of $r$ is necessary to calibrate a tuning parameter which plays a crucial role in the determination of fast rates of convergence. For this reason, an alternative is given by a nonlinear wavelet estimator. Using a thresholding rule of wavelet coefficients estimators, we develop a nonlinear wavelet estimator. To reach the goal of mildness assumptions on the model, we use a truncation rule in the definition of the wavelet coefficients estimators. This technique was introduced by Delyon and Juditsky (1996) in the nonparametric regression estimation setting, and recently improved in Chesneau (2013) (in a multidimensional regression function under mixing dependence framework) and Chaubey et al. (2015) (for a density estimation under a multiplicative censoring problem). The construction of the hard wavelet estimator does not depend on the smoothness of $r$ and we prove that, from a global point of view, it achieves a better rate of convergence under the MISE. In practice, the empirical performance of the estimators developed in this paper depends on the choice of several parameters, the truncation level of the linear estimator as well as the threshold parameter of the non-linear estimator. We propose here a method of automatic selection of these two parameters based on the 2-fold cross-validation method (2FCV) introduced by Nason (1996). A numerical study, in a context similar to stochastic frontier estimation, is being carried out to demonstrate the applicability of this approach.

The rest of the paper is organized as follows. Preliminaries on wavelets are described in Section 2. Section 3 specifies some assumptions on the model, presents our wavelet estimators and the main results on their performances. Numerical experiments are presented in Section 4. Section 5 is devoted to the proofs of the main result.

2 Preliminaries on wavelets

2.1 Wavelet bases on $[0,1]^{d}$

We begin with a classical notation in wavelet analysis. A multiresolution analysis (MRA) is a sequence of closed subspaces $\{V_{j}\}_{j\in\mathbb{Z}}$ of the square integrable function space $L^{2}(\mathbb{R})$ satisfying the following properties:

1. (i)

$V_{j}\subseteq V_{j+1}$, $j\in\mathbb{Z}$. $\mathbb{Z}$ denotes the integer set and $\mathbb{N}:=\{n\in\mathbb{Z},n\geq 0\};$ 2. (ii)

$\overline{\bigcup\limits_{j\in\mathbb{Z}}V_{j}}=L^{2}(\mathbb{R})$ (the space $\bigcup\limits_{j\in\mathbb{Z}}V_{j}$ is dense in $L^{2}(\mathbb{R})$); 3. (iii)

$f(2\cdot)\in V_{j+1}$ if and only if $f(\cdot)\in V_{j}$ for each $j\in\mathbb{Z}$; 4. (iv)

There exists $\phi\in L^{2}(\mathbb{R})$ (scaling function) such that $\{\phi(\cdot-k),\ k\in\mathbb{Z}\}$ forms an orthonormal basis of $V_{0}=\overline{\operatorname*{span}}\{\phi(\cdot-k)\}$.

See Meyer (1992) for further details. For the purpose of this paper, we use the compactly supported scaling function $\phi$ of the Daubechies family, and the associated compactly supported wavelet function $\psi$ (see Daubechies (1992)). Then we consider the wavelet tensor product bases on $[0,1]^{d}$ as described in Cohen et al. (1993). The main lines and notations are described below. We set $\Phi(\boldsymbol{x})=\prod_{v=1}^{d}\phi(x_{v})$ and wavelet functions: $\Psi_{u}(\boldsymbol{x})=\psi(x_{u})\prod_{\underset{v\not=u}{v=1}}^{d}\phi(x_{v})$ when $u\in\{1,\ldots,d\}$, and $\Psi_{u}(\boldsymbol{x})=\prod_{v\in A_{u}}\psi(x_{v})\prod_{v\not\in A_{u}}\phi(x_{v})$ when $u\in\{d+1,\ldots,2^{d}-1\}$, where $(A_{u})_{u\in\{d+1,\ldots,2^{d}-1\}}$ forms the set of all the non-void subsets of $\{1,\ldots,d\}$ of cardinal superior or equal to $2$. For any integer $j$ and any $\boldsymbol{k}=(k_{1},\ldots,k_{d})$, we set $\Phi_{j,\boldsymbol{k}}(\boldsymbol{x})=2^{jd/2}\Phi(2^{j}x_{1}-k_{1},\ldots,2^{j}x_{d}-k_{d})$, for any $u\in\{1,\ldots,2^{d}-1\}$, $\Psi_{j,\boldsymbol{k},u}(\boldsymbol{x})=2^{jd/2}\Psi_{u}(2^{j}x_{1}-k_{1},\ldots,2^{j}x_{d}-k_{d})$. Now, let us set $\Lambda_{j}=\{0,\ldots,2^{j}-1\}^{d}$. Then, with an appropriate treatment on the elements which step on the boundaries [math] and $1$, there exists an integer $\tau$ such that the system $\mathcal{S}=\{\Phi_{\tau,\boldsymbol{k}},\boldsymbol{k}\in\Lambda_{\tau};\ (\Psi_{j,\boldsymbol{k},u})_{u\in\{1,\ldots,2^{d}-1\}},\ j\geq\tau,\ \boldsymbol{k}\in\Lambda_{j}\}$ forms an orthonormal basis of $L^{2}([0,1]^{d})$. For any integer $j_{*}\geq\tau$, a function $h\in L^{2}([0,1]^{d})$ can be expressed via $\mathcal{S}$ by the following wavelet series:

[TABLE]

where $\alpha_{j,\boldsymbol{k}}=\langle h,\Phi_{j,\boldsymbol{k}}\rangle_{[0,1]^{d}}$ and $\beta_{j,\boldsymbol{k},u}=\langle h,\Psi_{j,\boldsymbol{k},u}\rangle_{[0,1]^{d}}$.

Also, let us mention that, by construction, $\int_{[0,1]^{d}}\Phi_{j,\boldsymbol{k}}(\boldsymbol{x})d\boldsymbol{x}=2^{-jd/2}$ and $\int_{[0,1]^{d}}\Psi_{j,\boldsymbol{k},u}(\boldsymbol{x})d\boldsymbol{x}=0$.

Let $P_{j}$ be the orthogonal projection operator from $L^{2}([0,1]^{d})$ onto the space $V_{j}$ with the orthonormal basis $\{\Phi_{j,\boldsymbol{k}}(\cdot)=2^{jd/2}\Phi(2^{j}\cdot-\boldsymbol{k}),\boldsymbol{k}\in\Lambda_{j}\}$. Then, for any $h\in L^{2}([0,1]^{d})$,

[TABLE]

2.2 Besov spaces

Besov spaces are important in theory and applications. They have the features to a wide variety of function spaces as Hölder and $L^{2}$ Sobolev spaces. Definitions of those spaces are given below. Suppose that $\phi$ is $m$ regular (i.e., $\phi\in C^{m}$ and $|D^{\alpha}\phi(x)|\leq c(1+|x|^{2})^{-l}$ for each $l\in\mathbb{Z}$, with $\alpha=0,1,\ldots,m$) and consider the wavelet framework defined in Subsection 2.1. Let $h\in L^{p}([0,1]^{d})$, $p,q\in[1,\infty]$ and $0<s<m$. Then the following assertions are equivalent:

(1) $h\in B^{s}_{p,q}([0,1]^{d})$;  (2) $\left\{2^{js}\|P_{j+1}h-P_{j}h\|_{p}\right\}\in l_{q};$  (3) $\{2^{j(s-\frac{d}{p}+\frac{d}{2})}\|\beta_{j,.,.}\|_{p}\}\in l_{q}.$

The Besov norm of $h$ can be defined by

[TABLE]

Further details on Besov spaces are given in Meyer (1992), Triebel (1994) and Härdle et al. (2012).

3 Assumptions, estimators and main result

We consider the model (1) with $\Delta=[0,1]^{d}$ for the sake of simplicity. Additional technical assumptions are formulated below.

1. A.1

We suppose that $f:[0,1]^{d}\rightarrow\mathbb{R}$ is bounded from above. 2. A.2

We suppose that $\boldsymbol{X}_{1}\sim\mathcal{U}([0,1]^{d})$. 3. A.3

We suppose that $U_{1}$ is reduced (mainly for the sake of simplicity in exposition) and has a moment of order $4$. 4. A.4

We suppose that $V_{1}$ has a moment of order $4$.

The two following assumptions involving $V_{i}$ and $\boldsymbol{X}_{i}$ are complementary and will be considered separately in the study:

1. A.5

We suppose that $\boldsymbol{X}_{i}$ and $V_{i}$ are independent for any $i\in\{1,\ldots,n\}$, and $U_{1}$ is centered or $V_{1}$ is centered. 2. A.6

We suppose that $V_{i}=g(\boldsymbol{X}_{i})$ where $g:[0,1]^{d}\rightarrow\mathbb{R}$ is known and bounded from above, and $U_{1}$ is centered.

These assumptions will be discussed later; some of them can be relaxed. In our main results, we will consider the two following sets of assumptions:

1. H.1

$=\{\ref{hyp:A1},\ref{hyp:A2},\ref{hyp:A3},\ref{hyp:A4},\ref{hyp:A5}\}$,  $=\{\ref{hyp:A1},\ref{hyp:A2},\ref{hyp:A3},\ref{hyp:A4},\ref{hyp:A6}\}$.

As usual in wavelet methods, the first step towards the estimation of $r$ is to consider its wavelet series given by (2). Then we aim to estimate the unknown wavelet coefficients $\alpha_{j,\boldsymbol{k}}=\langle r,\Phi_{j,\boldsymbol{k}}\rangle_{[0,1]^{d}}$ and $\beta_{j,\boldsymbol{k},u}=\langle r,\Psi_{j,\boldsymbol{k},u}\rangle_{[0,1]^{d}}$ by efficient estimators. In this study, we propose to estimate $\alpha_{j,\boldsymbol{k}}$ by

[TABLE]

where

[TABLE]

This is an unbiased estimator of $\alpha_{j,\boldsymbol{k}}$ and it converges to $\alpha_{j,\boldsymbol{k}}$ in $L^{2}$ (see Lemmas 5.1 and 5.2 in Section 5). On the other side, we propose to estimate $\beta_{j,\boldsymbol{k},u}$ by

[TABLE]

where $\mathbbm{1}_{A}$ denotes the indicator function over an event $A$, $\rho_{n}:=\sqrt{n/\ln n}$ and

[TABLE]

Due to the thresholding in its definition, this estimator is not unbiased of $\beta_{j,\boldsymbol{k},u}$ but it converges to $\beta_{j,\boldsymbol{k},u}$ in $L^{2}$ (see Lemmas 5.1 and 5.2 in Section 5). The role of the thresholding is to relax assumptions on $U_{1}$ and $V_{1}$; note that only moments of order $4$ is required in A.3 and A.4 including uniform or Gaussian distribution. This selection rule has been introduced in a wavelet setting in Delyon and Juditsky (1996). It has been recently improved in Chesneau (2013) (in a multidimensional regression function under mixing dependence framework) and Chaubey et al. (2015) (in a density estimation under multiplicative censoring setting). In this study, we adapt it to the general nonparametric regression model (1).

The next step in the construction of our wavelet estimators for $r$ is to expand the most informative of the wavelet coefficients estimators using the initial wavelet basis. We then define the linear wavelet estimator by

[TABLE]

We thus have projected the $\hat{\alpha}_{j,\boldsymbol{k}}$’s on the father wavelet basis at a certain level $j_{*}$. Despite the simplicity of its construction, this estimator has a serious drawback: its performance highly depends on the choice of the level $j_{*}$. A suitable choice of $j_{*}$, but depending on the smoothness of $r$, will be specified in our main result. To address this problem, an alternative is proposed by using a hard-thresholding rule that performs a term-by-term selection of the wavelet coefficient estimators $\hat{\beta}_{j,\boldsymbol{k},u}$ and to project them on the original wavelet basis. We define the nonlinear wavelet estimator by

[TABLE]

where $t_{n}:=\sqrt{\ln n/n}=\rho_{n}^{-1}$. The positive integer $j_{1}$ is specified in our main result, while the constant $\kappa$ will be chosen in its proof (see the proof of Lemma 5.3). The idea of keeping the estimators $\hat{\beta}_{j,\boldsymbol{k},u}$ with magnitude greater to $t_{n}$ is not new; it is a well-known wavelet techniques with strong mathematical and practical results for numerous nonparametric problems; $t_{n}$ is so-called “universal threshold”. We refer to Donoho et al. (1995), Delyon and Juditsky (1996) and Härdle et al. (2012). In this study, we describe how to calibrate such estimator when we deal with the general model (1).

In the sequel, we adopt the following notations: $x_{+}:=\max\{x,0\}$. $A\lesssim B$ denotes $A\leq cB$ for some constant $c>0$; $A\gtrsim B$ means $B\lesssim A$; $A\thicksim B$ stands for both $A\lesssim B$ and $B\lesssim A$.

Theorem 3.1 below determines the rates of convergence attained by $\hat{r}^{\mathrm{lin}}_{n}$ and $\hat{r}^{\mathrm{non}}_{n}$ over the MISE.

Theorem 3.1**.**

Consider the problem defined by (1) under the assumptions H.1 or H.1, let $r\in B^{s}_{p,q}([0,1]^{d})$ with $p,q\in[1,\infty)$, $s>d/p$. Then

• •

the linear wavelet estimator $\hat{r}^{\mathrm{lin}}_{n}$ with $2^{j_{*}}\thicksim n^{\frac{1}{2s^{\prime}+d}}$ and $s^{\prime}=s-d(1/p-1/2)_{+}$ satisfies

[TABLE]

• •

the nonlinear estimator with $2^{j_{*}}\sim n^{\frac{1}{2m+d}}~{}(m>s)$, $2^{j_{1}}\sim(n/\ln n)^{\frac{1}{d}}$ satisfies

[TABLE]

The obtained rates of convergence are those obtained in the standard density estimation problem or the regression function estimation problem under the MISE over Besov spaces (see Härdle et al. (2012)). Under some strong conditions on the model as $U_{i}:=1$, the rate of convergence $n^{-\frac{2s}{2s+d}}$ is proved to be optimal in the minimax sense (see Härdle et al. (2012) and Tsybakov (2009)). So our nonlinear wavelet estimator can be optimal in the minimax sense up to a $\ln n$. However, in full generality, without specifying the distributions of $U_{i}$ and $V_{i}$, the optimal lower bounds for the MISE are difficult to determine via standard techniques (Fano’s lemma, …) and the optimality of our estimators remains an open question.

Remark 3.1**.**

Some assumptions used in Theorem 3.1 can be relaxed without changing the result. In particular, one can consider the domain $\Delta=[a,b]^{d}$ with $(a,b)\in\mathbb{R}^{2}$ and $a<b$ with an adaptation of the wavelet basis. In this case, we can also replace A.2 by $\boldsymbol{X}_{1}$ with density function $h:[a,b]^{d}\rightarrow\mathbb{R}$ bounded from below, with the following wavelet coefficient estimators:

[TABLE]

and

[TABLE]

Finally, note that, in A.6 can be improved by assuming $g$ unknown. To the best of our knowledge, only Cai et al. (2008) have developed wavelet methods in this case for $d=1$, deterministic design $(\boldsymbol{X}_{i}:=i/n)$ and infinite moments for $U_{i}$. Extension of these methods in the general setting of (1) needs further developments that we leave for a future work.

4 Numerical Experiments

To illustrate the empirical performance of the estimators proposed in this work, we carried out a simulation study. The objective is to highlight some of the theoretical findings using numerical examples. We begin by giving some details about the specificities inherent in wavelet estimators in a non-deterministic design framework. We also try to propose a realistic simulation setting using an adaptive selection method to select both the truncation parameter of the linear estimator and the threshold parameter of the non-linear estimator. In this context, we compare their empirical performances in the model with both multiplicative and additive noise. Simulations were performed using R and in particular the rwavelet package Navarro and Chesneau (2019) (available from https://github.com/fabnavarro/rwavelet).

4.1 Computational aspect of wavelets and parameters selection

For fixed design, thanks to Mallat’s pyramidal algorithm (Mallat, 2008), the computation of wavelet-based estimators is simple and fast. When considering uniform random design, the implementation requires some changes and several strategies have been developed in the literature (see, e.g., Cai et al. (1998); Hall et al. (1997)). For uniform design regression, Cai and Brown (1999) has proposed to use an approach in which the wavelet coefficients are computed by a simple application of Mallat’s algorithm using the ordered $Y_{i}$’s as input variables. We have followed this approach because it preserves the simplicity of calculation and the efficiency of the equispaced algorithm. In the context of wavelet regression in random design with heteroscedastic noise, Kulik et al. (2009) and Navarro and Saumard (2017) also adopted this approach.

Nason successfully adjusted the standard two-Fold Cross Validation (2FCV) method to select the threshold parameter in wavelet shrinkage (see, Nason (1996)). For the calibration of linear wavelet estimators, his strategy was used by Navarro and Saumard (2017). We have chosen to apply this approach to select both the threshold and truncation parameter of linear and non-linear estimators. More precisely, in the linear case, we built a collection of linear estimators $\hat{r}_{j_{*},n}^{\mathrm{lin}},j_{*}=0,1,\ldots,\log 2(n)-1$ (by successively adding whole resolution levels of wavelet coefficients), and select the best among this collection by minimizing a 2FCV criterion denoted by $\mathrm{2FCV}(j_{*})$. The resulting estimator of the truncation level is denoted by $\hat{j}_{*}$ and the corresponding estimator of $r$ by $\hat{r}_{\hat{j}_{*},n}^{\mathrm{lin}}$ (see, Navarro and Saumard (2017, 2018) for more details). For the nonlinear estimator, the same estimator $\hat{j}_{*}$ of the truncation parameter $j_{*}$ obtained for the linear is used. The estimator of the thresholding parameter is obtained using the 2FCV method developed in Nason (1996). The parameter $j_{1}$ is fixed a priori as the maximum level allowed by the wavelet decomposition (i.e., $j_{1}=\log 2(n)-1$). It is a classic choice that allows the coefficients to be selected down to the smallest scale. In addition, in order to facilitate and not to overburden the implementation of the nonlinear estimator, we perform a standard hard thresholding of the wavelet coefficient estimators (rather than the double threshold used in its definition). In order to be able to evaluate the performance of these two criteria, the mean square error (MSE) is used (i.e., $\mathrm{MSE}(\hat{r}_{j_{*},n},r)=\frac{1}{n}\sum_{i=1}^{n}(r(X_{i})-\hat{r}_{j_{*},n}(X_{i}))^{2})$). We consider three test functions for $r$ (see Figure 1), commonly used in the wavelet literature, Parabolas, Ramp and Blip (see, e.g., Donoho et al. (1995)). In all simulations, we examine the case $d=1$, the design is chosen to satisfy A.2 (i.e., $\mathcal{U}([0,1])$) and the choice of the wavelet family used is also fixed (i.e., Daubechies compactly supported wavelet with 8 vanishing moments).

4.2 Additive-multiplicative regression

This subsection examines the behaviour and performance of linear and non-linear estimators in the context of additive and multiplicative regression by considering $V_{1}\sim\mathcal{N}(0,\sigma^{2})$, where $\sigma^{2}=0.01$ and $U_{1}\sim U([-1,1])$. Thus, the goal is to estimate the frontier $r$ from $(X_{i},Y_{i})$ sample simulated from one of the test functions. By applying one of the linear or nonlinear methods developed above to the estimation of $r$, one can construct an estimator whose rate of convergence is given by (5a) and (5b) respectively. Note that here, the nature of the frontier function is not necessarily the same as that commonly found in the literature on stochastic boundary estimation. Indeed, here $r$ is not necessarily a production function (e.g., $r$ is concave), the only assumption we make is given by A.1. Thus, the application here can be seen as the estimation of the boundary or frontier of a sample affected by some additive (positive) noise (see Jirak et al. (2014)).

A typical example of estimation for the Blip function, with $n=4096$ is given in Figure 2. It can be seen that the minimum of $\mathrm{2FCV}(j_{*})$ criteria coincides with that of unknown risk (i.e., $\hat{j}_{*}=4$) and therefore provides the best possible linear estimator for the collection under consideration (i.e., $\mathrm{MSE}(r,\hat{r}_{\hat{j}_{*},n}^{\mathrm{lin}})=0.0027$). We have not included the results of the non-linear estimator in Figure 2. Indeed, in this case, the value of the threshold obtained by minimizing the cross validated criterion leads to the elimination of all the thresholded coefficients (i.e., going from $\hat{j}_{*}$ to $j_{1}$) and therefore leads to the same estimate and the same risk as the linear estimator. We can see (Figure 2(c)) that the unknown risk behaves in the same way here. This is partly because the amplitude of the coefficients at the fine scales is so large and variable from one scale to another that it is not possible to obtain an overall optimal threshold value that makes it possible to maintain certain important coefficients and that, on the contrary, keeps coefficients associated with noise, with this specific thresholding policy (i.e., a ‘keep’ or ‘kill’ rule). In particular the important coefficients located on scales larger than $\hat{j}_{*}$ (especially those encoding the discontinuity of $r$) are too small in amplitude to be maintained by a global threshold.

In order to determine whether this phenomenon observed for a single function and a single realization is confirmed in a more general context, we compare the performance in terms of MSE (computed on the functions after reconstruction) for both estimators and for the three-test functions. For each function, a sample of $N=100$ is generated and we compare the average behavior of the MSE for parameters selected with the oracle obtained by minimizing the MSE using the original signal $r$ (denoted by $\mathrm{MSE}^{\mathrm{lin}}$ and $\mathrm{MSE}^{\mathrm{non}}$ respectively), the linear $\mathrm{2FCV}^{\mathrm{lin}}$ strategy and the non-linear $\mathrm{2FCV}^{\mathrm{non}}$ (i.e., calculated from $\hat{j}_{*}$ and the threshold that minimizes the $\mathrm{2FCV}(\lambda)$) strategy. Figure 3 presents the results in the form of boxplots, one for each function. On the one hand, for all three functions, we can see that the performance of $\mathrm{2FCV}^{\mathrm{lin}}$ is at the $\mathrm{MSE}^{\mathrm{lin}}$ level. This procedure therefore provides a remarkable surrogate of the unknown risk. On the other hand, the non-linear $\mathrm{2FCV}^{\mathrm{non}}$ oracle is similar to $\mathrm{MSE}^{\mathrm{lin}}$, which means that the optimal threshold here leads systematically to the suppression of all threshold coefficients — which corresponds to the selection of values of the threshold parameter which is greater than the largest noisy wavelet coefficient in absolute value. Finally, the variability of the $\mathrm{2FCV}^{\mathrm{non}}$ is high as a result of selected threshold values that are sometimes too small, resulting in the conservation of unnecessary coefficients in the reconstruction. This is because the curves associated with non-linear criteria do not generally allow a single global minimum, but the minimum is reached in the form of a plateau (see Figure 2(c)). In practice, when the minimum is reached on such a plateau, the first element that constitutes it is selected first. This has no influence on $\mathrm{MSE}^{\mathrm{non}}$ but generates this variability of $\mathrm{2FCV}^{\mathrm{non}}$, i.e., when the abscissa of the first point constituting the plateau associated with $\mathrm{2FCV}^{\mathrm{non}}$ is lower than that of $\mathrm{MSE}^{\mathrm{non}}$) Note that to overcome this problem, in the presence of a plateau, we could for example select a threshold value in the middle of it. We have not done so here to emphasize the fact that a cross validation strategy of the global threshold seems ineffective in this setting. It should also be noted that in our simulations, this finding is also verified for other noise levels or sample sizes (the results are generally very similar, so we give only another example by considering a lower number of samples, $n=2048$ and a lower additive noise level $\sigma^{2}=0.025$).

In conclusion, the linear approach seems more appropriate than the non-linear approach in the context of the simulations considered in this study. One way to fully benefit from the non-linear approach would be to consider an optimal threshold selection strategy on a scale by scale basis. The selection procedure used here, based on an interpolation performed in the original domain, does not facilitate this extension. For this purpose it would be necessary, for example, to define an interpolated version of the cross validation method in the wavelet coefficients domain.

5 Auxiliary results and proof of the main result

5.1 Auxiliary results

In this section, we provide some lemmas for the proof of the main Theorem.

Lemma 5.1**.**

Let $j\geq\tau$, $\boldsymbol{k}\in\Lambda_{j}$, $\hat{\alpha}_{j,\boldsymbol{k}}$ be (3). Then, under H.1 or H.1, we have

[TABLE]

Proof of Lemma 5.1. Using the independence assumptions on the random variables, H.1 or H.1, observe that

[TABLE]

and

[TABLE]

Therefore

[TABLE]

Using similar mathematical arguments, since $\int_{[0,1]^{d}}\Psi_{j,\boldsymbol{k},u}(\boldsymbol{x})d\boldsymbol{x}=0$, we have

[TABLE]

We prove the second equality. The proof of Lemma 5.1 is complete. $\hfill\Box$

Lemma 5.2**.**

Let $j\geq\tau$ such that $2^{j}\leq n$, $\boldsymbol{k}\in\Lambda_{j}$, $\hat{\alpha}_{j,\boldsymbol{k}}$ and $\hat{\beta}_{j,\boldsymbol{k},u}$ be (3) and (4) respectively. Then, under H.1 or H.1,

[TABLE]

Proof of Lemma 5.2. Owing to Lemma 5.1 we have $\mathbf{E}[\hat{\alpha}_{j,\boldsymbol{k}}]={\alpha}_{j,\boldsymbol{k}}$. Therefore

[TABLE]

By A.1 and $\mathbf{E}\left[\Phi^{2}_{j,\boldsymbol{k}}(\boldsymbol{X}_{1})\right]=\int_{[0,1]^{d}}\left(\Phi_{j,\boldsymbol{k}}(\boldsymbol{x})\right)^{2}d\boldsymbol{x}=1$, we have $\mathbf{E}\left[f^{4}(\boldsymbol{X}_{1})\Phi^{2}_{j,\boldsymbol{k}}(\boldsymbol{X}_{1})\right]\lesssim 1$. On the other hand, we have

[TABLE]

Thus all the terms in the brackets of (5.1) are bounded from above. The first inequality in Lemma 5.2 is proved.

Now, by the definition of $\hat{\beta}_{j,\boldsymbol{k},u}$, taking $K_{i}:=Y_{i}^{2}\Psi_{j,\boldsymbol{k},u}(\boldsymbol{X}_{i})-w_{j,\boldsymbol{k},u}$ and $D_{i}:=K_{i}\mathbbm{1}_{\{|K_{i}|\leq\rho_{n}\}}-\mathbf{E}\left[K_{i}\mathbbm{1}_{\{|K_{i}|\leq\rho_{n}\}}\right]$ for the sake of simplicity, the second equation in Lemma 5.1 yields

[TABLE]

Hence, using $\mathbf{E}\left[\left((1/n)\sum\limits_{i=1}^{n}D_{i}\right)^{2}\right]=({1}/{n^{2}})\mathbf{V}\left[\sum\limits_{i=1}^{n}D_{i}\right]=\mathbf{V}[D_{1}]/n\leq\mathbf{E}\left[D_{1}^{2}\right]/n$, we have

[TABLE]

Proceeding as for the proof of the first inequality, using the assumptions H.1 or H.1, note that $\mathbf{E}\left[K^{2}_{1}\right]\lesssim\mathbf{E}\left[Y^{4}_{1}\Psi^{2}_{j,\boldsymbol{k},u}(\boldsymbol{X}_{1})\right]+w^{2}_{j,\boldsymbol{k},u}\lesssim 1$ and

[TABLE]

Therefore,

[TABLE]

The second inequality in Lemma 5.2 is proved. This ends the proof of Lemma 5.2. $\hfill\Box$

Lemma 5.3**.**

Let $j\geq\tau$ such that $2^{jd}\lesssim n/\ln n$, $\boldsymbol{k}\in\Lambda_{j}$, $\hat{\beta}_{j,\boldsymbol{k},u}$ be (4). Then there exists a constant $\kappa>1$ such that

[TABLE]

Proof of Lemma 5.3. By the definition of $\hat{\beta}_{j,\boldsymbol{k},u}$, taking $K_{i}:=Y_{i}^{2}\Psi_{j,\boldsymbol{k},u}(\boldsymbol{X}_{i})-w_{j,\boldsymbol{k},u}$ and $D_{i}:=K_{i}\mathbbm{1}_{\{|K_{i}|\leq\rho_{n}\}}-E\left[K_{i}\mathbbm{1}_{\{|K_{i}|\leq\rho_{n}\}}\right]$ for the sake of simplicity, the second equation in Lemma 5.1 yields

[TABLE]

Using (7), there exists $c>0$ such that $\mathbf{E}\left[|K_{1}|\mathbbm{1}_{\{|K_{1}|>\rho_{n}\}}\right]\leq c\sqrt{\ln n/n}$. Then

[TABLE]

Note that $\mathbf{E}[D_{i}]=0$ thanks to Lemma 5.1. According to the proof of Lemma 5.2, $\mathbf{E}[D_{i}^{2}]:=\delta^{2}\lesssim 1$. This with $|D_{i}|\lesssim\sqrt{n/\ln n}$ and Bernstein inequality shows

[TABLE]

Then one choose large enough $\kappa$ such that

[TABLE]

This is the desired conclusion.$\hfill\Box$

5.2 Proof of the main result

This section is devoted to the proof of Theorem 3.1. We prove (5a) and (5b) in turn.

Proof of (5a)  Note that

[TABLE]

It is easy to see that

[TABLE]

According to Lemma 5.2, $|\Lambda_{j_{*}}|\thicksim 2^{j_{*}d}$ and $2^{j_{*}}\thicksim n^{\frac{1}{2s^{\prime}+d}}$,

[TABLE]

When $p\geq 2$, $s^{\prime}=s$. By Hölder inequality and $r\in B_{p,q}^{s}([0,1]^{d})$,

[TABLE]

When $1\leq p<2$ and $s>d/p$, $B_{p,q}^{s}([0,1]^{d})\subseteq B_{2,\infty}^{s^{\prime}}([0,1]^{d})$

[TABLE]

Therefore, in both cases,

[TABLE]

By (8), (9) and (10),

[TABLE]

Proof of (5b) We now follow the lines of (Delyon and Juditsky, 1996, Theorem 2) with adaptation to our statistical setting, by using the definitions of our estimators and the auxiliary results of Section 5.1. By the definitions of $\hat{r}^{\mathrm{lin}}_{n}$ and $\hat{r}^{\mathrm{non}}_{n}$, we have

[TABLE]

Hence,

[TABLE]

where T_{1}:=\mathbf{E}\left[\Big{\|}\hat{r}^{\mathrm{lin}}_{n}-P_{j_{*}}r\Big{\|}^{2}_{2}\right],~{}~{}T_{2}:=\Big{\|}r-P_{j_{1}+1}r\Big{\|}^{2}_{2} and

[TABLE]

According to (9) and $2^{j_{*}}\sim n^{\frac{1}{2m+d}}~{}(m>s)$,

[TABLE]

When $p\geq 2$, by the same arguments as (10) shows T_{2}=\Big{\|}r-P_{j_{1}+1}r\Big{\|}^{2}_{2}\lesssim 2^{-2j_{1}s}. This with $2^{j_{1}}\sim(n/\ln n)^{\frac{1}{d}}$ leads to

[TABLE]

On the other hand, $B_{p,q}^{s}([0,1]^{d})\subseteq B_{2,\infty}^{s+d/2-d/p}([0,1]^{d})$ when $1\leq p<2$ and $s>d/p$. Then

[TABLE]

Hence,

[TABLE]

for each $1\leq p<+\infty$.

The main work for the proof of (5b) is to show

[TABLE]

Note that

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

For $Q_{1}$, one observes that

[TABLE]

thanks to Hölder inequality. By Lemma 5.2, Lemma 5.3 and $|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}|^{2}\lesssim n/\ln n$,

[TABLE]

Then $Q_{1}\lesssim\sum\limits_{j=j_{*}}^{j_{1}}2^{jd}/n^{2}\lesssim 2^{j_{1}d}/n^{2}\lesssim 1/n\leq n^{-\frac{2s}{2s+d}}$, where one uses the choice $2^{j_{1}}\sim(n/\ln n)^{\frac{1}{d}}$. Hence,

[TABLE]

To estimate $Q_{2}$, one defines

[TABLE]

It is easy to see that $2^{j_{*}}\sim n^{\frac{1}{2m+d}}\leq 2^{j^{\prime}}\sim n^{\frac{1}{2s+d}}\leq 2^{j_{1}}\sim(n/\ln n)^{\frac{1}{d}}$. Furthermore, one rewrites

[TABLE]

By Lemma 5.2 and $2^{j^{\prime}}\sim n^{\frac{1}{2s+d}}$,

[TABLE]

On the other hand, it follows from Lemma 5.2 that

[TABLE]

When $p\geq 2$, since $r\in B_{p,q}^{s}([0,1]^{d})$, Lemma 5.2 and $t_{n}=\sqrt{\ln n/n}$,

[TABLE]

When $1\leq p<2$ and $s>d/p$, $B_{p,q}^{s}([0,1]^{d})\subseteq B_{2,\infty}^{s+d/2-d/p}([0,1]^{d})$. Then

[TABLE]

It follows from the upper bounds above that

[TABLE]

Finally, one evaluates $Q_{3}$. Clearly,

[TABLE]

This with the choice of $2^{j^{\prime}}$ shows

[TABLE]

On the other hand, $Q_{32}:=\sum\limits_{j=j^{\prime}+1}^{j_{1}}\sum\limits_{u=1}^{2^{d}-1}\sum\limits_{\boldsymbol{k}\in\Lambda_{j}}\left|\beta_{j,\boldsymbol{k},u}\right|^{2}\mathbbm{1}_{\{|\beta_{j,\boldsymbol{k},u}|\leq 2\kappa t_{n}\}}$. According to the arguments of (5.2), for $p\geq 2$,

[TABLE]

When $1\leq p<2$, $\left|\beta_{j,\boldsymbol{k},u}\right|^{2}\mathbbm{1}_{\{|\beta_{j,\boldsymbol{k},u}|\leq 2\kappa t_{n}\}}\leq\left|\beta_{j,\boldsymbol{k},u}\right|^{p}\left|2\kappa t_{n}\right|^{2-p}$. Then similar to the arguments of (5.2),

[TABLE]

It follows from the inequalities above that

[TABLE]

in both cases. Owing to (11), (12), (15), and (16), we prove that

[TABLE]

which is the desired conclusion. $\hfill\Box$

Acknowledgements

We would like to thank the reviewers for their thoughtful comments that have improved the manuscript.