Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
59,425 characters · 12 sections · 50 citation commands
Nonparametric estimation in a regression model with additive and multiplicative noise
{Keywords:} Nonparametric regression, multiplicative regression models, nonparametric frontier, rates of convergence, wavelets.
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
$f$ is an unknown regression function defined on a subset $\Delta$ of $\mathbb{R}^d$, with $d\ge 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 \lbrace 1, \ldots, n \rbrace$. 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 (ref) since besides the theoretical interest it embodies several potential applications. For example, for $U_i=1$, (ref) 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 hardle, tsybakov and comte, and the references therein. For $V_i=0$, (ref) becomes the standard nonparametric regression model with multiplicative noise. Recent studies can be found in chi,comte 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 chesneau2018linear. In the case where $V_i$ is a function of $\boldsymbol{X}_i$, (ref) becomes the popular nonparametric regression model with multiplicative heteroscedastic noises:
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 muller1987estimation,hall1990variance,hardletsy,wang, brown, cai for fixed design and recently kulik2011nonparametric,verzelen2018adaptive,shen2019optimal 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 (ref), 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 farrell1957measurement,de1984measuring,gijbels1999estimation,daouia2005robust for deterministic frontiers models and fan1996semiparametric,kumbhakar2007nonparametric,simar2011stochastic for stochastic frontier models. For general regression function and general nonparametric setting, we refer to girard2008frontier,girard2013frontier and jirak2014adaptive 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 as well as in speckle noise reduction encounter in particular in synthetic-aperture radar images kuan1985adaptive or in medical ultrasound images rabbani2008speckle,mateo2009finding), where noise sources can be both additive and multiplicative. In this context, one can also cite korostelev2012minimax 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 (ref), 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 anto2 and hardle, and the references therein. To the best of our knowledge, their development for (ref) 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 in the nonparametric regression estimation setting, and recently improved in chesneau (in a multidimensional regression function under mixing dependence framework) and chaubey (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:96. 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 (ref). Section (ref) specifies some assumptions on the model, presents our wavelet estimators and the main results on their performances. Numerical experiments are presented in Section (ref). Section (ref) is devoted to the proofs of the main result.
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:
See meyer 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 daub). Then we consider the wavelet tensor product bases on $[0,1]^d$ as described in cohen. 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^jx_1-k_1, \ldots,2^jx_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^jx_1-k_1, \ldots,2^jx_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 $0$ 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\ge \tau , \ \boldsymbol{k}\in \Lambda_{j}\}$ forms an orthonormal basis of $L^2(\lbrack 0,1 \rbrack^d)$. For any integer $j_* \ge \tau$, a function $h \in L^2(\lbrack 0,1 \rbrack^d)$ can be expressed via $\mathcal{S}$ by the following wavelet series:
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})$,
Besov spaces are important in theory and applications. They have the features to a wide variety of function spaces as H\"{o}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 (ref). 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})$; \quad (2) $\left\lbrace 2^{js}\|P_{j+1}h-P_{j}h\|_{p}\right\rbrace\in l_{q};$ \quad (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
Further details on Besov spaces are given in meyer, triebel and hardle.
We consider the model (ref) with $\Delta=[0,1]^d$ for the sake of simplicity. Additional technical assumptions are formulated below.
The two following assumptions involving $V_i$ and $\boldsymbol{X}_i$ are complementary and will be considered separately in the study:
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:
As usual in wavelet methods, the first step towards the estimation of $r$ is to consider its wavelet series given by (ref). 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
where \[ v_{j,\boldsymbol{k}}:=\left\{
\right. \] This is an unbiased estimator of $\alpha_{j,\boldsymbol{k}}$ and it converges to $\alpha_{j,\boldsymbol{k}}$ in $L^2$ (see Lemmas (ref) and (ref) in Section (ref)). On the other side, we propose to estimate $\beta_{j,\boldsymbol{k},u}$ by
where $\ensuremath{\mathbbm{1}}_A$ denotes the indicator function over an event $A$, $\rho_n:=\sqrt{n/\ln n}$ and \[ w_{j,\boldsymbol{k},u}:=\left\{
\right. \] 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 (ref) and (ref) in Section (ref)). 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 (ref) and (ref) including uniform or Gaussian distribution. This selection rule has been introduced in a wavelet setting in delyon. It has been recently improved in chesneau (in a multidimensional regression function under mixing dependence framework) and chaubey (in a density estimation under multiplicative censoring setting). In this study, we adapt it to the general nonparametric regression model (ref).
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
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
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 (ref)). 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 donoho3, delyon and hardle. In this study, we describe how to calibrate such estimator when we deal with the general model (ref).
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 (ref) below determines the rates of convergence attained by $\hat{r}^{\mathrm{lin}}_{n}$ and $\hat{r}^{\mathrm{non}}_{n}$ over the MISE.
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 hardle). 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 hardle and tsybakov). 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, \ldots) and the optimality of our estimators remains an open question.
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 \proglang{R} and in particular the {\fontseries{b}\selectfont rwavelet} package chesnav (available from \url{https://github.com/fabnavarro/rwavelet}).
For fixed design, thanks to Mallat's pyramidal algorithm mallat:08, 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:98,hall:97). For uniform design regression, cai:99 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:09 and navarro17 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:96). For the calibration of linear wavelet estimators, his strategy was used by navarro17. 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,\log2(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, navarro17,navarro2 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:96. The parameter $j_1$ is fixed a priori as the maximum level allowed by the wavelet decomposition (i.e., $j_1=\log2(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 (ref)), commonly used in the wavelet literature, Parabolas, Ramp and Blip (see, e.g., donoho3). In all simulations, we examine the case $d=1$, the design is chosen to satisfy (ref) (\emph{i.e.}, $\mathcal{U}([0,1])$) and the choice of the wavelet family used is also fixed (\emph{i.e.}, Daubechies compactly supported wavelet with 8 vanishing moments).
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 (ref) and (ref) 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 (ref). 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 jirak2014adaptive).
A typical example of estimation for the Blip function, with $n=4096$ is given in Figure (ref). 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 (ref). 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 (ref)(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 (ref) 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 (ref)(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.
In this section, we provide some lemmas for the proof of the main Theorem.
{\bf Proof of Lemma (ref).} Using the independence assumptions on the random variables, (ref) or (ref), observe that
and
Therefore
Using similar mathematical arguments, since $\int_{[0,1]^d} \Psi_{j,\boldsymbol{k},u}(\boldsymbol{x})d\boldsymbol{x}=0$, we have
We prove the second equality. The proof of Lemma (ref) is complete. $\hfill \Box$
{\bf Proof of Lemma (ref).} Owing to Lemma (ref) we have $\mathbf{E}[\hat{\alpha}_{j,\boldsymbol{k}}]={\alpha}_{j,\boldsymbol{k}}$. Therefore
By (ref) 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)^2d\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
Thus all the terms in the brackets of (ref) are bounded from above. The first inequality in Lemma (ref) 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\ensuremath{\mathbbm{1}}_{\{|K_i|\leq\rho_n\}}-\mathbf{E}\left[ K_i\ensuremath{\mathbbm{1}}_{\{|K_i|\leq\rho_n\}}\right]$ for the sake of simplicity, the second equation in Lemma (ref) yields
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\le \mathbf{E}\left[D_1^2\right]/n$, we have
Proceeding as for the proof of the first inequality, using the assumptions (ref) or (ref), 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}\lesssim1$ and
Therefore, \[ \mathbf{E}\left[(\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u})^{2}\right]\lesssim \frac{1}{n}+\frac{\ln n}{n}\lesssim\frac{\ln n}{n}. \] The second inequality in Lemma (ref) is proved. This ends the proof of Lemma (ref). $\hfill \Box$
{\bf Proof of Lemma (ref).} 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\ensuremath{\mathbbm{1}}_{\{|K_i|\leq\rho_n\}}-E\left[ K_i\ensuremath{\mathbbm{1}}_{\{|K_i|\leq\rho_n\}}\right]$ for the sake of simplicity, the second equation in Lemma (ref) yields
Using (ref), there exists $c>0$ such that $\mathbf{E}\left[ |K_1|\ensuremath{\mathbbm{1}}_{\{|K_1|>\rho_n\}}\right]\le c \sqrt{\ln n/n}$. Then
Note that $\mathbf{E}[D_{i}]=0$ thanks to Lemma (ref). According to the proof of Lemma (ref), $\mathbf{E}[D_{i}^{2}]:=\delta^{2}\lesssim1$. This with $|D_{i}|\lesssim \sqrt{n/\ln n}$ and Bernstein inequality shows
Then one choose large enough $\kappa$ such that
This is the desired conclusion.$\hfill \Box$
This section is devoted to the proof of Theorem (ref). We prove (ref) and (ref) in turn.
Proof of (ref) Note that
It is easy to see that \[ \mathbf{E}\left[\left\|\hat{r}^{\mathrm{lin}}_{n}-P_{j_{*}}r\right\|^{2}_{2}\right]=\mathbf{E}\left[\left\|\sum\limits_{\boldsymbol{k}\in\Lambda_{j_{*}}}(\hat{\alpha}_{j_{*},\boldsymbol{k}}-\alpha_{j_{*},\boldsymbol{k}})\Phi_{j_{*},\boldsymbol{k}}\right\|^{2}_{2}\right] =\sum\limits_{\boldsymbol{k}\in\Lambda_{j_{*}}} \mathbf{E}\left[\Big|\hat{\alpha}_{j_{*},\boldsymbol{k}}-\alpha_{j_{*},\boldsymbol{k}}\Big|^{2}\right]. \] According to Lemma (ref), $|\Lambda_{j_{*}}|\thicksim2^{j_{*}d}$ and $2^{j_{*}}\thicksim n^{\frac{1}{2s'+d}}$,
When $p\geq2$, $s'=s$. By H\"{o}lder inequality and $r\in B_{p,q}^{s}([0,1]^{d})$,
When $1\leq p<2$ and $s>d/p$, $B_{p,q}^{s}([0,1]^{d})\subseteq B_{2,\infty}^{s'}([0,1]^{d})$
Therefore, in both cases,
By (ref), (ref) and (ref),
Proof of (ref) We now follow the lines of delyon with adaptation to our statistical setting, by using the definitions of our estimators and the auxiliary results of Section (ref). By the definitions of $\hat{r}^{\mathrm{lin}}_{n}$ and $\hat{r}^{\mathrm{non}}_{n}$, we have
Hence,
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 \[ Q:=\mathbf{E}\left[\left\|\sum\limits_{j=j_{*}}^{j_{1}} \sum\limits_{u=1}^{2^{d}-1}\sum\limits_{\boldsymbol{k}\in\Lambda_j}\left(\hat{\beta}_{j,\boldsymbol{k},u}\ensuremath{\mathbbm{1}}_{\{|\hat{\beta}_{j,\boldsymbol{k},u}|\geq\kappa t_{n}\}}-\beta_{j,\boldsymbol{k},u}\right)\Psi_{j,\boldsymbol{k},u}\right\|^{2}_{2}\right]. \] According to (ref) and $2^{j_{*}}\sim n^{\frac{1}{2m+d}}~(m>s)$, \[ T_{1}=\mathbf{E}\left[\Big\|\hat{r}^{\mathrm{lin}}_{n}-P_{j_{*}}r\Big\|_{2}^{2}\right]\lesssim \frac{2^{j_{*}d}}{n}\thicksim n^{-\frac{2m}{2m+d}}<n^{-\frac{2s}{2s+d}}. \] When $p\geq2$, by the same arguments as (ref) shows $T_{2}=\Big\|r-P_{j_{1}+1}r\Big\|^{2}_{2}\lesssim2^{-2j_{1}s}.$ This with $2^{j_{1}}\sim(n/\ln n)^{\frac{1}{d}}$ leads to
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
Hence, \[ T_{2}\lesssim(\ln n)n^{-\frac{2s}{2s+d}}, \] for each $1\leq p<+\infty$.
The main work for the proof of (ref) is to show
Note that
where \[ Q_{1}=\sum\limits_{j=j_{*}}^{j_{1}}\sum\limits_{u=1}^{2^{d}-1}\sum\limits_{\boldsymbol{k}\in\Lambda_{j}}\mathbf{E}\left[\left|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}\right|^{2}\ensuremath{\mathbbm{1}}_{\{|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}|>\frac{\kappa t_{n}}{2}\}}\right], \] \[ Q_{2}=\sum\limits_{j=j_{*}}^{j_{1}}\sum\limits_{u=1}^{2^{d}-1}\sum\limits_{\boldsymbol{k}\in\Lambda_{j}}\mathbf{E}\left[\left|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}\right|^{2}\ensuremath{\mathbbm{1}}_{\{|\beta_{j,\boldsymbol{k},u}|\geq\frac{\kappa t_{n}}{2}\}}\right], \] \[ Q_{3}=\sum\limits_{j=j_{*}}^{j_{1}}\sum\limits_{u=1}^{2^{d}-1}\sum\limits_{\boldsymbol{k}\in\Lambda_{j}}\left|\beta_{j,\boldsymbol{k},u}\right|^{2}\ensuremath{\mathbbm{1}}_{\{|\beta_{j,\boldsymbol{k},u}|\leq2\kappa t_{n}\}}. \] For $Q_{1}$, one observes that \[ \mathbf{E}\left[\left|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}\right|^{2}\ensuremath{\mathbbm{1}}_{\{|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}|>\frac{\kappa t_{n}}{2}\}}\right]\leq\left(\mathbf{E}\left[\left|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}\right|^{4}\right]\right)^{\frac{1}{2}}\left(\mathbb{P} \left(|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}|>\frac{\kappa t_{n}}{2}\right)\right)^{\frac{1}{2}} \] thanks to H\"{o}lder inequality. By Lemma (ref), Lemma (ref) and $|\hat{\beta}_{j,\boldsymbol{k},u}-\beta_{j,\boldsymbol{k},u}|^{2}\lesssim n/\ln n$,
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,
To estimate $Q_{2}$, one defines \[ 2^{j'}\sim n^{\frac{1}{2s+d}}. \] It is easy to see that $2^{j_{*}}\sim n^{\frac{1}{2m+d}}\leq2^{j'}\sim n^{\frac{1}{2s+d}}\leq2^{j_{1}}\sim(n/\ln n)^{\frac{1}{d}}$. Furthermore, one rewrites
By Lemma (ref) and $2^{j'}\sim n^{\frac{1}{2s+d}}$,
On the other hand, it follows from Lemma (ref) that
When $p\geq2$, since $r\in B_{p,q}^{s}([0,1]^{d})$, Lemma (ref) and $t_{n}=\sqrt{\ln n / n}$,
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
It follows from the upper bounds above that
Finally, one evaluates $Q_{3}$. Clearly,
This with the choice of $2^{j'}$ shows
On the other hand, $Q_{32}:=\sum\limits_{j=j'+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}\ensuremath{\mathbbm{1}}_{\{|\beta_{j,\boldsymbol{k},u}|\leq2\kappa t_{n}\}}$. According to the arguments of (ref), for $p\geq2$,
When $1\leq p<2$, $\left|\beta_{j,\boldsymbol{k},u}\right|^{2}\ensuremath{\mathbbm{1}}_{\{|\beta_{j,\boldsymbol{k},u}|\leq2\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 (ref),
It follows from the inequalities above that
in both cases. Owing to (ref), (ref), (ref), and (ref), we prove that
which is the desired conclusion. $\hfill \Box$
We would like to thank the reviewers for their thoughtful comments that have improved the manuscript.