Non-Recursive estimator $\tilde{\pi}_n$:

To estimate the optimal bandwidth (2.24), we need to estimate $I_1$ and $I_2$. Therefore, we can state

$$ \tilde{a}_n(x, y) = \frac{1}{nh_n^d} \sum_{k=1}^{n} \chi_k(y) \mathbf{K} \left( \frac{x - X_k}{h_n} \right) = \frac{1}{nh_n^d} \sum_{k=1}^{n} \prod_{i=1}^{d} K \left( \frac{x - X_{ki}}{h_k} \right) \chi_{ki}(y) $$

and

$$ \tilde{f}_n(x) = \frac{1}{nh_n^d} \sum_{k=1}^{n} \mathbf{K} \left( \frac{x - X_k}{h_n} \right) = \frac{1}{nh_n^d} \sum_{k=1}^{n} \prod_{i=1}^{d} K \left( \frac{x_i - X_{ki}}{h_k} \right). $$

Estimation of $I_1$:

$$ \tilde{J}_1 = \frac{1}{n^3 b_n^6} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k}}^n \left[ \sum_{v=1}^d K_{b'}^{(2)} \left( \frac{X_{iv} - X_{jv}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq v}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{v=1}^d K_{b'}^{(2)} \left( \frac{X_{iv} - X_{kv}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq v}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] \prod_{s=1}^q \chi_{js}(Y_{is}) \chi_{ks}(Y_{is}), $$

$$ \tilde{J}_2 = \frac{1}{n^3 b_n^6} \sum_{\substack{i,j,k,u=1 \\ i \neq j \neq k \neq u \neq m}}^n \left[ \sum_{v=1}^d K_{b'}^{(2)} \left( \frac{X_{iv} - X_{jv}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq v}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{v=1}^d K_{b'}^{(2)} \left( \frac{X_{iv} - X_{kv}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq v}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] \prod_{s=1}^q \chi_{us}(Y_{is}) \chi_{js}(Y_{is}), $$

$$ \tilde{J}_3 = \frac{1}{n^4 b_n^6} \sum_{\substack{i,j,k,u,m=1 \\ i \neq j \neq k \neq u \neq m}}^n \left[ \sum_{v=1}^d K_{b'}^{(2)} \left( \frac{X_{iv} - X_{jv}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq v}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{v=1}^d K_{b'}^{(2)} \left( \frac{X_{iv} - X_{kv}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq v}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] \prod_{s=1}^q \chi_{us}(Y_{is}) \chi_{ms}(Y_{is}), $$

Then, we obtain

$$ \tilde{I}_1 = \mu^2(K) (\tilde{J}_1 - 2\tilde{J}_2 + \tilde{J}_3). $$

Estimation of $I_2$:

$$ \tilde{I}_2 = \frac{1}{n^2 b_n} \sum_{\substack{i,k,u=1 \\ i \neq k \neq u}}^n \prod_{l=1}^d K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \prod_{s=1}^q \chi_{us}(Y_{is}) (1 - \chi_{ks}(Y_{is})), $$

As a result, the plug-in estimator of (2.24) is denoted by

$$ (h_n) = \left( \left( \frac{\tilde{I}_2}{\tilde{I}_1} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{1}{d+4}} \right), \qquad (2.26) $$
Finally, a non-recursive estimator of MWISE[$\pi_n$] is provided by

$$ \widetilde{MWISE}[\tilde{\pi}_n] = \frac{5}{4} (\tilde{I}_2)^{\frac{4}{d+4}} (\tilde{I}_1)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{4}{d+4}} + \mathbf{o}(n^{-\frac{4}{d+4}}). $$

The major aim of our next section lies in comparing the performance of our recursive estimator (2.4) with that of non-recursive Nadaraya-Watson one (2.6).

## 2.6 Numerical applications

Let's start our numerical studies with some simulations with different dimensions Models.

### 2.6.1 Simulation studies

In order to compare the proposed recursive estimator with the Nadaraya-Watson non-recursive one, we consider three sample sizes: n=100, 200 and 500, a fixed number of simulations: N=500 and four distribution models:

* Model 1: $(X, Y) \in \mathbb{R} \times \mathbb{R}$:
  $Y = 2\sin(\pi X) + \epsilon$, where X follows the binomial distribution $\mathcal{B}(2, 1/3)$ and $\epsilon$ follows the normal distribution $\mathcal{N}(0, 1)$.

* Model 2: $(X, Y) \in \mathbb{R}^2 \times \mathbb{R}$:
  $Y = \exp(-X/2) + \epsilon$, where X follows the poisson distribution $\mathbb{P}(1/2)$ and $\epsilon$ follows the normal distribution $\mathcal{N}(0, 1/2)$.

* Model 3: $(X, Y) \in \mathbb{R}^3 \times \mathbb{R}^2$:
  $Y = AX + \epsilon$ with $A = \begin{pmatrix} 1 & 1 & 1 \\ 1 & -2 & 1 \end{pmatrix}$, $X = 0 \times 1_{Z \le 0.5} + 1_{Z > 0.5}$ where Z follows the uniform distribution $\mathcal{U}(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \end{pmatrix})$ and $\epsilon$ follows the normal distribution $\mathcal{N}(0, 1/2)$.

* Model 4: $(X, Y) \in \mathbb{R}^{10} \times \mathbb{R}^{10}$:
  $Y = \exp(X) + \epsilon$, with $X = \lfloor Z \rfloor$ where Z follows the 10-dimensional normal distribution $\mathcal{N}(0_{10}, I_{10})$ and $\epsilon$ follows the normal distribution $\mathcal{N}(0, 1)$.

We denote by $\pi_i^*$ the reference CCDF and by $\pi_i$ the test CCDF. Then, we calculate the following two measures:

* Mean squared error: $MSE = \frac{1}{n} \sum_{i=1}^{n} (\pi_i - \pi_i^*)^2$.

* The linear correlation: $Cor = \frac{Cov(\pi_i, \pi_i^*)}{\sigma(\pi_i)\sigma(\pi_i^*)}$.

In what follows, we portray the different steps of the simulation algorithm in the multivariate case.
2.6.2 Simulation Algorithm

**Algorithm 1** K is the Gaussian kernel, *d* the dimension size, *n* the simple size, *Np* the number of observations and *N* the number of iterations.

**Input:** K, d, n, Np and N.

1: A random initialization of $\hat{\Pi}^{(0)}$ (resp. $\tilde{\Pi}^{(0)}$).

2: **for** *l* = 1, ..., *N* **do**

3: A random sample vectors $X_1, \dots, X_d$ and $Y$ of length $n$.

4: A choice value for the recursive bandwidth vectors $h_1, \dots, h_n$. (resp. the non-recursive bandwidth values $h_n$) using the plug-in approach given in (2.25) (resp. (2.26)).

5: The choice of the stepsize ($\gamma_n$) = ($n^{-1}$).

6: We fix $x_1, \dots, x_d$ and consider an arbitrary sampling vector $T$ of $y$ of length $N_p$.

7: $\hat{\pi}_l(y|x) = \frac{\sum_{k=1}^{n} k\gamma_k 1_{\{Y_k \le y\}} h_k^{-d} \prod_{i=1}^{d} K\left(\frac{x_i - X_{k_i}}{h_k}\right)}{\sum_{k=1}^{n} k\gamma_k h_k^{-d} \prod_{i=1}^{d} K\left(\frac{x_i - X_{k_i}}{h_k}\right)}$ for the multivariate recursive CCDF estimator. (resp. $\tilde{\pi}_l(y|x) = \frac{\sum_{k=1}^{n} 1_{\{Y_k \le y\}} \prod_{i=1}^{d} K\left(\frac{x_i - X_{k_i}}{h_k}\right)}{\sum_{k=1}^{n} \prod_{i=1}^{d} K\left(\frac{x_i - X_{k_i}}{h_k}\right)}$ for the multivariate non-recursive CCDF estimator).

$$ \hat{\Pi}^{(l)} = \hat{\pi}_l(T|x). \quad (\text{resp. } \tilde{\Pi}^{(l)} = \tilde{\pi}_l(T|x).) $$

8: **end for**

9: $\hat{\pi} = N^{-1} \sum_{l=1}^{N} \hat{\Pi}^{(l)}$ (resp. $\tilde{\pi} = N^{-1} \sum_{l=1}^{N} \tilde{\Pi}^{(l)}$).

**output:** The vectors $\hat{\pi}$ and $\tilde{\pi}$.
<table><caption>Table 2.1: Quantitative comparison between Nadaraya-Watson estimator and the proposed estimator with stepsizes ($\gamma_n$) = ($n^{-1}$) through a plug-in method for Model 1.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">MSE/Cor</th><th rowspan="2">n</th><th colspan="2">x = 0</th><th colspan="2">x = 2</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="6">Model 1</td><td rowspan="3">MSE</td><td>100</td><td>3.671766e-07</td><td><b>1.333452e-07</b></td><td>1.342202e-05</td><td><b>8.228683e-06</b></td></tr><tr><td>200</td><td>3.377698e-07</td><td><b>1.185438e-07</b></td><td>3.434738e-06</td><td><b>2.050298e-06</b></td></tr><tr><td>500</td><td>1.748553e-07</td><td><b>4.915888e-08</b></td><td>1.407363e-06</td><td><b>6.641233e-07</b></td></tr><tr><td rowspan="3">Cor</td><td>100</td><td>9.999992e-01</td><td><b>9.999997e-01</b></td><td>9.999773e-01</td><td><b>9.999849e-01</b></td></tr><tr><td>200</td><td>9.999992e-01</td><td><b>9.999997e-01</b></td><td>9.999914e-01</td><td><b>9.999947e-01</b></td></tr><tr><td>500</td><td>9.999997e-01</td><td><b>9.999999e-01</b></td><td>9.999962e-01</td><td><b>9.999981e-01</b></td></tr></tbody></table>

<table><caption>Table 2.2: Quantitative comparison between the recursive estimator and the non-recursive one with stepsizes ($\gamma_n$) = ($n^{-1}$) through a plug-in method for Model 2.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">MSE/Cor</th><th rowspan="2">n</th><th colspan="2">x = (0, 0)</th><th colspan="2">x = (1, 1)</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="6">Model 2</td><td rowspan="3">MSE</td><td>100</td><td>0.005299303</td><td><b>0.002345599</b></td><td>0.011464770</td><td><b>0.005503603</b></td></tr><tr><td>200</td><td>0.004444445</td><td><b>0.001976603</b></td><td>0.009766589</td><td><b>0.004850061</b></td></tr><tr><td>500</td><td>0.003609582</td><td><b>0.001795659</b></td><td>0.006988903</td><td><b>0.004007790</b></td></tr><tr><td rowspan="3">Cor</td><td>100</td><td>0.993278442</td><td><b>0.997042336</b></td><td>0.985998000</td><td><b>0.993297472</b></td></tr><tr><td>200</td><td>0.993932848</td><td><b>0.997313718</b></td><td>0.987516900</td><td><b>0.993740343</b></td></tr><tr><td>500</td><td>0.994773451</td><td><b>0.997402840</b></td><td>0.990016792</td><td><b>0.994209320</b></td></tr></tbody></table>

<table><caption>Table 2.3: Quantitative comparison between the recursive estimator and the non-recursive one with stepsizes ($\gamma_n$) = ($n^{-1}$) through a plug-in method for Model 3.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">MSE/Cor</th><th rowspan="2">n</th><th colspan="2">x = (0, 0, 0)</th><th colspan="2">x = (1, 1, 1)</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="6">Model 3</td><td rowspan="3">MSE</td><td>100</td><td>0.005170397</td><td><b>0.003118756</b></td><td>0.007481887</td><td><b>0.002756530</b></td></tr><tr><td>200</td><td>0.004968340</td><td><b>0.002747727</b></td><td>0.007215537</td><td><b>0.002449973</b></td></tr><tr><td>500</td><td>0.004639918</td><td><b>0.002444590</b></td><td>0.007200206</td><td><b>0.002260629</b></td></tr><tr><td rowspan="3">Cor</td><td>100</td><td>0.989849852</td><td><b>0.993964835</b></td><td>0.989210840</td><td><b>0.996122260</b></td></tr><tr><td>200</td><td>0.990247670</td><td><b>0.994651835</b></td><td>0.989322443</td><td><b>0.996398260</b></td></tr><tr><td>500</td><td>0.990928518</td><td><b>0.995271000</b></td><td>0.989458730</td><td><b>0.996602105</b></td></tr></tbody></table>

<table><caption>Table 2.4: Quantitative comparison between the recursive estimator and the non-recursive one with stepsizes ($\gamma_n$) = ($n^{-1}$) through a plug-in method for Model 4.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">MSE/Cor</th><th rowspan="2">n</th><th colspan="2">x = (0, 0, 0, 0, 0, 0, 0, 0)</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="6">Model 4</td><td rowspan="3">MSE</td><td>100</td><td>0.012132600</td><td><b>0.008651937</b></td></tr><tr><td>200</td><td>0.012042460</td><td><b>0.007971624</b></td></tr><tr><td>500</td><td>0.009422018</td><td><b>0.006973041</b></td></tr><tr><td rowspan="3">Cor</td><td>100</td><td>0.985115900</td><td><b>0.989107240</b></td></tr><tr><td>200</td><td>0.982886350</td><td><b>0.987784007</b></td></tr><tr><td>500</td><td>0.984891065</td><td><b>0.987888565</b></td></tr></tbody></table>
Departing from Tables 4.1, 2.2 and 2.4, we conclude that:

1. The MSE of the proposed recursive estimator with stepsize $(\gamma_n) = (n^{-1})$ through a plug-in method is smaller than that of Nadaraya-Watson's non-recursive estimator.

2. The estimators get closer to the true CCDF function as sample size increases, i.e., the *MSE* decreases as the simple size increases and therefore the *Cor* increases as the sample size increases.

<figure><img src="image_22.png" /><figcaption>Figure 2.1: The reference CCDF for Model 1 for one simple simulation.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 2.2: The recursive CCDF estimator for Model 1 for one simple simulation.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 2.3: The non-recursive CCDF estimator for Model 1 for one simple simulation.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 2.4: Qualitative comparison between the recursive estimator and the non-recursive one for Model 1 with $n = 200$, $N = 500$ and $x = 0$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 2.5: Qualitative comparison between the recursive estimator and the non-recursive one for Model 1 with $n = 500$, $N = 500$ and $x = 0$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 2.6: Qualitative comparison between the recursive estimator and the non-recursive one for Model 2 with $n = 100$, $N = 500$ and $x = (0,0)$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 2.7: Qualitative comparison between the recursive estimator and the non-recursive one for Model 2 with $n = 500$, $N = 500$ and $x = (0,0)$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 2.8: The reference CCDF for Model 3 for one simple simulation with $n = 500$ and $x = (1,1,1)$.</figcaption></figure>

<caption>Figure 2.9: The recursive CCDF estimator for Model 3 for one simple simulation with $n = 500$ and $x = (1,1,1)$.</caption>

<caption>Figure 2.10: The non-recursive CCDF estimator for Model 3 for one simple simulation with $n = 500$ and $x = (1,1,1)$.</caption>
### 2.6.3 Real Datasets:

In this section, our focal point is to examine two real datasets Models, namely the Insurance Company Benchmark (COIL 2000) dataset as well as the French Hospital Data of COVID-19.

#### Application 1: The Insurance Company Benchmark (COIL 2000) dataset

The (COIL 2000) dataset is found in data.world website

https://data.world/uci/insurance-company-benchmark-coil-2000).

Information about customers consists of 86 variables and includes product usage data and sociodemographic data derived from zip area codes. The data are supplied by the Dutch data mining company Sentient Machine Research and rest on a real world business problem. The training set involves over 5000 descriptions of customers, including the information of whether or not they have a caravan insurance policy. A test set includes 4000 customers whom only the organizers know if they have a caravan insurance policy. This corresponds to a Dataset to train and validate prediction models and build up a description (5822 customer records). Each record consists of 86 attributes, incorporating sociodemographic data (attribute 1-43) and product ownership (attributes 44-86).

The sociodemographic data are derived from zip codes. All customers living in areas with the same zip code have the same sociodemographic attributes.

As far as our application is concerned, we shall consider the following two models:

* Model 1: X corresponds to the sociodemographic data attribute number 16 and Y stands for the whole 5822 observations of customer records.

* Model 2: X corresponds to the sociodemographic 5-dimensional data attributes number 6,8,11,12 and 13 and Y stands for the whole 5822 observations of customer records.

<table><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">x=0</th><th colspan="2">x=1</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="2">Model 1</td><td>MSE</td><td>0.006979541</td><td><strong>0.004748803</strong></td><td>0.001999743</td><td><strong>0.001390595</strong></td></tr><tr><td>Cor</td><td>0.987112335</td><td><strong>0.990653245</strong></td><td>0.995770041</td><td><strong>0.996574025</strong></td></tr></tbody></table>

<table><caption>Table 2.5: Quantitative comparison between Nadaraya-Watson estimator and the proposed estimator with stepsizes ($\gamma_n$) = ($n^{-1}$) through a plug-in method for the Insurance Company Benchmark (COIL 2000) dataset case.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">x=(0,0,0,0,0)</th><th colspan="2">x=(1,1,1,2,2)</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="2">Model 2</td><td>MSE</td><td>0.000708417</td><td><strong>0.0005969373</strong></td><td>0.00414241</td><td><strong>0.003705582</strong></td></tr><tr><td>Cor</td><td>0.996664739</td><td><strong>0.9979940979</strong></td><td>0.97540614</td><td><strong>0.976614363</strong></td></tr></tbody></table>
<figure><img src="image_30.png" /><figcaption>Figure 2.11: Qualitative comparison between the recursive estimator and the non-recursive one for the dataset Model 1 with $x = 1$.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 2.12: Qualitative comparison between the recursive estimator and the non-recursive one for the dataset Model 2 with $x = (0, 0, 0, 0, 0)$.</figcaption></figure>

## Application 2: French Hospital Data of COVID19

The French Hospital data of the COVID-19 epidemic are extracted from https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/.

The **Santé publique France**'s mission is to improve and protect the health of populations. During the health crisis linked to the COVID-19 epidemic, **Santé publique France** has taken in charge monitoring and understanding the dynamics of the epidemic, anticipating the different scenarios and implementing actions to prevent and limit the transmission of this virus on the national territory.

### Description of the dataset

This dataset provides information on the hospital situation regarding the COVID-19 epidemic. We have opted for the first proposed file:

Hospital data related to the COVID-19 epidemic by department (dep), sex of the patient (sex), number of hospitalized patients (hosp), number of persons currently in intensive care or resuscitation (rea), number of persons currently in follow-up and rehabilitation care (SSR) or long-term care units (USLD), number of persons currently in conventional hospitalization (HospConv), number of persons currently hospitalized in another type of service (autres), cumulative number of persons returning home (rad) or cumulative number of dead persons (dc).

The data have been daily updated. For the current application, we considered the data of 28/07/2021, with a total of 150894 observations. For simplicity reason, we have chosen to just study the department of 'Vienne' database.

Therefore, for our application, we served of a dataframe of 1494 observations and 6 variables. Hence, we shall consider the following three models:

* Model 1: $X = dc$ and $Y = hosp$.

* Model 2: $X1 = sex, X2 = rea, X3 = dc$ and $Y = hosp$.

* Model 3: $X1 = sex, X2 = rea, X3 = dc, Y1 = hosp$ and $Y2 = rad$.

<table><thead><tr><th rowspan="2"></th><th colspan="2">x=119</th><th colspan="2">x=17</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td>Model 1</td><td>MSE</td><td>0.02538461</td><td><strong>0.01900265</strong></td><td>0.008079359</td><td><strong>0.007509241</strong></td></tr><tr><td></td><td>Cor</td><td>0.86514927</td><td><strong>0.88894736</strong></td><td>0.781752320</td><td><strong>0.789593904</strong></td></tr></tbody></table>
<table><caption>Table 2.6: Quantitative comparison between Nadaraya-Watson estimator and the proposed estimator with stepsizes ($\gamma_n$) = ($n^{-1}$) through a plug-in method for the COVID-19 epidemic dataset case.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">x = (2, 0, 17)</th><th colspan="2">x = (1, 0, 0)</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="2">Model 2</td><td>MSE</td><td>0.02176777</td><td><strong>0.01416798</strong></td><td>0.01026319</td><td><strong>0.005681196</strong></td></tr><tr><td>Cor</td><td>0.55644221</td><td><strong>0.63197054</strong></td><td>0.74621610</td><td><strong>0.819493783</strong></td></tr><tr><td rowspan="2">Model 3</td><td>MSE</td><td>0.02107916</td><td><strong>0.009991914</strong></td><td>0.04130757</td><td><strong>0.03504694</strong></td></tr><tr><td>Cor</td><td>0.85362436</td><td><strong>0.920385555</strong></td><td>0.61896604</td><td><strong>0.64659332</strong></td></tr></tbody></table>

<figure><img src="image_32.png" /><figcaption>Figure 2.13: Qualitative comparison between the recursive estimator and the non-recursive one for the COVID-19 epidemic dataset Model 1 with $x = 17$.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 2.14: Qualitative comparison between the recursive estimator and the non-recursive one for the COVID-19 dataset Model 2 with $x = (2, 0, 17)$.</figcaption></figure>

### Data interpretation:

Referring to Tables 2.5, 2.6 and Figures 2.11, 2.12, 2.13 and 2.14, we conclude that:

1. For all considered Models, the proposed recursive estimator with stepsize ($\gamma_n$) = ($n^{-1}$) through a plug-in method outperformed the non-recursive one in terms of estimation error MSE and Cor.

2. The proposed recursive estimator is closer to the true CCDF function, compared with Nadaraya-Watson's non-recursive estimator.

Concerning the COVID-19 epidemic Model 1, we can infer that, for a fixed number of deaths $x = 17$, the proportion of hospitalized cases less than 20 is 50% and the proportion of hospitalized cases less than 50 is 85%. Moreover, 99% of the population have a number of hospitalized cases less than 100.

Likewise, the COVID-19 epidemic Model 2 yielded the same results as Model 1. Indeed, we recorded for the women gender a fixed value of sex $x_1 = 2$, a fixed number of REA persons $x_2 = 0$ and deaths $x_3 = 17$.

## 2.7 Conclusion

In this work, we elaborated a multivariate recursive CCDF estimator. We tackled the asymptotic properties of the proposed estimator by providing the bias as well as the variance in order to demonstrate that our estimator asymptotically follows a normal distribution. Subsequently, we revealed that the use of our recursive estimator with an appropriate choice of the bandwidth and the stepsize enables us to get closer to the true conditional cumulative distribution function rather than non-recursive one.
## 2.8 Proofs

Throughout this section, devoted to the proofs of our main results, it is noteworthy that

For all $x \in \mathbb{R}^d, y \in \mathbb{R}^q$,

$$ \mathcal{Z}_n(x, y) := h_n^{-d} \chi_n(y) \mathbf{K} \left( \frac{x - X_n}{h_n} \right) \quad \text{and} \quad W_n(x) := h_n^{-d} \mathbf{K} \left( \frac{x - X_n}{h_n} \right). $$

**Proof of Proposition 2.2.**

Our proof starts with the observation that, based on the expression (2.1) without assuming $a_0(x, y) = 0$ and with recurrence on $n$, we get

$$ \begin{aligned} a_n(x, y) - a(x, y) &= (1 - \gamma_n)a_{n-1}(x, y) + \gamma_n \mathcal{Z}_n(x, y) - a(x, y) \\ &= (1 - \gamma_n)[a_{n-1}(x, y) - a(x, y)] + \gamma_n[\mathcal{Z}_n(x, y) - a(x, y)] \\ &= \prod_{i=1}^{n} (1 - \gamma_i)[a_0(x, y) - a(x, y)] + \sum_{k=1}^{n-1} \prod_{i=k+1}^{n} (1 - \gamma_i)\gamma_k(\mathcal{Z}_k(x, y) - a(x, y)) \\ &\quad + \gamma_n(\mathcal{Z}_n(x, y) - a(x, y)). \end{aligned} $$

Grounded on the fact that, $\Pi_n = \prod_{i=1}^{n} (1 - \gamma_i)$ and $\prod_{i=k+1}^{n} (1 - \gamma_i) = \Pi_n \Pi_k^{-1}$, we deduce

$$ \forall x \in \mathbb{R}^d, y \in \mathbb{R}^q, \quad a_n(x, y) - a(x, y) = \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k (\mathcal{Z}_k(x, y) - a(x, y)) + \Pi_n [a_0(x, y) - a(x, y)]. $$

Therefore,

$$ \mathbb{E}[a_n(x, y)] - a(x, y) = \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k (\mathbb{E}[\mathcal{Z}_k(x, y)] - a(x, y)) + \Pi_n [a_0(x, y) - a(x, y)]. $$

In order to determine the bias of $a_n$, we need to simply focus on the quantity $\mathbb{E}[\mathcal{Z}_k(x, y)] - a(x, y)$. Relying upon the assumption (A1) and (A4), it follows that

$$ \begin{aligned} \mathbb{E}[\mathcal{Z}_k(x, y)] &= \int_{\mathbb{R}^{d+q}} h_k^{-d} \mathbf{K} \left( \frac{x-t}{h_k} \right) \chi_k(y) f_{(X,Y)}(t, u) dt du \\ &= \int_{\mathbb{R}^d} h_k^{-d} \mathbf{K} \left( \frac{x-t}{h_k} \right) \mathbb{E}[\chi_k(y) | X=t] f_X(t) dt \\ &= \int_{\mathbb{R}^d} h_k^{-d} \mathbf{K} \left( \frac{x-t}{h_k} \right) \pi(y|t) f_X(t) dt \\ &= \int_{\mathbb{R}^d} h_k^{-d} \mathbf{K} \left( \frac{x-t}{h_k} \right) a(t, y) dt \\ &= \int_{\mathbb{R}^d} \mathbf{K}(z) a(x - zh_k, y) dz. \end{aligned} $$

Moreover, Taylor's expansion with integral remainder ensures that

$$ \begin{aligned} \mathbb{E}[\mathcal{Z}_k(x, y)] - a(x, y) &= \int_{\mathbb{R}^d} \mathbf{K}(z) [a(x - zh_k, y) - a(x, y)] dz \\ &= \int_{\mathbb{R}^d} \mathbf{K}(z) \left[ \sum_{i=1}^d \frac{\partial a}{\partial x_i}(x, y) z_i h_k + \int_0^1 (1-t) \sum_{i,j=1}^d \frac{\partial^2 a}{\partial x_i \partial x_j} (x - tzh_k, y) z_i z_j h_k^2 dt \right] dz \\ &= h_k \sum_{i=1}^d a_i^{(1)}(x, y) \int_{\mathbb{R}^d} \mathbf{K}(z) z_i dz + h_k^2 \sum_{i,j=1}^d \int_{\mathbb{R}^d} \int_0^1 (1-t) a_{ij}^{(2)}(x - tzh_k, y) z_i z_j \mathbf{K}(z) dt dz. \end{aligned} $$
Referring to the assumption $(A_1)$, $\int_{\mathbb{R}^d} \mathbf{K}(z)z_i dz = 0$ and $\mu_i(\mathbf{K}) = \int_{\mathbb{R}^d} z_i^2 \mathbf{K}(z)dz$, it can be inferred that

$$ \mathbb{E}[\mathcal{Z}_k(x, y)] - a(x, y) = \frac{h_k^2}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{jj}^{(2)}(x, y) + h_k^2 \eta_k(x), $$

$$ \text{where, } \eta_k(x) := \sum_{i,j=1}^{d} \int_{\mathbb{R}^d} \int_{0}^{1} (1-t) [a_{ij}^{(2)}(x-tzh_k, y) - a_{ij}^{(2)}(x,y)] z_i z_j \mathbf{K}(z) dt dz. $$

Since we have $a_{ij}^{(2)}$ which is bounded and continuous at $x$ for all $i, j \in \{1, \dots, d\}$, we obtain $\lim_{k \to +\infty} \eta_k(x) = 0$, which ensures that, $\eta_k(x) = \mathbf{o}(1)$.

Hence,

$$ \begin{aligned} \mathbb{E}[a_n(x, y)] - a(x, y) &= \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k (\mathbb{E}[\mathcal{Z}_k(x, y)] - a(x, y)) + \Pi_n [a_0(x, y) - a(x, y)] \\ &= \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k \left( \frac{h_k^2}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{jj}^{(2)}(x, y) + h_k^2 \eta_k(x) \right) + \Pi_n [a_0(x, y) - a(x, y)] \\ &= \frac{1}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{jj}^{(2)}(x, y) \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k h_k^2 + \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k h_k^2 \eta_k(x) \\ &\quad + \Pi_n [a_0(x, y) - a(x, y)]. \end{aligned} $$

• For the case $a \le \alpha/(d+4)$, we have $\lim_{n \to +\infty} (n\gamma_n) > 2a$ and then $1 - 2a\xi > 0$. The application of lemma 1.2 enables us to write

$$ \begin{aligned} \mathbb{E}[a_n(x, y)] - a(x, y) &= \frac{1}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{jj}^{(2)}(x, y) \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k h_k^2 + \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k \mathbf{o}(h_k^2) + \mathbf{o}(\Pi_n) \\ &= \frac{1}{2(1 - 2a\xi)} \left( \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{jj}^{(2)}(x, y) \right) h_n^2 + \mathbf{o}(h_n^2) + \mathbf{o}(1) + \mathbf{o}(\Pi_n). \end{aligned} $$

We thus obtain the following desired result

$$ \mathbb{E}[a_n(x, y)] - a(x, y) = \frac{1}{2(1 - 2a\xi)} \left( \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{jj}^{(2)}(x, y) \right) h_n^2 + \mathbf{o}(h_n^2). $$

• For the case $a > \alpha/(d+4)$, we have $\lim_{n \to +\infty} (n\gamma_n) > \frac{\alpha - ad}{2}$ which yields that $h_n^2 = \mathbf{o}(\sqrt{\gamma_n h_n^{-d}})$. Then, the use of lemma 1.2 leads to

$$ \begin{aligned} \mathbb{E}[a_n(x, y)] - a(x, y) &= \frac{1}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{jj}^{(2)}(x, y) \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k \mathbf{o}(\sqrt{\gamma_k h_k^{-d}}) + \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k \mathbf{o}(\sqrt{\gamma_k h_k^{-d}}) \\ &\quad + \mathbf{o}(\Pi_n) \\ &= \mathbf{o}(\sqrt{\gamma_n h_n^{-d}}). \end{aligned} $$

Therefore, the claimed result (2.12) is established.
For the variance, and owing to the independence of $X_i$, for $i = 1, \dots, n$, it's obvious that

$$
\begin{align*}
\mathrm{Var}[a_n(x, y)] &= \mathrm{Var}[\Pi_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k \mathcal{Z}_k(x, y)] \\
&= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 \mathrm{Var}[\mathcal{Z}_k(x, y)] \\
&= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 (\mathbb{E}[\mathcal{Z}_k^2(x, y)] - \mathbb{E}[\mathcal{Z}_k(x, y)]^2).
\end{align*}
$$

$$
\begin{align*}
\mathbb{E}[\mathcal{Z}_k^2(x, y)] - \mathbb{E}[\mathcal{Z}_k(x, y)]^2 &= \int_{\mathbb{R}^d} h_k^{-2d} \mathbf{K}^2 \left(\frac{x-t}{h_k}\right) \mathbb{E}[\chi_k(y)^2 | X=t] f_X(t) dy \\
&\quad - \left(\int_{\mathbb{R}^d} h_k^{-d} \mathbf{K} \left(\frac{x-t}{h_k}\right) \pi(y|t) f_X(t) dt\right)^2 \\
&= \int_{\mathbb{R}^d} h_k^{-d} \mathbf{K}^2(z) a(x-zh_k, y) dz - \left(\int_{\mathbb{R}^d} \mathbf{K}(z) a(x-zh_k, y) dz\right)^2.
\end{align*}
$$

As matter of fact, the Taylor's expansions theorem ensures that

$$
\mathrm{Var}[a_n(x, y)] = \Pi_n^2 \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 h_k^{-d} \left[ a(x, y) \int_{\mathbb{R}^d} \mathbf{K}^2(z) dz + \nu_k(x) - h_k^d \eta_k(x) \right],
$$

where,

$$
\nu_k(x) = \int_{\mathbb{R}^d} \mathbf{K}^2(z) [a(x - zh_k, y) - a(x, y)] dz \text{ and } \eta_k(x) = \left( \int_{\mathbb{R}^d} \mathbf{K}(z) a(x - zh_k, y) dz \right)^2.
$$

• For the case $a < \alpha/(d+4)$, we have $\lim_{n \to +\infty} (n\gamma_n) > 2a$ which provides that $\gamma_n h_n^{-d} = \mathbf{o}(h_n^4)$.
By applying lemma 1.2, we infer that

$$
\begin{align*}
\mathrm{Var}[a_n(x, y)] &= \Pi_n^2 \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 h_k^{-d} [a(x, y) R(\mathbf{K}) + \mathbf{o}(1)] \\
&= \Pi_n^2 \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k \mathbf{o}(h_k^4).
\end{align*}
$$

Hence, we obtain the result

$$
\mathrm{Var}[a_n(x, y)] = \mathbf{o}(h_n^4).
$$

• For the case $a \ge \alpha/(d+4)$, we have $\lim_{n \to +\infty} (n\gamma_n) > \frac{\alpha-\alpha d}{2}$ and then $2-(\alpha-\alpha d)\xi > 0$. Since we have $\lim_{k \to +\infty} \nu_k(x) = 0$ and $\lim_{k \to +\infty} h_k \eta_k(x) = 0$, then the application of lemma 1.2 ensures that

$$
\begin{align*}
\mathrm{Var}[a_n(x, y)] &= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 h_k^{-d} \left[ a(x, y) R(\mathbf{K}) + \nu_k(x) - h_k^d \eta_k(x) \right] \\
&= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 h_k^{-d} \left[ a(x, y) R(\mathbf{K}) + \mathbf{o}(1) \right] \\
&= \frac{1}{2 - (\alpha - ad)\xi} \frac{\gamma_n}{h_n^d} \left[ a(x, y) R(\mathbf{K}) + \mathbf{o}(1) \right].
\end{align*}
$$

Thus, this leads to the result displayed in (2.14).

<figure><img src="image_34.png" /></figure>
Proof of Theorem 2.3.

Our proof rests upon the following decomposition, for $f_n \neq 0$, $n \ge 0$

$$ \pi_n(y|x) - \pi(y|x) = A_n(x, y) \frac{f_X(x)}{f_n(x)}, \qquad (2.27) $$

with

$$ A_n(x, y) = \frac{1}{f_X(x)} (a_n(x, y) - a(x, y)) - \frac{\pi(y|x)}{f_X(x)} (f_n(x) - f_X(x)). $$

It follows from (2.27) that the asymptotic behavior of $\pi_n(y|x) - \pi(y|x)$ can be deduced from the one of $A_n(x, y)$.

1. **Bias of $\pi_n$:** Here, we can state

$$ \mathbb{E}[A_n(x, y)] = \frac{1}{f_X(x)} (\mathbb{E}[a_n(x, y)] - a(x, y)) - \frac{\pi(y|x)}{f_X(x)} (\mathbb{E}[f_n(x)] - f_X(x)). $$

Now, using the first bias part of proposition 2.1 and proposition 2.2 and considering the fact that $a(x, y) = \pi(y|x)f_X(x)$, then by combining the assertions (2.7), (2.8), (2.11) and (2.12); we obtain the relations (2.15) and (2.16).

2. **Variance of $\pi_n$:** In order to confirm this statement, we have

$$ \begin{aligned} \mathrm{Var}[A_n(x, y)] = & \frac{1}{(f_X(x))^2} \mathrm{Var}[a_n(x, y)] + \frac{(\pi(y|x))^2}{(f_X(x))^2} \mathrm{Var}[f_n(x)] \\ & - 2 \frac{\pi(y|x)}{(f_X(x))^2} \mathrm{Cov}(a_n(x, y), f_n(x)). \end{aligned} \qquad (2.28) $$

Given that the $X_k$'s are independent, then for all $i = 1, \dots, n$, $k = 1, \dots, n$ and $i \neq k$, we have $\mathrm{Cov}(\mathcal{Z}_k(x, y), \mathcal{W}_i(x)) = 0$. Using lemma 1.2, classical computations entail

$$ \mathrm{Cov}(a_n(x, y), f_n(x)) = \frac{R(\mathbf{K})}{2 - (\alpha - ad)\xi} \pi(y|x)f_X(x) \frac{\gamma_n}{h_n^d} + \mathbf{o}\left(\frac{\gamma_n}{h_n^d}\right). \qquad (2.29) $$

In fact, we have

$$ \begin{aligned} \mathrm{Cov}(a_n(x, y), f_n(x)) &= \mathrm{Cov} \left( \Pi_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k \mathcal{Z}_k(x, y), \Pi_n \sum_{i=1}^n \Pi_i^{-1} \gamma_i \mathcal{W}_i(x) \right) \\ &= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 \mathrm{Cov}(\mathcal{Z}_k(x, y), \mathcal{W}_k(x)) \\ &= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 (\mathbb{E}[\mathcal{Z}_k(x, y) \mathcal{W}_k(x)] - \mathbb{E}[\mathcal{Z}_k(x, y)] \mathbb{E}[\mathcal{W}_k(x)]). \end{aligned} $$

Hence, the use of Taylor's expansion with integral remainder ensures that

$$ \mathbb{E}[\mathcal{Z}_k(x, y) \mathcal{W}_k(x)] = \mathbb{E} \left[ \chi_k(y) h_k^{-2d} \mathbf{K}^2 \left( \frac{x - X_k}{h_k} \right) \right] = R(\mathbf{K}) \pi(y|x) f_X(x) h_k^{-d} + \mathbf{o} \left( h_k^{-d} \right), $$

$$ \mathbb{E}[\mathcal{Z}_k(x, y)] = \mathbb{E} \left[ \chi_k(y) h_k^{-d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right) \right] = \pi(y|x) f_X(x) + \mathbf{o}(1), $$

$$ \text{and} \quad \mathbb{E}[\mathcal{W}_k(x)] = \mathbb{E} \left[ h_k^{-d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right) \right] = f_X(x) + \mathbf{o}(1). $$
Thus, by applying lemma 1.2, we can assert

$$
\begin{align*}
\mathrm{Cov}(a_n(x, y), f_n(x)) &= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 (R(\mathbf{K})\pi(y|x)f_X(x)h_k^{-d} - \pi(y|x)f_X^2(x) + \mathbf{o}(h_k^{-d}) + \mathbf{o}(1)) \\
&= \Pi_n^2 \sum_{k=1}^n \Pi_k^{-2} \gamma_k^2 h_k^{-d} (R(\mathbf{K})\pi(y|x)f_X(x) + \mathbf{o}(1)) \\
&= \frac{R(\mathbf{K})}{2 - (\alpha - ad)\xi} \pi(y|x)f_X(x) \frac{\gamma_n}{h_n^d} + \mathbf{o}\left(\frac{\gamma_n}{h_n^d}\right).
\end{align*}
$$

Consequently, relations (2.17) and (2.18) follow from the combination of assertions (2.9), (2.10), (2.13), (2.14) and (2.29).

For the case $a \ge \alpha/(d+4)$, we deduce with (2.28) that

$$
\mathrm{Var}[\pi_n(y|x)] = \frac{R(\mathbf{K})}{2 - (\alpha - ad)\xi} \frac{\pi(y|x)(1 - \pi(y|x))}{f_X(x)} \frac{\gamma_n}{h_n^d} + \mathbf{o}\left(\frac{\gamma_n}{h_n^d}\right).
$$

Proceeding with the same reasoning applied for the case $a < \alpha/(d+4)$, we obtain the desired result (2.17). $\square$

Proof of Theorem 2.4.

We can write, for all $n \ge 0$, $x \in \mathbb{R}^d$, $y \in \mathbb{R}^q$,

$$
\begin{align*}
A_n(x, y) - \mathbb{E}[A_n(x, y)]
&= \frac{1}{f_X(x)} [a_n(x, y) - \mathbb{E}[a_n(x, y)]] - \frac{\pi(y|x)}{f_X(x)} [f_n(x) - \mathbb{E}[f_n(x)]] \\
&= \frac{1}{f_X(x)} \Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k (T_k(x, y) - \mathbb{E}[T_k(x, y)]),
\end{align*}
$$

where

$$
T_k(x, y) = \mathcal{Z}_k(x, y) - \pi(y|x)W_k(x).
$$

Here and subsequently, it is worth noting

$$
S_k(x, y) := \Pi_k^{-1} \gamma_k (T_k(x, y) - \mathbb{E}[T_k(x, y)]).
$$

Hence, we state

$$
A_n(x, y) - \mathbb{E}[A_n(x, y)] = \frac{1}{f_X(x)} \Pi_n \sum_{k=1}^{n} S_k(x, y). \quad (2.30)
$$

This proof falls naturally into the application of Lyapunov’s theorem **1.14** for $S_k(x, y)$.
On the one hand, we can write

$$
\begin{align*}
u_n^2 &:= \sum_{k=1}^{n} \mathrm{Var}[S_k(x, y)] \\
&= \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 \mathrm{Var}[T_k(x, y)] \\
&= \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 \mathrm{Var}[\mathcal{Z}_k(x, y) - \pi(y|x)\mathcal{W}_k(x)] \\
&= \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 (\mathrm{Var}[\mathcal{Z}_k(x, y)] + \pi(y|x)^2 \mathrm{Var}[\mathcal{W}_k(x)]) \\
&\quad - 2\pi(y|x) \mathrm{Cov}(\mathcal{Z}_k(x, y), \mathcal{W}_k(x)). \tag*{[}
\end{align*}
$$
Moreover, we have

$$ \mathrm{Var}[\mathcal{Z}_k(x, y)] = h_k^{-d}(R(\mathbf{K}) \pi(y|x)f_X(x) + \mathbf{o}(1)), $$

$$ \mathrm{Var}[W_k(x)] = h_k^{-d}(R(\mathbf{K}) f_X(x) + \mathbf{o}(1)), $$

$$ \text{and} \quad \mathrm{Cov}(\mathcal{Z}_k(x, y), W_k(x)) = h_k^{-d}(R(\mathbf{K}) \pi(y|x)f_X(x) + \mathbf{o}(1)). $$

Therefore, by applying lemma 1.2, it can be deduced that

$$ u_n^2 = \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 h_k^{-d} (R(\mathbf{K}) f_X(x) \pi(y|x) (1 - \pi(y|x)) + \mathbf{o}(1)) = \frac{f_X(x)^2 \gamma_n}{\Pi_n^2 h_n^d} [\sigma^2(x, y) + \mathbf{o}(1)]. \quad (2.31) $$

On the other hand, we have

$$ \mathbb{E}[\mathcal{Z}_k^{2+p}(x, y)] = h_k^{-d(1+p)} \left[ \pi(y|x) f_X(x) \int_{\mathbb{R}^d} \mathbf{K}^{2+p}(z) dz + \mathbf{o}(1) \right] $$

$$ \text{and} \quad \mathbb{E}[W_k^{2+p}(x, y)] = h_k^{-d(1+p)} \left[ f_X(x) \int_{\mathbb{R}^d} \mathbf{K}^{2+p}(z) dz + \mathbf{o}(1) \right]. $$

Then, we can write

$$ \mathbb{E}[|T_k(x, y)|^{2+p}] = \mathbf{O} \left( \frac{1}{h_k^{d(1+p)}} \right), \quad \forall p > 0. $$

Therefore,

$$ \begin{aligned} \mathbb{E}[|S_k(x, y)|^{2+p}] &= \Pi_k^{-(2+p)} \gamma_k^{2+p} \mathbb{E}[|T_k(x, y) - \mathbb{E}[T_k(x, y)]|^{2+p}] \\ &= \Pi_k^{-(2+p)} \gamma_k^{2+p} \mathbb{E}[|T_k(x, y) + \mathbf{o}(1)|^{2+p}]. \end{aligned} $$

Hence,

$$ \mathbb{E}[|S_k(x, y)|^{2+p}] = \mathbf{O} \left( \Pi_k^{-(2+p)} \gamma_k^{2+p} \mathbb{E}[|T_k(x, y)|^{2+p}] \right). $$

This yields,

$$ \begin{aligned} \sum_{k=1}^{n} \mathbb{E}[|S_k(x, y)|^{2+p}] &= \mathbf{O} \left( \sum_{k=1}^{n} \Pi_k^{-(2+p)} \gamma_k^{2+p} \mathbb{E}[|T_k(x, y)|^{2+p}] \right) \\ &= \mathbf{O} \left( \sum_{k=1}^{n} \Pi_k^{-(2+p)} \gamma_k^{2+p} \frac{1}{h_k^{d(1+p)}} \right). \end{aligned} $$

For the application of lemma 1.2, let us assume that there exists a positive real $p$ such that

$$ \lim_{n \to +\infty} (n\gamma_n) > \frac{1+p}{2+p}(\alpha - ad). $$

Then, we obtain

$$ \sum_{k=1}^{n} \mathbb{E}[|S_k(x, y)|^{2+p}] = \mathbf{O} \left( \frac{\gamma_n^{1+p}}{\Pi_n^{2+p} h_n^{d(1+p)}} \right). $$

It follows that

$$ \frac{1}{u_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x, y)|^{2+p}] = \mathbf{O} \left( \frac{\gamma_n^{1+p}}{u_n^{2+p} \Pi_n^{2+p} h_n^{d(1+p)}} \right). $$
As a sequel, with the assertion (2.31), we can write

$$ \frac{1}{u_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x, y)|^{2+p}] = \mathbf{O}\left( \left( \frac{\gamma_n}{h_n^d} \right)^{p/2} \right) = \mathbf{o}(1). $$

Moreover, since we have

$$ \lim_{n \to +\infty} \frac{1}{u_n^{2+p}} \sum_{k=1}^{n} \mathbb{E} \left[ |S_k(x, y) - \mathbb{E}[S_k(x, y)]|^{2+p} \right] = \lim_{n \to +\infty} \frac{1}{u_n^{2+p}} \sum_{k=1}^{n} \mathbb{E} [ |S_k(x, y)|^{2+p} ] = 0, $$

then, by applying Lyapunov theorem 1.14, we get

$$ \frac{1}{u_n} \sum_{k=1}^{n} (S_k(x, y) - \mathbb{E}[S_k(x, y)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1). $$

This implies

$$ \frac{1}{u_n} \sum_{k=1}^{n} S_k(x, y) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1). $$

Additionally, the relations (2.27) and (2.30) ensure that

$$ \frac{1}{u_n \Pi_n} f_X(x) (\pi_n(y|x) - \mathbb{E}[\pi_n(y|x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1). \quad (2.32) $$

Given that

$$ u_n^2 = \frac{f_X(x)^2}{\Pi_n^2} \frac{\gamma_n}{h_n^d} [\sigma^2(x, y) + \mathbf{o}(1)], $$

and by replacing $u_n$ with its value in relation (2.32), we deduce that

$$ \sqrt{\gamma_n^{-1} h_n^{d+4}} (\pi_n(y|x) - \mathbb{E}[\pi_n(y|x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, \sigma^2(x, y)). \quad (2.33) $$

Since we have $\sqrt{\gamma_n^{-1} h_n^{d+4}} \xrightarrow{n \to +\infty} \sqrt{c}$, then the convergence (2.19) follows from the combination of relations (2.15), (2.16) and convergence (2.33). $\square$

### Proof of Proposition 2.5.

Based on the relation (2.20) and by distinguishing the different possible cases according to the expressions of the Bias ((2.15) and (2.16)) as well as the Variance ((2.17) and (2.18)), one can prove this proposition and find the required result.

In fact, we first note

$$ C_1 = \frac{I_1}{4(1 - 2a\xi)^2}, \quad C_2 = \frac{I_2}{2 - (\alpha - ad)\xi} R(\mathbf{K}). $$

$$ \text{MWISE}[\pi_n] = \begin{cases} C_1 h_n^4 + o(h_n^4) & \text{if } a \in (0, \frac{\alpha}{d+4}) \\ C_2 \gamma_n h_n^{-d} + C_1 h_n^4 + o(h_n^4) & \text{if } a = \frac{\alpha}{d+4} \\ C_2 \gamma_n h_n^{-d} + o(\gamma_n h_n^{-d}) & \text{if } a \in (\frac{\alpha}{d+4}, 1). \end{cases} $$

Set $\alpha \in ]1/2, 1]$.

* If $a < \frac{\alpha}{d+4}$, $C_1 h_n^4 \in \text{GS}(-4a)$ with $-\frac{4\alpha}{d+4} < -4a$.

* If $a = \frac{\alpha}{d+4}$, $C_2 \gamma_n h_n^{-d} + C_1 h_n^4 \in \text{GS}(-\frac{4\alpha}{d+4})$ with $-\frac{4\alpha}{d+4} = -4a$.

* If $a > \frac{\alpha}{d+4}$, $C_2 \gamma_n h_n^{-d} \in \text{GS}(-\alpha + ad)$ with $-\frac{4\alpha}{d+4} < -\alpha + ad$.

It follows that, for a given $\alpha$, to minimize the MWISE[$\pi_n$], we select the smallest quantity which is $-\frac{4\alpha}{d+4} = -4a$. Therefore, the parameter $a$ must be chosen equal to $\frac{\alpha}{d+4}$. $\square$
# Chapter 3

## The stochastic approximation method for semi-recursive multivariate kernel-type regression estimation

### 3.1 Introduction

Let $(X, Y)$ be a random vector with values in $\mathbb{R}^d \times \mathbb{R}$ with a joint density function $g(x, y)$ and let $f$ denote the probability density of $X$. Moreover, let $(X_1, Y_1), \dots, (X_n, Y_n)$ be independent random vectors identically distributed as $(X, Y)$. For a chosen measurable function $\varphi$ and $x \in \mathbb{R}^d$ the regression function, whenever it exists, is defined by

$$r_{\varphi}(x) = \mathbb{E}[\varphi(Y)|X = x] = \frac{1}{f(x)} \int_{\mathbb{R}} \varphi(y)g(x, y)dy.$$

Regression analysis stands for the study of how a response variable depends on one or more predictors. In fact, it's a reliable method of identifying which variables have impact on a topic of interest. The process of performing a regression allows us to confidently determine which factors matter most, which factors can be ignored, and how these factors influence each other. Regression problems can be usefully solved using nonparametric regression methods, which correspond to a category of regression analysis where the predictor does not take a predetermined form but is constructed according to information derived from the data. From this perspective, we have multiple methods of nonparametric estimation, such as Gaussian process regression (Kriging), kernel regression and regression trees.

An estimator of the $r_{\varphi}$ regression was first developed by Roussas (1990) and improved by Ein-mahl and Mason (2000) to determine exact rates of uniform strong consistency of kernel-type function estimators. Later, Deheuvels and Mason (2004) established uniform and non-uniform asymptotic simultaneous confidence bands for functionals of the distribution based on kernel-type estimators.

The classical recursive regression estimator was addressed in Mokaddem et al. (2009b) for univariate framework and a multivariate extension of this estimator was carried out by Mokaddem and Pelletier (2016). Subsequently, Slaoui (2016) established the semi-recursive case and introduced a new estimator which is the fraction of a recursive regression by a recursive density function. As far as this research is concerned, our basic objective is to extend this estimator for kernel-type estimation with large choice of parameters and properties in a multivariate case. Note that, recently, Bouzebda and Slaoui (2020) explored general kernel type estimators for censored data defined by the stochastic approximation algorithm.
### 3.1.1 Presentation of the method

Let us start with the presentation of our stochastic approximation method. The stochastic approximation algorithm, which estimates recursively the regression function

$$a_{\varphi}: x \mapsto r_{\varphi}(x)f(x) = \int_{\mathbb{R}} \varphi(y)g(x, y)dy$$

at a given vector $x$, can be expressed as follows:

$$a_{\varphi_n}(x) = (1 - \beta_n)a_{\varphi_{n-1}}(x) + \beta_n \varphi(Y_n) h_n^{-d} \mathbf{K} \left( \frac{x - X_n}{h_n} \right),$$

where $(\beta_n)$ is a positive sequences of real numbers decreasing towards zero verifying (18).

Here, we consider that $a_0(x) = 0$, then by a recurrence, we get

$$a_{\varphi_n}(x) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k \varphi(Y_k) h_k^{-d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right), \quad Q_n = \prod_{j=1}^{n} (1 - \beta_j).$$

Through this chapter, we consider the general multivariate kernel-type estimator for the regression function $r: x \mapsto \mathbb{E}[\varphi(Y)|X = x]$ at the vector $x$

$$r_{\varphi_n}(x) = \begin{cases} \frac{a_{\varphi_n}(x)}{f_n(x)} & \text{if } f_n(x) \neq 0, \\ 0 & \text{otherwise} \end{cases}, \qquad (3.1)$$

with $f_n$ stands for the recursive density estimator given in (2.3).

Our first aim is to examine the asymptotic properties of our proposed semi-recursive estimator of a multivariate regression function. Then, we prove its performance.

By introducing the non-recursive estimator of $a_{\varphi}$ given by

$$\tilde{a}_{\varphi_n}(x) = \frac{1}{nh_n^d} \sum_{k=1}^{n} \varphi(Y_k) \mathbf{K} \left( \frac{x - X_k}{h_n} \right),$$

we shall compare our estimator to the generalized non-recursive kernel regression estimator of Nadaraya-Watson Nadaraya (1964) and Watson (1964) $\tilde{r}_{\varphi_n}$ indicated by

$$\tilde{r}_{\varphi_n}(x) = \begin{cases} \frac{\tilde{a}_{\varphi_n}(x)}{\tilde{f}_n(x)} & \text{if } \tilde{f}_n(x) \neq 0, \\ 0 & \text{otherwise} \end{cases}, \qquad (3.2)$$

where $\tilde{f}_n$ was given in (2.5).

#### Particular cases:

1. For $\varphi(y) := \mathbb{I}(y) = y$, we have the classical regression function

$$r_{\mathbb{I}}(x) = \mathbb{E}[Y|X = x].$$

A recursive estimator of $r_{\mathbb{I}}$ was reported in Slaoui (2015).

2. For $\varphi(y) := \mathbb{I}(y) = y^m$, $m \in \mathbb{N}$, we have the conditional moments

$$r_{\mathbb{I}}(x) = \mathbb{E}[Y^m|X = x].$$

3. For $\varphi(y) := \chi_t(y) = \mathbb{1}_{\{y \le t\}}$, $t \in \mathbb{R}$, we have the conditional cumulative distribution function

$$r_{\chi_t}(y) = \pi(t|x) = \mathbb{P}[Y \le t | X = x].$$

A recursive estimator of $r_{\chi_t}$ was identified in Slama et al. (2021).
### 3.1.2 Notations and assumptions

For this section and under ($A_1$) and ($A_2$) given in 1.1.1, we provide the following notations and assumptions which will be intensively used for our theoretical main results.

(A5) (i) $(\beta_n)_{n \ge 1} \in \mathcal{GS}(-\beta)$, with $\beta \in (\frac{1}{2}, 1]$.

(ii) $(h_n)_{n \ge 1} \in \mathcal{GS}(-a)$, with $a \in (0, 1)$.

(iii) $\lim_{n \to +\infty} n\beta_n \in \left( \min\{2a, \frac{\beta - ad}{2} \}, +\infty \right]$.

(A6) (i) The functions $f$ and $a_\varphi$ are bounded and twice differentiable.

(ii) For all $i, j \in \{1, \dots, d\}$, $f_{ij}^{(2)} := \frac{\partial^2 f}{\partial x_i \partial x_j}$ and $a_{\varphi ij}^{(2)} := \frac{\partial^2 a_\varphi}{\partial x_i \partial x_j}$ are bounded and continuous at $x$.

(iii) For all $p > 0$, $s \mapsto \int_{\mathbb{R}} |\varphi(t)|^{2+p} g(s, t) dt$ is a bounded function.

(iv) The function $s \mapsto \int_{\mathbb{R}} \varphi(t)^2 g(s, t) dt$ is bounded and continuous at $s = x$.

#### Discussion of the assumptions:

It is to be noted that the assumption (A5)(iii) with regard to the limit of $(n\beta_n)$ as $n$ goes to infinity is quite common within the context of stochastic approximation algorithms. More specifically, the limit $\xi_\beta := \lim_{n \to +\infty} (n\beta_n)^{-1}$ is implied to be finite. Moreover, throughout the rest of this thesis manuscript, we shall use an other notation for $\lim_{n \to +\infty} (n\gamma_n)^{-1}$ denoted $\xi_\alpha$.

Furthermore, it is noteworthy that the assumptions given in (A6) are fulfilled under the following conditions:

- For every $t \in \mathbb{R}$, the function $g(., t)$ is twice continuously differentiable in $\mathbb{R}^d$.

- The functions $s \mapsto \int_{\mathbb{R}} g(s, t) dt$, $s \mapsto \int_{\mathbb{R}} \varphi(t) g(s, t) dt$ and $s \mapsto \int_{\mathbb{R}} \varphi(t)^2 g(s, t) dt$ are bounded and continuous at $s = x$.

- For $p > 0$, $s \mapsto \int_{\mathbb{R}} |\varphi(t)|^{2+p} g(s, t) dt$ is a bounded function.

- For all $i, j \in \{1, \dots, d\}$, $\int_{\mathbb{R}} \left| \frac{\partial g}{\partial s_i} (x, t) \right| dt < +\infty$ and $\int_{\mathbb{R}} |\varphi(t)| \left| \frac{\partial g}{\partial s_i} (x, t) \right| dt < +\infty$.

- For all $i, j \in \{1, \dots, d\}$, $s \mapsto \int_{\mathbb{R}} \frac{\partial^2 g}{\partial s_i \partial s_j} (s, t) dt$ and $s \mapsto \int_{\mathbb{R}} \varphi(t) \frac{\partial^2 g}{\partial s_i \partial s_j} (s, t) dt$ are bounded functions continuous at $s = x$.

Those conditions on the density of the couple $(X, Y)$ were applied in the non-recursive framework for the estimation of the regression function Nadaraya (1964), Watson (1964) and in the recursive framework Mokaddem et al. (2009b), Slaoui (2015).

Throughout this chapter, the following notations are used :

$$ \xi_{\alpha, \beta} = \xi_{\beta} \xi_{\alpha}^{-1} := \lim_{n \to +\infty} (\gamma_n \beta_n^{-1}) \quad \text{and} \quad \xi_{\beta, \alpha} = (\xi_{\alpha, \beta})^{-1} := \lim_{n \to +\infty} (\beta_n \gamma_n^{-1}). $$

## 3.2 Main results

In order to investigate the asymptotic properties of our estimator $r_{\varphi n}$, we need to first introduce the following proposition which provide the bias and the variance of $a_{\varphi n}$.
### 3.2.1 Bias and variance of $a_{\varphi_n}$

**Proposition 3.1.** Under the assumptions (A₁), (A₂), (A₅) and (A₆), we obtain

$$1. \text{ If } a \in (0, \frac{\beta}{d+4}], \text{ then} \\ \mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x) = \frac{h_n^2}{2(1-2a\xi_\beta)} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x) + o(h_n^2). \quad (3.3)$$

$$\text{If } a \in \left(\frac{\beta}{d+4}, 1\right), \text{ then} \\ \mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x) = o\left(\sqrt{\beta_n h_n^{-d}}\right). \quad (3.4)$$

$$2. \text{ If } a \in (0, \frac{\beta}{d+4}), \text{ then} \\ \mathrm{Var}[a_{\varphi_n}(x)] = o(h_n^4). \quad (3.5)$$

$$\text{If } a \in \left[\frac{\beta}{d+4}, 1\right), \text{ then} \\ \mathrm{Var}[a_{\varphi_n}(x)] = \frac{\beta_n \mathbb{E}[\varphi(Y)^2 | X=x]}{h_n^d (2-(\beta-ad)\xi_\beta)} f(x) R(\mathbf{K}) + o\left(\beta_n h_n^{-d}\right). \quad (3.6)$$

### 3.2.2 Bias and variance of $r_{\varphi_n}$

Our main result rests on the following theorem, which yields the bias and the variance of $r_{\varphi_n}$.

**Theorem 3.2.** Under the assumptions (A₁), (A₂), (A₅) and (A₆), we obtain

$$1. \text{ If } a \in \left(0, \frac{\min(\beta, \alpha)}{d+4}\right], \text{ then} \\ \mathbb{E}[r_{\varphi_n}(x)] - r_{\varphi}(x) = \frac{h_n^2}{f(x)} \left( \frac{\sum_{j=1}^d \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x)}{2(1-2a\xi_\beta)} - \frac{r_{\varphi}(x) \sum_{j=1}^d \mu_j(\mathbf{K}) f_{jj}^{(2)}(x)}{2(1-2a\xi_\alpha)} \right) + o(h_n^2). \quad (3.7)$$

$$\text{If } a \in \left(\frac{\min(\beta, \alpha)}{d+4}, 1\right), \text{ then} \\ \mathbb{E}[r_{\varphi_n}(x)] - r_{\varphi}(x) = o\left(\sqrt{\beta_n h_n^{-d}}\right) \mathbf{1}_{\{\beta \le \alpha\}} + o\left(\sqrt{\gamma_n h_n^{-d}}\right) \mathbf{1}_{\{\alpha < \beta\}}. \quad (3.8)$$

$$2. \text{ If } a \in (0, \frac{\min(\beta, \alpha)}{d+4}), \text{ then} \\ \mathrm{Var}[r_{\varphi_n}(x)] = o(h_n^4). \quad (3.9)$$

$$\text{If } a \in \left[\frac{\min(\beta, \alpha)}{d+4}, 1\right), \text{ then} \\ \mathrm{Var}[r_{\varphi_n}(x)] = \frac{\beta_n R(\mathbf{K})}{h_n^d f(x)} \left[ \frac{\mathbb{E}[\varphi(Y)^2 | X=x]}{2-(\beta-ad)\xi_\beta} \right. \\ \left. -r_{\varphi}^2(x) \left( \frac{2}{1-(\beta-ad-\xi_\beta^{-1})\xi_\alpha} - \frac{\xi_{\alpha,\beta}}{2-(\alpha-ad)\xi_\alpha} \right) + o\left(\frac{\beta_n}{h_n^d}\right) \mathbf{1}_{\{\beta \le \alpha\}} \right. \\ \left. + \frac{\gamma_n R(\mathbf{K})}{h_n^d f(x)} \left[ \frac{\mathbb{E}[\varphi(Y)^2 | X=x]\xi_{\beta,\alpha}}{2-(\beta-ad)\xi_\beta} \right. \right. \\ \left. \left. -r_{\varphi}^2(x) \left( \frac{2}{1-(\alpha-ad-\xi_\alpha^{-1})\xi_\beta} - \frac{1}{2-(\alpha-ad)\xi_\alpha} \right) + o\left(\frac{\gamma_n}{h_n^d}\right) \mathbf{1}_{\{\alpha < \beta\}} \right] \right). \quad (3.10)$$
Therefore, the bias and the variance of the estimator $r_{\varphi_n}$ defined by the stochastic approximation algorithm (3.1) depend heavily on the choice of the stepsizes ($\beta_n$) and ($\gamma_n$).

**Remark 3.3.** Notice that, for the case where ($\gamma_n$) = ($\beta_n$) and then $\alpha = \beta$, the expression (3.10) will be written as follows

$$ \mathrm{Var}[r_{\varphi_n}(x)] = \frac{\beta_n R(\mathbf{K})}{h_n^d f(x)} \frac{\mathrm{Var}[\varphi(Y)|X=x]}{2 - (\beta - ad)\xi_\beta} + o(\beta_n h_n^{-d}). $$

The asymptotic normality of the generalized semi-recursive estimator $r_{\varphi_n}$ is indicated by the following theorem.

### 3.2.3 Weak pointwise convergence rate of $r_{\varphi_n}$

**Theorem 3.4.** Under the assumptions (A$_1$), (A$_2$), (A$_5$) and (A$_6$), we obtain:

1. For the case $\beta \le \alpha$:

(a) If there exists $c \ge 0$ such that $\beta_n^{-1}h_n^{d+4} \xrightarrow{n \to +\infty} c$, then

$$ \sqrt{\beta_n^{-1} h_n^d} (r_{\varphi_n}(x) - r_{\varphi}(x)) \xrightarrow{n \to +\infty} \mathcal{N}\left(\sqrt{c}\mathbf{M}_{\beta}(x), \mathbf{\Sigma}_{\beta}(x)\right), \quad (3.11) $$

with

$$ \mathbf{M}_{\beta}(x) = \frac{1}{2f(x)} \left( \frac{\sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x)}{(1 - 2a\xi_{\beta})} - \frac{r_{\varphi}(x) \sum_{j=1}^{d} \mu_j(\mathbf{K}) f_{jj}^{(2)}(x)}{(1 - 2a\xi_{\alpha})} \right) \quad (3.12) $$

and

$$ \mathbf{\Sigma}_{\beta}(x) = \frac{R(\mathbf{K})}{f(x)} \left[ \frac{\mathbb{E}[\varphi(Y)^2 | X=x]}{2 - (\beta - ad)\xi_{\beta}} - r_{\varphi}^2(x) \left( \frac{2}{1 - (\beta - ad - \xi_{\beta}^{-1})\xi_{\alpha}} - \frac{\xi_{\alpha, \beta}}{2 - (\alpha - ad)\xi_{\alpha}} \right) \right]. \quad (3.13) $$

(b) If $\beta_n^{-1}h_n^{d+4} \xrightarrow{n \to +\infty} +\infty$, then

$$ \frac{1}{h_n^2} (r_{\varphi_n}(x) - r_{\varphi}(x)) \xrightarrow{n \to +\infty} \mathbf{M}_{\beta}(x). $$

2. For the case $\beta > \alpha$:

(a) If there exists $c \ge 0$ such that $\gamma_n^{-1}h_n^{d+4} \xrightarrow{n \to +\infty} c$, then

$$ \sqrt{\gamma_n^{-1} h_n^d} (r_{\gamma_n}(x) - r_{\gamma}(x)) \xrightarrow{n \to +\infty} \mathcal{N}\left(\sqrt{c}\mathbf{M}_{\gamma}(x), \mathbf{\Sigma}_{\gamma}(x)\right), \quad (3.14) $$

with

$$ \mathbf{M}_{\gamma}(x) = \frac{1}{2f(x)} \left( \frac{\sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x)}{(1 - 2a\xi_{\beta})} - \frac{r_{\varphi}(x) \sum_{j=1}^{d} \mu_j(\mathbf{K}) f_{jj}^{(2)}(x)}{(1 - 2a\xi_{\alpha})} \right) \quad (3.15) $$
and

$$ \Sigma_{\gamma}(x) = \frac{R(\mathbf{K})}{f(x)} \left[ \frac{\mathbb{E}[\varphi(Y)^2 | X = x] \xi_{\beta, \alpha}}{2 - (\alpha - ad)\xi_{\alpha}} - r_{\varphi}^2(x) \left( \frac{2}{1 - (\alpha - ad - \xi_{\alpha}^{-1})\xi_{\beta}} - \frac{1}{2 - (\alpha - ad)\xi_{\alpha}} \right) \right]. \quad (3.16) $$

(b) If $\gamma_n^{-1}h_n^{d+4} \xrightarrow{n \to +\infty} +\infty$, then

$$ \frac{1}{h_n^2} (r_{\varphi_n}(x) - r_{\varphi}(x)) \xrightarrow{n \to +\infty} \mathbf{M}_{\gamma}(x). $$

The following theorem demonstrates the strong pointwise convergence rate of our estimator $r_{\varphi_n}$.

### 3.2.4 Strong pointwise convergence rate of $r_{\varphi_n}$

**Theorem 3.5.** Under the assumptions (A$_1$), (A$_2$), (A$_5$) and (A$_6$), we get:

1. For the case $\beta \le \alpha$:

(a) If there exists $b \ge 0$ such that $\frac{\beta_n^{-1}h_n^{d+4}}{\ln\left(\sum_{i=1}^n \beta_i\right)} \xrightarrow{n \to +\infty} b$, then with probability one, the sequence

$$ \left( \sqrt{\frac{\beta_n^{-1}h_n^d}{2 \ln\left(\sum_{i=1}^n \beta_i\right)}} (r_{\varphi_n}(x) - r_{\varphi}(x)) \right) $$

is relatively compact and its limit set is the interval

$$ \left[ \sqrt{\frac{b}{2}} \mathbf{M}_{\beta}(x) - \sqrt{\mathbf{\Sigma}_{\beta}(x)}, \sqrt{\frac{b}{2}} \mathbf{M}_{\beta}(x) + \sqrt{\mathbf{\Sigma}_{\beta}(x)} \right]. $$

(b) If $\frac{\beta_n^{-1}h_n^{d+4}}{\ln\left(\sum_{i=1}^n \beta_i\right)} \xrightarrow{n \to +\infty} +\infty$, then, with probability one,

$$ \lim_{n \to +\infty} \frac{1}{h_n^2} (r_{\varphi_n}(x) - r_{\varphi}(x)) = \mathbf{M}_{\beta}(x). $$

2. For the case $\beta > \alpha$:

(a) If there exists $b \ge 0$ such that $\frac{\gamma_n^{-1}h_n^{d+4}}{\ln\left(\sum_{i=1}^n \gamma_i\right)} \xrightarrow{n \to +\infty} b$, then with probability one, the sequence

$$ \left( \sqrt{\frac{\gamma_n^{-1}h_n^d}{2 \ln\left(\sum_{i=1}^n \gamma_i\right)}} (r_{\varphi_n}(x) - r_{\varphi}(x)) \right) $$

is relatively compact and its limit set is the interval

$$ \left[ \sqrt{\frac{b}{2}} \mathbf{M}_{\gamma}(x) - \sqrt{\mathbf{\Sigma}_{\gamma}(x)}, \sqrt{\frac{b}{2}} \mathbf{M}_{\gamma}(x) + \sqrt{\mathbf{\Sigma}_{\gamma}(x)} \right]. $$
$$
(b) \text{ If } \frac{\gamma_n^{-1} h_n^{d+4}}{\ln\left(\sum_{i=1}^{n} \gamma_i\right)} \xrightarrow{n \to +\infty} +\infty, \text{ then, with probability one,}
$$

$$
\lim_{n \to +\infty} \frac{1}{h_n^2} (r_{\varphi_n}(x) - r_{\varphi}(x)) = \mathbf{M}_{\gamma}(x).
$$

In what follows, we clarify the choices of the stepsizes (β_n) as well as (γ_n) and the bandwidth (h_n)
based on the *MWISE* of the recursive estimator minimization, and then enact a comparison
with Nadaraya Watson's estimator.

3.3 Optimal choice of the stepsizes

In order to measure the optimal choice of the couple of stepsizes (β_n, γ_n), we need to minimize
the Mean Weighted Integrated Squared Error (MWISE) of the semi-recursive estimator r_φn.
For the sequel, we will need the following notations. We first note,

$$
I_1 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^d \mu_j(\mathbf{K}) a_{\varphi jj}^{(2)}(x) \right)^2 f(x)dx, \quad I_4 = \int_{\mathbb{R}^d} \mathbb{E}[\varphi(Y)^2 | X=x] f^2(x)dx,
$$

$$
I_2 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi jj}^{(2)}(x) \right) \left( \sum_{j=1}^{d} \mu_j(\mathbf{K}) f_{jj}^{(2)}(x) \right) r_{\varphi}(x) f(x) dx,
$$

$$
I_3 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^{d} \mu_j(\mathbf{K}) f_{jj}^{(2)}(x) \right)^2 r_{\varphi}^2(x) f(x) dx, \quad I_5 = \int_{\mathbb{R}^d} r_{\varphi}^2(x) f^2(x) dx.
$$

3.3.1 Asymptotic expressions of MWISE [$r_{\varphi_n}$]

The *MWISE* of the estimator $r_{\varphi_n}$ is determined by the following expression

$$
\text{MWISE}[r_{\varphi_n}] = \int_{\mathbb{R}^d} (\mathbb{E}[r_{\varphi_n}(x)] - r_{\varphi}(x))^2 f^3(x)dx + \int_{\mathbb{R}^d} \text{Var}[r_{\varphi_n}(x)]f^3(x)dx.
$$

**Proposition 3.6.** We first note,

$$
C_1 = \frac{I_1}{(1 - 2a\xi_\beta)^2} - \frac{2I_2}{(1 - 2a\xi_\beta)(1 - 2a\xi_\alpha)} + \frac{I_3}{(1 - 2a\xi_\alpha)^2}, \quad C_2 = \frac{I_4}{2 - (\beta - ad)\xi_\beta}
$$

$$
C_3 = \frac{I_4 \xi_{\beta, \alpha}}{2 - (\beta - ad)\xi_{\beta}}, \quad C_4 = I_5 \left( \frac{2}{1 - (\beta - ad - \xi_{\beta}^{-1})\xi_{\alpha}} - \frac{\xi_{\alpha, \beta}}{2 - (\alpha - ad)\xi_{\alpha}} \right),
$$

$$
C_5 = I_5 \left( \frac{2}{1 - (\alpha - ad - \xi_{\alpha}^{-1})\xi_{\beta}} - \frac{1}{2 - (\alpha - ad)\xi_{\alpha}} \right).
$$

1. For the case $\beta \le \alpha$:

$$
MWISE[r_{\varphi_n}] = \begin{cases} \frac{1}{4}C_1 h_n^4 + o(h_n^4) & \text{if } a \in (0, \frac{\beta}{d+4}) \\ (C_2 - C_4)R(\mathbf{K})\beta_n h_n^{-d} + \frac{1}{4}C_1 h_n^4 + o(h_n^4) & \text{if } a = \frac{\beta}{d+4} \\ (C_2 - C_4)R(\mathbf{K})\beta_n h_n^{-d} + o(\beta_n h_n^{-d}) & \text{if } a \in (\frac{\beta}{d+4}, 1) \end{cases} .
$$
2. For the case $\beta > \alpha$:

$$ MWISE[r_{\varphi_n}] = \begin{cases} \frac{1}{4}C_1h_n^4 + o(h_n^4) & \text{if } a \in (0, \frac{\alpha}{d+4}) \\ (C_3 - C_5)R(\mathbf{K})\gamma_n h_n^{-d} + \frac{1}{4}C_1h_n^4 + o(h_n^4) & \text{if } a = \frac{\alpha}{d+4} \\ (C_3 - C_5)R(\mathbf{K})\gamma_n h_n^{-d} + o(\gamma_n h_n^{-d}) & \text{if } a \in (\frac{\alpha}{d+4}, 1) \end{cases}. $$

The following corollary ensures that the bandwidth which minimizes the MWISE of $r_{\varphi_n}$ depends on the choice of the stepsizes $(\beta_n)$ and $(\gamma_n)$ and then the corresponding MWISE depends in turn on $(\beta_n)$ and $(\gamma_n)$.

**Corollary 3.7.** Let assumptions (A$_1$), (A$_2$), (A$_5$) and (A$_6$) hold. To minimize the MWISE of $r_{\varphi_n}$, the bandwidth ($h_n$) needs to be equal to the following expressions.

1. For the case $\beta \le \alpha$:

$$ h_n = d^{\frac{1}{d+4}} \left( \frac{C_2 - C_4}{C_1} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} \beta_n^{\frac{1}{d+4}}. $$

Hence, the corresponding MWISE is specified by

$$ MWISE[r_{\varphi_n}] = \frac{(d+4)}{4d^{\frac{d}{d+4}}} C_1^{\frac{d}{d+4}} (C_2 - C_4)^{\frac{4}{d+4}} R(\mathbf{K})^{\frac{4}{d+4}} \beta_n^{\frac{4}{d+4}} + o\left(\beta_n^{\frac{4}{d+4}}\right). $$

2. For the case $\alpha < \beta$:

$$ h_n = d^{\frac{1}{d+4}} \left( \frac{C_3 - C_5}{C_1} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} \gamma_n^{\frac{1}{d+4}}. $$

Thus, the corresponding MWISE is expressed by

$$ MWISE[r_{\varphi_n}] = \frac{(d+4)}{4d^{\frac{d}{d+4}}} C_1^{\frac{d}{d+4}} (C_3 - C_5)^{\frac{4}{d+4}} R(\mathbf{K})^{\frac{4}{d+4}} \gamma_n^{\frac{4}{d+4}} + o\left(\gamma_n^{\frac{4}{d+4}}\right). $$

The following corollary is provided in the special cases, where $(\beta_n)$ is chosen as $(\beta_n) = (\beta_0 n^{-1})$ in order to minimize the $MWISE[a_{\varphi_n}]$ and $(\gamma_n)$ is chosen as $(\gamma_n) = (\gamma_0 n^{-1})$ in order to minimize the $MWISE[f_n]$.

**Proposition 3.8.** Let assumptions (A$_1$), (A$_2$), (A$_5$) and (A$_6$) hold. It is worth noting,

$$ \Theta_1 = \frac{\beta_0 I_4}{(d+4)\beta_0 - 2} - \frac{2\gamma_0 I_5}{(d+4)(\frac{\gamma_0+\beta_0}{2}) - 2} + \frac{\beta_0^{-1}\gamma_0^2 I_5}{(d+4)\gamma_0 - 2}, $$

$$ \Theta_2 = \frac{\gamma_0 I_5}{(d+4)\gamma_0 - 2} - \frac{2\beta_0 I_5}{(d+4)(\frac{\gamma_0+\beta_0}{2}) - 2}, $$

$$ \Theta_3 = \frac{\beta_0^2 I_1}{((d+4)\beta_0 - 2)^2} - \frac{2\beta_0\gamma_0 I_2}{((d+4)\beta_0 - 2)((d+4)\gamma_0 - 2)} + \frac{\gamma_0^2 I_3}{((d+4)\gamma_0 - 2)^2}. $$

To minimize the MWISE of $r_{\varphi_n}$, we need to choose the stepsize $(\gamma_n)$ in $GS(-1)$ such that $\lim_{n \to +\infty} n\gamma_n = \gamma_0$ and the stepsize $(\beta_n)$ in $GS(-1)$ such that $\lim_{n \to +\infty} n\beta_n = \beta_0$. As a matter of fact,
1. For the case $\beta \le \alpha$:

The bandwidth ($h_n$) needs to be equal to

$$ \left( \beta_0^{\frac{1}{d+4}} \left( \frac{d}{2(d+4)} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-1}{d+4}} \left( \frac{\Theta_1}{\Theta_3} \right)^{\frac{1}{d+4}} \right). $$

Consequently, the corresponding MWISE is determined by

$$ MWISE[r_{\varphi_n}] = \frac{(d+4)^{\frac{3d+8}{d+4}}}{4^{\frac{d+6}{d+4}} d^{\frac{d}{d+4}}} \beta_0^{\frac{4}{d+4}} \Theta_1^{\frac{4}{d+4}} \Theta_3^{\frac{d}{d+4}} R(\mathbf{K})^{\frac{4}{d+4}} n^{\frac{-4}{d+4}} + o(n^{\frac{-4}{d+4}}). $$

2. For the case $\alpha < \beta$:

The bandwidth ($h_n$) needs to be equal to

$$ \left( \gamma_0^{\frac{1}{d+4}} \left( \frac{d}{2(d+4)} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-1}{d+4}} \left( \frac{\Theta_2}{\Theta_3} \right)^{\frac{1}{d+4}} \right). $$

Therefore, the corresponding MWISE is specified by

$$ MWISE[r_{\varphi_n}] = \frac{(d+4)^{\frac{3d+8}{d+4}}}{4^{\frac{d+6}{d+4}} d^{\frac{d}{d+4}}} \gamma_0^{\frac{4}{d+4}} \Theta_2^{\frac{4}{d+4}} \Theta_3^{\frac{d}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-4}{d+4}} + o(n^{\frac{-4}{d+4}}). $$

Additionally, the minimum of $MWISE[r_{\varphi_n}]$ is achieved at $(\beta_0, \gamma_0) = (1, 1)$. From this perspective, the optimal bandwidth ($h_n$) must be equal to

$$ \left( \left( \frac{d(d+2)}{2(d+4)} \right)^{\frac{1}{d+4}} \left( \frac{I_4 - I_5}{I_1 - 2I_2 + I_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-1}{d+4}} \right). \quad (3.17) $$

Thus, the corresponding MWISE is indicated by

$$ MWISE[r_{\varphi_n}] = \frac{(d+4)^{\frac{3d+8}{d+4}}}{4^{\frac{d+6}{d+4}} d^{\frac{d}{d+4}} (d+2)^{\frac{d+6}{d+4}}} (I_1 - 2I_2 + I_3)^{\frac{d}{d+4}} (I_4 - I_5)^{\frac{4}{d+4}} R(\mathbf{K})^{\frac{4}{d+4}} n^{\frac{-4}{d+4}} + o(n^{\frac{-4}{d+4}}). $$

**Remark 3.9.** Note that for the particular case where the stepsize $(\beta_n)$ is in $GS(-1)$ such that $\lim_{n \to +\infty} n\beta_n = 1$, $(\gamma_n)$ is in $GS(-1)$ such that $\lim_{n \to +\infty} n\gamma_n = 1$ and the bandwidth ($h_n$) is chosen such that $\lim_{n \to +\infty} nh_n^{d+4} = 0$ (which corresponds to undersmoothing), the asymptotic normality of the proposed estimator is represented as follows

$$ \sqrt{nh_n^d} (r_{\varphi_n}(x) - r_{\varphi}(x)) \xrightarrow{n \to +\infty} \mathcal{N} \left( 0, \frac{1}{a+d} R(\mathbf{K}) \frac{\mathrm{Var}[\varphi(Y)|X=x]}{f(x)} \right). $$

## 3.4 Asymptotic properties of $\tilde{r}_{\varphi_n}$

The main properties of the generalized non-recursive regression function estimator $\tilde{r}_{\varphi_n}$ are displayed in the following proposition.

**Proposition 3.10.** Let assumptions (A1) and (A6) hold. Then, the asymptotic properties of Nadaraya-Watson's estimator are denoted as follows.

* The bias of $\tilde{r}_{\varphi_n}$:

$$ \mathbb{E}[\tilde{r}_{\varphi_n}(x)] - r_{\varphi}(x) = \frac{1}{2f(x)} h_n^2 \left( \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi jj}^{(2)}(x) - r_{\varphi}(x) \sum_{j=1}^{d} \mu_j(\mathbf{K}) f_{jj}^{(2)}(x) \right) + o(h_n^2). $$
• The variance of $\tilde{r}_{\varphi_n}$:

$$Var[\tilde{r}_{\varphi_n}(x)] = \frac{1}{nh_n^d f(x)} Var[\varphi(Y)|X=x] R(\mathbf{K}) + o\left(\frac{1}{nh_n^d}\right).$$

• The MWISE of $\tilde{r}_{\varphi_n}$:

$$MWISE[\tilde{r}_{\varphi_n}] = \frac{1}{4}(I_1 - 2I_2 + I_3)h_n^4 + \frac{1}{nh_n^d}(I_4 - I_5)R(\mathbf{K}) + o\left(h_n^4 + \frac{1}{nh_n^d}\right).$$

To minimize the MWISE of $\tilde{r}_{\varphi_n}$, the bandwidth ($h_n$) must be equal to

$$\left(d^{\frac{1}{d+4}} \left( \frac{I_4 - I_5}{I_1 - 2I_2 + I_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{1}{d+4}}\right). \quad (3.18)$$

Therefore, the corresponding MWISE is expressed by

$$MWISE[\tilde{r}_{\varphi_n}] = \frac{(d+4)}{4d^{\frac{d}{d+4}}} (I_4 - I_5)^{\frac{4}{d+4}} (I_1 - 2I_2 + I_3)^{\frac{d}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{4}{d+4}} + o\left(n^{-\frac{4}{d+4}}\right).$$

• The asymptotic normality of $\tilde{r}_{\varphi_n}$:

Suppose that $nh_n^{d+4} \xrightarrow{n \to +\infty} 0$. Thus,

$$\sqrt{nh_n^d} (\tilde{r}_{\varphi_n}(x) - r_{\varphi}(x)) \xrightarrow{n \to +\infty} \mathcal{N}\left(0, R(\mathbf{K}) \frac{Var[\varphi(Y)|X=x]}{f(x)}\right).$$

• The weak pointwise convergence rate of $\tilde{r}_{\varphi_n}$:

If $nh_n^{d+4} \xrightarrow{n \to +\infty} +\infty$, then

$$\frac{1}{h_n^2} (\tilde{r}_{\varphi_n}(x) - r_{\varphi}(x)) \xrightarrow{n \to +\infty} 0.$$

## 3.5 Bandwidth selection

Kernel smoothing in non-parametric statistics requires the choice of a bandwidth parameter. There are numerous methods for bandwidth selection, namely the cross-validation method, the bootstrap procedure and the second generation plug-in approach.

First of all, we adopt (1.14) and the kernel $\mathbf{K}$ we shall use is considered as a product of univariate kernels. Hence, we have

$$r_{\varphi_n}(x) = \frac{Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k \varphi(Y_k) h_k^{-d} \prod_{i=1}^{d} K\left(\frac{x_i - X_{ki}}{h_k}\right)}{\Pi_n \sum_{k=1}^{n} \Pi_k^{-1} \gamma_k h_k^{-d} \prod_{i=1}^{d} K\left(\frac{x_i - X_{ki}}{h_k}\right)}$$

and

$$\tilde{r}_{\varphi_n}(x) = \frac{(nh_n^d)^{-1} \sum_{k=1}^{n} \varphi(Y_k) \prod_{i=1}^{d} K\left(\frac{x_i - X_{ki}}{h_n}\right)}{(nh_n^d)^{-1} \sum_{k=1}^{n} \prod_{i=1}^{d} K\left(\frac{x_i - X_{ki}}{h_n}\right)}.$$

Let us start by introducing our bandwidth selection methods.
### 3.5.1 Plug-in method

In statistics, Altman and Leger (1995) set forward an efficient method of bandwidth selection, a plug-in estimate which minimizes an estimate of the mean weighted integrated squared error, using the density function as a weight function. Since the MWISE depends on the unknown quantities $I_j$, $j = 1, \dots, 5$, we attempt to construct an asymptotic unbiased estimator of those quantities.

For this purpose, we let $\mu(K) = \int_{\mathbb{R}} z^2 K(z) dz$ and

$$I_i = \mu^2(K) I'_i, i = 1, 2, 3,$$

where

$$I'_1 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^d a_{\varphi_{jj}}^{(2)}(x) \right)^2 f(x) dx.$$

$$I'_2 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^d a_{\varphi_{jj}}^{(2)}(x) \right) \left( \sum_{j=1}^d f_{jj}^{(2)}(x) \right) r_{\varphi}(x) f(x) dx.$$

$$I'_3 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^d f_{jj}^{(2)}(x) \right)^2 r_{\varphi}(x) f(x) dx.$$

At this stage of analysis, in order to estimate the optimal bandwidth (3.17), we need to estimate $I_j$, $j = 1, \dots, 5$. For this purpose, we recall that $K_b$ is a kernel and $b_n$ is the associated bandwidth, such that $\delta = 2/5$, and $K_{b'}^{(2)}$ is the second derivative of a kernel $K_{b'}$ with the associated bandwidth $b'_n$ such that $\delta = 3/14$.

#### Semi-Recursive estimator $r_{\varphi_n}$:

To estimate the optimal bandwidth (3.17), we need to estimate $I_j$, $j = 1, \dots, 5$.

**Estimation of $I_1$, $I_2$ and $I_3$:** Here, the plug-in estimate gives

$$\hat{I}'_1 = \frac{Q_n^2}{n} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k}}^n Q_j^{-1} Q_k^{-1} \beta_j \beta_k b_j^{-(d+2)} b_k^{-(d+2)} \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_j'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_j} \right) \right] \\ \times \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_k'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \right] \varphi(Y_j) \varphi(Y_k),$$

$$\hat{I}'_2 = \frac{Q_n \Pi_n}{n} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k}}^n Q_j^{-1} \Pi_k^{-1} \beta_j \gamma_k b_j^{-(d+2)} b_k^{-(d+2)} \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_j'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_j} \right) \right] \\ \times \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_k'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \right] \varphi(Y_j) \varphi(Y_i),$$

$$\hat{I}'_3 = \frac{\Pi_n^2}{n} \sum_{\substack{i,j,k,m=1 \\ i \neq j \neq k \neq m}}^n \Pi_j^{-1} \Pi_k^{-1} \gamma_j \gamma_k b_j^{-(d+2)} b_k^{-(d+2)} \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_j'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_j} \right) \right] \\ \times \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_k'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \right] \varphi(Y_i) \varphi(Y_m),$$
Therefore, we obtain

$$ \hat{I}_i = \mu^2(\mathbf{K}) \hat{I}'_i, \quad i = 1 \dots 3. $$

**Estimation of $I_4$ and $I_5$:**

$$ \hat{I}_4 = \frac{\Pi_n}{n} \sum_{\substack{i,k=1 \\ i \neq k}}^{n} \Pi_k^{-1} \gamma_k b_k^{-d} \prod_{l=1}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \varphi(Y_i)^2 $$

and

$$ \hat{I}_5 = \frac{Q_n}{n} \sum_{\substack{i,k=1 \\ i \neq k}}^{n} Q_k^{-1} \beta_k b_k^{-d} \prod_{l=1}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \varphi(Y_i) \varphi(Y_k), $$

As a result, the plug-in estimator of (3.17) is denoted in terms of:

$$ h_n = \left( \left( \frac{d(d+2)}{2(d+4)} \right)^{\frac{1}{d+4}} \left( \frac{\hat{I}_4 - \hat{I}_5}{\hat{I}_1 - 2\hat{I}_2 + \hat{I}_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-1}{d+4}} \right), \quad (3.19) $$

Finally, an estimator of $MWISE[r_{\varphi_n}]$ is expressed as

$$ \widehat{MWISE}[r_{\varphi_n}] = \frac{(d+4)^{\frac{3d+8}{d+4}}}{4^{\frac{d+6}{d+4}} d^{\frac{d+6}{d+4}} (d+2)^{\frac{d}{d+4}}} (\hat{I}_1 - 2\hat{I}_2 + \hat{I}_3)^{\frac{d}{d+4}} (\hat{I}_4 - \hat{I}_5)^{\frac{4}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-4}{d+4}} + o(n^{\frac{-4}{d+4}}). $$

**Non-Recursive estimator $\tilde{r}_{\varphi_n}$:**

To estimate the optimal bandwidth (3.18), we need to estimate $I_j, j = 1, \dots, 5$.

**Estimation of $I_1, I_2$ and $I_3$:** For the non-recursive case, the plug-in estimate yields

$$ \tilde{I}'_1 = \frac{1}{n^3 b_n^{2(d+2)}} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k}}^{n} \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] \varphi(Y_j)\varphi(Y_k), $$

$$ \tilde{I}'_2 = \frac{1}{n^3 b_n^{2(d+2)}} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k \neq m}}^{n} \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] \varphi(Y_j)\varphi(Y_i), $$

$$ \tilde{I}'_3 = \frac{1}{n^4 b_n^{2(d+2)}} \sum_{\substack{i,j,k,m=1 \\ i \neq j \neq k \neq m}}^{n} \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] \varphi(Y_i)\varphi(Y_m), $$

Therefore, we obtain

$$ \tilde{I}_i = \mu^2(\mathbf{K}) \tilde{I}'_i, \quad i = 1 \dots 3. $$
Estimation of $I_4$ and $I_5$:

$$ \tilde{I}_4 = \frac{1}{n^2 b_n^d} \sum_{i,k=1}^{n} \prod_{l=1}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \varphi(Y_i)^2 $$

and

$$ \tilde{I}_5 = \frac{1}{n^2 b_n^d} \sum_{\substack{i,k=1 \\ i \neq k}}^{n} \prod_{l=1}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \varphi(Y_i) \varphi(Y_k), $$

As a consequence, the plug-in estimator of (3.18) is indicated by

$$ h_n = \left( \left( \frac{\tilde{I}_4 - \tilde{I}_5}{\tilde{I}_1 - 2\tilde{I}_2 + \tilde{I}_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{1}{d+4}} \right), \quad (3.20) $$

Finally, a non-recursive estimator of $MWISE[r_{\varphi_n}]$ is determined by

$$ \widetilde{MWISE}[\tilde{r}_{\varphi_n}] = \frac{5}{4} (\tilde{I}_4 - \tilde{I}_5)^{\frac{4}{d+4}} (\tilde{I}_1 - 2\tilde{I}_2 + \tilde{I}_3)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{4}{d+4}} + o(n^{-\frac{4}{d+4}}). $$

### 3.5.2 Wild Bootstrap approach

The basic idea of the wild bootstrap introduced in Hardle and Marron (1991) lies in resampling from the estimated residuals

$$ \varepsilon_i = \varphi(Y_i) - r_n(X_i) $$

instead of resampling from the pairs $(Y_i, X_i)_{i=1}^n$ and then investing the obtained data to construct an estimator whose distribution will approximate the distribution of the original estimator. Notice that each bootstrapped residual $\varepsilon_i$ is drawn from a two-point distribution, such that

$$ \mathbb{E}(\varepsilon_i^*) = 0, \quad \mathbb{E}(\varepsilon_i^{*2}) = \hat{\varepsilon}_i^2 \quad \text{and} \quad \mathbb{E}(\varepsilon_i^{*3}) = \hat{\varepsilon}_i^3. $$

Such distribution is expressed by

$$ G_i^* = \left( \frac{5 + \sqrt{5}}{10} \right) \delta_{\hat{\varepsilon}_i \frac{(1-\sqrt{5})}{2}} + \left( \frac{5 - \sqrt{5}}{10} \right) \delta_{\hat{\varepsilon}_i \frac{(1+\sqrt{5})}{2}}. $$

Our adapted procedure for bandwidth selection to estimate the operator $r_{\varphi}$ recursively relies on three steps:

1. Giving the bootstrapped residuals $\varepsilon_i^*$ drawn from the distribution $G_i^*$.

2. Resampling new observations $\varphi(Y_i^*) = r_n(X_i, g) + \varepsilon_i^*$ such that $g$ should be oversmoothed ($g$ needs to be larger than $h$).

3. Computing the kernel regression estimator $r_n^*(X_i, h)$, based on the bootstrapped data $(X_i, Y_i^*)_{i=1}^n$.

The bootstrapped bandwidth $h^*$ is then indicated by:

$$ h^* = \underset{h \in H}{\operatorname{argmin}} \left( \frac{1}{NB} \sum_{i=1}^{NB} (r_n^*(X_i, h) - r_n(X_i, g))^2 \right), \quad (3.21) $$

where $H$ is a fixed set of bandwidths and $NB$ is the number of replications.

In order to ameliorate the performance of the bootstrap procedure over the plug-in method, we set $H = ]h_n - \epsilon, h_n + \epsilon[$, where $h_n$ is the plug-in bandwidth and $\epsilon$ is quite close to zero.
## 3.6 Confidence intervals

Now, let $\phi$ denote the distribution function of the standard normal distribution, and let $t_{\lambda/2}$ be such that $\phi(t_{\lambda/2}) = 1 - \frac{\lambda}{2}$ with $\lambda \in (0, 1)$. We set

$$I_{r_n} = [r_n(x) - t_{\frac{\lambda}{2}}\Lambda, r_n(x) + t_{\frac{\lambda}{2}}\Lambda],$$

with

$$\Lambda = \sqrt{C_f(r_n)[C_\sigma(r_n)\sigma_n^2(x) - C_r(r_n)r_n^2(x)]}, \quad \sigma_n^2(x) = \frac{1}{n}\sum_{i=1}^n (\varphi(Y_i) - r_n(X_i))^2$$

and

<table><thead><tr><th>r<sub>n</sub></th><th>case</th><th>C<sub>f</sub>(r<sub>n</sub>)</th><th>C<sub>σ</sub>(r<sub>n</sub>)</th><th>C<sub>r</sub>(r<sub>n</sub>)</th></tr></thead><tbody><tr><td>r<sub>φn</sub></td><td>α > β</td><td>β<sub>n</sub>R(K)<br>h<sup>d</sup><sub>n</sub>f<sub>n</sub>(x)</td><td>1<br>2 - (β - ad)ξ<sub>β</sub></td><td>2<br>1 - (β - ad - ξ<sup>-1</sup><sub>β</sub>)ξ<sub>α</sub> - ξ<sub>α,β</sub><br>2 - (α - ad)ξ<sub>α</sub> - 1<br>2 - (β - ad)ξ<sub>β</sub></td></tr><tr><td>r<sub>φn</sub></td><td>α < β</td><td>γ<sub>n</sub>R(K)<br>h<sup>d</sup><sub>n</sub>f<sub>n</sub>(x)</td><td>ξ<sub>β,α</sub><br>2 - (β - ad)ξ<sub>β</sub></td><td>2<br>1 - (α - ad - ξ<sup>-1</sup><sub>α</sub>)ξ<sub>β</sub> - 1<br>2 - (α - ad)ξ<sub>α</sub> - ξ<sub>β,α</sub><br>2 - (β - ad)ξ<sub>β</sub></td></tr><tr><td>r<sub>φn</sub></td><td>γ<sub>n</sub> = β<sub>n</sub> = 1/n</td><td>R(K)<br>nh<sup>d</sup><sub>n</sub>f<sub>n</sub>(x)</td><td>1<br>1 + ad</td><td>0</td></tr><tr><td>r̃<sub>φn</sub></td><td>1/n</td><td>R(K)<br>nh<sup>d</sup><sub>n</sub>f̃<sub>n</sub>(x)</td><td>1</td><td>0</td></tr></tbody></table>

In fact, since we have (3.28), the Confidence Intervals for means with unknown standard deviation approach ensure

$$\mathbb{P}\left[-t_{\frac{\lambda}{2}} < \sqrt{\beta_n^{-1} h_n^d} \left(\frac{r_n(x) - \mathbb{E}[r_n(x)]}{\sqrt{\Sigma_{\beta,n}(x)}}\right) < t_{\frac{\lambda}{2}}\right] = 1 - \lambda,$$

with

$$\Sigma_{\beta,n}(x) = \frac{R(K)}{f_n(x)} \left[ \frac{\sigma_n^2(x)}{2 - (\beta - ad)\xi_\beta} - r_n^2(x) \left( \frac{2}{1 - (\beta - ad - \xi_\beta^{-1})\xi_\alpha} - \frac{\xi_{\alpha,\beta}}{2 - (\alpha - ad)\xi_\alpha} - \frac{1}{2 - (\beta - ad)\xi_\beta} \right) \right]$$

is an estimator of (3.16).

Therefore, a confidence interval for the coverage error is given by

$$I_{r_n} = \left[ r_n(x) - t_{\frac{\lambda}{2}} \sqrt{\frac{\Sigma_{\beta,n}(x)}{\beta_n^{-1} h_n^d}} , r_n(x) + t_{\frac{\lambda}{2}} \sqrt{\frac{\Sigma_{\beta,n}(x)}{\beta_n^{-1} h_n^d}} \right].$$

## 3.7 Numerical applications

The main target of this section is to perform a simulation study comparing the performance of our semi-recursive estimator (3.1) to that of Nadaraya-Watson (3.2) from confidence interval point of view. Throughout this section, we consider the regression model defined as

$$\varphi(Y) = r_{\varphi}(X) + \varepsilon,$$

where $X$ follows the multivariate normal distribution $N(0_d, \sigma I_d)$ and $\varepsilon$ follows the normal distribution $N(0, \sigma_\varepsilon)$, with $\sigma$ and $\sigma_\varepsilon$ are two positive constants.
### 3.7.1 Simulation studies

We shall start by specifying our kernel function $\mathbf{K}$ choice which is not carried out at random but according to several criteria. The Gaussian kernel has as an expression

$$K(x) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x^2}{2}\right), \text{ for all } x \in \mathbb{R}.$$

Moreover, we shall consider the case where $\lambda = 0.05$ which yield $t_{\lambda/2} = 1.96$.

When applying our estimator $r_{\varphi_n}$, we must choose two quantities:

* The stepsizes $(\beta_n, \gamma_n) = (\beta_0 n^{-1}, \gamma_0 n^{-1})$, where $\beta_0 = 1$ and $\gamma_0 = 1$.

* The bandwidth $(h_n)$ which is chosen to be equal to (3.19) for plug-in recursive estimator (resp. (3.21) for bootstrapped recursive one).

For this special case, we set

$$\hat{I}_n = \left[ r_{\varphi_n}(x) - 1.96 \sqrt{ \frac{R(\mathbf{K}) \sum_{i=1}^{n} (\varphi(Y_i) - r_{\varphi_n}(X_i))^2}{(1+ad)n^2 h_n^d f_n(x)} } , r_{\varphi_n}(x) + 1.96 \sqrt{ \frac{R(\mathbf{K}) \sum_{i=1}^{n} (\varphi(Y_i) - r_{\varphi_n}(X_i))^2}{(1+ad)n^2 h_n^d f_n(x)} } \right].$$

When applying our estimator $\tilde{r}_{\varphi_n}$, we have to opt for the following quantity:

* The bandwidth $(h_n)$ which is chosen to be equal to (3.20) for plug-in non-recursive estimator (resp. (3.21) for bootstrapped non-recursive one).

For this special case, we set

$$\tilde{I}_n = \left[ \tilde{r}_{\varphi_n}(x) - 1.96 \sqrt{ \frac{R(\mathbf{K}) \sum_{i=1}^{n} (\varphi(Y_i) - \tilde{r}_{\varphi_n}(X_i))^2}{n^2 h_n^d \tilde{f}_n(x)} } , \tilde{r}_{\varphi_n}(x) + 1.96 \sqrt{ \frac{R(\mathbf{K}) \sum_{i=1}^{n} (\varphi(Y_i) - \tilde{r}_{\varphi_n}(X_i))^2}{n^2 h_n^d \tilde{f}_n(x)} } \right].$$

In what follows, we denote by $r_i^*$ the reference regression, by $r_i$ the test regression and by $L_i$ the average length of the test confidence interval, then we compute the following measures:

* Mean squared error: $MSE = \frac{1}{n} \sum_{i=1}^{n} (r_i - r_i^*)^2$.

* The linear correlation: $Cor = Cov(r_i, r_i^*) \sigma(r_i)^{-1} \sigma(r_i^*)^{-1}$.

* Mean amplitude of the confidence interval: $MAIC = \frac{1}{Np} \sum_{i=1}^{Np} L_i$.

Aiming to compare the proposed semi-recursive estimator to the non-recursive Nadaraya-Watson one, we consider four sample sizes: $n = 50, 100, 200$ and $500$, a fixed number of simulations: $N=500$ and three models:

* Model 1: $X$ follows the normal distribution $N(0, 5)$ and $r_{\varphi}(x) = \frac{1}{1 + \exp(-x)}$.

* Model 2: $X$ follows the standard normal distribution $N(0, 1)$ and $r_{\varphi}(x) = \cos(x)$.

* Model 3: $X$ follows the standard bivariate normal distribution $N((0^0), (1^0))$ and $r_{\varphi}(x_1, x_2) = \exp(-x_1^2) + \sin(x_2)$.
<table><caption>Table 3.1: Quantitative comparison between Nadaraya-Watson estimator and the proposed estimator with stepsizes $(\beta_n, \gamma_n) = (n^{-1}, n^{-1})$ through a plug-in method and a bootstrap one in the unidimensional case.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">σ<sub>ε</sub></th><th rowspan="2"></th><th rowspan="2">n</th><th colspan="2">Plug-in</th><th colspan="2">Bootstrap</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td rowspan="12">Model 1</td><td rowspan="12">0.01</td><td rowspan="4">MSE</td><td>50</td><td>0.00316476</td><td>0.00282563</td><td>0.00307640</td><td><b>0.00273586</b></td></tr><tr><td>100</td><td>0.00176918</td><td>0.00165111</td><td>0.00171129</td><td><b>0.00159965</b></td></tr><tr><td>200</td><td>0.00111536</td><td>0.00103373</td><td>0.00108654</td><td><b>0.00099738</b></td></tr><tr><td>500</td><td>0.00090141</td><td>0.00085514</td><td>0.00088022</td><td><b>0.00083644</b></td></tr><tr><td rowspan="4">Cor</td><td>50</td><td>0.99499632</td><td>0.99573655</td><td>0.99512281</td><td><b>0.99585851</b></td></tr><tr><td>100</td><td>0.99690637</td><td>0.99723254</td><td>0.99699979</td><td><b>0.99731362</b></td></tr><tr><td>200</td><td>0.99797460</td><td>0.99819921</td><td>0.99802550</td><td><b>0.99825906</b></td></tr><tr><td>500</td><td>0.99828980</td><td>0.99841832</td><td>0.99832927</td><td><b>0.99845168</b></td></tr><tr><td rowspan="4">MAIC</td><td>50</td><td>0.05291563</td><td>0.05051224</td><td>0.04785009</td><td><b>0.04556449</b></td></tr><tr><td>100</td><td>0.03647338</td><td>0.03533060</td><td>0.02995457</td><td><b>0.02911132</b></td></tr><tr><td>200</td><td>0.02599057</td><td>0.02479223</td><td>0.02181338</td><td><b>0.02103482</b></td></tr><tr><td>500</td><td>0.02089086</td><td>0.02039699</td><td>0.01663979</td><td><b>0.01614454</b></td></tr><tr><td rowspan="12">Model 2</td><td rowspan="12">0.01</td><td rowspan="4">MSE</td><td>50</td><td>0.04793966</td><td>0.04395630</td><td>0.04683956</td><td><b>0.04317905</b></td></tr><tr><td>100</td><td>0.04366443</td><td>0.03475473</td><td>0.04086334</td><td><b>0.03142012</b></td></tr><tr><td>200</td><td>0.01696163</td><td>0.00593385</td><td>0.01698463</td><td><b>0.00592159</b></td></tr><tr><td>500</td><td>0.01596187</td><td>0.00438506</td><td>0.01534741</td><td><b>0.00436068</b></td></tr><tr><td rowspan="4">Cor</td><td>50</td><td>0.98491608</td><td>0.98493530</td><td>0.98295729</td><td><b>0.98600537</b></td></tr><tr><td>100</td><td>0.98738420</td><td>0.98756152</td><td>0.98746333</td><td><b>0.98772430</b></td></tr><tr><td>200</td><td>0.98387705</td><td>0.99434872</td><td>0.98379161</td><td><b>0.99423778</b></td></tr><tr><td>500</td><td>0.99639055</td><td>0.99739599</td><td>0.99640605</td><td><b>0.99739693</b></td></tr><tr><td rowspan="4">MAIC</td><td>50</td><td>0.53178680</td><td>0.50550970</td><td>0.52328180</td><td><b>0.47993540</b></td></tr><tr><td>100</td><td>0.43686240</td><td>0.41610057</td><td>0.42394440</td><td><b>0.37903800</b></td></tr><tr><td>200</td><td>0.45451390</td><td>0.36258430</td><td>0.42906110</td><td><b>0.35669820</b></td></tr><tr><td>500</td><td>0.36321757</td><td>0.32708432</td><td>0.34213073</td><td><b>0.30261850</b></td></tr></tbody></table>

<figure><img src="image_35.png" /><figcaption>Figure 3.1: Qualitative comparison between the Nadaraya-Watson estimator and the recursive estimator for Model 1 with n=50 and $\sigma_{\epsilon} = 0.01$.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Figure 3.2: Qualitative comparison between the Nadaraya-Watson estimator and the recursive estimator for Model 1 with n=500 and $\sigma_{\epsilon} = 0.01$.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Figure 3.3: Qualitative comparison between the Nadaraya-Watson estimator and the recursive estimator for Model 2 with n=100 and $\sigma_{\epsilon} = 0.1$.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 3.4: Qualitative comparison between the Nadaraya-Watson estimator and the recursive estimator for Model 2 with n=200 and $\sigma_{\epsilon} = 0.1$.</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>Figure 3.5: The reference regression function for Model 3 for one simple simulation with $n = 500$.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Figure 3.6: The recursive regression estimator for Model 3 for one simple simulation with $n = 500$.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Figure 3.7: The non-recursive regression estimator for Model 3 for one simple simulation with $n = 500$.</figcaption></figure>
### 3.7.2 Real Datasets

#### Application 1: French Hospital Data of COVID19

The French Hospital data of the COVID-19 epidemic are found in

https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/.

The **Santé publique France**'s mission is devoted to improve and protect the health of population. During the health crisis related to the COVID-19 epidemic, **Santé publique France** has been in charge of monitoring and understanding the dynamics of the epidemic, anticipating the different scenarios and implementing actions so as to prevent and limit the spread of this virus on the national territory.

#### Description of the dataset

This dataset provides information on the hospital situation regarding the COVID-19 epidemic. We have chosen the first proposed file:

Hospital data related to the COVID-19 epidemic by department (dep) and sex (sex) of the patient: number of hospitalized patients (hosp), number of persons currently in intensive care or resuscitation (rea), number of persons currently in follow-up and rehabilitation care (SSR) or long-term care units (USLD), number of persons currently in conventional hospitalization (HospConv), number of persons currently hospitalized in another type of service (autres) or cumulative number of persons having returned home (rad), cumulative number of persons who died (dc).

The data are daily updated. For the current application, we have selected the data of 28/07/2021, with a total of 150894 observations. For simplicity reasons, we opted for focusing just on the department of 'Paris' database.

As a matter of fact, our application rests upon a dataframe of 1494 observations and 6 variables. The following two models are considered :

* Model 1: $X = \text{rea}, Y = \text{hosp}$ and $\varphi: y \mapsto y$.

* Model 2: $X_1 = \text{rea}, X_2 = \text{dc}, Y = \text{hosp}$ and $\varphi: y \mapsto y$.

<table><caption>Table 3.2: Quantitative comparison between Nadaraya-Watson estimator and the proposed one with stepsizes $(\beta_n, \gamma_n) = (n^{-1}, n^{-1})$ through plug-in method and the bootstrap one.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">Plug-in</th><th colspan="2">Bootstrap</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td>Model 1</td><td>3.648179</td><td>3.647376</td><td>3.648197</td><td><strong>3.647362</strong></td></tr><tr><td>Model 2</td><td>2.852123</td><td>2.514630</td><td>2.852662</td><td><strong>2.514177</strong></td></tr></tbody></table>
<figure><img src="image_42.png" /><figcaption>Figure 3.8: Box-plot of the relative error estimation of the four considered estimators for the bivariate COVID-19 application Model 1.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Figure 3.9: Box-plot of the relative error estimation of the four considered estimators for the bivariate COVID-19 application Model 1.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Figure 3.10: Box-plot of the relative error estimation of the four considered estimators for the bivariate COVID-19 application Model 2.</figcaption></figure>
<figure><img src="image_45.png" /><figcaption>Figure 3.11: Box-plot of the relative error estimation of the four considered estimators for the bivariate COVID-19 application Model 2.</figcaption></figure>

## Application 2: Plasmodium falciparum Parasite Load

As far as our application is concerned, we considered a dataset of 176 families belonging to Senegal, living in two villages of Niakhar (Diohine and Toucar), with 505 children aged between 2 and 19 years old. The total number of observations was 6986. We measured Plasmodium falciparum Parasite Load (PL) from thick blood smears obtained by finger-prick during two different seasons and regularly over a three-year observation period (2001–2003). The number of measurements per child ranged from 1 to 15. For more details about the data, we refer the reader to consult Milet et al. (2010). This application relies upon the following variables:

* PL : Parasite Load, as our response variable $Y$.

* malariae : The presence of co-infection with P. malariae, a factor with two levels (infected: 1 or not infected: 0).

* sex : A factor with two levels (a boy: 0 or a girl: 1).

* age : Age of the child in years between 2 and 19.

* season : A factor with two levels (July-October and October-March).

Therefore, for our selection we have a dataframe of 500 observations and 3 variables. The following two models are considered :

* Model 3: $X_1$ = sex, $X_2$ = age, $Y$ = PL and $\phi: y \mapsto \log(y+1)$.

* Model 4: $X_1$ = age, $X_2$ = malariae, $X_3$ = season, $Y$ = PL and $\phi: y \mapsto \log(y+1)$.

<table><caption>Table 3.3: Quantitative comparison between Nadaraya-Watson estimator and the proposed one with stepsizes $(\beta_n, \gamma_n) = (n^{-1}, n^{-1})$ through plug-in method and the bootstrap one.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">Plug-in</th><th colspan="2">Bootstrap</th></tr><tr><th>Nadaraya's estimator</th><th>Recursive estimator</th><th>Nadaraya's estimator</th><th>Recursive estimator</th></tr></thead><tbody><tr><td>Model 3</td><td>0.8057840</td><td><strong>0.8015399</strong></td><td>0.8058742</td><td>0.8019406</td></tr><tr><td>Model 4</td><td>0.8019984</td><td>0.7962954</td><td>0.8025314</td><td><strong>0.7954731</strong></td></tr></tbody></table>
<figure><img src="image_46.png" /><figcaption>Figure 3.12: Box-plot of the relative error estimation of the four considered estimators for the multivariate PL application Model 3.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>Figure 3.13: Box-plot of the relative error estimation of the four considered estimators for the multivariate PL application Model 3.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>Figure 3.14: Box-plot of the relative error estimation of the four considered estimators for the multivariate PL application Model 4.</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>Figure 3.15: Box-plot of the relative error estimation of the four considered estimators for the multivariate PL application Model 4.</figcaption></figure>

## 3.8 Conclusion

This chapter reports an extension of the semi-recursive regression function estimator. Initially, we tackled the asymptotic properties of the proposed estimator in order to demonstrate that our estimator asymptotically follows a normal distribution. The proposed estimator was compared to the non-recursive multivariate Nadaraya Watson regression estimator. Basically, we revealed that using a specific bandwidth selection, the plug-in approach as well as the bootstrap procedure, and particular stepsizes couple $(\gamma_n, \beta_n) = (n^{-1}, n^{-1})$; the proposed estimator (3.1) often provides better results compared to the non-recursive Nadaraya Watson’s one in terms of estimation error. The simulation studies and real datasets illustrate our findings. Even if the bootstrap approach outperforms the plug-in method, it’s not quite accurate to assert that one method is better than the other. They are indistinguishable and it has been widely proven that they behave similarly. We recommend the reader to consult Delaigle and Gijbels (2004) for a detailed comparison of practical bandwidth selection procedures. In conclusion, the use of our recursive estimator, with an appropriate choice of the bandwidth, enables us to get closer to the true regression function rather than non-recursive one.

## 3.9 Proofs

Throughout this section, we will need the following notations:

$$ \mathcal{Z}_n(x) = h_n^{-d} \varphi(Y_n) \mathbf{K} \left( \frac{x - X_n}{h_n} \right) \quad \text{and} \quad W_n(x) = h_n^{-d} \mathbf{K} \left( \frac{x - X_n}{h_n} \right), \text{ for all } x \in \mathbb{R}^d. $$

### Proof of Proposition 3.1.

This proof is mainly based on the same concept as the second chapter, by assuming $a := a_\varphi$. To this extent, we just briefly outline the proof. We have

$$ a_{\varphi_n}(x) - a_{\varphi}(x) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k (\mathcal{Z}_k(x) - a_{\varphi}(x)) + Q_n [a_0(x) - a_{\varphi}(x)]. $$

Hence,

$$ \mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k (\mathbb{E}[\mathcal{Z}_k(x)] - a_{\varphi}(x)) + Q_n [a_0(x) - a_{\varphi}(x)]. $$
**Bias of $a_{\varphi_n}$**: Resting upon the assumptions $(A_1)$ and $(A_6)$ and by applying Taylor's development formula for $a_{\varphi}$, we deduce that

$$
\begin{align*}
\mathbb{E}[Z_k(x)] - a_{\varphi}(x) &= \int_{\mathbb{R}^{d+1}} h_k^{-d} \mathbf{K} \left( \frac{x-y}{h_k} \right) \varphi(t) g(y,t) dy dt - \int_{\mathbb{R}^d} \mathbf{K}(y) a_{\varphi}(x) dy \\
&= \int_{\mathbb{R}^d} h_k^{-d} \mathbf{K} \left( \frac{x-y}{h_k} \right) \mathbb{E}[\varphi(Y) | X=y] f(y) dy - \int_{\mathbb{R}^d} \mathbf{K}(y) a_{\varphi}(x) dy \\
&= \int_{\mathbb{R}^d} \mathbf{K}(z) [a_{\varphi}(x-zh_k) - a_{\varphi}(x)] dz \\
&= \int_{\mathbb{R}^d} \mathbf{K}(z) \left[ \sum_{i=1}^d \frac{\partial a_{\varphi}}{\partial x_i}(x) z_i h_k + \int_0^1 (1-t) \sum_{i,j=1}^d \frac{\partial^2 a_{\varphi}}{\partial x_i \partial x_j}(x-tzh_k) z_i z_j h_k^2 dt \right] dz \\
&= \frac{h_k^2}{2} \sum_{j=1}^d \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x) + h_k^2 \eta_k(x),
\end{align*}
$$

where $\eta_k(x) = \sum_{i,j=1}^{d} \int_{\mathbb{R}^d} \int_0^1 (1-t) [a_{\varphi_{ij}}^{(2)}(x - tzh_k) - a_{\varphi_{ij}}^{(2)}(x)] z_i z_j \mathbf{K}(z) dt dz$.

We thus get

$$
\mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x) = \frac{1}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x) Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k h_k^2 + Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k h_k^2 \eta_k(x) \\
+ Q_n [a_0(x) - a_{\varphi}(x)].
$$

Since $a_{\varphi_{ij}}^{(2)}$ is bounded and continuous at $x$, we deduce that $\lim_{k \to +\infty} \eta_k(x) = 0$.

For the case $a \le \beta/(d+4)$, we have $\lim_{+\infty} (n\beta_n) > 2a$ and then $1-2a\xi_\beta > 0$. Hence, the application of lemma 1.2 provides

$$
\begin{align*}
\mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x) &= \frac{1}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x) Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k h_k^2 + Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k o(h_k^2) + O(Q_n) \\
&= \frac{h_n^2}{2(1-2a\xi_\beta)} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x) + o(h_n^2) + o(1) + O(Q_n).
\end{align*}
$$

Thus the result can be written as

$$
\mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x) = \frac{h_n^2}{2(1 - 2a\xi_\beta)} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x) + o(h_n^2).
$$

For the case $a > \beta/(d+4)$, we have $\lim_{+\infty} (n\beta_n) > \frac{\beta-a}{2}$, which ensures that $h_n^2 = o(\sqrt{\beta_n h_n^{-d}})$.
Therefore, the application of lemma 1.2 entails

$$
\begin{align*}
\mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x) &= \frac{1}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) a_{\varphi_{jj}}^{(2)}(x) Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k o \left( \sqrt{\beta_k h_k^{-d}} \right) \\
&\quad + Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k o \left( \sqrt{\beta_k h_k^{-d}} \right) + O(Q_n) \\
&= o \left( \sqrt{\beta_n h_n^{-d}} \right).
\end{align*}
$$
**Variance of $a_{\varphi_n}$:** For the variance, we infer that

$$ 
\begin{aligned}
\mathrm{Var}[a_{\varphi_n}(x)] &= \mathrm{Var}\left[Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k \mathcal{Z}_k(x)\right] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 (\mathbb{E}[\mathcal{Z}_k^2(x)] - \mathbb{E}[\mathcal{Z}_k(x)]^2).
\end{aligned}
$$

We have

$$ 
\begin{aligned}
\mathbb{E}[\mathcal{Z}_k^2(x)] &= \int_{\mathbb{R}^d} h_k^{-2d} \mathbb{E}[\varphi(Y)^2 | X = y] \mathbf{K}^2 \left(\frac{x-y}{h_k}\right) f(y) dy \\
&= \int_{\mathbb{R}^d} h_k^{-d} \mathbf{K}^2(z) \mathbb{E}[\varphi(Y)^2 | X = x - zh_k] f(x - zh_k) dz \\
&= h_k^{-d} \left[ \mathbb{E}[\varphi(Y)^2 | X = x] f(x) \int_{\mathbb{R}^d} \mathbf{K}^2(z) dz + \nu_k(x) \right],
\end{aligned}
$$

with

$$ \nu_k(x) = \int_{\mathbb{R}^d} \mathbf{K}^2(z) \left[ \mathbb{E}[\varphi(Y)^2 | X = x - zh_k] f(x - zh_k) - \mathbb{E}[\varphi(Y)^2 | X = x] f(x) \right] dz. $$

Thus,

$$ \mathrm{Var}[a_{\varphi_n}(x)] = Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \left[ \mathbb{E}[\varphi(Y)^2 | X = x] f(x) \int_{\mathbb{R}^d} \mathbf{K}^2(z) dz + \nu_k(x) - h_k^d \kappa_k(x) \right], $$

$$ \text{where } \kappa_k(x) = \left( \int_{\mathbb{R}^d} \mathbf{K}(z) a_{\varphi}(x - zh_k) dz \right)^2. $$

Resting upon $(A_6)$, we have that the function $s \mapsto \mathbb{E}[\varphi(Y)^2 | X = s] f(s) = \int_{\mathbb{R}} \varphi(y)^2 g(s, y) dy$ is bounded and continuous at $s = x$ which ensures that $\lim_{k \to +\infty} \nu_k(x) = 0$.

For the case $a \ge \beta/(d+4)$, we have $\lim_{+\infty} (n\beta_n) > \frac{\beta - ad}{2}$ and therefore $2 - (\beta - ad)\xi_\beta > 0$. Since we have $h_k\kappa_k(x) = o(1)$ and $\nu_k(x) = o(1)$, then the application of lemma 1.2 yields

$$ 
\begin{aligned}
\mathrm{Var}[a_{\varphi_n}(x)] &= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \left[ \mathbb{E}[\varphi(Y)^2 | X = x] f(x) R(\mathbf{K}) + \nu_k(x) - h_k^d \kappa_k(x) \right] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \left[ \mathbb{E}[\varphi(Y)^2 | X = x] f(x) R(\mathbf{K}) + o(1) \right] \\
&= \frac{\mathbb{E}[\varphi(Y)^2 | X = x] \beta_n}{2 - (\alpha - ad)\xi_\beta} \frac{h_n}{h_n} [f(x) R(\mathbf{K}) + o(1)].
\end{aligned}
$$

Thus, the result is indicated in terms of

$$ \mathrm{Var}[a_{\varphi_n}(x)] = \frac{\mathbb{E}[\varphi(Y)^2 | X = x] \beta_n}{2 - (\alpha - a)\xi_\beta} \frac{h_n}{h_n} f(x) R(\mathbf{K}) + o\left(\frac{\beta_n}{h_n}\right). $$

For the case $a < \beta/(d+4)$, we have $\lim_{+\infty} (n\beta_n) > 2a$, which ensures that $\beta_n h_n^{-d} = o(h_n^4)$. By applying lemma 1.2, we obtain

$$ 
\begin{aligned}
\mathrm{Var}[a_{\varphi_n}(x)] &= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \left[ \mathbb{E}[\varphi(Y)^2 | X = x] f(x) R(\mathbf{K}) + o(1) \right] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k o(h_k^4) \\
&= o(h_n^4).
\end{aligned}
$$
### Proof of Theorem 3.2.

This proof is based on the following observation

$$r_{\varphi_n}(x) - r_{\varphi}(x) = D_n(x) \frac{f(x)}{f_n(x)}, \quad f_n \neq 0 \qquad (3.22)$$

with

$$D_n(x) = \frac{1}{f(x)} (a_{\varphi_n}(x) - a_{\varphi}(x)) - \frac{r_{\varphi}(x)}{f(x)} (f_n(x) - f(x)).$$

The only remaining point concerns the asymptotic behaviour of $r_{\varphi_n}(x) - r_{\varphi}(x)$, which can be deduced from that of $D_n(x)$. Hence, we can state

$$\mathbb{E}[D_n(x)] = \frac{1}{f(x)} (\mathbb{E}[a_{\varphi_n}(x)] - a_{\varphi}(x)) - \frac{r_{\varphi}(x)}{f(x)} (\mathbb{E}[f_n(x)] - f(x)).$$

Combining the bias of $a_{\varphi_n}(x)$ (3.3) and (3.4)) as well as that of $f_n(x)$ (2.7) and (2.8)) yields the desired results (3.7) and (3.8).

For the variance, we get

$$\mathrm{Var}[D_n(x)] = \frac{1}{(f(x))^2} \mathrm{Var}[a_{\varphi_n}(x)] - \frac{(r_{\varphi}(x))^2}{(f(x))^2} \mathrm{Var}[f_n(x)] - 2 \frac{r_{\varphi}(x)}{(f(x))^2} \mathrm{Cov}(a_{\varphi_n}(x), f_n(x)).$$

### 1. For the case $\beta \le \alpha$:

Since $X_k$'s are independent, then for all $i \ne k$, $\mathrm{Cov}(Z_k(x), W_i(x)) = 0$ and by applying lemma 1.2, classical computations entail

$$\mathrm{Cov}(a_{\varphi_n}(x), f_n(x)) = r_{\varphi}(x)f(x)R(\mathbf{K})\beta_n h_n^{-d} \left( \frac{1}{1 - (\beta - ad - \xi_{\beta}^{-1})\xi_{\alpha}} + o(1) \right). \quad (3.23)$$

In fact, we have

$$
\begin{align*}
\mathrm{Cov}(a_{\varphi_n}(x), f_n(x)) &= \mathrm{Cov} \left( Q_n \sum_{k=1}^n Q_k^{-1} \beta_k \varphi(Y_k) h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right), \Pi_n \sum_{i=1}^n \Pi_i^{-1} \gamma_i h_i^{-d} \mathbf{K} \left( \frac{x-X_i}{h_i} \right) \right) \\
&= Q_n \sum_{k=1}^n Q_k^{-1} \beta_k \Pi_n \sum_{i=1}^n \Pi_i^{-1} \gamma_i \mathrm{Cov} \left( \varphi(Y_k) h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right), h_i^{-d} \mathbf{K} \left( \frac{x-X_i}{h_i} \right) \right) \\
&= Q_n \Pi_n \sum_{k=1}^n \Pi_k^{-1} Q_k^{-1} \gamma_k \beta_k \mathrm{Cov} \left( \varphi(Y_k) h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right), h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right) \right) \\
&= Q_n \Pi_n \sum_{k=1}^n \Pi_k^{-1} Q_k^{-1} \gamma_k \beta_k \left( \mathbb{E} \left[ \varphi(Y_k) h_k^{-2d} \mathbf{K}^2 \left( \frac{x-X_k}{h_k} \right) \right] \right. \\
&\qquad \left. - \mathbb{E} \left[ \varphi(Y_k) h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right) \right] \mathbb{E} \left[ h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right) \right] \right) \\
&= Q_n \Pi_n \sum_{k=1}^n \Pi_k^{-1} Q_k^{-1} \gamma_k \beta_k \left( \mathbf{E}[\varphi(Y)|X=x] f(x) R(\mathbf{K}) h_k^{-d} + o(h_k^{-d}) \right. \\
&\qquad \left. - \mathbf{E}[\varphi(Y)|X=x] f^2(x) + o(1) \right) \\
&= Q_n \Pi_n \sum_{k=1}^n \Pi_k^{-1} Q_k^{-1} \gamma_k \beta_k h_k^{-d} \left( r_{\varphi}(x) f(x) R(\mathbf{K}) + o(1) \right) \\
&= \frac{\beta_n h_n^{-d}}{1 - (\beta - ad - \xi_{\beta}^{-1})\xi_{\alpha}} r_{\varphi}(x) f(x) R(\mathbf{K}) + o(\beta_n h_n^{-d}).
\end{align*}
$$
Consequently, (3.9) and (3.10) follow from the combination of the variance of $a_{\varphi_n}(x)$ ((3.5) and (3.6)), as well as from that of $f_n(x)$ ((2.9) and (2.10)) and the covariance expression (3.23).

It is noteworthy that, for the case $a \ge \beta/(d+4)$, we deduce

$$
\begin{align*}
\mathrm{Var}[r_{\varphi_n}(x)] &= \frac{1}{f(x)} \frac{\beta_n \mathbb{E}[\varphi(Y)^2 | X=x]}{h_n^d (2-(\beta-ad)\xi_\beta)} R(\mathbf{K}) + o\left(\beta_n h_n^{-d}\right) + \frac{r_\varphi(x)^2 \gamma_n}{f(x)} \frac{1}{h_n^d (2-(\alpha-ad)\xi_\alpha)} R(\mathbf{K}) \\
&\quad + o\left(\gamma_n h_n^{-d}\right) - 2 \frac{r_\varphi(x)}{f(x)^2} \frac{\beta_n h_n^{-d}}{1-(\beta-ad-\xi_\beta^{-1})\xi_\alpha} r_\varphi(x) f(x) R(\mathbf{K}) + o\left(\beta_n h_n^{-d}\right) \\
&= \frac{1}{f(x)} \frac{\beta_n \mathbb{E}[\varphi(Y)^2 | X=x]}{h_n^d (2-(\beta-ad)\xi_\beta)} R(\mathbf{K}) + o\left(\beta_n h_n^{-d}\right) + \frac{r_\varphi(x)^2 \beta_n}{f(x)} \frac{\xi_{\alpha,\beta}}{h_n^d (2-(\alpha-ad)\xi_\alpha)} R(\mathbf{K}) \\
&\quad - 2 \frac{r_\varphi(x)^2}{f(x)} \frac{\beta_n h_n^{-d}}{1-(\beta-ad-\xi_\beta^{-1})\xi_\alpha} R(\mathbf{K}) + o\left(\beta_n h_n^{-d}\right) \\
&= \frac{\beta_n R(\mathbf{K})}{h_n^d f(x)} \left[ \frac{\mathbb{E}[\varphi(Y)^2 | X=x]}{2-(\beta-ad)\xi_\beta} - r_\varphi(x)^2 \left( \frac{2}{1-(\beta-ad-\xi_\beta^{-1})\xi_\alpha} - \frac{\xi_{\alpha,\beta}}{2-(\alpha-ad)\xi_\alpha} \right) \right] \\
&\quad + o\left(\beta_n h_n^{-d}\right).
\end{align*}
$$

2. For the case $\alpha < \beta$:

Similarly to the first case, and taking the stepsize $(\gamma_n)$ as a reference, we infer the result.

□

Proof of Theorem 3.4.

We have

$$
\begin{align*}
D_n(x) - \mathbb{E}[D_n(x)] &= \frac{1}{f(x)}[a_{\varphi_n}(x) - \mathbb{E}[a_{\varphi_n}(x)]] - \frac{r_{\varphi}(x)}{f(x)}[f_n(x) - \mathbb{E}[f_n(x)]] \\
&= \frac{1}{f(x)}Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k (T_k(x) - \mathbb{E}[T_k(x)]),
\end{align*}
$$

with

$$
T_k(x) = \mathcal{Z}_k(x) - r_\varphi(x) Q_n^{-1} \Pi_n \Pi_k^{-1} Q_k \beta_k^{-1} \gamma_k \mathcal{W}_k(x).
$$

We note

$$
S_k(x) = Q_k^{-1} \beta_k (T_k(x) - \mathbb{E}[T_k(x)]). \qquad (3.24)
$$

Hence, we can write

$$
D_n(x) - \mathbb{E}[D_n(x)] = \frac{1}{f(x)} Q_n \sum_{k=1}^{n} S_k(x). \qquad (3.25)
$$

Now, we are trying to apply Lyapunov's theorem 1.14 for $S_k(x)$. For this reason, we assume

$$
\begin{align*}
\upsilon_n^2 &= \sum_{k=1}^{n} \mathrm{Var}[S_k(x)] \\
&= \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 \mathrm{Var}[T_k(x)] \\
&= \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 \mathrm{Var}[\mathcal{Z}_k(x)] + r_\varphi(x)^2 Q_n^{-2} \Pi_n^2 \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 \mathrm{Var}[\mathcal{W}_k(x)] \\
&\quad - 2 r_\varphi(x) Q_n^{-1} \Pi_n \sum_{k=1}^{n} Q_k^{-1} \beta_k \Pi_k^{-1} \gamma_k \mathrm{Cov}(\mathcal{Z}_k(x), \mathcal{W}_k(x)).
\end{align*}
$$
Here, we consider the case $\beta \le \alpha$. Since we have

$$ \mathrm{Var}[\mathcal{Z}_k(x)] = h_k^{-d} (\mathbb{E}[\varphi(Y)^2 | X = x] f(x) R(\mathbf{K}) + o(1)), $$

$$ \mathrm{Var}[\mathcal{W}_k(x)] = h_k^{-d} (f(x) R(\mathbf{K}) + o(1)), $$

$$ \mathrm{Cov}(\mathcal{Z}_k(x), \mathcal{W}_k(x)) = h_k^{-d} (r_{\varphi}(x) f(x) R(\mathbf{K}) + o(1)), $$

then, the application of lemma 1.2 ensures that

$$ \begin{aligned} v_n^2 &= \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \left( \mathbb{E}[\varphi(Y)^2 | X = x] f(x) R(\mathbf{K}) + o(1) \right) \\ &\quad + r_{\varphi}(x)^2 Q_n^{-2} \Pi_n^2 \sum_{k=1}^{n} \Pi_k^{-2} \gamma_k^2 h_k^{-d} \left( f(x) R(\mathbf{K}) + o(1) \right) \\ &\quad - 2r_{\varphi}(x) Q_n^{-1} \Pi_n \sum_{k=1}^{n} Q_k^{-1} \beta_k \Pi_k^{-1} \gamma_k h_k^{-d} \left( r_{\varphi}(x) f(x) R(\mathbf{K}) + o(1) \right) \\ &= \frac{\beta_n}{h_n^d} \frac{f(x)^2}{Q_n^2} [\Sigma_{\beta}(x) + o(1)]. \end{aligned} $$

On the other side, we have

$$ \forall p > 0, \quad \mathbb{E}[|T_k(x)|^{2+p}] = O\left(\frac{1}{h_k^{d(1+p)}}\right). $$

Therefore,

$$ \begin{aligned} \mathbb{E}[|S_k(x)|^{2+p}] &= Q_k^{-2-p} \beta_k^{2+p} \mathbb{E}[|T_k(x) - \mathbb{E}[T_k(x)]|^{2+p}] \\ &\le 2Q_k^{-2-p} \beta_k^{2+p} \mathbb{E}[|T_k(x)|^{2+p}]. \end{aligned} $$

Hence,

$$ \mathbb{E}[|S_k(x)|^{2+p}] = O\left(Q_k^{-2-p} \beta_k^{2+p} \frac{1}{h_k^{d(1+p)}}\right). \quad (3.26) $$

As a consequence,

$$ \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = O\left(\sum_{k=1}^{n} Q_k^{-2-p} \beta_k^{2+p} \frac{1}{h_k^{d(1+p)}}\right). $$

In what follows, let us suppose that $a \ge \beta/(d+4)$ and assume that there is $p > 0$, such that

$$ \lim_{n \to +\infty} (n\beta_n) > \frac{1+p}{2+p}(\beta - ad). $$

The application of lemma 1.2 yields

$$ \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = O\left(\frac{\beta_n^{1+p}}{Q_n^{2+p} h_k^{d(1+p)}}\right). $$

Hence,

$$ \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = O\left(\frac{\beta_n^{1+p}}{v_n^{2+p} Q_n^{2+p} h_n^{d(1+p)}}\right). $$
Thus, we deduce

$$ \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = O\left(\left(\frac{\beta_n}{h_n^d}\right)^{p/2}\right) = o(1). $$

In addition, since we have

$$ \lim_{n \to +\infty} \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}\left[|S_k(x) - \mathbb{E}[S_k(x)]|^{2+p}\right] = \lim_{n \to +\infty} \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = 0, $$

therefore, by applying the Lyapunov theorem, we get

$$ \frac{1}{\sqrt{v_n^2}} \sum_{k=1}^{n} (S_k(x) - \mathbb{E}[S_k(x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1), $$

which implies

$$ \frac{1}{v_n} \sum_{k=1}^{n} S_k(x) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1). $$

Moreover, (3.22) and (3.25) ensure that

$$ f(x)Q_n^{-1}v_n^{-1}(r_{\varphi_n}(x) - \mathbb{E}[r_{\varphi_n}(x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1). \quad (3.27) $$

Given that $v_n^2 = \frac{\beta_n f(x)^2}{h_n^d Q_n^2}$ [$\Sigma_\beta(x) + o(1)$], where $\Sigma_\beta(x)$ is defined in (3.12) and by replacing $v_n$ with its value in (3.27), we conclude that

$$ \sqrt{\beta_n^{-1} h_n^d} (r_{\varphi_n}(x) - \mathbb{E}[r_{\varphi_n}(x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, \Sigma_\beta(x)). \quad (3.28) $$

The convergence in (3.11) then follows from the application of Lyapounov's Theorem and the combination between (3.7), (3.8) and (3.28).

The case $a < \beta/(d+4)$ is fulfilled in the convergence in probability. By applying the Bien-aymé-Chebyshev inequality, we get

$$ \mathbb{P}\left[\left|\frac{r_{\varphi_n}(x) - r_{\varphi}(x)}{h_n^2}\right| - \mathbb{E}\left[\left|\frac{r_{\varphi_n}(x) - r_{\varphi}(x)}{h_n^2}\right|\right| \ge \epsilon\right] \le \frac{\mathrm{Var}[r_{\varphi_n}(x)]}{h_n^4 \epsilon^2}. $$

Since we have $\beta_n^{-1}h_n^{d+4} \xrightarrow{n \to +\infty} +\infty$, then we deduce that

$$ \frac{1}{h_n^2} (r_{\varphi_n}(x) - r_{\varphi}(x)) \xrightarrow{n \to +\infty} \mathbf{M}_{\beta}(x), $$

with $\mathbf{M}_{\beta}(x)$ is provided in (3.12). $\square$

For the next proof of the strong convergence rate, we first need to introduce a real LIL version of the Mokaddem and Pelletier (2007b) theorem 1.

Let $(X_n)_{n \ge 1}$ be a sequence of independent random real variables with $\mathbb{E}[X_n] = 0$. Set

$$ \zeta_n = \sum_{k=1}^{n} X_k \quad \text{and} \quad B_n = \sum_{k=1}^{n} \mathbb{E}[X_k^2]. $$

We suppose that we have the following assumptions
(AS₁) (i) $\mathbb{E}[|X_n|^3] < +\infty$.

(ii) $\frac{1}{n\sqrt{(n)}} \sum_{k=1}^{n} \mathbb{E}[|B_n^{-1/2}X_k|^3] \le \frac{1}{\ln(B_n)}$.

(AS₂) (i) There exists a positive quantity $\Gamma$ such that $\lim_{n \to +\infty} H_n^2 B_n = \Gamma$.

(ii) $\lim_{n \to +\infty} H_n^{-2} = +\infty$ and $\lim_{n \to +\infty} \frac{H_n^2}{H_{n-1}^2} = 1$.

Note that the assumption (AS₁)(ii) is the condition which (together with the Lyapunov condition) ensures that $\zeta_n$ satisfies the central limit theorem

$$H_n \zeta_n \xrightarrow{n \to +\infty} \mathcal{N}(0, \Gamma).$$

**Theorem 3.11.** *Laws of the iterated logarithm (LIL) for $\zeta_n$, under assumptions (AS₁) and (AS₂), ensure that, with probability one, the sequence*

$$\left( \frac{H_n \zeta_n(x)}{\sqrt{2 \ln \ln(H_n^{-2})}} \right)$$

is relatively compact and its limit set is the interval

$$[-\Gamma, \Gamma].$$

*Proof of Theorem 3.5.*

For this proof, we state

$$\zeta_n(x) := \sum_{k=1}^{n} S_k(x) = f(x)Q_n^{-1}(D_n(x) - \mathbb{E}[D_n(x)]),$$

where $S_k$ is given in (3.24).

Here, we consider the case $\beta \le \alpha$. We suppose that $a \ge \frac{\beta}{d+4}$ and set $\beta_0 = h_0 = 1$ and $H_n^2 = Q_n^2 \beta_n^{-1} h_n^d$, then we get

$$
\begin{align*}
\ln(H_n^{-2}) &= \ln(Q_n^{-2}) + \ln(\beta_n h_n^{-d}) \\
&= -2 \ln(Q_n) + \ln\left(\prod_{k=1}^{n} \frac{\beta_{k-1}^{-1} h_{k-1}^d}{\beta_k^{-1} h_k^d}\right) \\
&= -2 \sum_{k=1}^{n} \ln(1 - \beta_k) + \sum_{k=1}^{n} \ln\left(1 - \frac{\beta - ad}{k} + o\left(\frac{1}{k}\right)\right) \\
&= -2 \sum_{k=1}^{n} (-\beta_k + o(\beta_k)) + \sum_{k=1}^{n} (-\beta - ad) \beta_k \xi + o(\beta_k) \\
&= \sum_{k=1}^{n} (2\beta_k - (\beta - ad) \beta_k \xi + o(\beta_k)).
\end{align*}
$$

Hence, using the notation $s_n = \sum_{k=1}^{n} \beta_k$, we can write

$$\ln(H_n^{-2}) = (2 - (\beta - ad)\xi) s_n + o(s_n). \quad (3.29)$$
Since $2-(\beta-ad)\xi > 0$ and $s_\infty$ diverges, then we deduce that $\lim_{n \to +\infty} H_n^{-2} = +\infty$ and $\lim_{n \to +\infty} \frac{H_{n-1}^{-2}}{H_n^{-2}} = 1$.

Moreover, we have $\sum_{k=1}^{n} \mathrm{Var}[S_k(x)] = \frac{\beta_n}{h_n^d} \frac{f(x)^2}{Q_n^2} [\Sigma_\beta(x) + o(1)]$, where $\Sigma_\beta(x)$ is defined in (3.16). From this perspective, it's obvious that

$$ \lim_{n \to +\infty} H_n^2 \sum_{k=1}^{n} \mathrm{Var}[S_k(x)] = f(x)^2 \Sigma_\beta(x). $$

Considering the particular case of $p=1$ in (3.26), we have $\mathbb{E}[|S_k(x)|^3] = O(Q_k^{-3}\beta_k^3 h_k^{-2d})$ and then we deduce that

$$ \begin{aligned} \frac{1}{n\sqrt{n}} \sum_{k=1}^{n} \mathbb{E}[|H_n S_k(x)|^3] &= O\left(\frac{H_n^3}{n\sqrt{n}} \sum_{k=1}^{n} Q_k^{-3} \beta_k^3 h_k^{-2d}\right) \\ &= O\left(\frac{H_n^3}{n\sqrt{n}} \sum_{k=1}^{n} Q_k^{-3} \beta_k o\left([\beta_k h_k^{-d}]^{\frac{3}{2}}\right)\right) \\ &= O\left(\frac{H_n^3}{n\sqrt{n}} Q_n^{-3} o\left([\beta_n h_n^{-d}]^{\frac{3}{2}}\right)\right) \\ &= o\left(\frac{H_n^3}{n\sqrt{n}} Q_n^{-3} [\beta_n h_n^{-d}]^{\frac{3}{2}}\right) \\ &= o\left(\frac{1}{n\sqrt{n}}\right) \\ &= o([\ln(H_n^{-2})]^{-1}). \end{aligned} $$

The application of Theorem (3.11) then ensures that, with probability one, the sequence

$$ \left( \frac{H_n \zeta_n(x)}{\sqrt{2 \ln \ln(H_n^{-2})}} \right) = \left( \frac{\sqrt{\beta_n^{-1} h_n^d} f(x) (D_n(x) - \mathbb{E}[D_n(x)])}{\sqrt{2 \ln \ln(H_n^{-2})}} \right) $$

is relatively compact and its limit set is the interval

$$ \left[-f(x)\sqrt{\Sigma_{\beta}(x)}, f(x)\sqrt{\Sigma_{\beta}(x)}\right]. $$

On account of (3.29), we have $\lim_{n \to +\infty} \frac{\ln \ln(H_n^{-2})}{\ln(s_n)} = 1$, and referring to (3.22) and (3.25), we deduce that

$$ \left( \frac{\sqrt{\beta_n^{-1} h_n^d} (r_{\varphi_n}(x) - \mathbb{E}[r_{\varphi_n}(x)])}{\sqrt{2 \ln s_n}} \right) $$

is relatively compact and its limit set is the interval

$$ \left[-\sqrt{\Sigma_{\beta}(x)}, \sqrt{\Sigma_{\beta}(x)}\right]. $$

The combination between (3.7) and (3.8) then entails

$$ \left( \sqrt{\frac{\beta_n^{-1} h_n^d}{2 \ln(s_n)}} (r_{\varphi_n}(x) - r_{\varphi}(x)) \right) $$
is relatively compact and its limit set is the interval

$$ \left[ \sqrt{\frac{b}{2}} \mathbf{M}_{\beta}(x) - \sqrt{\mathbf{\Sigma}_{\beta}(x)}, \sqrt{\frac{b}{2}} \mathbf{M}_{\beta}(x) + \sqrt{\mathbf{\Sigma}_{\beta}(x)} \right], $$

where $\mathbf{M}_{\beta}(x)$ is defined in (3.12) and $\mathbf{\Sigma}_{\beta}(x)$ is provided in (3.16).

Now we suppose that $a < \frac{\beta}{d+4}$. Set $H_n^{-2} = Q_n^{-2}h_n^4(\ln \ln(Q_n^{-2}h_n^4)^{-1})$, then we get

$$ \begin{aligned} \ln(Q_n^{-2}h_n^4) &= \ln(Q_n^{-2}) + \ln(h_n^4) \\ &= -2 \ln(Q_n) + \ln\left(\prod_{k=1}^{n} \frac{h_{k-1}^{-4}}{h_k^{-4}}\right) \\ &= -2 \sum_{k=1}^{n} \ln(1 - \beta_k) + \sum_{k=1}^{n} \ln\left(1 - \frac{4a}{k} + o\left(\frac{1}{k}\right)\right) \\ &= -2 \sum_{k=1}^{n} (-\beta_k + o(\beta_k)) + \sum_{k=1}^{n} (-4a\beta_k\xi + o(\beta_k)). \end{aligned} $$

Hence, using the notation $s_n = \sum_{k=1}^{n} \beta_k$, we can write

$$ \ln(H_n^{-2}) = (2 - 4a\xi) s_n + o(s_n). \quad (3.30) $$

Since $2 - 4a\xi > 0$ and $s_{\infty}$ diverges, then we deduce that $\lim_{n \to +\infty} H_n^{-2} = +\infty$ and $\lim_{n \to +\infty} \frac{H_{n-1}^{-2}}{H_n^{-2}} = 1$. Moreover, we have

$$ \begin{aligned} H_n^2 \sum_{k=1}^{n} \mathrm{Var}[S_k(x)] &= O \left( Q_n^2 h_n^{-4} \ln \ln(Q_n^{-2}h_n^4) \frac{\beta_n}{h_n^d Q_n^2} \right) \\ &= o(1). \end{aligned} $$

Considering the particular case of $p = 1$ in (3.26), we have $\mathbb{E}[|S_k(x)|^3] = O(Q_k^{-3}\beta_k^3 h_k^{-2d})$ and then we deduce that

$$ \begin{aligned} \frac{1}{n\sqrt{n}} \sum_{k=1}^{n} \mathbb{E}[|H_n S_k(x)|^3] &= O\left(\frac{H_n^3}{n\sqrt{n}} \sum_{k=1}^{n} Q_k^{-3} \beta_k^3 h_k^{-2d}\right) \\ &= O\left(\frac{H_n^3}{n\sqrt{n}} \sum_{k=1}^{n} Q_k^{-3} \beta_k o(h_k^6)\right) \\ &= o\left(\frac{H_n^3}{n\sqrt{n}} Q_n^{-3} h_n^6\right) \\ &= o([\ln(H_n^{-2})]^{-1}). \end{aligned} $$

The application of Theorem (3.11) then ensures that, with probability one, the sequence

$$ \left( \frac{H_n \zeta_n(x)}{\sqrt{2 \ln \ln(H_n^{-2})}} \right) = \left( \frac{h_n^{-2} \sqrt{\ln \ln(Q_n^{-2}h_n^4)} f(x) (D_n(x) - \mathbb{E}[D_n(x)])}{\sqrt{2 \ln \ln(H_n^{-2})}} \right) $$

is relatively compact and its limit set is 0.

On account of (3.29), we have $\lim_{n \to +\infty} \frac{\ln \ln(H_n^{-2})}{\ln \ln(Q_n^{-2}h_n^4)} = 1$, and referring to (3.22) and (3.25), we deduce that

$$ \lim_{n \to +\infty} h_n^{-2} (r_{\varphi_n}(x) - \mathbb{E}[r_{\varphi_n}(x)]) = 0. $$
The combination between (3.7) and (3.8) then entails

$$ \lim_{n \to +\infty} h_n^{-2} (r_{\varphi_n}(x) - r_{\varphi}(x)) = \mathbf{M}_{\beta}(x), $$

where $\mathbf{M}_{\beta}(x)$ is defined in (3.12). $\square$
# Chapter 4

# Non-parametric multivariate kernel regression estimation to describe cognitive processes and mental representations

## 4.1 Introduction

Research on the handwritten word production aims to describe the cognitive processes and mental representations mobilized when a human being prepares to hand-write a word from an idea (see Perret and Olive(2019)). The most frequently used method to explore this issue relies on relating a behavioral variable, reaction time and a set of factors aiming at predicting different cognitive treatments (e.g., Bonin et al. (2002), Perret and Laganaro (2013)). It is possible to imagine some variations in the cognitive treatments performed by participants. This could result in variations in the relationship between the behavioral variables and the explanatory factors. The intrinsic target lies in being able to group participants with similar degrees of variation.

In order to achieve our purpose, we resort to regression analysis, which corresponds to the study of how a response variable depends on one or more predictors. In fact, it is a reliable method for identifying which variables have impact on a topic of interest. The process of performing a regression allows us to confidently determine which factors matter most, which factors can be ignored, and how these factors influence each other. Regression problems can be usefully summarized using non-parametric regression methods which represent a category of regression analysis in which the predictor does not take a predetermined form but is constructed according to information derived from the data. Since we ignore the behavior of our data, and we don't have the normality (see McGill(1963), McCormack and Wright (1964) and Luce(1986)), we resorted to non-parametric approach. In this task, we shall focus on Kernel regression which is a non-parametric technique in statistics to estimate the conditional expectation of a random variable. The main objective is to find a non-linear relation between a pair of random variables X and T. In addition to the non-parametric fact, we rely on the recursive approach of estimation using stochastic approximation method. The missing data question is a former problem in psychology, which can contaminate the results and disrupt them. In order to settle the missing data problem, multiple 'naive' methods have been incorporated to solve this problem, such as the replacement of the missing value by the mean/median or complete outliers detection and treatment (see Cousineau and Chartier (2010)). Afterword, Slaoui (2017) used the propensity score probability technique and constructed an estimator of the density function under missing data. Our central focus resides in building up a multivariate kernel regression estimator under missing data.
### 4.1.1 Presentation of the method

Let $(X, T)$ be a random vector with values in $\mathbb{R}^d \times \mathbb{R}$ with a joint density function $h(x, t)$ and let $f$ denote the probability density of $X$. Moreover, let $(X_1, T_1), \dots, (X_n, T_n)$ be independent random vectors identically distributed as $(X, T)$. Assuming that $T_1, \dots, T_n$ are subjects to missing data, the observed random variables are then $Y_i$ and $\delta_i$, where

$$ \delta_i = \mathbb{1}_{\{T_i \text{ is observed}\}} \text{ and } Y_i = T_i \times \delta_i, \text{ for all } i \in \{1, \dots, n\}. $$

Accordingly, when some $T_i$ are missing, we introduce the propensity score, a probability elabo- rated by Rosenbaum and Rubin (1983) and defined as follows

$$ \psi_i := \mathbb{P}[\delta_i = 1 | T_i], \text{ for all } i \in \{1, \dots, n\} \text{ and } \psi = \lim_{n \to +\infty} \psi_n. $$

In the remainder, $Y$ is considered as the response variable of interest and $X$ its associated regressor vector variable.

Our basic purpose in this chapter is to propose a recursive estimator to estimate recursively the regression function $p(x) = \mathbb{E}[T|X=x]$ under censoring data.

Our aim then resides in building up a stochastic algorithm, which approaches the regression function

$$ m : x \mapsto \mathbb{E}[T|X=x]f(x) = \int_{\mathbb{R}} t h(x,t)dt $$

at a given vector $x$. For this reason, we define an algorithm of search of the zero function $\phi : y \mapsto m(x) - y$. We therefore proceed as follows, we fix $m_0(x) \in \mathbb{R}$, and then we set for all $n \ge 1$, $m_n(x) = m_{n-1}(x) + \beta_n U_n(x)$, where $U_n(x)$ is an observation of the function $\phi$ at the point $m_{n-1}(x)$. By choosing

$$ U_n(x) = Y_n \psi_n^{-1} h_n^{-d} \mathbf{K} \left( \frac{x - X_n}{h_n} \right) - m_{n-1}(x), $$

the stochastic approximation algorithm that we consider to estimate recursively the regression function $m$ at a vector $x$ can be expressed by :

$$ m_n(x) = (1 - \beta_n)m_{n-1}(x) + \beta_n Y_n \psi_n^{-1} h_n^{-d} \mathbf{K} \left( \frac{x - X_n}{h_n} \right). \quad (4.1) $$

Throughout this section, we consider that $m_0(x) = 0$. It follows that

$$ m_n(x) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k Y_k \psi_k^{-1} h_k^{-d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right). \quad (4.2) $$

In this chapter, we consider the following recursive estimator of the regression function $p$ at the vector $x$

$$ p_n(x) = \begin{cases} \frac{m_n(x)}{f_n(x)} & \text{if } f_n(x) \neq 0, \\ 0 & \text{otherwise} \end{cases}, \quad (4.3) $$

with $f_n$ stands for the recursive density estimator given in (2.3).

We explore the asymptotic properties of our proposed multivariate recursive kernel regression estimator. Afterwards, by introducing the non-recursive estimator of $m$ given by

$$ \tilde{m}_n(x) = \frac{1}{nh_n^d} \sum_{k=1}^{n} Y_k \psi_k^{-1} \mathbf{K} \left( \frac{x - X_k}{h_n} \right), $$
we compare our proposed estimator to the multivariate non-recursive generalized Nadaraya-Watson's regression estimator indicated by

$$ \tilde{p}_n(x) = \begin{cases} \frac{\tilde{m}_n(x)}{\tilde{f}_n(x)} & \text{if } \tilde{f}_n(x) \neq 0 \\ 0 & \text{otherwise} \end{cases}, \quad (4.4) $$

with $\tilde{f}_n$ stands for the non-recursive density estimator given in (2.5).

### 4.1.2 Notations and assumptions

For this section and under (A1) given in 1.1.1 and (A5) provided in 3.1.2, the assumptions upon which we shall rely in this chapter are the following.

**Assumptions:**

(A7) (i) The functions $f$ and $m$ are bounded and twice differentiable.

(ii) For all $i, j \in \{1, \dots, d\}$, $f_{ij}^{(2)} := \frac{\partial^2 f}{\partial x_i \partial x_j}$ and $m_{ij}^{(2)} := \frac{\partial^2 m}{\partial x_i \partial x_j}$ are bounded and continuous at $x$.

(iii) For all $p > 0$, $s \mapsto \int_{\mathbb{R}} |t|^{2+p} h(s, t) dt$ is a bounded function.

(iv) The function $s \mapsto \int_{\mathbb{R}} t^2 h(s, t) dt$ is bounded and continuous at $s = x$.

## 4.2 Main results

In order to investigate the asymptotic properties of our estimator $p_n$, we first need to introduce the following two propositions which yield the bias and the variance of $m_n$.

### 4.2.1 Bias and variance of $m_n$

**Proposition 4.1.** Let assumptions (A1), (A5) and (A7) hold. Hence, we obtain

$$ \text{1. If } a \in (0, \frac{\beta}{d+4}], \text{ then} \qquad \mathbb{E}[m_n(x)] - m(x) = \frac{h_n^2}{2(1 - 2a\xi_\beta)} \sum_{j=1}^{d} \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) + o(h_n^2). \quad (4.5) $$

$$ \text{If } a \in (\frac{\beta}{d+4}, 1), \text{ then} \qquad \mathbb{E}[m_n(x)] - m(x) = o\left(\sqrt{\beta_n h_n^{-d}}\right). \quad (4.6) $$

$$ \text{2. If } a \in (0, \frac{\beta}{d+4}), \text{ then} \qquad \mathrm{Var}[m_n(x)] = o(h_n^4). \quad (4.7) $$

$$ \text{If } a \in [\frac{\beta}{d+4}, 1), \text{ then} \qquad \mathrm{Var}[m_n(x)] = \frac{\beta_n}{h_n^d} \psi_n^{-1} \frac{\mathbb{E}[T^2 | X = x]}{2 - (\beta - ad)\xi_\beta} f(x) R(\mathbf{K}) + o(\beta_n h_n^{-d}). \quad (4.8) $$

Our main result rests upon the following theorem, which provides us the bias and the variance of $p_n$.
### 4.2.2 Bias and variance of $p_n$

**Theorem 4.2.** Let assumptions (A1), (A5) and (A7) hold. Hence we obtain

$$1. \text{ If } a \in (0, \frac{\beta}{d+4}], \text{ then} \\ \mathbb{E}[p_n(x)] - p(x) = \frac{1}{2(1 - 2a\xi_\beta)} \frac{h_n^2}{f(x)} \sum_{j=1}^{d} \mu_j(\mathbf{K}) (m_{jj}^{(2)}(x) - p(x)f_{jj}^{(2)}(x)) + o(h_n^2). \quad (4.9)$$

$$\text{If } a \in \left(\frac{\beta}{d+4}, 1\right), \text{ then} \\ \mathbb{E}[p_n(x)] - p(x) = o\left(\sqrt{\beta_n h_n^{-d}}\right). \quad (4.10)$$

$$2. \text{ If } a \in (0, \frac{\beta}{d+4}), \text{ then} \\ \mathrm{Var}[p_n(x)] = o(h_n^4). \quad (4.11)$$

$$\text{If } a \in \left[\frac{\beta}{d+4}, 1\right), \text{ then} \\ \mathrm{Var}[p_n(x)] = \frac{\beta_n}{h_n^d} \frac{\psi_n^{-1}}{2 - (\beta - ad)\xi_\beta} \frac{R(\mathbf{K})}{f(x)} \left(\mathbb{E}[T^2|X=x] - \psi p^2(x)\right) + o\left(\frac{\beta_n}{h_n^d}\right). \quad (4.12)$$

The bias and the variance of the estimator $p_n$ defined by the stochastic approximation algorithm (4.3) then heavily depend on the choice of the stepsizes $(\beta_n)$.

Now, let us state the following theorem which yields the asymptotic normality of the proposed multivariate recursive regression estimator under missing data $p_n$ denoted in (4.3).

### 4.2.3 Asymptotic normality of $p_n$

**Theorem 4.3.** Let assumptions (A1), (A5) and (A7) hold. We therefore have

$$\text{If there exists } c \ge 0 \text{ such that } \beta_n^{-1} h_n^{d+4} \xrightarrow{n \to +\infty} c, \text{ then} \\ \sqrt{\beta_n^{-1} h_n^d \psi_n} (p_n(x) - p(x)) \xrightarrow{n \to +\infty} \mathcal{N}\left(\sqrt{cM}(x), \Sigma(x)\right), \quad (4.13)$$

with

$$M(x) = \frac{1}{2(1 - 2a\xi_\beta)f(x)} \sum_{j=1}^{d} \mu_j(\mathbf{K}) (m_{jj}^{(2)}(x) - p(x)f_{jj}^{(2)}(x))$$

and

$$\Sigma(x) = \frac{1}{2 - (\beta - ad)\xi_\beta} \frac{R(\mathbf{K})}{f(x)} \left( \mathbb{E}[T^2 | X = x] - \psi p^2(x) \right),$$

## 4.3 Optimal choice of the stepsizes

In order to measure the asymptotic performance of the proposed recursive kernel regression estimator under missing data $p_n$ and to be able to use a data-driven bandwidth selection procedure, through proposing an asymptotic unbiased estimators of the unknown quantities, we consider the Mean Weighted Integrated Squared Error (MWISE), where the weight function is selected to be equal to $f^3(x)$. This choice was motivated by the fact that we can propose an asymptotic unbiased kernel estimator for the unknown quantities, which will appear in the MWISE as reported previously in Slaoui (2016), and which shall be detailed later.
### 4.3.1 Asymptotic expressions of MWISE of $p_n$

The MWISE of the estimator $p_n$ is determined by,

$$MWISE[p_n] = \int_{\mathbb{R}^d} (\mathbb{E}[p_n(x)] - p(x))^2 f^3(x) dx + \int_{\mathbb{R}^d} \mathrm{Var}[p_n(x)] f^3(x) dx. \quad (4.14)$$

In the sequel, we will need the following notations

$$I_1 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) \right)^2 f(x) dx, \quad I_2 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) \right) \left( \sum_{j=1}^d \mu_j(\mathbf{K}) f_{jj}^{(2)}(x) \right) p(x) f(x) dx, \\ I_3 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^d \mu_j(\mathbf{K}) f_{jj}^{(2)}(x) \right)^2 p^2(x) f(x) dx, \quad I_4 = \int_{\mathbb{R}^d} \mathbb{E}[T^2 | X = x] f^2(x) dx, \quad I_5 = \int_{\mathbb{R}^d} p^2(x) f^2(x) dx.$$

**Proposition 4.4.** For simplicity, we first set

$$C_1 = \frac{I_1 - 2I_2 + I_3}{(1 - 2a\xi_\beta)^2} \quad \text{and} \quad C_2 = \frac{I_4 - \psi I_5}{2 - (\beta - ad)\xi_\beta}.$$

It follows that

$$MWISE[p_n] = \begin{cases} \frac{1}{4}C_1h_n^4 + o(h_n^4) & \text{if } a \in (0, \frac{\beta}{d+4}) \\ C_2R(\mathbf{K})\beta_n h_n^{-d}\psi_n^{-1} + \frac{1}{4}C_1h_n^4 + o(h_n^4) & \text{if } a = \frac{\beta}{d+4} \\ C_2R(\mathbf{K})\beta_n h_n^{-d}\psi_n^{-1} + o(\beta_n h_n^{-d}) & \text{if } a \in (\frac{\beta}{d+4}, 1) \end{cases}.$$

The corollary below ensures that the bandwidth which minimizes the MWISE of $p_n$ depends on the choice of the stepsizes $(\beta_n)$ and then the corresponding MWISE depends in turn on $(\beta_n)$.

**Corollary 4.5.** Let assumptions (A1), (A5) and (A7) hold. To minimize the MWISE of $p_n$, the bandwidth $(h_n)$ must be equal to

$$\left( d^{\frac{1}{d+4}} \left( \frac{C_2}{C_1} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} \beta_n^{\frac{1}{d+4}} \psi_n^{\frac{-1}{d+4}} \right).$$

Then, the corresponding MWISE is estimated in terms of

$$MWISE[p_n] = \frac{(d+4)}{4d^{\frac{d}{d+4}}} C_1^{\frac{d}{d+4}} C_2^{\frac{4}{d+4}} R(\mathbf{K})^{\frac{4}{d+4}} \beta_n^{\frac{4}{d+4}} \psi_n^{\frac{-4}{d+4}} + o\left(\beta_n^{\frac{4}{d+4}}\right).$$

The following corollary is presented in the special case, where $(\beta_n)$ is chosen as $(\beta_n) = (\beta_0 n^{-1})$. We can check easily that the optimal choice of $\beta_0$ is obtained by getting $\beta_0$ equal to 1.

**Corollary 4.6.** Let assumptions (A1), (A5) and (A7) hold. To minimize the MWISE of $p_n$, we must choose the stepsize $(\beta_n)$ in $GS(-1)$ such that $\lim_{n \to +\infty} n\beta_n = 1$. Consequently, the optimal bandwidth $(h_n)$ must be equal to

$$\left( \left( \frac{d(d+2)}{2(d+4)} \right)^{\frac{1}{d+4}} \left( \frac{I_4 - \psi I_5}{I_1 - 2I_2 + I_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-1}{d+4}} \psi_n^{\frac{-1}{d+4}} \right). \quad (4.15)$$

Thus, the corresponding MWISE is provided by

$$MWISE[p_n] = \frac{(d+4)^{\frac{3d+8}{d+4}}}{4^{\frac{d+6}{d+4}} d^{\frac{d}{d+4}} (d+2)^{\frac{d+6}{d+4}}} (I_1 - 2I_2 + I_3)^{\frac{d}{d+4}} (I_4 - \psi I_5)^{\frac{4}{d+4}} R(\mathbf{K})^{\frac{4}{d+4}} n^{-\frac{4}{d+4}} \psi_n^{-\frac{4}{d+4}} \\ + o\left(n^{-\frac{4}{d+4}} \psi_n^{-\frac{4}{d+4}}\right).$$
## 4.4 Asymptotic properties of $\tilde{p}_n$

The main properties of the generalized non-recursive regression function estimator $\tilde{p}_n$ are displayed in the following proposition.

**Proposition 4.7.** Let assumptions (A₁) and (A₇) hold. Therefore, the bias and variance of Nadaraya-Watson's regression estimator are equal to:

$$ \mathbb{E}[\tilde{p}_n(x)] - p(x) = \frac{1}{2f(x)}h_n^2 \left( \sum_{j=1}^{d} \mu_j(\mathbf{K})m_{jj}^{(2)}(x) - p(x) \sum_{j=1}^{d} \mu_j(\mathbf{K})f_{jj}^{(2)}(x) \right) + o(h_n^2). $$

$$ \mathrm{Var}[\tilde{p}_n(x)] = \frac{1}{nh_n^d} \psi_n^{-1} \frac{1}{f(x)} \mathrm{Var}[T|X = x] R(\mathbf{K}) + o\left(\frac{1}{nh_n^d}\right). $$

It is inferred that

$$ \mathrm{MWISE}[\tilde{p}_n] = \frac{1}{4}(I_1 - 2I_2 + I_3)h_n^4 + \frac{1}{nh_n^d}\psi_n^{-1}(I_4 - \psi I_5)R(\mathbf{K}) + o\left(h_n^4 + \frac{1}{nh_n^d\psi_n}\right). $$

**Corollary 4.8.** Let assumptions (A₁) and (A₇) hold. To minimize the MWISE of $\tilde{p}_n$, the bandwidth ($h_n$) must be equal to

$$ \left( d^{\frac{1}{d+4}} \left( \frac{I_4 - \psi I_5}{I_1 - 2I_2 + I_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{1}{d+4}} \psi_n^{-\frac{1}{d+4}} \right). \quad (4.16) $$

Then, the corresponding MWISE is specified by

$$ \mathrm{MWISE}[\tilde{p}_n] = \frac{(d+4)}{4d^{\frac{d}{d+4}}} (I_4 - \psi I_5)^{\frac{4}{d+4}} (I_1 - 2I_2 + I_3)^{\frac{d}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{4}{d+4}} \psi_n^{-\frac{4}{d+4}} + o\left(n^{-\frac{4}{d+4}} \psi_n^{-\frac{4}{d+4}}\right). $$

It is obvious that, the use of such bandwidth (4.16), is not possible when we use real data. From this perspective, the next section is devoted to build up a data-driven bandwidth procedure, which will be helpful in practice.

## 4.5 Bandwidth selection

Within the framework of non-parametric kernel estimation, the choice of the smoothing parameter is crucial for the effective performance of the estimators. There are a myriad of data-driven bandwidth selection methods recorded in literature which can be divided into three broad classes: cross-validation techniques, plug-in methods, and the bootstrap approach. A detailed comparison of the three techniques is exhibited in Delaigle and Gijbels (2004). In this chapter, based on the previous work conducted by Slaoui (2014a), Slaoui (2014b), Slaoui (2016) for unidimensional data, we propose a second generation Plug-in bandwidth data-driven procedures in the multivariate data for regression estimation.

### 4.5.1 Plug-in bandwidth selection method

A widely used criterion stands for selecting a bandwidth that minimizes the estimate of the mean squared error, using the density function as a weight function. Altman and Leger (1995) proposed an efficient method of bandwidth selection, a plug-in estimate. Since the MWISE depends on unknown quantities $I_j, j = 1...5$, we suggest elaborating an asymptotic unbiased estimator of those quantities. For this purpose, we adopt the relation (1.14) and let $\mu(K) = \int_{\mathbb{R}} z^2 K(z) dz$ and

$$ I_j = \mu^2(K) I'_j, \quad j = 1, 2, 3, $$
where

$$I'_1 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^{d} m_{jj}^{(2)}(x) \right)^2 f(x) dx,$$

$$I'_2 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^{d} m_{jj}^{(2)}(x) \right) \left( \sum_{j=1}^{d} f_{jj}^{(2)}(x) \right) p(x) f(x) dx,$$

$$I'_3 = \int_{\mathbb{R}^d} \left( \sum_{j=1}^{d} f_{jj}^{(2)}(x) \right)^2 p^2(x) f(x) dx.$$

At this stage of analysis, in order to estimate the optimal bandwidth (4.15), we need to estimate $I_j$, $j = 1, \dots, 5$. For this purpose, we assume that $K_b$ is a kernel and $b_n$ is the associated bandwidth, such that $\delta = 2/5$, and $K_{b'}^{(2)}$ is the second derivative of a kernel $K_{b'}$ with the associated bandwidth $b'_n$ such that $\delta = 3/14$.

### Multivariate recursive kernel regression estimator under missing data $p_n$:

Here, we can state

$$m_n(x) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k h_k^{-d} Y_k \psi_k^{-1} \mathbf{K} \left( \frac{x - X_k}{h_k} \right) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k Y_k \psi_k^{-1} h_k^{-d} \prod_{i=1}^{d} K \left( \frac{x_i - X_{ki}}{h_k} \right)$$

and

$$f_n(x) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k h_k^{-d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k h_k^{-d} \prod_{i=1}^{d} K \left( \frac{x_i - X_{ki}}{h_k} \right).$$

### Estimation of $I_1$, $I_2$ and $I_3$:

We consider the following kernel estimators to estimate respectively $I_1$, $I_2$ and $I_3$:

$$\hat{I}'_1 = \frac{Q_n^2}{n} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k}}^n Q_j^{-1} Q_k^{-1} \beta_j \beta_k b_j^{-(d+2)} b_k^{-(d+2)} \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_j'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_j} \right) \right] \\ \times \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_k'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \right] Y_j \psi_j^{-1} Y_k \psi_k^{-1},$$

$$\hat{I}'_2 = \frac{Q_n^2}{n} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k}}^n Q_j^{-1} Q_k^{-1} \beta_j \beta_k b_j^{-(d+2)} b_k^{-(d+2)} \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_j'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_j} \right) \right] \\ \times \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_k'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \right] Y_j \psi_j^{-1} Y_i \psi_i^{-1},$$
$$ \hat{I}'_3 = \frac{Q_n^2}{n} \sum_{\substack{i,j,k,m=1 \\ i \neq j \neq k \neq m}}^n Q_j^{-1} Q_k^{-1} \beta_j \beta_k b_j^{-(d+2)} b_k^{-(d+2)} \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_j'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{jl}}{b_j} \right) \right] \\ \times \left[ \sum_{t=1}^d K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_k'} \right) \prod_{\substack{l=1 \\ l \neq t}}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) \right] Y_i \psi_i^{-1} Y_m \psi_m^{-1}, $$

Therefore, we obtain

$$ \hat{I}_i = \mu^2(\mathbf{K}) \hat{I}'_i, \quad i = 1 \dots 3. $$

### Estimation of $I_4$ and $I_5$:

We consider the following kernel estimators to estimate respectively $I_4$ and $I_5$:

$$ \hat{I}_4 = \frac{\Pi_n}{n} \sum_{\substack{i,k=1 \\ i \neq k}}^n \Pi_k^{-1} \gamma_k b_k^{-d} \prod_{l=1}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) Y_i^2 \psi_i^{-2}, $$

and

$$ \hat{I}_5 = \frac{Q_n}{n} \sum_{\substack{i,k=1 \\ i \neq k}}^n Q_k^{-1} \beta_k b_k^{-d} \prod_{l=1}^d K_b \left( \frac{X_{il} - X_{kl}}{b_k} \right) Y_i \psi_i^{-1} Y_k \psi_k^{-1}. $$

It follows that, the plug-in bandwidth selection estimator of (4.15) is expressed by

$$ (h_n) = \left( \left( \frac{d(d+2)}{2(d+4)} \right)^{\frac{1}{d+4}} \left( \frac{\hat{I}_4 - \psi \hat{I}_5}{\hat{I}_1 - 2\hat{I}_2 + \hat{I}_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-1}{d+4}} \psi_n^{-\frac{1}{d+4}} \right), \quad (4.17) $$

Then, the plug-in estimator of $MWISE[p_n]$ is equal to

$$ \widehat{MWISE}[p_n] = \frac{(d+4)^{\frac{3d+8}{d+4}}}{4^{\frac{d+6}{d+4}} d^{\frac{d}{d+4}} (d+2)^{\frac{d+6}{d+4}}} (\hat{I}_1 - 2\hat{I}_2 + \hat{I}_3)^{\frac{d}{d+4}} (\hat{I}_4 - \psi \hat{I}_5)^{\frac{4}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{\frac{-4}{d+4}} \psi_n^{-\frac{4}{d+4}} \\ + o\left(n^{\frac{-4}{d+4}} \psi_n^{-\frac{4}{d+4}}\right). $$

### Multivariate non-recursive kernel regression estimator under missing data $\tilde{p}_n$:

Here, we can state

$$ \tilde{m}_n(x) = \frac{1}{nh_n^d} \sum_{k=1}^n Y_k \psi_k^{-1} \mathbf{K} \left( \frac{x - X_k}{h_n} \right) = \frac{1}{nh_n^d} \sum_{k=1}^n Y_k \psi_k^{-1} \prod_{i=1}^d K \left( \frac{x_i - X_{ki}}{h_k} \right) $$

and

$$ \tilde{f}_n(x) = \frac{1}{nh_n^d} \sum_{k=1}^n \mathbf{K} \left( \frac{x - X_k}{h_n} \right) = \frac{1}{nh_n^d} \sum_{k=1}^n \prod_{i=1}^d K \left( \frac{x_i - X_{ki}}{h_k} \right). $$

In order to estimate the optimal bandwidth (4.16), we need to estimate $I_j$, $j = 1 \dots 5$.

### Estimation of $I_1$, $I_2$ and $I_3$:

We consider the following kernel estimators to estimate respectively $I_1$, $I_2$ and $I_3$:
$$ \tilde{I}'_1 = \frac{1}{n^3 b_n^{2(d+2)}} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k}}^n \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] Y_j \psi_j^{-1} Y_k \psi_k^{-1}, $$

$$ \tilde{I}'_2 = \frac{1}{n^3 b_n^{2(d+2)}} \sum_{\substack{i,j,k=1 \\ i \neq j \neq k \neq m}}^n \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] Y_j \psi_j^{-1} Y_i \psi_i^{-1}, $$

$$ \tilde{I}'_3 = \frac{1}{n^4 b_n^{2(d+2)}} \sum_{\substack{i,j,k,m=1 \\ i \neq j \neq k \neq m}}^n \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{jt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{jl}}{b_n} \right) \right] \\ \times \left[ \sum_{t=1}^{d} K_{b'}^{(2)} \left( \frac{X_{it} - X_{kt}}{b_n'} \right) \prod_{\substack{l=1 \\ l \neq t}}^{d} K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) \right] Y_i \psi_i^{-1} Y_m \psi_m^{-1}, $$

Therefore, we obtain

$$ \tilde{I}_i = \mu^2(\mathbf{K}) \tilde{I}'_i, \quad i = 1 \dots 3. $$

### Estimation of $I_4$ and $I_5$:

We consider the following kernel estimators to estimate respectively $I_4$ and $I_5$:

$$ \tilde{I}_4 = \frac{1}{n^2 b_n^d} \sum_{\substack{i,k=1 \\ i \neq k}}^n \prod_{l=1}^d K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) Y_i^2 \psi_i^{-2}, $$

and

$$ \tilde{I}_5 = \frac{1}{n^2 b_n^d} \sum_{\substack{i,k=1 \\ i \neq k}}^n \prod_{l=1}^d K_b \left( \frac{X_{il} - X_{kl}}{b_n} \right) Y_i \psi_i^{-1} Y_k \psi_k^{-1}, $$

Hence, the plug-in bandwidth selection estimator of (4.16) is indicated by

$$ (h_n) = \left( d^{\frac{1}{d+4}} \left( \frac{\tilde{I}_4 - \psi \tilde{I}_5}{\tilde{I}_1 - 2\tilde{I}_2 + \tilde{I}_3} \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{1}{d+4}} \psi_n^{-\frac{1}{d+4}} \right), \quad (4.18) $$

It follows that, the plug-in non-recursive estimator of MWISE[$p_n$] is equal to

$$ \widetilde{MWISE}[\tilde{p}_n] = \frac{(d+4) 5}{4 d^{\frac{d}{d+4}}} \left( \tilde{I}_4 - \psi \tilde{I}_5 \right)^{\frac{4}{d+4}} \left( \tilde{I}_1 - 2 \tilde{I}_2 + \tilde{I}_3 \right)^{\frac{1}{d+4}} R(\mathbf{K})^{\frac{1}{d+4}} n^{-\frac{4}{d+4}} \psi_n^{-\frac{4}{d+4}} \\ + o \left( n^{-\frac{4}{d+4}} \psi_n^{-\frac{4}{d+4}} \right). $$
## 4.6 Application to the handwritten word production

Research on the handwritten word production aims to describe the cognitive processes and mental representations mobilized when a human being prepares to handwrite a word from an idea of Perret and Olive (2019). One of the most widely used tasks to experimentally explore these issues is object naming. Participants have to produce words corresponding to the names of a set of drawings in handwriting as quickly as possible. It is generally accepted that the handwritten objects naming involves four levels of processing Perret and Bonin (2019). First, a perceptual analysis of the visual input is performed, which results in activation of stored structural knowledge about the object. A second processing level corresponds to the retrieval of semantic/conceptual information. The lexical selection level makes orthographic word form information available. Eventually, the motoric programming level allows the access to motoric codes corresponding to each produced letter.

These theoretical propositions concerning the cognitive processes and representations involved in the handwritten object naming stem from studies aiming at finding predictors of reaction times (RTs hereafter), i.e., the time between the presentation of the image and the first graphic movement (e.g., Bonin et al. (2002); Perret and Laganaro (2013)).

Four factors have been reported to significantly influence RTs, each of which allows indexing a specific processing level.

**Image Agreement (IA)** captures the similarity between structural representations stored in memory and the visual characteristics of an object's drawing. This factor has extensive influence in terms of the perceptual analysis. The IA is measured on a Like rt scale, generally in five points, from '1 - weakly similar' to '5 - strongly similar'. A negative linear relationship is observed between this variable and the RTs (see Bonin et al. (2002); Perret and Laganaro (2013)).

**Image variability (Ivar or Image ability)** is designed to index the 'richness' of semantic representations. Like AI, it is rated on a 5-point scale, from 1 = few images to 5 = many images. A negative linear relationship is reported between handwritten RTs and Ivar (see Bonin et al. (2002); Perret and Laganaro (2013)).

**Name agreement (NA)** refers to the degrees of agreement on the use of a specific label for an image, measured using an entropy measure (h-index). A positive linear relationship is reported between RTs and the h-index (see Bonin et al. (2002); Perret and Laganaro (2013)). NA indexes the influence of the number of correct alternative names existing for an image (e.g., couch = > sofa). Latencies would be more or less impacted by the time needed to manage the competition between the higher or lower number of alternatives during lexical access.

Finally, the influence of age-limited learning (Age of Acquisition, AoA) has been systematically emphasized in studies on the predictors of handwritten RTs (see Bonin et al. (2002); Perret and Laganaro (2013)). AoA is usually measured using a Like rt scale (from 1 = learned at 0-3 years to 5 = learned at age +12, with 3-year bands in between), with a population of young adults who are asked to estimate the age at which they learned the proposed word. A positive linear relationship is observed between the RTs and the rated values of AoA (see Bonin et al. (2002); Perret and Laganaro (2013)). Experimental work Perret et al. (2014) suggests that this variable influences the orthographic wordform encoding processes.

### Problematic:

The major target of this work is to classify the participants in groups of clusters. From this perspective, we first have to predict the regression function, i.e the relation between the variable T = RTs and the four covariates $X_1 = H$, $X_2 = IA$, $X_3 = Ivar$ and $X_4 = AoA$. Since the response variable RTs is subject of missing data, we should introduce a correction variable Y := CRTs defined as followed:

$$Y_i = T_i \times 1\{T_i \text{ is observed}\}.$$
Here, we have Np individual estimators of each participant $\hat{Y}_1, \dots, \hat{Y}_{Np}$ ($\tilde{Y}_1, \dots, \tilde{Y}_{Np}$) and a general estimator $\hat{Y}^g$ ($\tilde{Y}^g$) which estimates the whole database of Np participants.

It's worth noting that, for each participant/covariate behavior test, we invested a different method for bandwidth selection, namely the plug-in univariate selection for multivariate data. This implies that, instead of opting for a single value of bandwidth $h_n$, we considered a vector $h_{n1}, \dots, h_{nd}$, an individual choice of bandwidth for each covariate. Then, for the recursive case, we have a matrix of bandwidths:

$$H = \begin{pmatrix} h_{11} & \dots & h_{1d} \\ \vdots & \ddots & \vdots \\ h_{n1} & \dots & h_{nd} \end{pmatrix}.$$

We denote by $p_i^*$ the reference regression vector and by $p_i$ the test regression. Thus, we calculate the two following measures:

* The Mean Squared Error: $MSE = \frac{1}{n} \sum_{i=1}^{n} (p_i - p_i^*)^2$.

* The Mean Relative Error: $MRE = \frac{1}{n} \sum_{i=1}^{n} \left| \frac{p_i - p_i^*}{p_i^*} \right|$.

**Algorithm 2** $X_1, \dots, X_4$ are the covariates such that $X_1 = H$, $X_2 = IA$, $X_3 = Ivar$ and $X_4 = AoA$, $Y$ is the response variable with $Y = CRTs$, $K$ is the Gaussian kernel, $n$ the number of items and $Np$ is the number of participants.

**Input:** $Y$, $X_1, \dots, X_4$, $K$, $n$ and $Np$.

1: A choice value for the recursive bandwidth vectors $h_1, \dots, h_n$. (resp. the non-recursive bandwidth values $h_n$) using the plug-in approach provided in (4.15) (resp. (4.18)).

2: The choice of the stepsize $(\beta_n) = (n^{-1})$ (then, $(Q_n) = (n^{-1})$).

3: An arbitrary sampling vectors $x_1, \dots, x_4$.

4: The estimation of $\psi$ is carried out according to the algorithm proposed in Slaoui (2017).

5: **for** $l = 1, \dots, Np$ **do**

6: $$ \hat{Y}_l = \frac{\sum_{k=1+(l-1)n}^{ln} k\beta_k Y_k \psi_k^{-1} \prod_{i=1}^4 h_{ki}^{-1} \prod_{i=1}^4 K\left(\frac{x_i - X_{ki}}{h_{ki}}\right)}{\sum_{k=1+(l-1)n}^{ln} k\beta_k \prod_{i=1}^4 h_{ki}^{-1} \prod_{i=1}^4 K\left(\frac{x_i - X_{ki}}{h_{ki}}\right)} $$ for the multivariate recursive kernel

regression estimator (resp. $\tilde{Y}_l = \frac{\sum_{k=1+(l-1)n}^{ln} Y_k \psi_k^{-1} \prod_{i=1}^4 K\left(\frac{x_i - X_{ki}}{h_{ki}}\right)}{\sum_{k=1+(l-1)n}^{ln} \prod_{i=1}^4 K\left(\frac{x_i - X_{ki}}{h_{ki}}\right)}$ for the multivariate non-recursive kernel regression estimator).

7: **end for**

**output:** $\hat{Y}_1, \dots, \hat{Y}_{Np}$ and $\tilde{Y}_1, \dots, \tilde{Y}_{Np}$.
<figure><img src="image_50.png" /><figcaption>Figure 4.1: Participants’ behavior representations, the regression between the reactions time variable $Y = CRTs$ and each covariate ($X_1 = H$, $X_2 = IA$, $X_3 = Ivar$ and $X_4 = AoA$) with the entire database (a total of 137 participants).</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure 4.2: Box-plot of the relative error estimation of both considered estimators, the recursive one on the left and the non-recursive one on the right.</figcaption></figure>

<table><caption>Table 4.1: Quantitative comparison between the mean relative error of the multivariate non-recursive Nadaraya-Watson’s regression estimator (Non-recursive) and the proposed multivariate recursive kernel regression estimator (Recursive) with stepsize ($\beta_n$) = ($n^{-1}$) through a plug-in method.</caption><thead><tr><th>Relative Error</th><th>Min</th><th>1st Qu</th><th>Median</th><th>Mean</th><th>3rd Qu</th><th>Max</th></tr></thead><tbody><tr><td>Recursive</td><td>0.0000014</td><td>0.1208883</td><td>0.2643889</td><td>0.3837469</td><td>0.5322441</td><td>2.4902917</td></tr><tr><td>Non-recursive</td><td>0.0000103</td><td>0.1412407</td><td>0.3301375</td><td>0.4748705</td><td>0.9999976</td><td>2.4525833</td></tr></tbody></table>

Let us underline that in order to classify participants in groups, we shall use the *MSE* as a reference vector. Thus, we resort to the k-means method to specify the maximum number of needed clusters.
**Algorithm 3** Participants classification algorithm:

Y is the response variable with $Y = CRTs$, $\hat{Y}_1, \dots, \hat{Y}_{Np}$ are the predicted multivariate recursive kernel regression estimators and $\tilde{Y}_1, \dots, \tilde{Y}_{Np}$ are the predicted multivariate non-recursive kernel regression estimators.

**Input:** $Y, \hat{Y}_1, \dots, \hat{Y}_{Np}$ and $\tilde{Y}_1, \dots, \tilde{Y}_{Np}$.

1. Start with writing $Y$ in a matrix form participant per participant.

2. **for** $l = 1, \dots, Np$ **do**

3. $MSER_l = \frac{1}{n} \sum_{i=1}^{n} (\hat{Y}_i - Y_i)^2$. for the recursive estimation (resp. $MSET_l = \frac{1}{n} \sum_{i=1}^{n} (\tilde{Y}_i - Y_i)^2$. for the non-recursive estimation).

4. **end for**

5. A classification of the remote distance through kmeans package in R.

**output:** The classification list using both considered estimators.

<figure><img src="image_52.png" /><figcaption>Figure 4.3: The elbow method of selecting the optimal number of clusters ($k = 3$) for K-means clustering on the MSE vector.</figcaption></figure>

**Result Analysis:**

Departing from figure 4.2 and table 4.1, we deduce that the proposed recursive estimator outperforms the non-recursive one in terms of mean relative error estimation. Meanwhile, figures 4.1 and 4.3 indicate that it is advisable to consider three clusters. As far as written production behavior is concerned, this implies that the classification procedure suggests three clusters to measure the distance of each participant from the reference. In other words, three forms of variation can be observed when participants have to write the label of a drawing. Further exploration of the available characteristics of the participants suggests that such anthropological factors as the age and gender do not account for the result of clustering. Descriptive analysis of executive function task data suggests that there are differences between the three groups of participants. This indicates that the variations would be interpreted in part by the participants' cognitive processing ability and by differences in the mobilization of participants' executive functions. Studies based upon procedures for fitting reaction time distributions with ex-Gaussian-type probability density distributions (convolution of a normal and exponential law) have corroborated the role of these executive functions in simple tasks (e.g., Schmiedek et al. (2007); Unsworth et al. (2010)). Our analyses yield that this result can be extended to more complex activities such as written production. Eventually, this work confirms the significance of the use of non-parametric regressions for modeling behavior in experimental psychology area.
## 4.7 Conclusion

In this research work, we elaborated a multivariate recursive regression estimator under missing data. We first investigated the asymptotic properties of the proposed estimator by providing the bias as well as the variance in order to demonstrate that our estimator asymptotically follows a normal distribution. Subsequently, we compared our recursive estimator with the non-recursive multivariate Nadaraya-Watson's regression estimator using the plug-in bandwidth selection approach. In our application of real dataset, and for all the cases, the proposed estimator (4.3) with stepsize ($\beta_n$) = ($n^{-1}$) yielded smaller MSE and MRE compared to the non-recursive Nadaraya-Watson's estimator.

As part of the application, it was possible to estimate the response variable RTs (Reaction Times) according to the other covariates through classifying the participants into clusters of membership according to their approximation to the real value of RTs.

To conclude, the use of the multivariate recursive kernel regression estimator under missing data enabled us to obtain better results compared to the multivariate non-recursive kernel regression estimator under missing data. With an appropriate choice of the bandwidth, we depicted that our proposed estimator is closer to the true regression function than the non-recursive one.

## 4.8 Proofs

Throughout this section, we use the following notations:

$$ \mathcal{Z}_n(x) = h_n^{-d} Y_n \psi_n^{-1} \mathbf{K} \left( \frac{x - X_n}{h_n} \right) \quad \text{and} \quad \mathcal{W}_n(x) = h_n^{-d} \mathbf{K} \left( \frac{x - X_n}{h_n} \right), \text{ for all } x \in \mathbb{R}^d. $$

*Proof of Proposition 4.1.*

We have

$$ \begin{aligned} m_n(x) - m(x) &= (1 - \beta_n)m_{n-1}(x) + \beta_n \mathcal{Z}_n(x) - m(x) \\ &= (1 - \beta_n)[m_{n-1}(x) - m(x)] + \beta_n[\mathcal{Z}_n(x) - m(x)] \\ &= \prod_{i=1}^{n} (1 - \beta_i)[m_0(x) - m(x)] + \sum_{k=1}^{n-1} \prod_{i=k+1}^{n} (1 - \beta_i)\beta_k(\mathcal{Z}_k(x) - m(x)) + \beta_n(\mathcal{Z}_n(x) - m(x)) \\ &= Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k(\mathcal{Z}_k(x) - m(x)) + Q_n[m_0(x) - m(x)]. \end{aligned} $$

It follows that,

$$ \mathbb{E}[m_n(x)] - m(x) = Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k (\mathbb{E}[\mathcal{Z}_k(x)] - m(x)) + Q_n[m_0(x) - m(x)]. \quad (4.19) $$

Moreover, we have

$$ \begin{aligned} \mathbb{E}[\mathcal{Z}_k(x)] &= h_k^{-d} \psi_k^{-1} \mathbb{E} \left[ Y_k \mathbf{K} \left( \frac{x - X_k}{h_k} \right) \right] \\ &= h_k^{-d} \psi_k^{-1} \mathbb{E} \left[ T_k \mathbf{1}_{\{T_k = Y_k\}} \mathbf{K} \left( \frac{x - X_k}{h_k} \right) \right] \\ &= h_k^{-d} \psi_k^{-1} \mathbb{E} \left[ \mathbf{1}_{\{T_k = Y_k\}} \right] \int_{\mathbb{R}^d} \mathbb{E}[T | X = y] \mathbf{K} \left( \frac{x - y}{h_k} \right) f(y) dy \\ &= h_k^{-d} \int_{\mathbb{R}^d} \mathbf{K} \left( \frac{x - y}{h_k} \right) m(y) dy. \end{aligned} $$
Since we have $\int_{\mathbb{R}^d} \mathbf{K}(z)dz = 1$, we infer that

$$ 
\begin{aligned}
\mathbb{E}[\mathcal{Z}_k(x)] - m(x) &= \int_{\mathbb{R}^d} h_k^{-d} \mathbf{K} \left( \frac{x-y}{h_k} \right) m(y) dy - \int_{\mathbb{R}^d} \mathbf{K}(y) m(x) dy \\
&= \int_{\mathbb{R}^d} \mathbf{K}(z) [m(x-zh_k) - m(x)] dz.
\end{aligned}
$$

A Taylor expansion of $m$ around $x$ ensures that

$$ 
\begin{aligned}
\mathbb{E}[\mathcal{Z}_k(x)] - m(x) &= \int_{\mathbb{R}^d} \mathbf{K}(z) [m(x-zh_k) - m(x)] dz \\
&= \int_{\mathbb{R}^d} \mathbf{K}(z) \left[ \sum_{i=1}^d \frac{\partial m}{\partial x_i}(x)z_i h_k + \int_0^1 (1-t) \sum_{i,j=1}^d \frac{\partial^2 m}{\partial x_i \partial x_j}(x-tzh_k)z_i z_j h_k^2 dt \right] dz \\
&= h_k \sum_{i=1}^d \frac{\partial m}{\partial x_i}(x) \int_{\mathbb{R}^d} \mathbf{K}(z) z_i dz + h_k^2 \sum_{i,j=1}^d \int_{\mathbb{R}^d} \int_0^1 (1-t) \frac{\partial^2 m}{\partial x_i \partial x_j}(x-tzh_k)z_i z_j \mathbf{K}(z) dt dz \\
&= \frac{h_k^2}{2} \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) + h_k^2 \eta_k(x).
\end{aligned}
$$

$$ \text{where } \eta_k(x) = \sum_{i,j=1}^d \int_{\mathbb{R}^d} \int_0^1 (1-t) [m_{ij}^{(2)}(x-tzh_k) - m_{ij}^{(2)}(x)] z_i z_j \mathbf{K}(z) dt dz. $$

Owing to the fact that $m_{ij}^{(2)}$ is bounded and continuous at $x$ for all $i, j \in \{1, \dots, d\}$, we thus get

$$ 
\begin{aligned}
\mathbb{E}[m_n(x)] - m(x) &= Q_n \sum_{k=1}^n Q_k^{-1} \beta_k (\mathbb{E}[\mathcal{Z}_k(x)] - m(x)) + Q_n[m_0(x) - m(x)] \\
&= Q_n \sum_{k=1}^n Q_k^{-1} \beta_k \left( \frac{h_k^2}{2} \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) + h_k^2 \eta_k(x) \right) + Q_n[m_0(x) - m(x)] \\
&= \frac{1}{2} \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) Q_n \sum_{k=1}^n Q_k^{-1} \beta_k h_k^2 + Q_n \sum_{k=1}^n Q_k^{-1} \beta_k h_k^2 \eta_k(x) + Q_n[m_0(x) - m(x)].
\end{aligned}
$$

For the case $a \le \beta/(d+4)$, we have $\lim_{n \to +\infty} (n\beta_n) > 2a$ and then $1-2a\xi_\beta > 0$. The application of lemma 1.2 ensures that

$$ 
\begin{aligned}
\mathbb{E}[m_n(x)] - m(x) &= \frac{1}{2} \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) Q_n \sum_{k=1}^n Q_k^{-1} \beta_k h_k^2 + Q_n \sum_{k=1}^n Q_k^{-1} \beta_k o(h_k^2) + O(Q_n) \\
&= \frac{h_n^2}{2(1-2a\xi_\beta)} \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) + o(h_n^2).
\end{aligned}
$$

We infer that

$$ \mathbb{E}[m_n(x)] - m(x) = \frac{1}{2(1-2a\xi_\beta)} h_n^2 \sum_{j=1}^d \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) + o(h_n^2). $$

For the case $a > \beta/(d+4)$, we have $\lim_{n \to +\infty} (n\beta_n) > \frac{\beta-a}{2}$, which yields $h_n^2 = o(\sqrt{\beta_n h_n^{-d}})$. Hence,
the application of lemma 1.2 ensures that

$$
\begin{align*}
\mathbb{E}[m_n(x)] - m(x) &= \frac{1}{2} \sum_{j=1}^{d} \mu_j(\mathbf{K}) m_{jj}^{(2)}(x) Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k o\left(\sqrt{\beta_k h_k^{-d}}\right) + Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k o\left(\sqrt{\beta_k h_k^{-d}}\right) \\
&= o\left(\sqrt{\beta_n h_n^{-d}}\right).
\end{align*}
$$

As a matter of fact, the result can be expressed as

$$
\mathbb{E}[m_n(x)] - m(x) = o\left(\sqrt{\beta_n h_n^{-d}}\right).
$$

Let us now compute the variance of $m_n(x)$. We state

$$
\begin{align*}
\mathrm{Var}[m_n(x)] &= \mathrm{Var}[Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k \mathcal{Z}_k(x)] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 \mathrm{Var}[\mathcal{Z}_k(x)] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 (\mathbb{E}[\mathcal{Z}_k^2(x)] - \mathbb{E}[\mathcal{Z}_k(x)]^2).
\end{align*}
$$

Moreover, we have

$$
\begin{align*}
\mathbb{E}[\mathcal{Z}_k^2(x)] &= \int_{\mathbb{R}^d} h_k^{-2d} \psi_k^{-2} \mathbb{E}[T^2 | X = y] \psi_k \mathbf{K}^2 \left( \frac{x-y}{h_k} \right) f(y) dy \\
&= \int_{\mathbb{R}^d} h_k^{-d} \psi_k^{-1} \mathbf{K}^2 (z) \mathbb{E}[T^2 | X = x - zh_k] f(x - zh_k) dz.
\end{align*}
$$

Hence, the Taylor's expansion for the function $x \mapsto \mathbb{E}[T^2|X=x]f(x) = \int_{\mathbb{R}} y^2h(x,y)dy$ ensures
that

$$
\mathbb{E}[\mathcal{Z}_k^2(x)] = h_k^{-d} \psi_k^{-1} \left[ \mathbb{E}[T^2 | X = x] f(x) \int_{\mathbb{R}^d} \mathbf{K}^2(z) dz + \nu_k(x) \right].
$$

Thus,

$$
\begin{align*}
\mathrm{Var}[m_n(x)] &= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 \left[ \mathbb{E}[T^2 | X = x] \int_{\mathbb{R}^d} h_k^{-d} \psi_k^{-1} \mathbf{K}^2(z) f(x - zh_k) dz \right. \\
&\qquad \left. - \left( \int_{\mathbb{R}^d} \mathbf{K}(z) m(x - zh_k) dz \right)^2 \right] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \psi_k^{-1} \left[ \mathbb{E}[T^2 | X = x] f(x) \int_{\mathbb{R}^d} \mathbf{K}^2(z) dz + \nu_k(x) - h_k^d \psi_k \eta_k(x) \right],
\end{align*}
$$

$$
\text{where } \nu_k(x) = \int_{\mathbb{R}^d} \mathbf{K}^2(z) [\mathbb{E}[T^2 | X = x - zh_k] f(x - zh_k) - \mathbb{E}[T^2 | X = x] f(x)] dz
$$

$$
\text{and } \eta_k(x) = \left( \int_{\mathbb{R}^d} \mathbf{K}(z) m(x - zh_k) dz \right)^2
$$
For the case $a \ge \beta/(d+4)$, we have $\lim_{n \to +\infty} (n\beta_n) > \frac{\beta - ad}{2}$ and then $1 - 2a\xi_\beta > 0$. Since we have $\lim_{k \to +\infty} \nu_k(x) = 0$ and $\lim_{k \to +\infty} h_k \eta_k(x) = 0$, then the application of lemma 1.2 ensures that

$$
\begin{align*}
\mathrm{Var}[m_n(x)] &= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \psi_k^{-1} \left[ \mathbb{E}[T^2 | X = x] f(x) R(\mathbf{K}) + \nu_k(x) - h_k^d \eta_k(x) \right] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \psi_k^{-1} \left[ \mathbb{E}[T^2 | X = x] f(x) R(\mathbf{K}) + o(1) \right] \\
&= \frac{\mathbb{E}[T^2 | X = x]}{2 - (\alpha - ad)\xi_\beta} \frac{\beta_n}{h_n} \psi_n^{-1} [f(x) R(\mathbf{K}) + o(1)].
\end{align*}
$$

Therefore, the result is provided by

$$
\mathrm{Var}[m_n(x)] = \frac{\mathbb{E}[T^2 | X = x]}{2 - (\alpha - a)\xi_\beta} \frac{\beta_n}{h_n} \psi_n^{-1} f(x) R(\mathbf{K}) + o\left(\frac{\beta_n}{h_n}\right).
$$

For the case $a < \beta/(d+4)$, we have $\lim_{n \to +\infty} (n\beta_n) > 2a$ which yields $\beta_n h_n^{-d} = o(h_n^4)$. Then, the
application of lemma 1.2 ensures that

$$
\begin{align*}
\mathrm{Var}[m_n(x)] &= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \psi_k^{-1} \left[ \mathbb{E}[T^2 | X = x] f(x) R(\mathbf{K}) + o(1) \right] \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k o(h_k^4) \\
&= o(h_n^4).
\end{align*}
$$

□

Proof of Theorem 4.2.

For this proof, let us note that for $f_n(x) \neq 0$, we have

$$
p_n(x) - p(x) = A_n(x) \frac{f(x)}{f_n(x)}, \qquad (4.20)
$$

with

$$
A_n(x) = \frac{1}{f(x)} (m_n(x) - m(x)) - \frac{p(x)}{f(x)} (f_n(x) - f(x)). \quad (4.21)
$$

It follows from (4.20) that the asymptotic behavior of $p_n(x) - p(x)$ can be deduced from the one
of $A_n(x)$. Hence, we can state

$$
\mathbb{E}[A_n(x)] = \frac{1}{f(x)} (\mathbb{E}[m_n(x)] - m(x)) - \frac{p(x)}{f(x)} (\mathbb{E}[f_n(x)] - f(x)).
$$

Since we already have the bias of $m_n(x)$ as well as that of $f_n(x)$, and considering the fact that
$m(x) = p(x)f(x)$, then we just need to combine the results (4.5), (4.6), (2.7) and (2.8) in order
to obtain (4.9) and (4.10). Now, we have

$$
\mathrm{Var}[A_n(x)] = \frac{1}{(f(x))^2} \mathrm{Var}[m_n(x)] - \frac{(p(x))^2}{(f(x))^2} \mathrm{Var}[f_n(x)] - 2 \frac{p(x)}{(f(x))^2} \mathrm{Cov}(m_n(x), f_n(x)).
$$
Let us now compute the covariance between $m_n(x)$ and $f_n(x)$. Indeed, we have

$$
\begin{align*}
\mathrm{Cov}(m_n(x), f_n(x)) &= \mathrm{Cov} \left( Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k Y_k \psi_k^{-1} h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right), Q_n \sum_{i=1}^{n} Q_i^{-1} \beta_i h_i^{-d} \mathbf{K} \left( \frac{x-X_i}{h_i} \right) \right) \\
&= Q_n \sum_{k=1}^{n} Q_k^{-1} \beta_k Q_n \sum_{i=1}^{n} Q_i^{-1} \beta_i \mathrm{Cov} \left( Y_k \psi_k^{-1} h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right), h_i^{-d} \mathbf{K} \left( \frac{x-X_i}{h_i} \right) \right) \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 \mathrm{Cov} \left( Y_k \psi_k^{-1} h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right), h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right) \right) \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 \left( \mathrm{E} \left[ Y_k \psi_k^{-1} h_k^{-2d} \mathbf{K}^2 \left( \frac{x-X_k}{h_k} \right) \right] \right. \\
&\qquad \left. - \mathrm{E} \left[ Y_k \psi_k^{-1} h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right) \right] \mathrm{E} \left[ h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right) \right] \right) \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 \left( \mathrm{E}[T|X=x] f(x) R(\mathbf{K}) h_k^{-d} - \mathrm{E}[T|X=x] f^2(x) \right) + o(h_k^{-d}) \\
&= Q_n^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} (p(x) f(x) R(\mathbf{K}) + o(1)) \\
&= \frac{\beta_n h_n^{-d}}{2 - (\beta - ad)\xi_\beta} p(x) f(x) R(\mathbf{K}) + o(\beta_n h_n^{-d}). \tag{4.22}
\end{align*}
$$

Consequently, (4.11) and (4.12) follow from the combination of (4.7), (4.8), (2.9), (2.10) and (4.22). For the case $a \ge \beta/(d+4)$, we can deduce

$$
\begin{align*}
\mathrm{Var}[p_n(x)] &= \frac{1}{f(x)} \frac{\beta_n}{h_n^d} \psi_n^{-1} \frac{\mathrm{E}[T^2|X=x]}{2-(\beta-ad)\xi_\beta} R(\mathbf{K}) + \frac{p(x)^2}{f(x)} \frac{\beta_n}{h_n^d} \frac{1}{2-(\beta-ad)\xi_\beta} R(\mathbf{K}) \\
&\quad - 2 \frac{p(x)}{f(x)^2} \frac{\beta_n h_n^{-d}}{2-(\beta-ad)\xi_\beta} p(x) f(x) R(\mathbf{K}) + o(\beta_n h_n^{-d}) \\
&= \frac{\beta_n}{h_n^d} \frac{\psi_n^{-1}}{2-(\beta-ad)\xi_\beta} \frac{R(\mathbf{K})}{f(x)} \left( \mathrm{E}[T^2|X=x] - \psi p(x)^2 \right) + o(\beta_n h_n^{-d}).
\end{align*}
$$

□

### Proof of Theorem 4.3.

We have

$$
\begin{align*}
A_n(x) - \mathbb{E}[A_n(x)] &= \frac{1}{f(x)} [m_n(x) - \mathbb{E}[m_n(x)]] - \frac{p(x)}{f(x)} [f_n - \mathbb{E}[f_n]] \\
&= \frac{1}{f(x)} Q_n \sum_{k=1}^{n} (L_k(x) - \mathbb{E}[L_k(x)]),
\end{align*}
$$

with

$$L_k(x) = Q_k^{-1} \beta_k (\mathcal{Z}_k(x) - p(x) \mathcal{W}_k(x)).$$

In this proof, we note

$$S_k(x) = L_k(x) - \mathbb{E}[L_k(x)].$$

On the one hand, it's obvious that

$$p_n(x) - \mathbb{E}[p_n(x)] = \frac{1}{f(x)} Q_n \sum_{k=1}^{n} S_k(x). \tag{4.23}$$
On the other hand, we attempt to apply Lyapunov's theorem 1.14 for $S_k(x)$. For this reason, we consider

$$
\begin{align*}
\nu_n^2 &= \sum_{k=1}^{n} \mathrm{Var}[S_k(x)] \\
&= \sum_{k=1}^{n} \mathrm{Var}[L_k(x)] \\
&= \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 (\mathrm{Var}[\mathcal{Z}_k(x)] + p(x)^2 \mathrm{Var}[\mathcal{W}_k(x)] - 2p(x) \mathrm{cov}(\mathcal{Z}_k(x), \mathcal{W}_k(x))).
\end{align*}
$$

Moreover, we have

$$
\mathrm{Var}[\mathcal{Z}_k(x)] = h_k^{-d} \psi_k^{-1} \left( \mathbb{E}[T^2 | X = x] f(x) R(\mathbf{K}) + o(1) \right).
$$

$$
\mathrm{Var}[\mathcal{W}_k(x)] = h_k^{-d} (f(x)R(\mathbf{K}) + o(1)).
$$

$$
\mathrm{cov}(\mathcal{Z}_k(x), \mathcal{W}_k(x)) = h_k^{-d} (p(x)f(x)R(\mathbf{K}) + o(1)).
$$

Hence, by applying lemma 1.2, it can be inferred that

$$
\begin{align}
\nu_n^2 &= \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \psi_k^{-1} \left( \mathbb{E}[T^2 | X=x] f(x) R(\mathbf{K}) + o(1) \right) \nonumber \\
&\quad + p(x)^2 \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \left( f(x) R(\mathbf{K}) + o(1) \right) - 2p(x) \sum_{k=1}^{n} Q_k^{-2} \beta_k^2 h_k^{-d} \left( p(x) f(x) R(\mathbf{K}) + o(1) \right) \nonumber \\
&= \frac{\beta_n}{h_n^d} \psi_n^{-1} \frac{f(x)^2}{Q_n^2} [\Sigma + o(1)]. \tag{4.24}
\end{align}
$$

In addition, we have

$$
\forall p > 0, \quad \mathbb{E}[|L_k(x)|^{2+p}] = O \left( \frac{1}{h_k^{d(1+p)}} \right).
$$

Therefore,

$$
\begin{align*}
\mathbb{E}[|S_k(x)|^{2+p}] &= \mathbb{E}[|L_k(x) - \mathbb{E}[L_k(x)]|^{2+p}] \\
&\le 2Q_k^{-2-p}\beta_k^{2+p}\mathbb{E}[|L_k(x)|^{2+p}].
\end{align*}
$$

Hence,

$$
\mathbb{E} [|S_k(x)|^{2+p}] = O \left( Q_k^{-2-p} \beta_k^{2+p} \mathbb{E} [|L_k(x)|^{2+p}] \right).
$$

We then deduce that

$$
\begin{align*}
\sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] &= O\left(\sum_{k=1}^{n} Q_k^{-2-p} \beta_k^{2+p} \mathbb{E}[|L_k(x)|^{2+p}]\right) \\
&= O\left(\sum_{k=1}^{n} Q_k^{-2-p} \beta_k^{2+p} h_k^{-d(1+p)}\right).
\end{align*}
$$

In the following, let us assume that there is $p > 0$, such that

$$
\lim_{n \to +\infty} n \beta_n > \frac{1+p}{2+p} (\beta - ad).
$$
By applying lemma 1.2, we obtain

$$ \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = O\left(\frac{\beta_n^{1+p}}{Q_n^{2+p} h_k^{d(1+p)}}\right). $$

Thus,

$$ \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = O\left(\frac{\beta_n^{1+p}}{v_n^{2+p} Q_n^{2+p} h_n^{d(1+p)}}\right). $$

Then, it follows that

$$ \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = O\left(\left(\frac{\beta_n}{h_n^d}\right)^{p/2}\right) = o(1). $$

Moreover, since we have

$$ \lim_{n \to +\infty} \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}\left[|S_k(x) - \mathbb{E}[S_k(x)]|^{2+p}\right] = \lim_{n \to +\infty} \frac{1}{v_n^{2+p}} \sum_{k=1}^{n} \mathbb{E}[|S_k(x)|^{2+p}] = 0, $$

by applying the Lyapunov theorem, we get

$$ \frac{1}{\sqrt{v_n^2}} \sum_{k=1}^{n} (S_k(x) - \mathbb{E}[S_k(x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1), $$

which implies

$$ \frac{1}{v_n} \sum_{k=1}^{n} S_k(x) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1). $$

Moreover, (4.23) ensures that

$$ f(x) Q_n^{-1} v_n^{-1} (p_n(x) - \mathbb{E}[p_n(x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, 1). \quad (4.25) $$

Then, the combination of (4.24) and (4.25) ensures that

$$ \sqrt{\beta_n^{-1} h_n^d \psi_n} (p_n(x) - \mathbb{E}[p_n(x)]) \xrightarrow{n \to +\infty} \mathcal{N}(0, \Sigma). \quad (4.26) $$

Hence, the application of Lyapounov's Theorem coupled with the combination of (4.9), (4.10) and (4.26) ensures the convergence in (4.13). □
# Conclusion and perspectives

In this research work, we elaborated a multivariate recursive functional estimators. The basic merit of recursive estimators resides in the fact that one can update the estimation with each additional new observation. Therefore, instead of re-running the data each time, it is possible to rewrite our considered estimator as a combination of two (or more) estimators, where each estimator is based on separate dataset, which can be very interesting to keep the computational cost reasonably low. For this purpose, we investigated the stochastic approximation method in order to set forward four basic recursive estimators, namely the multivariate distribution function estimator, the multivariate conditional cumulative distribution function estimator, the multivariate kernel-type regression estimator and the propensity score regression estimator. We first tackled the asymptotic properties of each proposed estimator by providing the bias as well as the variance in order to demonstrate that our estimator asymptotically follows a normal distribution.

Subsequently, we revealed that, using a specific bandwidth selection, namely the cross-validation method, the plug-in procedure as well as the bootstrap technique, and a particular choice of the stepsize; for all the cases, each proposed recursive estimator yielded better results compared to the corresponding non-recursive one in terms of estimation error. The simulation studies and real datasets illustrate our findings. Following our numerical applications, we can demonstrate that our multivariate functional recursive estimators are closer to the true functional than the multivariate non-recursive ones.

At this stage of reflection, it would be appropriate to assert that the present thesis would be worthwhile in terms of opening up further promising directions of investigation and for providing valuable perspectives for future researches. First of all, we shall state some ongoing research works in progress following the same spirit of multivariate functional estimation in order to construct some new classes of recursive estimators.

## The conditional hazard function

The conditional hazard function $H(y|x)$ is defined by, for all real $y$ and $x$,

$$H(y|x) := \frac{f(y|x)}{S(y|x)}$$

where

* $f(y|x)$ denotes the conditional density of $Y$ given $X = x$.

* $S(y|x)$ denotes the conditional survival function defined by, for all real $y$ and $x$, $S(y|x) := 1 - \pi(y|x)$.

Then, an estimator of the hazard function is expressed by

$$H_n(y|x) = \frac{\hat{f}_n(y|x)}{S_n(y|x)},$$
with

$$ \hat{f}_n(y|x) = \frac{\prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k h_k^{-2d} \mathbf{K}_1 \left( \frac{y-X_k}{h_k} \right) \mathbf{K} \left( \frac{x-X_k}{h_k} \right)}{\prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k h_k^d \mathbf{K} \left( \frac{x-X_k}{h_k} \right)}, $$

and

$$ S_n(y|x) = 1 - \pi_n(y|x) = \frac{\prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k (1 - \chi_k(y)) h_k^{-d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right)}{\prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k h_k^{-d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right)}. $$

## The conditional variance

An estimator of the conditional variance $Var(Y|X)$ is determined by

$$ \sigma_n^2(x) = \sigma_{\varphi_n}^2(x) := \frac{\prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k (\varphi(Y_k) - r_{\varphi_n}(x))^2 h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right)}{\prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k h_k^{-d} \mathbf{K} \left( \frac{x-X_k}{h_k} \right)}. $$

It is noteworthy that the empirical estimator of the conditional variance is specified by

$$ \sigma_n^2(x) = \frac{1}{n} \sum_{i=1}^{n} (\varphi(Y_i) - r_n(X_i))^2. $$

An other perspective approach concerning the estimation of stochastic process is given in the following.

## Distribution function under Markov Renewal Process

A Markov Renewal Process (MRP) $(J_n, S_n)_{n \ge 0}$ is a random process that generalizes the notion of Markov jump processes, where $(J_n)_{n \ge 0}$ is a Markov chain and the process $(S_n)_{n \ge 0}$ are the jump times. Here, we may consider the distribution function associated with the sojourn time in state $i$ before going to state $j$,

$$ F_{ij}(x) = \mathbb{P}(X_{n+1} \le x | J_n = i, J_{n+1} = j), $$

where, for $n \ge 0$, $J_0, J_1, \dots, J_n$ are the consecutive states to be visited by the (MRP) and $(X_n)_{n \ge 0}$ defined by $X_0 = S_0 = 0$ and $X_n = S_n - S_{n-1}$ for $n \ge 1$, are the sojourn times in these states.

Then, an estimator of Jumps process with two conditions is identified by

$$ \hat{F}_n(x, y, t) = \frac{a_n(x, y, t)}{f_n(x, y)} $$

with

$$ a_n(x, y, t) = \prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k \chi_k(y) h_k^{-2d} \mathbf{K} \left( \frac{y-X_k}{h_k} \right) \mathbf{K} \left( \frac{x-X_k}{h_k} \right). $$

and

$$ f_n(x, y) = \prod_n \sum_{k=1}^n \Pi_k^{-1} \gamma_k h_k^{-2d} \mathbf{K} \left( \frac{x - X_k}{h_k} \right) \mathbf{K} \left( \frac{y - X_k}{h_k} \right). $$
It's noteworthy that, even if our proposed recursive estimators outperformed the non-recursive ones, we still have the possibility to ameliorate their performances in order to minimize the risk and guaranty a much faster convergence. For this purpose, we will attempt to propose other efficient estimators using new techniques of risk optimization.

Further future research directions would be to extend our findings to the following topics.

1. Extend our findings to the setting of serially dependent observations, α-mixing framework like in Khardani and Slaoui (2019).

2. Extend our findings to the case of functional data like in Slaoui (2019) and Slaoui (2020) as well as spatial data like in Bouzebda and Slaoui (2020).

3. Consider the k-nearest neighbours smoothing with functional regressor, see Almanjahie et al. (2020) in finite dimensional data and Kara et al. (2017) in the case of functional data.

4. Explore the idea developed in the recent work of Boukabour and Masmoudi (2020) through considering some semi-parametric Bayesian networks approaches based on the idea developed in our last chapter.
Bibliography

Abdelkader, B. (2014). Asymptotic normality of the recursive kernel estimate of conditional cumulative distribution function.

Altman, N. and Leger, C. (1995). Bandwidth Selection for Kernel Distribution Function Estimation. *J. Statist. Plann. Inference, **46***, 195-214.

Al-Awadhi, F., Kaid, Z., Laksaci, A., Ouassou, I. and Rachdi, M. (2019). Functional data analysis: local linear estimation of the $L_1$-conditional quantiles. *Stat. Methods Appl.*, **28**, 217-240.

Almanjahie, I M., Chikr Elmezouar, Z., Laksaci, A. and Rachdi, M. (2018). kNN local linear estimation of the conditional cumulative distribution function: Dependent functional data case. *C. R. Math. Acad. Sci. Paris*, **356**, 1036-1039.

Almanjahie, I. M, Aissiri, K. A., Laksaci, A. and Chikr Elmezouar, Z. (2020). The k nearest neighbors smoothing of the relative-error regression with functional regressor *Comm. Statist. Theory Methods*.

Amiri A. (2010). Estimateurs fonctionnels récursifs et leurs applications à la prévision. *Thèse en Mathématiques générales [math.GM]. Université d'Avignon*.

Benziadi, F., Laksaci, A. and Tebboune, F. (2016). Recursive kernel estimate of the conditional quantile for functional ergodic data. *Comm. Statist. Theory Methods, **45***, 3097-3113.

Berlinet, A., Gannoun, A. and Matzner-Lø ber, E. (1998). Normalité asymptotique d'estimateurs convergents du mode conditionnel. *Canad. J. Statist.*, **26**, 365-380.

Berlinet, A., Gannoun, A. and Matzner-Lø ber, E. (1998). Propriétés asymptotiques d'estimateurs convergents des quantiles conditionnels. *C. R. Acad. Sci. Paris Sér. I Math.*, **326**, 611-614.

Blum, J.R. (1954) Multidimensional stochastic approximation methods. *Ann. Math. Statist.*, **25**, 737-744.

Bonin, P., Chalard, M., Meot, A. and Fayol, M. (2002). The determinants of spoken and written picture naming latencies. *British Journal of Psychology, **93***, 89-114. *JPSS J. Probab. Stat. Sci.*, **12**, 117-126.

Bouanani, O., Rahmani, S. and Ait-Hennani, L. (2020). Local linear conditional cumulative distribution function with mixing data. *Arab. J. Math. (Springer)*, **9**, 289-307.

Boukabour, S. and Masmoudi, A. (2020). Semiparametric Bayesian networks for continuous data. *Comm. Statist. Theory Methods* doi:10.1080/03610926.2020.1738486.

Bouzebda, S. and Slaoui, Y. (2020). *Nonparametric Recursive Method for Kernel-Type Function Estimators for Censored Data, J. Stoch. Anal.* **3**, 19.
Brunel, E., Comte, F. and Lacour, C. (2010). Minimax estimation of the conditional cumulative distribution function. *Sankhya A*, **72**, 293–330.

Chikr-Elmezouar, Z., Almanjahie, I M., Laksaci, A. and Rachdi, M. (2019). FDA: strong consistency of the kNN local linear estimation of the functional conditional density and mode. *J. Nonparametr. Stat.*, **31**, 175–195.

Chilinski, P. and Silva, R. (2020). Neural likelihoods via cumulative distribution function. arXiv:1811.00974v2[stat.ML].

Chowdhury, M., Wu, C. and Modarres, R. (2018). Nonparametric estimation of conditional distribution functions with longitudinal data and time-varying parametric models. *Metrika*, **81**, 61–83.

Cousineau, D. and Chartier, S. (2010). Outliers detection and treatment: a review. *International Journal of Psychological Research* **3**, 58–67.

Deheuvels, P. and Mason, D. M. (2004). *General asymptotic confidence bands based on kernel-type function estimators*, Stat. Inference Stoch. Process. **7**, 225–277.

Deheuvels, P. (1973). *Sur l'estimation séquentielle de la densité*. C. R. Acad. Sci. Paris Sér. A-B **276**, A1119–A1121.

Delaigle, A and Gijbels, I. (2004). Practical bandwidth selection in deconvolution kernel density estimation, *Comput. Statist. Data Anal.* **45**, 249–267.

Dippon, J. and Renz, J. (1997). Weighted means in stochastic approximation of minima. *SIAM J. Control Optim.*, **35**, 1811–1827.

Dippon, J. (2003). Accelerated randomized stochastic optimization. *Ann. Statist.*, **31**, 1260–1281.

Duflo, M. (1997). *Random iterative models*, volume 34 of *Applications of Mathematics (New York)*. Springer-Verlag, Berlin. Translated from the 1990 French original by Stephen S. Wilson and revised by the author.

Einmahl, U. and Mason, D. M. (2000). *An empirical process approach to the uniform consistency of kernel-type function estimators*, *J. Theoret. Probab.* **13**, 1–37.

Fan, J., Yao, Q. and Tong, H. (1996). Estimation of conditional densities and sensitivity measures in nonlinear dynamical systems. *Biometrika*, **83**, 189–206.

Ferraty, F. and Vieu, P. (2002). The functional nonparametric model and application to spectrometric data. *Comput. Statist.*, **17**, 545–564.

Ferrigno, S. Foliguet, B. Myriam, M. and Muller-Gueudin, A. (2014). Certainty bands for the conditional cumulative distribution function and applications

Galambos, J. and Seneta E. (1973). Regularly Varying Sequences. *Proc. Amer. Math. Soc.* **41**, 110–116.

Gannoun, A., Girard, S., Guinot, C. and Saracco, J. (2002). Reference curves based on non-parametric quantile regression. *Statistics in Medicine*, **21**, 3119–3135.

Hall, P., Wolff, R C. L. and Yao, Q. (1999). Methods for estimating a conditional distribution function. *J. Amer. Statist. Assoc.*, **94**, 154–163.
Hardle, W. and Marron, J.S. (1985). Optimal Bandwidth Selection in Nonparametric Regression Function Estimation. *The Ann. Statist.* **13**, 1465–1481.

Hardle, W. and Marron, J. S. (1991). *Bootstrap Simultaneous Error Bars for Nonparametric Regression*, *The Ann. Statist.* **19**, 778–796.

Hill, P.D. (1985). *Kernel Estimation of a Distribution Function*. *Comm. Statist. Theory Methods*, **14**, 605–620.

Honda, T. (2000). *Nonparametric estimation of a conditional quantile for $\alpha$-mixing processes*. *Ann. Inst. Statist. Math.*, **52**, 459–470.

Honda, T. (2000). *Nonparametric estimation of the conditional median function for long-range dependent processes*. *J. Japan Statist. Soc.*, **30**, 129–142.

Jmaei, A. Slaoui Y. and Dellagi, W. (2017). Recursive distribution estimators defined by stochastic approximation method using Bernstein polynomials *J. Nonparametr. Stat.*, **29**, 792–805.

Kara, L.Z., Laksaci, A., Rachdi, M. and Vieu, P. (2017) Data-driven kNN estimation in non-parametric functional data analysis, *J. Multivariate Anal.* **153**, 176–188.

S. Khardani and Y. Slaoui. Recursive kernel density estimation and optimal bandwidth selection under alpha-mixing data, *J. Stat. Theory Pract.*, 13(36), 2019.

Kiefer, J. and Wolfowitz, J. (1952). Stochastic Estimation of the Maximum of a Regression Function. *Ann. Math. Statist.*, **23**, 462–466.

Kiwitt, S. and Neumeyer, N. (2012). Estimating the conditional error distribution in non-parametric regression. *Scand. J. Stat.*, **39**, 259–281.

Kushner, H.J. and Clark, D.S. (1978). Stochastic approximation methods for constrained and unconstrained systems. Springer, New York.

Kara, L.Z., Laksaci, A., Rachdi, M. and Vieu, P. (2017) Data-driven kNN estimation in non-parametric functional data analysis, *J. Multivariate Anal.* **153**, 176–188.

Laksaci, A. and Hachemi, N. (2012). Note on the functional linear estimate of conditional cumulative distribution function. *JPSS J. Probab. Stat. Sci.*, **10**, 153–160.

Laksaci, A. and Maref, F. (2009). Conditional cumulative distribution estimation and its applications. *JPSS J. Probab. Stat. Sci.*, **7**, 57–69.

Li, Q., Lin, J. and Racine, J S. (2013). Optimal bandwidth selection for nonparametric conditional distribution and quantile functions. *J. Bus. Econom. Statist.*, **31**, 57–65.

Luce, R.D. (1986). Response times: their role in inferring elementary mental organization. Handbook of Mathematical Psychology New York: Oxford.

McCormack, P.F. and Wright, N.M. (1964). The positive skew observed in reaction time distributions. *Canadian Journal of Psychology* **18**, 43–51.

McGill, W.J. (1963). Stochastic latency mechanisms. *Handbook of Mathematical Psychology* 1, 309–360. New York: John Wiley and Sons, Inc.

Milet, J. AND Nuel, G. AND Watier, L. AND Courtin, D. AND Slaoui, Y. AND Senghor, P. AND Migot-Nabias, F. AND Gaye, O. AND Garcia, A. *Genome Wide Linkage Study, Using a 250K SNP Map, of Plasmodium falciparum Infection and Mild Malaria Attack in a Senegalese Population*, PLOS ONE 5 (2010), 1–11.
Mokkadem, A. and Pelletier, M. (2007). *A companion for the Kiefer-Wolfowitz-Blum stochastic approximation algorithm*, Ann.Statist. **35**,1749–1772.

Mokkadem, A. and Pelletier, M. (2007). *Compact law of the iterated logarithm for matrix-normalized sums of random vectors*, Teor. Veroyatn. Primen. **52**, 752–767.

Mokkadem, A. Pelletier, M. and Slaoui, Y. (2009a). The stochastic approximation method for the estimation of a multivariate probability density. *J. Statist. Plann. Inference.* **139**, 2459–2478.

Mokkadem, A. Pelletier, M. and Slaoui, Y. (2009b). Revisiting Revesz's Stochastic Approximation Method for the Estimation of a Regression Function. *ALEA Lat. Amer.J. Probab.Math. Statist.* **6**, 63–114.

Mokkadem, A. and Pelletier, M. (2016). The Multivariate Revesz's Online Estimator of a Regression Function and Its Averaging. *M. Math. Meth. Stat.* **25**, 151–167.

Nadaraya, E.A. (1964). Some New Estimates for Distribution Functions. *Theory Probab. Appl.*, **9**, 497–500.

Ould-Saïd, E. (2006). A strong uniform convergence rate of kernel conditional quantile estimator under random censorship. *Statist. Probab. Lett.*, **76**, 579–586.

Parzen, E. (1962). On Estimation of a Probability Density Function and Mode. *Ann. Math. Statist.* **33**, 1065–1076.

Perret, C. and Bonin, P. (2019). Which variables should be controlled for to investigate picture naming in adults? A Bayesian meta-analysis. *Behavior Research Methods* **51**, 2533–2545.

Perret, C. and Laganaro, M. (2013). Why are written naming latencies (not) longer than spoken naming? *Reading and Writing An Interdisciplinary Journal* **26**, 225–239.

Perret, C. and Olive, T. (2019). Spelling and Writing Words: Theoretical and Methodological Advances. *Brills Edition*.

Perret, C., Bonin, P., and Laganaro, M. (2014). Exploring the multiple-level hypothesis of AoA effects in spoken and written picture naming using a topographic ERP analysis. *Brain and Language* **135**, 20–31.

Sandra Plancade. (2013) Adaptive estimation of the conditional cumulative distribution function from current status data. *J. Statist. Plann. Inference*, **143**, 1466–1485.

Polyak, B.T. and Tsybakov, A.B. (1990). Optimal orders of accuracy for search algorithms of stochastic optimization. *Problems Inform. Transmission.*, **26**, 126–133.

Reiss, R.D. (1981). Non parametric Estimation of Smooth Distribution Function. *Scand. J. Statist.*, **8**, 116–119.

Révész, P. (1973). Robbins-Monro procedure in a Hilbert space and its application in the theory of learning processes. *I. Studia Sci. Math. Hung.*, **8**, 391–398.

Révész, P. (1977). How to Apply the Method of Stochastic Approximation in the Non parametric Estimation of a Regression Function. *Math. Operationsforsch. Statist. Ser. Statistics.*, **8**, 119–126.

Robbins, H. and Monro, S. (1951). A Stochastic Approximation Method. *Ann. Statist.*, **22**, 400–407.
Rosenbaum, P. R. and Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. *Biometrika* **70**, 41–55.

Rosenblatt, M. (1956). Remarks on Some Nonparametric Estimates of a Density Function. *Ann. Math. Statist.*, **27**, 832–837.

Roussas, G. (1990). *Nonparametric regression estimation under mixing conditions*, *Stochastic Process. Appl.* **36**, 107–116.

Ruppert, D. (1982). Almost Sure Approximations to the Robbins–Monro and Kiefer–Wolfowitz Processes with Dependent Noise. *Ann. Probab.*, **10**, 178–187.

Sarda, P. (1993). Smoothing parameter selection for smooth distribution functions. *J. Statist Plann. Inference*, **35**, 65–75.

Schmiedek, F., Oberauer, K., Wilhelm, O., Süß, H.-M. and Wittmann, W.W. (2007). Individual differences in components of reaction time distributions and their relations to working memory and intelligence. *Journal of Experimental Psychology: General*, **136**, 414–429.

Slama, S. and Slaoui, Y. (2021). Multivariate distribution function estimation using stochastic approximation method, *International Journal of Mathematics and Statistics*, vol. 22, no. 2, pp. 31–59.

Slama, S. Slaoui, Y. and Fathallah, H. (2022). *Statistical inference for multivariate conditional cumulative distribution function estimation by stochastic approximation method*. *Statistics, Optimization & Information Computing*. 10(3), 789-814, https://doi.org/10.19139/soic-2310-5070-1416.

Slama, S., Slaoui, Y. Le Du, G. and Perret, C. (2022). Non parametric multivariate kernel regression estimation to describe cognitive processes and mental representations. *Statistics, Optimization & Information Computing*. 10(4), 1021-1043. https://doi.org/10.19139/soic-2310-5070-1507.

Slama, S. Slaoui, Y. and Fathallah, H. (2022). *The stochastic approximation method for semi-recursive multivariate kernel-type regression estimation*. *Theory of Stochastic Processes*.

Slaoui, Y. (2006). Application des méthodes d'approximations stochastiques à l'estimation de la densité et de la régression. Mathématiques [math]. Université de Versailles-Saint Quentin en Yvelines, HAL Id : tel-00131964, version 1.

Slaoui, Y. (2013). Large and Moderate Principles for Recursive Kernel Density Estimators Defined by Stochastic Approximation Method. *Serdica Math. J.*, **39**, 53–82.

Slaoui, Y. (2014a). Bandwidth selection for recursive kernel density estimators defined by stochastic approximation method. *Journal of Probability and Statistics*, ID 739640, doi:10.1155/2014/739640.

Slaoui, Y. (2014b). The stochastic approximation method for the estimation of a distribution function. *Math. Methods Statist.* **23**, 306–325.

Slaoui, Y. (2015). Plug-in bandwidth selector for recursive kernel regression estimators defined by stochastic approximation method. *Statistica Neerlandica*. **69**, 483–509.

Slaoui, Y. (2016). Optimal bandwidth selection for semi-recursive kernel regression estimators. *Statistics and Its Interface*. **9**, 375–388.

Slaoui, Y. (2017). Recursive kernel density estimators under missing data. *Comm. Statist. Theory Methods*. **18**, 9101–9125.
Slaoui, Y. (2018). Bias reduction in kernel density estimation. *J. Nonparametr. Stat.* **30**, 505–522.

Slaoui, Y. (2019). Wild Bootstrap Bandwidth Selection of Recursive Nonparametric Relative Regression for Independent Functional Data, *J. Multivariate Anal.*, **173**, 494–511.

Slaoui, Y. (2019). Large and moderate deviation principles for nonparametric recursive kernel distribution estimators defined by stochastic approximation method. *Opuscula Mathematica*, **39**, 733–746.

Slaoui, Y. and Jmaei, A. (2019). Recursive density estimators based on Robbins-Monro's scheme and using Bernstein polynomials. *Stat. Interface* **12**, 439–455.

Slaoui, Y. and Khardani, S. (2020). Adaptive recursive kernel conditional density estimators under censoring data. *ALEA Lat. Am. J. Probab. Math. Stat.* **17**, 389–417.

Slaoui, Y. (2020). Recursive non-parametric regression estimation for independent functional data, *Statist. Sinica*, **30**, 417–437.

Stute, W. (1986). Conditional Empirical Processes *Ann. Statist.*, **14**, 638–647.

Tsybakov, A.B. (1990). Recurrent Estimation of the Mode of a Multidimensional Distribution. *Probl. Inf. Transm.*, **8**, 119–126.

Unsworth, N., Redick, T.S., Lakey, C.E., and Young, D.L. (2010). Lapses in sustained attention and their relation to executive control and fluid abilities: An individual differences investigation. *Intelligence*, **38**, 111–122.

Veraverbeke, N., Gijbels, I. and Omelka, M. (2014) Preadjusted non-parametric estimation of a conditional distribution function. *J. R. Stat. Soc. Ser. B. Stat. Methodol.*, **76**, 399–438.

Watson, G.S. (1964). Smooth regression analysis. *Sankhya Ser. A.* **26**, 359–372.

Wegman, E.J. and Davies, H.I. (1979). Remarks on some recursive estimators of a probability density. *Ann. Statist.* **7**, 316–327.

Wolverton, C.T. and Wagner, T.J. (1969). Asymptotically optimal discriminant functions for pattern classification. *IEEE Trans. Inform. Theory.* **15**, 258–265.

Yu, K. and Jones, M. C. (1998) Local linear quantile regression. *J. Amer. Statist. Assoc.*, **93**, 228–237.