# Optimal Almost Sure Rate of Convergence for the Wavelets Estimator in the Partially Linear Additive Models

Khalid Chokri, Salim Bouzebda

► To cite this version:

Khalid Chokri, Salim Bouzebda. Optimal Almost Sure Rate of Convergence for the Wavelets Estimator in the Partially Linear Additive Models. Symmetry, 2025, 17 (3), pp.394. 10.3390/sym17030394. hal-04979046

HAL Id: hal-04979046

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

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal Almost Sure Rate of Convergence for the Wavelets Estimator in the Partially Linear Additive Models

Khalid Chokri 1,† and Salim Bouzebda 2,*,†

**Abstract:** In this article, we examine a class of partially linear additive models (PLAM) defined via a measurable mapping $\Psi: \mathbb{R}^q \to \mathbb{R}$. More precisely, we consider $\Psi(\mathbf{Y}_i) := \mathcal{Y}_i = \mathbf{Z}_i^\top \beta + \sum_{\ell=1}^d m_\ell(X_{\ell,i}) + \varepsilon_i$, $i = 1, \dots, n$, where $\mathbf{Z}_i = (Z_{i,1}, \dots, Z_{i,p})^\top$ and $\mathbf{X}_i = (X_{1,i}, \dots, X_{d,i})^\top$ denote vectors of explanatory variables. The unknown parameter vector is $\beta = (\beta_1, \dots, \beta_p)^\top$, and $m_1, \dots, m_d$ are real-valued functions of a single variable whose forms are not specified. The error terms $\varepsilon_1, \dots, \varepsilon_n$ are identically distributed with mean zero and finite variance $\sigma_\varepsilon$, and they fulfill the condition $\mathbb{E}(\varepsilon | \mathbf{X}, \mathbf{Z}) = 0$ almost surely. These models are broadly applicable in finance, biology, and engineering, where capturing intricate nonlinear effects is essential. We propose an estimation method that leverages marginal integration in conjunction with linear wavelet-based techniques to obtain estimators for the unknown components $m_1, \dots, m_d$. Under suitable regularity conditions, we establish strong uniform convergence of these estimators, demonstrating that they achieve practically relevant convergence rates. Our theoretical results indicate that these estimators converge uniformly at rates that are favorable for practical applications, underscoring the adaptability and scope of this partially linear additive model.

**Keywords:** additive model; regression function; marginal integration; non-parametric estimation; density estimation; regression estimation; wavelet estimation; linear estimators; linear models

**MSC:** 62G20; 62G32; 62J05; 60G07; 60F15; 62G08

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

check for
updates

Academic Editor: Junesang Choi

Received: 10 February 2025

Revised: 28 February 2025

Accepted: 4 March 2025

Published: 5 March 2025

**Citation:** Chokri, K.; Bouzebda, S. Optimal Almost Sure Rate of Convergence for the Wavelets Estimator in the Partially Linear Additive Models. *Symmetry* **2025**, 17, 394. https://doi.org/10.3390/sym17030394

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Consider a random vector $(\mathbf{X}, \mathbf{Y}, \mathbf{Z}) \in \mathbb{R}^d \times \mathbb{R}^q \times \mathbb{R}^p$, and let $\Psi: \mathbb{R}^q \to \mathbb{R}$ be a measurable mapping from some function class $\mathcal{F}$. We investigate a modeling structure where the relation between the response variables and the covariates is captured by a partially linear specification of the form

$$ \Psi(\mathbf{Y}) = \mathcal{Y} = m_0 + \mathbf{Z}^\top \beta + m(\mathbf{X}) + \varepsilon, \quad (1) $$

where $m_0$ is a constant, $\beta \in \mathbb{R}^p$ is a vector of unknown coefficients, and $m(\cdot)$ is a nonparametric function designed to capture potential nonlinear effects of the variables in $\mathbf{X}$. The random error term $\varepsilon$ represents noise. Throughout, the transpose of any vector $\mathbf{V}$ will be written as $\mathbf{V}^\top$. Model (1) is notably flexible, accommodating a variety of established statistical frameworks as particular cases (see, for example, [1]). In practice, it is common

<footnote n="1">¹ Modeling and Complex Systems Laboratory-M.C.S.L., Cadi Ayyad University, P.B. 549, Marrakech 40001, Morocco</footnote>

<footnote n="2">² Laboratoire de Mathématiques Appliquées de Compiègne-L.M.A.C., Université de Technologie de Compiègne, 60205 Compiègne Cedex, France</footnote>

<footnote>* Correspondence: salim.bouzebda@utc.fr</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
to collect a large set of variables, which then requires the elimination of those that are not influential before finalizing the model. Hence, the selection of significant variables is a key consideration in both parametric and nonparametric regression contexts; for a comprehensive discussion, see [2]. One illustrative application concerns the study of factors influencing East Germans' decisions to migrate to West Germany one year after reunification [3,4], where certain covariates (e.g., discrete variables) may exert linear effects on the response, whereas others act in a nonlinear manner. Similar scenarios arise in various domains, such as civil engineering materials [5], the Pima Indian diabetes dataset [6,7], and investigations related to guidelines for urinary incontinence discussion and evaluation [8].

Introducing a linear component $Z^\top \beta$ alongside a nonparametric function $m(X)$ in the model (1) offers an appealing compromise between purely parametric and fully nonparametric approaches. Parametric models are easy to use but can be severely biased if their assumptions are violated, whereas nonparametric techniques can accommodate complex data structures but often need large sample sizes for reliable inference. By incorporating both a linear and a nonparametric term, the partially linear formulation in (1) combines the interpretability and efficiency of parametric methods with the flexibility to capture nonlinear relationships via $m(X)$. A salient special case is isotonic regression, which arises if $\Psi(y) = y$ and the linear term is omitted. This isotonic setting dates back to the pioneering work of [9] and remains a rich area of investigation in statistics.

The partially linear framework, initially proposed by [10] for the classical situation $\Psi(y) = y$, has been extensively used in fields such as economics, epidemiology, engineering, and the life sciences. Foundational contributions by [11–13] deepened the theoretical underpinnings for estimation and inference, motivating a large body of ensuing research. In recent decades, numerous efforts have extended (1) to the case $\Psi(y) = y$ and refined inference for both the parametric and nonparametric components; see [14–22], as well as the comprehensive volume [23]. More contemporary developments appear in [24–33].

When parametric models are correctly specified, they can yield a strong performance but may introduce significant biases if crucial assumptions are violated. In contrast, fully nonparametric smoothing does not rely on strict parametric constraints, though it can be hindered by the “curse of dimensionality” in moderate or high dimensions. To address this challenge, refs. [34,35] introduced additive models, which decompose multivariate nonparametric relationships into sums of lower-dimensional functions. This decomposition promotes interpretability, reduces computational overhead, and can lead to performance gains over fully nonparametric techniques. Detailed expositions of additive models and their development appear in [36,37], and related sources.

In our specific context, we handle the potentially high dimension of $m(X)$ in (1) by imposing an additive structure. Concretely, we replace $m(X)$ with

$$ \Psi(\mathbf{Y}) = m_0 + Z^\top \beta + \sum_{\ell=1}^{d} m_\ell(X_\ell) + \varepsilon, \quad (2) $$

where $X_\ell$ denotes the $\ell$-th component of $X$ and $m_\ell(\cdot)$ captures its influence on $\Psi(\mathbf{Y})$. This partially linear additive design retains a transparent linear part $Z^\top \beta$ while permitting each covariate $X_\ell$ to appear in a potentially nonlinear form. Empirical results such as those in [38] suggest that exploiting this additive framework can substantially improve efficiency. In related work, refs. [39,40] illustrate how additional constraints in semiparametric models can also boost estimation precision.

Researchers have long taken an interest in partially linear additive models across diverse domains. Kernel-based techniques are frequently employed to estimate (2), partly because of their relative ease of implementation and strong theoretical foundations. In-depth reviews of kernel-based methods can be found in [41–43]. In parallel, ref. [30]
introduced a fusion-penalized inverse probability weighted least squares method designed for missing responses in partially linear additive setups with latent subgroups. Meanwhile, ref. [32] proposed a smooth-threshold generalized method of moments estimator and a variable selection procedure for partially linear additive spatial autoregressive models, utilizing basis expansions for the nonparametric components.

Recent work by [29] examined partially linear additive models in pooled data settings, supporting their approach with simulations and an application based on environmental health information from the National Health and Nutrition Examination Survey. Independently, ref. [26] introduced a weighted quantile regression methodology with inverse probability weighting to address selective data retention, establishing large-sample behavior for both parametric and nonparametric components. Other significant studies include [44], which proposed a scalar-on-function paradigm incorporating nonparametric functional predictors and semiparametric covariates, and [45], which developed generalized partial functional linear additive models (GPFLAM) and employed a quasi-likelihood approach. Additionally, ref. [46] investigated partially linear additive dynamic panels with fixed effects and autocorrelation, and [47] examined transformation models for right-censored survival data under the same partially linear additive framework. Wavelet-based methods have also been explored; for example, ref. [48] advanced a spline-backfitted kernel smoothing algorithm displaying oracle-type properties under stationarity and geometric mixing conditions.

Although kernel-based estimators are effective, they may face complications near boundaries or with discontinuous target functions. Wavelet methods inherently address such localized irregularities by allowing for abrupt changes and flexible boundary treatments. Foundational contributions to wavelet theory—see [49–52] and the summary in [53]—established a robust mathematical basis for wavelet-based estimation. Important theoretical expansions of this framework can be found in [54–57].

Wavelet techniques have seen extensive application in density estimation and regression. For instance, ref. [58] analyzed a linear wavelet estimator, proving its strong consistency under $L_2$-type smoothing with the help of sophisticated empirical process tools. Subsequently, ref. [59] introduced statistical inference procedures for wavelet estimators of symmetric positive-definite function curves, employing a log-Euclidean metric to maintain key properties such as positive-definiteness and permutation invariance. Wavelet-based methods have also been used to detect and categorize local symmetry in images [60], to characterize refinement masks for symmetric wavelets meeting particular sum rules [61], and to examine DNA sequence symmetries [62]. Further advances include [63], which developed nearly frame-like wavelets with specified approximation properties, and [64], focusing on block thresholding for nonparametric regression. Broader contexts of wavelet transforms involving symmetry groups or tube domains are discussed in [65,66], while [67] investigated perfect reconstruction of two-dimensional signals on hexagonal grids with enhanced regularity. In [68,69], a methodology based on the $k$-th correlation coefficient was introduced to verify whether the true model error in partial linear settings follows a symmetric distribution.

A notable takeaway from this body of literature is that linear wavelet estimation can asymptotically achieve performance comparable to the well-known nonlinear hard-threshold wavelet estimator developed by [56], reinforcing the practical advantages of wavelet-based approaches. Under $\alpha$-mixing processes, ref. [70] showed that wavelet estimators for densities on bounded domains in $\mathbb{R}^d$ attain optimal or near-optimal almost sure convergence rates, mirroring those in the i.i.d. context. Extending this line of inquiry, ref. [71] demonstrated strong uniform convergence for wavelet regression estimators on compact subsets of $\mathbb{R}^d$. Discussions of linear wavelet estimators appear in [70–74].
That said, the distributional properties of wavelet estimators are comparatively less well-established. Early progress was made by [75], who proved almost sure uniform convergence and asymptotic normality for wavelet estimators of density and regression functions under ergodic $\mathbb{R}^d$-valued processes, later extending these insights to continuous-time ergodic processes in [76]. Subsequent work has refined and broadened these findings; see, for example, [77–80]. Collectively, this evolving literature underscores the versatility and strength of wavelet-based methods in a wide variety of applications.

**Contributions of this paper.** In this study, we seek to establish the optimal almost sure convergence rate of wavelet estimators for each component function in partially linear additive regression models (PLAM). Our work substantially extends and refines earlier theoretical findings, in particular those reported by [43,80–82]. Rather than merely synthesizing existing methodologies, we rigorously examine the dual challenges posed by mixing dependence structures—commonly arising in both time series and spatial data—and the complexities introduced by wavelet expansions. To meet these challenges, we employ advanced large-sample theory, the mixing conditions proposed by [71], and the sophisticated wavelet expansion techniques discussed in [79]. By integrating these elements, we can derive optimal almost sure convergence results that remain valid beyond the standard assumption of independent and identically distributed observations. Crucially, even in the simpler i.i.d. setting, the optimal almost sure convergence rate of wavelet estimators for component functions in additive regression models has not been fully determined. Addressing this unresolved issue constitutes a key contribution to our work and highlights the significant theoretical and technical developments required to advance the field. In previous works [71, 75, 76, 78, 79, 83], wavelet-based techniques have been employed to address multivariate regression. However, in these studies, the convergence rate was shown to depend on the dimension of the covariates, causing performance to deteriorate as the number of covariates increases. As a result, larger sample sizes become necessary to draw robust conclusions. In contrast, the principal contribution of the present paper lies in achieving a dimension-free convergence rate, which offers a significant advantage in practical applications.

**Paper organization.** The subsequent sections provide the essential background, introduce the primary model, and present the theoretical findings. In Section 2, we offer a concise overview of wavelets, Besov spaces, and wavelet estimators, which collectively serve as the foundational framework for our later developments. Section 3 then elaborates on the partially linear additive model and specifies the estimators that constitute the core of our analysis. Next, Section 4 outlines the principal assumptions and formally states our main theoretical contributions, including detailed assertions regarding the asymptotic normality of each additive component. Section 5 concludes by underscoring the key results and proposing directions for further study. To maintain a coherent narrative, the complete proofs of the central theorems are deferred to Section 6, while the Appendix A consolidates additional lemmas and other technical derivations, thereby preserving the continuity of the main exposition.

## 2. Wavelets and Besov Space

In this section, we introduce the notation associated with wavelets and Besov spaces. We begin by providing a brief overview of wavelet multiresolution theory as developed in [49], and further elaborated in [70–72]. Consider a multiresolution analysis $\{V_j\}_{j=1}^\infty$ of the space $L^2(\mathbb{R}^d)$, where $L^p(\mathbb{R}^d)$ (for $1 \le p < \infty$) is the set of all measurable functions $f: \mathbb{R}^d \to \mathbb{R}$ satisfying

$$ \int_{\mathbb{R}^d} |f(x)|^p dx < \infty $$
and is equipped with the norm

$$ \|f\|_{\mathbf{L}^p} = \left( \int_{\mathbb{R}^d} |f(x)|^p dx \right)^{1/p}. $$

A multiresolution analysis is defined as a decomposition of $\mathbf{L}^2(\mathbb{R}^d)$ into an increasing sequence of closed subspaces $\{V_j : j \in \mathbb{Z}\}$, satisfying the following properties:

1. For all integers $j$,

$$ V_j \subseteq V_{j+1}. $$

2.

$$ \bigcap_j V_j = \{0\}, \quad \overline{\bigcup_j V_j} = \mathbf{L}^2(\mathbb{R}^d). $$

3. For any $f(\cdot) \in V_j$,

$$ f(2\cdot) \in V_{j+1}, \quad \text{and} \quad f(\cdot + k) \in V_j \text{ for all } k \in \mathbb{Z}^d. $$

4. There exists a scaling function $\phi \in \mathbf{L}^2(\mathbb{R}^d)$, normalized such that

$$ \int_{\mathbb{R}^d} \phi(x) dx = 1, $$

and the set

$$ \{\phi_k(x) = \phi(x - k) : k \in \mathbb{Z}^d\} $$

forms an orthonormal basis of $V_0$. Hence,

$$ V_0 = \overline{\operatorname{span}}\{\phi(\cdot - k) : k \in \mathbb{Z}^d\}. $$

By scaling, we obtain an orthonormal basis for $V_j$:

$$ \{\phi_{j,k}(x) = 2^{jd/2}\phi(2^j x - k) : k \in \mathbb{Z}^d\}. $$

Throughout this paper, we assume that $\phi(\cdot)$ is $r$-regular ($r \ge 1$), meaning that $\phi \in C^r(\mathbb{R}^d)$ and satisfies, for each $i > 0$ and all multi-indices $\beta = (\beta_1, \dots, \beta_d)$ with

$$ |\beta| = \sum_{i=1}^{d} \beta_i \le r, $$

the bound

$$ |D^\beta \phi(x)| \le \frac{C_i}{(1 + \|x\|)^i}, \quad (3) $$

where

$$ D^\beta \phi(x) = \frac{\partial^{|\beta|} \phi}{\partial x_1^{\beta_1} \dots \partial x_d^{\beta_d}}(x), $$

and $C_i$ is a constant depending only on $i$.

Define the wavelet space $W_j$ by

$$ V_j \oplus W_j = V_{j+1}. $$

Following the construction in [50], one obtains $N = 2^d - 1$ wavelet functions

$$ \{\psi_i : i = 1, \dots, N\} $$
associated with the scaling function $\phi$, such that the following applies:

1. Orthonormal Basis for $W_0$:

$$ \{\psi_i(x - k) : k \in \mathbb{Z}^d, i = 1, \dots, N\} $$

forms an orthonormal basis for $W_0$.

2. Orthonormal Basis for $\mathbf{L}^2(\mathbb{R}^d)$:

$$ \{\psi_{i,j,k}(x) = 2^{jd/2}\psi_i(2^jx - k) : j \in \mathbb{Z}, k \in \mathbb{Z}^d, i = 1, \dots, N\} $$

forms an orthonormal basis for $\mathbf{L}^2(\mathbb{R}^d)$.

3. Regularity and Compact Support: Each $\psi_i(\cdot)$ shares the same regularity as $\phi(\cdot)$, and both $\phi$ and $\psi_i$ have compact support contained within $[-L, L]^d$ for some $L > 0$.

Starting from some integer $j_0$, any $f(\cdot) \in \mathbf{L}^2(\mathbb{R}^d)$ can be expanded as

$$ f(x) = \sum_{k \in \mathbb{Z}^d} a_{j_0,k} \phi_{j_0,k}(x) + \sum_{i=1}^{2^d-1} \sum_{j=j_0}^{\infty} \sum_{k \in \mathbb{Z}^d} b_{i,j,k} \psi_{i,j,k}(x), \quad (4) $$

where

$$ a_{j_0,k} = \int_{\mathbb{R}^d} f(u) \phi_{j_0,k}(u) du, $$

and

$$ b_{i,j,k} = \int_{\mathbb{R}^d} f(u) \psi_{i,j,k}(u) du. $$

For any $\ell \ge j_0$, the orthogonal projection of $f$ onto $V_\ell$ can be written in two equivalent forms:

$$ \begin{aligned} (P_{V_\ell} f)(x) &= \sum_{k \in \mathbb{Z}^d} a_{\ell,k} \phi_{\ell,k}(x) \\ &= \sum_{k \in \mathbb{Z}^d} a_{j_0,k} \phi_{j_0,k}(x) + \sum_{j=j_0}^{\ell} \sum_{i=1}^{N} \sum_{k \in \mathbb{Z}^d} b_{i,j,k} \psi_{i,j,k}(x). \end{aligned} $$

Besov spaces $B_{s,p,q}$ are of central importance in statistical estimation and approximation theory as they provide a flexible framework for characterizing smoothness properties; for instance, see [84–88]. They include a wide variety of both homogeneous and inhomogeneous function spaces frequently employed in statistical analysis. Using wavelet coefficients, ref. [49] characterizes the Besov space as follows. Consider a real-valued smoothness parameter $0 < s < r$ and let $f \in \mathbf{L}^p(\mathbb{R}^d)$. Then $f \in B_{s,p,q}$ if and only if

(B.1)

$$ J_{s,p,q}(f) = \|P_{V_0}f\|_{\mathbf{L}^p} + \left( \sum_{j \ge 0} (2^{js} \|P_{W_j}f\|_{\mathbf{L}^p})^q \right)^{1/q} < \infty, $$

(B.2)

$$ J'_{s,p,q}(f) = \|a_0 \cdot \|_{l_p} + \left( \sum_{j \ge 0} (2^{j(s+d(1/2-1/p))} \|b_j \cdot \|_{l_p})^q \right)^{1/q} < \infty, $$

where

$$ \|a_0 \cdot \|_{l_p} = \left( \sum_{k \in \mathbb{Z}^d} |a_{0,k}|^p \right)^{1/p}, $$

and

$$ \|b_j \cdot \|_{l_p} = \left( \sum_{i=1}^{N} \sum_{k \in \mathbb{Z}^d} |b_{i,j,k}|^p \right)^{1/p}. $$
For $q = \infty$, the usual sup-norm modification applies. Additional equivalent characterizations and advantages of Besov spaces in approximation theory and statistics can be found in [56,89–91] and in the Appendix A.

## Linear Wavelets Regression Estimators

We begin by introducing the notation and definitions that will be used in the following sections. We start with a reminder of the definition of strong mixing. To this end, we introduce the following notation.

Let $\mathcal{F}_i^k(Z)$ denote the $\sigma$-algebra generated by the random variables $\{Z_j : i \le j \le k\}$. The strong mixing coefficient as introduced by [92] is defined as follows. For additional details, the reader may consult [93].

**Definition 1.** Let $Z = \{Z_i : i = 1, 2, \dots\}$ be a strictly stationary sequence of random variables. For a positive integer $n$, we set

$$ \alpha(n) = \sup \{ |\mathbb{P}(A \cap B) - \mathbb{P}(A)\mathbb{P}(B)| : A \in \mathcal{F}_1^k(Z), B \in \mathcal{F}_{k+n}^\infty(Z), k \in \mathbb{N}^* \}. $$

The sequence $Z$ is said to be $\alpha$-mixing (or strongly mixing) if the mixing coefficient $\alpha(n)$ converges to 0 as $n \to \infty$.

In the study of mixing conditions, $\alpha$-mixing emerges as a relatively weak yet broadly applicable criterion, encompassing a wide range of stochastic processes, including many time series models. The seminal contributions of [94,95] established sufficient conditions for linear processes to be $\alpha$-mixing. These works further demonstrated that both linear autoregressive and bilinear models possess strong mixing properties, with mixing coefficients that decay at an exponential rate. In addition, ref. [96] highlighted the critical role of $\alpha$-mixing in understanding nonlinear time series models, including aspects of geometric ergodicity, a concept further explored and extended in [75]. Along similar lines, ref. [97] showed that functional autoregressive processes can achieve geometric ergodicity under appropriate conditions.

Moreover, the studies [98,99] demonstrated that, even with minimal assumptions, autoregressive conditional heteroscedastic (ARCH) processes and nonlinear additive autoregressive models incorporating exogenous variables remain both stationary and $\alpha$-mixing.

Next, let $\{X_i, Y_i\}$ be jointly stationary processes, and let $\Psi(\cdot)$ be a Borel measurable function on the real line. Define

$$ m(\mathbf{x}, \Psi) = \mathbb{E}[\Psi(Y_1) | \mathbf{X}_1 = \mathbf{x}], \qquad (5) $$

for $\mathbf{x} \in \mathbb{R}^d$, whenever it exists. That is, we assume

$$ \mathbb{E}(|\Psi(Y_1)|) < \infty. $$

Our framework includes several special cases of the regression model (5), depending on the choice of the function $\Psi(\cdot)$. Notable examples include the following:

* $\Psi(Y) = 1\{Y \le y\}$, which yields the conditional distribution function of $Y_1$ given $X_1 = \mathbf{x}$.

* $\Psi(Y) = Y^k$, which provides the $k$-th conditional moment of $Y_1$ given $X_1 = \mathbf{x}$.
Throughout, we shall assume

$$ \mathbb{E}(|\Psi(Y_1)|^p) < \infty, \quad p \ge 2. \qquad (6) $$

Following the approach of previous works, our estimator of $m(\Psi, \mathbf{x})$ will be obtained as the ratio of wavelet estimators of

$$ g(\Psi, \mathbf{x}) = m(\mathbf{x}, \Psi) f_X(\mathbf{x}) $$

and $f_X(\mathbf{x})$, where $f_X(\mathbf{x})$ is the unknown density function of $\mathbf{X}$. First, we recall that from (4), the linear wavelet estimator of $f_X(\mathbf{x}) \in L^2(\mathbb{R}^d)$ is given by

$$ \hat{f}_{\mathbf{X}}(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}_{\tau, \mathbf{k}} \phi_{\tau, \mathbf{k}}(\mathbf{x}), \qquad (7) $$

or equivalently,

$$ \hat{f}_{\mathbf{X}}(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}_{j_0, \mathbf{k}} \phi_{j_0, \mathbf{k}}(\mathbf{x}) + \sum_{j=j_0}^{\tau} \sum_{i=1}^{N} \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{b}_{ij\mathbf{k}} \psi_{i,j,\mathbf{k}}(\mathbf{x}), \qquad (8) $$

where the estimators of the coefficients $\{a_{\tau, \mathbf{k}}\}$ and $\{b_{ij\mathbf{k}}\}$ are given by

$$ \hat{a}_{\tau, \mathbf{k}} = \frac{1}{n} \sum_{i=1}^{n} \phi_{\tau, \mathbf{k}}(\mathbf{X}_i), \quad \text{and} \quad \hat{b}_{ij\mathbf{k}} = \frac{1}{n} \sum_{i=1}^{n} \psi_{i,j,\mathbf{k}}(\mathbf{X}_i), \qquad (9) $$

for any fixed $j_0 \le \tau$, where $\tau = \lambda(n)$ represents the resolution level, a strictly positive integer depending on $n$ and tending to infinity at a specified rate.

By the regularity and compact support conditions on $\phi(\cdot)$ and $\psi(\cdot)$, the above sums are finite for any fixed $\mathbf{x}$, which is crucial for practical implementation. The supports of $\phi(\cdot)$ and $\psi_i(\cdot)$ grow only as a function of their degree of differentiability; see [50]. Next, since $g(\Psi, \cdot) \in L^2(\mathbb{R}^d)$, it admits a wavelet expansion as well. Thus, a linear wavelet estimator for $g(\Psi, \mathbf{x})$ can be written as

$$ \hat{g}_n(\Psi, \mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}'_{\tau, \mathbf{k}} \phi_{\tau, \mathbf{k}}(\mathbf{x}), \qquad (10) $$

or equivalently,

$$ \hat{g}_n(\Psi, \mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}'_{j_0, \mathbf{k}} \phi_{j_0, \mathbf{k}}(\mathbf{x}) + \sum_{j=j_0}^{\tau} \sum_{i=1}^{N} \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{b}'_{ij\mathbf{k}} \psi_{i,j,\mathbf{k}}(\mathbf{x}), \qquad (11) $$

where the coefficient estimators are

$$ \hat{a}'_{\tau, \mathbf{k}} = \frac{1}{n} \sum_{i=1}^{n} \Psi(Y_i) \phi_{\tau, \mathbf{k}}(\mathbf{X}_i), $$

and

$$ \hat{b}'_{ij\mathbf{k}} = \frac{1}{n} \sum_{i=1}^{n} \Psi(Y_i) \psi_{i,j,\mathbf{k}}(\mathbf{X}_i), $$

for any $j_0 \le \tau$. From a practical perspective, the multiresolution properties of wavelets facilitate computationally efficient estimation algorithms that are memory-efficient. In particular, the bandwidth is chosen in the form $2^j$, where $j$ can be selected from a small set of values (e.g., three or four); see [57].
Returning to the wavelet-based estimation, for $\mathbf{u}, \mathbf{v} \in \mathbb{R}^d$, we define the kernel $K(\mathbf{u}, \mathbf{v})$ by

$$K(\mathbf{u}, \mathbf{v}) := \sum_{\mathbf{k} \in \mathbb{Z}^d} \phi(\mathbf{u} - \mathbf{k})\phi(\mathbf{v} - \mathbf{k}). \quad (12)$$

Since

$$|\phi(\mathbf{x})| \le \frac{A_{d+1}}{(1 + \|\mathbf{x}\|)^{d+1}},$$

the kernel function $K(\cdot, \cdot)$ defined in (12) converges uniformly in $\mathbf{u}, \mathbf{v} \in \mathbb{R}^d$. In particular,

$$\sup_{\mathbf{u}, \mathbf{v} \in \mathbb{R}^d} \left| \sum_{\mathbf{k} \in \mathbb{Z}^d} \phi(\mathbf{u} - \mathbf{k})\phi(\mathbf{v} - \mathbf{k}) \right| < \infty.$$

Moreover, for any $j \ge 1$, there exists a constant $C_j$ such that ([49], p. 33)

$$|K(\mathbf{v}, \mathbf{u})| \le \frac{C_j}{(1 + \|\mathbf{v} - \mathbf{u}\|)^j}. \quad (13)$$

From (14), it follows that

$$\int_{\mathbb{R}^d} |K(\mathbf{v}, \mathbf{u})|^j d\mathbf{v} \le G_j(d), \quad (14)$$

where

$$G_j(d) = 2\pi^{d/2} \frac{\Gamma(d)\Gamma(j+d(j-1))}{\Gamma(d/2)\Gamma((d+1)j)} C_{d+1}^j,$$

and $\Gamma(t)$ denotes the Gamma function,

$$\Gamma(t) := \int_{0}^{\infty} y^{t-1}e^{-y} dy.$$

Furthermore, under the assumptions on $\phi(\cdot)$, for any partial derivative of order $|\beta| = 1$, we have ([49], p. 33)

$$\left| \frac{\partial K(\mathbf{u}, \mathbf{v})}{\partial \mathbf{u}_i} \right| \le \frac{C_2}{(1 + \|\mathbf{u} - \mathbf{v}\|)^2}, \quad i = 1, \dots, d. \quad (15)$$

Combining (7), (9), and (12), we see that the linear estimates of the regression function $\hat{m}_n(\varphi, \mathbf{x})$ can be expressed as an extended kernel estimator:

$$K_{h_n, \mathbf{x}}(\mathbf{X}_i) = \frac{1}{h_n^d} K\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right), \quad \text{with } h_n = 2^{-\lambda(n)}.$$

Thus, we have

$$\hat{m}_n(\Psi, \mathbf{x}) = \begin{cases} \frac{\hat{g}_n(\Psi, \mathbf{x})}{\hat{f}_\chi(\mathbf{x})} = \frac{\hat{m}_{2,n}(\Psi, \mathbf{x})}{\hat{m}_{1,n}(\mathbf{x})}, & \text{if } \hat{f}_\chi(\mathbf{x}) \neq 0, \\ \frac{1}{n} \sum_{i=1}^n \varphi(Y_i), & \text{otherwise.} \end{cases}$$

Here,

$$\hat{m}_{1,n}(\Psi, \mathbf{x}) = \frac{\sum_{i=1}^{n} K_{h_n, \mathbf{x}}(\mathbf{X}_i)}{n \mathbb{E}[K_{h_n, \mathbf{x}}(\mathbf{X}_1)]}$$
and

$$ \hat{m}_{2,n}(\Psi, \mathbf{x}) = \frac{\sum_{i=1}^{n} \Psi(Y_i) K_{h_n, \mathbf{x}}(\mathbf{x}_i)}{n \mathbb{E}[K_{h_n, \mathbf{x}}(\mathbf{x}_1)]} $$

This formulation neatly integrates wavelet-based estimation of regression functions with kernel-like estimation procedures.

### 3. Presentation of Estimators

Consider a stationary sequence $\{(\mathbf{X}_i, \mathbf{Y}_i, \mathbf{Z}_i)\}_{i \ge 1}$ of random replicates from the random vector $(\mathbf{X}, \mathbf{Y}, \mathbf{Z})$. Let $f_{\mathbf{X}\mathbf{Y}\mathbf{Z}}(\cdot)$ be the joint density of $(\mathbf{X}, \mathbf{Y}, \mathbf{Z})$ with respect to the Lebesgue measure, and let $f_{\mathbf{X}}(\cdot), f_{\mathbf{Y}}(\cdot)$, and $f_{\mathbf{Z}}(\cdot)$ denote the corresponding marginal densities. For any $\mathbf{x} \in \mathbb{R}^d$, we define the kernel estimator of $f_{\mathbf{X}}(\cdot)$ by

$$ f_n(\mathbf{x}) = \frac{2^{\lambda(n)d}}{n} \sum_{i=1}^{n} K(2^{\lambda(n)}\mathbf{x}, 2^{\lambda(n)}\mathbf{x}_i), $$

where $K(\cdot)$ is a bounded, non-negative kernel function integrating to 1, and $\lambda(n)$ controls the bandwidth growth rate.

For simplicity in what follows, we assume $m_0 = 0$. Observing that model (1) can be recast for a measurable function $\Psi(\cdot) \in \mathcal{F}$ as

$$ \Psi(\mathbf{Y}) - \mathbf{Z}^{\top}\beta = m(\Psi, \mathbf{X}) + \varepsilon, \qquad (16) $$

we propose an internal estimator of the nonparametric component in (16). For any $\mathbf{x} \in \mathbb{R}^d$, this estimator is given by

$$ \hat{m}_{n}^{\beta}(\Psi, \mathbf{x}, h) = \sum_{i=1}^{n} \frac{2^{\lambda(n)d} (\Psi(\mathbf{Y}_i) - \mathbf{Z}_i^{\top} \beta)}{n f_n(\mathbf{X}_i)} K(2^{\lambda(n)}\mathbf{x}, 2^{\lambda(n)}\mathbf{x}_i), \qquad (17) $$

where $h = 2^{-\lambda(n)}$ is a smoothing parameter tending to zero as specified by Condition (H.2).

Since $\hat{m}_{n}^{\beta}(\cdot)$ depends on the unknown $\beta$, we recognize that the true function $m(\cdot)$ in (16) may also naturally depend on $\beta$. In particular, when $m(\cdot)$ has an additive form, we can write

$$ m_{add}^{\beta}(\Psi, \mathbf{x}) = \sum_{\ell=1}^{d} m_{\ell}^{\beta}(\Psi, x_{\ell}), $$

where $x_{\ell}$ is the $\ell$-th component of $\mathbf{x}$. To ensure identifiability of each additive component $m_{\ell}(\cdot)$, we impose

$$ \mathbb{E}[m_{\ell}^{\beta}(\Psi, X_{\ell})] = 0 \quad \text{for} \quad 1 \le \ell \le d. $$

There are various approaches to estimating the additive components in regression models. Two notable examples are the marginal integration technique (see, for instance, [100–104]) and back-fitting algorithms (cf. [105,106], and references therein). While back-fitting poses more intricate asymptotic challenges (though progress is reported in [107,108]), marginal integration methods tend to exhibit more transparent asymptotic properties, making them particularly appealing in statistical and econometric contexts.

Guided by these advantages, we adopt a marginal integration framework to estimate the linear parameter $\beta$ through a least squares criterion. Crucially, our estimators require only kernel functions, smoothing parameters, and the density estimators, avoiding the need for complex optimization techniques or weighting schemes that depend on unknown quantities (as is sometimes required in other methods; see, e.g., [109]).
### 3.1. The Marginal Integration Method

In this work, we rely on marginal integration to build our estimators. To outline this approach, let us first establish some notation. For $1 \le \ell \le d$, let

$$ \mathbf{x}_{-\ell} = (x_1, \dots, x_{\ell-1}, x_{\ell+1}, \dots, x_d). $$

We then define the functions

$$ q_{-\ell}(\mathbf{x}_{-\ell}) = \prod_{\substack{j=1 \\ j \ne \ell}}^{d} q_j(x_j) \quad \text{and} \quad q(\mathbf{x}) = \prod_{\ell=1}^{d} q_{\ell}(x_{\ell}), $$

where $q_{\ell}$ is a known univariate density for each $\ell \in \{1, \dots, d\}$. Throughout, all integrals over continuous variables are assumed to be taken with respect to the Lebesgue measure. By applying marginal integration, the additive regression function can be represented as

$$ m_{add}^{\beta}(\Psi, \mathbf{x}) = \sum_{\ell=1}^{d} \eta_{\ell}^{\beta}(\Psi, x_{\ell}) + \int_{\mathbb{R}^d} m_{add}^{\beta}(\Psi, \mathbf{z}) q(\mathbf{z}) d\mathbf{z}, \quad (18) $$

where each component $\eta_{\ell}^{\beta}$ is given by

$$ \eta_{\ell}^{\beta}(\Psi, x_{\ell}) = \int_{\mathbb{R}^{d-1}} m_{add}^{\beta}(\Psi, \mathbf{x}) q_{-\ell}(\mathbf{x}_{-\ell}) d\mathbf{x}_{-\ell} - \int_{\mathbb{R}^d} m_{add}^{\beta}(\Psi, \mathbf{x}) q(\mathbf{x}) d\mathbf{x}. \quad (19) $$

For notational convenience, we will often write $m_{add}(\cdot)$ instead of $m_{add}^{\beta}(\cdot)$ and $\eta_{\ell}(\cdot)$ instead of $\eta_{\ell}^{\beta}(\cdot)$.

Given the parameter $\beta$, we define an estimator of $m_{add}(\cdot)$ via marginal integration as follows:

$$ \hat{m}_{add}^{\beta}(\Psi, \mathbf{x}) = \sum_{\ell=1}^{d} \hat{\eta}_{\ell}^{\beta}(\Psi, x_{\ell}) + \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\Psi, \mathbf{z}) q(\mathbf{z}) d\mathbf{z}, \quad (20) $$

where

$$ \hat{\eta}_{\ell}^{\beta}(\Psi, x_{\ell}) = \int_{\mathbb{R}^{d-1}} \hat{m}_{n}^{\beta}(\Psi, \mathbf{x}) q_{-\ell}(\mathbf{x}_{-\ell}) d\mathbf{x}_{-\ell} - \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\Psi, \mathbf{x}) q(\mathbf{x}) d\mathbf{x}. \quad (21) $$

Throughout the remainder of the paper, we will simplify the notation by writing $\hat{m}_{add}(\cdot)$ for $\hat{m}_{add}^{\beta}(\cdot)$, $\hat{m}_{n}(\cdot)$ for $\hat{m}_{n}^{\beta}(\cdot)$, and $\hat{\eta}_{\ell}(\cdot)$ for $\hat{\eta}_{\ell}^{\beta}(\cdot)$. Consequently, once $\beta$ has been estimated, we obtain a fully specified estimator $\hat{m}_{add}(\cdot)$ for the additive regression function.

### 3.2. Estimation of the Model's Parameter

Consider the following partially linear additive regression model:

$$ \Psi(\mathbf{Y}) = \mathbf{Z}^{\top}\beta + m_{add}(\mathbf{X}) + \varepsilon. $$

Throughout the subsequent discussion, for $\ell = 1, \dots, d$, we define

$$ K_{\ell}(u_{\ell}) = \int_{\mathbb{R}^{d-1}} K(\mathbf{u}) d\mathbf{u}_{-\ell} \quad \text{and} \quad K_{-\ell}(\mathbf{u}_{-\ell}) = \int_{\mathbb{R}} K(\mathbf{u}) du_{\ell}. $$

Following the methodology outlined in [42], we estimate $\beta$ using

$$ \hat{\beta} = [\tilde{Z}\tilde{Z}^{\top}]^{-1}\tilde{Z}\tilde{Y}, \quad (22) $$
where

$$ \tilde{Y} = \left[ \psi(\mathbf{Y}_i) - \sum_{j=1}^{n} W_{nj}(\mathbf{X}_i) \Psi(\mathbf{Y}_j) \right]_{1 \le i \le n}^{\top}, \qquad (23) $$

$$ \tilde{Z} = \left[ Z_i - \sum_{j=1}^{n} W_{nj}(X_i) Z_j \right]_{1 \le i \le n}, \qquad (24) $$

$$ W_{nj}(\mathbf{X}_i) = \frac{U_{nj}(\mathbf{X}_i)}{nf_n(\mathbf{X}_j)}, $$

and

$$ U_{nj}(\mathbf{X}_i) = \sum_{\ell=1}^{d} 2^{\lambda(n)} K_{\ell}(2^{\lambda(n)} X_{\ell,i}, 2^{\lambda(n)} X_{j,\ell}) D_{-\ell} \\ -(d-1) \int_{\mathbb{R}^d} 2^{\lambda(n)d} K(2^{\lambda(n)} \mathbf{z}, 2^{\lambda(n)} \mathbf{X}_j) q(\mathbf{z}) d\mathbf{z}. $$

Here,

$$ D_{-\ell} = \int_{\mathbb{R}^{d-1}} 2^{\lambda(n)(d-1)} K(2^{\lambda(n)} \mathbf{z}_{-\ell}, 2^{\lambda(n)} \mathbf{X}_{-\ell,j}) q_{-\ell}(\mathbf{z}_{-\ell}) d\mathbf{z}_{-\ell}. $$

For additional information regarding the estimation of $\beta$, we refer the reader to [16]. Finally, to estimate the regression function and the additive components, we use

$$ \hat{m}_{add}(\Psi, \mathbf{x}, h) = \sum_{\ell=1}^{d} \hat{\eta}_{\ell}(\mathbf{x}_{\ell}) + \int_{\mathbb{R}^d} \hat{m}_n(\mathbf{z}) q(\mathbf{z}) d\mathbf{z}, \qquad (25) $$

and

$$ \hat{\eta}_{\ell}(\Psi, x_{\ell}, h) = \int_{\mathbb{R}^{d-1}} \hat{m}_n(\mathbf{x}) q_{-\ell}(\mathbf{x}_{-\ell}) d\mathbf{x}_{-\ell} - \int_{\mathbb{R}^d} \hat{m}_n(\mathbf{x}) q(\mathbf{x}) d\mathbf{x}. \qquad (26) $$

## 4. Main Results

To present our findings, we introduce an additional structural assumption on the model. Specifically, we assume that

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} \tilde{Z}_{i}^{\top} \tilde{Z}_{i} = B \quad \text{a.s.,} \qquad (27) $$

where $B$ is a $p \times p$ positive-definite matrix. Let

$$ I^d = \prod_{\ell=1}^{d} I_{\ell} = \prod_{\ell=1}^{d} [a_{\ell}, c_{\ell}], $$

and

$$ J^d = \prod_{\ell=1}^{d} J_{\ell} = \prod_{\ell=1}^{d} [a'_{\ell}, c'_{\ell}], $$

be two fixed partitions of $\mathbb{R}^d$ such that

$$ a'_{\ell} < a_{\ell} < c_{\ell} < c'_{\ell}, \quad \text{for all } 1 \le \ell \le d. $$

Additional assumptions regarding the density of $\mathbf{X}$, the regression function $m(\cdot)$, the kernel functions $K(\cdot)$, $K(\cdot)$, and the smoothing parameters are collected below for ease of reference
and the reader's convenience. The first set of conditions concerns the density function $f_X(\cdot)$ and the regression function $m(\cdot)$.

(G.1) The density $f_X(\cdot)$ is continuous and bounded away from 0 on $\mathbb{R}^d$.

(G.2) One of the following holds:

(i) $f_X(\cdot)$ is bounded and uniformly continuous on $\mathbb{R}^d$, or

(ii) $f_X(\cdot)$ belongs to a Besov space $B_{s,p,q}$ for some $0 < s < r, 1 \le p, q \le \infty$.

(G.3) For all $\mathbf{x}, \mathbf{y}$, we have

$$|f_j(\mathbf{x}, \mathbf{y}) - f_X(\mathbf{x})f_X(\mathbf{y})| \le M < \infty \quad \text{for all } j \ge 1,$$

where $f_j(\mathbf{x}, \mathbf{y})$ is the joint probability density function of the vectors $(\mathbf{X}_0, \mathbf{X}_j)$, assumed to exist.

(G.4) (i) For all $1 \le \ell \le d$, the functions $m_{\ell}^{\beta}(\Psi, \cdot)$ are Lipschitz continuous of order $(s - \frac{1}{p})$.

(ii) $\sup_{\mathbf{x} \in I^d} g(\varphi, \mathbf{x}) < \infty$.

(iii) The function $H(\cdot)$ belongs to the Besov space $B_{s,p,q}$ for some $0 < s < r, 1 \le p, q \le \infty$.

(G.5) The process $\{\mathbf{X}_i, \mathbf{Y}_i, \mathbf{Z}_i\}$ is $\alpha$-mixing with a mixing coefficient $\alpha(s)$ satisfying

$$\alpha(s) = \mathcal{O}(s^{-\nu}), \text{ where } \nu > 2d + 3.$$

In what follows, let $\lambda(n)$ denote a sequence of positive constants satisfying the following conditions:

$$\text{(H.1) As } n \to \infty, \qquad n^{1/2\lambda(n)} \to \infty, \quad 2^{\lambda(n)} \to \infty, \quad n2^{-\frac{\lambda(n)}{a}} \to 0, \\ \text{and} \\ \qquad (n2^{-\lambda(n)})^{1/2} (n^{-1} \log n)^{\frac{s}{d+2s}} \to 0.$$

$$\text{(H.2) As } n \to \infty, \qquad n^{-1}2^{\lambda(n)} \log n \log_2(n) \to 0.$$

$$\text{(H.3) As } n \to \infty, \qquad \sqrt{\frac{2^{3\lambda(n)}}{n}} \left(\frac{n}{\log n}\right)^{a\left(s-\frac{1}{p}\right)} \to 0.$$

Next, we state assumptions on the random variables $\mathbf{Z}$ and $\mathbf{Y}$.

(M.1) $\Psi(\mathbf{Y})$ and $\mathbf{Z}$ are bounded.

To present our main results, we impose conditions on the density functions $q_\ell(\cdot)$ for $1 \le \ell \le d$.

(Q.1) Each $q_\ell(\cdot)$ is bounded and continuous for all $1 \le \ell \le d$. Moreover, $q_\ell(\cdot)$ is bounded and uniformly continuous on $\mathbb{R}$ for all $1 \le \ell \le d$.

Alternatively, we assume the following:

(Q.2) The function $q_\ell(\cdot)$ belongs to a Besov space $B_{s,p,q}$ for some $0 < s < r, 1 \le p, q \le \infty$, for all $1 \le \ell \le d$.

(R.1) The multiresolution analysis is $r$-regular.

Our principal finding is encapsulated in the following result.
**Theorem 1.** Suppose that all the conditions (G.1)–(G.5), (H.1)–(H.3), (M.1), (Q.1), (Q.2), and (R.1) are satisfied. Then, for a given measurable function $\Psi(\cdot) \in \mathcal{F}$, it holds that

$$ \sup_{x \in \mathcal{C}} |\hat{m}_{\text{add}}^{\beta}(\Psi, x) - m(x)| = \mathcal{O}\left(\left[\frac{\log n}{n}\right]^{a\left(s-\frac{1}{p}\right)}\right). \quad (28) $$

**Remark 1.** It is noteworthy that the quantity $2^{-\lambda(n)}$ appearing in the preceding theorems corresponds directly to the bandwidth $h_n$ utilized in Parzen–Rosenblatt convolution kernel density estimators. From a practical perspective, selecting the multiresolution level $\tau_n$ associated with wavelets is often more straightforward than choosing the bandwidth $h_n$ itself. Indeed, only a relatively small number of candidate values for $\lambda(n)$ (e.g., three or four) need to be considered in practice, thus rendering the procedure more accessible and computationally feasible.

**Remark 2.** It is important to emphasize that the precise logarithmic convergence rates are influenced by both the chosen resolution level and the underlying smoothness parameters $s$ associated with $f(\cdot)$, as posited in (G.2). The function $f(\cdot)$ lies within a Besov-type space $B_{s,p,q}$, which allows for more general and subtle smoothness characterizations than classical integer-order differentiability conditions. Such complexity in smoothness assumptions is well-known in nonparametric estimation and aligns with established results in the literature. While $s$ may be unknown in practice, these convergence rates are standard and commonly observed across nonparametric procedures. To render the estimator fully data-driven, various adaptive techniques have been proposed for selecting the resolution level $\tau$, including methods inspired by Stein's unbiased risk estimation, heuristic "rules of thumb," and cross-validation approaches. For detailed discussions on the design and analysis of asymptotically optimal bandwidth selection rules, the interested reader may refer to [53,57].

**Remark 3.** The moment condition (M.1) can be generalized substantially. Instead of requiring boundedness conditions, one may impose more general assumptions on the moments of $\mathbf{Y}$, as outlined in [110]. Specifically, condition (M.1) can be replaced by one of the following stronger but more flexible alternatives:

(M.1)' There exists a measurable envelope function $F: \mathbb{R}^q \to [0, \infty)$ such that for all $\mathbf{y} \in \mathbb{R}^q$,

$$ F(\mathbf{y}) \geq \sup_{\Psi \in \mathcal{F}} |\Psi(\mathbf{y})|, $$

and for some $s > 2$,

$$ \sup_{\mathbf{x} \in I^d} \mathbb{E}(|F(\mathbf{Y})|^s | \mathbf{X} = \mathbf{x}) < \infty. $$

(M.1)'' Consider a nonnegative, continuous, and nondecreasing function $\{\mathcal{M}(x) : x \geq 0\}$ defined on $[0, \infty)$. Assume that for some $s > 2$, as $x \to \infty$:

$$ (i) \quad x^{-s} \mathcal{M}(x) \text{ is eventually decreasing, and } (ii) \quad x^{-1} \log \mathcal{M}(x) \text{ is eventually increasing.} \quad (29) $$

For any $t \geq \mathcal{M}(0)$, define $\mathcal{M}^{\text{inv}}(t) \geq 0$ by

$$ \mathcal{M}(\mathcal{M}^{\text{inv}}(t)) = t. $$

We further require

$$ \sup_{\mathbf{x} \in I^d} \mathbb{E}(\mathcal{M}(|F(\mathbf{Y})|) | \mathbf{X} = \mathbf{x}) < \infty. $$

Notable examples of functions $\mathcal{M}(\cdot)$ that satisfy these criteria include the following:

1. $\mathcal{M}(x) = x^p$ for some $p > 2$;

2. $\mathcal{M}(x) = \exp(sx)$ for some $s > 0$.
In other words, the assumption (M1) may be replaced by weaker, but more intricate, finite-moment-type conditions. This change, however, significantly complicates the technical details of the proofs. For further insights and a more comprehensive discussion, the reader may consult [111]. Finally, the introduction of the function $\Psi(\cdot)$ into our framework is motivated by findings and observations stated in Remark 1.2 of [111].

### 4.1. The Conditional Distribution Function

By incorporating $\Psi(\mathbf{Y}) = 1\{\mathbf{Y} \le \mathbf{t}\}$, for $\mathbf{t} \in \mathbb{R}^q$, into (17) under the assumption $\beta = 0$, we obtain the following kernel estimator of the conditional distribution function

$$F(\mathbf{t} | \mathbf{x}) := \mathbb{P}(\mathbf{Y} \le \mathbf{t} | \mathbf{X} = \mathbf{x}),$$

namely,

$$\tilde{m}_n(\mathbf{t} | \mathbf{x}) := \sum_{i=1}^{n} \frac{2^{\lambda(n)d} \mathbb{1}(\mathbf{Y}_i \le \mathbf{t})}{n f_n(\mathbf{X}_i)} K(2^{\lambda(n)}\mathbf{x}, 2^{\lambda(n)}\mathbf{X}_i). \quad (30)$$

ref. [112] refers to this as the conditional empirical distribution function and was the first to establish uniform consistency results for it. Next, consider the following additive structure:

$$F_{add}(\mathbf{t} | \mathbf{x}) = \sum_{\ell=1}^{d} m_{\ell}(\mathbf{t}, \mathbf{x}_{\ell}).$$

Following the marginal integration procedure, for any $\mathbf{x} \in \mathbb{R}^d$, the corresponding additive regression function takes the form

$$F_{add}(\mathbf{t} | \mathbf{x}) = \sum_{\ell=1}^{d} \eta_{\ell}(\mathbf{t}, \mathbf{x}_{\ell}) + \int_{\mathbb{R}^d} F_{add}(\mathbf{t} | \mathbf{z}) q(\mathbf{z}) d\mathbf{z},$$

where

$$\eta_{\ell}(\mathbf{t}, \mathbf{x}_{\ell}) = \int_{\mathbb{R}^{d-1}} F_{add}(\mathbf{t} | \mathbf{x}) q_{-\ell}(\mathbf{x}_{-\ell}) d\mathbf{x}_{-\ell} - \int_{\mathbb{R}^d} F_{add}(\mathbf{t} | \mathbf{x}) q(\mathbf{x}) d\mathbf{x}. \quad (31)$$

To estimate both the conditional empirical distribution and the additive components, we define

$$\tilde{F}_{add}(\mathbf{t}, \mathbf{x}, h_n) = \sum_{\ell=1}^{d} \tilde{\eta}_{\ell}(\mathbf{t}, \mathbf{x}_{\ell}, h_n) + \int_{\mathbb{R}^d} \tilde{m}_n(\mathbf{t} | \mathbf{z}) q(\mathbf{z}) d\mathbf{z}, \quad (32)$$

and

$$\tilde{\eta}_{\ell}(\mathbf{t}, \mathbf{x}_{\ell}, h_n) = \int_{\mathbb{R}^{d-1}} \tilde{m}_n(\mathbf{t} | \mathbf{x}) q_{-\ell}(\mathbf{x}_{-\ell}) d\mathbf{x}_{-\ell} - \int_{\mathbb{R}^d} \tilde{m}_n(\mathbf{t} | \mathbf{x}) q(\mathbf{x}) d\mathbf{x}. \quad (33)$$

By applying Theorem 1, we obtain the uniform convergence rate

$$\sup_{\mathbf{x} \in \mathcal{C}} |\tilde{F}_{add}(\mathbf{t}, \mathbf{x}, h) - F_{add}(\mathbf{t} | \mathbf{x})| = \mathcal{O}\left(\left[\frac{\log n}{n}\right]^{a\left(s-\frac{1}{p}\right)}\right). \quad (34)$$

### 4.2. Relative-Error Prediction

Classical regression typically estimates the operator $r(\cdot)$ by minimizing the conditional expectation of the squared error, $\mathbb{E}[(Y - r(X))^2 | X]$. While widely used, this approach imposes equal emphasis on all observations, rendering it sensitive to outliers and potentially unsuitable in situations involving large variability in the response values.
To address these drawbacks, we focus on minimizing the mean squared relative error (MSRE) defined by

$$ \mathbb{E}\left[\left(\frac{Y - r(X)}{Y}\right)^2 \mid X\right], \quad \text{with } Y > 0 \text{ almost surely.} $$

This criterion offers a more robust assessment of predictive performance when the range of predicted values is extensive. Moreover, the unique feature of the MSRE-based approach is that its solution can be expressed as a ratio of the first two conditional inverse moments of Y. Specifically, letting

$$ g_{\gamma}(x) = \mathbb{E}[Y^{-\gamma} \mid X = x], \quad \gamma = 1, 2, $$

and assuming these moments exist and are finite almost surely, one can show, following [113–115], that the optimal MSRE predictor of Y given X is

$$ \check{r}(x) = \frac{\mathbb{E}[Y^{-1} \mid X = x]}{\mathbb{E}[Y^{-2} \mid X = x]}, \quad \text{a.s.} $$

By incorporating $\Psi(Y) = Y^{-\gamma}$, into (17) under the assumption $\beta = 0$, we obtain the following kernel estimator

$$ \check{r}_n(x) = \frac{\hat{m}_{\text{add}}^0(y^{-1}, \mathbf{x})}{\hat{m}_{\text{add}}^0(y^{-2}, \mathbf{x})} $$

By leveraging Theorem 1 for the particular cases $\Psi(y) = y^{-1}$ and $\Psi(y) = y^{-2}$, we establish the following uniform convergence rate:

$$ \sup_{\mathbf{x} \in C} |\check{r}_n(\mathbf{x}) - \check{r}(\mathbf{x})| = \mathcal{O}\left(\left[\frac{\log n}{n}\right]^{a\left(s - \frac{1}{p}\right)}\right). $$

This result not only corroborates previous findings but also extends and enhances the analyses presented in [116].

**Remark 4.** It has long been recognized that kernel estimators become less reliable when the dimensionality of the data increases. This phenomenon, commonly known as the curse of dimensionality [117], emerges because, in high-dimensional settings, local regions typically require exceedingly large sample sizes to accumulate enough observations. As a result, unless sample sizes are enormous, bandwidths must be made so broad that the notion of “local” averaging is compromised. A thorough exploration of these issues, with numerical illustrations, can be found in [118], and more recent references include [119,120]. Despite their popularity, however, there remains significant uncertainty regarding the asymptotic properties of penalized splines, even in standard nonparametric settings, and only a few theoretical investigations address their performance. In parallel, it should be noted that most functional regression techniques rely on minimizing an $L_2$-norm and are thus sensitive to outliers. Other less conventional but potentially valuable tools include methods using delta sequences and wavelet-based approaches.

**Remark 5.** In statistical learning for time series, it is often assumed that the underlying data-generating process exhibits some form of asymptotic independence, commonly referred to as “mixing”. While many standard processes are presumed to follow known types of mixing rates, the precise coefficients that govern these rates generally remain unspecified. Moreover, such mixing assumptions are seldom examined empirically, and there are currently no established methods to
estimate mixing coefficients directly from observed data. To address this gap, ref. [121] introduces an estimator for beta-mixing coefficients based on a single stationary realization. Subsequently, ref. [122] proposes strongly consistent estimators for the $ℓ_1$-norm of both $\alpha$- and $\beta$-mixing coefficients in a stationary ergodic process. These estimators then serve as the foundation for strongly consistent hypothesis tests for goodness-of-fit. In particular, ref. [122] develops tests to determine whether the $\alpha$-mixing or $\beta$-mixing coefficients of a process remain bounded above by a specified rate function, under the same summability assumptions.

**Remark 6.** Nonlinear thresholding approaches yield the following estimators for $g(\Psi, \mathbf{x})$ and $f_X(\mathbf{x})$. Specifically, we define

$$ \hat{f}_{\mathbf{X}}(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}_{j_0, \mathbf{k}} \phi_{j_0, \mathbf{k}}(\mathbf{x}) + \sum_{j=j_0}^{\tau} \sum_{i=1}^{N} \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{b}_{ijk} \mathbf{1}\{ |\hat{b}_{ijk}| > \delta_{j,n} \} \psi_{i,j,\mathbf{k}}(\mathbf{x}), \quad (35) $$

and

$$ \hat{g}_n(\Psi, \mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}'_{j_0, \mathbf{k}} \phi_{j_0, \mathbf{k}}(\mathbf{x}) + \sum_{j=j_0}^{\tau} \sum_{i=1}^{N} \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{b}'_{ijk} \mathbf{1}\{ |\hat{b}'_{ijk}| > \delta_{j,n} \} \psi_{i,j,\mathbf{k}}(\mathbf{x}), \quad (36) $$

where $\delta_{j,n}$ is an appropriately selected threshold. In the one-dimensional setting ($d = 1$), the estimator $\hat{f}_{\mathbf{X}}$ was originally proposed by [56]. Exploring the applicability and performance of these estimators for higher-dimensional cases would be a valuable direction for future research.

**Remark 7.** Numerous investigations, such as [54–56], report that thresholded nonlinear wavelet-based methods frequently outperform linear alternatives when assessed by mean integrated error. However, Theorem 3 and its corollary reveal that, for any fixed function $g$ in the Besov space $B_{s,p,q}$, linear wavelet estimators already attain the classic optimal convergence rates under the almost sure sup-norm criterion in the standard i.i.d. setting. This implies that, in terms of convergence rates, the nonlinear thresholded estimators in (35)–(36) confer no additional advantage over the linear estimators in (8)–(11) once the almost sure sup-norm is adopted as the performance measure.

This result naturally raises the question of why the sup-norm metric is preferable to the mean integrated error. Many functions in the broad Besov space $B_{s,p,q}$ exhibit significant local oscillations, which can cause large estimation errors in small subsets of $\mathbb{R}^d$. Such localized errors have a negligible impact on mean integrated metrics yet are readily detected by the sup-norm. This issue is especially acute for discontinuous functions, as substantial misestimation around a discontinuity may remain hidden under mean integrated error. Empirical findings by [55] corroborate this, illustrating how pointwise bias and variance in wavelet-based estimators can oscillate considerably while appearing largely diminished in the mean integrated squared error.

Moreover, in most practical scenarios, we typically observe only one realization of the underlying process. An almost sure sup-norm analysis reveals how estimators perform for almost every realization, in contrast to metrics such as mean-square error or mean integrated squared error, which capture only average behavior across multiple realizations. Therefore, we propose that the almost sure sup-norm constitutes a more fitting benchmark for density and regression estimation in $B_{s,p,q}$. Additional insights on this perspective can be found in the introduction of [71].

**Remark 8.** Our approach, which integrates marginal integration with linear wavelet-based techniques, facilitates the effective separation of linear influences from complex nonlinear patterns. This capability is especially valuable in managerial settings, where comprehending both direct and subtle impacts on outcomes is crucial. For instance:

1. **Finance:** Managers can utilize this model to enhance risk assessment and optimize portfolio strategies. By capturing nonlinear dependencies among market variables, the model yields more
precise forecasts and deeper insights into asset behavior, thereby supporting more informed investment decisions.

2. Biology: In research and development, the ability to detect and quantify nonlinear effects assists in guiding resource allocation and experimental design, ultimately leading to more efficient innovation and product development strategies.

3. Engineering: In engineering contexts, our method's robust estimation of critical performance indicators—despite the presence of complex nonlinear interactions—substantially enhances operational control and system optimization. This increased precision enables managers to implement quality control measures that are both accurate and reliable, ensuring that engineering systems consistently perform within established standards.

**5. Concluding Remarks**

In this work, we establish strong uniform convergence results for estimators of nonlinear additive components derived from integrating marginal methods with linear wavelet-based estimation strategies. These findings significantly advance the theoretical underpinnings of additive modeling under mixing conditions, which characterizes the dependence structure of the underlying random variable sequences. In particular, wavelet estimators excel in detecting discontinuities, offering a distinct advantage over kernel-based techniques. This benefit is most pronounced when target functions exhibit spatial irregularities, especially for $1 \le p < 2$, or when function smoothness is not specified. In such scenarios, purely linear smoothing strategies cannot attain the minimax optimal convergence rate, making nonlinear wavelet estimators (e.g., soft and hard thresholding) indispensable. Nevertheless, constructing and analyzing these nonlinear approaches introduces considerable mathematical complexity, which remains an important avenue for future research to broaden the applicability and bolster the performance of wavelet-based estimators.

A promising direction for further exploration involves the application of wavelet techniques to conditional *U*-statistics. Integrating wavelet-based methods into this setting may open up new theoretical insights and spur the development of enhanced inference procedures. Advancing this line of research is therefore highly encouraged as it has the potential to substantially enrich methodological progress in the analysis of complex stochastic processes.

The advent of high-throughput genomic technologies—covering genome-wide single nucleotide polymorphisms and RNA sequencing—has generated extensive datasets that require robust statistical models to quantify and interpret how genetic variation and gene expression influence quantitative traits [123]. In contrast to conventional data structures, where the number of variables is typically much smaller than the sample size, modern high-dimensional data often have as many or more variables ($p$) than observations ($n$). This imbalance poses serious challenges to traditional statistical methods, which generally assume $n \gg p$.

Consequently, it is critical to examine the viability of partial linear additive models (PLAM) in high-dimensional frameworks, particularly where both linear and nonlinear components can exceed the sample size, as emphasized in [123]. Such investigations promise to refine our analytical toolkit for complex genomic data and deepen our understanding of the genetic mechanisms that underlie quantitative traits.

Finally, to reinforce the practical significance of this work, comprehensive simulation studies and applications to real datasets would be highly advantageous. We also intend to perform a comprehensive sensitivity analysis to assess the robustness of our findings under various conditions. Such empirical investigations will offer more concrete guidelines for the effective use of these methods, which will underscore their potential to advance contemporary statistical analysis.
## 6. Proofs

This section is dedicated to proving our main result. Throughout this proof, we will continue to use the notation introduced earlier. To unburden our notation a bit, $\hat{m}_{add}^{\hat{\beta}}(\Psi, \cdot)$ will be denoted by $\hat{m}_{add}^{\hat{\beta}}(\cdot)$. We begin by establishing the result under the assumption that the marginal density $f(\cdot)$ is known. Subsequently, we provide a brief outline of the proof in which the marginal density is estimated.

**Proof of Theorem 1.** Since $C$ is compact, we can cover it with a finite collection of balls $B_p(t_p, r_p)$ such that

$$C \subset \bigcup_{p=1}^{r(n)} B_p,$$

where $r(n)$ will be specified later. Therefore, for any point $x \in C$, there exists a ball centered at $t(x)$ that contains $x$. Consequently, the approximation error can be decomposed as follows:

$$ \begin{aligned} & \sup_{x \in C} |\hat{m}_{add}^{\hat{\beta}}(\mathbf{x}) - m(\mathbf{x})| \\ & \le \sup_{x \in C} |\hat{m}_{add}^{\hat{\beta}}(\mathbf{x}) - \hat{m}_{add}^{\hat{\beta}}(\mathbf{x})| + \sup_{x \in C} |\hat{m}_{add}^{\hat{\beta}}(\mathbf{x}) - \hat{m}_{add}^{\hat{\beta}}(t(\mathbf{x}))| \\ & \quad + \sup_{x \in C} |\hat{m}_{add}^{\hat{\beta}}(t(\mathbf{x})) - \mathbb{E}[\hat{m}_{add}^{\hat{\beta}}(t(\mathbf{x}))]| + \sup_{x \in C} |\mathbb{E}[\hat{m}_{add}^{\hat{\beta}}(t(\mathbf{x}))] - \mathbb{E}[\hat{m}_{add}^{\hat{\beta}}(\mathbf{x})]| \\ & \quad + \sup_{x \in C} |\mathbb{E}[\hat{m}_{add}^{\hat{\beta}}(\mathbf{x})] - m(\mathbf{x})|. \end{aligned} \quad (37) $$

This decomposition shows how each term contributes to the overall approximation error. Recall that

$$\hat{m}_{add}^{\hat{\beta}}(\mathbf{x}) = \sum_{\ell=1}^{d} \hat{\eta}_{\ell}^{\beta}(x_{\ell}) + \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) d\mathbf{z}.$$

Hence,

$$ \begin{aligned} & |\mathbb{E}[\hat{m}_{add}^{\hat{\beta}}(\mathbf{x})] - m(\mathbf{x})| \\ &= \left| \mathbb{E} \left[ \sum_{\ell=1}^{d} \hat{\eta}_{\ell}^{\beta}(x_{\ell}) + \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) d\mathbf{z} \right] - \sum_{\ell=1}^{d} \eta_{\ell}(x_{\ell}) - \int_{\mathbb{R}^d} m(\mathbf{z}) q(\mathbf{z}) d\mathbf{z} \right| \\ &= \left| \sum_{\ell=1}^{d} \left( \mathbb{E}[\hat{\eta}_{\ell}^{\beta}(x_{\ell})] - \eta_{\ell}(x_{\ell}) \right) + \int_{\mathbb{R}^d} \left( \mathbb{E}[\hat{m}_{n}^{\beta}(\mathbf{z})] - m(\mathbf{z}) \right) q(\mathbf{z}) d\mathbf{z} \right| \\ &=: |\mathcal{B}_1 + \mathcal{B}_2|. \end{aligned} \quad (38) $$

Here, $\mathcal{B}_1$ captures the main bias for the additive components $\hat{\eta}_{\ell}^{\beta}(\cdot)$, while $\mathcal{B}_2$ corresponds to the integral term involving $\hat{m}_{n}^{\beta}(\mathbf{z})$. Next, recall that

$$\hat{\eta}_{\ell}^{\beta}(x_{\ell}) = \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) - \int_{\mathbb{R}} \hat{\zeta}_{n,\ell}^{\beta}(u_{\ell}) q_{\ell}(u_{\ell}) du_{\ell},$$

where

$$\hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) = \frac{2^{\lambda(n)}}{n} \sum_{i=1}^{n} \frac{\tilde{Y}_{i,n}^{\beta}}{f_{\ell}(X_{\ell,i})} K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}X_{\ell,i}),$$
and

$$ \tilde{Y}_{i,n}^{\beta} = (\Psi(\mathbf{Y}_i) - \mathbf{Z}_i^{\top} \boldsymbol{\beta}) \int_{\mathbb{R}^{d-1}} 2^{(d-1)\lambda(n)} K_{-\ell}(2^{\lambda(n)}\mathbf{x}_{-\ell}, 2^{\lambda(n)}\mathbf{x}_{-\ell,i}) \frac{q_{-\ell}(\mathbf{x}_{-\ell})}{f(\mathbf{x}_{-\ell,i} | \mathbf{X}_{\ell,i})} d\mathbf{x}_{-\ell}. $$

The term $\hat{\zeta}_{n,\ell}^{\beta}(\cdot)$ can be viewed as a univariate regression estimator of

$$ \tilde{m}_n(\cdot) = \mathbb{E}[\tilde{Y}_{i,n}^{\beta} | \mathbf{X}_{\ell,i} = \cdot]. $$

In accordance with the approach of [103], let

$$ C_n := \mu + \int_{\mathbb{R}^{d-1}} \sum_{j=1, j \neq \ell}^{d} m_j(z_j) g_n(\mathbf{z}_{-\ell}) d\mathbf{z}_{-\ell}, $$

where

$$ g_n(\mathbf{z}_{-\ell}) := \int_{\mathbb{R}^{d-1}} 2^{(d-1)\lambda(n)} K_{-\ell}(2^{\lambda(n)}\mathbf{u}_{-\ell}, 2^{\lambda(n)}\mathbf{z}_{-\ell}) q_{-\ell}(\mathbf{u}_{-\ell}) d\mathbf{u}_{-\ell}. $$

Then

$$ \tilde{m}_n(x_\ell) = m_\ell(x_\ell) + C_n. \quad (39) $$

Furthermore,

$$ \eta_{\ell}(x_{\ell}) = m_{\ell}(x_{\ell}) - \int_{\mathbb{R}} m_{\ell}(u_{\ell}) q_{\ell}(u_{\ell}) du_{\ell}. \quad (40) $$

Combining (39) and (40) leads to

$$ \eta_{\ell}(x_{\ell}) = \tilde{m}_{n}(x_{\ell}) - C_{n} - \int_{\mathbb{R}} m_{\ell}(u_{\ell}) q_{\ell}(u_{\ell}) du_{\ell}. \quad (41) $$

Additionally, for $1 \le \ell \le d$, we know

$$ \hat{\eta}_{\ell}^{\beta}(x_{\ell}) = \int_{\mathbb{R}^{d-1}} \hat{m}_{n}^{\Psi, \beta}(\mathbf{x}) q_{-\ell}(\mathbf{x}_{-\ell}) d\mathbf{x}_{-\ell} - \int_{\mathbb{R}^{d}} \hat{m}_{n}^{\Psi, \beta}(\mathbf{x}) q(\mathbf{x}) d\mathbf{x}, $$

which implies

$$ \hat{\eta}_{\ell}^{\beta}(x_{\ell}) = \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) - \int_{\mathbb{R}^d} \hat{m}_{n}^{\Psi, \beta}(\mathbf{x}) q(\mathbf{x}) d\mathbf{x}. \quad (42) $$

Therefore, exploiting (41) and (42), we can write the first bias term $\mathcal{B}_1$ in (38) as

$$ \mathcal{B}_1(x_\ell) = \mathbb{E}[\hat{\zeta}_{n,\ell}^{\beta}(x_\ell)] - \tilde{m}_n(x_\ell) - \left(\mathbb{E}\left[\int_{\mathbb{R}^d} \hat{m}_n^{\Psi, \beta}(\mathbf{x}) q(\mathbf{x}) d\mathbf{x}\right] - \int_{\mathbb{R}} m_\ell(u_\ell) q_\ell(u_\ell) du_\ell - C_n\right). \quad (43) $$

To further estimate $|\mathbb{E}[\hat{\zeta}_{n,\ell}^{\beta}(x_\ell)] - \tilde{m}_n(x_\ell)|$, observe:

$$ \begin{aligned} & |\mathbb{E}[\hat{\zeta}_{n,\ell}^{\beta}(x_\ell)] - \tilde{m}_n(x_\ell)| \\ &= \left|\mathbb{E}\left[\frac{2^{\lambda(n)}\tilde{Y}_{i,n}^{\beta}}{f_\ell(\mathbf{X}_{\ell,i})} K_\ell(2^{\lambda(n)}x_\ell, 2^{\lambda(n)}\mathbf{X}_{\ell,i})\right] - \tilde{m}_n(x_\ell)\right| \\ &= \left|\int_{\mathbb{R}} 2^{\lambda(n)} \tilde{m}_n(u_\ell) K_\ell(2^{\lambda(n)}x_\ell, 2^{\lambda(n)}u_\ell) du_\ell - \tilde{m}_n(x_\ell)\right| \\ &= \left|\int_{\mathbb{R}} (\tilde{m}_n(x_\ell + 2^{-\lambda(n)}v_\ell) - \tilde{m}_n(x_\ell)) K_\ell(2^{\lambda(n)}x_\ell, 2^{\lambda(n)}x_\ell + v_\ell) dv_\ell\right| \\ &\le \int_{\mathbb{R}} |\tilde{m}_n(x_\ell + 2^{-\lambda(n)}v_\ell) - \tilde{m}_n(x_\ell)| |K_\ell(2^{\lambda(n)}x_\ell, 2^{\lambda(n)}x_\ell + v_\ell)| dv_\ell \end{aligned} $$
$$ \leq 2^{-\lambda(n)(s-\frac{1}{p})} C \int_{\mathbb{R}} |v_{\ell}|^{(s-\frac{1}{p})} |K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell} + v_{\ell})| dv_{\ell}, $$

where in the final step we use the fact that $\tilde{m}_n(\cdot)$ is Lipschitz of order $(s - \frac{1}{p})$. Note additionally that each kernel $K_\ell(\cdot)$ has compact support, which simplifies bounding the above integrals. We begin by noting that under the assumption

$$ 2^{-\lambda(n)} \sim \left( \frac{\log n}{n} \right)^a, $$

it follows that

$$ \left| \mathbb{E}_{\tilde{\zeta}_{n,\ell}^{\beta}}(x_{\ell}) - \tilde{m}_n(x_{\ell}) \right| = \mathcal{O}\left( \left( \frac{\log n}{n} \right)^{a(s-\frac{1}{p})} \right). \quad (44) $$

To proceed, we next examine

$$
\begin{align*}
& \mathbb{E} \left[ \int_{\mathbb{R}^d} \hat{m}_n^{\Psi, \beta}(\mathbf{x}) q(\mathbf{x}) \, d\mathbf{x} - C_n \right] \\
&= \int_{\mathbb{R}^d} \mathbb{E} \left[ \hat{m}_n^{\Psi, \beta}(\mathbf{x}) \right] q(\mathbf{x}) \, d\mathbf{x} - \int_{\mathbb{R}^{d-1}} \sum_{\substack{j=1 \\ j \neq \ell}}^d m_j(z_j) g_n(\mathbf{z}_{-\ell}) \, d\mathbf{z}_{-\ell} - \mu \\
&= 2^{\lambda(n)d} \int_{\mathbb{R}^d} \mathbb{E} \left[ \frac{(\Psi(Y_i) - \mathbf{Z}_i^{\top} \boldsymbol{\beta}) K(2^{\lambda(n)}\mathbf{x}, 2^{\lambda(n)}\mathbf{X}_i)}{f(\mathbf{X}_i)} \right] q(\mathbf{x}) \, d\mathbf{x} \\
&\quad - \int_{\mathbb{R}^{d-1}} \sum_{\substack{j=1 \\ j \neq \ell}}^d m_j(z_j) g_n(\mathbf{z}_{-\ell}) \, d\mathbf{z}_{-\ell} - \mu \\
&= 2^{\lambda(n)d} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} m_{add}(\mathbf{u}) K(2^{\lambda(n)}\mathbf{x}, 2^{\lambda(n)}\mathbf{u}) q(\mathbf{x}) \, d\mathbf{x} \, d\mathbf{u} \\
&\quad - \int_{\mathbb{R}^{d-1}} \sum_{\substack{j=1 \\ j \neq \ell}}^d m_j(z_j) g_n(\mathbf{z}_{-\ell}) \, d\mathbf{z}_{-\ell} - \mu \\
&= 2^{\lambda(n)d} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \sum_{j=1}^n m_j(u_j) K(2^{\lambda(n)}\mathbf{x}, 2^{\lambda(n)}\mathbf{u}) q(\mathbf{x}) \, d\mathbf{x} \, d\mathbf{u} \\
&\quad - \int_{\mathbb{R}^{d-1}} \sum_{\substack{j=1 \\ j \neq \ell}}^d m_j(z_j) g_n(\mathbf{z}_{-\ell}) \, d\mathbf{z}_{-\ell} \\
&= 2^{\lambda(n)} \int_{\mathbb{R}} \int_{\mathbb{R}} m_{\ell}(u_{\ell}) K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}u_{\ell}) q_{\ell}(x_{\ell}) \, dx_{\ell} \, du_{\ell} \\
&\leq \int_{\mathbb{R}} m_{\ell}(x_{\ell}) q_{\ell}(x_{\ell}) \, dx_{\ell} + \int_{\mathbb{R}} q_{\ell}(x_{\ell}) \int_{\mathbb{R}} \left[ m_{\ell}(x_{\ell} + 2^{-\lambda(n)}v_{\ell}) - m_{\ell}(x_{\ell}) \right] \\
&\quad \times K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell} + v_{\ell}) \, dx_{\ell} \, dv_{\ell} \\
&\leq \int_{\mathbb{R}} m_{\ell}(x_{\ell}) q_{\ell}(x_{\ell}) \, dx_{\ell} \\
&\quad + \int_{\mathbb{R}} q_{\ell}(x_{\ell}) \int_{\text{Supp}(K_{\ell})} 2^{-\lambda(n)(s-\frac{1}{p})} |v_{\ell}|^{(s-\frac{1}{p})} K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell} + v_{\ell}) \, dx_{\ell} \, dv_{\ell} \\
&= \int_{\mathbb{R}} m_{\ell}(x_{\ell}) q_{\ell}(x_{\ell}) \, dx_{\ell} + \mathcal{O}\left(\left(\frac{\log n}{n}\right)^{a(s-\frac{1}{p})}\right). \tag{45}
\end{align*}
$$
Combining (44) and (45) immediately implies

$$ B_1 = O\left(\left(\frac{\log n}{n}\right)^a\left(s-\frac{1}{p}\right)\right). \quad (46) $$

We now consider the additional term:

$$ \begin{aligned} & |\mathbb{E}[\hat{m}_n^\beta(\mathbf{z})] - m(\mathbf{z})| \\ &= |\mathbb{E}\left(2^{\lambda(n)d} \sum_{i=1}^n \frac{\Psi(\mathbf{Y}_i) - \mathbf{Z}_i^\top \beta}{n f_{\mathbf{X}}(\mathbf{X}_i)} K(2^{\lambda(n)}\mathbf{z}, 2^{\lambda(n)}\mathbf{X}_i)\right) - m(\mathbf{z})| \\ &= \left|\int_{\mathbb{R}^d} 2^{\lambda(n)d} m(\mathbf{u}) K(2^{\lambda(n)}\mathbf{z}, 2^{\lambda(n)}\mathbf{u}) d\mathbf{u} - m(\mathbf{z})\right| \\ &= \left|\int_{\mathbb{R}^d} (m(\mathbf{z} + 2^{-\lambda(n)}\mathbf{v}) - m(\mathbf{z})) K(2^{\lambda(n)}\mathbf{z}, 2^{\lambda(n)}\mathbf{z} + \mathbf{v}) d\mathbf{v}\right| \\ &\le 2^{-\lambda(n)} \left(s - \frac{1}{p}\right) \left|\int_{\mathbb{R}^d} \|\mathbf{v}\|^{(s-\frac{1}{p})} K(2^{\lambda(n)}\mathbf{z}, 2^{\lambda(n)}\mathbf{z} + \mathbf{v}) d\mathbf{v}\right|. \end{aligned} $$

Hence, we conclude that

$$ B_2 = O\left(\left(\frac{\log n}{n}\right)^a\left(s-\frac{1}{p}\right)\right). \quad (47) $$

From (38), (46), and (47), we deduce the following:

$$ \sup_{\mathbf{x} \in C} |\mathbb{E}[\hat{m}_{add}^{\beta}(\mathbf{x})] - m(\mathbf{x})| = O\left(\left(\frac{\log n}{n}\right)^a\left(s-\frac{1}{p}\right)\right). \quad (48) $$

Finally, we consider the difference

$$ \left|\hat{m}_{add}^{\beta}(\mathbf{x}) - \hat{m}_{add}^{\beta}(t(\mathbf{x}))\right| = \left|\sum_{\ell=1}^{d} \left(\hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) - \hat{\zeta}_{n,\ell}^{\beta}(t(x_{\ell}))\right)\right|, $$

where

$$ \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) = \frac{2^{\lambda(n)}}{n} \sum_{i=1}^{n} \frac{\tilde{Y}_{i,n}^{\beta}}{f_{\ell}(X_{\ell,i})} K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell,i}), $$

and

$$ \tilde{Y}_{i,n}^{\beta} = \left(\Psi(\mathbf{Y}_i) - \mathbf{Z}_i^{\top} \beta\right) \times \int_{\mathbb{R}^{d-1}} 2^{(d-1)\lambda(n)} K_{-\ell}(2^{\lambda(n)}\mathbf{x}_{-\ell}, 2^{\lambda(n)}\mathbf{x}_{-\ell,i}) \frac{q_{-\ell}(\mathbf{x}_{-\ell})}{f(\mathbf{x}_{-\ell,i} | \mathbf{x}_{-\ell,i})} d\mathbf{x}_{-\ell}. $$

Since $K_\ell$ is a Lipschitz function, there exists a suitable Lipschitz constant $L$; thus,

$$ \left|\hat{m}_{add}^{\beta}(\mathbf{x}) - \hat{m}_{add}^{\beta}(t(\mathbf{x}))\right| = O(2^{2\lambda(n)} r_p), $$

where $r_p$ is the radius of $B_p(t_p, r_p)$. For

$$ r_p = \left(\frac{\varepsilon(n)}{2^{\lambda(n)}}\right)^{\gamma} \quad \text{with} \quad \gamma \ge 2, $$
we obtain

$$ \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\hat{m}_{\text{add}}^{\beta}(\mathbf{x}) - \hat{m}_{\text{add}}^{\beta}(t(\mathbf{x}))\right| \ge \varepsilon(n)/4\right) = 0, \quad (49) $$

provided $\varepsilon(n)$ is chosen sufficiently small. Similarly, one can establish that

$$ \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\mathbb{E}[\hat{m}_{\text{add}}^{\beta}(\mathbf{x})] - \mathbb{E}[\hat{m}_{\text{add}}^{\beta}(t(\mathbf{x}))]\right| \ge \varepsilon(n)/4\right) = 0. \quad (50) $$

Next, let $\{(p(n), q(n))\}_{n \in \mathbb{N}^*}$ be sequences of positive integers such that

$$ 2p(n)q(n) \le n < 2p(n)(q(n) + 1). $$

We then obtain

$$ \begin{aligned} & \hat{m}_{\text{add}}^{\beta}(t(\mathbf{x})) - \mathbb{E}[\hat{m}_{\text{add}}^{\beta}(t(\mathbf{x}))] \\ &= \sum_{\ell=1}^{d} \hat{\eta}_{\ell}^{\beta}(t(x_{\ell})) + \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) \, d\mathbf{z} \\ & \quad - \mathbb{E}\left[\sum_{\ell=1}^{d} \hat{\eta}_{\ell}^{\beta}(t(x_{\ell})) + \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) \, d\mathbf{z}\right] \\ &= \frac{1}{n} \sum_{j=1}^{q(n)} \xi_j^1(t(\mathbf{x})) + \frac{1}{n} \sum_{j=1}^{q(n)} \xi_j^2(t(\mathbf{x})) + \frac{1}{n} \delta_n(t(\mathbf{x})), \end{aligned} $$

where

$$ \xi_j^1 = \sum_{i=2(\ell-1)p(n)+1}^{(2j-1)p(n)} \xi_i - \mathbb{E}[\xi_i], $$

$$ \xi_j^2 = \sum_{i=(2j-1)p(n)+1}^{2jp(n)} \xi_i - \mathbb{E}[\xi_i], $$

$$ \delta_n = \sum_{i=2q(n)p(n)+1}^{n} \xi_i - \mathbb{E}[\xi_i], $$

$$ \begin{aligned} \xi_i = & \sum_{\ell=1}^{d} \frac{2^{\lambda(n)}}{n} \frac{\tilde{Y}_{i,n}^{\beta}}{f_{\ell}(X_{\ell,i})} K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell,i}) \\ & - \int_{\mathbb{R}} \frac{\tilde{Y}_{i,n}^{\beta}}{f_{\ell}(X_{\ell,i})} K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell,i}) \\ & + 2^{\lambda(n)d} \frac{\Psi(Y_i) - \mathbf{Z}_i^{\top}\boldsymbol{\beta}}{f(\mathbf{X}_i)} \int_{\mathbb{R}^d} K(2^{\lambda(n)}\mathbf{x}, 2^{\lambda(n)}\mathbf{x}_i) q(\mathbf{x}) \, d\mathbf{x}. \end{aligned} $$

It is straightforward to show that for $1 \le j \le q(n)$, there exist constants $m$ and $M$ such that

$$ m p(n) \le \xi_j^i \le M p(n) \quad \text{for } i = 1, 2. $$

Consequently, we have

$$ \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\hat{m}_{\text{add}}^{\beta}(t(\mathbf{x})) - \mathbb{E}[\hat{m}_{\text{add}}^{\beta}(t(\mathbf{x}))]\right| \ge \varepsilon(n)\right) $$
$$ \leq \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\sum_{j=1}^{q(n)} \xi_j^1(t(\mathbf{x}))\right| \geq \frac{2n \varepsilon(n)}{5}\right) + \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\sum_{j=1}^{q(n)} \xi_j^2(t(\mathbf{x}))\right| \geq \frac{2n \varepsilon(n)}{5}\right) \\ + \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\delta_n(t(\mathbf{x}))\right| \geq \frac{n \varepsilon(n)}{5}\right). \quad (51) $$

We begin by examining the leading term of (51). Specifically, we have

$$ \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\sum_{j=1}^{q(n)} \xi_j^1(t(\mathbf{x}))\right| \geq \frac{2n \varepsilon(n)}{5}\right) \leq r(n) \max_{p=1, \dots, r(n)} \mathbb{P}\left(\left|\sum_{j=1}^{q(n)} \xi_j^1(t_p)\right| \geq \frac{2n \varepsilon(n)}{5}\right). \quad (52) $$

Next, by choosing

$$ \zeta = \min\left(\frac{n\varepsilon(n)}{5q(n)}, \frac{m p(n)}{2}\right), $$

we obtain

$$ \left\{ \left| \sum_{j=1}^{q(n)} \xi_j^1 \right| \geq \frac{2n\varepsilon(n)}{5} \right\} \subset \left\{ \begin{aligned} & \left| \sum_{j=1}^{q(n)} \xi_j^1 \right| \geq \frac{2n\varepsilon(n)}{5} \quad \text{and} \quad \left| \xi_j^1 - \xi_j^{1*} \right|_{1 \leq j \leq q(n)} \leq \zeta \\ & \cup \left\{ \bigcup_{j=1}^{q(n)} \left| \xi_j^1 - \xi_j^{1*} \right| > \zeta \right\} \end{aligned} \right\} $$

where $\xi_j^{1*}$ are defined as in Lemma A1. Consequently,

$$ \mathbb{P}\left(\left|\sum_{j=1}^{q(n)} \xi_j^1\right| \geq \frac{2n \varepsilon(n)}{5}\right) \leq \mathbb{P}\left(\left|\sum_{j=1}^{q(n)} \xi_j^{1*}\right| \geq \frac{2n \varepsilon(n)}{5} - q(n) \zeta\right) \\ + \sum_{j=1}^{q(n)} \mathbb{P}\left(\left|\xi_j^1 - \xi_j^{1*}\right| \geq \zeta\right). \quad (53) $$

By applying Hoeffding's inequality [124], we deduce

$$ \begin{aligned} \mathbb{P}\left(\left|\sum_{j=1}^{q(n)} \xi_j^{1*}\right| \geq \frac{2n \varepsilon(n)}{5} - q(n) \zeta\right) & \leq \mathbb{P}\left(\left|\sum_{j=1}^{q(n)} \xi_j^{1*}\right| \geq \frac{n \varepsilon(n)}{5}\right) \\ & \leq 2 \exp\left(-\frac{2 n^2 \varepsilon^2(n) / 5^2}{q(n) 4 M^2 p(n)^2}\right) \quad (54) \\ & \leq 2 \exp\left(-\frac{n \varepsilon^2(n)}{5^2 M^2 p(n)}\right). \end{aligned} $$

Next, applying Bradley's Lemma A1 with $\gamma = +\infty$ gives

$$ \begin{aligned} \mathbb{P}\left(\left|\xi_j - \xi_j^*\right| \geq \zeta\right) & \leq 11 \sqrt{\frac{M p(n)}{\min\left(\frac{n \varepsilon(n)}{5 p(n)}, \frac{m p(n)}{2}\right)}} \alpha(p) \\ & \leq 11 \left(\frac{5}{2}\right)^{\frac{1}{2}} \left(\min\left(\frac{\varepsilon(n)}{M}, \frac{5m}{4M}\right)\right)^{-\frac{1}{2}} \alpha(p), \end{aligned} $$

since $n \geq 2p(n)q(n)$. If we additionally require $\varepsilon(n) \leq M$, then

$$ \min\left(\frac{\varepsilon(n)}{M}, \frac{5m}{4M}\right) \geq \frac{\varepsilon(n) M}{m^2}. $$
Hence,

$$ \mathbb{P}(|\xi_j - \xi_j^*| \ge \zeta) \le 11 \left(\frac{5}{2}\right)^{\frac{1}{2}} \sqrt{\frac{M^2}{m\epsilon(n)}} \alpha(p). \quad (55) $$

Combining (53), (54), and (55) then yields

$$ \mathbb{P}\left(\left|\sum_{j=1}^{q(n)} \xi_j^1\right| \ge \frac{2n\epsilon(n)}{5}\right) \le 2 \exp\left(-\frac{n\epsilon^2(n)}{5^2 M^2 p(n)}\right) \\ + 11 \left(\frac{5}{2}\right)^{\frac{1}{2}} \sqrt{\frac{M^2}{m\epsilon(n)}} q(n) \alpha(p). \quad (56) $$

Finally, by combining (52) and (56), we obtain

$$ \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\sum_{j=1}^{q(n)} \xi_j^1(t(\mathbf{x}))\right| \ge \frac{2n\epsilon(n)}{5}\right) \le r(n) \left(2 \exp\left(-\frac{n\epsilon^2(n)}{5^2 M^2 p(n)}\right) \right. \\ \left. + 11 \left(\frac{5}{2}\right)^{\frac{1}{2}} \sqrt{\frac{M^2}{m\epsilon(n)}} q(n) \alpha(p)\right). \quad (57) $$

By the same reasoning, we can also prove

$$ \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\sum_{j=1}^{q(n)} \xi_j^2(t(\mathbf{x}))\right| \ge \frac{2n\epsilon(n)}{5}\right) \le r(n) \left(2 \exp\left(-\frac{n\epsilon^2(n)}{5^2 M^2 p(n)}\right) \right. \\ \left. + 11 \left(\frac{5}{2}\right)^{\frac{1}{2}} \sqrt{\frac{M^2}{m\epsilon(n)}} q(n) \alpha(p)\right). \quad (58) $$

Now, let $v > 0$. Then we have

$$ \begin{aligned} \mathbb{P}\left(\left|\delta_n\right| \ge \frac{n\epsilon(n)}{5}\right) &\le \exp\left(-\frac{vn\epsilon(n)}{5}\right) \mathbb{E}\exp(v|\delta_n|) \\ &\le \exp\left(-v\left(\frac{n\epsilon(n)}{5} - M(n - 2p(n)q(n))\right)\right). \end{aligned} $$

Since $Mv(n - 2p(n)q(n)) = 1$ and $n < 2p(n)(q(n) + 1)$, we obtain

$$ \begin{aligned} \mathbb{P}\left(\left|\delta_n\right| \ge \frac{n\epsilon(n)}{5}\right) &\le e \exp\left(-\frac{n\epsilon(n)}{5M(n-2p(n)q(n))}\right) \\ &\le e \exp\left(-\frac{\epsilon(n)}{5M} q(n)\right). \end{aligned} $$

In addition, if $\epsilon(n) \le M$, then

$$ \mathbb{P}\left(\left|\delta_n\right| \ge \frac{n\epsilon(n)}{5}\right) \le e \exp\left(-\frac{\epsilon(n)^2}{5M^2} q(n)\right). \quad (59) $$

Hence,

$$ \begin{aligned} \mathbb{P}\left(\sup_{\mathbf{x} \in C} \left|\delta_n(t(\mathbf{x}))\right| \ge \frac{n\epsilon(n)}{5}\right) &\le r(n) \max_{p=1,\dots,r(n)} \mathbb{P}\left(\left|\delta_n(t_p)\right| \ge \frac{n\epsilon(n)}{5}\right) \\ &\le e r(n) \exp\left(-\frac{\epsilon(n)^2}{5M^2} q(n)\right). \end{aligned} \quad (60) $$

Recall that

$$ 2p(n)q(n) \le n < 2p(n)(q(n) + 1), $$
which implies

$$q(n) = \frac{n}{2p(n)} - \theta, \quad 0 \le \theta < 1.$$

Then, we have

$$ \begin{aligned} \mathbb{P}\left(\sup_{\mathbf{x} \in C} |\delta_n(t(\mathbf{x}))| \ge \frac{n\epsilon(n)}{5}\right) &\le e r(n) \exp\left(-\frac{n\epsilon(n)^2}{10M^2p(n)}\right) \exp\left(\frac{\epsilon(n)^2}{5M^2}\right) \\ &\le e^{6/5} r(n) \exp\left(-\frac{n\epsilon(n)^2}{10M^2p(n)}\right), \end{aligned} \quad (61) $$

since $\epsilon(n) \le M$.

It follows from (51), (57), (58), and (61) that

$$ \begin{aligned} & \mathbb{P}\left(\sup_{\mathbf{x} \in C} |\hat{m}_{\text{add}}^{\beta}(t(\mathbf{x})) - E[\hat{m}_{\text{add}}^{\beta}(t(\mathbf{x}))]| \ge \epsilon(n)\right) \\ &\le 2\left(r(n)\left(2\exp\left(-\frac{n\epsilon^2(n)}{25M^2p(n)}\right) + 11\left(\frac{5}{2}\right)^{\frac{1}{2}}\sqrt{\frac{M^2}{m\epsilon(n)}}q(n)\alpha(p)\right)\right. \\ &\qquad \left. + e^{6/5}r(n)\exp\left(-\frac{n\epsilon(n)^2}{10M^2p(n)}\right)\right) \\ &\le r(n)\left(4 + e^{6/5}\right)\exp\left(-\frac{n\epsilon(n)^2}{25M^2p(n)}\right) + 22r(n)\left(\frac{5}{2}\right)^{\frac{1}{2}}\sqrt{\frac{M^2}{m\epsilon(n)}}q(n)\alpha(p). \end{aligned} \quad (62) $$

Therefore, by the Borel–Cantelli lemma, we obtain

$$ \sum_{n \ge 1} \mathbb{P} \left( \sup_{\mathbf{x} \in C} \left| \hat{m}_{\text{add}}^{\beta}(t(\mathbf{x})) - \mathbb{E}[\hat{m}_{\text{add}}^{\beta}(t(\mathbf{x}))] \right| \ge \epsilon(n) \right) < \infty. \quad (63) $$

We begin by observing the following relationship:

$$ \begin{aligned} & \left| \hat{m}_{\text{add}}^{\beta}(\mathbf{x}) - \hat{m}_{\text{add}}^{\beta}(\mathbf{x}) \right| \\ &= \sum_{\ell=1}^{d} \hat{\eta}_{\ell}^{\beta}(x_{\ell}) + \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) \, d\mathbf{z} - \sum_{\ell=1}^{d} \hat{\eta}_{\ell}^{\beta}(x_{\ell}) - \int_{\mathbb{R}^d} \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) \, d\mathbf{z} \\ &= \sum_{\ell=1}^{d} \left( \hat{\eta}_{\ell}^{\beta}(x_{\ell}) - \hat{\eta}_{\ell}^{\beta}(x_{\ell}) \right) + \int_{\mathbb{R}^d} \left( \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) \, d\mathbf{z} - \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) \, d\mathbf{z} \right). \end{aligned} \quad (64) $$

Recall that

$$ \hat{\eta}_{\ell}^{\beta}(x_{\ell}) = \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) - \int_{\mathbb{R}} \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) q_{\ell}(\mathbf{x}_{\ell}) \, d\mathbf{x}_{\ell}, $$

where

$$ \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) = \frac{2^{\lambda(n)}}{n} \sum_{i=1}^{n} \frac{\tilde{Y}_{i,n}^{\beta}}{f_{\ell}(X_{\ell,i})} K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}X_{\ell,i}), $$

and

$$ \tilde{Y}_{i,n}^{\beta} = (\Psi(\mathbf{Y}_i) - \mathbf{Z}_i^{\top} \boldsymbol{\beta}) \int_{\mathbb{R}^{d-1}} 2^{(d-1)\lambda(n)} K_{-\ell}(2^{\lambda(n)}\mathbf{x}_{-\ell}, 2^{\lambda(n)}\mathbf{X}_{-\ell,i}) \frac{q_{-\ell}(\mathbf{x}_{-\ell})}{f(\mathbf{X}_{-\ell,i} | \mathbf{X}_{\ell,i})} d\mathbf{x}_{-\ell}. $$

Consequently,

$$ \begin{aligned} & \hat{\eta}_{\ell}^{\beta}(x_{\ell}) - \hat{\eta}_{\ell}^{\beta}(x_{\ell}) \\ &= \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) - \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) - \int_{\mathbb{R}} (\hat{\zeta}_{n,\ell}^{\beta}(x_{\ell}) - \hat{\zeta}_{n,\ell}^{\beta}(x_{\ell})) q_{\ell}(\mathbf{x}_{\ell}) \, d\mathbf{x}_{\ell} \end{aligned} $$
$$ = \sum_{i=1}^{n} \frac{2^{\lambda(n)} (\tilde{Y}_{i,n}^{\hat{\beta}} - \tilde{Y}_{i,n}^{\beta})}{n f_{\ell}(X_{\ell,i})} \left[ K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell,i}) - \int_{\mathbb{R}} K_{\ell}(2^{\lambda(n)}x_{\ell}, 2^{\lambda(n)}x_{\ell,i}) q_{\ell}(x_{\ell}) dx_{\ell} \right]. \quad (65) $$

It is immediately apparent that

$$ \begin{aligned} & \tilde{Y}_{i,n}^{\hat{\beta}} - \tilde{Y}_{i,n}^{\beta} \\ &= Z_{i}^{\top}(\hat{\beta} - \beta) \int_{\mathbb{R}^{d-1}} 2^{(d-1)\lambda(n)} K_{-\ell}(2^{\lambda(n)}x_{-\ell}, 2^{\lambda(n)}x_{-\ell,i}) \frac{q_{-\ell}(x_{-\ell})}{f(x_{-\ell,i} | X_{\ell,i})} dx_{-\ell}. \end{aligned} \quad (66) $$

Observe that

$$ \begin{aligned} & \sqrt{\frac{n}{2^{\lambda(n)}}} (\hat{\beta} - \beta) \\ &= (n^{-1} \tilde{Z}^{\top} \tilde{Z})^{-1} \frac{1}{\sqrt{n} 2^{\lambda(n)}} \left[ \sum_{i=1}^{n} \tilde{Z}_{i} \tilde{m}_{add}(X_{i}) - \sum_{i=1}^{n} \tilde{Z}_{i} \sum_{j=1}^{n} W_{nj}(X_{j}) \varepsilon_{j} + \sum_{i=1}^{n} \tilde{Z}_{i} \varepsilon_{i} \right]. \end{aligned} \quad (67) $$

Then,

$$ \sum_{i=1}^{n} \tilde{Z}_{i} \tilde{m}_{add}(X_{i}) = \sum_{i=1}^{n} Z_{i} \tilde{m}_{add}(X_{i}) - \sum_{i=1}^{n} \sum_{j=1}^{n} W_{nj}(X_{i}) Z_{j} \tilde{m}_{add}(X_{i}), \quad (68) $$

where

$$ \begin{aligned} \max_{1 \le i \le n} |\tilde{m}_{add}(X_i)| &= \max_{1 \le i \le n} |m_{add}(X_i) - \sum_{j=1}^{n} m_{add}(X_j) W_{nj}(X_i)| \\ &\le \max_{1 \le i \le n} |m_{add}(X_i) - \tilde{m}_{add}^{\Psi, \beta}(X_i)| + \max_{1 \le i \le n} \left| \sum_{j=1}^{n} W_{nj}(X_i) \varepsilon_j \right|. \end{aligned} \quad (69) $$

By Lemma 1 of [125], under hypotheses (G.1-3), (K.1-3), (M.1), and (Q.1-2), it follows that

$$ \max_{1 \le i \le n} |m_{add}(X_i) - \tilde{m}_{add}^{\Psi, \beta}(X_i)| = O\left(\sqrt{\frac{2^{\lambda(n)} \log n}{n}}\right). \quad (70) $$

Moreover, it is straightforward to see that

$$ \max_{1 \le i \le n} \left| \sum_{j=1}^{n} W_{nj}(X_i) \right| = O(1). $$

By Theorem 1.1 of [126], choosing $c_n = \sqrt{2n \log_2(n)}$ yields

$$ \sum_{i=1}^{n} Z_{i} = O\left(\sqrt{2n \log_{2}(n)}\right). \quad (71) $$

Hence, combining (69), (70), and (71) with hypothesis (H.2), we obtain

$$ \sum_{i=1}^{n} Z_{i} \tilde{m}_{add}(X_{i}) = o(\sqrt{n 2^{\lambda(n)}}). \quad (72) $$

Next, consider

$$ \sum_{i=1}^{n} \sum_{j=1}^{n} W_{nj}(X_i) Z_j \tilde{m}_{add}(X_i) $$
$$ \leq n \max_{1 \leq i \leq n} \hat{m}_{add}(\mathbf{X}_i) \max_{1 \leq i \leq n} \max_{1 \leq j \leq n} W_{nj}(\mathbf{X}_i) \sum_{j=1}^{n} \mathbf{Z}_j. \quad (73) $$

A similar argument, using

$$ \max_{1 \leq i \leq n} \max_{1 \leq j \leq n} W_{nj}(\mathbf{X}_i) = O(n^{-1} 2^{\lambda(n)}), $$

implies

$$ \sum_{i=1}^{n} \sum_{j=1}^{n} W_{nj}(\mathbf{X}_i) \mathbf{Z}_j \hat{m}_{add}(\mathbf{X}_i) = o(\sqrt{n 2^{\lambda(n)}}). \quad (74) $$

Therefore, by (68), (72), and (74), it follows that

$$ \sum_{i=1}^{n} \tilde{\mathbf{Z}}_i \hat{m}_{add}(\mathbf{X}_i) = o(\sqrt{n 2^{\lambda(n)}}). \quad (75) $$

Applying the same arguments, we also deduce

$$ \sum_{i=1}^{n} \tilde{\mathbf{Z}}_i \sum_{j=1}^{n} W_{nj}(\mathbf{X}_i) \varepsilon_j = o(\sqrt{n 2^{\lambda(n)}}), \quad (76) $$

and

$$ \sum_{i=1}^{n} \tilde{\mathbf{Z}}_i \varepsilon_i = o(\sqrt{n 2^{\lambda(n)}}). \quad (77) $$

Finally, combining (67), (75), (76), and (77) with the fact that

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} \tilde{\mathbf{Z}}_{i}^{\top} \tilde{\mathbf{Z}}_{i} = B, $$

we conclude that

$$ \sqrt{\frac{n}{2^{\lambda(n)}}} (\hat{\beta} - \beta) = o(1). \quad (78) $$

Hence, from (65), (66), (67), (78), and hypothesis (H.3), we obtain

$$ \sum_{\ell=1}^{d} (\hat{\eta}_{\ell}^{\hat{\beta}}(x_{\ell}) - \hat{\eta}_{\ell}^{\beta}(x_{\ell})) = o(\frac{\log n}{n})^{a(s-\frac{1}{p})}. \quad (79) $$

Next, observe that

$$ \begin{aligned} & \hat{m}_{n}^{\hat{\beta}}(\mathbf{z}) q(\mathbf{z}) d\mathbf{z} - \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) d\mathbf{z} \\ & = 2^{\lambda(n)d} \sum_{i=1}^{n} \frac{\mathbf{Z}_{i}^{\top} (\hat{\beta} - \beta)}{n f_{\chi}(\mathbf{X}_{i})} K(2^{\lambda(n)}\mathbf{z}, 2^{\lambda(n)}\mathbf{X}_{i}). \end{aligned} \quad (80) $$

By (78) and hypothesis (H.3), it follows that

$$ \int_{\mathbb{R}^d} (\hat{m}_{n}^{\hat{\beta}}(\mathbf{z}) q(\mathbf{z}) d\mathbf{z} - \hat{m}_{n}^{\beta}(\mathbf{z}) q(\mathbf{z}) d\mathbf{z}) = o(\frac{\log n}{n})^{a(s-\frac{1}{p})}. \quad (81) $$
Therefore, by (64), (79), and (81), we conclude that

$$ \sup_{x \in C} |\hat{m}_{add}^{\beta}(\mathbf{x}) - \hat{m}_{add}^{\beta}(\mathbf{x})| = o\left(\left(\frac{\log n}{n}\right)^a \left(s - \frac{1}{p}\right)\right). \quad (82) $$

Finally, combining (37), (48), (49), (50), 63, and (82) completes the proof of the theorem. Observe that, up to this point, we have focused on the scenario in which the density function $f_X(\cdot)$ is known. When $f_X(\cdot)$ is unknown, however, it can be estimated by the empirical function $f_n(\cdot)$. In this case, we use the following decomposition:

$$ \frac{1}{f_n} = \frac{1}{f_X} + \frac{f_X - f_n}{f_X f_n}. $$

According to [70], under conditions (G.1)–(G.3), it holds on every compact subset $I^d \subset \mathbb{R}^d$ that

$$ \sup_{x \in I^d} |f_n(\mathbf{x}) - f_X(\mathbf{x})| = O\left(\left(\frac{\log n}{n}\right)^{\frac{s}{d+2s}}\right). \quad (83) $$

By combining this result with the preceding findings, we can now complete the proof of the theorem. □

**Author Contributions:** K.C.: Conceptualization, Methodology, Investigation, Writing—Original Draft, Writing—Review and Editing; S.B.: Conceptualization, Methodology, Investigation, Writing—Original Draft, Writing—Review and Editing. Both authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** No new data were created or analyzed in this study.

**Acknowledgments:** The authors wish to extend their sincere gratitude to the Editor-in-Chief, the Associate Editor, and the six referees for their exceptionally thoughtful and constructive feedback. Their rigorous and detailed review process not only illuminated several insights in our initial submission but also provided insightful guidance that has significantly enhanced the depth, clarity, and focus of our work. We are deeply appreciative of the academic rigor and kindness embodied in their suggestions, which have been instrumental in refining our presentation and elevating the overall quality of this research.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

# Appendix A

**Lemma A1.** Let $\zeta$ and $\gamma$ be positive constants, and consider random variables $\zeta_1, \dots, \zeta_q$ such that

$$ \zeta \le \|\zeta_j\|_{\gamma} < \infty $$

for all $j = 1, \dots, q$. Then there exist independent random variables $\zeta_1^*, \dots, \zeta_q^*$ with $\mathbb{P}_{\zeta_j^*} = \mathbb{P}_{\zeta_j}$ for each $j$, and

$$ \mathbb{P}(|\zeta_j^* - \zeta_j| \ge \zeta) \le 11 \left( \frac{\|\zeta_j\|_{\gamma}}{\zeta} \right)^{\frac{\gamma}{2\gamma+1}} \alpha(p)^{\frac{2\gamma}{2\gamma+1}}, \quad (A1) $$

where $\alpha(p)$ denotes the strong mixing coefficients given by

$$ \alpha(p) = \sup_{t} \sup_{A \in \mathcal{F}_{-\infty}^{t}, B \in \mathcal{F}_{t+p}^{+\infty}} |\mathbb{P}(A \cap B) - \mathbb{P}(A)\mathbb{P}(B)|, $$
for $p = 1, 2, \dots$, and $\mathfrak{F}_{-\infty}^t$ and $\mathfrak{F}_{t+p}^{+\infty}$ are the $\sigma$-fields generated by the centered real-valued stochastic process $(X_s)_{s \in \mathbb{Z}}$, with $X_s$ for $s \le t$ and $s \ge t + p$, respectively.

**Proof.** See Theorem 3 of [127] and p. 202 of [128]. □

Besov Spaces

Following [71], suppose $1 \le p, q \le \infty$. Define the shift operator

$$(S_{\tau}f)(x) = f(x - \tau).$$

For $0 < s < 1$, set

$$\gamma_{s,p,q}(f) = \left( \int_{\mathbb{R}^d} \left( \frac{\|S_{\tau}f - f\|_{L_p}}{\|\tau\|^s} \right)^q \frac{d\tau}{\|\tau\|^d} \right)^{1/q},$$

and

$$\gamma_{s,p,\infty}(f) = \sup_{\tau \in \mathbb{R}^d} \frac{\|S_{\tau}f - f\|_{L_p}}{\|\tau\|^s}.$$

If $s = 1$, define

$$\gamma_{1,p,q}(f) = \left( \int_{\mathbb{R}^d} \left( \frac{\|S_{\tau}f + S_{-\tau}f - 2f\|_{L_p}}{\|\tau\|} \right)^q \frac{d\tau}{\|\tau\|^d} \right)^{1/q},$$

and

$$\gamma_{1,p,\infty}(f) = \sup_{\tau \in \mathbb{R}^d} \frac{\|S_{\tau}f + S_{-\tau}f - 2f\|_{L_p}}{\|\tau\|}.$$

For $0 < s < 1$ and $1 \le p, q \le \infty$, define

$$\mathbf{B}_{s,p,q} = \{f \in L_p(\mathbb{R}^d) : \gamma_{s,p,q}(f) < \infty\}.$$

For $s > 1$, write

$$s = [s]^{-} + \{s\}^{+},$$

where $[s]^{-}$ is an integer and $0 < \{s\}^{+} \le 1$. Then, $\mathbf{B}_{s,p,q}$ consists of all functions $f \in L_p(\mathbb{R}^d)$ such that $D^j f \in \mathbf{B}_{\{s\}^{+},p,q}$ for every multi-index $j$ with $|j| \le [s]^{-}$. The norm is given by

$$\|f\|_{B_{s,p,q}} = \|f\|_{L_p} + \sum_{|j| \le [s]^{-}} \gamma_{\{s\}^{+},p,q}(D^j f).$$

Notable examples of Besov spaces include the Sobolev space $H_2^s = \mathbf{B}_{s,2,2}$ and the space of bounded s-Lipschitz functions $\mathbf{B}_{s,\infty,\infty}$.

References

1. Liu, Y.; Lu, M.; McMahan, C.S. A penalized likelihood approach for efficiently estimating a partially linear additive transformation model with current status data. *Electron. J. Stat.* **2021**, 15, 2247–2287. [CrossRef]

2. Fan, J.; Li, R. Variable selection via nonconcave penalized likelihood and its oracle properties. *J. Am. Stat. Assoc.* **2001**, 96, 1348–1360. [CrossRef]

3. Härdle, W.; Mammen, E.; Müller, M. Testing parametric versus semiparametric modeling in generalized linear models. *J. Am. Stat. Assoc.* **1998**, 93, 1461–1474. [CrossRef]

4. Müller, M.; Rönz, B. Credit scoring using semiparametric methods. In *Measuring Risk in Complex Stochastic Systems*; Springer: New York, NY, USA, 2000; pp. 83–97.

5. Zhang, J.; Zhou, Y.; Lin, B.; Yu, Y. Estimation and hypothesis test on partial linear models with additive distortion measurement errors. *Comput. Stat. Data Anal.* **2017**, 112, 114–128. [CrossRef]
6. Wang, L.; Liu, X.; Liang, H.; Carroll, R.J. Estimation and variable selection for generalized additive partial linear models. *Ann. Stat.* **2011**, *39*, 1827–1851. [CrossRef]

7. Zhang, X.; Liang, H. Focused information criterion and model averaging for generalized additive partial linear models. *Ann. Stat.* **2011**, *39*, 174–200. [CrossRef]

8. Tian, Y.; Jiang, B. Equalities for estimators of partial parameters under linear model with restrictions. *J. Multivar. Anal.* **2016**, *143*, 299–313. [CrossRef]

9. Brunk, H.D. On the estimation of parameters restricted by inequalities. *Ann. Math. Stat.* **1958**, *29*, 437–454. [CrossRef]

10. Engle, R.F.; Granger, C.W.J.; Rice, J.; Weiss, A. Semiparametric Estimates of the Relation Between Weather and Electricity Sales. *J. Am. Stat. Assoc.* **1986**, *81*, 310–320. [CrossRef]

11. Speckman, P. Kernel smoothing in partial linear models. *J. R. Stat. Soc. Ser. B* **1988**, *50*, 413–436. [CrossRef]

12. Liang, H.; Härdle, W.; Carroll, R.J. Estimation in a semiparametric partially linear errors-in-variables model. *Ann. Stat.* **1999**, *27*, 1519–1535. [CrossRef]

13. Severini, T.A.; Staniswalis, J.G. Quasi-likelihood estimation in semiparametric models. *J. Am. Stat. Assoc.* **1994**, *89*, 501–511. [CrossRef]

14. Rice, J. Convergence rates for partially splined models. *Stat. Probab. Lett.* **1986**, *4*, 203–208. [CrossRef]

15. Chen, H. Convergence rates for parametric components in a partly linear model. *Ann. Stat.* **1988**, *16*, 136–146. [CrossRef]

16. Robinson, P.M. Root-N-consistent semiparametric regression. *Econometrica* **1988**, *56*, 931–954. [CrossRef]

17. Chen, H.; Shiau, J.J.H. A two-stage spline smoothing method for partially linear models. *J. Stat. Plann. Inference* **1991**, *27*, 187–201. [CrossRef]

18. Eubank, R.L.; Speckman, P. Curve fitting by polynomial-trigonometric regression. *Biometrika* **1990**, *77*, 1–9. [CrossRef]

19. Donald, S.G.; Newey, W.K. Series estimation of semilinear models. *J. Multivar. Anal.* **1994**, *50*, 30–40. [CrossRef]

20. Shi, P.D.; Li, G.Y. Asymptotic normality of the M-estimators for parametric components in partly linear models. *Northeast. Math. J.* **1995**, *11*, 127–138.

21. Shi, P.D.; Li, G.Y. A note on the convergent rates of M-estimates for a partly linear model. *Statistics* **1995**, *26*, 27–47. [CrossRef]

22. Hamilton, S.A.; Truong, Y.K. Local linear estimation in partly linear models. *J. Multivar. Anal.* **1997**, *60*, 1–19. [CrossRef]

23. Härdle, W.; Liang, H.; Gao, J. *Partially Linear Models*; Contributions to Statistics; Physica-Verlag: Heidelberg, Germany, 2000; p. x+203. [CrossRef]

24. Lee, E.R.; Han, K.; Park, B.U. Estimation of errors-in-variables partially linear additive models. *Stat. Sin.* **2018**, *28*, 2353–2373.

25. Boente, G.; Martínez, A.M. Robust variable selection for partially linear additive models. *Stat. Comput.* **2024**, *34*, 201. [CrossRef]

26. Maidman, A.; Wang, L.; Zhou, X.H.; Sherwood, B. Quantile partially linear additive model for data with dropouts and an application to modeling cognitive decline. *Stat. Med.* **2023**, *42*, 2729–2745. [CrossRef] [PubMed]

27. Yang, X.; Li, L.; Wu, H.; Xu, W. Composite quantile regression for partially linear additive model with censored responses and its application. *Chin. J. Appl. Probab. Stat.* **2023**, *39*, 604–622.

28. Li, Z.; Song, Y. Two-stage Walsh-average-based robust estimation and variable selection for partially linear additive spatial autoregressive models. *Braz. J. Probab. Stat.* **2023**, *37*, 667–692. [CrossRef]

29. Mou, X.; Wang, D. Additive partially linear model for pooled biomonitoring data. *Comput. Stat. Data Anal.* **2024**, *190*, 107862. [CrossRef]

30. Cai, T.; Li, J.; Zhou, Q.; Yin, S.; Zhang, R. Subgroup detection based on partially linear additive individualized model with missing data in response. *Comput. Stat. Data Anal.* **2024**, *192*, 107910. [CrossRef]

31. Zhao, W.; Li, R.; Lian, H. Estimation and variable selection of quantile partially linear additive models for correlated data. *J. Stat. Comput. Simul.* **2024**, *94*, 315–345. [CrossRef]

32. Lu, F.; Tian, G.; Yang, J. GMM estimation and variable selection of partially linear additive spatial autoregressive model. *Stat. Pap.* **2024**, *65*, 2253–2288. [CrossRef]

33. Lu, M.; Li, C.S.; Wagner, K.D. Penalised estimation of partially linear additive zero-inflated Bernoulli regression models. *J. Nonparametr. Stat.* **2024**, *36*, 863–890. [CrossRef]

34. Stone, C.J. Additive regression and other nonparametric models. *Ann. Stat.* **1985**, *13*, 689–705. [CrossRef]

35. Stone, C.J. The dimensionality reduction principle for generalized additive models. *Ann. Stat.* **1986**, *14*, 590–606. [CrossRef]

36. Fan, J.; Gijbels, I. Local polynomial modelling and its applications. In *Monographs on Statistics and Applied Probability; Chapman & Hall: London, UK*, 1996; Volume 66, p. xvi+341.

37. Härdle, W. Applied nonparametric regression. In *Econometric Society Monographs*; Cambridge University Press: Cambridge, UK, 1990; Volume 19, p. xvi+333. [CrossRef]

38. Manzan, S.; Zerom, D. Kernel estimation of a partially linear additive model. *Stat. Probab. Lett.* **2005**, *72*, 313–322. [CrossRef]

39. Chamberlain, G. Efficiency bounds for semiparametric regression. *Econometrica* **1992**, *60*, 567–596. [CrossRef]

40. Yu, K.; Mammen, E.; Park, B.U. Semi-parametric regression: Efficiency gains from modeling the nonparametric part. *Bernoulli* **2011**, *17*, 736–748. [CrossRef]
41. Bouzebda, S.; Didi, S. Additive regression model for stationary and ergodic continuous time processes. *Comm. Stat. Theory Methods* **2017**, *46*, 2454–2493. [CrossRef]

42. Bouzebda, S.; Chokri, K. Statistical tests in the partially linear additive regression models. *Stat. Methodol.* **2014**, *19*, 4–24. [CrossRef]

43. Chokri, K.; Bouzebda, S. Uniform-in-bandwidth consistency results in the partially linear additive model components estimation. *Commun. Stat. Theory Methods* **2024**, *53*, 3383–3424. [CrossRef]

44. Tang, Q.; Tu, W.; Kong, L. Estimation for partial functional partially linear additive model. *Comput. Stat. Data Anal.* **2023**, *177*, 107584. [CrossRef]

45. Du, J.; Cao, R.; Kwessi, E.; Zhang, Z. Estimation for generalized partially functional linear additive regression model. *J. Appl. Stat.* **2019**, *46*, 914–925. [CrossRef]

46. Ding, F.; Chen, J. Penalized quadratic inference estimation for partially linear additive dynamic panel model with fixed effects. *Appl. Math. Ser. A (Chin. Ed.)* **2018**, *33*, 21–35.

47. Liu, L.; Li, J.; Zhang, R. General partially linear additive transformation model with right-censored data. *J. Appl. Stat.* **2014**, *41*, 2257–2269. [CrossRef]

48. Ma, S.; Yang, L. Spline-backfitted kernel smoothing of partially linear additive model. *J. Stat. Plann. Inference* **2011**, *141*, 204–219. [CrossRef]

49. Meyer, Y. Wavelets and operators. In *Cambridge Studies in Advanced Mathematics*; Cambridge University Press: Cambridge, UK, 1992; Volume 37, p. xvi+224.

50. Daubechies, I. Ten lectures on wavelets. In *CBMS-NSF Regional Conference Series in Applied Mathematics*; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, USA, 1992; Volume 61, p. xx+357. [CrossRef]

51. Mallat, S. *A Wavelet Tour of Signal Processing*, 3rd ed.; Elsevier/Academic Press: Amsterdam, The Netherlands, 2009; p. xxii+805.

52. Vidakovic, B. *Statistical Modeling by Wavelets*; Wiley Series in Probability and Statistics: Applied Probability and Statistics; John Wiley & Sons, Inc.: New York, NY, USA, 1999; p. xiv+382. [CrossRef]

53. Härdle, W.; Kerkyacharian, G.; Picard, D.; Tsybakov, A. *Wavelets, Approximation, and Statistical Applications*; Lecture Notes in Statistics; Springer: New York, NY, USA, 1998; Volume 129, p. xviii+265. [CrossRef]

54. Hall, P.; Patil, P. Formulae for mean integrated squared error of nonlinear wavelet-based density estimators. *Ann. Stat.* **1995**, *23*, 905–928. [CrossRef]

55. Hall, P.; Patil, P. On wavelet methods for estimating smooth functions. *Bernoulli* **1995**, *1*, 41–58. [CrossRef]

56. Donoho, D.L.; Johnstone, I.M.; Kerkyacharian, G.; Picard, D. Density estimation by wavelet thresholding. *Ann. Stat.* **1996**, *24*, 508–539. [CrossRef]

57. Hall, P.; Penev, S. Cross-validation for choosing resolution level for nonlinear wavelet curve estimators. *Bernoulli* **2001**, *7*, 317–341. [CrossRef]

58. Giné, E.; Nickl, R. Uniform limit theorems for wavelet density estimators. *Ann. Probab.* **2009**, *37*, 1605–1646. [CrossRef]

59. Rademacher, D.; Krebs, J.; von Sachs, R. Statistical inference for wavelet curve estimators of symmetric positive definite matrices. *J. Stat. Plann. Inference* **2024**, *231*, 106140. [CrossRef]

60. Püspöki, Z.; Unser, M. Template-free wavelet-based detection of local symmetries. *IEEE Trans. Image Process.* **2015**, *24*, 3009–3018. [CrossRef]

61. Krivoshein, A.V. On construction of multivariate symmetric MRA-based wavelets. *Appl. Comput. Harmon. Anal.* **2014**, *36*, 215–238. [CrossRef]

62. Cattani, C. On the existence of wavelet symmetries in archaea DNA. *Comput. Math. Methods Med.* **2012**, *2012*, 673934. [CrossRef] [PubMed]

63. Krivoshein, A.V. From frame-like wavelets to wavelet frames keeping approximation properties and symmetry. *Appl. Math. Comput.* **2019**, *344/345*, 204–218. [CrossRef]

64. Doosti, H.; Iranmanesh, A.; Arashi, M.; Tabatabaey, S.M.M. On minimaxity of block thresholded wavelets under elliptical symmetry. *J. Stat. Plann. Inference* **2011**, *141*, 1526–1534. [CrossRef]

65. Rauhut, H. Wavelet transforms associated to group representations and functions invariant under symmetry groups. *Int. J. Wavelets Multiresolut. Inf. Process.* **2005**, *3*, 167–187. [CrossRef]

66. Liu, H. Wavelet transforms and symmetric tube domains. *J. Lie Theory* **1998**, *8*, 351–366.

67. Cohen, A.; Schlenker, J.M. Compactly supported bidimensional wavelet bases with hexagonal symmetry. *Constr. Approx.* **1993**, *9*, 209–236. [CrossRef]

68. Gai, Y.; Zhang, J. Detection the symmetry or asymmetry of model errors in partial linear models. *Comm. Stat. Simul. Comput.* **2022**, *51*, 2217–2234. [CrossRef]

69. Gai, Y.; Zhang, J. Detection of the symmetry of model errors for partial linear single-index models. *Comm. Stat. Simul. Comput.* **2022**, *51*, 3410–3427. [CrossRef]
70. Masry, E. Multivariate probability density estimation by wavelet methods: Strong consistency and rates for stationary time series. *Stoch. Process. Appl.* **1997**, *67*, 177–193. [CrossRef]

71. Masry, E. Wavelet-based estimation of multivariate regression functions in Besov spaces. *J. Nonparametr. Stat.* **2000**, *12*, 283–308. [CrossRef]

72. Masry, E. Probability density estimation from dependent observations using wavelets orthonormal bases. *Stat. Probab. Lett.* **1994**, *21*, 181–194. [CrossRef]

73. Chesneau, C.; Doosti, H. A note on the adaptive estimation of a conditional continuous-discrete multivariate density by wavelet methods. *Chin. J. Math.* **2016**, *2016*, 6204874. [CrossRef]

74. Chesneau, C.; Doosti, H.; Stone, L. Adaptive wavelet estimation of a function from an *m*-dependent process with possibly unbounded *m*. *Commun. Stat. Theory Methods* **2019**, *48*, 1123–1135. [CrossRef]

75. Bouzebda, S.; Didi, S. Multivariate wavelet density and regression estimators for stationary and ergodic discrete time processes: Asymptotic results. *Commun. Stat. Theory Methods* **2017**, *46*, 1367–1406. [CrossRef]

76. Bouzebda, S.; Didi, S.; El Hajj, L. Multivariate wavelet density and regression estimators for stationary and ergodic continuous time processes: Asymptotic results. *Math. Methods Stat.* **2015**, *24*, 163–199. [CrossRef]

77. Allaoui, S.; Bouzebda, S.; Chesneau, C.; Liu, J. Uniform almost sure convergence and asymptotic distribution of the wavelet-based estimators of partial derivatives of multivariate density function under weak dependence. *J. Nonparametr. Stat.* **2021**, *33*, 170–196. [CrossRef]

78. Allaoui, S.; Bouzebda, S.; Liu, J. Multivariate wavelet estimators for weakly dependent processes: Strong consistency rate. *Comm. Stat. Theory Methods* **2023**, *52*, 8317–8350. [CrossRef]

79. Allaoui, S.; Bouzebda, S.; Liu, J. Asymptotic Distribution of the Wavelet-based estimators of multivariate regression functions in Besov spaces under weak dependence. *J. Math. Inequal.* **2023**, *17*, 481–515. [CrossRef]

80. Chokri, K.; Bouzebda, S. Asymptotic normality for the wavelet partially linear additive model components estimation. *Comm. Stat. Theory Methods* **2024**, *53*, 8376–8411. [CrossRef]

81. Debbarh, M. Normalité asymptotique de l'estimateur par ondelettes des composantes d'un modèle additif de régression. *C. R. Math. Acad. Sci. Paris* **2006**, *343*, 601–606. [CrossRef]

82. Bouzebda, S.; Chokri, K.; Louani, D. Some uniform consistency results in the partially linear additive model components estimation. *Commun. Stat. Theory Methods* **2016**, *45*, 1278–1310. [CrossRef]

83. Doosti, H.; Niroumand, H.A. Multivariate stochastic regression estimation by wavelet methods for stationary time series. *Pak. J. Stat.* **2009**, *25*, 37–46.

84. Donoho, D.L.; Johnstone, I.M. Minimax estimation via wavelet shrinkage. *Ann. Stat.* **1998**, *26*, 879–921. [CrossRef]

85. Donoho, D.L.; Vetterli, M.; DeVore, R.A.; Daubechies, I. Data compression and harmonic analysis. *IEEE Trans. Inf. Theory* **1998**, *44*, 2435–2476. [CrossRef]

86. Birgé, L.; Massart, P. An adaptive compression algorithm in Besov spaces. *Constr. Approx.* **2000**, *16*, 1–36. [CrossRef]

87. Nickl, R.; Pötscher, B.M. Bracketing metric entropy rates and empirical central limit theorems for function classes of Besov- and Sobolev-type. *J. Theoret. Probab.* **2007**, *20*, 177–199. [CrossRef]

88. Nickl, R. Empirical and Gaussian processes on Besov classes. In *High Dimensional Probability*; IMS Lecture Notes Monograph Series; Institute of Mathematical Statistics: Beachwood, OH, USA, 2006; Volume 51, pp. 185–195. [CrossRef]

89. Schneider, C. *Beyond Sobolev and Besov—Regularity of Solutions of PDEs and Their Traces in Function Spaces*; Lecture Notes in Mathematics; Springer: Cham, Switzerland, 2021; Volume 2291, p. xviii+327. [CrossRef]

90. Sawano, Y. Theory of Besov spaces. In *Developments in Mathematics*; Springer: Singapore, 2018; Volume 56, p. xxiii+945. [CrossRef]

91. Peetre, J. *New Thoughts on Besov Spaces*; Duke University Mathematics Series, No. 1; Duke University, Mathematics Department: Durham, NC, USA, 1976; p. vi+305.

92. Rosenblatt, M. A central limit theorem and a strong mixing condition. *Proc. Natl. Acad. Sci. USA* **1956**, *42*, 43–47. [CrossRef]

93. Bradley, R.C. *Introduction to Strong Mixing Conditions*; Kendrick Press: Heber City, UT, USA, 2007; Volume 3, p. xii+597.

94. Gorodecki, V.V. The strong mixing property for linearly generated sequences. *Teor. Verojatnost. Primenen.* **1977**, *22*, 421–423.

95. Withers, C.S. Conditions for linear processes to be strong-mixing. *Z. Wahrsch. Verw. Geb.* **1981**, *57*, 477–480. [CrossRef]

96. Auestad, B.; Tjøstheim, D. Identification of nonlinear time series: First order characterization and order determination. *Biometrika* **1990**, *77*, 669–687. [CrossRef]

97. Chen, R.; Tsay, R.S. Functional-coefficient autoregressive models. *J. Am. Stat. Assoc.* **1993**, *88*, 298–308. [CrossRef]

98. Masry, E.; Tjøstheim, D. Nonparametric estimation and identification of nonlinear ARCH time series. *Econom. Theory* **1995**, *11*, 258–289. [CrossRef]

99. Masry, E.; Tjøstheim, D. Additive nonlinear ARX time series and projection estimates. *Econom. Theory* **1997**, *13*, 214–252. [CrossRef]

100. Newey, W.K. Kernel estimation of partial means and a general variance estimator. *Econom. Theory* **1994**, *10*, 233–253. [CrossRef]
101. Tjøstheim, D.; Auestad, B.r.H. Nonparametric identification of nonlinear time series: Projections. *J. Am. Stat. Assoc.* **1994**, *89*, 1398–1409.

102. Linton, O.; Nielsen, J.P. A kernel method of estimating structured nonparametric regression based on marginal integration. *Biometrika* **1995**, *82*, 93–100. [CrossRef]

103. Camlong-Viot, C.; Sarda, P.; Vieu, P. Additive time series: The kernel integration method. *Math. Methods Stat.* **2000**, *9*, 358–375.

104. Camlong-Viot, C.; Rodríguez-Póo, J.M.; Vieu, P. Nonparametric and semiparametric estimation of additive models with both discrete and continuous variables under dependence. In *The Art of Semiparametrics*; Contributions to Statistics; Physica-Verlag/Springer: Heidelberg, Germany, 2006; pp. 155–178. [CrossRef]

105. Hastie, T.; Tibshirani, R. Generalized additive models. *Stat. Sci.* **1986**, *1*, 297–318. [CrossRef]

106. Sperlich, S.; Linton, O.B.; Härdle, W. Integration and backfitting methods in additive models–Finite sample properties and comparison. *Test* **1999**, *8*, 419–458. [CrossRef]

107. Opsomer, J.D.; Ruppert, D. Fitting a bivariate additive model by local polynomial regression. *Ann. Stat.* **1997**, *25*, 186–211. [CrossRef]

108. Mammen, E.; Linton, O.; Nielsen, J. The existence and asymptotic properties of a backfitting projection algorithm under weak conditions. *Ann. Stat.* **1999**, *27*, 1443–1490. [CrossRef]

109. Fan, J.; Härdle, W.; Mammen, E. Direct estimation of low-dimensional components in additive models. *Ann. Stat.* **1998**, *26*, 943–971. [CrossRef]

110. Einmahl, U.; Mason, D.M. Uniform in bandwidth consistency of kernel-type function estimators. *Ann. Stat.* **2005**, *33*, 1380–1403. [CrossRef]

111. Deheuvels, P.; Mason, D.M. General asymptotic confidence bands based on kernel-type function estimators. *Stat. Inference Stoch. Process.* **2004**, *7*, 225–277. [CrossRef]

112. Stute, W. Conditional empirical processes. *Ann. Stat.* **1986**, *14*, 638–647. [CrossRef]

113. Park, H.; Stefanski, L.A. Relative-error prediction. *Stat. Probab. Lett.* **1998**, *40*, 227–236. [CrossRef]

114. Jones, M.C.; Park, H.; Shin, K.I.; Vines, S.K.; Jeong, S.O. Relative error prediction via kernel regression smoothers. *J. Stat. Plann. Inference* **2008**, *138*, 2887–2898. [CrossRef]

115. Demongeot, J.; Hamie, A.; Laksaci, A.; Rachdi, M. Relative-error prediction in nonparametric functional statistics: Theory and practice. *J. Multivar. Anal.* **2016**, *146*, 261–268. [CrossRef]

116. Bouzebda, S.; Taachouche, N. Oracle inequalities and upper bounds for kernel conditional *U*-statistics estimators on manifolds and more general metric spaces associated with operators. *Stochastics* **2024**, *96*, 2135–2198. [CrossRef]

117. Bellman, R. *Adaptive Control Processes: A Guided Tour*; Princeton University Press: Princeton, NJ, USA, 1961; p. xvi+255.

118. Scott, D.W.; Wand, M.P. Feasibility of multivariate density estimates. *Biometrika* **1991**, *78*, 197–205. [CrossRef]

119. Bouzebda, S. Weak convergence of the conditional single index *U*-statistics for locally stationary functional time series. *AIMS Math.* **2024**, *9*, 14807–14898. [CrossRef]

120. Attaoui, S.; Bentata, B.; Bouzebda, S.; Laksaci, A. The strong consistency and asymptotic normality of the kernel estimator type in functional single index model in presence of censored data. *AIMS Math.* **2024**, *9*, 7340–7371. [CrossRef]

121. McDonald, D.J.; Shalizi, C.R.; Schervish, M. Estimating beta-mixing coefficients via histograms. *Electron. J. Stat.* **2015**, *9*, 2855–2883. [CrossRef]

122. Khaleghi, A.; Lugosi, G. Inferring the mixing properties of an ergodic process. *arXiv* **2021**, arXiv:2106.07054.

123. Li, X.; Wang, L.; Nettleton, D. Additive partially linear models for ultra-high-dimensional regression. *Stat* **2019**, *8*, e223. [CrossRef]

124. Hoeffding, W. Probability inequalities for sums of bounded random variables. *J. Am. Stat. Assoc.* **1963**, *58*, 13–30. [CrossRef]

125. Camlong-Viot, C. Vers un test d'additivité en régression non paramétrique sous des conditions de mélange. *C. R. Acad. Sci. Paris Sér. I Math.* **2001**, *333*, 877–880. [CrossRef]

126. Fu, K.; Zhang, L. LIL behavior for *B*-valued strong mixing random variables. *Sci. China Math.* **2011**, *54*, 785–792. [CrossRef]

127. Bradley, R.C. Approximation theorems for strongly mixing random variables. *Mich. Math. J.* **1983**, *30*, 69–81. [CrossRef]

128. Tran, L.T. The $L_1$ convergence of kernel density estimates under dependence. *Canad. J. Stat.* **1989**, *17*, 197–208. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.