Pierre Ailliot

Julie Bessac

Valérie Monbet

Françoise Pène

Non-homogeneous hidden Markov-switching models for wind time series

Keywords: Wind time series, Markov-switching autoregressive process, non-homogeneous hidden Markov process, linear-circular time series, consistency

In this paper we propose various Markov-switching autoregressive models for bivariate time series which describe wind conditions at a single location. The main originality of the proposed models is that the hidden Markov chain is not homogeneous, its evolution depending on the past wind conditions. It is shown that they have good theoretical properties and permit to reproduce complex features of wind time series such as non-linear dynamics and the multimodal marginal distributions.

Introduction

Meteorological time series are a key input in many risk forecasting and impact studies applications and historical data are often available over periods of time that are not long enough to get reliable estimates of the quantities of interest. Stochastic weather generators have been developed to overcome this insufficiency by simulating artificial sequences of unlimited number of meteorological variables with statistical properties similar to those of the observations (see [START_REF] Srikanthan | Stochastic generation of annual, monthly and daily climate data: A review[END_REF] and references therein). These generators can also be useful for downscaling global climate models (see e.g. [START_REF] Maraun | Precipitation downscaling under climate change: Recent developments to bridge the gap between dynamical models and the end user[END_REF] and references therein) or infilling missing values by conditional simulation. In this paper we focus on wind time series. Various approaches have been proposed in the literature for modeling time series of wind speed (see e.g. [START_REF] Monbet | Survey of stochastic models for wind and sea state time series[END_REF] for a review). In comparison, there exists only very few models for circular time series of wind direction or for bivariate time series describing simultaneously the evolution of the wind speed and the wind direction. There is thus a need for models which can reproduce the specificities of such time series and this paper aims at filling this gap.

In this work we consider a wind time series for a location off the French Brittany coast which bivariate marginal distribution has complex features (see Figure 1). In particular it clearly exhibits two modes, each one corresponding to a different meteorological regime or 'weather type': the prevailing mode (westerly winds) corresponds to cyclonic conditions, i.e. low pressure systems (e.g. storms) coming from the North-Atlantic ocean, whereas the second mode (easterly winds) corresponds to anticyclonic conditions which can temporarily deviate or block the westerly flow. The alternation of such weather regimes is a well-known characteristic of the North-Atlantic/European area, and after [START_REF] Vautard | Multiple weather regimes over the north atlantic: Analysis of precursors and successors[END_REF] brought some evidence for quasi-stationary solutions in the equations of the atmospheric flow, thus giving a physical meaning to the statistically-derived regimes, they have been broadly used in climate studies. More generally, the presence of regimes with distinct weather conditions is a usual feature of meteorological time series and a classical approach for modeling these meteorological regimes consists in introducing a hidden (or latent) variable. This idea goes back to [START_REF] Zucchini | A hidden Markov model for space-time precipitation[END_REF] where Hidden Markov Models (HMMs) were proposed for modeling the space-time evolution of daily rainfall. HMMs have also been proposed for modeling time series of wind direction in [START_REF] Zucchini | Hidden Markov Models for Time Series: An Introduction Using R[END_REF]. However HMMs assume that the successive observations are conditionally independent given the latent weather type and cannot the strong relationship which exists between the wind conditions at successive time steps for the dataset considered in this work. Markov-Switching AutoRegressive (MS-AR) have been proposed in this context to model time series of wind speed in [START_REF] Monbet | Survey of stochastic models for wind and sea state time series[END_REF][START_REF] Pinson | Regime-switching modelling of the fluctuations of offshore wind generation[END_REF][START_REF] Ailliot | Markov-switching autoregressive models for wind time series[END_REF]. MS-AR models extend both the usual HMMs, by adding dynamics in the regimes, and AR models, which are often used to model wind time series (see e.g. [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF]), by introducing regime switchings through a latent variable.

In HMMs or MS-AR models, the evolution of the weather type is independent of the past weather conditions. For our particular example, it would imply for example that the probability of switching from the cyclonic conditions to the anticyclonic conditions between time t and time t + 1 does not depend on the wind conditions observed at time t whereas we know that these switchings generally occur when the wind is blowing from the north and is very unlikely to occur when the wind is blowing from the south. One originality of the models proposed in this paper is that the evolution of the latent weather type depends on past wind direction leading to non-homogeneous MS-AR (NHMS-AR) models. We show that NHMS-AR models lead to a better description of important characteristics of the data considered in this work, such as multimodality and non-linear dynamics, compared to MS-AR models.

The wind condition at a single location at time t can be described using the polar coordinates {U t , Φ t }, where U t denotes the wind speed with values in R + and Φ t the wind direction with values in T = R/2πZ or the Cartesian coordinates {u t , v t } where u t and v t denote respectively the zonal and meridional components with values in R. The polar coordinates are generally used by meteorologists, probably because they are easier to interpret. However, from a statistical point of view, it is probably more straightforward to model the time series of Cartesian components since many models, such as Gaussian vector AR models, have been proposed for bivariate time series with values in R 2 whereas the process {U t , Φ t } is a linear-circular process with values in R + × T and very few models have been proposed for such variables. Both representations are considered in this paper and a discussion of their respective advantages is given.

The paper is organized as follows. NHMS-AR models are introduced in Section 2 with specific parameterizations proposed when considering Cartesian and polar coordinates. In Section 3, we briefly describe the EM algorithm which has been used to maximize the likelihood function and discuss the asymptotic properties of the maximum likelihood estimates (MLE). Then the performances of the models are discussed and compared in Section 4. The data used in this paper are also introduced at the beginning of this section. At last, we make a synthesis of the obtained results and we give some perspectives in Section 5.

(x t |x t-1 , y t-1 ) = P (X t = x t |X t-1 = x t-1 , Y t-1 = y t-1 ),
• the conditional distribution of Y t given the values of {Y Let us write p(.|x t-1 t-u , y t-1 t-u ) for the conditional p.d.f. of (X t , Y t ) given

(X t-1 t-u = x t-1 t-u , Y t-1 t-u = y t-1 t-u ). Hypothesis 1 implies that for u ≥ s p(x t , y t |x t-1 t-u , y t-1 t-u ) = p 1 (x t |x t-1 , y t-1 )p 2 (y t |x, y t-1 t-s ). (1)
The various conditional independence assumptions are summarized by the directed graph below for s = 1.

• • • → X t-1 → X t → X t+1 → • • • ↓ ր ↓ ր ↓ • • • → Y t-1 → Y t → Y t+1 → • • •
NHMS-AR models define a quite general family of models:

• If p 1 (x t |x t-1 , y t-1 ) does not depend on y t-1 , we retrieve the usual MS-AR models which include the HMMs as a particular case (s = 0).

• If M = 1, {Y t } is an autoregressive process of order s.

• If p 1,θ (x k |x k-1 , y k-1 k-s
) does not depend on x k-1 and is parametrized using indicator functions, we obtain the Threshold AutoRegressive (TAR) models which is another important family of models with regime-switching in the literature (see e.g. [START_REF] Tong | Non-Linear Time Series: A Dynamical System Approach[END_REF]).

The following sections propose specific parametric models for p 1 (see Section 2.2) and p 2 when using Cartesian coordinates (see 2.3) or polar coordinates (see Sections 2.4 and 2.5).

Non-homogeneous Markov model for the weather type

As mentioned earlier, we introduce the latent process {X t } to describe the weather type which evolution may depend on previous wind direction. For example, we expect that the probability of switching from the cyclonic to the anticyclonic conditions generally is higher when the wind is blowing from the north that when it is blowing from the South. Such features can be modeled through the transition kernel p 1 . Hereafter we assume that

p 1 (x t |x t-1 , φ t-1 ) ∝ q x t-1 ,xt f V M (φ t-1 ; λ x t-1 ,xt , ψ x t-1 ,xt ), (2)
where f V M (.; κ, φ) is the p.d.f. of the von Mises distribution, φ t-1 is the wind direction at time t -1, Q = (q x,x ′ ) x,x ′ ∈{1,...,M } is a stochastic matrix and, for x, x ′ ∈ {1, ..., M }, λ x,x ′ ≥ 0 and ψ x,x ′ ∈ T are unknown parameters. The von Mises distribution is a natural distribution for circular variables (see [START_REF] Mardia | Statistics of directional data[END_REF]) which p.d.f. with respect to the Lebesgue measure on T, is given by

∀z ∈ T, f γ (z) := f V M (z; κ, φ) = 1 2πI 0 (κ) exp (κ cos(z -φ)) = 1 2πI 0 (|γ|) e γe -iz , (3)
with γ := κe iφ a complex parameter. I 0 denotes the modified Bessel function of order 0 defined as

I 0 (κ) := 1 2π T exp(κ cos(z)) dz.
φ ∈ T corresponds to the circular mean of the distribution and κ ≥ 0 describes the concentration of the distribution around φ: when κ = 0 we get the uniform distribution whereas the larger is κ the more concentrated around φ the distribution is. This distribution is denoted by V M (γ) hereafter.

According to [START_REF] Ailliot | An autoregressive model with time-varying coefficients for wind fields[END_REF], the probability that the hidden Markov chain {X t } switches from x t-1 to x t will increase when the wind direction Φ t-1 is close to ψ x t-1 ,xt and λ x t-1 ,xt models the directional spreading in which this transition is likely to occur. When λ x,x ′ = 0 for all x, x ′ ∈ {1, ..., M } we obtain again the homogeneous MS-AR models. Observe that (2) can be rewritten

p 1 (x t |x t-1 , φ t-1 ) = q x t-1 ,xt exp λx t-1 ,xt e -iφ t-1 M x ′ =1 q x t-1 ,x ′ exp λx t-1 ,x ′ e -iφ t-1 , (4)
with λx,x ′ ∈ C (by taking λx,x ′ = λ x,x ′ e iψ x,x ′ ). With this expression, it can be easily seen that replacing ( λx,x ′ ) x,x ′ by ( λx,x ′ -a x ) x,x ′ (for any choice of (a x ) x ) does not change p 1 and thus that identifiability constraints are needed.

In order to reduce the number of unknown parameters we add the following constraints for the non-homogeneous models developed in the sequel

λx,x ′ = λx ′ (5)
for all x, x ′ ∈ {1, ..., M } such that x = x ′ with the identifiability constraint

M x ′ =1 λx ′ = 0 (6)
We have also fitted the model without the constraint (5) and found that the likelihood of these models is similar to the one of the models with the constraint (5) whereas they have a significantly larger number of parameters. Even when assuming (5), we found that the parameter (λ x ′ ) is sometimes hard to fit in practice and that fixing its values to e.g. the concentration parameter of the von-Mises distribution fitted to the time series of wind direction leads to satisfactory models. The results obtained with these alternative strategies are not further discussed below.

Modeling the Cartesian coordinates conditionally to the weather type

In this section we propose a model for the bivariate process {Y t } = {u t , v t } conditionally to the weather type {X t }. This process has values in R 2 and the most classical autoregressive model for such process is the linear Gaussian vector autoregressive (VAR) model of order s. With this model, if

X t = x t then Y t = A (xt) 0 + A (xt) 1 Y t-1 + ... + A (xt) s Y t-s + Σ (xt) 1 2 ǫ t (7)
where

A (x) 0 ∈ R 2 , A (x) l
∈ R 2×2 for l ∈ {1, ..., s} and x ∈ {1, ..., M }, Σ (x) ∈ R 2×2 are symmetric positive matrices for x ∈ {1, ..., M } and {ǫ t } is a bivariate white noise sequence.

VAR models have been proposed for wind fields in a space-time context in [START_REF] Ailliot | An autoregressive model with time-varying coefficients for wind fields[END_REF][START_REF] Wikle | Spatiotemporal hierarchical bayesian modeling tropical ocean surface winds[END_REF][START_REF] Fuentes | Modeling and predicting complex space-time structures and patterns of coastal wind fields[END_REF]. On our particular dataset, we found that this model was not appropriate to reproduce the 'hole' around the origin which can be seen on the joint distribution on Figure 1. This hole corresponds to a low probability of observing low wind speed. We can get around this issue by applying a power transformation as follows ũt = U α t cos(Φ t ) ṽt = U α t sin(Φ t ) and fit the MS-AR model to {ũ t , ṽt } instead of {u t , v t }. The value α = 1.5 was chosen experimentally to remove the 'hole' close to the origin in the original distribution. The model with homogeneous hidden Markov chain is denoted HMS-AR (u,v) and the non-homogeneous model, where p 1 is given by ( 2), is denoted NHMS-AR (u,v) .

Modeling the wind direction conditionally to the weather type

In this section we propose a model for the circular process {Y t } = {Φ t }. The inclusion of {U t } in this model is discussed in the next section. Several autoregressive models have been proposed in the literature for modeling directional time series (see [START_REF] Breckling | The Analysis of Directional Time Series: Applications to Wind Speed and Direction[END_REF], [START_REF] Fisher | Time series analysis of circular data[END_REF], [START_REF] Holzmann | Hidden markov models for circular and linear-circular time series[END_REF], [START_REF] Kato | A markov process for circular data[END_REF]). We have chosen to focus on the von Mises process initially introduced in [START_REF] Breckling | The Analysis of Directional Time Series: Applications to Wind Speed and Direction[END_REF] and assume that the conditional distribution

of Y t given X t = x t , Y t-1 t-s = y t-1 t-s is V M γ (xt) 0 + s ℓ=1 γ (xt) ℓ e iy t-ℓ with γ (x) ℓ = κ (x) ℓ e iφ (x)
ℓ ∈ C for x ∈ {1, ..., M } and ℓ ∈ {0, ..., s}. This can be rewritten

p 2 (y t |x t , y t-1 t-s ) = 1 b(x t , y t-1 t-s ) exp κ (xt) 0 cos(y t -φ (xt) 0 ) + s ℓ=1 κ (xt) ℓ cos(y t -y t-ℓ -φ (x) ℓ ) = 1 b(x t , y t-1 t-s ) exp [γ (xt) 0 + s ℓ=1 γ (xt) ℓ e iy t-ℓ ]e -iyt (8) with b(x t , y t-1 t-s ) := T exp κ (xt) 0 cos(y -φ (xt) 0 ) + s ℓ=1 κ (xt) ℓ cos(y -y t-ℓ ) dy = I 0 γ (xt) 0 + s ℓ=1 γ (xt) ℓ e iy t-ℓ .
In [START_REF] Breckling | The Analysis of Directional Time Series: Applications to Wind Speed and Direction[END_REF], it was assumed that γ

(x) ℓ ∈ R for ℓ ∈ {1, ..., s}.
We have chosen to extend it to a model with complex parameters in order to be able to reproduce the prevailing rotation of the wind direction in the clockwise direction (see Section 4).

The model with homogeneous hidden Markov chain is denoted HMS-EVM and the nonhomogeneous model, where p 1 is given by (2), is denoted NHMS-EVM.

A conditional model for the wind speed given the wind direction

In [START_REF] Ailliot | Markov-switching autoregressive models for wind time series[END_REF] it was proposed to model the wind speed {U t } using a homogeneous MS-AR model with M = 3 regimes and Gaussian linear AR models (see ( 7)) of order s = 2. Figure 2 shows typical examples of wind speed and wind direction time series together with the regimes identified by the fitted MS-AR models. These regimes basically correspond to periods with different temporal variability and there seems to be no simple relation between the regimes identified on the wind speed and the wind direction. In this context, it does not seem relevant to use the same weather type for the two time series. We thus propose to introduce a different weather type X (U ) t for the wind speed and X (Φ) t for the wind direction.

In order to explore the link between X (U ) t and Φ, we have computed the most likely values of X (U ) t

given the observed time series of wind speed {U t } and produced rose plots of the wind direction in the different weather types which were identified for the wind speed. We got plots very similar to the ones shown on Figure 3 (right panel). The first regime, which corresponds to periods with low temporal variability for the wind speed (anticyclonic conditions), can occur in any wind direction whereas the more variable regimes 2 and 3 are mainly associated to southwesterlies (cyclonic conditions). Stated otherwise the wind direction provides information on the synoptic weather conditions which control the intensity and the variability of the wind speed. It suggests the use of a non-homogeneous MS-AR model for the wind speed where the transition probabilities depend on the wind direction. Hereafter NHMS-AR (U,Φ) denotes the model for

{U t , Φ t } such that -{Φ t } is modeled by the NHMS-EVM model for {x (Φ) t , Φ t } with M = 4 and s = 2,
-{U t } is modeled conditionally to {Φ t } by a NHMS-AR model with p 1 given by ( 2) and p 2 by a linear Gaussian AR model [START_REF] Cappé | Inference in hidden Markov models[END_REF].

The structure of the model, with two layers of hidden variables, one for the wind speed and one for the wind direction, is shown on the directed graph below.

• • • / / X (Φ) t-1 / / X (Φ) t / / X (Φ) t+1 / / • • • • • • / / = = ! ! Φ t-1 / / < < " " Φ t / / < < " " Φ t+1 / / = = ! ! • • • • • • / / X (U ) t-1 / / X (U ) t / / X (U ) t+1 / / • • • • • • / / U t-1 / / U t / / U t+1 / / • • • 3 Parameter estimation 3.

Numerical computation of the MLE

The parameter vector of NH-MSAR models is composed of the parameters θ Q of the transition probabilities p 1 (x t |x t-1 , y t-1 ), the parameters θ (x) of the transition kernel p 2 (y t |x, y t-1 t-s ) for each regime x ∈ {1, ..., M }.

They are estimated by maximizing the likelihood function using a generalized EM algorithm. This algorithm was initially introduced in [START_REF] Baum | A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov chains[END_REF] for HMMs and then generalized to models with latent variables in [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. This recursive algorithm computes successive approximations θi of the maximum likelihood estimate (MLE) θ by cycling through the following steps.

E-step: Compute Q θ| θi = E θi (log(p θ (X T 1 , Y T 1 ))|y T -s+1
) as a function of θ. M-step: Determine the updated parameter estimate θi+1 = arg max θ Q θ| θi .

The conditional probabilities involved in the computation of Q(θ| θn ) are computed using the so-called forward-backward recursions (see e.g. [START_REF] Cappé | Inference in hidden Markov models[END_REF] and references therein). The particular implementation of these recursions for homogeneous MS-AR models is discussed in [START_REF] Hamilton | Analysis of time series subject to changes in regime[END_REF] and [START_REF] Hughes | A non-homogeneous hidden markov model for precipitation occurrence[END_REF] discusses it for non-homogeneous HMMs. It can be easily generalized to the models considered in this paper. The M-step requires numerical optimization leading to the so-called Generalized EM (GEM) algorithm. In order to get an efficient EM algorithm, it is important to implement carefully the numerical optimization procedure. In practice, the function Q(.| θi ) which has to be maximized in the M-step can be written as the sum of M + 1 functions as follows

Q θ| θi = Q X θ Q | θi + M x=1 Q Y θ (x) | θi .
This leads to solving M + 1 separate optimization problems on spaces with reduced dimension which is far more efficient than maximizing directly Q(.| θi ) over all the parameters. Note that analytical expressions are available for arg max

θ (x) Q Y θ (x)
| θi when Gaussian linear AR models [START_REF] Cappé | Inference in hidden Markov models[END_REF] are used and for arg max

θ Q Q X θ Q | θi when homogeneous MS-AR models are considered.
In order to avoid convergence to non-interesting maxima and save computational time, a proper initialization of this algorithm with realistic parameter values θ0 is needed. In practice, we have used the nested nature of the models. We have first fitted homogeneous models and then use the estimated parameters as a starting point for the corresponding non-homogeneous models. In the same spirit, the results obtained for the model of order s have been used to initialize the EM for the models of order s + 1.

Ergodicity of the models and asymptotic properties of the MLE

Conditions on the parameters which warrant the existence of a unique stationary ergodic solution with finite moments of order 2 are given in [START_REF] Ailliot | Consistency of the maximum likelihood estimate for nonhomogeneous markov-switching models[END_REF] for the models with linear Gaussian AR models [START_REF] Cappé | Inference in hidden Markov models[END_REF]. These are desirable properties for the fitted models since we expect that the wind process satisfies such properties and imply that the MLE are consistent (see [START_REF] Ailliot | Consistency of the maximum likelihood estimate for nonhomogeneous markov-switching models[END_REF]). These conditions apply to the models with linear Gaussian AR discussed in Section 4 and thus the MLE are consistent for these models. In this section we show that similar results hold true for the model with von Mises distribution introduced in Section 2.4 and prove more precisely ψ-irreducibility, aperiodicity, Harris-recurrence, identifiability and consistency for the NHMS-EVM model.

In the rest of this section, we assume Hypothesis 1 with p 1 and p 2 given by ( 4) and (8) respectively. Let Θ ′ be the set of parameters θ = (q x,x ′ , λx,x ′ , γ [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF]. Let us now state our main result.

(x) ℓ ) ℓ,x,x ′ such that γ (x) ℓ ∈ C, q x,x ′ > 0 such that x ′ q x,x ′ = 1 and λx,x ′ ∈ C satisfying
Theorem 2 (Consistency for NHMS-EVM). Assume that Θ is a compact subset of Θ ′ and that the coordinates of the true parameter θ * satisfy

x = x ′ ⇒ (γ (x) 0, * , ..., γ (x) s, * ) = (γ (x ′ ) 0, * , ..., γ (x ′ ) s, * ). (9)
Then, for every x 0 ∈ {1, ..., M } and any initial measure ν on {1, ..., M }×T, on a set of probability one, the limit values θ = (γ, Q, λ) of the sequence of MLE ( θn,x 0 ) n are equal to θ * = (γ * , Q * , λ * ) up to a permutation of indices, i.e. for any such limit value θ, there exists a permutation σ of {1, ..., M } such that, for every x, x ′ ∈ {1, ..., M }, for every j = 0, ..., s, the following relations hold true γ

(x) j = γ (σ(x)) j, *
, q x,x ′ = q σ(x),σ(x ′ ), * and λx,x ′ = λσ(x),σ(x ′ ), * .

One can notice that (9) just means that there is no couple of regimes (x, x ′ ) with x = x ′ in which the behavior of the process {Y t } is the same.

The proof of Theorem 2 is based on two ingredients: a general consistency result established in [START_REF] Ailliot | Consistency of the maximum likelihood estimate for nonhomogeneous markov-switching models[END_REF][Thm 2] and the proof of the identifiability up to a permutation of indices (see Proposition 4). Let us first check that the conditions of [START_REF] Ailliot | Consistency of the maximum likelihood estimate for nonhomogeneous markov-switching models[END_REF][Thm 2] apply for the NHMS-EVM model.

Since for every (θ, x, y) ∈ Θ × {1, ..., M } × T, q θ (x, y|•, •) is continuous on the compact set {1, ..., M } × T s , we have

α = E×K γ(x, y) dµ 0 (x, y) > 0, with γ(x, y) := inf x ′ ,y -1 -s q θ (x, y|x ′ , y -1 -s ).
Now we consider the probability density function (w.r.t. µ 0 ) β given by β(x, y) := γ(x, y) α .

For every x 0 , x -1 ∈ E and every y 0 -s , we have

q θ (x 0 , y 0 |x -1 , y -1 -s ) ≥ αβ(x 0 , y 0 ).
This implies the ψ-irreducibility, the strong aperiodicity (the ν 1 -small set being the whole space), the Harris recurrence (since we can decompose the whole set in a union of uniformly accessible sets from the whole set) and the positiveness (the invariant measure being unique and finite).

In particular, this gives Assumption (5

) of [3][Thm 2].
Moreover, since p 1,θ (x 1 |x 0 , y 0 ) and p 2,θ (y 0 |x 0 , y -1 ) are continuous in (θ, x 1 , x 0 , y 0 ) and in (θ, x 0 , y 0 , y -1 ) (respectively), all the other assumptions of [START_REF] Ailliot | Consistency of the maximum likelihood estimate for nonhomogeneous markov-switching models[END_REF][Thm 2] are satisfied for any compact subset of Θ ′ . Hence, we have Corollary 3. Assume that Θ is a compact subset of Θ ′ . Then, for all θ ∈ Θ, there exists a unique invariant probability and, for every x 0 ∈ E and every initial probability ν, the limit values of ( θn,x 0 ) n are Pθ * -almost surely contained in {θ ∈ Θ : Pθ = Pθ * }. Now, Theorem 2 will follow from the following result giving the identifiability of the parameter (up to permutation of indices) .

Proposition 4 (Identifiability). Let θ 1 and θ 2 belong to Θ ′ :

θ i = (γ (x) j,(i) ) j,x , (q x,x ′ ,(i) ) x,x ′ , ( λx,x ′ ,(i) ) x,x ′ .
Assume that the parameters (γ (x) j,(1) ) j,x which model the evolution of the wind direction in the different regimes through [START_REF] Cattiaux | European temperatures in cmip5: origins of present-day biases and future uncertainties[END_REF] for θ 1 are such that

x = x ′ ⇒ (γ (x) 0,(1) , ..., γ (x) s,(1) ) = (γ (x ′ ) 0,(1) , ..., γ (x ′ ) s,(1) ). ( 10
)
Then PY θ 1 = PY θ 2 if and only θ 1 is equal to θ 2 up to a permutation of indices.

Proof of Proposition 4. We write Leb for the Lebesgue measure on T. Assume that PY θ 1 = PY θ 2 . In particular, we have Since pθ 1 (y t t-s ) > 0 (since the invariant density h 1 satisfies h 1 > 0 since α > 0 and the transition density p satisfies p > 0 by construction) and since (4) holds true, we deduce that, for Leb ⊗(s+1)a.e. y t t-s , we have

pθ 1 (Y t = y t |Y t-1 t-s = y t-1 t-s ) = pθ 2 (Y t = y t |Y t-1 t-s = y t-1 t-s ),
M x=1 Pθ 1 (X t = x|y t-1 t-s )f γ (x) 0,(1) + s ℓ=1 γ (x) ℓ,(1) e iy t-ℓ (y t ) = M x=1 Pθ 2 (X t = x|y t-1 t-s )f γ (x) 0,(2) + s ℓ=1 γ (x)
ℓ,(2) e iy t-ℓ (y t ) with f γ defined by (3). Due to [START_REF] Fraser | Identifiability of finite mixtures of von mises distributions[END_REF], finite mixtures of von Mises distributions are identifiable. Hence if .

M x=1 π (x) 1 f γ (x) 1 (y) = M x=1 π (x) 2 f γ (x) 2 (
Recall that we have assumed that θ

(x) Y,1 = θ (x ′ ) Y,1 if x = x ′ , which implies that f or Leb ⊗s -a.e. y t-1 t-s , γ (x) 0,1 + s ℓ=1 γ (x)
ℓ,1 e iy t-ℓ = γ

(x ′ ) 0,1 + s ℓ=1 γ (x ′ )
ℓ,1 e iy t-ℓ .

Therefore, since for every x ∈ {1, ..., M } and for Leb ⊗s -almost every y t-1 t-s , Pθ 1 (X t = x|y t-1 t-s ) > 0 (since h θ 1 > 0), for Leb ⊗s -almost every y t-1 t-s there exists a permutation

σ y t-1 t-s of {1, ..., M } such that, ∀x ∈ {1, ..., M }, γ (x) 0,(1) + s ℓ=1 γ (x) ℓ,(1) e iy t-ℓ = γ (σ y t-1 t-s (x)) 0,(2) + s ℓ=1 γ (σ y t-1 t-s (x)) ℓ,(2) e iy t-ℓ .
Since the set of permutations of {1, ..., M } is finite, there exists a positive Lebesgue measure subset of T s on which the permutation is the same permutation σ. From this, we deduce that ∀x ∈ {1, ..., M }, ∀j ∈ {0, ..., s}, γ [START_REF] Ailliot | An autoregressive model with time-varying coefficients for wind fields[END_REF] and that, for Lebesgue almost every y t+1 t-s , the following holds true ∀x ∈ {1, ..., M }, Pθ 1 (X t = x|y t-1 t-s ) = Pθ 2 (X t = σ(x)|y t-1 t-s ).

(x) j,(1) = γ (σ(x)) j,
Let us now discuss the identifiability of the other components of θ 1 and θ 2 . If PY

θ 1 = PY θ 2 then pθ 1 (Y t = y t , Y t+1 = y t+1 |Y t-1 t-s = y t-1 t-s ) = pθ 2 (Y t = y t , Y t+1 = y t+1 |Y t-1 t-s = y t-1 t-s ) PY θ 1 a.s.
and thus, for Lebesgue almost every y t+1 t-s , we have

M x,x ′ =1 Pθ 1 (X t = x|y t-1 t-s )p 1,θ 1 (x ′ |x, y t )f γ (x) 0,(1) + s ℓ=1 γ (x)
ℓ,(1) e iy t-ℓ (y t )f γ (x ′ ) 0,(1) + s ℓ=1 γ

(x ′ ) ℓ,(1) e iy t-ℓ+1 (y t+1 ) = M x,x ′ =1 Pθ 2 (X t = x|y t-1 t-s )p 1,θ 2 (x ′ |x, y t )f γ (x) 0,(2) + s ℓ=1 γ (x) ℓ,(2) e iy t-ℓ (y t )f γ (x ′ ) 0,(2) + s ℓ=1 γ (x ′ )
ℓ,(2) e iy t-ℓ+1 (y t+1 ).

Using again the identifiability of von Mises distribution, we obtain ∀x, x ′ , p 1,θ 1 (x ′ |x, y t ) = p 1,θ 2 (σ(x ′ )|σ(x), y t ) for Leb -a.e. y t .

Now, due to the special form of p 1,θ specified in (4), we get ∀x, x ′ , Leb -a.e. y t , q x,x ′ ,(1) exp λx,x ′ ,(1) e -iyt

M

x ′′ =1 q x,x ′′ ,(1) exp λx,x ′′ ,(1) e -iyt = q σ(x),σ(x ′ ),(2) exp λσ(x),σ(x ′ ),(2) e -iyt

M

x ′′ =1 q σ(x),x ′′ ,(2) exp λσ(x),x ′′ ,(2) e -iyt .

(11) Let x ∈ {1, ..., M } be fixed. Applying (11) a first time with x ′ = x and a second time with any x ′ , we get ∀x ′ , for Leb -a.e. y t , q x,x ′ ,(1) exp λx,x ′ ,(1) e -iyt q x,x,(1) exp λx,x,(1) e -iyt = q σ(x),σ(x ′ ),(2) exp λσ(x),σ(x ′ ),(2) e -iyt q σ(x),σ(x),(2) exp λσ(x),σ(x),(2) e -iyt and so ∀x ′ , q x,x ′ ,(1) q x,x,(1) = q σ(x),σ(x ′ ),(2) q σ(x),σ(x),(2) [START_REF] Fraser | Identifiability of finite mixtures of von mises distributions[END_REF] and ∀x ′ , λx,x ′ ,(1) -λx,x,(1) = λσ(x),σ(x ′ ),(2) -λσ(x),σ(x),(2) .

Now, since x ′ q x,x ′ ,(1) = 1 = x ′ q σ(x),σ(x ′ ),(2) , due to [START_REF] Fraser | Identifiability of finite mixtures of von mises distributions[END_REF], it comes q x,x,(1) = q σ(x),σ(x),(2) and so ∀x ′ ∈ E, q x,x ′ ,(1) = q σ(x),σ(x ′ ), [START_REF] Ailliot | An autoregressive model with time-varying coefficients for wind fields[END_REF] .

Since θ 1 and θ 2 are in Θ ′ , x ′ λx,x ′ ,(1) = 0 = x ′ λσ(x),σ(x ′ ), [START_REF] Ailliot | An autoregressive model with time-varying coefficients for wind fields[END_REF] , and due to (13), we get λx,x,(1) = λσ(x),σ(x),( 2) and, applying again (13), we conclude that

∀x ′ ∈ E, λx,x ′ ,(1) = λσ(x),σ(x ′ ),(2) .
4 Numerical results

Data

In this paper, we focus on a wind time series extracted from the ERA-40 data set which consists in a global reanalysis with 6-hourly data covering the period from 1958 to 2001. This reanalysis was carried out by the European Centre for Medium-range Weather Forecast (ECMWF) and can be freely downloaded and used for scientific purposes at the URL: http://data.ecmwf.int/data. We consider the surface wind data for the point with geographical coordinates (47.5 0 N, 5 0 W) which is located off the Brittany coast (northwest of France). We have performed a comparison with in-situ data which indicates that this reanalysis data provides an accurate description of the wind condition observed at this location with the advantage of being easy to use in a statistical study (long time series with no missing data). The resulting time series is non-stationary since it exhibits an important seasonal component but also diurnal and inter-annual components. A classical approach for treating seasonality in meteorological time series consists in blocking the data, typically by period ranging from a month to a trimester depending on the amount of data available, and in fitting a separate model for each period in the year. This approach has been used in this paper and we have chosen to focus on the months of January. It leads to 44 time series of length 124 (31 days with 4 observations each day), each time series describing the wind conditions during the month of January for a particular year. In the sequel, we assume that these time series are independent realizations of a stationary process. It seems realistic according to the results given in [START_REF] Ailliot | Markov-switching autoregressive models for wind time series[END_REF] for the wind speed at the same location since the diurnal components can be neglected during the winter season. Following [START_REF] Ailliot | Markov-switching autoregressive models for wind time series[END_REF], another approach would consist in letting some of the coefficients of the models introduced hereafter to evolve in time with periodic functions for the diurnal and seasonal components and eventually a trend.

Model selection

Before analyzing any numerical result, one has to discuss the choice of the number of regimes and of the order of the AR models. In practice we found that the BIC criterion generally permits to identify parsimonious models which fit well the data. It is defined as

BIC = -2 log L + k log(N )
and L is the likelihood of the data, k is the number of parameters and N is the number of observations. In order to make the final selection among the best models identified by BIC, we have compared their abilities to generate realistic wind time series since this is the main motivation for this work. For this, a large number of realizations of the various models under competition have been simulated and various statistics of these synthetic sequences have been compared with the ones of the original data.

The models were fitted with a number of regimes M varying from 1 to 6 at the most and the BIC values suggest selecting models with M = 3 or M = 4 regimes (see Tables 1 and2). For the wind direction the model with M = 4 regimes tends to better reproduce the marginal distribution of the process compared to the models with M = 3 regimes and we thus have chosen to select this model. The NHMS-EVM model with M = 4 and s = 2, which is used in the NHMS-AR (U,Φ) model, has 43 parameters. For the Cartesian coordinates {u t , v t } and for the wind intensity {U t } the models with M = 3 and M = 4 regimes lead to similar results and we have thus chosen to keep the simplest model with M = 3.

We also varied the order s of the autoregressive models from s = 0 (y t is independent of y t-1 0 given x t ) to s = 5 and the BIC values are generally decreasing with s suggesting that a model of order s ≥ 5 may be needed. Notice however that there is generally a big improvement in the BIC values when s increases from 0 to 1 and from 1 to 2 whereas the difference is much smaller when comparing the models of order s = 2 and s ≥ 3 (not shown). We will focus on models of order s = 2 in the sequel. We believe that models of reduced order are more realistic from a physical point of view and we get similar simulation results with s = 2 compared to the models with s ≥ 3.

The BIC of {u t , v t } models are generally smaller than the ones of {U t , Φ t } models except for s = 0. It may be due to the higher number of parameters involved in the NHMS-AR (U,Φ) model (66 parameters when M = 3), which has two layers of hidden variables whereas the NHMS-AR (u,v) model has one common weather type for {u t } and {v t } and only 44 parameters when M = 3. Note however that BIC does not permit to make a clear distinction between both parameterizations (polar or Cartesian) since the differences in BIC values are relatively small. The bold values correspond to the selected models.

Regimes can be interpreted as weather types

An important benefit of using weather type models for meteorological variables is that they generally lead to interpretable models. This is illustrated in this section on NHMS-AR (u,v) and NHMS-AR (U,Φ) models. In order to compare the regimes of these two models, they have been ordered according to the variance of the innovation of the autoregressive processes Σ (s) . Figure 3 shows that the distributions of the wind direction in the different regimes are broadly similar for both models. The first regime corresponds mainly to anticyclonic conditions with easterly wind and a slow varying intensity (the variance of the innovation of the AR model is lower than in the other regimes and the first AR coefficient is larger). This regime is also the most likely (probability of occurrence of about 46%). The two other regimes correspond to cyclonic conditions with westerly wind and higher temporal variability in the intensity. These two regimes are discriminated mainly by the temporal variability, which is higher in the third regime, and the wind direction with the second regime corresponding mainly to south-westerlies and the third regime corresponding mainly to north-westerlies (see Figure 3). In order to assess the physical consistency of these three local regimes, we confronted them to the large-scale regimes commonly used in climate studies. More precisely, we considered the four wintertime regimes of [START_REF] Cattiaux | European temperatures in cmip5: origins of present-day biases and future uncertainties[END_REF], obtained over the North-Atlantic/European sector (90W-30E / 20-80N) by a kmeans-clustering of 3607 daily maps of 500 mb geopotential anomalies (days of December, January and February 1981-2010). The classification of each day of January 1979-2001 into these four types was used to compute the conditional probabilities given in Table 3. It shows a clear link between the regimes identified by the NHMS-AR (u,v) model and the largescale regimes. For example the first regime of the NHMS-AR (u,v) model, with low temporal variability, is more likely to occur when an anticyclone, generally located over Scandinavia, blocks the westerly flow (large scale regime denoted BL). At the opposite the more variable third regime of the NHMS-AR (u,v) model is generally associated with the large scale regime NAO+ (positive phase of the North Atlantic Oscillation). The regimes AR (Atlantic Ridge) and NAO-(negative phases of the North Atlantic Oscillation) have intermediate temporal variability with AR being more stable than NAO-. These results are coherent with the climatology of the area. 3: Joint probability of occurrence of the three regimes identified by the the NHMS-AR (u,v) model (lines) and the large-scale regimes provided by J. Cattiaux (see [START_REF] Cattiaux | European temperatures in cmip5: origins of present-day biases and future uncertainties[END_REF]) in columns.

Regime

Both NHMS-AR (u,v) and NHMS-AR (U,Φ) models have similar transition probabilities (see Figure 4) with a more pronounced dependence on the wind direction for the NHMS-AR (U,Φ) model. The more persistent regime is clearly the first one (mean duration of about 3.4 days) with a high probability of staying in this regime in any wind direction. The probability of switching directly from regime 1 to regime 3 is very small and the Markov chain will generally transit through the regime 2. Transitions from regime 1 to regime 2 are more likely when the wind is blowing from the west and transitions from regime 2 to regime 3 generally occur when the wind is from south. Regime 3 is persistent only when the wind is from south-west. If the wind blows from other directions, the weather type will quickly switch to regime 1 or 2.

Marginal distributions

According to Figure 5, the models with non-homogeneous transition probabilities provide a better description of the marginal distribution of the wind direction than homogeneous models which have difficulties in reproducing the second mode of the distribution (associated to easterlies). The NHMS-AR (U,Φ) model seems to perform slightly better than the NHMS-AR (u,v) model. The joint distribution of {u t , v t } is globally well reproduced by both non-homogeneous models (see Figure 6). Simulated data exhibit two modes like in the original data. The modes seem to be slightly better located with the NHMS-AR (U,Φ) model. It may be due to the small differences in the non-homogeneous transition probabilities shown in Figure 4 with the NHMS-AR (U,Φ) model having a slightly higher probability of staying in regime 3 when the wind is blowing from the south-west. It may help to create two distinct modes at the correct locations.

Both models generate too much low wind and as a consequence fail to reproduce accurately the hole at origin. Similar lack of fit for the marginal distribution was observed on other datasets with MS-AR models. It seems to be especially sensible when the model is miss-specified. This discrepancy may be reduced by developing alternative estimation methods which would give more importance to the stationary distribution of the process. This will be the topic of future research.

Dependence structure

All the models reproduce approximatively the first lags of autocorrelation function of {U t } (see Figure 7) and the circular autocorrelation of Φ (not shown) defined as (see [START_REF] Fisher | A correlation coefficient for circular data[END_REF])

ρ(h) = E[cos(Y 0 ) cos(Y h )] + E[sin(Y 0 ) sin(Y h )] -E[sin(Y 0 ) cos(Y h )] -E[cos(Y 0 ) sin(Y h )] E[cos(Y 0 ) 2 ]E[sin(Y 0 ) 2 ] -E[sin(Y 0 ) cos(Y 0 )] 2 (14)
for any positive integer h. To further validate the models, we have also plotted the various terms which appear in [START_REF] Hamilton | Analysis of time series subject to changes in regime[END_REF]. According to Figure 8, the autocorrelation function of {cos(Φ t )} is generally better reproduced than the one of {sin(Φ t )}. The empirical autocorrelation of {sin(Φ t )} has a more complex shape, with a quick decrease close to the origin and a bump around 4 days, than the one of {cos(Φ t )} which exhibits a more monotonic decrease. Figure 7 shows the cross-correlation function between the time series {cos(Y t )} and {sin(Y t )}. The sample crosscorrelation function is at its maximum value for a lag between 18 hours and 24 hours, with the time series {sin(Y t )} being in advance of the time series {cos(Y t )}. This may be related to the fact that, for the location of interest, the wind direction tends to rotate more often clockwise than anti-clockwise between two successive time steps (see Figure 9). Note that the complex parametrization of the von Mises autoregressive models (see Section 2.4) permits to model rotation in a prevailing direction and significantly improves the boxplot shown on Figure 9 compared to models with real parametrization (not shown). One can also remark that the first order autoregressive matrices of the NHMS-AR (u,v) model have diagonal coefficients which are close to each other and out-diagonal coefficients which are almost opposed and thus may be interpreted as the product of rotation and dilatation matrices. Figure 9 shows however that they do not generate enough anticlockwise rotations.

The non-homogeneous models generally lead to a better description of the correlation functions compared to the homogeneous models. All the models lead to an underestimation of the empirical autocorrelations functions of the time series {cos(Y t )} and {sin(Y t )}. Increasing the order s of the autoregressive models leads to a better description of the second order structure of the process but models of order s ≥ 3 cannot reproduce the second mode of the marginal distribution and thus models of order s = 2 seem to provide a good compromise.

Conclusion

In this paper we propose to model bivariate wind time series considering Cartesian coordinates on one hand and polar coordinates on the other hand. Both approaches have advantages. The {u t , v t } model is easier to write and to fit since it is based on Gaussian distributions whereas a linear-circular model is required when considering polar coordinates. The {u t , v t } model permits to globally well restore the second order structure observed on the data while the {U t , Φ t } model seems to give a better description of the marginal distributions. However, the differences between both models are slight.

Models with homogeneous and non-homogeneous latent Markov chains are compared. In nonhomogeneous models, the transitions depend on the wind direction at the previous time. At the location of interest, wind is rotating more often clockwise but wind direction may also oscillate around two prevailing directions (northeast for anti-cyclonic conditions and southwest for cyclonic conditions). These features induce respectively some cycles which can be seen in the second order structure and modes in the marginal distribution. In broad outline, nonhomogeneous transitions help the process to stay in the same weather regime when the wind direction is close to the prevailing directions and lead to sojourn duration in the regimes which are not geometric. In order to generate the prevailing rotation, it is necessary to command the wind direction to turn in the right direction inside the regimes. In {u t , v t } models the rotations are reproduced by the autoregressive A matrices, but they are specified more naturally in NHMS-EVM model by using a complex parametrization of the von Mises autoregressive models.

The proposed models allow to generate wind time series with features very close to the main features of the observed time series. The introduction of the latent state allows to simulate the different time scale which are present in the data, with the autoregressive part describing the short-term fluctuations whereas the weather type, which lasts typically a few days, describes longer-term fluctuations and is related to large scale circulation. Another layer could be added to simulate shorter time scales for very local features. The model could also be extended to a space-time model in several ways. For this, it will probably be easier to work with the {u t , v t } model based on Gaussian distributions which can naturally handle a space-time information. Then several strategies could be considered for the weather type process which could be local, with a different weather type at each site, or regional with a common weather type for the different locations. With the first strategy one has to deal with a space-time process of latent discrete variables and this is challenging from both a modeling and computational point of view. The second strategy is probably simpler to implement but requires some space-time homogeneity in the data. These and other related modeling issues are currently investigated.

Figure 1 :

 1 Figure 1: Histogram of {U t } (top left), rose plot of {Φ t } (bottom left) and histograms of {u t } and {v t } and joint distribution of {u t , v t } (right). The lines on the scatter plots are levels of a non-parametric kernel estimate of the bivariate density. Results for the months of January.

Figure 2 :

 2 Figure 2: Example of time series of wind direction (top plot) and wind speed (bottom plot). The colors indicate the most likely regimes for the fitted NHMS-EVM model with 4 regimes (top plot) and Gaussian homogeneous MS-AR model for wind speed with 3 regimes (bottom plot). The regimes have been ordered according to the time variability (the darker the more variability).

1 for x = x ′ and π (x) 1 >

 11 y) for Leb -a.e. y with γ 0 for x ∈ {1, ..., M } then there exists a permutation τ : {1, ..., M } → {1, ..., M } such that γ

  NHMS-EVM 2 7568 5882 5872 5882 5968 6075 M(M+1)-1+6M

Figure 3 :

 3 Figure 3: Rose plot of the wind direction in the three regimes identified on wind speed by the NHMS-AR (U,Φ) model (top) and the NHMS-AR (u,v) model (bottom) .

Figure 4 :

 4 Figure 4: Matrix of non-homogeneous transitions of NHMS-AR (u,v) model (plain line) and NHMS-AR (U,Φ) model (dashed line)

Figure 5 :

 5 Figure5: Rose plot of the marginal distribution of wind direction for the various models with M = 4 (resp. M = 3) regimes for Φ t (resp.{u t , v t }) and order s = 2. From left to right: HMS-AR (u,v) , HMS-AR (U,Φ) , NHMS-AR (u,v) , NHMS-AR (U,Φ)

Figure 6 :

 6 Figure 6:Joint distributions of {u t , v t } for observed wind (left) and wind simulated with NHMS-AR (u,v) model (middle) and NHMS-AR (U,Φ) model (right)

Figure 7 :

 7 Figure 7: Correlation function of {U t } and cross-correlation functions between the time series {cos(Y t )} and {sin(Y t )} for the various models. The full grey line corresponds to the sample functions and the dashed line to the fitted model with a 95% prediction intervals (dotted line). The distributions for the fitted model was obtained by simulation. Time on the x-axis is expressed in days. (a): HMS-AR (u,v) , (b): HMS-AR (U,Φ) , (c): NHMS-AR (u,v) , (d): NHMS-AR (U,Φ)

Figure 8 :

 8 Figure 8: Autocorrelation functions of the time series {cos(Y t )} and {sin(Y t )} for the various models. The full grey line corresponds to the sample functions and the dashed line to the fitted model with a 95% prediction intervals (dotted line). The distributions for the fitted model was obtained by simulation. Time on the x-axis is expressed in days. (a): HMS-AR (u,v) , (b): HMS-AR (U,Φ) , (c): NHMS-AR (u,v) , (d): NHMS-AR (U,Φ)

Figure 9 :

 9 Figure 9: Frequency of anticlockwise rotations between successive observations for the various models with M = 4 regimes and autoregressive model of order s = 2. The grey line corresponds to the value obtained on the data (45.4 % of anticlockwise rotations against 54.6 % of clockwise rotations). The boxplots show the distributions for the fitted models. They were obtained by simulation (results based on 4400 time series of length 124). (a): HMS-AR (u,v) , (b): HMS-AR (U,Φ) , (c): NHMS-AR (u,v) , (d): NHMS-AR (U,Φ)

  Let s, M ≥ 1 be some integers. The sequence (X t , Y t t-s+1 ) t∈Z follows a NHMS-AR model if it is a Markov chain with values in {1, ..., M } × E such that • the conditional distribution of X t given the values of {X t ′ } t ′ <t and {Y t ′ } t ′ <t only depends on X t-1 and Y t-1 and we denote p 1

t } (top left), rose plot of {Φ t } (bottom left) and histograms of {u t } and {v t } and joint distribution of {u t , v t } (right). The lines on the scatter plots are levels of a non-parametric kernel estimate of the bivariate density. Results for the months of January.

2 Models

2.1 Non-homogeneous Markov-switching autoregressive models

Let X t ∈ {1, ..., M } represents the latent weather type and Y t denotes the observed wind conditions at time t. Throughout the article {Y t } will represent successively the bivariate process of Cartesian coordinates of wind in Section 2.3, the wind direction in Section 2.4, and finally {Y t } stands for the wind speed in Section 2.5. Let us write E for the space in which Y t takes values (E will respectively refer to R 2 , T and R + in the following sections). It will be useful to introduce notation Y t+u t := (Y t , ..., Y t+u ), y t+u t := (y t , ..., y t+u ) (as well as X t+u t , x t+u t ) for t > 0 and u > 0.

Hypothesis 1.

  t ′ } t ′ <t and {X t ′ } t ′ ≤t only depends on X

t and Y t-1 , . . . , Y t-s and this conditional distribution has a probability density function (p.d.f.) p 2 y t |x t , y t-1 t-s .

Table 1 :

 1 BIC values for the various fitted wind direction models. The last column gives the number of parameters. The bold value corresponds to the selected model.

Table 2 :

 2 BIC values for the various bivariate models. The last column gives the number of parameters.

Table

		BL	AR	NA0-NAO+ Total
	1	0.1832 0.0793 0.0621 0.0548 0.3794
	2	0.0897 0.0919 0.0996 0.1728 0.4540
	3	0.0080 0.0166 0.0316 0.1103 0.1666
	Total	0.2810 0.1878 0.1933 0.3380	1

Acknowledgments: We thank Julien Cattiaux (CNRS) for providing data and very helpfull comments on this work.


Fabienne Comte

Johanna Kappus

Density deconvolution from repeated measurements without symmetry assumption on the errors

Keywords: Nonparametric estimation, Density deconvolution, Repeated measurements, Panel data AMS Subject Classification: 62G05 62G07 62G20

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Density deconvolution is one of the classical topics in nonparametric statistics and has been extensively studied during the past decades. The aim is to identify the density of some random variable X, which cannot be observed directly, but is contaminated by some additional additive error ε, independent of X. A large amount of literature is available on the case where the distribution of the errors is perfectly known. To mention only a few of the various publications on this subject, we refer to [START_REF] Carroll | Optimal rates of convergence for deconvolving a density[END_REF], [START_REF] Stefanski | Rates of convergence of some estimators in a class of deconvolution problems[END_REF], [START_REF] Stefanski | Deconvoluting kernel density estimators[END_REF], [START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF], [START_REF] Efromovich | Density estimation for the case of supersmooth measurement errors[END_REF], [START_REF] Pensky | Adaptive wavelet estimator for nonparametric density deconvolution[END_REF], [START_REF] Comte | Penalized contrast estimator for adaptive density deconvolution[END_REF]. However, perfect knowledge of the error distribution is hardly ever realistic in applications. For this reason, the interest in deconvolution problems with unknown error distribution has grown. [START_REF] Meister | On the effect of misspecifying the error density in a deconvolution problem[END_REF] has investigated deconvolution with misspecified error distribution. [START_REF] Diggle | A Fourier Approach to Nonparametric Deconvolution of a Density Estimate[END_REF] replace the unknown characteristic function of the errors by its empirical counterpart and then apply standard kernel deconvolution techniques. The effect of estimating the characteristic function of the errors is then systematically studied by [START_REF] Michael | On the effect of estimating the error density in nonparametric deconvolution[END_REF]. Let us also mention [START_REF] Johannes | Deconvolution with unknown error distribution[END_REF] for deconvolution problems with unknown errors. The last mentioned publications have been working under the standing assumption that an additional sample of the pure noise is available. This is realistic in some practical examples. For example, if the noise is due to some measurement error, it is possible to carry out additional measurements in absence of any signal. However, in many fields of applications it is not realistic to assume that an additional training set is available. It is clear that, to make the problem identifiable, some additional information on the noise is required. In the present work, we are interested in the case where information can be drawn from repeated measurements of X, perturbed by independent errors. This framework is known as model of repeated measurements or panel data model. The observations are of the type

Y j,k = X j + ε j,k ; j = 1, • • • , n; k = 1, • • • , N,
where all X j and ε j,k are independent. This problem is relatively well-studied under the assumption that the distribution of the errors is symmetric. We refer to [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF], [START_REF] Comte | Deconvolution Estimation of Onset of Pregnancy with Replicate Observations[END_REF] and [START_REF] Kappus | Adaptive density estimation in deconvolution problems with unknown error distribution[END_REF]. In the present paper, we consider deconvolution from repeated observations when the symmetry assumption on the errors is no longer satisfied. The estimation strategies which have been developed for the symmetric error case cannot be generalized to this framework and a completely different approach is in order. The same problem has been investigated in earlier publications by [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF] and by [START_REF] Neumann | Deconvolution from panel data with unknown error distribution[END_REF]. The paper by Li and Vuong has two major drawbacks. On one hand, the rates of convergence presented therein are extremely slow, in comparison to the rate results which are usually found in deconvolution problems. On the other hand, the mentioned authors impose extremely restrictive assumptions on the target density and on the distribution of the noise, which are only met in some exceptional cases. Neumann succeeds in overcoming this second drawback and constructing consistent estimators under most general assumptions. However, rate results are not given in that paper so the question whether the convergence rates found by Li and Vuong can be improved has so far remained unanswered. Moreover, the estimator proposed by Neumann is only implicitly given and nonconstructive, so it is difficult to investigate the practical performance. In the present work, we study a fully constructive estimator, which is based on a modification of the original procedure by Li and Vuong. It is interesting to note that we are able to improve substantially upon the rates of convergence found by Li and Voung and, at the same time, dispose of most of their restrictive assumptions. Surprisingly, it can also be shown that our estimator outperforms, in some cases, the estimators which have been studied for the structurally simpler case of repeated observations with symmetric errors. This paper is organized as follows: In Section 2, we introduce the statistical model and define estimators for the target density, as well as for the residuals. In Section 3, we provide upper risk bounds and derive rates of convergence. In Section 4, we present some data examples. All proofs are postponed to Section 5.

Statistical model an estimation procedure

Let ε 1 and ε 2 be independent copies of a random variable ε and let X be independent of ε 1 and ε 2 . By Y , we denote the random vector

Y = (Y 1 , Y 2 ) = (X + ε 1 , X + ε 2 ). We observe n independent copies Y j = (Y j,1 , Y j,2 ), j = 1, • • • , n of Y .
The following assumptions are imposed on X and ε:

(A1) X and ε have a square integrable Lebesgue densities f X and f ε .

(A2) The characteristic functions ϕ ε (•) = E[e i•ε ] and ϕ X (•) = E[e i•X ] vanish nowhere. (A3) E[ε] = 0.
Our objective is to estimate f X and f ε . This statistical framework allows a straightforward generalization to the case where more than two observations of the noisy random variable X are feasible. However, for sake of simplicity and clarity, we content ourselves with considering the two dimensional case.

In the sequel, we denote by ψ the characteristic function of the two dimensional random vector Y , u1Y1+u2Y2) ].

ψ(u 1 , u 2 ) = E[e i(
By independence of X, ε 1 an ε 2 , the following holds for ψ:

ψ(u 1 , u 2 ) = E[e i(u1+u2)X e iu1ε1 e iu2ε2 ] = ϕ X (u 1 + u 2 )ϕ ε (u 1 )ϕ ε (u 2 ). (2.1)
From formula (2.1) one derives the following Lemma, which has been formulated and proved in [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF]. Lemma 2.1 is then the key to the construction of the estimator.

2.1 Lemma. Assume that E[|Y 1 |] < ∞ and E[ε] = 0.
Then ϕ X is determined by ψ via the following formula:

ϕ X (u) = exp u 0 ∂ ∂u1 ψ(0, u 2 ) ψ(0, u 2 ) du 2 .
Li and Vuong propose the following estimator of ϕ X : u1Yj,1+u2Yj,2) and u1Yj,1+u2Yj,2) denoting the empirical version of ψ and its first partial derivative. Given a kernel K and bandwidth h, the corresponding estimator of f X is

ϕ LV X (u) := exp u 0 ∂ ∂u1 ψ(0, u 2 ) ψ(0, u 2 ) du 2 , with ψ(u 1 , u 2 ) = 1 n n j=1 e i(
∂ ∂u 1 ψ(u 1 , u 2 ) = 1 n n j=1 iY j,1 e i(
f LV X h (x) = 1 2π e -iux ϕ LV X (u)F K h (u) du, with K h (•) = 1/h K(•/h) and with F K h (u) = e iux K h (x) dx denoting the Fourier transformation.
We propose a modified version of f LV X h . First of all, it is well known that small values of the denominator lead to unfavorable effects in the estimation procedure, so it is preferable to consider some regularized version of ψ. One possible approach is to replace ψ in the denominator by ψ + ρ with some Ridge-parameter to be appropriately chosen. See, for example, [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF]. However, following ideas in [START_REF] Michael | On the effect of estimating the error density in nonparametric deconvolution[END_REF], we prefer to define

ψ(0, u 2 ) = ψ(0, u 2 ) min{n 1/2 | ψ(0, u 2 )|, 1} .
and use 1/ ψ(0, u 2 ) as an estimator of 1/ψ(0, u 2 ). This leads to defining the following modified estimator of ϕ X :

ϕ mod X (u) = exp u 0 ∂ ∂u1 ψ(0, u 2 ) ψ(0, u 2 ) du 2 .
Next, we have to pay attention to the fact that, by definition, neither ϕ LV X nor ϕ mod X need to be characteristic functions and they may take values larger than one. Indeed, much of the complexity in the proofs presented in [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF] and many of the restrictive assumptions imposed therein are a consequence of the fact that there appears an unbounded exponential term in the definition of ϕ LV X which has to be controlled, leading to some Bernstein-type arguments and hence to the assumption that the supports are bounded. However, the quantity to be estimated is, in any case, a characteristic function, so the quality of the estimator can be improved by bounding the absolute value of ϕ mod X . These considerations lead to defining our final estimator of the characteristic function of X,

ϕ X (u) := ϕ mod X (u) max{1, | ϕ mod X (u)|} . (2.2)
Sometimes, one may not only be interested in the estimation of the target density itself, but also in the distribution of the residuals. The following holds true for the characteristic function of ε:

ϕ ε (u) = ψ(0, u) ϕ X (u) .
This quantity can hence be recovered, using a plug-in estimator. We set

ϕ X (u) = ϕ X (u) min{n 1/2 | ϕ X (u)|, 1}
and then

ϕ ε (u) := ψ(0, u) ϕ X (u) . (2.3)
Given a kernel K and bandwidth h > 0, the kernel estimators of f X and f ε corresponding to formula (2.2) and (2.3) are

f X h (x) = 1 2π e -iux ϕ X (u)F K h (u) du and f ε h (x) = 1 2π e -iux ϕ ε (u)F K h (u) du.
3 Risk bounds and rates of convergence

Non-asymptotic risk bounds

We start by analyzing the performance of f X . It is important to stress that we can dispose of most of the assumptions which have been imposed in earlier publications on the subject. Indeed the conditions on X and ε which are imposed in [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF], namely boundedness of the support of f X and f ε and nowhere vanishing characteristic functions, are violated for any distribution which is commonly studied in probability theory. In [START_REF] Bonhomme | Generalized nonparametric deconvolution with an application to earning dynamics[END_REF] an estimator is constructed under weaker assumptions on the distributions. But still, it is required in that paper that X have moments of all orders, which is certainly quite restrictive. Moreover, the rates which are found by those authors turn out to be even slower than the rate results presented in [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF]. It is interesting to note that we can substantially improve upon these results, even though our assumptions are much weaker.

In [START_REF] Neumann | Deconvolution from panel data with unknown error distribution[END_REF], an implicit estimator of f X is proposed. The strength of this approach lies in the fact that it is fully general. However, the price one has to pay is the lack of constructivity. The estimator is found as the solution to an abstract minimization problem, so the practical computation is not clear. Moreover, consistency of the estimator is shown, but rate results are not given, so nothing can be said about the quality of the procedure. Finally, [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF] and [START_REF] Comte | Deconvolution Estimation of Onset of Pregnancy with Replicate Observations[END_REF] have studied estimators in a repeated measurement model, but it is assumed in both papers that the distribution of the noise is symmetric. It is the main concern of the present publication to be able to dispose of the symmetry.

In the sequel, we impose the following mild regularity assumption on the characteristic function of X:

(A4) For some positive constant C X , ∀u, v ∈ R + : (v ≤ u) ⇒ (|ϕ X (u)| ≤ C X |ϕ X (v)|).
The following bound can be given on the mean integrated squared error:

3.1 Theorem. Let K be supported on [-1, 1]. Assume that (A1)-( A4) are satisfied and that for some positive integer p,

E[|Y 1 | 2p ] < ∞. Assume, moreover, that ϕ ′′ X ϕ ε is integrable. Then for some positive constant C depending only on p, E f X -f X h 2 L 2 ≤ 2 f X -K h * f X 2 L 2 + CC X G(X, ε, 1, 1/h) n 1/h -1/h |u| 0 1 |ϕ ε (z)| 2 dz du + CG(X, ε, p, 1/h) 1/h -1/h 1 n |u| 0 1 |ϕ Y (z)| 2 dz p du, with ϕ Y (z) = E[e izY2 ] = ψ(0, z) = ϕ X (z)ϕ ε (z)
and

G(X, ε, p, u) :=( ϕ ′′ X ϕ ε L 1 + E[ε 2 ] ϕ X ϕ ε L 1 + ϕ ′ X ϕ ε 2 L 2 ) p + |u| 0 | ∂ ∂u 1 log ψ(0, x)| 2 dx p + E[|Y 1 | 2p ]1 {p≥2} u p n p-1 + E 1 2 [|Y 1 | 2p ].
In analogy with (A4), we impose the following assumption on the characteristic function of the errors:

(A5) For some positive constant C ε , ∀u, v ∈ R + : (v ≤ u) ⇒ (|ϕ ε (u)| ≤ C ε |ϕ ε (v)|).
The following bound can then be given on the mean integrated squared error of f ε :

3.2 Theorem. Assume that K is supported on [-1, 1] and the assumptions (A1)-(A5) are met. Assume, moreover, that for some positive integer q ≥ 2, E[|Y 1 | 4q ] is finite. Then for some positive constant C depending only on q,

E f ε -f ε h 2 L 2 ≤ f ε -K h * f ε 2 L 2 + CC ε G(X, ε, 1, 1/h) n 1/h -1/h |u| 0 1 |ϕ X (z)| 2 dz du + G(X, ε, q, 1/h) n q-1 1/h -1/h 1 |ϕ X (u)| 2 |u| 0 1 |ϕ X (z)| 2 dz |u| 0 1 |ϕ Y (z)| 2 dz q-1 du + G(X, ε, 2, 1/h) 1/2 n 2 1/h -1/h 1 |ϕ X (u)| 2 |u| 0 1 |ϕ Y (z)| 2 dz du + G(X, ε, 2q, 1/h) 1/2 n q 1/h -1/h 1 |ϕ X (u)| 4 |u| 0 1 |ϕ Y (z)| 2 dz q du + 1 n 2 1/h -1/h 1 |ϕ X (u)| 4 du .
Discussion It is easily seen that the assumptions (A4) and (A5) are not very restrictive. They are met, for example, for normal or mixed normal distributions, Gamma distributions, bilateral Gamma distributions and many others. By a location shift, one can always ensure that E[ε] = 0. The integrability condition on ϕ ′′ X ϕ ε is also very mild. Under the above assumptions, it is automatically met if ϕ ε is integrable but can also be checked in most other cases.

The upper bound in Theorem 3.1 differs from the bounds which are commonly found in deconvolution problems in two ways: For one thing, there appears an additional inner integral in the variance term. This could be a consequence of the two-dimensional nature of the underlying problem. On the other hand, it is completely unexpected to find, in the second variance term, the characteristic function of the target density appearing in the denominator. On an intuitive level, this phenomenon could be understood as follows: To draw inference on X some information on the noise is required. However, in comparison to standard deconvolution problems, ε is itself an unobservable quantity and is contaminated by X. Consequently, X does not only play the role of a random variable of interest but also, with respect to the error term, the role of a contamination. This might explain the occurrence of ϕ X in the denominator.

Rates of convergence

In what follows, we derive rates of convergence under regularity assumptions on the target density f X and on the density f ε of the noise. For sake of simplicity, we assume in this section that K is the sinc-kernel,

F K = 1 [-1,1] . Let us introduce some notation: For ρ, C 1 > 0, β, c ≥ 0, C 2 ≥ 1, we denote by F u (C 1 , C 2 , c, β, ρ) the class of square integrable densities f such that the characteristic function ϕ(•) = e i•x f (x) dx satisfies ∀u, v ∈ R + : (u ≥ v) ⇒ (|ϕ(u)| ≤ C 2 |ϕ(v)|) (3.1) and ∀u ∈ R : |ϕ(u)| ≤ (1 + C 1 |u| 2 ) -β 2 e -c|u| ρ . If c = 0, the functions collected in F u (C 1 , C 2 , c, β, ρ) are called ordinary smooth. For c > 0, they are called supersmooth. By F ℓ (C 1 , C 2 , c, β, ρ),
we denote the class of square integrable densities for which (3.1) holds and, in addition,

∀u ∈ R : |ϕ(u)| ≥ (1 + C 1 |u| 2 ) -β 2 e -c|u| ρ .
For C 3 > 0, we denote by G(C 3 , p) the class of pairs (f X , f ε ) of square integrable densities for which the following conditions are met: For the characteristic function ϕ X of f X and ϕ ε of f ε ,

( ϕ ′′ X ϕ ε + E[ε 2 ]ϕ X ϕ ε L 1 + ϕ ′ X ϕ ε 2 L 2 ) p + E[|X + ε| 2p ] ≤ C 3 holds and moreover, (log ϕ X+ε ) ′ is square integrable, with (log ϕ X+ε ) ′ 2p L 2 ≤ C 3 . Finally, we use the short notation F u,ℓ (X, ε, p) = F u (C 1,X , C 2,X , c X , β X , ρ X ) × F ℓ (C 1,ε , C 2,ε , c ε , β ε , ρ ε ) ∩ G(C 3 , p) and F ℓ,u (X, ε, p) = F ℓ (C 1,X , C 2,X , c X , β X , ρ X ) × F u (C 1,ε , C 2,ε , c ε , β ε , ρ ε ) ∩ G(C 3 , p).

Estimation of the target density

We start by providing rates of convergence for the estimation of f X . Let p ≥ 2. We may limit the considerations to bandwidths h ≥ n -1/2 , so the term u p n p-1 appearing in the definition of G(X, ε, p, u) is readily negligible. We consider three different cases:

Case I: Ordinary smooth density with ordinary smooth errors , c

X = c ε = 0, β X > 1/2, β ε > 1/2.
Then the choice of the kernel, Theorem 3.1 and the definition of F u,ℓ (X, ε, p) give sup

(f X ,fε)∈F u,ℓ (X,ε,p) E f X -f X,h 2 L 2 = O(r n,h ) := O (1/h) γ1 + 1 n (1/h) γ2 + 1 n p (1/h) γ3 with γ 1 = -2β X + 1, γ 2 = 2β ε + 2, γ 3 = p(2β X + 2βε + 1) + 1.
Minimizing r n,h with respect to h yields for the optimal bandwidth h * ,

1/h * ≍ n 1 2βε +2(1+1/p)β X +1 . Plugging h * in gives sup (f X ,fε)∈F u,ℓ (X,ε,p) E f X -f X,h * 2 L 2 = O n - (2β X -1) 2βε +2(1+1/p)β X +1
.

Case II Ordinary smooth density with supersmooth errors,

β X > 1/2, c X = 0, c ε > 0. Then sup (f X ,fε)∈F u,ℓ (X,ε,p) E f X -f X,h 2 L 2 = O(r n,h ) := O (1/h) γ1 + 1 n (1/h) γ2 exp(2c ε (1/h) ρε ) + 1 n p (1/h) γ3 exp(2pc ε (1/h) ρε ) , with γ 1 = -2β X + 1, γ 2 = [(2β ε + 1 -ρ ε ) + + 1 -ρ ε ] + , γ 3 = [p(2β ε + 2β X + 1 -ρ ε ) + + 1 -ρ ε ] + .
Selecting h * as the minimizer of r n,h gives

1/h * = 1 2c ε (log n) - 1 2c ε log(log n) γ + O(1) 1/ρε . with γ = max{ γ 2 ρ ε + 2β X -1, 1/p( γ 3 ρ ε + 2β X -1)}.
From this we derive that sup

(f X ,fε)∈F u,ℓ (X,ε,p) E f X -f X,h * 2 L 2 = O (log n) -2β X -1 ρε .
Case III Supersmooth density with ordinary smooth errors, c X > 0, c ε = 0, β ε > 1/2. In this case, sup

(f X ,fε)∈F u,ℓ (X,ε,p) E f X -f X,h * 2 L 2 = O(r n,h ) := O (1/h) γ1 exp(-2c X (1/h) ρ X ) + 1 n (1/h) γ2 + 1 n p (1/h) γ3 exp(2pc X (1/h) ρ X , with γ 1 = -2β X + 1 -ρ X , γ 2 = 2β ε + 2, γ 3 = [p(2β X + 2β ε + 1 -ρ X ) + + 1 -ρ X ] + . Minimizing r n,h yields 1/h * = p 2c X (p + 1) (log n) - 1 2c X log(log n) γ + O(1) 1/ρ X with γ = γ 3 -γ 1 (p + 1)ρ X Then it holds that sup (f X ,fε)∈F u,ℓ (X,ε,p) E f X -f X,h * 2 L 2 = O n -p p+1 (log n) p/(p+1)γ 1 +1/(p+1)γ 3 ρ X
.

Estimation of the residuals

In analogy with the rates for the estimation of f X , we consider the following different cases:

Case I Both, f ε and f X are ordinary smooth, c X = c ε = 0, β X > 1/2, β ε > 1/2.
Then by Theorem 3.2 and by the definition of F ℓ,u (X, ε, p),

sup (f X ,fε)∈F ℓ,u (X,ε,p) E f ε -f ε,h 2 L 2 = O(r n,h ) := O (1/h) γ1 + 1 n (1/h) γ2 + 1 n 2 (1/h) γ3 + 1 n p-1 (1/h) γ4 + 1 n p (1/h) γ5 , with γ 1 = -2β ε + 1, γ 2 = 2β X + 2, γ 3 = 4β X + 2β ε + 2, γ 4 = 2(p + 1)β X + 2(p -1)β ε + p + 1, γ 5 = 2(p + 2)β X + 2pβ ε + p + 1. Minimizing r n,h gives 1/h * ≍ n 1 2(1+2/(p-1))β X +2(1+1/(p-1))βε +1+1/(p-1) . Consequently, sup (f X ,fε)∈F ℓ,u (X,ε,p) E f ε -f ε,h * 2 L 2 = O n - 2βε -1 2(1+2/(p-1))β X +2(1+1/(p-1))βε +1+1/(p-1) . Case II Ordinary smooth f ε and supersmooth f X , c X > 0, c ε = 0, β ε > 1/2. In this case, sup (f X ,fε)∈F ℓ,u (X,ε,p) E f ε -f ε,h 2 L 2 = O(r n,h ) :=O (1/h) γ1 + 1 n (1/h) γ2 exp(2c X (1/h) ρ X ) + 1 n 2 (1/h) γ3 exp(4c X (1/h) ρ X ) + 1 n p-1 (1/h) γ4 exp(2c X (p + 1)(1/h) ρ X ) + 1 n p (1/h) γ5 exp((2c X (p + 2)(1/h) ρ X ) , with γ 1 = -2β ε + 1, γ 2 = [(2β X + 1 -ρ X ) + + 1 -ρ X ] + , γ 3 = [(2β X + 2β ε + 1 -ρ X ) + + 2β X + 1 -ρ X ] + , γ 4 = [(p -1)(2β X + 2β ε + 1 -ρ X ) + + (2β X + 1 -ρ X ) + + 2β X + 1 -ρ X ] + , γ 5 = [p(2β X + 2β ε + 1 -ρ X ) + 4β X + 1 -ρ X ] + . Then 1/h * = (p -1) 2c X (p + 1) log n - 1 2c X log (log n) γ + O(1) 1/ρ X with γ = 1/(p + 1)( γ 4 ρ X + 2β ε -1). This implies sup (f X ,fε)∈F ℓ,u (X,ε,p) E f ε -f ε,h * 2 L 2 = O (log n) -2βε -1 ρ X .
Case III f ε is supersmooth and f X is ordinary smooth. Then sup

(f X ,fε)∈F ℓ,u (X,ε,p) E f ε -f ε,h 2 L 2 = O(r n,h ) := O (1/h) γ1 exp(-2c ε (1/h) ρε ) + 1 n (1/h) γ2 + 1 n 2 (1/h) γ3 exp(2c ε (1/h) ρε ) + 1 n p-1 (1/h) γ4 exp(2c ε (p -1)(1/h) ρε ) + (1/h) γ5 n p exp(2c ε q(1/h) ρε ) ,
with

γ 1 = -2β ε + 1 -ρ ε , γ 2 = 2β X + 2, γ 3 = [2β X + 2β ε + 1 -ρ ε ) + + 2β X + 1 -ρ ε ] + , γ 4 = [(p -1)(2β X + 2β ε + 1 -ρ ε ) + + 4β X + 2 -ρ ε ] + , γ 5 = [p(2β X + 2β ε + 1 -ρ ε ) + + 4β X + 1 -ρ ε ] + .
We arrive at

1/h * = (p -1) p2c ε log n - 1 2c ε log(log n) γ 1/ρε , with γ = 1/p(γ 3 -γ 1 ) ρ ε which, in turn, implies sup (f X ,fε)∈F ℓ,u (X,ε,p) E f ε -f ε,h * 2 L 2 = O n -p-1 p (log n) p-1 p γ 1 + 1 p γ 3 ρε .
Discussion: We have not considered the case where both, the target density and the error density, are supersmooth. Deriving rates of convergence in this framework requires the consideration of various different subcases, leading to rather tedious and cumbersome calculations. We omit the details and refer to [START_REF] Lacour | Rates of convergence for nonparametric deconvolution[END_REF] for a detailed discussion on the subject.

Comparison to earlier results

We have mentioned that the rates of convergence derived above differ substantially from the rate results given in [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF]. To illustrate this point, the rates are listed in the table below.

c X = 0, c ε = 0 c X = 0, c ε > 0 c X > 0, c ε = 0 f X n - 2β X -1 2(1+1/p)β X +2βε +1 (log n) -2β X -1 ρε (log n) γ n -p p+1 f LV X n - 2β X -1 4β X +6βε +4 (log n) -2β X -1 ρε (log n) γ n -1 3

Table 1: Rates of convergence for estimating the target density

We need to be careful about the fact that the rates of convergence given in [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF] are derived under the assumption the moments of all orders and even all exponential moments are finite, which compares to p = ∞. There is no difference in the rate when an ordinary smooth target with supersmooth noise is being considered. In any other case, the gap in the rate is striking. It is interesting to note that the rates of convergence found in the present publication do also differ from the rates which have been found for estimators in the structurally simpler case of panel data

c X = 0, c ε = 0 c X = 0, c ε > 0 c X > 0, c ε = 0 f ε n - 2βε -1 2(p+1) p-1 β X + 2(p+1) p-1 βε + p+1 p-1 (log n) γ n -p-1 p (log n) -2βε -1 ρ X f LV ε n - 2β X -1 6β X +6βε +4 (log n) γ n -1 3 (log n) -2βε -1 ρ X
Table 2: Rates of convergence for estimating the noise density

c X = 0, c ε = 0 c X = 0, c ε > 0 c X > 0, c ε = 0 f X n - 2β X -1 2(1+1/p)β X +2βε +1 (log n) -2β X -1 ρε (log n) γ n -p p+1 f sym X n - 2β X -1 2(β X ∨βε )+2βε (log n) -2β X -1 ρε (log n) γ n -1+ε
Table 3: Rates of convergence for estimating f X , symmetric vs. non-symmetric error case.

with symmetric errors, see [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF] and [START_REF] Comte | Deconvolution Estimation of Onset of Pregnancy with Replicate Observations[END_REF]. In the table below, f sym X is understood to be the estimator for the symmetric error case, defined according to [START_REF] Comte | Deconvolution Estimation of Onset of Pregnancy with Replicate Observations[END_REF] and ε ∈ (0, 1/2) is arbitrary. The convergence rates coincide if an ordinary smooth target density with supersmooth errors is being considered. When both, f X and f ε are ordinary smooth and β X ≥ β ε holds, f sym attains the rate n -2β X -1 2β X +2βε , which is known to be optimal in deconvolution problems. In this situation, f X shows a slightly worse performance than f sym X , which is not surprising in light of the fact that the model with non-symmetric errors has a more complicated structure. However, it is certainly surprising to notice that for β ε >> (1 + 1/p)β X + 1/2, the rates for f X are substantially better than the rates for f sym X .

4 Simulation studies

Some data examples

For the practical choice of the smoothing parameter, we use a leave-p-out cross validation strategy. We consider the parameter set M = {1, • • • , √ n} with each parameter m corresponding to the bandwidth 1/m. Given any subset

N := {n 1 , • • • , n p } ⊆ {1, • • • , n} of size p, we build an estimator ϕ N X of ϕ X based on the subsample (Y k ) k∈N , as well as an estimator ϕ -N X based on (Y k ) k ∈N . For m ∈ M , we may use ℓ(ϕ X , ϕ X,1/m ) := 1 n p N ={n1,••• ,np} ϕ N X F K 1/ √ n -ϕ -N X F K 1/m 2 L 2
as an empirical approximation to the loss function

ℓ(ϕ X , ϕ X,1/m ) = ϕ X -ϕ X,1/m 2 L 2 .
Minimizing the empirical loss leads to selecting

m = argmin{m ∈ M : ℓ(ϕ X , ϕ X,1/m )},
and working with the bandwidth h = 1/ m, thus defining the adaptive estimator ϕ ad X = ϕ X, h . Simulation experiments indicate that the procedure works reasonably well with p = 10. However, it is evident that even for small sample sizes, the complexity of the algorithm explodes and the procedure is numerically intractable. To deal with this problem, we use a modified algorithm. We subdivide {1, • • • , n} into n/5 disjoint blocks B 1 , • • • , B n/5 of size 5 and build our leave-10-out estimators, based on the subsets

N k = B k ∪ B k+1 , k = 1, • • • , n/5 -1.
We work with a Gaussian kernel and try the procedure for the following target densities and errors:

(i) X has a Γ(4, 2) distribution and ε has a bilateral Gamma distribution with parameters 2, 2, 3, 3, that is, the corresponding density is the convolution of a Γ(2, 2)-density, supported on R + and a Γ(3, 3)-density, supported on R -. In the sequel, we abbreviate this type of distributions by bΓ(2, 2, 3, 3).

(ii) X has a bΓ(1, 1, 2, 2)-distribution and the errors are, up to a location shift, Γ(4, 2)-distributed, that is, ε + 2 ∼ Γ(4, 2). In the sequel, we write ε ∼ Γ(4, 2) -2. (The location shift is necessary to ensure that E[ε] = 0 holds true.)

(iii) X has a standard normal distribution, X ∼ N(0, 1) and ε ∼ bΓ(2, 2, 3, 3).

(iv) X has again a standard normal distribution. ε is a mixture of two normal distributions with parameters -2, 1 and 2, 2. We use the notation ε ∼ mN(-2, 1, 2, 2).

We use a Gaussian kernel and run the procedure for n = 100, 1000, 10000 observations. Based on 500 repetitions of the adaptive procedure, we calculate the empirical risk r ad and compare this quantity to the empirical risk r or of the "estimator" with oracle choice of the bandwidth. The values are summarized in the table below.

X ∼ Γ(4, 2), ε ∼ bΓ(2, 2, 3, 3) X ∼ N(0, 1), ε ∼ bΓ(2,
X ∼ N(0, 1), ε ∼ mN(-2, 1, 2, 2) X ∼ bΓ(1, 1, 2, 2), ε ∼ Γ(4,

Comparison to the symmetric error case

We have mentioned that so far, the model of repeated observations has mainly been studied under the additional assumption that the error terms are symmetric. In Section 3, it turned out that our rates of convergence are, in some cases, better than the rate results presented in [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF] or [START_REF] Comte | Deconvolution Estimation of Onset of Pregnancy with Replicate Observations[END_REF]. So far, it is not clear if this gap in the rate is due to a sub-optimal upper bound in the mentioned papers or to a different performance of the estimators themselves.

Simulation studies indicate that the estimator which has been designed to handle the case of skew errors does indeed outperform, in some cases, the standard estimator for the symmetric error case.

Before having a look at some data examples, let us give a brief outline on the estimation strategy for the symmetric error case: In the panel data model, suppose that ε has a symmetric distribution. In this case,

Yj,1 -Yj,2 = εj,1 -εj,2 d = εj,1 + εj,2, j = 1, • • • , n.
Consequently, ϕY 1 -Y 2 = ϕ 2 ε . An unbiased estimator of ϕ 2 ε can then be built from the data set (Yj,1 -Yj,2)j=1,••• ,n. Taking square roots gives an estimator ϕε of ϕε and a regularized version of this estimator is plugged in the denominator. Again ϕY can be estimated directly from the data. For the details, we refer to [START_REF] Comte | Deconvolution Estimation of Onset of Pregnancy with Replicate Observations[END_REF]. In the sequel, we denote by ϕ sym X the estimator for the symmetric error case. As indicated by the theory, it turns out that ϕX performs substantially better than ϕ sym X if the error density is very smooth, in comparison to the target density. To illustrate this phenomenon, we have a look at the following examples:

(i) X has a Γ(2, 4)-distribution and ε ∼ bΓ(3, 2, 3, 2). (ii) X ∼ bΓ(1, 2, 1, 2) an ε ∼ bΓ(4, 3, 4, 3).
When we consider target densities which are very smooth, in comparison to the error density, the estimator discussed in the present paper does, on small or medium sample sizes, still perform slightly better than the estimator for the symmetric error case. However, this difference in the performance is small and vanishes completely as the sample size increases. For illustration, we consider the following examples:

(iii) X ∼ bΓ(4, 3, 4, 3) and ε ∼ bΓ(1, 2, 1, 2).

(iv) X ∼ N(0, 1) and ε ∼ bΓ(3, 5, 3, 5).

In the table below, based on 500 repetitions of the estimation procedure (with oracle choice of the bandwidth), we compare the empirical risk r or of ϕX to the empirical risk r sym,or of ϕ or X .

X ∼ Γ(2, 4), ε ∼ bΓ(3, 5, 3, 5) X ∼ bΓ(1, 2, 1, 2) ε ∼ bΓ(4,
X ∼ bΓ(4, 3, 4, 3), ε ∼ bΓ(1, 2, 1, 2) X ∼ N(0, 1), ε ∼ bΓ(3,

Conclusion:

Our simulation studies indicate that our estimator is, in some cases, preferable to the estimation procedures designed for the symmetric error case. In other cases, the performance of both procedures is practically identical. However, if the errors are unknown it is clear that in practical applications, one cannot be sure if the symmetry assumption on the errors is satisfied, so we conclude that it is preferable, in either case, to work with the procedure which is designed for the non-symmetric case.

Proofs

Proof of Theorem 3.1

We start by providing some auxiliary results to prepare the proof of Theorem 3.1. In the sequel, we use the following short notation.

R(u) := 1 ψ(u) - 1 ψ(u) ; c(u) := ∂ ∂u1 ψ(0, u) - ∂ ∂u1 ψ(0, u); cj(u) := iYj,1e iuY j,2 - ∂ ∂u1 ψ(0, u); b(u) := ψ(0, u) -ψ(0, u) and Ψ ′ (u2) := ∂ ∂u1 log ψ(0, u2).
Moreover,

∆(u) := u 0 ∂ ∂u 1 ψ(0, u2) ψ(0, u2) - ∂ ∂u 1 ψ(0, u2) ψ(0, u2) du2.
First, we consider the deviation of 1/ ψ from its target:

5.1 Lemma. It holds that for some positive constant C depending on p,

E 1 ψ(0, u) - 1 ψ(0, u) 2p ≤ C min n -p |ψ(0, u)| 4p , 1 |ψ(0, u)| 2p .
Proof. Consider first the case where |ψ(0, u)| ≥ n -1/2 . We start by observing that

E | b(u)| 2p ≤ 4 p E |ψ(0, u) -ψ(0, u)| 2p + E | ψ(0, u) -ψ(0, u)| 2p .
By Rosenthal's inequality, for some constant C depending on p,

E |ψ(0, u) -ψ(0, u)| 2p ≤ C n p .
Moreover, by definition of ψ,

E | ψ(0, u) -ψ(0, u)| 2p ≤ E | ψ(0, u)| + n -1/2 2p 1 {| ψ(0,u)|≤n -1/2 } ≤ 4 p n -p . Consequently, E | b(u)| 2p ≤ Cn -p . Now, 1 ψ(0, u) - 1 ψ(0, u) 2p = b(u) ψ(0, u) ψ(0, u) 2p ≤ 4 p | b(u)| 2p |ψ(0, u)| 4p + | b(u)| 4p |ψ(0, u)| 4p | ψ(0, u)| 2p .
We have

E | b(u)| 2p |ψ(0, u)| 4p ≤ C n -p |ψ(0, u)| 4p and, since 1/|ψ(0, u)| ≤ √ n by definition, E | b(u)| 4p |ψ(0, u)| 4p | ψ(0, u)| 2p ≤ C n -p |ψ(0, u)| 4p .
On the other hand, for |ψ(0, u)| ≤ n -1/2 , we have the series of inequalities

E 1 ψ(0, u) - 1 ψ(0, u) 2p ≤ 4 p 1 |ψ(0, u)| 2p + E 1 | ψ(0, u)| 2p ≤ 4 p 1 |ψ(0, u)| 2p + n p ≤4 p 2 |ψ(0, u)| 2p .
This completes the proof.

The following result gives control on ∆:

5.2 Lemma. Assume that E[|Y1| 2p ] < ∞. Then for some positive constant C, E |∆(u)|1 {|∆(u)|>1} ≤ CG(X, ε, u, p) 1 n p |u| 0 1 |ψ(0, u2)| 2 du2 p with G(X, ε, u, p) =( ϕ ′′ X ϕε + E[ε 2 ]ϕX ϕε L 1 + ϕ ′ X ϕε 2 L 2 ) p + |u| 0 |Ψ ′ (x)| 2 dx p + u p 1 {p≥2} E[|Y1| 2p ] n p-1 + E 1/2 [|Y1| 2p ]. Moreover E |∆(u)| 2p 1 {|∆(u)|≤1} ≤ CG(X, ε, u, p) 1 n p |u| 0 1 |ψ(0, u2)| 2 du2 p . Proof. We can estimate ∆(u) = u 0 ∂ ∂u 1 ψ(0, u2) ψ(0, u2) - ∂ ∂u 1 ψ(0, u2) ψ(0, u2) du2 ≤ u 0 ∂ ∂u 1 ψ(0, u2) -∂ ∂u 1 ψ(0, u2) ψ(0, u2) du2 + u 0 ∂ ∂u1 ψ(0, u2) R(u2) du2 + u 0 ∂ ∂u1 ψ(0, u2) - ∂ ∂u1 ψ(0, u2) R(u2) du2 =:∆1(u) + ∆2(u) + ∆3(u).
Rosenthal's inequality (see, for example, [START_REF] Ibragimov | The exact constant in the Rosenthal inequality for random variables with mean zero[END_REF] ) gives for some constant C depending only on p:

E ∆1(u) 2p =E u 0 ∂ ∂u 1 ψ(0, u2) -∂ ∂u 1 ψ(0, u2) ψ(0, u2) du2 2p = E 1 n n j=1 u 0 cj(u2) ψ(0, u2) du2 2p ≤ C( 1 n p E u 0 cj(u2) ψ(0, u2) du2 2 p + 1 n 2p-1 E u 0 cj(u2) ψ(0, u2) du2 2p
).

Using Fubini's theorem, the Cauchy-Schwarz inequality and Lemma 6.1, we derive that

E u 0 cj(u2) ψ(0, u2) du2 2 = u 0 u 0 Cov( cj(x), cj(y)) ψ(0, x)ψ(0, -y) dx dy = u 0 u 0 E[(iY1) 2 e i(x-y)Y 2 ] ψ(0, x)ψ(0, -y) dx dy - u 0 u 0 E[iY1e ixY 2 ] E[iY1e -iyY 2 ] ψ(0, x)ψ(0, -y) dx dy ≤ u 0 u 0 | E[(iY1) 2 e i(x-y)Y 2 ]| |ψ(0, x)| 2 dx dy ≤ sup x∈[0,u] u 0 | E[(iY1) 2 e i(x-y)Y 2 | dy u 0 1 |ψ(0, y)| 2 dy ≤ ϕ ′′ X ϕε + E[ε 2 ]ϕX ϕε L 1 u 0 1 |ψ(0, y)| 2 dy.
For p ≥ 2, the Cauchy-Schwarz inequality gives

1 n 2p-1 E u 0 cj(u2) ψ(0, u2) du2 2p ≤ 1 n 2p-1 u 0 1 |ψ(0, u2)| 2 du2 p E   u 0 | cj| 2 du2   p ≤ 4 p E[|Yj| 2p ] n 2p-1 u p u 0 1 |ψ(0, u2)| 2 du2 p .
We have thus shown

E ∆1(u) 2p ≤ CG(X, ε, p, u) 1 n |u| 0 1 |ψ(0, u2)| 2 du2 p .
Next, thanks to Lemma 5.1 and the Hölder inequality,

E ∆2(u) 2p =E u 0 ∂ ∂u1 ψ(0, u2) R(u2) du2 2p = E u 0 Ψ ′ (u2)ψ(0, u2) R(u2) 2p ≤ |u| 0 |Ψ ′ (x)| 2 dx p E |u| 0 |ψ(0, u2)| 2 | R(u2)| 2 du2 p ≤ |u| 0 |Ψ ′ (x)| 2 dx p |u| 0 1 |ψ(0, u2)| 2 du2 p-1 |u| 0 1 |ψ(0, u2)| 2 |ψ(0, u2)| 4p E | R(u2)| 2p du2 ≤ C n p |u| 0 |Ψ ′ (x)| 2 dx p |u| 0 1 |ψ(0, u2)| 2 du2 p ≤ CG(X, ε, p, u) 1 n |u| 0 1 |ψ(0, u2)| 2 du2 p .
Finally, another application of Lemma 5.1 and the Hölder inequality gives

E ∆3(u) p =E u 0 ∂ ∂u1 ψ(0, u2) - ∂ ∂u1 ψ(0, u2) R(u2) du2 p = E u 0 c(u2) R(u2) du2 p ≤ |u| 0 1 |ψ(0, u2)| 2 du2 p-1 |u| 0 |ψ(0, u2)| 2p E | c(u2) R(u2)| p |ψ(0, u2)| 2 du2 ≤ |u| 0 1 |ψ(0, u2)| 2 du2 p-1 |u| 0 |ψ(0, u2)| 2p E 1 2 | c(u2)| 2p E 1 2 | R(u2)| 2p |ψ(0, u2)| 2 du2 ≤C E 1 2 [|Y1| 2p ] n p |u| 0 1 |ψ(0, u2)| 2 du2 p ≤ CG(X, ε, p, u) 1 n |u| 0 1 |ψ(0, u2)| 2 du2 p . We set Aj = Aj(u) := {|∆(u)| > 1} ∩ {argmax k=1,2,3 |∆ k (u)| = j}.
We may use the fact that on Aj, ∆(u) ≤ 3∆j(u) as well as ∆j(u) > 1/3, to conclude that

E |∆(u)|1 {|∆(u)|>1} ≤ 3 E |∆1(u)|1A 1 + E |∆2(u)|1A 2 + E |∆3(u)|1A 3 ≤ 3 2p E |∆1(u)| 2p + E |∆2(u)| 2p + E |∆3(u)| p .
Combining this inequality with the moment bounds on the ∆j(u), we have shown that for a constant C depending only on p,

E |∆(u)|1 {|∆(u)|>1} ≤ CG(X, ε, p, u) 1 n p 1 |ψ(0, u2)| 2 du2 p . Next, we define Bj := {|∆(u)| ≤ 1} ∩ {max{∆ k (u)|k = 1, 2, 3} = j}, j = 1, 2, 3. It holds that E |∆(u)| 2p 1 {|∆(u)|≤1} ≤ 9 p E[∆1(u) 2p 1B 1 ] + E[∆2(u) 2p 1B 2 ] + E |∆(u)| p 1B 3 ≤ 9 p E[∆1(u) 2p 1B 1 ] + E[∆2(u) 2p 1B 2 ] + E[∆3(u) p 1B 3 ] .
This implies, using again the moment bounds on the ∆j,

E |∆(u)| 2p 1 {|∆(u)|≤1} ≤ CG(X, ε, p, u) 1 n |u| 0 1 |ψ(0, u2)| 2 du2 p .
We can now prove the upper bounds for f X,h : Proof of Theorem 3.1. Parseval's identity gives

E f -f h 2 L 2 ≤ 2 f -K h * f 2 L 2 + 1 π |F K h (u)| 2 E |ϕX (u) -ϕX (u)| 2 du.
We use the trivial observation that |ϕX (u) -ϕX (u)| ≤ |ϕX (u) -ϕ mod X (u)|, as well as the fact that for z

∈ C with |z| ≤ 1, |1 -exp(z)| ≤ 2|z| holds, to derive that |ϕX (u) -ϕX (u)| 2 1 {|∆(u)|≤1} ≤ |ϕX (u) -ϕ mod X (u)| 2 1 {|∆(u)|≤1} =|ϕX (u)(1 -ϕ mod X (u)/ϕX (u))| 2 1 {|∆(u)|≤1} = |ϕX (u)(1 -exp(∆(u))| 2 1 {|∆(u)|≤1} ≤2|ϕX (u)∆(u)| 2 1 {|∆(u)|≤1} .
On the other hand, using the fact that |ϕX (u) -ϕX (u)| ≤ 2, as well as the Markov-inequality, we can estimate

|ϕX (u) -ϕX (u)| 2p 1 {|∆(u)|>1} ≤ 4 p |∆(u)|1 {|∆(u)|>1} .
Lemma 5.1, Lemma 5.2 and (A4) thus give

E |ϕX (u) -ϕX (u)| 2 ≤ 2|ϕX (u)| 2 E |∆(u)| 2 1 {|∆(u)|≤1} + 4E |∆(u)|1 {|∆(u)|>1} ≤ CG(X, ε, 1, u) n |ϕX (u)| 2 |u| 0 1 |ψ(0, u2)| 2 du2 + CG(X, ε, p, u) 1 n |u| 0 1 |ψ(0, u2)| du2 p ≤ CCX G(X, ε, 1, u) n |u| 0 1 |ϕε(u2)| 2 du2 + CG(X, ε, p, u) 1 n |u| 0 1 |ψ(0, u2)| du2 p .
Hence, by assumption on the support of K,

1/h -1/h | K h (u)| 2 E |ϕX (u) -ϕX (u)| 2 du ≤ CCX G(X, ε, 1, 1/h) n 1/h -1/h |u| 0 1 |ϕε(z)| 2 dz du + CG(X, ε, p, 1/h) n p 1/h -1/h |u| 0 1 |ψ(0, u2)| du2 p du.
This completes the proof.

Proof of Theorem

3.2 5.3 Lemma. Let q ≥ p. Assume that E[|Y1| 2q ] < ∞.
Then for some positive constant C depending only on p and q,

E 1 ϕX (u) - 1 ϕX (u) 2p ≤C G(X, ε, p, u) |ϕX (u)| 2p 1 n u 0 1 |ψ(0, u2)| 2 du2 p + G(X, ε, q, u) |ϕX (u)| 4p n q-p u 0 1 |ψ(0, u2)| 2 du2 q + 1 n p |ϕX (u)| 4p . Proof. We have E 1 ϕX (u) - 1 ϕX (u) 2p = E |ϕX (u) -ϕX (u)| 2p |ϕX (u) ϕX (u)| 2p .
Using the definition of ϕ mod X , as well as the fact that | exp(z)| ≥ 1/e holds for z ∈ C, |z| ≤ 1, we derive that

| ϕ mod X (u)|1 {|∆(u)|≤1} = |ϕX (u)|| exp(∆(u))|1 {|∆(u)|≤1} ≥ 1/e|ϕX (u)|1 {|∆(u)|≤1} .
Consequently, by definition of ϕX and ϕX ,

| ϕX (u)|1 {|∆(u)|≤1} ≥ | ϕX (u)|1 {|∆(u)|≤1} = | ϕ mod X (u)|1 {| ϕ mod X (u)|≤1} + 1 {| ϕ mod X (u)|≥1} 1 {|∆(u)|≤1} ≥ 1 e |ϕX (u)|1 {|∆(u)|≤1} . Next, |ϕX (u) -ϕX (u)| 2p ≤ 4 p |ϕX (u) -ϕX (u)| 2p + | ϕX (u) -ϕX (u)| 2p
and it holds that

E | ϕX (u) -ϕX (u)| 2p = E | ϕX (u) -ϕX (u)| 2p 1 {| ϕ X (u)|≤n -1/2 } ≤ E (| ϕX (u)| + n -1/2 ) 2p 1 {| ϕ X (u)|≤n -1/2 } ≤ 4 p n -p .
We use Lemma 5.2 to conclude that

E |ϕX (u) -ϕX (u)| 2p |ϕX (u) ϕX (u)| 2p 1 {|∆(u)|≤1} ≤ E |ϕX (u) -ϕX (u)| 2p 1 {|∆(u)|≤1}
1/e 2p |ϕX (u)| 4p This completes the proof of the lemma.

Proof of Theorem 3.2 . By Parseval's inequality and by assumption on the support of K,

E fε -f ε,h 2 L 2 ≤ 2 fε -K h * fε 2 L 2 + 1 π 1/h -1/h E |ϕε(u) -ϕε(u)| 2 du.
+ n -1 |ϕX (u)| 4 ≤C|ϕε(u)| 2 G(X, ε, 1, u) n u 0 1 |ψ(0, z)| 2 dz + G(X, ε, q, u) |ϕX (u)| 2 n q-1 u 0 1 |ψ(0, z)| 2 dz q + 1 n|ϕX (u)| 2 ≤CCε G(X, ε, 1, u) n u 0 1 |ϕX (z)| 2 dz + G(X, ε, 1, u) |ϕX (u)| 2 n q-1 u 0 1 |ϕX (z)| 2 dz u 0 1 |ψ(0, z)| 2 dz q-1 + 1 n|ϕX (u)| 2 .
Finally, using the Cauchy-Schwarz inequality and again Lemma 5.3, we derive that Putting the above together, we have shown that for some positive constant C,

E |ψ(0, u) -ψ(0, u)| 2 1 ϕX (u) - 1 ϕX (u) 2 ≤ E 1 2 |ψ(0, u) -ψ(0, u)| 4 E 1 2 1 ϕX (u) - 1 ϕX (u)
1/h -1/h E ϕε(u) -ϕε(u) 2 du ≤ CCε G(X, ε, 1, 1/h) n 1/h -1/h |u| 0 1 |ϕX (z)| 2 dz du + G(X, ε, q, 1/h) n q-1 1/h -1/h 1 |ϕX (u)| 2 |u| 0 1 |ϕX (z)| 2 dz |u| 0 1 |ψ(0, z)| 2 dz q-1 du + G(X, ε, 2, 1/h) 1/2 n 2 1/h -1/h 1 |ϕX (u)| 2 |u| 0 1 |ψ(0, z)| 2 dz du + G(X, ε, 2q, 1/h) 1/2 n q 1/h -1/h 1 |ϕX (u)| 4 |u| 0 1 |ψ(0, z)| 2 dz q + 1 n 2 1/h -1/h 1 |ϕX (u)| 4 du ,
which gives the statement of the theorem.

Appendix

6.1 Lemma. The following holds for the partial derivatives of ψ:

∂ k ∂u k 1 ψ(0, u2) = E (iY1) k e iu 2 Y 2 = k m=0 k m E[(iε) k-m ]ϕε(u2)ϕ (k) X (u2).
Proof. By definition of ψ and by independence of X, ε1 and ε2,

∂ k ∂u k 1 ψ(0, u2) = E ∂ k ∂u k 1 e iu 1 Y 1 +iu 2 Y 2 u 1 =0 = E (iY1) k e iu 2 Y 2 = k m=0 k m E (iX) m (iε1) k-m e iu 2 X e iu 2 ε 2 = k m=0 k m E[(iε) k-m ] E[(iX) m e iu 2 X ] E[e iu 2 ε ] = k m=0 k m E[(iε) k-m ]ϕε(u2)ϕ (m)
X (u2).

≤

  2|ϕX (u)| 2p E[|∆(u)| 2p 1 {|∆(u)|≤1} ] + n -p 1/e 2p |ϕX (u)| 4p (u)| 4p .Next, using the fact that by definition of ϕX , |1/ ϕX (u)| ≤ √ n holds, as well as the fact (u) -ϕX (u)| |ϕX (u) ϕX (u)| , we conclude that E |ϕX (u) -ϕX (u)| 2p |ϕX (u) ϕX (u)| 2p 1 {|∆(u)|>1} ≤4 p E |ϕX (u) -ϕX (u)| 2p 1 {|∆(u)|>1} |ϕX (u)| 4p + n p E |ϕX (u) -ϕX (u)| 4p 1 {|∆(u)|>1} |ϕX (u)| 4p ≤8 p E |∆(u)|1 {|∆(u)|>1} + n -p |ϕX (u)| 4p + n p E |∆(u)|1 {|∆(u)|>1} + n -2p

.

  Since ψ(0, u) is a characteristic function and ψ(0, u) its empirical counterpart,E[|ψ(0, u) -ψ(0, u)| 2 ]|ϕX (