# Wavelet Estimation of Partial Derivatives in Multivariate Regression Under Discrete-Time Stationary Ergodic Processes

Sultana Didi, Salim Bouzebda

► To cite this version:

Sultana Didi, Salim Bouzebda. Wavelet Estimation of Partial Derivatives in Multivariate Regression Under Discrete-Time Stationary Ergodic Processes. Mathematics, 2025, 13 (10), pp.1587-10.3390/math13101587. [https://hal.science/hal-05064958v1](https://hal.science/hal-05064958v1)

HAL Id: hal-05064958

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
Article

# Wavelet Estimation of Partial Derivatives in Multivariate Regression Under Discrete-Time Stationary Ergodic Processes

Sultana Didi <fnref n="1" />,† ID and Salim Bouzebda <fnref n="2" />,*,† ID

**Abstract:** This study introduces a wavelet-based framework for estimating derivatives of a general regression function within discrete-time, stationary ergodic processes. The analysis focuses on deriving the integrated mean squared error (IMSE) over compact subsets of $\mathbb{R}^d$, while also establishing rates of uniform convergence and the asymptotic normality of the proposed estimators. To investigate their asymptotic behavior, we adopt a martingale-based approach specifically adapted to the ergodic nature of the data-generating process. Importantly, the framework imposes no structural assumptions beyond ergodicity, thereby circumventing restrictive dependence conditions. By establishing the limiting behavior of the wavelet estimators under these minimal assumptions, the results extend existing findings for independent data and highlight the flexibility of wavelet methods in more general stochastic settings.

**Keywords:** regression estimation; stationarity; ergodicity; rates of strong convergence; wavelet-based estimators; martingale differences; discrete time; stochastic processes; time series

**MSC:** 62G07; 62G08; 62G05; 62G20; 62H05; 60G42; 60G46

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

check for
updates

Academic Editors: Hamdy F.F. Mahmoud, Jiancheng Jiang and Stefano Bonnini

Received: 11 April 2025

Revised: 6 May 2025

Accepted: 9 May 2025

Published: 12 May 2025

**Citation:** Didi, S.; Bouzebda, S. Wavelet Estimation of Partial Derivatives in Multivariate Regression Under Discrete-Time Stationary Ergodic Processes. Mathematics **2025**, *13*, 1587. https://doi.org/10.3390/math13101587

**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

Nonparametric estimation has attracted sustained attention over many years, leading to a broad spectrum of methodological developments. Given its extensive range of applications and significant role within mathematical statistics, the estimation of both density and regression functions has become a central research theme. Among the most frequently employed approaches are kernel-type estimators, recognized for their flexibility and robust performance. For comprehensive treatments of these methods and their applications, the reader is referred to [1,2], and the references therein. Estimating derivatives of a function—whether a density or a regression function—serves as a powerful technique in statistical data analysis. Despite their importance, however, the nonparametric estimation of higher-order density derivatives has not been extensively explored. A key objective of this work is thus to investigate wavelet-based nonparametric estimators for the partial derivatives of multivariate densities. Derivative estimation arises in many disciplines, including economics and industry, where complex systems must often be modeled under limited prior knowledge. For instance, the second-order derivative of a density can support statistical tests for mode detection [3] and also guides bandwidth selection in kernel density

<footnote n="1">¹ Department of Statistics and Operations Research, College of Sciences, Qassim University, P.O. Box 6688, Buraydah 51452, Saudi Arabia; s.biha@qu.edu.sa</footnote>

<footnote n="2">² Université de Technologie de Compiègne, LMAC (Laboratory of Applied Mathematics of Compiègne), CS 60 319-60 203 Compiègne, France</footnote>

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

<footnote>† These authors contributed equally to this work.</footnote>
estimation [4]. In nonparametric signal estimation, the logarithmic derivative of the den-
sity—defined by the ratio of the derivative to the density itself—is instrumental for optimal
filtering and interpolation procedures [5], thereby making precise estimation crucial for
accurate signal processing. Moreover, gradient estimation is essential for filament detection
in point cloud data, a task of broad relevance in medical imaging, remote sensing, seis-
mology, and cosmology [6]. Additional motivations and applications appear in regression
analysis, Fisher information estimation, parameter estimation, and hypothesis testing [7].
Seminal investigations into density derivative estimation include [8–10], among others.
Further context for related statistical challenges, such as regression, Fisher information
estimation, parameter estimation, and hypothesis testing, may be found in [7].

It is likewise well documented that estimating first- and higher-order regression deriva-
tives holds substantial practical importance. Examples include modeling human growth
processes [11], assessing kidney function in lupus nephritis patients [12], and analyzing
Raman spectra of bulk materials [13]. In nonparametric regression, derivative estimation is
indispensable for constructing confidence intervals [14], selecting kernel bandwidths [15],
and comparing regression curves [16]. In the setting of a homoscedastic regression model,
ref. [17] proposed kernel *M*-estimators to estimate the first derivative of the regression
function in a nonparametric fashion, and extended these ideas heuristically to higher-order
derivatives. Derivative information also underpins modal regression, an alternative to
traditional regression methods for investigating the relationship between a response vari-
able *Y* and a predictor variable *X*; see [18]. Foundational studies of regression function
estimation can be found in [19,20]. In contrast, fewer works address the estimation of
derivatives for stationary densities or regression functions, with most established results
pertaining to independent and identically distributed data. Wavelet-based techniques have
recently attracted growing interest in statistical estimation, largely owing to their adaptabil-
ity to varying function regularity and their capacity to handle discontinuities effectively.
Furthermore, wavelet methods often lead to computationally efficient algorithms that
require relatively modest memory. For a detailed overview of wavelet methodologies in
nonparametric functional estimation, see [21]. Examples of wavelet applications include es-
timating the integrated squared derivative of a univariate density for independent data [22]
and for negatively or positively associated sequences [23]. The work of [24] extended these
methods to partial derivatives of multivariate densities under independence, while [25,26]
addressed the mixing scenario. More recently, ref. [27] studied wavelet estimators for the
partial derivatives of multivariate densities under additive noise.

To the best of our knowledge, methods for wavelet-based estimation of partial deriva-
tives of multivariate densities have not yet been extended to more general dependence
structures beyond the strong mixing framework. Addressing this gap forms the central mo-
tivation of our work. In particular, we draw on a collection of martingale-based techniques
that markedly differ from the tools typically employed under strong mixing conditions.
Nevertheless, as the subsequent sections will illustrate, bridging this gap entails much more
than merely combining pre-existing methodologies: it requires advanced mathematical
developments tailored to wavelet-based estimation in an ergodic setting.

The remainder of this article is organised as follows. Section 2 reviews the necessary
mathematical foundations and introduces the proposed class of linear wavelet estima-
tors. Section 3 specifies the underlying assumptions and presents the principal theoretical
results—uniform convergence rates and asymptotic normality established under weak-
dependence conditions. Section 4 applies the methodology to regression-derivative estima-
tion, whereas Section 5 explores its extension to mode regression. Concluding remarks and
avenues for further investigation are given in Section 6. To maintain the flow of exposition,
all proofs are collected in Section 7.
## Notation

Unless otherwise specified, $C$ denotes a positive constant whose value may change from line to line. We write $1_A(\cdot)$ for the indicator function of a set $A$. For two sequences of positive real numbers $\{a_n\}_{n \ge 1}$ and $\{b_n\}_{n \ge 1}$, we use the Landau notation

$$a_n = O(b_n)$$

when there exists a constant $C > 0$ such that $a_n \le C b_n$ for all sufficiently large $n$, and

$$a_n = o(b_n)$$

when $\lim_{n \to \infty} \frac{a_n}{b_n} = 0$.

# 2. Mathematical Backgrounds

## 2.1. Besov Spaces

Throughout this study we adopt the Besov scale $\mathbf{B}_{s,p,q}$ with parameters $s > 0$ and $1 \le p, q \le \infty$. Meyer's wavelet characterization [28] shows that, for $0 < s < r$ and any $f \in L_p(\mathbb{R}^d)$, the inclusion $f \in \mathbf{B}_{s,p,q}$ is equivalent to the finiteness of either of the two seminorms

$$ (B.1) \qquad J_{s,p,q}(f) = \|P_{V_0}f\|_{L_p} + \left(\sum_{j>0} (2^{js} \|P_{W_j}f\|_{L_p})^q\right)^{1/q} < \infty, $$

$$ (B.2) \qquad J'_{s,p,q}(f) = \|a_{0,\cdot}\|_{l_p} + \left(\sum_{j>0} (2^{j[s+d(\frac{1}{2}-\frac{1}{p})]} \|b_{j,\cdot}\|_{l_p})^q\right)^{1/q} < \infty, $$

where the scaling and wavelet coefficients are defined by

$$ a_{0,\mathbf{k}} = \int_{\mathbb{R}^d} f(\mathbf{u}) \phi_{0,\mathbf{k}}(\mathbf{u}) \, d\mathbf{u}, \qquad b_{i,j,\mathbf{k}} = \int_{\mathbb{R}^d} f(\mathbf{u}) \psi_{i,j,\mathbf{k}}(\mathbf{u}) \, d\mathbf{u}, $$

and

$$ \|b_{j,\cdot}\|_{l_p} = \left(\sum_{i=1}^{d} \sum_{\mathbf{k} \in \mathbb{Z}^d} |b_{i,j,\mathbf{k}}|^p\right)^{1/p}. $$

When $q = \infty$, the summations in (B.1) and (B.2) are replaced by the essential supremum. The Besov framework subsumes many classical smoothness classes that are ubiquitous in statistics and machine learning—for instance, the Sobolev (Hilbert) space $H^s$ corresponds to $\mathbf{B}_{s,2,2}$, while the (non-integer) Hölder-Lipschitz-Zygmund class $C^s$ coincides with $\mathbf{B}_{s,\infty,\infty}$. Additional equivalent formulations, together with their advantages in approximation theory and statistical methodology, are detailed in [29–32] and in Appendix A.

## 2.2. Linear Wavelets Estimator

This section opens with a concise review of the essentials of wavelet theory, following the notation of [28]. Let $\{V_j\}_{j \in \mathbb{Z}}$ be a multiresolution analysis of the Hilbert space $L_2(\mathbb{R}^d)$. Denote by $\phi$ the scaling function and by $\psi$ the associated orthogonal wavelet, both assumed to be $r$-regular ($r \ge 1$) and compactly supported in the hyper-cube $[-L, L]^d$ for some $L > 0$. For every integer $j$ and every index vector $\mathbf{k} \in \mathbb{Z}^d$, define

$$ \phi_{j,\mathbf{k}}(\mathbf{x}) = 2^{jd/2} \phi(2^j \mathbf{x} - \mathbf{k}). $$
The family {$\phi_{j,k}$}$_{k \in \mathbb{Z}^d}$ forms an orthonormal basis of $V_j$; moreover, each partial derivative of $\phi$ up to total order $r$ is rapidly decreasing. Specifically, for every integer $i > 0$ there exists a constant $A_i > 0$ such that (see [28], p. 29, Thm.2)

$$ |(\partial^\beta \phi)(\mathbf{x})| \le \frac{A_i}{(1 + \|\mathbf{x}\|)^i}, \quad \text{for all } |\beta| \le r. \quad (1) $$

In addition, there exist precisely $2^d - 1$ companion wavelets $\psi_1, \dots, \psi_{2^d - 1}$. Together they generate the collection

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

which constitutes an orthonormal basis of $L_2(\mathbb{R}^d)$.

We next recall the notions of strong mixing and ergodicity and the relationship between them. Let {$X_n$}$_{n \in \mathbb{Z}}$ be a strictly stationary sequence, and set

$$ \mathcal{A}_n = \sigma(X_k : k \le n) $$

and

$$ \mathcal{B}_m = \sigma(X_k : k \ge m). $$

The sequence is (a-)strongly mixing if

$$ \alpha(n) = \sup_{A \in \mathcal{A}_0, B \in \mathcal{B}_n} |\mathbb{P}(A \cap B) - \mathbb{P}(A)\mathbb{P}(B)| \to 0 \quad (n \to \infty). $$

It is ergodic when

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{k=0}^{n-1} |\mathbb{P}(A \cap \tau^{-k}B) - \mathbb{P}(A)\mathbb{P}(B)| = 0, $$

where $\tau$ denotes the left-shift operator. The foregoing definition of strong mixing is stricter than the one often used for measure-preserving transformations—namely, $\lim_{n \to \infty} \mathbb{P}(A \cap \tau^{-n}B) = \mathbb{P}(A)\mathbb{P}(B)$ for all measurable sets $A, B$ [33]. While strong mixing implies ergodicity, the converse need not hold (cf. Remark 2.6, p. 50, and Proposition 2.8, p. 51, of [34]). A number of authors have argued that an ergodic, rather than a strongly mixing, dependence structure is often preferable; see, for example, the discussion and illustrative examples in [35].

Finally, let $(\mathbf{X}, \mathbf{Y})$ be a random vector with $\mathbf{X} = (X_1, \dots, X_d) \in \mathbb{R}^d$ and $\mathbf{Y} = (Y_1, \dots, Y_q) \in \mathbb{R}^q$. Define the joint distribution function (df) of $(\mathbf{X}, \mathbf{Y})$ by

$$ F(\mathbf{x}, \mathbf{y}) := \mathbb{P}(\mathbf{X} \le \mathbf{x}, \mathbf{Y} \le \mathbf{y}), \quad \mathbf{x} \in \mathbb{R}^d, \mathbf{y} \in \mathbb{R}^q. $$

Henceforth, for vectors $\mathbf{v}' = (v'_1, \dots, v'_r) \in \mathbb{R}^r$ and $\mathbf{v}'' = (v''_1, \dots, v''_r) \in \mathbb{R}^r$, we write $\mathbf{v}' \le \mathbf{v}''$ if and only if $v'_j \le v''_j$ for all $j = 1, \dots, r$. We let $\mathbf{I}$ and $\mathbf{J}$ be two fixed subsets of $\mathbb{R}^d$ such that

$$ \mathbf{I} = \prod_{j=1}^{d} [a_j, b_j] \subset \mathbf{J} = \prod_{j=1}^{d} [c_j, d_j] \subset \mathbb{R}^d, $$

with

$$ -\infty < c_j < a_j < b_j < d_j < \infty \quad \text{for each } j = 1, \dots, d. $$
Suppose $(\mathbf{X}, \mathbf{Y})$ has a joint density function

$$f_{\mathbf{X},\mathbf{Y}}(\mathbf{x}, \mathbf{y}) := \frac{\partial^{d+q}}{\partial x_1 \dots \partial x_d \partial y_1 \dots \partial y_q} F(\mathbf{x}, \mathbf{y}) \quad \text{on} \quad \mathbf{J} \times \mathbb{R}^q,$$

with respect to the Lebesgue measure $d\mathbf{x}d\mathbf{y}$. Denote

$$f_{\mathbf{X}}(\mathbf{x}) = \int_{\mathbb{R}^q} f(\mathbf{x}, \mathbf{y}) d\mathbf{y}, \quad \mathbf{x} \in \mathbf{J},$$

as the marginal density of $\mathbf{X}$ on $\mathbf{J}$. Let $\rho: \mathbb{R}^q \to \mathbb{R}$ be a measurable function. We are chiefly interested in estimating derivatives of the form

$$\begin{aligned} (\partial^{\beta} r)(\rho; \mathbf{x}) &:= \partial^{\beta} [\mathbb{E}[\rho(\mathbf{Y}) | \mathbf{X} = \mathbf{x}]] f_{\mathbf{X}}(\mathbf{x}) \\ &= \int_{\mathbb{R}^q} \rho(\mathbf{y}) (\partial^{\beta} f)(\mathbf{x}, \mathbf{y}) d\mathbf{y}, \end{aligned}$$

where

$$(\partial^{\beta} f)(\mathbf{x}, \mathbf{y}) = \frac{\partial^{|\beta|} f(\mathbf{x}, \mathbf{y})}{\partial x_1^{\beta_1} \dots \partial x_d^{\beta_d}}, \quad \beta = (\beta_1, \dots, \beta_d), \quad |\beta| = \sum_{i=1}^d \beta_i.$$

When the regression function

$$m_{\rho}(\mathbf{x}) = \mathbb{E}[\rho(\mathbf{Y}) | \mathbf{X} = \mathbf{x}]$$

is sufficiently smooth, we may also consider its derivative

$$\partial_{\mathbf{x}}^{\beta} m_{\rho}(\mathbf{x}) = \frac{\partial^{|\beta|} m_{\rho}(\mathbf{x})}{\partial x_1^{\beta_1} \dots \partial x_d^{\beta_d}}.$$

Assume now that $(\partial^{\beta}r)(\rho; \cdot) \in L_2(\mathbb{R}^d)$. According to [36], for any integer $\tau$, the function $(\partial^{\beta}r)(\rho; \cdot)$ admits an expansion in the subspace $V_{\tau}$ of the multiresolution analysis:

$$(\partial^{\beta} r)(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} a_{\tau, \mathbf{k}} \phi_{\tau, \mathbf{k}}(\mathbf{x}), \quad (2)$$

where each wavelet coefficient $a_{\tau, \mathbf{k}}$ can be represented (via repeated integration by parts) as

$$\begin{aligned} a_{\tau, \mathbf{k}} &= \int_{\mathbb{R}^d} (\partial^{\beta} r)(\rho; \mathbf{u}) \phi_{\tau, \mathbf{k}}(\mathbf{u}) d\mathbf{u} \\ &= (-1)^{|\beta|} \int_{\mathbb{R}^d} r(\rho; \mathbf{u}) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{u}) d\mathbf{u}. \end{aligned}$$

Let $\partial^{\beta}\phi_{\tau, \mathbf{k}}(\cdot)$ denote the $\beta$-order partial derivative of $\phi_{\tau, \mathbf{k}}(\cdot)$. We define the linear estimator of $(\partial^{\beta}r)(\rho; \cdot)$ at the resolution level $\tau = m(n) \to \infty$ (the precise divergence rate of $m(n)$ will be specified below) by

$$(\widehat{\partial^{\beta} r})_n(\rho; \mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}_{\tau, \mathbf{k}} \phi_{\tau, \mathbf{k}}(\mathbf{x}), \quad (3)$$

where $\hat{a}_{\tau, \mathbf{k}}$ is the unbiased empirical estimator of $a_{\tau, \mathbf{k}}$:

$$\hat{a}_{\tau, \mathbf{k}} = \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \rho(\mathbf{Y}_i) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{X}_i). \quad (4)$$

**Remark 1.** Since $\phi(\cdot)$ and $\psi_i(\cdot)$ are bounded and compactly supported—where the support grows in a controlled manner with increasing differentiability—it follows that, for any fixed point $\mathbf{x}$, only
finitely many terms in the sums over $\mathbf{k} \in \mathbb{Z}^d$ contribute to the value of the wavelet expansions (see [36] for further details). This implies pointwise convergence of the expansions.

### 3. Assumptions and Main Results

To facilitate the presentation of our main results, we introduce the following notation. Let us denote by $g(\cdot)$ the density function of the random variable $Y$. Let $\mathcal{F}_n$ be the $\sigma$-field generated by $\{(\mathbf{X}_i, Y_i) : 1 \le i \le n\}$, let $f_{(\mathbf{X},Y)}^{\mathcal{F}_{i-1}}(\cdot)$ and $g_{(\mathbf{X},Y)}^{\mathcal{F}_{i-1}}(\cdot)$ and be the conditional densities of $(\mathbf{X}, Y)$ and $Y$ respectively, given the $\sigma$-field $\mathcal{F}_{i-1}$. Define the $\sigma$-field

$$ \mathcal{G}_n = \sigma((\mathbf{X}_k, Y_k); (\mathbf{X}_{n+1}) : 1 \le k \le n). $$

The following assumptions are imposed throughout the paper.

(C.1) For every $\mathbf{x} \in \mathbf{S}$, the sequence

$$ \frac{1}{n} \sum_{i=1}^{n} f_{(\mathbf{X},Y)}^{\mathcal{F}_{i-1}}(\mathbf{x}) $$

converges to $f(\mathbf{x})$ as $n \to \infty$, both almost surely (a.s.) and in the $L^2$ sense.

(C.2) Moreover, for every $\mathbf{x} \in \mathbf{S}$,

$$ \lim_{n \to \infty} \sup_{\mathbf{x} \in \mathbb{R}^d} \left| \frac{1}{n} \sum_{i=1}^{n} f_{(\mathbf{X},Y)}^{\mathcal{F}_{i-1}}(\mathbf{x}) - f(\mathbf{x}) \right| = 0, $$

again in both the almost sure and $L^2$ senses.

(C.3) (i) Assume that the multiresolution analysis is $r$-regular.

(ii) The density $(\partial^\beta r)(\rho; \cdot) \in B_{spq}$ for some $0 < s < r, 1 \le p, q \le \infty$.

(C.4) (i) The density $f(\cdot) \in B_{spq}$ for some $0 < s < r, 1 \le p, q \le \infty$.

(ii) The conditional density $f_{(\mathbf{X},Y)}^{\mathcal{F}_{i-1}}(\cdot) \in B_{spq}$ for some $0 < s < r, 1 \le p, q \le \infty$.

(N.0) $\mathbb{E}[|\rho(Y_1)|^\nu] < \infty$ for some $\nu > 2$.

(N.1) For any $\mathbf{y} \in \mathbb{R}$

$$ \lim_{n \to \infty} \sup_{\mathbf{y} \in \mathbb{R}} \left| \frac{1}{n g(\mathbf{y})} \sum_{i=1}^{n} g_{(\mathbf{X},Y)}^{\mathcal{F}_{i-1}}(\mathbf{y}) - 1 \right| = 0, \quad \text{in the a.s. and } L^2 \text{ sense.} $$

We may refer to [37] for further details.

(N.2) (i) The conditional mean of $Y_i$ given the $\sigma$-field $\mathcal{S}_i$ depends only on $\mathbf{X}_i$, i.e., for any $i \ge 1$,

$$ \mathbb{E}[\rho(Y_i) | \mathcal{G}_{i-1}] = \mathbb{E}[\rho(Y_i) | \mathbf{X}_i] = \mathbf{m}_{\rho}(\mathbf{x}). $$

(ii) for all $\eta \le 2$,

$$ \mathbb{E}[|\rho(Y_i)|^\eta | \mathcal{G}_{i-1}] = \mathbb{E}[|\rho(Y_i)|^\eta | \mathbf{X}_i] = \mathbf{m}_{\rho,\eta}(\mathbf{x}). $$

and the function $\mathbf{m}_{\rho,\eta}(\mathbf{x})$ is continuous at $\mathbf{x}$

(N.3) (i) We shall suppose that there exist two known constants $C_m > 0$ such that

$$ \sup_{\mathbf{x} \in \mathbf{S}} \mathbf{m}_{\rho}(\mathbf{x}) \le C_m. $$

(ii) The function regression satisfies a Hölder condition, that is, there exist constants $\beta > 0$ and $c_3 > 0$ such that, for any $(\mathbf{u}, \mathbf{v}) \in \mathbb{R}^{2d}$

$$ |\mathbf{m}_{\rho}(\mathbf{u}) - \mathbf{m}_{\rho}(\mathbf{v})| \le c_3 \|\mathbf{u} - \mathbf{v}\|^{\beta}. $$
(iii) The function regression $m_{\rho,\eta}$ satisfy a Hölder condition that is, there exist constants $\beta_1 > 0$ and $c_4 > 0$ such that, for any $(\mathbf{u}, \mathbf{v}) \in \mathbb{R}^{2d}$

$$|\mathbf{m}_{\rho,\eta}(\mathbf{u}) - \mathbf{m}_{\rho,\eta}(\mathbf{v})| \le c_4 \|\mathbf{u} - \mathbf{v}\|^{\beta_1}.$$

**Lemma 1.** Under condition (C.3), for $s > d/p$, we have

$$\sup_{\mathbf{x} \in \mathbb{R}^d} |\mathbb{E}[(\widehat{\partial^{\beta} r})_n(\rho; \mathbf{x})] - (\partial^{\beta} r)(\rho; \mathbf{x})| = \sup_{\mathbf{x} \in \mathbb{R}^d} |(\partial^{\beta} r)(\rho; \mathbf{x}) - (P_{V_m}(\partial^{\beta} r)(\rho; \mathbf{x}))(\mathbf{x})| \\ \le \text{const. } 2^{-(s-d/p)m(n)} J_{spq}((\partial^{\beta} r)(\rho; \cdot)).$$

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 \text{and} \quad h_n = 2^{-m(n)}. \qquad (5)$$

consequently

$$K^{(\beta)}(\mathbf{u}, \mathbf{v}) := \sum_{\mathbf{k} \in \mathbb{Z}^d} \phi(\mathbf{u} - \mathbf{k}) (\partial_{\mathbf{v}}^{\beta} \phi)(\mathbf{v} - \mathbf{k}). \qquad (6)$$

where $K^{(\beta)}(\mathbf{u}, \mathbf{v}) := (\partial_{\mathbf{v}}^{\beta} k)(\mathbf{u}, \mathbf{v})$ denotes the $\beta$-th partial derivative of $K(\mathbf{u}, \mathbf{v})$ with respect to $\mathbf{v}$. We infer that the derivative kernel function $K^{(\beta)}(\cdot, \cdot)$ defined in (6) converges uniformly and satisfies ([28], p. 33), i.e., for $|\alpha| \le r$, $|\beta| \le r$, and for some constant $C_m > 0$ with $m \ge 1$, the following inequalities holds:

$$| (\partial_{\mathbf{u}}^{\alpha} \partial_{\mathbf{v}}^{\beta} K)(\mathbf{u}, \mathbf{v}) | \le \frac{C_m}{(1 + \|\mathbf{v} - \mathbf{u}\|_2^m)^m},$$

where, for any $\mathbf{x} = (x_1, \dots, x_d) \in \mathbb{R}^d$,

$$\|\mathbf{x}\|_2 = \sqrt{\sum_{i=1}^{d} x_i^2}.$$

For $\alpha = 0$ and $m = d + |\beta| + 1$, we obtain

$$|K^{(\beta)}(\mathbf{u}, \mathbf{v})| \le \frac{C_{d+|\beta|+1}}{(1 + \|\mathbf{v} - \mathbf{u}\|_2^{d+|\beta|+1})}. \qquad (7)$$

That is, for any $j \ge 1$, we have

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

The function $G_j(d)$ is defined as:

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

where $\ell(t)$ represents the ell function, given by

$$\ell(t) := \int_{0}^{\infty} y^{t-1} e^{-y} dy.$$
Furthermore, by assuming $|\alpha| = 1$ and $m = 2$, we deduce that

$$ \left| \frac{\partial_{\mathbf{u}} \partial_{\mathbf{v}}^{\beta} K(\mathbf{u}, \mathbf{y})}{\partial u_i} \right| \le \frac{C_m}{(1 + \|\mathbf{u} - \mathbf{y}\|_2)^2} \le C_2, \quad i = 1, \dots, d. \quad (8) $$

This, in turn, implies that

$$ \left| K^{\beta}(\mathbf{u}, \mathbf{y}) - K^{\beta}(\mathbf{v}, \mathbf{y}) \right| \le C_2 \sum_{i=1}^{d} |u_i - v_i| \le d^{1/2} C_2 \|\mathbf{u} - \mathbf{v}\|_2. \quad (9) $$

By incorporating Equations (4), (3), and (6), the estimator $\widehat{(\partial^{\beta}r)}_n(\mathbf{x})$ can be expressed in an extended kernel estimation framework as follows:

$$ \widehat{(\partial^{\beta}r)}_n(\mathbf{x}) = \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^{n} \rho(Y_i) K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right), \quad \text{where } h_n = 2^{-m(n)}. \quad (10) $$

We refer the reader to [24] for further details.

**Theorem 1.** Under the stated assumptions (C.1), (C.4)(i), (N.0) and (N.2)(i), let $r$ be an element of the Besov space $B_{s,p,q}$ with $s > 1/p$, and $1 \le p, q \le \infty$. In this setting, the linear wavelet estimator $\widehat{(\partial^{\beta}r)}_n$ satisfies

$$ \mathbb{E} \left\| \widehat{(\partial^{\beta}r)}_n(\rho; \mathbf{x}) - (\partial^{\beta}r)(\rho; \mathbf{x}) \right\|_2^2 = O\left(n^{-\frac{2(s-|\beta|-d/4)}{2s+1}}\right). $$

**Theorem 2.** Assume that assumptions (C.2), (N.1), (N.2), (N.3)(i) and

$$ m(n) = m \to \infty \quad \text{and} \quad \frac{2^{dm(n)} \log n}{n} \to 0 \quad \text{as} \quad n \to \infty. $$

Define $T_n$ and assume that for $m(n) \to \infty$ slowly enough such that

$$ \left( \frac{2^{m(n)d} \log n}{n} \right)^{1/2} T_n = \left( \frac{\log n}{nh_n^d} \right)^{1/2} T_n \to 0, \quad \text{as} \quad n \to \infty. \quad (11) $$

For every compact subset $D \subset \mathbb{R}^d$, under assumptions (C.1) and (C.3), we have almost surely,

$$ \sup_{\mathbf{x} \in D} \left| \widehat{(\partial^{\beta}r)}_n(\rho; \mathbf{x}) - \mathbb{E}[\widehat{(\partial^{\beta}r)}_n(\rho; \mathbf{x})] \right| = O\left( \left( \frac{\log n}{nh_n^{d+2|\beta|}} \right)^{1/2} \right) + O(h_n^{-1/2}). $$

**Remark 2.** The factor $2^{-m(n)}$ that appears in the preceding theorems corresponds directly to the bandwidth $h_n$ in Parzen-Rosenblatt kernel density estimation. In practice, however, selecting the multiresolution level $m(n)$ within the wavelet framework is typically much simpler than choosing $h_n$. Because only a small, discrete set of candidate levels—usually three or four—needs to be examined, the procedure remains both conceptually transparent and computationally inexpensive. Adaptive choice of $m(n)$ most commonly relies on Stein's unbiased risk estimate, the classical rule of thumb, and cross-validation. Comprehensive treatments of these methods, as well as their use in establishing asymptotically optimal data-driven bandwidth selection rules, are provided by [38] and [21]. For the univariate setting ($d = 1$), the cross-validation criterion at resolution level $j$ is given by

$$ CV(j) = \sum_k \left[ \frac{2}{n(n-1)} \sum_{i=1}^{n} (\phi_{j,k}(X_i))^2 - \frac{n+1}{n^2(n-1)} \left( \sum_{i=1}^{n} \phi_{j,k}(X_i) \right)^2 \right]. $$
The statistic CV(j) depends solely on the observations X₁, ..., Xₙ and the resolution index j. The optimal level is selected by

$$j_0 = \arg \min_j \mathrm{CV}(j).$$

**Remark 3.** Kernel estimators suffer a well-documented decline in accuracy as the dimension of the covariate space grows—a manifestation of the curse of dimensionality [39]. In high dimensions, the volume of a local neighbourhood expands so rapidly that an impractically large sample is required to collect even a modest number of observations inside it. Unless the sample size is very large, one must therefore adopt bandwidths so wide that the resulting estimate is no longer truly local. A thorough exposition, complete with numerical illustrations, appears in [40], while more recent analyses are given in [41,42]. Although penalised splines enjoy considerable empirical popularity, their asymptotic behaviour is still not fully resolved even in classical non-parametric settings, and only a handful of theoretical investigations address the issue. In addition, most functional-regression techniques minimise an $L_2$ criterion and thus remain sensitive to outliers. Promising, if less conventional, alternatives include methods based on delta sequences and wavelet representations.

### 3.1. Asymptotic Normality Results

This section establishes a central limit theorem for the estimator defined in Equation (3). This problem was considered in several papers. In [10,43], the authors analyzed a non-parametric regression model with repeated measurements corrupted by mixing errors and, under mild moment and mixing-coefficient assumptions, established the r-th mean consistency, almost-sure consistency, complete consistency, optimal almost-sure convergence rates, and asymptotic normality of the associated wavelet estimator. For deconvolution density estimation with moderately ill-posed noise, ref. [44] proves the asymptotic normality of wavelet estimators when the target density lies in suitable Besov spaces. The study in [45] revisits non-parametric regression, deriving asymptotic normality of the wavelet estimator when the errors constitute an asymptotically negatively associated sequence. In a censored-data setting, ref. [46] treats density estimation where survival and censoring times form a stationary α-mixing sequence, again demonstrating asymptotic normality of the wavelet estimator. A fixed-design non-parametric regression model driven by strongly mixing errors is considered in [47], which likewise confirms asymptotic normality for the wavelet estimator of the regression function. Wavelet methods for linear processes—Gaussian or otherwise—with long, short, or negative memory are examined in [48]; both the log-regression and Whittle wavelet estimators are shown to be asymptotically normal, and explicit expressions for the limiting variance are derived via a general result for the suitably centred and normalised scalogram. Finally, ref. [49] proposes a wavelet estimator for fixed-design non-parametric regression with strictly stationary, associated errors, proving pointwise weak consistency and uniform asymptotic normality. For additional treatments of asymptotic normality of wavelet estimators in alternative frameworks, see [10,50,51]. Our results are obtained under mild regularity conditions on the estimator and only minimal bandwidth assumptions. We write

$$Z_n \xrightarrow{D} \mathcal{N}(0, \Sigma^2(\mathbf{x}))$$

to indicate that the sequence of random variables $(Z_n)_{n \ge 1}$ converges in distribution to a mean-zero normal distribution with covariance matrix $\Sigma^2(\mathbf{x})$.

**Theorem 3.** Assume that the hypotheses (C.1)–(C.4), (N.0), (N.2)(ii) and (N.3)(iii) are satisfied. Additionally, suppose that

$$nh_n^{d+2|\beta|+2\delta} \to 0 \quad \text{as} \quad n \to \infty. \qquad (12)$$
Then, the following convergence holds:

$$ \sqrt{n h_n^{d+2|\beta|}} (\widehat{(\partial^\beta r)}_n(\rho; \mathbf{x}) - (\partial^\beta r)(\rho; \mathbf{x})) \xrightarrow{D} \mathcal{N}(0, \Sigma_{(\beta)}^2(\mathbf{x})) $$

where

$$ \Sigma_{(\beta)}^2(\mathbf{x}) := \mathbf{m}_2(\rho, \mathbf{x}) f_X(\mathbf{x}) \int_{\mathbb{R}^d} (K^{(\beta)})^2(\mathbf{0}, \mathbf{u}) \, d\mathbf{u}, $$

and

$$ \mathbf{m}_2(\rho, \mathbf{x}) = \mathbb{E}[\rho^2(\mathbf{Y}) | \mathbf{X} = \mathbf{x}] = \frac{1}{f_X(\mathbf{x})} \int_{\mathbb{R}^d} \rho^2(\mathbf{y}) f_{\mathbf{X},\mathbf{Y}}(\mathbf{x}, \mathbf{y}) d\mathbf{y}. $$

The proof of Theorem 3 is detailed in Section 7.

**Remark 4.** Invoking Theorem 1, we obtain a refined bound for the mean-squared error of the wavelet-based estimator of multivariate density derivatives. Under the usual smoothness and moment assumptions,

$$ \mathbb{E} \| (\widehat{\partial^\beta f})_n - \partial^\beta f \|_2^2 = O(n^{-\frac{2(s-|\beta|)}{2s+1}}). $$

A parallel analysis, this time relying on Theorem 3, yields an asymptotic normality result. Specifically, for any fixed $\mathbf{x} \in \mathbb{R}^d$,

$$ \sqrt{n h_n^{d+2|\beta|}} (\widehat{(\partial^\beta f)}_n(\mathbf{x}) - \partial^\beta f(\mathbf{x})) \xrightarrow{D} \mathcal{N}(0, \Sigma_{(\beta)}^2(\mathbf{x})) $$

with asymptotic variance

$$ \Sigma_{(\beta)}^2(\mathbf{x}) = f_X(\mathbf{x}) \int_{\mathbb{R}^d} (K^{(\beta)})^2(\mathbf{0}, \mathbf{u}) \, d\mathbf{u}. $$

Both results are consistent with the theoretical developments in Theorem 4.8 of [52], further demonstrating the efficacy of wavelet methods for non-parametric estimation in multivariate contexts.

**Remark 5.** A noteworthy case of $(\partial^\beta r)(\rho, \cdot)$ arises when $\rho(\mathbf{y}) = 1_{\{\mathbf{y} \le \mathbf{t}\}}$. In this setting, from (3) we obtain

$$ \widehat{(\partial^\beta r)}_n(1_{\{\cdot \le \mathbf{t}\}}; \mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \tilde{a}_{\tau, \mathbf{k}} \phi_{\tau, \mathbf{k}}(\mathbf{x}), $$

where $\tilde{a}_{\tau, \mathbf{k}}$ is the unbiased estimator given by

$$ \tilde{a}_{\tau, \mathbf{k}} = \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} 1_{\{\mathbf{Y}_i \le \mathbf{t}\}} (\partial^\beta \phi_{\tau, \mathbf{k}})(\mathbf{X}_i). $$

Since this is a specific instance of the estimator in (3), Theorem 3 implies that

$$ \sqrt{n h_n^{d+2\beta}} (\widehat{(\partial^\beta r)}_n(1_{\{\cdot \le \mathbf{t}\}}; \mathbf{x}) - \partial^\beta r(1_{\{\cdot \le \mathbf{t}\}}; \mathbf{x})) \xrightarrow{D} \mathcal{N}(0, \Sigma^2(\mathbf{x})) $$

where

$$ \Sigma^2(\mathbf{x}) := m_2(1_{\{\cdot \le \mathbf{t}\}}, \mathbf{x}) f_X(\mathbf{x}) \int_{\mathbb{R}^d} (K^{(\beta)})^2(\mathbf{0}, \mathbf{u}) \, d\mathbf{u}. $$

This result is pivotal for determining the asymptotic behavior of wavelet estimators derived from the conditional distribution.

## 3.2. Confidence Interval

The asymptotic variance $\Sigma_{(\beta)}^2(\mathbf{x})$ that features in the central limit theorem depends on two unknown components—the conditional variance of $\mathbf{Y}$ given $\mathbf{X}$ and the density $f_X(\cdot)$ of $\mathbf{X}$—and therefore must be estimated in practice. We obtain consistent estimates
by employing a bounded, compactly supported wavelet basis drawn from the literature, such as the widely used Daubechies family [36]. Using a sufficiently large multiresolution level $\tau_n$ and an adaptively chosen initial level $j_0$, we recover the requisite nuisance parameters via wavelet-based estimation and subsequently insert them through a plug-in scheme. Using the estimators (3) we establish $\Sigma_{n,(\beta)}^2(\mathbf{x})$ a consistent estimate of the variance $\Sigma_{(\beta)}^2(\mathbf{x})$ given by

$$ \Sigma_{n,(\beta)}^2(\mathbf{x}) := \mathbf{m}_{n,2}(\rho, \mathbf{x}) f_{n,\mathbf{x}}(\mathbf{x}) \int_{\mathbb{R}^d} (K^{(\beta)})^2(\mathbf{0}, \mathbf{u}) \, d\mathbf{u}, \quad (13) $$

where $f_{n,\mathbf{x}}(\cdot)$ is the wavelet estimator of $f_{\mathbf{x}}(\cdot)$ and

$$ \mathbf{m}_{n,2}(\rho, \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}). \quad (14) $$

The considered coefficients estimators are given, respectively, for any $j_0 \le \tau$, by

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

and

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

The approximate confidence intervals of $(\partial^{\beta}r)(\rho; \mathbf{x})$ can be obtained as

$$ (\partial^{\beta}r)(\rho; \mathbf{x}) \in \left[ \widehat{(\partial^{\beta}r)}_n(\rho; \mathbf{x}) \pm c_{\alpha} \frac{\Sigma_{n,(\beta)}^2(\mathbf{x})}{\sqrt{n h_n^{d+2|\beta|}}} \right], $$

where $c_\alpha$ denotes the $(1 - \frac{\alpha}{2})$-quantile of the standard normal condition.

## 4. Application to the Regression Derivatives

In this section, we will follow the same notation as in [20]. We will consider especially the conditional expectation of $\rho(Y)$ given $X = x$, for $p = q = 1$. Recall that

$$ \begin{aligned} \mathbf{m}(x, \rho) &= \mathbb{E}(\rho(Y) | X = x) = \frac{1}{f_X(x)} \int_{\mathbb{R}} \rho(y) f_{X,Y}(x, y) dy \\ &= \frac{r(\rho, x)}{f_X(x)}. \end{aligned} $$

Recall the following derivatives

$$ \mathbf{m}_{\rho}^{\prime}(x) = \frac{r'(\rho, x)}{f_X(x)} - \frac{r(\rho, x) f_X'(x)}{f_X^2(x)}, \quad (15) $$

and

$$ \mathbf{m}_{\rho}''(x) = \frac{r''(\rho, x)}{f_X(x)} - \frac{2r'(\psi, x) f_X'(x)}{f_X^2(x)} + \frac{r(\rho, x)\{2(f_X'(x))^2 - f_X(x)f_X''(x)\}}{f_X^3(x)}. \quad (16) $$

In order to estimate the derivatives of $\mathbf{m}_{\rho}^{\prime}(x)$ in (15) and $\mathbf{m}_{\rho}''(x)$ in (16) by replacing $f_X(\cdot)$, $f_X'(\cdot)$, $f_X''(\cdot)$, $r(\rho; \cdot)$, $r'(\rho; \cdot)$ and $r''(\rho; \cdot)$ by $f_{X;n}(\cdot, h_n)$, $f_{X;n}'(\cdot, h_n)$, $f_{X;n}''(\cdot, h_n)$,
$r_n(\rho; \cdot, h_n)$, $r'_n(\rho; \cdot, h_n)$ and $r''_n(\rho; \cdot, h_n)$. We so define $\mathbf{m}'_{\rho,n}(x; h_n)$ and $\mathbf{m}''_{\rho,n}(x; h_n)$ when $f_{X;n}(x, h_n) \neq 0$. The definition of $\mathbf{m}'_{\rho,n}(x; h_n)$ and $\mathbf{m}''_{\rho,n}(x; h_n)$ is completed by setting $\mathbf{m}'_{\rho,n}(x; h_n) = \mathbf{m}''_{\rho,n}(x; h_n) = 0$ when $f_{X;n}(x, h_n) = 0$.

The following theorem is more or less a straightforward consequence of Theorem 2.

**Corollary 1.** Under the assumptions of Theorem 2, we have

$$ \sup_{x \in J} \left| \mathbf{m}'_{\rho,n}(x; h_n) - \mathbf{m}'_{\rho}(x) \right| = o_P(1), $$

$$ \sup_{x \in J} \left| \mathbf{m}''_{\rho,n}(x; h_n) - \mathbf{m}''_{\rho}(x) \right| = o_P(1). $$

## 5. Mode Regression

Now the location $\theta$ and size $\mathbf{m}_{\rho}(\theta)$ of the mode of $\mathbf{m}_{\rho}(\cdot)$ are estimated from the respective functionals $\hat{\theta}_n$ and $\hat{m}_{\rho,n}(\hat{\theta}_n)$ pertaining to $\hat{m}_{\rho,n}(\cdot)$, i.e., $\hat{\theta}_n$ is chosen through the equation

$$ \hat{m}_{\rho,n}(\hat{\theta}_n) = \max_{x \in I} \hat{m}_{\rho,n}(x). \qquad (17) $$

Then, a natural estimator of $\theta$ is

$$ \hat{\theta}_n := \arg \max_{x \in I} \hat{m}_{\rho,n}(x). \qquad (18) $$

Note that the estimate $\hat{\theta}_n$ is not necessarily unique and our results are valid for any chosen value satisfying (18). We point out that we can specify our choice by taking

$$ \hat{\theta}_n = \inf \left\{ y \in \mathbb{R} \text{ such that } \hat{m}_{\rho,n}(y) = \sup_{x \in \mathbb{R}} \hat{m}_{\rho,n}(x) \right\}. $$

It is known that kernel estimators tend to produce some additional and superfluous modality. However, this has no bearing on the asymptotic theory; our results are valid for any choice of $\hat{\theta}_n$ satisfying (18). Following the paper [53], we have the convergence result given in the following corollary.

**Corollary 2.** Under the assumptions of Theorem 2, we have

$$ |\hat{\theta}_n - \theta| = o_P(1). $$

**Remark 6.** Using weak-convergence techniques in the space of continuous functions, ref. [54] strengthened Parzen's distributional result under broadly weaker assumptions. He showed that an appropriately rescaled kernel-density estimator converges weakly to a randomly translated parabola that passes through the origin and possesses a fixed second derivative. In [55], the key insight is that, on a shrinking neighbourhood around the mode, the estimator $f$ converges weakly—under very general conditions—to a Gaussian process whose mean function is a parabola through the origin. This mean depends on the kernel moments and on the unknown density and its derivatives evaluated at the mode. More precisely, define

$$ Z_n(t) = h_n^{-2} [\hat{f}_n(\theta + h_n t) - \hat{f}_n(\theta)], \quad t \in \mathbb{R}, $$
where $\tilde{\theta}$ denotes the mode of $f$. The process $Z_n$ is thus a normalised version of the estimator in an interval centred at the mode. The limiting Gaussian process $Z(t)$ is given by

$$Z(t) = Y_{\alpha}(t) + (-1)^{p+1} \frac{d}{c p!} f^{(p+1)}(\tilde{\theta}) B_p t + \frac{1}{2} f^{(2)}(\tilde{\theta}) t^2, \quad t \in \mathbb{R},$$

where $c > 0, d \ge 0, p \ge 2$ is fixed,

$$B_p = \int x^p K(x) dx,$$

and $Y_\alpha$ is a mean-zero Gaussian process with covariance

$$R(s,t) = \frac{f(\tilde{\theta})}{c^{2+\alpha}} V_{\alpha}(s/t) s t^{1-\alpha}, \quad 0 \le \alpha \le 1, |s| \le |t|.$$

The function $V_\alpha$ is defined by

$$V_{\alpha}(\gamma) = \lim_{\delta \to 0} \delta^{\alpha} \int \left[ \frac{K(\gamma\delta - x) - K(-x)}{\gamma\delta} \right] \left[ \frac{K(\delta - x) - K(-x)}{\delta} \right] dx.$$

Ref. [55] establishes weak convergence of $\{Z_n\}$ to $Z$ for independent data. A central step in his Theorem 2.1 is to express $Z_n(t)$ as an average of i.i.d. variables, facilitating covariance computations and convergence of finite-dimensional distributions. Tightness is obtained via Theorem 15.7 of [56], which reduces the argument to evaluating the variance of independent variables (relation 2.9). The present study tackles the same problem via alternative methods, extending the framework of [18] to censored, dependent observations. Adapting Eddy's approach to this setting would be an attractive direction for future work. Such an extension, however, requires new probabilistic results analogous to those in [55] but tailored to dependent—e.g., mixing—samples, a significant undertaking that we leave for subsequent research.

**Remark 7.** We note that, when $|s| \ge 2$, $m_{\rho,n}^{(|s|)}(x, h_n) = D^{|s|}(m_{\rho,n}(x, h_n))$ may be obtained likewise through the usual Leibniz expansion of derivatives of products given by

$$m_{\rho,n}^{(|s|)}(x, h_n) = \sum_{|j|=0}^{|s|} C_{|s|}^{j} r_n^{(|j|)}(\rho; x, h_n) \{f_{X,n}^{-1}(x, h_n)\}^{(|s|-|j|)}, \quad f_{X,n}(x, h_n) \ne 0.$$

**Remark 8.** It is well known that the accuracy of kernel estimators deteriorates as the dimension of the predictor space grows. This phenomenon—popularly termed the curse of dimensionality [39]—stems from the fact that, in high dimensions, an impractically large number of observations is needed within any local neighborhood to secure a dependable estimate. When the available sample size is modest, practitioners are forced to adopt a bandwidth so wide that the very notion of “local” averaging effectively collapses. A thorough treatment of these difficulties, complemented by illustrative numerical experiments, can be found in [40]; see also the more recent contributions [41,42,57] for additional insights. Despite their popularity, penalized splines remain theoretically under-explored: rigorous results describing their asymptotic behaviour are sparse even in classical non-parametric settings. Moreover, a large proportion of functional-regression techniques are predicated on minimising the $L_2$ loss, rendering them particularly susceptible to outlying observations. More robust alternatives—such as methods based on delta sequences or on wavelet decompositions—offer promising directions for mitigating this sensitivity.

**Remark 9.** Let us recall some integral functions of the density function

$$T_1(F) = \int_{\mathbb{R}} (f'(x))^2 dx, \quad T_2(F) = - \int_{\mathbb{R}} a(F(x)) f'(x) dx$$
and

$$T_3(F) = \int_{\mathbb{R}} (f(x))^2 dx.$$

Notice that the functional $T_3(F)$ is a special case of $T_2(F)$. The functionals $T_1(F)$ and $T_3(F)$ appear in plug-in data-driven bandwidth selection procedures in density estimation (refer to [58] and the references therein) and the functional $T_2(F)$ arises as part of the variance in non-parametric location and regression estimation based on linear rank statistics (refer especially to [59]). Consider the following general class of integral functionals of the density:

$$T(F) = \int_{\mathbb{R}} \rho(x, F(x), F^{(1)}(x), \dots, F^{(r)}(x)) dF(x),$$

where $F$ is a cumulative distribution function on $\mathbb{R}$ with $r \ge 1$ derivatives $F^{(m)}$, refer to the [60,61] for more details. One can estimate $T(F)$ the plug-in methods making use of the wavelet estimate of the density function and its derivatives. The proof of such a statement, however, should require a different methodology than that used in the present paper, and we leave this problem open for future research.

**Remark 10.** The nonlinear thresholding framework furnishes an alternative class of estimators for the unknown density $f_X(x)$. Specifically,

$$\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}} \mathbf{1}\{ |\hat{b}_{ij\mathbf{k}}| > \delta_{j,n} \} \psi_{i,j,\mathbf{k}}(\mathbf{x}), \quad (19)$$

where $\delta_{j,n}$ denotes an appropriately selected threshold. In the univariate case ($d=1$), this estimator was first proposed by [29]. Extending its theoretical guarantees and empirical performance to higher-dimensional settings remains an intriguing direction for future research.

**Remark 11.** The sup-norm convergence rates proved in our theorems coincide with the optimal rates reported by [62–64]. The exact logarithmic factors depend on the resolution level, which is itself governed by the smoothness index $s$ of the target function $f$ in the Besov space $B_{s,p,q}$. Such a dependence on $s$ is a hallmark of non-parametric estimation and is widely documented in the literature. By allowing $f$ to possess general Besov regularity, our analysis relaxes the customary requirement of an integer-order derivative that underpins classical convolution-kernel methods, even though $s$ is typically unknown in practice.

**Remark 12.** To demonstrate that the assumptions stated in [65] are attainable, we provide three illustrative examples in which they are fulfilled:

1. Long-memory discrete-time processes: Let $(\epsilon_t)_{t \in \mathbb{Z}}$ be white noise with variance $\sigma^2$; denote the identity and back-shift operators by $I$ and $B$, respectively. According to [66] (Theorem 1, p. 55), the $k$-factor Gegenbauer process satisfies

$$\prod_{i \le i \le k} (I - 2v_i B + B^2)^{d_i} X_t = \epsilon_t,$$

with $0 < d_i < 1/2$ when $|v_i| < 1$ and $0 < d_i < 1/4$ when $|v_i| = 1$ for $i = 1, \dots, k$. This specification yields a stationary, causal, invertible series that exhibits long-range dependence. Moreover, it admits the moving-average representation

$$X_t = \sum_{j \ge 0} \psi_j(d, v) \epsilon_{t-j},$$
and the condition

$$ \sum_{j=0}^{\infty} \psi_j^2(d, v) < \infty $$

guarantees asymptotic stability. Nevertheless, ref. [67] shows that if $(\epsilon_t)_{t \in \mathbb{Z}}$ is Gaussian, the process is not strongly mixing. Even so, the moving-average form secures stationarity, Gaussianity, and ergodicity, clarifying the subtle influence of mixing conditions and emphasizing the interpretive value of the moving-average representation for long-memory dynamics.

2. Stationary solution of a linear Markov AR(1) process: Consider

$$ X_i = \frac{1}{2} X_{i-1} + \epsilon_i, $$

where $(\epsilon_i)$ are independent symmetric Bernoulli variables taking the values -1 and 1. As shown in [68], this model is not $\alpha$-mixing because of its dependence structure. It nevertheless remains stationary, Markov, and ergodic, illustrating that strong mixing is not necessary for either Markovianity or ergodicity—a point of direct relevance to statistical inference for time series and functional data.

3. A stationary process with an AR(1) representation: Let $(u_i)$ be an independent and identically distributed sequence uniformly distributed on $\{1, \dots, 9\}$, and define

$$ X_t := \sum_{i=0}^{\infty} 10^{-i-1} u_{t-i}, $$

so that $u_t, u_{t-1}, \dots$ constitute the decimal expansion of $X_t$. The series is stationary and can be written in AR(1) form:

$$ X_t = \frac{1}{10} X_{t-1} + \frac{1}{10} u_t = \frac{1}{10} X_{t-1} + \frac{1}{2} + \epsilon_t, $$

where

$$ \epsilon_t = \frac{1}{10} u_t - \frac{1}{2} $$

is strong white noise. Although it fails the $\alpha$-mixing criterion [69] (Example A.3, p. 349), the process is ergodic. This confirms that ergodicity may persist even in the absence of strong mixing, underscoring its suitability for non-parametric functional data analysis.

**Remark 13.** The paper already includes several well-known examples from the literature of processes that are ergodic yet not mixing. Broadening this to encompass more general linear Markov models—such as ARMA($p, q$) or ARIMA($p, d, q$) with $q \neq 0$—would certainly be of considerable interest. Nonetheless, pursuing this extension would necessitate a meticulous and comprehensive study to rigorously delineate the circumstances under which such processes remain ergodic while failing to satisfy mixing conditions.

## 6. Concluding Remarks

This study tackles the estimation of partial derivatives of multivariate regression functions. We introduce a family of non-parametric estimators based on linear wavelet methods and provide a rigorous theoretical analysis. Specifically, we establish strong uniform consistency on compact subsets of $\mathbb{R}^d$ and derive the corresponding convergence rates. We also prove the estimators' asymptotic normality and extend these results to ergodic processes, thereby broadening the scope of existing theory.

A key open problem is how to select smoothing parameters optimally so as to minimise the mean-squared error; addressing this question will be the focus of future work. Other promising directions include extending the methodology to functional ergodic data—a task that presents substantial mathematical challenges—and adapting it to settings with
incomplete observations, such as data missing at random or subject to various censoring mechanisms, particularly in spatially dependent contexts.

Further progress could come from relaxing the stationarity assumption to accommodate locally stationary processes and from developing comparable results for stationary continuous-time models. Both extensions would require fundamentally different analytical tools.

Finally, extensive numerical experiments on simulated and real datasets would enhance the practical relevance of our procedures, and the development of weighted bootstrap techniques—building on recent contributions such as [70,71]—offers another fruitful avenue for investigation.

## 7. Proofs

We derive an upper bound for the partial sums of unbounded martingale differences—a key ingredient in analyzing the asymptotic properties of the wavelet estimator based on strictly stationary, ergodic observations. Throughout the paper, $C$ denotes a generic positive constant that may change from one occurrence to the next. The necessary inequality is presented in the ensuing lemmas.

**Lemma 2 (Burkholder-Rosenthal inequality).** Following Notation 1 in [72].

Let $(X_i)_{i \ge 1}$ be a stationary martingale adapted to the filtration $(\mathcal{F}_i)_{i \ge 1}$, define $(d_i)_{i \ge 1}$ is the sequence of martingale differences adapted to $(\mathcal{F}_i)_{i \ge 1}$ and

$$S_n = \sum_{i=1}^{n} d_i,$$

then for any positive integer $n$,

$$ \| \max_{1 \le j \le n} |S_j| \|_p \ll n^{1/p} \|d_1\|_p + \left\| \sum_{k=1}^{n} \mathbb{E}(d_k^2 / \mathcal{F}_{k-1}) \right\|_{p/2}^{1/2}, \quad \text{for any } p \ge 2; \quad (20) $$

where as usual the norm

$$ \| \cdot \|_p = (\mathbb{E}[| \cdot |^p])^{1/p}. $$

**Lemma 3.** Let $(Z_n)_{n \ge 1}$ be a sequence of real martingale differences with respect to the sequence of $\sigma$-fields $(\mathcal{F}_n = \sigma(Z_1, \dots, Z_n))_{n \ge 1}$, where is the $\sigma$-field generated by the random variables $Z_1, \dots, Z_n$. Set

$$S_n = \sum_{i=1}^{n} Z_i.$$

For any $p \ge 2$ and any $n \ge 1$, assume that there exists some nonnegative constants $C$ and $d_n$ such that

$$ \mathbb{E}[Z_n^p | \mathcal{F}_{n-1}] \le C^{p-1} p! d_n^2, \quad \text{almost surely.} $$

Then, for any $\epsilon > 0$, we have

$$ \mathbb{P}(|S_n| > \epsilon) \le 2 \exp \left\{ - \frac{\epsilon^2}{2(D_n + C\epsilon)} \right\}, $$

where

$$ D_n = \sum_{i=1}^{n} d_i^2. $$

**Proof of Lemma 3.** The proof follows as a particular case of Theorem 8.2.2 due to [73]. □
To prove Theorem 1, we utilize the following two lemmas. Define the conditional expectation of $a_{\tau, \mathbf{k}}$

$$ \tilde{a}_{\tau, \mathbf{k}} = \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \mathbb{E} \left[ \rho(\mathbf{Y}_i) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) \mid \mathcal{F}_{i-1} \right]. \quad (21) $$

**Lemma 4.** For any $k \in \mathbb{Z}^d$, under the assumptions (C.1) and (N.2)(i), the following holds

$$ \tilde{a}_{\tau, \mathbf{k}} = a_{\tau, \mathbf{k}}, \quad \text{as} \quad n \to \infty. \quad (22) $$

**Proof of Lemma 4.** Observe that by invoking assumptions (N.2)(i) and (C.1), and the fact that $\mathcal{F}_{i-1} \subset \mathcal{G}_{i-1}$, we can further expand $\tilde{a}_{\tau, \mathbf{k}}$ as follows:

$$ \begin{aligned} \tilde{a}_{\tau, \mathbf{k}} &= \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \mathbb{E} \left[ \mathbb{E} \left[ \rho(\mathbf{Y}_i) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) \mid \mathcal{G}_{i-1} \right] \mid \mathcal{F}_{i-1} \right] \\ &= \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \mathbb{E} \left[ \mathbb{E} \left[ \rho(\mathbf{Y}_i) \mid \mathbf{X}_i \right] (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) \mid \mathcal{F}_{i-1} \right] \\ &= \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \int_{\mathbb{R}^d} m_{\rho}(\mathbf{x}) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{x}) f^{\mathcal{F}_{i-1}}(\mathbf{x}) d\mathbf{x} \\ &= (-1)^{|\beta|} \int_{\mathbb{R}^d} m_{\rho}(\mathbf{x}) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{x}) \left( \frac{1}{n} \sum_{i=1}^{n} f^{\mathcal{F}_{i-1}}(\mathbf{x}) \right) d\mathbf{x} \\ &= (-1)^{|\beta|} \int_{\mathbb{R}^d} m_{\rho}(\mathbf{x}) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{x}) (f(\mathbf{x}) + o(1)) d\mathbf{x} \\ &= (-1)^{|\beta|} \int_{\mathbb{R}^d} r(\rho, \mathbf{x}) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{x}) d\mathbf{x} + o(1) \\ &= a_{\tau, \mathbf{k}} + o(1). \end{aligned} \quad (23) $$

Consequently,

$$ \tilde{a}_{\tau, \mathbf{k}} = a_{\tau, \mathbf{k}} \quad \text{as} \quad n \to \infty, \quad (24) $$

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

**Lemma 5.** For any $k \in \mathbb{Z}^d$, under the assumption (N.2)(i), the following holds:

$$ \mathbb{E}[\hat{a}_{\tau, \mathbf{k}}] = a_{\tau, \mathbf{k}}, \quad \text{as} \quad n \to \infty \quad \text{a.s.} \quad (25) $$

**Proof of Lemma 5.** Following a similar argument as for (23) and using the property $\mathbb{E}[\mathbb{E}[\cdot|\mathcal{F}]] = \mathbb{E}[\cdot]$ where $\mathcal{F}$ is a sigma field, under the assumption (N.2)(i), we have

$$ \mathbb{E}[\hat{a}_{\tau, \mathbf{k}}] = \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \mathbb{E} \left[ \rho(\mathbf{Y}_i) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) \right] \quad (26) $$

$$ = \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \mathbb{E} \left[ \mathbb{E} \left[ \rho(\mathbf{Y}_i) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) \mid \mathbf{X}_i \right] \right] $$

$$ = \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \mathbb{E} \left[ m_{\rho}(\mathbf{X}_i) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) \right] $$

$$ = \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \int_{\mathbb{R}^d} m_{\rho}(\mathbf{x}) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{x}, \mathbf{y}) f(\mathbf{x}) d\mathbf{x} d\mathbf{y} $$

$$ = (-1)^{|\beta|} \int_{\mathbb{R}^d \times \mathbb{R}} \frac{r(\rho, \mathbf{x})}{f(\mathbf{x})} (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{x}) f(\mathbf{x}) d\mathbf{x} $$

$$ = (-1)^{|\beta|} \int_{\mathbb{R}^d \times \mathbb{R}} r(\rho, \mathbf{x}) (\partial^{\beta} \phi_{\tau, \mathbf{k}})(\mathbf{x}) d\mathbf{x} \quad (27) $$
which implies

$$ \mathbb{E}[\hat{a}_{\tau, \mathbf{k}}] = a_{\tau, \mathbf{k}}, \quad \text{as } n \to \infty. \quad \text{a.s.} \qquad (28) $$

This completes the proof of (25). □

**Lemma 6.** For any $k \in \mathbb{Z}^d$, under the assumptions (C.1), (C.4)(i), (N.0), and (N.2)(i), the following holds:

$$ \mathbb{E}[(\hat{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}})^2] = O\left(\frac{2^{2m(d/4+|\beta|)}}{n}\right), \quad \text{as } n \to \infty. \qquad (29) $$

**Proof of Lemma 6.** We present the following decomposition:

$$ \hat{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}} = (\hat{a}_{\tau, \mathbf{k}} - \tilde{a}_{\tau, \mathbf{k}}) + (\tilde{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}}). \qquad (30) $$

Therefore, we conclude that

$$ \begin{aligned} \mathbb{E}[(\hat{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}})^2] &= \mathbb{E}[(\hat{a}_{\tau, \mathbf{k}} - \tilde{a}_{\tau, \mathbf{k}})^2] + \mathbb{E}[(\tilde{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}})^2] \\ &\quad + 2\mathbb{E}[(\hat{a}_{\tau, \mathbf{k}} - \tilde{a}_{\tau, \mathbf{k}})(\tilde{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}})] \\ &= A_{m, \mathbf{k}, 1} + A_{m, \mathbf{k}, 2} + A_{m, \mathbf{k}, 3}. \end{aligned} \qquad (31) $$

By applying using Lemma 4 and its assumptions, especially statement (22), one obtains

$$ A_{m, \mathbf{k}, 2} = o(1). \qquad (32) $$

Moreover, using the same statement (22), we deduce

$$ \begin{aligned} A_{m, \mathbf{k}, 3} &= o(1) \mathbb{E}[\hat{a}_{\tau, \mathbf{k}} - \tilde{a}_{\tau, \mathbf{k}}] \\ &= o(1) \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \left( \mathbb{E}[\rho(Y_i)(\partial^\beta \phi_{\tau, \mathbf{k}})(\mathbf{X}_i)] - \mathbb{E}[\mathbb{E}[\rho(Y_i)(\partial^\beta \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) | \mathcal{F}_{i-1}]] \right) \\ &= 0. \end{aligned} \qquad (33) $$

We now direct our focus to the first term in the decomposition (31). Observe that

$$ \begin{aligned} \hat{a}_{\tau, \mathbf{k}} - \tilde{a}_{\tau, \mathbf{k}} &= \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} (\rho(Y_i)(\partial^\beta \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) - \mathbb{E}[\rho(Y_i)(\partial^\beta \phi_{\tau, \mathbf{k}})(\mathbf{X}_i) | \mathcal{F}_{i-1}]) \\ &= \frac{(-1)^{|\beta|}}{n} \sum_{i=1}^{n} \Phi_{\tau, \mathbf{k}}(\mathbf{X}_i, Y_i). \end{aligned} $$

Observe that $(\Phi_{\tau, \mathbf{k}}(\mathbf{X}_i))_{0 \le k \le m_n}$ forms a martingale difference sequence with respect to the filtration $(\mathcal{F}_i)_{0 \le k \le m_n}$. Applying Lemma 2, we directly conclude:

$$ A_{m, \mathbf{k}, 1} = \frac{(-1)^{2|\beta|}}{n^2} \mathbb{E}\left[\left|\sum_{i=1}^{n} \Phi_{\tau, \mathbf{k}}(\mathbf{X}_i, Y_i)\right|^2\right]. $$

Moreover, the following holds:

$$ \left( \mathbb{E} \left[ \left| \sum_{i=1}^{n} \Phi_{\tau, \mathbf{k}}(\mathbf{X}_i, Y_i) \right|^2 \right] \right)^{\frac{1}{2}} $$
$$ 
\begin{aligned}
&\leq \sqrt{n} \|\Phi_{\tau,k}(\mathbf{X}_1, Y_1)\|_2 + \left\|\sum_{i=1}^{n} \mathbb{E}[\Phi_{\tau,k}^2(\mathbf{X}_i, Y_i) | \mathcal{F}_{i-1}]\right\|_1^{\frac{1}{2}} \\
&= \Phi_{(1)} + \Phi_{(2)}. \qquad (34)
\end{aligned}
$$

To examine these components, we employ a standard decomposition approach and observe that $\mathcal{F}_0$ corresponds to the trivial $\sigma$-algebra. Consequently, the following result is derived:

$$ 
\begin{aligned}
\frac{1}{n}\Phi_{(1)}^2 &= \|\Phi_{\tau,k}(\mathbf{X}_1, Y_1)\|_2^2 \\
&= \mathbb{E}\left[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1) - \mathbb{E}[\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1) | \mathcal{F}_0]|^2\right] \\
&\leq \sum_{j=0}^2 \binom{2}{j} \mathbb{E}\left[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|^j \left(\mathbb{E}[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|]\right)^{2-j}\right] \\
&= \sum_{j=0}^2 \binom{2}{j} \mathbb{E}\left[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|^j\right] \left(\mathbb{E}[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|]\right)^{2-j} \\
&= \mathbb{E}\left[|\rho(\mathbf{Y}_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|^2\right] + 3 \left(\mathbb{E}[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|]\right)^2 \\
&\leq 4 \mathbb{E}\left[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|^2\right]. \qquad (35)
\end{aligned}
$$

By invoking the Cauchy-Schwarz inequality we obtain

$$ \mathbb{E}[|\rho(Y_1)(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|^2] \leq \mathbb{E}^{1/2}[|\rho(Y_1)|^4] \mathbb{E}^{1/2}[|(\partial^\beta\phi_{\tau,k})(\mathbf{X}_1)|^4] $$

Under assumption (N.0), there exists a constant $C_{\rho,4}$ such that

$$ \mathbb{E}[|\rho(Y_1)|^4] \leq C_{\rho,4} < \infty. $$

on the other hand, by employing a first-order Taylor expansion alongside Equation (1) and assumption (C.4)(i), we obtain

$$ 
\begin{aligned}
\mathbb{E}[|\partial^\beta \phi_{\tau,k}(\mathbf{X}_1)|^4] &= \int_{\mathbb{R}^d} (\partial^\beta \phi_{\tau,k})^4(\mathbf{u}) f(\mathbf{u}) \, d\mathbf{u} \\
&= 2^{m(2d+4|\beta|)} \int_{\mathbb{R}^d} (\partial^\beta \phi)^4 (2^m \mathbf{u} - \mathbf{k}) f(\mathbf{u}) \, d\mathbf{u} \\
&= 2^{m(d+4|\beta|)} \int_{\mathbb{R}^d} (\partial^\beta \phi)^4(\mathbf{v}) f\left(\frac{\mathbf{v}+\mathbf{k}}{2^m}\right) d\mathbf{v}. \\
&= 2^{m(d+4|\beta|)} \int_{\mathbb{R}^d} (\partial^\beta \phi)^4(\mathbf{v}) \left(f(\mathbf{v}) + O(2^{-md})\right) d\mathbf{v} \\
&\leq 2^{m(d+4|\beta|)} \left(\frac{A_i}{(1 + \|\mathbf{x}\|)^i}\right)^4 \\
&= O(2^{4m(d/4+|\beta|)}). \qquad (36)
\end{aligned}
$$

This leads to

$$ \Phi_{(1)} = O(2^{m(d/4+|\beta|)n^{1/2}}). \qquad (37) $$

Next, we analyze the second component of the decomposition in Equation (34). Specifically, we consider
$$
\begin{align*}
\Phi_{(2)} &= \left( \mathbb{E} \left( \sum_{i=1}^{n} \mathbb{E} \left[ \Phi_{\tau, \mathbf{k}}^2 (\mathbf{X}_i, Y_i) | \mathcal{F}_{i-1} \right] \right) \right)^{1/2} \\
&= \left( \sum_{i=1}^{n} \mathbb{E} \left( E \left[ \Phi_{\tau, \mathbf{k}}^2 (\mathbf{X}_i, Y_i) | \mathcal{F}_{i-1} \right] \right) \right)^{1/2} \\
&= \left( \sum_{i=1}^{n} \mathbb{E} \left[ \Phi_{\tau, \mathbf{k}}^2 (\mathbf{X}_i, Y_i) \right] \right)^{1/2}.
\end{align*}
$$

Applying a standard identity, we find

$$
\begin{align*}
\mathbb{E}[\Phi_{\tau,k}^2(\mathbf{X}_i, Y_i)] &= \mathbb{E}[\rho(Y_i) ((\partial^\beta \phi_{\tau,k})(\mathbf{X}_i) - \mathbb{E}[\rho(Y_i) (\partial^\beta \phi_{\tau,k})(\mathbf{X}_i) | \mathcal{F}_{i-1}])^2] \\
&\le 2\mathbb{E}[|\rho(Y_i) (\partial^\beta \phi_{\tau,k})(\mathbf{X}_i)|^2] + 2\mathbb{E}[\mathbb{E}[|\rho(Y_i) (\partial^\beta \phi_{\tau,k})(\mathbf{X}_i)|^2 | \mathcal{F}_{i-1}]] \\
&\le 4\mathbb{E}[|\rho(Y_i) (\partial^\beta \phi_{\tau,k})(\mathbf{X}_i)|^2].
\end{align*}
$$

Using Equation (36), it follows that

$$
\Phi_{(2)} = O\left(2^{m\left(\frac{d}{4} + \frac{|\beta|}{n}\right)} n^{\frac{1}{2}}\right). \quad (38)
$$

Combining the results from Equations (42) and (38), we derive

$$
\left( \mathbb{E} \left[ \left| \sum_{i=1}^{n} \Phi_{\tau, \mathbf{k}}(\mathbf{X}_i, Y_i) \right|^2 \right] \right)^{\frac{1}{2}} = O \left( 2^{m \left( \frac{d}{4} + \frac{|\beta|}{n} \right)} n^{\frac{1}{2}} \right).
$$

Consequently,

$$
\begin{align}
A_{m,\mathbf{k},1} &= \mathbb{E} \left[ (\hat{a}_{\tau,\mathbf{k}} - \tilde{a}_{\tau,\mathbf{k}})^2 \right] \nonumber \\
&= \frac{(-1)^{2|\beta|}}{n^2} \mathbb{E} \left[ \left| \sum_{i=1}^{n} \Phi_{\tau,\mathbf{k}}(\mathbf{X}_i, Y_i) \right|^2 \right] \nonumber \\
&= O\left(2^{2m\left(\frac{d}{4} + |\beta|\right)} n^{-1}\right). \tag{39}
\end{align}
$$

Finally, by integrating the findings from Equations (32), (33), and (39), we conclude that

$$
\mathbb{E}[(\hat{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}})^2] = O\left(\frac{2^{2m(d/4+|\beta|)}}{n}\right).
$$

This completes the proof. □

**Proof of Lemma 1.** The analysis of the bias term is purely analytical and follows from arguments analogous to those in [24], as it remains unaffected by the dependence structure. For brevity, we omit the detailed derivation. □

**Proof of Theorem 1.** Building on standard wavelet estimation techniques (see [21] for a detailed exposition), we derive the main result in the following manner. First, by applying the definition of projector normality, we obtain

$$
\mathbb{E} \| (\partial^{\beta} r)_n (\rho; \mathbf{x}) - (\partial^{\beta} r) (\rho; \mathbf{x}) \|_2^2
$$
$$ = \mathbb{E} \left\| \sum_{\mathbf{k} \in \mathbb{Z}^d} (\hat{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}}) \phi_{\tau, \mathbf{k}} \right\|_2^2 = \sum_{\mathbf{k} \in \mathbb{Z}^d} \mathbb{E} \left[ (\hat{a}_{\tau, \mathbf{k}} - a_{\tau, \mathbf{k}})^2 \right]. $$

Next, by invoking Lemma 6 and using the facts that $|V_\ell| \sim 2^m$ and $2^m \sim n^{\frac{1}{2s+1}}$, we deduce

$$ \begin{aligned} \mathbb{E} \left\| (\widehat{\partial^\beta r})_n (\rho; \mathbf{x}) - (\partial^\beta r) (\rho; \mathbf{x}) \right\|_2^2 &= \left( n^{\frac{1}{2s+1}} \right) \mathcal{O} \left( \frac{2^{2m(d/4+|\beta|)}}{n} \right) \\ &= \mathcal{O} \left( n^{-\frac{2(s-|\beta|-d/4)}{2s+1}} \right). \end{aligned} $$

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

**Proof of Theorem 2.** Now, we consider the following decomposition:

$$ \begin{aligned} & \sup_{\mathbf{x} \in \mathcal{D}} |(\widehat{\partial^\beta r})_n(\rho; \mathbf{x}) - \mathbb{E}[(\widehat{\partial^\beta r})_n(\rho; \mathbf{x})]| \\ & \le \sup_{\mathbf{x} \in \mathcal{D}} |(\widehat{\partial^\beta r})_n(\rho; \mathbf{x}) - \bar{r}_n(\rho; \mathbf{x})| + \sup_{\mathbf{x} \in \mathcal{D}} |\bar{r}_n(\rho; \mathbf{x}) - \mathbb{E}[(\widehat{\partial^\beta r})_n(\rho; \mathbf{x})]| \\ & = G_{n,1}(\rho; \mathbf{x}) + G_{n,2}(\rho; \mathbf{x}). \end{aligned} \quad (40) $$

To establish the uniform consistency of the above term we will follow the method used in ([74]). Let us introduce the truncated version of $(\widehat{\partial^\beta r})_n(\rho; \mathbf{x})$ as follows. Let

$$ \hat{a}_{\tau, \mathbf{k}}^T = \frac{1}{n} \sum_{i=1}^{n} \rho(\mathbf{Y}_i) \mathbf{1}\{|\rho(\mathbf{Y}_i)| \le T_n\} (\partial^\beta \phi_{\tau, \mathbf{k}})(\mathbf{X}_i), \quad (41) $$

and

$$ (\widehat{\partial^\beta r})_n^T (\rho; \mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} \hat{a}_{\tau, \mathbf{k}}^T \phi_{\tau, \mathbf{k}} (\mathbf{x}). \quad (42) $$

Throughout this work, the notation $\mathbf{1}\{A\}$ represents the indicator function for the set A. In a similar way as in the (10), we write $(\widehat{\partial^\beta r})_n^T(\rho; \mathbf{x})$ as an extended kernel estimator

$$ (\widehat{\partial^\beta r})_n^T (\rho; \mathbf{x}) = \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^{n} \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right). \quad (43) $$

Furthermore, we define:

$$ \bar{r}_n(\rho; \mathbf{x}) = \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^{n} \mathbb{E} \left[ \rho(Y_i) K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \mid \mathcal{F}_{i-1} \right]. $$

and

$$ \bar{r}_n^T(\mathbf{x}) = \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^{n} \mathbb{E} \left[ \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \mid \mathcal{F}_{i-1} \right]. \quad (44) $$

We begin by breaking down the initial term of the Equation (40), $G_{n,1}(\rho; \mathbf{x})$ into three distinct components, expressed as a sum through the following formulation:

$$ \begin{aligned} G_{n,1}(\rho; \mathbf{x}) &= \sup_{\mathbf{x} \in \mathcal{D}} |(\widehat{\partial^\beta r})_n(\rho; \mathbf{x}) - \bar{r}_n(\rho; \mathbf{x})| \\ &\le \sup_{\mathbf{x} \in \mathcal{D}} |(\widehat{\partial^\beta r})_n(\rho; \mathbf{x}) - (\widehat{\partial^\beta r})_n^T(\rho; \mathbf{x})| + \sup_{\mathbf{x} \in \mathcal{D}} |(\widehat{\partial^\beta r})_n^T(\rho; \mathbf{x}) - \bar{r}_n^T(\rho; \mathbf{x})| \\ &\quad + \sup_{\mathbf{x} \in \mathcal{D}} |\bar{r}_n^T(\rho; \mathbf{x}) - \bar{r}_n(\rho; \mathbf{x})| \end{aligned} $$
$$ = W_{n,1}(\rho; \mathbf{x}) + W_{n,2}(\rho; \mathbf{x}) + W_{n,3}(\rho; \mathbf{x}). \quad (45) $$

Recalling statement (7), we obtain readily that

$$ \sup_{\mathbf{x} \in D} |\widehat{(\partial^\beta \mathbf{r})}_n(\rho; \mathbf{x}) - (\widehat{\partial^\beta \mathbf{r}})_n^T(\rho; \mathbf{x})| \le \frac{(-1)^{|\beta|} h_n C_{d+|\beta|+1}}{n} \sum_{i=1}^n |\rho(Y_i) 1\{|\rho(Y_i)| > T_n\}|. \quad (46) $$

Markov inequality, in turn, implies that

$$ \mathbb{P}(|\rho(Y_n)| > T_n) \le T_n^{-\nu} \mathbb{E}(|\rho(Y_n)|^\nu), $$

under assumption (N.0) and using the fact that

$$ \sum_{n>1} T_n^{-\nu} < \infty, $$

By the Borel-Cantelli lemma,

$$ |\rho(Y_n)| \le T_n $$

almost surely for all sufficiently large $n$. Given the monotonicity of $T_n$, this implies $|\rho(Y_i)| \le T_n$ for all $i \le n$. Combining these results with the inequality (46) as in [74], we conclude that

$$ W_{n,1}(\rho; \mathbf{x}) = \sup_{\mathbf{x} \in D} |\widehat{(\partial^\beta \mathbf{r})}_n(\rho; \mathbf{x}) - (\widehat{\partial^\beta \mathbf{r}})_n^T(\rho; \mathbf{x})| = o(h_n), \quad \text{a.s.} \quad (47) $$

Once more, by (7), we infer that

$$ \begin{aligned} & \sup_{\mathbf{x} \in D} |\bar{r}_n^T(\rho; \mathbf{x}) - \bar{r}_n(\rho; \mathbf{x})| \\ & \le \frac{1}{nh_n^{d+|\beta|}} \sum_{i=1}^n |\mathbb{E}\left[\rho(Y_i) 1\{|\rho(Y_i)| \le T_n\} K^\beta\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) |\mathcal{F}_{i-1}\right] \\ & \phantom{\le} - \mathbb{E}\left[\rho(Y_i) K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) |\mathcal{F}_{i-1}\right]| \\ & \le \frac{1}{nh_n^{d+|\beta|}} \sum_{i=1}^n |\mathbb{E}\left[\rho(Y_i) 1\{|\rho(Y_i)| > T_n\} K^\beta\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) |\mathcal{F}_{i-1}\right] \\ & \le \frac{1}{n} \sum_{i=1}^n |\mathbb{E}\left[\rho(Y_i) 1\{|\rho(Y_i)| > T_n\} \frac{C_{d+|\beta|+1}}{h_n^{d+|\beta|} (1 + \|\mathbf{x} - \mathbf{X}_i\|_2 h_n^{-1})^{d+|\beta|+1}} |\mathcal{F}_{i-1}\right] \\ & \le \frac{h_n C_{d+|\beta|+1}}{n} \sum_{i=1}^n |\mathbb{E}[\rho(Y_i) 1\{|\rho(Y_i)| > T_n\} |\mathcal{F}_{i-1}]| \end{aligned} $$

Observe that the measurability of the functions $\rho(\cdot)$ and $v \mapsto v^q$, combined with the ergodicity of $(Y_n)_{n \in \mathbb{N}}$, ensures that the process $(g^{\mathcal{F}_{i-1}}(v))_{n \in \mathbb{N}}$ satisfies condition (N.1). Applying the Hölder and Markov inequalities, for any $\epsilon > 0$ and exponents $p, q$ fulfilling

$$ \frac{1}{p} + \frac{1}{q} = 1, $$

that

$$ \begin{aligned} & \mathbb{E}[|\rho(Y_i)| 1\{|\rho(Y_i)| > T_n\} |\mathcal{F}_{i-1}] \\ & \le (\mathbb{E}[|\rho(Y_i)|^q |\mathcal{F}_{i-1}])^{1/q} (\mathbb{P}\{|\rho(Y_i)| > T_n |\mathcal{F}_{i-1}\})^{1/p} \\ & \le T_n^{-q} \mathbb{E}[|\rho(Y_i)|^q |\mathcal{F}_{i-1}] \end{aligned} $$
$$ = T_n^{-q} \int_{\mathbb{R}} |\rho(v)|^q g^{\mathcal{F}_{i-1}}(v) dv. $$

This, in turn, readily implies that

$$ \begin{aligned} & \sup_{\mathbf{x} \in D} |\bar{r}_n^T(\rho; \mathbf{x}) - \bar{r}_n(\rho; \mathbf{x})| \\ & \le h_n C_{d+|\beta|+1} \sup_{y \in \mathbb{R}} \left\| \frac{1}{ng(y)} \sum_{i=1}^n g^{\mathcal{F}_{i-1}}(y) \right\| T_n^{-q} \int_{\mathbb{R}} |\rho(v)|^q g(v) dv \\ & = h_n \sup_{y \in \mathbb{R}} \left\| \frac{1}{ng(y)} \sum_{i=1}^n g^{\mathcal{F}_{i-1}}(y) \right\| T_n^{-q} \mathbb{E}[|\rho(Y_i)|^q], \end{aligned} $$

which gives

$$ W_{n,3}(\rho; \mathbf{x}) = \sup_{\mathbf{x} \in D} |\bar{r}_n^T(\rho; \mathbf{x}) - \bar{r}_n(\rho; \mathbf{x})| = O(h_n) \quad \text{a.s.} \quad (48) $$

Lastly, we turn our attention to second term of the decomposition (45), given the compactness of D, there exists a finite covering by L(n) cubes I$_{n,j}$, each centered at x$_{n,j}$ with side length ℓ$_{n}$, where j = 1, ..., L(n). It follows directly that

$$ \ell_n = \frac{\text{const}}{L^{1/d}(n)}. $$

where

$$ L(n) = \left( \frac{2^{(d+2)m(n)} n T_n^2}{\log n} \right)^{d/2} = \left( \frac{n T_n^2}{h_n^{d+2} \log n} \right)^{d/2}. $$

we readily infer that

$$ \begin{aligned} W_{n,2}(\rho; \mathbf{x}) & = \sup_{\mathbf{x} \in D} |\widehat{(\partial \beta r)}_n^T(\rho; \mathbf{x}) - \bar{r}_n^T(\rho; \mathbf{x})| \\ & = \max_{1 \le j \le L(n)} \sup_{\mathbf{x} \in D \cap I_j} |\widehat{(\partial \beta r)}_n^T(\rho; \mathbf{x}) - \bar{r}_n^T(\rho; \mathbf{x})| \\ & \le \max_{1 \le j \le L(n)} \sup_{\mathbf{x} \in D \cap I_j} |\widehat{(\partial \beta r)}_n^T(\rho; \mathbf{x}) - \widehat{(\partial \beta r)}_n^T(\rho; \mathbf{x}_j)| \\ & \quad + \max_{1 \le j \le L(n)} |\widehat{(\partial \beta r)}_n^T(\rho; \mathbf{x}_j) - \bar{r}_n^T(\rho; \mathbf{x}_j)| \\ & \quad + \max_{1 \le j \le L(n)} \sup_{\mathbf{x} \in D \cap I_j} |\bar{r}_n^T(\rho; \mathbf{x}_j) - \bar{r}_n^T(\rho; \mathbf{x})| \\ & = Q_1 + Q_2 + Q_3. \end{aligned} \quad (49) $$

The statement (9) allows us to infer that

$$ \begin{aligned} & |\widehat{(\partial \beta r)}_n^T(\rho; \mathbf{x}) - \widehat{(\partial \beta r)}_n^T(\rho; \mathbf{x}_j)| \\ & = \left| \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} \left( K\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) - K\left(\frac{\mathbf{x}_j}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right) \right| \\ & \le \frac{d^{1/2} C_2 T_n}{h_n^{d+|\beta|+1}} \|\mathbf{x} - \mathbf{x}_i\|, \end{aligned} $$
which implies that

$$
\begin{align}
Q_1 &= \max_{1 \le j \le L(n)} \sup_{\mathbf{x} \in D \cap I_j} \left| (\widehat{\partial^\beta r})_n^T (\rho; \mathbf{x}) - \bar{r}_n^T (\rho; \mathbf{x}) \right| \nonumber \\
&\le \frac{d^{1/2} C_2 T_n \ell_n}{h_n^{d+|\beta|+1}} = \frac{d^{1/2} C_2 T_n \text{const}}{L^{1/d}(n) h_n^{d+|\beta|+1}} \nonumber \\
&= O\left( \left( \frac{\log n}{nh_n^{d+2|\beta|}} \right)^{1/2} \right), \quad \text{a.s.} \tag{50}
\end{align}
$$

A similar argument shows, likewise, that

$$
\begin{align}
Q_3 &= \max_{1 \le j \le L(n)} \sup_{\mathbf{x} \in D \cap I_j} \left| \bar{r}_n^T(\rho; \mathbf{x}_j) - \bar{r}_n^T(\rho; \mathbf{x}) \right| \nonumber \\
&\le \frac{d^{1/2} C_2 T_n \ell_n}{h_n^{d+|\beta|+1}} = \frac{d^{1/2} C_2 T_n \text{const}}{L^{1/d}(n) h_n^{d+|\beta|+1}} \nonumber \\
&= O\left( \left( \frac{\log n}{nh_n^{d+2|\beta|}} \right)^{1/2} \right), \quad \text{a.s.} \tag{51}
\end{align}
$$

We now analyze the term $Q_2$ in the right side of (49) and show that

$$
Q_2 = O\left( \left( \frac{\log n}{nh_n^{d+2|\beta|}} \right)^{1/2} \right), \quad \text{a.s.} \tag{52}
$$

Observe that

$$
\begin{align*}
Q_2 &= \max_{1 \le j \le L(n)} \left| \widehat{(\partial^\beta r)}_n^T (\rho; \mathbf{x}_j) - \bar{r}_n^T (\rho; \mathbf{x}_j) \right| \\
&= \max_{1 \le j \le L(n)} \left| \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \left( \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}_i}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right. \right. \\
&\qquad \left. - \mathbb{E}\left[ \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}_i}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \middle| \mathcal{F}_{i-1} \right] \right| \\
&= \max_{1 \le j \le L(n)} \left| \frac{1}{nh_n^{d+|\beta|}} \sum_{i=1}^n Z(\rho, \mathbf{x}_j) \right|,
\end{align*}
$$

where, for $i = 1, \dots, n$,

$$
Z(\rho, \mathbf{x}_j) = (-1)^{|\beta|} \left( \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}_i}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right. \\
\left. - \mathbb{E}\left[ \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}_i}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \middle| \mathcal{F}_{i-1} \right] \right)
$$

forms a sequence of martingale-difference arrays relative to the $\sigma$-field $\mathcal{F}_i$. We now utilize
the Lemma 3 for the partial sums of unbounded of martingale differences to obtain an
upper bound. Observe that

$$
\begin{align*}
(Z(\rho, \mathbf{x}_j))^p \\
&= (-1)^{p|\beta|} \left( \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right. \\
&\qquad \left. - \mathbb{E}\left[ \rho(Y_i) \mathbf{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \middle| \mathcal{F}_{i-1} \right] \right)^p
\end{align*}
$$
$$ = (-1)^{p|\beta|} \sum_{k=0}^{p} (-1)^{p-k} C_p^k \left( \rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right)^k \left( \mathbb{E}\left[\rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) | \mathcal{F}_{i-1}\right] \right)^{p-k} . $$

By using the fact that

$$ \left( \mathbb{E}\left[\rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}_j}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) | \mathcal{F}_{i-1}\right] \right)^{p-k} $$

is $\mathcal{F}_{i-1}$-measurable, it follows then

$$ \begin{aligned} & |\mathbb{E}[Z(\rho, \mathbf{x}_j)^p | \mathcal{F}_{i-1}]| \\ & \le \sum_{k=0}^{p} C_p^k \mathbb{E} \left[ \left| \left( \rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}_j}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right)^k \right| |\mathcal{F}_{i-1}| \right] \\ & \times \mathbb{E} \left[ \left| \left( \rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}_j}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right)^{p-k} \right| |\mathcal{F}_{i-1}| \right]. \end{aligned} \quad (53) $$

recalling the fact that

$$ \mathcal{F}_{i-1} \subset \mathcal{G}_{i-1} $$

and applying the assumption (N.2)(ii), for any integer $\ell > 1$, we obtain

$$ \begin{aligned} & \mathbb{E}\left[\left|\rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right|^{\ell} | \mathcal{F}_{i-1}\right] \\ &= \mathbb{E}\left[\mathbb{E}\left[\left|\rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right|^{\ell} \middle| \mathcal{G}_{i-1}\right] | \mathcal{F}_{i-1}\right] \\ &\le \mathbb{E}\left[\mathbb{E}\left[\left|\rho(Y_i)\right|^{\ell} | \mathcal{G}_{i-1}\right] \middle| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right|^{\ell} | \mathcal{F}_{i-1}\right] \\ &= \mathbb{E}\left[m_{\rho, \ell}(\mathbf{X}_i) \middle| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right|^{\ell} | \mathcal{F}_{i-1}\right], \end{aligned} \quad (54) $$

where $m_{\rho,\ell}(\mathbf{X}_i)$ is a positive function as stipulated in assumption (N.2)(ii), ensures that

$$ \begin{aligned} & \mathbb{E}\left[\left|\rho(Y_i) \mathbb{1}\{|\rho(Y_i)| \le T_n\} K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right|^{\ell} | \mathcal{F}_{i-1}\right] \\ & \le \mathbb{E}\left[\left| m_{\rho,\ell}(\mathbf{X}_i) - m_{\rho,\ell}(\mathbf{x}) \right| \middle| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right|^{\ell} | \mathcal{F}_{i-1}\right] \\ & \quad + m_{\rho,\ell}(\mathbf{x}) \mathbb{E}\left[\left| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right|^{\ell} | \mathcal{F}_{i-1}\right] \\ &= \left[ \sup_{\|\mathbf{u}-\mathbf{x}\| \le h_n} \left| m_{\rho,\ell}(\mathbf{u}) - m_{\rho,\ell}(\mathbf{x}) \right| + m_{\rho,\ell}(\mathbf{x}) \right] \mathbb{E}\left[\left| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right|^{\ell} | \mathcal{F}_{i-1}\right] \\ &\le C_0 \mathbb{E}\left[\left| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right|^{\ell} | \mathcal{F}_{i-1}\right], \end{aligned} $$

where $C_0$ is a positive constant. It follows readily from (7),

$$ \mathbb{E}\left[\left| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right|^{\ell} | \mathcal{F}_{i-1}\right] \le \mathbb{E}\left[\left| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right|^{\ell} | \mathcal{F}_{i-1}\right] $$
$$
\begin{align*}
&= \int_{\mathbb{R}^d} \left| K^{(\boldsymbol{\beta})} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) \right|^{\ell} f^{\mathcal{F}_{i-1}}(\mathbf{y}) d\mathbf{y} \\
&\le \int_{\mathbb{R}^d} \left( \frac{C_{d+|\boldsymbol{\beta}|}}{(1 + h_n^{-1} \|\mathbf{x} - \mathbf{y}\|_2)^{d+|\boldsymbol{\beta}|}} \right)^{\ell} f^{\mathcal{F}_{i-1}}(\mathbf{y}) d\mathbf{y} \\
&\le h_n^{\ell(d+|\boldsymbol{\beta}|)} C_{d+|\boldsymbol{\beta}|}^{\ell} \tag{55}
\end{align*}
$$

statements (53), (54) and (55) give the following upper bound

$$
\begin{align*}
|\mathbb{E}[Z(\rho, \mathbf{x}_j)^p | \mathcal{F}_{i-1}]| &\le C_0^2 \sum_{k=0}^p C_p^k h_n^{p(d+|\boldsymbol{\beta}|)} C_{d+|\boldsymbol{\beta}|}^p \\
&\le C_0^2 2^p h_n^{d+|\boldsymbol{\beta}|} C_{d+|\boldsymbol{\beta}|}^p \le p! C^{p-1} d_i^2.
\end{align*}
$$

We apply Lemma 3 to the summation of the martingale differences $Z_n(\mathbf{x}_j, \mathbf{X}_i)$, where

$$
C = C_{d+|\boldsymbol{\beta}|}, \quad d_i^2 = C_0^2 C_{d+|\boldsymbol{\beta}|} |h_n^{d+|\boldsymbol{\beta}|}|, \quad D_n = \sum_{i=1}^n d_i^2 = O(nh_n^{d+|\boldsymbol{\beta}|})
$$

and

$$
\epsilon_n = \epsilon_0 (\log n / nh_n^{d+2|\beta|})^{1/2}.
$$

Consequently, there exists a positive constant $C_2 > 0$ such that the following inequalities hold:

$$
\begin{align*}
& \mathbb{P} \left( \max_{1 \le j \le L(n)} \left| \frac{1}{nh_n^{d+|\boldsymbol{\beta}|}} \sum_{i=1}^n Z(\rho, \mathbf{x}_j) \right| > \epsilon_n \right) \\
\le & \sum_{j=1}^{L(n)} \mathbb{P} \left( \left| \sum_{i=1}^n Z(\rho, \mathbf{x}_j) \right| > \epsilon_n nh_n^{d+|\boldsymbol{\beta}|} \right) \\
\le & 2L(n) \exp \left\{ - \frac{\epsilon_0^2 (nh_n^{d+|\boldsymbol{\beta}|})^2 (\log n / nh_n^{d+2|\boldsymbol{\beta}|})}{2(D_n + 2C_{d+1}nh_n^{d+|\boldsymbol{\beta}|} (\log n / nh_n^{d+2|\boldsymbol{\beta}|})^{1/2})} \right\} \\
\le & \left( \frac{nT_n^2}{h_n^{d+1} \log n} \right)^{d/2} \exp \left\{ - \frac{\epsilon_0^2 nh_n^{d+|\boldsymbol{\beta}|} (\log n / h_n^{|\boldsymbol{\beta}|})}{O(nh_n^{d+|\boldsymbol{\beta}|}) (1 + 2C_{d+1} (\log n / nh_n^{d+2|\boldsymbol{\beta}|})^{1/2})} \right\} \\
\le & \left( T_n \left( \frac{\log n}{nh_n^d} \right)^{1/2} \right)^d \left( \frac{1}{(nh_n)^{1/2} \log n} \right)^n \left( \frac{C_2 \epsilon_0^2}{h_n^{|\boldsymbol{\beta}|}} - \frac{3d}{2} \right).
\end{align*}
$$

Observe that

$$
nh_n \to \infty \text{ as } n \to \infty, \text{ and } \frac{1}{h_n^{|\boldsymbol{\beta}|}} \to \infty \text{ as } n \to \infty.
$$

this implies that when $n \to \infty$ we have

$$
\frac{C_2 \epsilon_0^2}{h_n^{|\boldsymbol{\beta}|}} - \frac{3d}{2} > 0,
$$

we readily obtain, under condition (11), that

$$
\sum_{n \ge 1} \mathbb{P} \left( \max_{1 \le j \le L(n)} \left| \frac{1}{nh_n^{d+|\boldsymbol{\beta}|}} \sum_{i=1}^n Z(\rho, \mathbf{x}_j) \right| > \epsilon_n \right) < \infty. \quad (56)
$$
The assertion (52) is established through a standard application of the Borel-Cantelli lemma. Consequently, the result stated in Theorem 2 follows by applying the decompositions (40) that resulted in two terms: the first one was rewritten in (45) as the sum of three terms, the statements (47) and (48) and the decomposition (49) in conjunction with estimates (50), (51), and (52). To complete the argument, we now analyze the second term on the right-hand side of (40). Employing the same reasoning as in (54), and in view of assumptions (N.2)(i) and (N.3)(i), we observe that

$$
\begin{align*}
& \sup_{x \in D} \left| \bar{r}_n(\rho; \mathbf{x}) - \mathbb{E}[(\widehat{\partial \beta r})_n(\rho; \mathbf{x})] \right| \\
&= \sup_{x \in D} \left| \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \left( \mathbb{E} \left[ \rho(Y_i) K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) | \mathcal{F}_{i-1} \right] - \mathbb{E} \left[ \rho(Y_i) K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right] \right) \right| \\
&= \sup_{x \in D} \left| \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \left( \mathbb{E} \left[ m_\rho(\mathbf{X}_i) K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) | \mathcal{F}_{i-1} \right] - \mathbb{E} \left[ m_\rho(\mathbf{X}_i) K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right] \right) \right| \\
&= \sup_{x \in D} \left| \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \int_{\mathbb{R}^d} m_\rho(\mathbf{y}) K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) (f^{\mathcal{F}_{i-1}}(\mathbf{y}) - f(\mathbf{y})) dy \right| \\
&= \sup_{x \in D} \left| \frac{(-1)^{|\beta|}}{h_n^{d+|\beta|}} \int_{\mathbb{R}^d} m_\rho(\mathbf{y}) K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) \left( \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}}(\mathbf{y}) - f(\mathbf{y}) \right) dy \right| \\
&\le C_m \sup_{x \in D} \left| \frac{(-1)^{|\beta|}}{h_n^{d+|\beta|}} \int_{\mathbb{R}^d} K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) \left( \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}}(\mathbf{y}) - f(\mathbf{y}) \right) dy \right|.
\end{align*}
$$

Making use of the Cauchy-Schwarz inequality and statement (7) when

$$m = d + 2|\beta| + 1,$$

we obtain readily that

$$
\begin{align*}
& \left| \frac{(-1)^{|\beta|}}{h_n^{d+|\beta|}} \int_{\mathbb{R}^d} K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) \left( \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}}(\mathbf{y}) - f(\mathbf{y}) \right) dy \right| \\
&\le \left( \int_{\mathbb{R}^d} \left| \frac{1}{h_n^{d+|\beta|}} K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) \right|^2 dy \right)^{1/2} \left( \int_{\mathbb{R}^d} \left| \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}}(\mathbf{y}) - f(\mathbf{y}) \right|^2 dy \right)^{1/2} \\
&\le \left\| \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}} - f \right\|_{L^2} \left( \int_{\mathbb{R}^d} \left( \frac{1}{h_n^{d+2|\beta|}} \frac{C_{d+2|\beta|+1}}{(1+h_n^{-1}\|\mathbf{x}-\mathbf{y}\|)^{d+2|\beta|+1}} \right) \left| \frac{1}{h_n^d} K^\beta \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) \right| dy \right)^{1/2} \\
&\le h_n^{1/2} C_{d+2|\beta|+1} \left\| \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}} - f \right\|_{L^2} \left( \int_{\mathbb{R}^d} K^\beta \left( \frac{\mathbf{x}}{h_n}, \mathbf{z} \right) |dz| \right)^{1/2} \\
&\le h_n^{1/2} C_{d+2|\beta|+1} \left\| \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}} - f \right\|_{L^2} G_1^{1/2}(d).
\end{align*}
$$

Under Assumption (C.2), we deduce that

$$
\begin{align}
G_{n,2}(\mathbf{x}) &= \sup_{x \in D} \left| \bar{r}_n(\rho; \mathbf{x}) - \mathbb{E}[(\widehat{\partial \beta r})_n(\rho; \mathbf{x})] \right| \nonumber \\
&\le C_m h_n^{1/2} C_{d+2|\beta|+1} G_1^{1/2}(d) \left\| \frac{1}{n} \sum_{i=1}^n f^{\mathcal{F}_{i-1}} - f \right\|_{L^2} \nonumber \\
&= O(h_n^{1/2}). \tag{57}
\end{align}
$$

Hence the proof is complete. □
**Proof of Theorem 3.** Recall the decomposition

$$
\begin{aligned}
& \sqrt{nh_n^{d+2|\beta|}} (\widehat{(\partial^\beta r)}_n(\rho; \mathbf{x}) - (\partial^\beta r)(\rho; \mathbf{x})) \\
& = \sqrt{nh_n^{d+2|\beta|}} \left( (\widehat{(\partial^\beta r)}_n(\rho; \mathbf{x}) - \bar{r}_n(\rho; \mathbf{x})) \right. \\
& \qquad \left. + (\bar{r}_n(\rho; \mathbf{x}) - (\widehat{\partial^\beta r})(\rho; \mathbf{x})) \right) \\
& = \sqrt{nh_n^{d+2|\beta|}} (Q_n(\rho; \mathbf{x}) + B_n(\rho; \mathbf{x})).
\end{aligned}
\tag{58}
$$

Observe that

$$
\begin{aligned}
B_n(\rho; \mathbf{x}) &= \left( \bar{r}_n(\rho; \mathbf{x}) - \mathbb{E}[\widehat{(\partial^\beta r)}_n(\rho; \mathbf{x})] \right) \\
&\quad + \left( \mathbb{E}[\widehat{(\partial^\beta r)}_n(\rho; \mathbf{x})] - \widehat{(\partial^\beta r)}(\rho; \mathbf{x}) \right) \\
&= B_{n,1}(\rho; \mathbf{x}) + B_{n,2}(\rho; \mathbf{x}).
\end{aligned}
\tag{59}
$$

Thus, employing reasoning analogous to that used for (57), under hypothesis (C.2), statement (7) with $m = 2(d + |\beta|)$, and Condition (12), we directly conclude that

$$
\begin{aligned}
& (nh_n^{d+2|\beta|})^{1/2} B_{n,1}(\mathbf{x}) \\
&= (nh_n^{d+2|\beta|})^{1/2} \left( \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \left( \mathbb{E}\left[\rho(Y_i) K^\beta\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) | \mathcal{F}_{i-1}\right] - \mathbb{E}\left[\rho(Y_i) K^\beta\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right] \right) \right) \\
&= O(h_n^{d/2}) (nh_n^{d+2|\beta|})^{1/2} \\
&= o(1).
\end{aligned}
\tag{60}
$$

We turn our attention to the term $B_{n,2}(\rho; \mathbf{x})$ and observe that under assumption (C.3), we have

$$
\begin{aligned}
\mathbb{E}[\widehat{(\partial^\beta r)}_n(\rho; \mathbf{x})] &= \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \mathbb{E}\left[\rho(Y_i) K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right] \\
&= \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \mathbb{E}\left[\rho(Y_i) |\mathbf{X}_i| K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right] \\
&= \frac{(-1)^{|\beta|}}{nh_n^{d+|\beta|}} \sum_{i=1}^n \int_{\mathbb{R}^d} m_\rho(\mathbf{u}) K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{u}}{h_n}\right) f(\mathbf{u}) d\mathbf{u} \\
&= \frac{1}{h_n^d} \int_{\mathbb{R}^d} (\partial^\beta r)(\rho; \mathbf{u}) K\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{u}}{h_n}\right) d\mathbf{u} \\
&= \int_{[-2L,2L]^d} (\partial^\beta r)(\rho; \mathbf{x} + h_n \mathbf{v}) K\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{x}}{h_n} + \mathbf{v}\right) d\mathbf{v} \\
&= (\partial^\beta r)(\rho; \mathbf{x}) + h_n (\partial^{\beta+1} r)(\rho; \mathbf{x}) \int_{[-2L,2L]^d} K(0, \mathbf{v}) d\mathbf{v}. \tag{61}
\end{aligned}
$$

On the other hand, by Condition (12), we have

$$
(nh_n^{d+2|\beta|})^{1/2} B_{n,2}(\rho; \mathbf{x}) = O(h_n^{\delta}(nh_n^{d+2|\beta|})^{1/2}) = o(1). \tag{62}
$$

Observe that

$$
\sqrt{nh_n^{d+2|\beta|}} Q_n(\mathbf{x}) = \sum_{i=1}^{n} \left\{ \left( \frac{(-1)^{|\beta|}}{\sqrt{nh_n^d}} \rho(Y_i) K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) \right) \right\}
$$
$$ -\mathbb{E}\left[\frac{(-1)^{|\beta|}}{\sqrt{n}h_n^d} \rho(Y_i) K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right) | \mathcal{F}_{i-1}\right] \quad (63) $$

$$ = \sum_{i=1}^{n} \{\tilde{\xi}_{ni}(\mathbf{x}_j, \mathbf{X}_i) - \mathbb{E}[\tilde{\xi}_{ni}(\mathbf{x}_j, \mathbf{X}_i) | \mathcal{F}_{i-1}]\} \\ = \sum_{i=1}^{n} \chi_{ni}(\mathbf{x}_j, \mathbf{X}_i), \quad (64) $$

Here, $\chi_{ni}(\mathbf{x}, \rho)$ constitutes a martingale difference array adapted to the filtration $\mathcal{F}_i$. This justifies applying the martingale central limit theorem for discrete-time arrays (cf. [75]) to demonstrate the asymptotic normality of $\sqrt{n}h_n^{d+2|\beta|}Q_n(\mathbf{x})$. To achieve this, it suffices to verify the following conditions:

(a) Lyapunov condition:

$$ \sum_{i=1}^{n} \mathbb{E} \left[ \chi_{ni}^2(\mathbf{x}_j, \mathbf{X}_i) \mid \mathcal{F}_{i-1} \right] \xrightarrow{\mathbb{P}} \Sigma^2(\mathbf{x}); $$

(b) Lindberg condition:

$$ n \mathbb{E} \left[ \chi_{ni}^2(\mathbf{x}_j, \mathbf{X}_i) 1_{\{|\chi_{ni}(\mathbf{x}_j, \mathbf{X}_i)| > \epsilon\}} \right] = o(1), \text{ holds true for any } \epsilon > 0. $$

**Proof of part (a).** Observe that

$$ \left| \sum_{i=1}^{n} \mathbb{E} \left[ \tilde{\xi}_{ni}^2(\mathbf{x}_j, \mathbf{X}_i) \mid \mathcal{F}_{i-1} \right] - \sum_{i=1}^{n} \mathbb{E} \left[ \chi_{ni}^2(\mathbf{x}_j, \mathbf{X}_i) \mid \mathcal{F}_{i-1} \right] \right| = \sum_{i=1}^{n} \left( \mathbb{E} \left[ \tilde{\xi}_{ni}(\mathbf{x}_j, \mathbf{X}_i) \mid \mathcal{F}_{i-1} \right] \right)^2. $$

Employing the same reasoning as in (54), and employing a first-order Taylor expansion alongside Equation (7), and in view of assumptions (N.2)(i), (N.3)(i), and (C.4)(ii), we obtain

$$ \begin{aligned} \mathbb{E}[\tilde{\xi}_{ni}(\mathbf{x}_j, \mathbf{X}_i) | \mathcal{F}_{i-1}] &= \frac{(-1)^{|\beta|}}{\sqrt{n}h_n^d} \mathbb{E} \left[ \rho(Y_i) K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) | \mathcal{F}_{i-1} \right] \\ &= \frac{(-1)^{|\beta|}}{\sqrt{n}h_n^d} \mathbb{E} \left[ \rho(Y_i) |G_{i-1}| K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) | \mathcal{F}_{i-1} \right] \\ &= \frac{(-1)^{|\beta|}}{\sqrt{n}h_n^d} \int_{\mathbb{R}^d} m_{\rho}(y) K^{\beta} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{y}}{h_n} \right) f^{\mathcal{F}_{i-1}}(\mathbf{y}) dy \\ &\le \frac{(-1)^{|\beta|} C_{\mathbf{m}} \sqrt{h_n^d}}{\sqrt{n}} \int_{\mathbb{R}^d} K^{\beta} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{x}}{h_n} + \mathbf{u} \right) f^{\mathcal{F}_{i-1}}(\mathbf{x} + h_n \mathbf{u}) d\mathbf{u} \\ &\le \left( \frac{(-1)^{|\beta|} C_{\mathbf{m}} \sqrt{h_n^d}}{\sqrt{n}} \right) \left( \frac{C_{d+|\beta|}}{(1 + \|\mathbf{u}\|_2)^{d+|\beta|}} \right) \left( f^{\mathcal{F}_{i-1}}(\mathbf{x}) + o(1) \right). \end{aligned} $$

Therefore, we infer

$$ \begin{aligned} & \sum_{i=1}^{n} \left( \mathbb{E} \left[ \tilde{\xi}_{ni}(\mathbf{x}_j, \mathbf{X}_i) \mid \mathcal{F}_{i-1} \right] \right)^2 \\ &= \left( h_n^d \right) \left( \frac{C_{d+|\beta|} C_{\mathbf{m}}}{(1 + \|\mathbf{u}\|_2)^{d+|\beta|}} \right)^2 \left( \frac{1}{n} \sum_{i=1}^{n} \left( f^{\mathcal{F}_{i-1}}(\mathbf{x}) \right)^2 + o(1) \right). \end{aligned} $$
The ergodicity of the process $(\mathbf{X}_n)_{n \ge 1}$ implies that the process defined by $((f^{\mathcal{F}_{i-1}}(\mathbf{x}))^2)_{n \ge 1}$ is ergodic too and verifies the condition (C.2), which means that

$$ \frac{1}{n} \sum_{i=1}^{n} (f^{\mathcal{F}_{i-1}}(\mathbf{x}))^2 = (f(\mathbf{x}))^2, $$

in both the almost sure and $L^2$ senses. This implies

$$ \sum_{i=1}^{n} (\mathbb{E}[\xi_{ni}(\mathbf{x}_j, \mathbf{X}_i) | \mathcal{F}_{i-1}])^2 = O(h_n^d). $$

The statement (a) follows then from

$$ \lim_{n \to \infty} \sum_{i=1}^{n} \mathbb{E} [\xi_{ni}^2(\mathbf{x}_j, \mathbf{X}_i) | \mathcal{F}_{i-1}] \stackrel{\mathbb{P}}{=} \Sigma_{\rho}^2(\mathbf{x}). $$

Observe that by the fact that (see [28])

$$ K^{(\beta)}(\mathbf{u}, \mathbf{v}) = K^{(\beta)}(\mathbf{u} + \mathbf{k}, \mathbf{v} + \mathbf{k}), \text{ for } \mathbf{k} \in \mathbb{Z}^d. $$

By assumption (N.2)(ii), we have

$$ \begin{aligned} & \sum_{i=1}^{n} \mathbb{E} [\xi_{ni}^2(\mathbf{x}_j, \mathbf{X}_i) | \mathcal{F}_{i-1}] \\ &= \sum_{i=1}^{n} \mathbb{E} \left[ \left( \frac{(-1)^{|\beta|}}{\sqrt{nh_n^d}} \rho(Y_i) K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right)^2 | \mathcal{F}_{i-1} \right] \\ &= \frac{1}{nh_n^d} \sum_{i=1}^{n} \mathbb{E} \left[ \mathbb{E} [\rho^2(Y_i) | \mathcal{G}_{i-1}] \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right)^2 | \mathcal{F}_{i-1} \right] \\ &= \frac{1}{nh_n^d} \sum_{i=1}^{n} \mathbb{E} \left[ m_{\rho,2}(\mathbf{X}_i) \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right)^2 | \mathcal{F}_{i-1} \right] \\ &\le \frac{1}{nh_n^d} \sum_{i=1}^{n} \mathbb{E} \left[ \sup_{\{\mathbf{u}: \|\mathbf{x}-\mathbf{u}\|<h\}} |m_{\rho,2}(\mathbf{X}_i) - m_{\rho,2}(\mathbf{x})| \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right)^2 | \mathcal{F}_{i-1} \right] \\ &\quad + \frac{1}{nh_n^d} \sum_{i=1}^{n} \mathbb{E} \left[ m_{\rho,2}(\mathbf{x}) \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right)^2 | \mathcal{F}_{i-1} \right] \\ &= I_{n,1}(\mathbf{x}) + I_{n,2}(\mathbf{x}). \end{aligned} \tag{65} $$

Under assumptions (C.1), (C.4)(ii), and (N.3)(iii) we have

$$ \begin{aligned} I_{n,1}(\mathbf{x}) &= \frac{1}{nh_n^d} \sum_{i=1}^{n} \mathbb{E} \left[ \sup_{\{\mathbf{u}: \|\mathbf{x}-\mathbf{u}\|<h\}} |m_{\rho,2}(\mathbf{X}_i) - m_{\rho,2}(\mathbf{x})| \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right)^2 | \mathcal{F}_{i-1} \right] \\ &= \frac{o(1)}{nh_n^d} \sum_{i=1}^{n} \mathbb{E} \left[ \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n} \right) \right)^2 | \mathcal{F}_{i-1} \right] \\ &= o(1) \frac{1}{n} \sum_{i=1}^{n} \int_{[-2L,2L]^d} \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{x}}{h_n} + \mathbf{v} \right) \right)^2 f^{\mathcal{F}_{i-1}}(\mathbf{x} + h_n \mathbf{v}) d\mathbf{v} \\ &= o(1) \left( \frac{1}{n} \sum_{i=1}^{n} f^{\mathcal{F}_{i-1}}(\mathbf{x}) + o(1) \right) \int_{[-2L,2L]^d} \left( K^{(\beta)} \left( \frac{\mathbf{x}}{h_n}, \frac{\mathbf{x}}{h_n} + \mathbf{v} \right) \right)^2 d\mathbf{v} \\ &= o(1)(f(\mathbf{x}) + o(1)) \int_{[-2L,2L]^d} \left( K^{(\beta)}(\mathbf{0}, \mathbf{v}) \right)^2 d\mathbf{v} \end{aligned} $$
$$ = o(1). \qquad (66) $$

Observe that

$$
\begin{align*}
I_{n,2}(\mathbf{x}) &= m_{\rho,2}(\mathbf{x}) \left(\frac{1}{nh_n^d}\right) \sum_{i=1}^{n} \mathbb{E}\left[\left(K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right)^2 | \mathcal{F}_{i-1}\right] \\
&= m_{\rho,2}(\mathbf{x}) \left(\frac{1}{nh_n^d}\right) \sum_{i=1}^{n} \int_{\mathbb{R}^d} \left(K^{(\beta)}\right)^2 \left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{u}}{h_n}\right) f_{\mathcal{F}_{i-1}}(\mathbf{u}) d\mathbf{u} \\
&= m_{\rho,2}(\mathbf{x}) \frac{1}{n} \sum_{i=1}^{n} \int_{\mathbb{R}^d} \left(K^{(\beta)}(\mathbf{0}, \mathbf{v})\right)^2 f_{\mathcal{F}_{i-1}}(\mathbf{x} + h_n \mathbf{v}) d\mathbf{v} \\
&= m_{\rho,2}(\mathbf{x}) \int_{\mathbb{R}^d} \left(K^{(\beta)}(\mathbf{0}, \mathbf{v})\right)^2 \left(\frac{1}{n} \sum_{i=1}^{n} f_{\mathcal{F}_{i-1}}(\mathbf{x}) + o(1)\right) d\mathbf{v}. \tag{67}
\end{align*}
$$

Combining Equations (66) and (67) we deduce

$$
\begin{align}
\lim_{n \to \infty} \sum_{i=1}^{n} \mathbb{E} \left[ \bar{\zeta}_{ni}^2 (\mathbf{x}_j, \mathbf{X}_i) \mid \mathcal{G}_{i-1} \right] &\stackrel{\mathbb{P}}{=} \Sigma_{\rho}^2(\mathbf{x}) \nonumber \\
&= m_{\rho,2}(\mathbf{x}) f(\mathbf{x}) \int_{\mathbb{R}^d} \left( K^{(\beta)}(\mathbf{0}, \mathbf{v}) \right)^2 d\mathbf{v}. \tag{68}
\end{align}
$$

**Proof of (b).** The Lindberg condition results from Corollary 9.5.2 in [76] which implies that

$$ n\mathbb{E}[\chi_{ni}^2(\mathbf{x}_j, \mathbf{X}_i)1_{\{|\chi_{ni}(\mathbf{x}_j, \mathbf{X}_i)| > \epsilon\}}] \le 4n\mathbb{E}[\bar{\zeta}_{ni}^2(\mathbf{x}_j, \mathbf{X}_i)1_{\{|\bar{\zeta}_{ni}(\mathbf{x}_j, \mathbf{X}_i)| > \epsilon/2\}}]. $$

Let $a > 1$ and $b > 1$ such that

$$ \frac{1}{a} + \frac{1}{b} = 1. $$

Making use of Hölder and Markov inequalities one can write, for all $\epsilon > 0$,

$$ \mathbb{E}[\bar{\zeta}_{ni}^2(\mathbf{x}_j, \mathbf{X}_i)1_{\{|\bar{\zeta}_{ni}(\mathbf{x}_j, \mathbf{X}_i)| > \epsilon/2\}}] \le \frac{\mathbb{E}|\bar{\zeta}_{ni}(\mathbf{x}_j, \mathbf{X}_i)|^{2a}}{(\epsilon/2)^{2a/b}}. $$

Therefore, by using the condition (7) when $m = d + 1$, we obtain

$$
\begin{align*}
& 4n\mathbb{E}[\bar{\zeta}_{ni}^2(\mathbf{x}_j, \mathbf{X}_i)1_{\{|\bar{\zeta}_{ni}(\mathbf{x}_j, \mathbf{X}_i)| > \epsilon/2\}}] \\
&\le \frac{4}{n^{a-1}h_n^{ad}(\epsilon/2)^{2a/b}} \mathbb{E}\left|\rho(Y_i)\left(K^{(\beta)}\left(\frac{\mathbf{x}}{h_n}, \frac{\mathbf{X}_i}{h_n}\right)\right)\right|^{2a} \\
&\le \frac{4}{n^{a-1}h_n^{ad}(\epsilon/2)^{2a/b}} \mathbb{E}\left|\rho(Y_i)\left(\frac{C_{d+1}}{(1+h_n^{-1}\|\mathbf{x}-\mathbf{X}_i\|_2^{d+1})}\right)\right|^{2a} \\
&= \frac{4C_{d+1}^{2a}}{n^{a-1}h_n^{(a-1)d}(\epsilon/2)^{2a/b}} \mathbb{E}|\rho(Y_i)|^{2a}.
\end{align*}
$$

Hence, under (N.0) we infer that

$$
\begin{align}
4n\mathbb{E}[\bar{\zeta}_{ni}^2(\mathbf{x}, \rho)1_{\{|\bar{\zeta}_{ni}(\mathbf{x}, \rho)| > \epsilon/2\}}] &\le \frac{4C_{d+1}^{2a}}{n^{a-1}h_n^{(a-1)d}(\epsilon/2)^{2a/b}} \nonumber \\
&= O\left(\left(\frac{1}{nh_n^d}\right)^{a-1}\right). \tag{69}
\end{align}
$$

Combining statements (68) and (69), we achieve the proof of the theorem. $\square$
**Author Contributions:** Conceptualization, S.D. and S.B.; methodology, S.D. and S.B.; validation, S.D. and S.B.; formal analysis, S.D. and S.B.; investigation, S.D. and S.B.; resources, S.D. and S.B.; writing—original draft preparation, S.D. and S.B.; writing—review and editing, S.D. and S.B. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research. The authors extend their sincere gratitude to the Editor-in-Chief, the Associate Editor, and the four reviewers for their invaluable feedback and for pointing out a number of oversights in the version initially submitted. Their insightful comments have greatly refined and focused the original work, resulting in markedly improved presentation.

**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. Besov Spaces

Following [74], let $1 \le p, q \le \infty$ and introduce the shift operator

$$ (S_{\tau}f)(x), f(x - \tau), \quad \tau \in \mathbb{R}^d. $$

Define

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

with the usual modification

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

For first-order regularity one replaces the first difference with a second difference:

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

and

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

For $0 < s < 1$ the Besov space

$$ B_{s,p,q} = \{f \in L_p(\mathbb{R}^d) : \ell_{s,p,q}(f) < \infty\}. $$

When $s > 1$ write

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

with $[s]^{-} \in \mathbb{N}$ and $0 < \{s\}^{+} \le 1$. Then $f \in B_{s,p,q}$ iff every weak derivative $D^j f$ with $|j| \le [s]^{-}$ satisfies

$$ D^j f \in B_{\{s\}^{+}, p,q}; $$

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

Prominent examples include the Sobolev space $H_2^s = B_{s,2,2}$ and the class of bounded s-Lipschitz functions $B_{s,\infty,\infty}$.

**Remark A1.** For density estimation over a Sobolev ball of smoothness s, the minimax $L^2$-risk is $O(n^{-2s/(2s+1)})$; see [77]. Comprehensive surveys of the relationships between classical function spaces and Besov spaces—including Fourier-analytic characterisations of Sobolev spaces when $p \ne 2$—are given in [78,79]. Connections between Besov spaces and the $V_p$ spaces of functions of bounded p-variation are developed in [80], relying on interpolation theory from [32]; a traditional exposition of p-variation may be found in [81]. For Besov spaces on more general geometric settings such as manifolds or Dirichlet spaces, consult [82].

## References

1. Wand, M.P.; Jones, M.C. Kernel smoothing. In *Monographs on Statistics and Applied Probability; Chapman and Hall, Ltd.: London*, UK, 1995; Volume 60, pp. xii+212. [CrossRef]

2. Eggermont, P.P.B.; LaRiccia, V.N. *Maximum Penalized Likelihood Estimation. Volume II. Regression*; Springer Series in Statistics; Springer: Dordrecht, The Netherlands, 2009; pp. xx+571. [CrossRef]

3. Genovese, C.R.; Perone-Pacifico, M.; Verdinelli, I.; Wasserman, L. Non-parametric inference for density modes. *J. R. Stat. Soc. Ser. B Stat. Methodol.* **2016**, *78*, 99–126. [CrossRef]

4. Noh, Y.K.; Sugiyama, M.; Liu, S.; du Plessis, M.C.; Park, F.C.; Lee, D.D. Bias reduction and metric learning for nearest-neighbor estimation of Kullback-Leibler divergence. *Neural Comput.* **2018**, *30*, 1930–1960. [CrossRef]

5. Dobrovidov, A.V.; Ruds'ko, I.M. Bandwidth selection in nonparametric estimator of density derivative by smoothed cross-validation method. *Autom. Remote Control* **2010**, *71*, 209–224. [CrossRef]

6. Genovese, C.R.; Perone-Pacifico, M.; Verdinelli, I.; Wasserman, L. On the path density of a gradient field. *Ann. Stat.* **2009**, *37*, 3236–3271. [CrossRef] [PubMed]

7. Singh, R.S. Applications of estimators of a density and its derivatives to certain statistical problems. *J. R. Stat. Soc. Ser. B* **1977**, *39*, 357–363. [CrossRef]

8. Meyer, T.G. Bounds for estimation of density functions and their derivatives. *Ann. Stat.* **1977**, *5*, 136–142. [CrossRef]

9. Silverman, B.W. Weak and strong uniform consistency of the kernel estimate of a density and its derivatives. *Ann. Stat.* **1978**, *6*, 177–184. [CrossRef]

10. 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]

11. Ramsay, J.O.; Silverman, B.W. *Applied Functional Data Analysis*; Springer Series in Statistics; Springer: New York, NY, USA, 2002; pp. x+190, Methods and case studies. [CrossRef]

12. Ramsay, J.O.; Silverman, B.W. *Functional Data Analysis*, 2nd ed.; Springer Series in Statistics; Springer: New York, NY, USA, 2005; pp. xx+426.

13. Liu, S.; Kong, X. A generalized correlated $C_p$ criterion for derivative estimation with dependent errors. *Comput. Statist. Data Anal.* **2022**, *171*, 107473. [CrossRef]

14. Eubank, R.L.; Speckman, P.L. Confidence bands in nonparametric regression. *J. Amer. Statist. Assoc.* **1993**, *88*, 1287–1301. [CrossRef]

15. Ruppert, D.; Sheather, S.J.; Wand, M.P. An effective bandwidth selector for local least squares regression. *J. Amer. Statist. Assoc.* **1995**, *90*, 1257–1270. [CrossRef]

16. Park, C.; Kang, K.H. SiZer analysis for the comparison of regression curves. *Comput. Statist. Data Anal.* **2008**, *52*, 3954–3970. [CrossRef]

17. Härdle, W.; Gasser, T. On robust kernel estimation of derivatives of regression functions. *Scand. J. Statist.* **1985**, *12*, 233–240.

18. Ziegler, K. On the asymptotic normality of kernel regression estimators of the mode in the nonparametric random design model. *J. Statist. Plann. Inference* **2003**, *115*, 123–144. [CrossRef]

19. Georgiev, A.A. Speed of convergence in nonparametric kernel estimation of a regression function and its derivatives. *Ann. Inst. Statist. Math.* **1984**, *36*, 455–462. [CrossRef]
20. Deheuvels, P.; Mason, D.M. General asymptotic confidence bands based on kernel-type function estimators. *Stat. Inference Stoch. Process.* **2004**, *7*, 225-277. [CrossRef]

21. Härdle, W.; Kerkyacharian, G.; Picard, D.; Tsybakov, A. Wavelets, approximation, and statistical applications. In *Lecture Notes in Statistics*, Springer: New York, NY, USA, 1998; Volume 129, pp. xviii+265. [CrossRef]

22. Prakasa Rao, B.L.S. Nonparametric estimation of the derivatives of a density by the method of wavelets. *Bull. Inf. Cybern.* **1996**, *28*, 91-100.

23. Chaubey, Y.P.; Doosti, H.; Prakasa Rao, B.L.S. Wavelet based estimation of the derivatives of a density with associated variables. *Int. J. Pure Appl. Math.* **2006**, *27*, 97-106.

24. Rao, B.L.S.P. Nonparametric Estimation of Partial Derivatives of a Multivariate Probability Density by the Method of Wavelets. In *Asymptotics in Statistics and Probability. Papers in Honor of George Gregory Roussas*; Puri, M.L., Ed.; De Gruyter: Berlin, Germany; Boston, MA, USA, 2000; pp. 321-330. [CrossRef]

25. Hosseinioun, N.; Doosti, H.; Niroumand, H.A. Nonparametric estimation of a multivariate probability density for mixing sequences by the method of wavelets. *Ital. J. Pure Appl. Math.* **2011**, *28*, 31-40.

26. Koshkin, G.; Vasil'iev, V. An estimation of a multivariate density and its derivatives by weakly dependent observations. In *Statistics and Control of Stochastic Processes. The Liptser Festschrift. Papers from the Steklov Seminar Held in Moscow, Russia, 1995-1996; World Scientific: Singapore, 1997*; pp. 229-241.

27. Prakasa Rao, B.L.S. Wavelet estimation for derivative of a density in the presence of additive noise. *Braz. J. Probab. Stat.* **2018**, *32*, 834-850. [CrossRef]

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

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

30. Schneider, C. Beyond Sobolev and Besov—Regularity of solutions of PDEs and their traces in function spaces. In *Lecture Notes in Mathematics*; Springer: Cham, Switzerland, 2021; Volume 2291, pp. xviii+327. [CrossRef]

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

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

33. Rosenblatt, M. Uniform ergodicity and strong mixing. *Z. Wahrscheinlichkeitstheorie Verw. Geb.* **1972**, *24*, 79-84. [CrossRef]

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

35. Didi, S.; Bouzebda, S. Wavelet Density and Regression Estimators for Continuous Time Functional Stationary and Ergodic Processes. *Mathematics* **2022**, *10*, 4356. [CrossRef]

36. 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, pp. xx+357. [CrossRef]

37. Peškir, G. The uniform mean-square ergodic theorem for wide sense stationary processes. *Stoch. Anal. Appl.* **1998**, *16*, 697-720. [CrossRef]

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

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

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

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

42. Bouzebda, S. Limit Theorems in the Nonparametric Conditional Single-Index U-Processes for Locally Stationary Functional Random Fields under Stochastic Sampling Design. *Mathematics* **2024**, *12*, 1996. [CrossRef]

43. Shen, A.; Li, X.; Zhang, Y.; Qiu, Y.; Wang, X. Consistency and asymptotic normality of wavelet estimator in a nonparametric regression model. *Stochastics* **2021**, *93*, 868-885. [CrossRef]

44. Liu, Y.; Zeng, X. Asymptotic normality for wavelet deconvolution density estimators. *Appl. Comput. Harmon. Anal.* **2020**, *48*, 321-342. [CrossRef]

45. Tang, X.; Xi, M.; Wu, Y.; Wang, X. Asymptotic normality of a wavelet estimator for asymptotically negatively associated errors. *Statist. Probab. Lett.* **2018**, *140*, 191-201. [CrossRef]

46. Niu, S.I. Asymptotic normality of wavelet density estimator under censored dependent observations. *Acta Math. Appl. Sin. Engl. Ser.* **2012**, *28*, 781-794. [CrossRef]

47. Li, Y.; Guo, j. Asymptotic normality of wavelet estimator for strong mixing errors. *J. Korean Statist. Soc.* **2009**, *38*, 383-390. [CrossRef]
48. Roueff, F.; Taqqu, M.S. Asymptotic normality of wavelet estimators of the memory parameter for linear processes. *J. Time Ser. Anal.* **2009**, *30*, 534–558. [CrossRef]

49. Li, Y.; Yang, S.; Zhou, Y. Consistency and uniformly asymptotic normality of wavelet estimator in regression model with associated samples. *Statist. Probab. Lett.* **2008**, *78*, 2947–2956. [CrossRef]

50. Debbah, 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]

51. Allaoui, S.; Bouzebda, S.; Liu, J. Asymptotic distribution of the wavelet-based estimators of multivariate regression functions under weak dependence. *J. Math. Inequal.* **2023**, *17*, 481–515. [CrossRef]

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

53. Bouzebda, S.; Didi, S. Some asymptotic properties of kernel regression estimators of the mode for stationary and ergodic continuous time processes. *Rev. Mat. Complut.* **2021**, *34*, 811–852. [CrossRef]

54. Eddy, W.F. Optimum kernel estimators of the mode. *Ann. Statist.* **1980**, *8*, 870–882. [CrossRef]

55. Eddy, W.F. The asymptotic distributions of kernel estimators of the mode. *Z. Wahrsch. Verw. Geb.* **1982**, *59*, 279–290. [CrossRef]

56. Billingsley, P. *Convergence of Probability Measures*; John Wiley & Sons, Inc.: New York, NY, USA; London, UK; Sydney, Australia, 1968; pp. xii+253.

57. Bouzebda, S.; Chaouch, M. Uniform limit theorems for a class of conditional Z-estimators when covariates are functions. *J. Multivar. Anal.* **2022**, *189*, 104872. [CrossRef]

58. Hall, P.; Marron, J.S. Estimation of integrated squared density derivatives. *Statist. Probab. Lett.* **1987**, *6*, 109–115. [CrossRef]

59. Jurečková, J. Asymptotic linearity of a rank statistic in regression parameter. *Ann. Math. Statist.* **1969**, *40*, 1889–1900. [CrossRef]

60. Giné, E.; Mason, D.M. Uniform in bandwidth estimation of integral functionals of the density function. *Scand. J. Statist.* **2008**, *35*, 739–761. [CrossRef]

61. Levit, B.Y. Asymptotically efficient estimation of nonlinear functionals. *Probl. Peredachi Informatsii* **1978**, *14*, 65–72.

62. 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]

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

64. Bouzebda, S. Limit theorems for wavelet conditional U-statistics for time series models. *Math. Methods Statist.* **2025**, *35*, 1–42.

65. Chaouch, M.; Laïb, N. Regression estimation for continuous-time functional data processes with missing at random response. *J. Nonparametr. Stat.* **2024**, *36*, 1–32. [CrossRef]

66. Giraitis, L.; Leipus, R. A generalized fractionally differencing approach in long-memory modeling. *Liet. Mat. Rink.* **1995**, *35*, 65–81. [CrossRef]

67. Guégan, D.; Ladoucette, S. Non-mixing properties of long memory processes. *C. R. Acad. Sci. Paris Sér. I Math.* **2001**, *333*, 373–376. [CrossRef]

68. Andrews, D.W.K. Non-strong mixing autoregressive processes. *J. Appl. Probab.* **1984**, *21*, 930–934. [CrossRef]

69. Francq, C.; Zakoïan, J.M. GARCH Models: Structure, Statistical Inference and Financial Applications; John Wiley & Sons, Ltd.: Chichester, UK, 2010; pp. xiv+489. [CrossRef]

70. Bouzebda, S.; Limnios, N. On general bootstrap of empirical estimator of a semi-Markov kernel with applications. *J. Multivar. Anal.* **2013**, *116*, 52–62. [CrossRef]

71. Bouzebda, S. On the weak convergence and the uniform-in-bandwidth consistency of the general conditional *U*-processes based on the copula representation: Multivariate setting. *Hacet. J. Math. Stat.* **2023**, *52*, 1303–1348. [CrossRef]

72. Burkholder, D.L. Distribution function inequalities for martingales. *Ann. Probab.* **1973**, *1*, 19–42. [CrossRef]

73. de la Peña, V.H.; Giné, E. *Decoupling*; Probability and Its Applications (New York); Springer: New York, NY, USA, 1999; pp. xvi+392. [CrossRef]

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

75. Hall, P.; Heyde, C.C. *Martingale Limit Theory and Its Application*; Probability and Mathematical Statistics; Academic Press, Inc.; Harcourt Brace Jovanovich, Publishers: New York, NY, USA; London, UK, 1980; pp. xii+308.

76. Chow, Y.S.; Teicher, H. *Probability Theory*; Springer: New York, NY, USA; Berlin/Heidelberg, Germany, 1978; pp. xv+455.

77. Efromovich, S. Lower bound for estimation of Sobolev densities of order less $\frac{1}{2}$. *J. Statist. Plann. Inference* **2009**, *139*, 2261–2268. [CrossRef]

78. Triebel, H. Theory of function spaces. In *Monographs in Mathematics*; Birkhäuser Verlag: Basel, Switzerland, 1983; Volume 78, p. 284. [CrossRef]

79. DeVore, R.A.; Lorentz, G.G. Constructive approximation. In *Grundlehren der Mathematischen Wissenschaften* [Fundamental Principles of Mathematical Sciences]; Springer: Berlin/Heidelberg, Germany, 1993; Volume 303, pp. x+449.
80. Bourdaud, G.; Lanza de Cristoforis, M.; Sickel, W. Superposition operators and functions of bounded *p*-variation. *Rev. Mat. Iberoam.* **2006**, *22*, 455–487. [CrossRef]

81. Dudley, R.M.; Norvaiša, R. Differentiability of six operators on nonsmooth functions and *p*-variation; In *Lecture Notes in Mathematics*, Springer: Berlin/Heidelberg, Germany, 1999; Volume 1703, pp. viii+277. [CrossRef]

82. Geller, D.; Pesenson, I.Z. Band-limited localized Parseval frames and Besov spaces on compact homogeneous manifolds. *J. Geom. Anal.* **2011**, *21*, 334–371. [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.