Jean-Pierre Florens
email: jean-pierre.florens@tse-fr.eu

Léopold Simar
email: leopold.simar@uclouvain.be

Ingrid Van Keilegom
email: ingrid.vankeilegom@kuleuven.be

Estimation of the Boundary of a Variable observed with Symmetric Error

Keywords: Characteristic function, cumulant function, flexible parametric family, frontier estimation, Laguerre polynomials

Consider the model Y = X + ε with X = τ + Z, where τ is an unknown constant (the boundary of X), Z is a random variable defined on R + , ε is a symmetric error, and ε and Z are independent. Based on a iid sample of Y we aim at identifying and estimating the boundary τ when the law of ε is unknown (apart from symmetry) and in particular its variance is unknown. We propose an estimation procedure based on a minimal distance approach and by making use of Laguerre polynomials. Asymptotic results as well as finite sample simulations are shown. The paper also proposes an extension to stochastic frontier analysis, where the model is conditional to observed variables. The model becomes Y = τ (w 1 , w 2 ) + Z + ε, where Y is a cost, w 1 are the observed outputs and w 2 represents the observed values of other conditioning variables, so Z is the cost inefficiency. Some simulations illustrate again how the approach works in finite samples, and the proposed procedure is illustrated with data coming from post offices in France.

Introduction

The objective of this paper is to identify and estimate the lower endpoint of the support of a random variable that is bounded from below and that is observed with error when 'minimal' assumptions are imposed on the error term, in the sense that the error is only assumed to be symmetric around zero with no distributional assumptions. The problem of estimating the minimum of a random variable that is observed with noise, has received a lot of attention in the literature. Consider the relation Y = X +ε = τ +Z +ε, where Y is the observed variable, X is the variable of interest, τ is the lower endpoint of the support of X, Z ≥ 0, ε is the noise, and Z and ε are independent. When the distribution of ε is known, a large body of papers exists on the estimation of the support of the variable X, most of which are based on deconvolution techniques. See e.g. [START_REF] Goldenshluger | Estimating the endpoint of a distribution in presence of additive observation errors[END_REF], [START_REF] Delaigle | Data-driven boundary estimation in deconvolution problems[END_REF], Meister (2006a), and Aarts, Groeneboom and Jongbloed (2007), among others. The assumption that the distribution of ε is given, is a strong hypothesis (see [START_REF] Daouia | Robust frontier estimation from noisy data: a Tikhonov regularization approach[END_REF]). However, it has been shown that thanks to the positivity of Z the variance of ε is identified if ε is normally distributed (see [START_REF] Schwarz | Consistent density deconvolution under partially known error distribution[END_REF], and in that case the minimum τ may be estimated nonparametrically (see Hall andSimar, 2002, and[START_REF] Kneip | Frontier estimation in the presence of measurement error with unknown variance[END_REF]. Also note that a different but related literature is on the identification and estimation of the density of X under certain smoothness assumptions on the density of X and ε. This has been studied in [START_REF] Butucea | Minimax estimation of the noise level and of the deconvolution density in a semiparametric convolution model[END_REF], [START_REF] Meister | Density estimation with normal measurement error with unknown variance[END_REF][START_REF] Meister | Deconvolving compactly supported densities[END_REF], [START_REF] Butucea | Adaptivity in convolution models with partially known noise distribution[END_REF], [START_REF] Schwarz | Consistent density deconvolution under partially known error distribution[END_REF] and [START_REF] Delaigle | Methodology for non-parametric deconvolution when the error distribution is unknown[END_REF], among others. This paper goes one step further in the sense that we weaken the hypotheses on the stochastic error ε: we only assume that ε is symmetric around zero. Our approach is in the line of [START_REF] Delaigle | Methodology for non-parametric deconvolution when the error distribution is unknown[END_REF], who estimate the density of X using the empirical phase function. See also [START_REF] Butucea | Semiparametric mixtures of symmetric distributions[END_REF] and [START_REF] Butucea | Semiparametric topographical mixture models with symmetric errors[END_REF] for somewhat related papers, in which ε has a symmetric density and X is discrete. The price of this weak assumption on ε is that the distribution of X is no longer identified. As we will show below, only the odd cumulants of the distribution of X are identified, which is not sufficient to identify in a unique way the distribution. However, a precise parametric but flexible approximation of the density of X may be estimated from the odd cumulants. This model may be extended to conditional models. Our main motivation for such an extension is to analyze the so-called stochastic frontier models. Let us assume that X is a cost generated conditionally on variables W = (W 1 , W 2 ). The distribution of X is conditional to W 1 ≥ w 1 in R p 1 and to W 2 = w 2 in R p 2 . In econometric applications W 1 represents the outputs and W 2 describes the conditions of the production process. We are interested in the minimal cost given that W 1 ≥ w 1 and W 2 = w 2 , which we denote by τ (w 1 , w 2 ). If we define Z = X -τ (w 1 , w 2 ), the usual expression of the stochastic frontier model is

Y = τ (w 1 , w 2 ) + Z + ε,
where Z is positive, the error ε has a symmetric but unknown density, and ε and Z are independent conditionally on W 1 ≥ w 1 and W 2 = w 2 . This model has been studied extensively in the case where Z, ε and τ have a parametric form (for example, a Cobb-Douglas model for τ , an exponential distribution for Z, and a normal distribution for ε), see Aigner et al. (1977), Meeusen and van den Broek (1977). Some semiparametric approaches have also been proposed where for instance the frontier function is unspecified but the stochastic part is fully parametric (see e.g. [START_REF] Fan | Semiparametric estimation of stochastic production frontier models[END_REF] or [START_REF] Kuosmanen | Stochastic non-smooth envelopment of data: Semiparametric frontier estimation subject to shape constraints[END_REF]. We extend this analysis to the case where the density of ε is only assumed to be symmetric.

The paper is organized as follows. In Section 2 we discuss the identifiability of the marginal model. Estimation in the case without conditioning variables is treated in Section 3. A tractable flexible parametric family of distributions based on Laguerre polynomials is considered in Section 4, together with finite sample simulations. Section 5 presents the extension to frontier estimation, with some simulations and an illustration based on data coming from post offices in France. Conclusions and further research are given in Section 6.

The proofs of the main results are gathered in Appendix A, whereas Appendix B collects some properties of Laguerre polynomials. Finally, Appendix C gives some details about the practical implementation of our estimators.

Identifiability of the model

We consider the following model. We observe the random variable Y defined as

Y = X + ε = τ + Z + ε, (2.1)
where Z is an unobservable random variable with support starting at 0 and with density g, τ is an unknown constant corresponding to the lower endpoint of the support of the random variable X = τ + Z, Y has density f , the error ε has a symmetric density h (around zero), and ε and Z are independent. The model parameters are τ , g and h, but we are mostly interested in estimating τ . Note that we do not specify the distribution of ε, nor its variance, we only assume symmetry. We suppose throughout the paper that all variables are continuous and have an analytic characteristic function, and in particular that all moments exist (see Lukacs, 1960, Chapter 7). The characteristic function of Y (and similarly for all other variables) will sometimes be written as ψ Y and sometimes as ψ f , depending on the context.

Let θ 0 = (τ 0 , g 0 , h 0 ) be the true parameter values, and let Θ = R × G × H, where G is the class of densities living on [0, ∞), and H is the class of symmetric densities around zero.

We suppose that θ 0 ∈ Θ.

We know that the true density f of Y satisfies

f (y) = +∞ -∞ g 0 (y -τ 0 -e)h 0 (e) de (2.2)
for all y. The question we like to answer is how to identify τ, g and h from (2.2) given that g lives on [0, ∞) and that h is symmetric, or in other words under which conditions does the following hold true : if for some (τ, g, h) ∈ Θ we have that f (y) = g(y -τ -e)h(e) de for all y, then necessarily τ = τ 0 , g = g 0 and h = h 0 ?

It is clear that additional assumptions are needed to identify the model. Consider for example model (2.1) with τ = 2, Z ∼ Exp(1) and ε distributed according to the independent sum of a Un[-1, 1] variable and a standard normal variable. Then, exactly the same density f will be obtained if τ = 1, Z is distributed according to the independent sum of a Exp [START_REF] Aarts | Estimating the upper support point in deconvolution[END_REF] and a Un[0, 2] variable, and ε is standard normal. More generally, if in model (2.1) we take

τ = τ -c, Z = c + Z + U, ε = ε -U,
with c ≥ 0, U independent of ( Z, ε), Z and ε independent, the density of U symmetric around zero and the support of U equal to [-c, c], then we can equally well write Y as Y = τ + Z + ε.

Note that this shows that the constraints on Z to make the model identifiable are related to the concept of indecomposability of Z, which means that Z cannot be decomposed as the sum of two independent non-constant random variables (see [START_REF] Lukacs | Characteristic Functions[END_REF][START_REF] Lukacs | Developments in Characteristic Function Theory[END_REF] for more details about indecomposability of a distribution). Here it is necessary that such a decomposition is not possible if one of the two variables is required to be symmetric, which we call symmetric indecomposability. The following lemma gives sufficient conditions on a positive random variable Z to be symmetric indecomposable.

Lemma 2.1 Let ψ Z be the characteristic function of a positive random variable Z, and suppose that ψ Z has at most a finite number of zeros. Then, Z is symmetric indecomposable.

Proof. Suppose that Z would be symmetric decomposable and could be written as Z = Z + U with U symmetric around zero and Z independent of U . Since Z > 0, we necessarily have that |U | ≤ a for some finite a > 0 and Z ≥ a, because min Z = min Z + min U = 0.

Then, by Theorem 7.2.3 in [START_REF] Lukacs | Characteristic Functions[END_REF], the characteristic function ψ U of U has an infinite number of zeros. Since ψ Z = ψ Z ψ U , this contradicts the assumption that ψ Z has at most a finite number of zeros. 2

So, for example, the exponential, the half-normal and the truncated normal distribution are symmetric indecomposable, since ψ Z (t) has no zeros in these cases. The same is true for the gamma distribution or for finite linear combinations of these, since ψ Z (t) is in these cases a polynomial of finite order in t with at most a finite number of zeros. So the above lemma indicates that the symmetric indecomposability of a distribution defined on the positive real line is rather common. But still this is not a universal property as the simple example above reminds us.

In order to formulate conditions that guarantee the identifiability of our model, we follow a somewhat similar approach as Delaigle and Hall (2016), who studied the identifiability of the model Y = X + ε, where ε has a symmetric density (as in our case), but they do not assume that the support of X has a finite lower bound. Assumption (2.3) in their paper is a crucial assumption to identify the law of X, and states that the distribution of X uniquely has least variance among all distributions sharing the phase function ψ Y / ψ Y , where • is the modulus. Note that this implies that X and hence Z must be symmetric indecomposable. 1One way to identify our model would therefore be to impose their assumption (2.3) restricted to the class F of densities living on [τ, ∞) with τ unknown, since τ is identified as soon as the law of X is identified. However, we are in a more particular framework than Delaigle and Hall (2016) and we feel that it is more natural to disentangle Z and τ in our model.

Therefore, we will look for alternative assumptions to identify the density g of Z, since once g is identified we know τ via the relation

E(Y ) = τ + E(Z).
Our approach is based on cumulants, where for = 1, 2, . . ., the -th cumulant of X is defined as κ (f ) = K ( ) f (0), and the cumulant generating function K f (t) is given by

K f (t) = log E[exp(tY )].
Using standard properties of cumulants, we have that

κ 1 (f ) = E(Y ) = τ + E(Z) + E(ε) = τ + κ 1 (g).
On the other hand, for p ≥ 1,

κ 2p+1 (f ) = κ 2p+1 (δ τ ) + κ 2p+1 (g) + κ 2p+1 (h) = κ 2p+1 (g),
with δ τ the Dirac measure at τ , since κ 2p+1 (δ τ ) = 0 and since odd cumulants of symmetric densities are zero (because of formula (2.3) below, and the fact that the characteristic function of a symmetric function is real). So the odd cumulants of g of order 3 or higher are identifiable. We also know that again for p ≥ 1,

κ 2p (f ) = κ 2p (g) + κ 2p (h),
and since the even cumulants of h do not vanish in general, the even cumulants of g are not identifiable. We therefore introduce a class of densities of Z denoted by G, that satisfies the following hypothesis :

(A1) The density g 0 belongs to G, and if κ 2p+1 (g) = κ 2p+1 (g 0 ) for p ≥ 1 and for some g ∈ G, then g = g 0 .

Then as κ 2p+1 (f ) = κ 2p+1 (g 0 ), the model restricted to the class G becomes identified.

In particular we will consider cases where G is parametrized by a vector of parameters λ ∈ Λ of finite dimension :

G = {g(•|λ) : λ ∈ Λ ⊂ R m , Λ compact}.
In [START_REF] Delaigle | Methodology for non-parametric deconvolution when the error distribution is unknown[END_REF] identification is obtained by the selection of the distribution of smallest variance in the class of equivalent densities (i.e. defining the same f ).

In order to formulate assumption (A1) in a different way, we introduce now the function η f (and similarly for η g ) :

η f (t) = ∂ ∂t log ψ f (t) ,
where (a + ib) = b for any complex number a + ib, and where for any complex number z = ρ exp(iθ) with θ ∈ (-π, π], the logarithm of z is defined as log z = log ρ + iθ. Note that

log ψ f (t) = ∞ p=1 (it) p p! κ p (f ) = ∞ p=1 (-t 2 ) p (2p)! κ 2p (f ) + i ∞ p=1 (-1) p+1 t 2p-1 (2p -1)! κ 2p-1 (f ), (2.3)
and hence

log ψ f (t) = ∞ p=1 (-1) p+1 t 2p-1 (2p -1)! κ 2p-1 (f ), or η f (t) = κ 1 (f ) + ∞ p=2 (-1) p+1 t 2p-2 (2p -2)! κ 2p-1 (f ) = κ 1 (f ) + ∞ p=1 (-1) p t 2p (2p)! κ 2p+1 (f ).
This shows that we can write

η f (t) -η f (0) = ∞ p=1 (-1) p t 2p (2p)! κ 2p+1 (f ),
which depends only on κ 2p+1 (f ) for p ≥ 1. The above derivations now lead to the following equivalent formulation of assumption (A1) :

(A1') The density g 0 belongs to G, and if η g (t) -η g (0) = η g 0 (t) -η g 0 (0) for all t and for some g ∈ G, then g = g 0 .

Note that the mean (first cumulant) κ 1 (f ) is not identified, as it contains the unknown parameter τ . By considering the function η g (t) -η g (0) instead of the function η g (t) the mean cancels out, and we end up with an identified function. If we denote by T λ (t) the function η g (t) -η g (0) with g ≡ g(•|λ), the identification condition means that for all λ ∈ Λ,

T λ ≡ T λ 0 ⇒ λ = λ 0 .
In a non-linear context local identification is a useful concept and is defined as follows : there exists a neighborhood V of λ 0 such that

T λ ≡ T λ 0 for some λ ∈ V ⇒ λ = λ 0 .
In this paper we identify V without loss of generality with Λ. We denote now by κ (λ) the -th cumulant of g(•|λ). Then, the equality T λ ≡ T λ 0 is equivalent to the property that if κ 2p+1 (λ) = κ 2p+1 (λ 0 ) for all p ≥ 1 , then λ = λ 0 . The implicit function theorem implies that if κ 2p+1 (λ 0 ) is differentiable with respect to λ, then there exists a P ≥ m such that the matrix ∂ ∂λ j κ 2p+1 (λ) λ=0 p=1,...,P j=1,...,m has rank equal to m.

To conclude, if λ is identified, the distribution of Z is also identified, and so E(Z) and E(Z 2 ) are identified as well as a function of λ. Since τ 0 = E(Y ) -E(Z) and σ 2 0 = Var(ε) = Var(Y ) -Var(Z), the identification of τ 0 and σ 2 0 is guaranteed as well.

Estimation of the model

We start by explaining the methodology, followed by the asymptotic properties of the proposed estimators. Throughout this section the parametric family G can be any parametric class which is such that it contains the true density g(•|λ 0 ) and such that the cumulants of odd order starting from order 3 identify the densities in the class. In Section 4 we will work with one particular class of densities determined by Laguerre polynomials.

Methodology

Recall that the class G is a parametric class of densities g(•|λ) defined on R + and indexed by λ ∈ Λ ⊂ R m . We equip Λ with the canonical inner product λ, λ R m = m j=1 λ j λj (we will suppress the subscript R m in •, • R m whenever there is no ambiguity). Let η λ (t) = ∂ ∂t log ψ λ (t) , where ψ λ is the characteristic function corresponding to g(•|λ). We suppose that the function

T λ (•) = η λ (•) -η λ (0) belongs to the space E of functions from R to C such that v 1 , v 2 E = v 1 (t)v 2 (t)w(t) dt < ∞ for v 1 , v 2 ∈ E,
for a certain measure w, and where v is the complex conjugate of v. (Note that T λ (•) takes values in R, which we embed in C.)

We suppose that the Fréchet derivative T λ : Λ → E exists and that T λ and T λ are injective.

The linear operator T λ is bounded, and hence the dual operator T * λ , defined by

T λ (λ), v E = λ, T * λ (v) R m
exists, for any v ∈ E. Note that T * λ T λ : R m → R m is regular (since T λ is injective), but often ill-conditioned, and hence we will use penalization techniques to estimate λ 0 .

Define T (t) = η(t) -η(0), where η(t) = ∂ ∂t log ψY (t) and ψY (t) = n -1 n j=1 exp(itY j ), where Y 1 , . . . , Y n is an i.i.d. sample of Y . Then, the estimator of λ 0 is defined as

λα = argmin λ∈Λ +∞ -∞ T (t) -T λ (t) 2 w(t) dt + αP (λ) (3.1) = argmin λ∈Λ T -T λ 2 E + α λ 2 R m ,
where P (λ) = λ 2 and α > 0 is a regularization parameter. We could also work with other penalties like sparsity penalties (P (λ) = λ ) or entropy penalties (P (λ) = m j=1 λ j log |λ j |), but we do not consider them here for reasons of brevity. Note that instead of working with P (λ) = λ 2 , which favors models for which λ is close to zero, we could also work with P (λ) = λ -λ * 2 , if we have an a priori guess of the value of λ. Properties of λα are not affected by the choice of λ * . Remark 3.1 We assume that there exists a parametric family Λ of dimension m such that g = g(•|λ 0 ) for some λ 0 ∈ Λ. This assumption is not restrictive if m may be large (but finite). In that case, the model remains identified but severely ill-conditioned (i.e. the ratio of the largest and smallest eigenvalue of T * λ 0 T λ 0 is very high). The penalization solves the numerical problems caused by this ill-conditioning.

Note that for all λ ∈ Λ, λα satisfies

-T λα (λ), T -T λα + α λ, λα = -λ, T * λα ( T -T λα ) + α λ, λα = 0,
and hence

T * λα ( T -T λα ) -α λα = 0. (3.2)
For what follows, we need to define the regularized version λ α of λ 0 :

λ α = argmin λ∈Λ T λ 0 -T λ 2 E + α λ 2 R m ,
which satisfies (using a similar derivation as for λα )

T * λα (T λ 0 -T λα ) -αλ α = 0. (3.3)
Finally, note that since E(Z) = a(λ) and Var(Z) = b(λ) for some known functions a and b, we may estimate τ and σ 2 by

τα = Y -a( λα ) and σ2 α = σ2 Y -b( λα ), (3.4)
where Y and σ2 Y are the empirical average and variance of Y , respectively. Appendix C gives some details about the practical implementation of the whole method.

Asymptotic properties

The proofs of the results of this subsection are deferred to Appendix A, as well as the conditions under which these results are valid. We start with the weak convergence of the

process ( T (•) -T λ 0 (•)) as a process in E, equipped with •, • E . Theorem 3.1 The process n 1/2 ( T (•)-T λ 0 (•)) converges weakly in E to a mean-zero Gaussian process W (•) with covariance operator Γ given by (Γψ)(t) = γ(s, t)ψ(s)w(s)ds with γ(s, t) = σ(s, t) -σ(s, 0) -σ(0, t) + σ(0, 0), where σ(s, t) = 1 4 ∂ 2 ∂s∂t ψ Y (t -s) ψ Y (t)ψ Y (s) + ψ Y (t -s) ψ Y (t)ψ Y (s) - ψ Y (t + s) ψ Y (t)ψ Y (s) - ψ Y (t + s) ψ Y (t) ψ Y (s) .
Next, we consider the weak consistency and asymptotic normality of λα for fixed α > 0.

Note that the term ∇(T * λ S)| S=T λα -T λ 0 ,λ=λα in the formula of the asymptotic variance is nonstandard and is caused by the non-linearity of the operator T λ , where ∇ denotes the gradient vector with respect to λ. Theorem 3.2 Assume (C1)-(C4). Then, for fixed α > 0, λα → P λ α , and

n 1/2 ( λα -λ α ) d → N (0, Ω),
where

Ω = ∆ -1 T * λα ΓT λα ∆ -1 , ∆ = αI + T * λα T λα + ∇(T * λ S)| S=T λα -T λ 0 ,λ=λα .
The variance matrix Ω depends on T λ and T * λ , which can be calculated analytically, although their formula might be complicated. Let us consider T * λ . Define β jλ = T λ e j , where e j is the vector in R m with j-th component equal to one and all other components equal to zero. Then, for all λ ∈ Λ, T λ λ = T λ m j=1 λj e j = m j=1 λj β jλ , and therefore

λ, T * λ a R m = T λ λ, a E = m j=1 λj β jλ , a E = λ, β λ , a L 2 R m .
Hence, T * λ a = β λ , a . Note that the choice of α is an important but delicate issue. We refer to Theorem 10.7 in [START_REF] Engl | Regularization of Inverse Problems[END_REF] for a result based on balancing the variance and the squared bias of λα .

Next, we consider the asymptotic normality of τα and σ2 α .

Theorem 3.3 Assume (C1)-(C5). Then, for fixed α > 0,

n 1/2 (τ α -τ α ) d → N (0, v τ ) and n 1/2 (σ 2 α -σ 2 α ) d → N (0, v σ 2 ),
where

τ α = E(Y ) -a(λ α ), σ 2 α = Var(Y ) -b(λ α )
, and where v τ and v σ 2 are given in the proof in the Appendix.

We finish this section with the weak consistency of λα , τα and σ2 α when α tends to zero.

Theorem 3.4 Assume (C1)-(C4) and assume that α → 0 and (nα) -1 = O(1). Then, there exists a subsequence λα,ss of λα that converges in probability to λ 0 . Similarly, there exist subsequences τα,ss of τα and σ2 α,ss of σ2 α that converge in probability to τ 0 and σ 2 0 , respectively.

4 A flexible parametric family of densities on R +

Laguerre approximation of densities

We now consider a particular parametric family of densities, that is rich and flexible and that is based on Laguerre polynomials. The idea is to extend the basic exponential density e -z by using polynomial expansions which are orthogonal to it. These polynomials are the Laguerre polynomials (see e.g. [START_REF] Rice | The Approximation of Functions[END_REF]. Similar ideas have already been used, e.g., by Geerdens

f m (z|θ) = e -z ||θ|| 2 m k=0 θ k v k (z) 2 , (4.1)
where θ = (θ 0 , θ 1 , . . . , θ m ) t , θ is the Euclidean norm of θ, and v k are polynomials orthonormal to the exponential density e -z . We fix θ 0 = 1 to fix the scale of θ. Note that if m = 0, we have f 0 (z|θ) = e -z , and we can check that for any m, f m (z|θ) is a density on R + . Note that in the notation of Section 3, the relation between θ and λ is given by

θ t = (1, λ t ).
The idea of orthogonal polynomials is as follows. The polynomials have to be orthogonal with respect to the inner product

ϕ 1 , ϕ 2 = ∞ 0 ϕ 1 (z)ϕ 2 (z)e -z
dz and they can be obtained by the Gram-Schmidt orthogonalization procedure. Next, they are divided by their norm to get an orthonormal sequence. As pointed out in [START_REF] Geerdens | Goodness-of-fit tests for the frailty distribution in proportional hazards models with shared frailty[END_REF], the interesting feature is that since ∞ 0 e cz e -z dz < ∞ for any 0 < c < 1, the system of polynomials

{v k : k = 0, 1, . . .} is closed to continuous functions ϕ on R + that satisfy ∞ 0 ϕ 2 (z)e -z dz < ∞. (4.2)
It follows thatby defining

c k = ϕ, v k ,we have that lim m→∞ ∞ 0 [ϕ(z) -m k=0 c k v k (z)]
2 e -z dz = 0 for any continuous function ϕ satisfying (4.2). Now consider any continuous density

f Z (•) of Z on R + . By defining ϕ(z) = [f Z (z)/e -z ]
1/2 , we see that (4.2) is verified, so with

c k = ∞ 0 [f Z (z)e -z ] 1/2 v k (z)dz we obtain lim m→∞ ∞ 0 f Z (z) - m k=0 c k √ e -z v k (z) 2 dz = 0,
and hence

lim m→∞ A f Z (z) -f m (z|θ)
potentially approximate any continuous density on R + with respect to the Hellinger distance restricted to the set A.

f m (z|θ) = 2m j=0 a j f γ (z|j + 1, 1), (4.3)
where 2m j=0 a j = 1 and f γ (z|j + 1, 1) =

z j e -z Γ(j + 1)
are the gamma densities with shape parameter j + 1 and scale equal to 1. The a j are real numbers defined as

a j = ||θ|| -2 Γ(j + 1) j k=0 c k c j-k for j = 0, . . . , m ||θ|| -2 Γ(j + 1) m k=j-m c k c j-k for j = m + 1, . . . , 2m , (4.4)
where the coefficients c k are given by

c k = m j=k j k (-1) j-k θ j Γ(k + 1) , for k = 0, . . . , m. (4.5)
Note that equation (4.3) does not define a mixture of Gamma densities, but a linear combination with coefficients a j summing to one but having possibly positive or negative values, as shown by equation (4.5). We refer to Appendix B for the calculation of the coefficients a j in a simulated example.

Thanks to this result we can provide explicit expressions for the moments and an analytical expression for ψ Z (t|θ), its derivative and finally for the function η Z (t|θ). For the moments we have

E(Z k |θ) = 2m j=0 a j (j + 1) k , for k = 1, 2, . . . , (4.6)
2 As pointed out in [START_REF] Geerdens | Goodness-of-fit tests for the frailty distribution in proportional hazards models with shared frailty[END_REF], with an additional assumption on the function ϕ, i.e. [START_REF] Nikiforov | Special Functions of Mathematical Physics[END_REF].

∞ 0 [ϕ (z)] 2 ze -z dz < ∞, the convergence is uniform on any interval [z 1 , z 2 ] ⊂ R + (see
where (j + 1) k is the Pochhammer symbol representing Γ(j + 1 + k)/Γ(j + 1). In the same way we have

ψ Z (t|θ) = 2m j=0 a j 1 1 -it j+1 = 2m j=0 a j 1 + it 1 + t 2 j+1 , (4.7)
therefore we can derive

ψ Z (t|θ) = (-2t + i(1 -t 2 )) (1 + t 2 ) 2 2m j=0 a j j + 1 (1 -it) j . (4.8)
Finally we can compute η Z (t|θ) = ψ Z (t|θ)/ψ Z (t|θ) .

Simulations

For exploring the behavior of our estimator in finite samples, we will compare its performance with that of the estimator suggested in Kneip Our results are displayed in Tables 1 to 4. First of all, a quick scan through the tables confirms that the results are going in the expected directions: as the sample size n increases the RM SE's decrease for a fixed level of ρ nts . This is not uniform in all the cases: see, e.g. Table 3, for ρ nts = 0.25, RM SE( τ α ) is the same for n = 50 and n = 100, but this is because the optimal pair (m, α) is selected for optimizing the sum of the two RM SE's (for estimating τ and σ). But in most cases, the behavior is as expected. Second, for a fixed sample size n, the RM SE increases with ρ nts as expected. Finally we observe that in the 4 scenarios, the level of the accuracy is quite comparable through the different scenarios for the density of Z. Now we will compare the results with the corresponding results from KSVK. We kept the same order for the tables and the same "presentation" to facilitate the comparison. Globally, at least for the 4 preceding cases where the normal assumption is correct, we expect to have less accurate performances, since we only use the symmetry of ε in our method. We summarize in what follows the global picture and the main differences.

Except for the no-noise case (ρ nts = 0), globally our results in terms of the RM SE of the estimators are of the same order as in KSVK. When the size of the noise increases, and when the sample size increases, our results are sometimes better than those obtained in KSVK, although the difference is small. But this is globally good news, because our method does not use the (correct) assumption of normality of ε. The case without noise is much harder to estimate with our approach, which only uses the information on the odd cumulants of order greater than 3. We should also stress that in Table 4, the truncated normal case, where the density of Z is bell-shaped with mode 0.7726 and truncated at zero, our approach based on Laguerre polynomial approximations, works much better than the histogram approach used in KSVK: here the RM SE's are substantially smaller in our approach, in particular for the estimation of τ , in all the cases where ρ nts > 0.

Next, in their Examples 3 and 4, KSVK illustrate how their approach is robust to the normality assumption, by generating samples with ε ∼ C 1 × Stud(4), where C 1 = σ/ √ 2 is a constant to adjust the same noise-to-signal ratios as in the preceding examples. 3 The signal Z is chosen as above as N + (0, 0.8 2 ). The second non-normal specification is ε ∼ C 2 ×Laplace, where again C 2 = σ/ √ 2 to tune correctly ρ nts . Our results are displayed in Tables 5 and6.

In the last two examples we may expect better results from our approach since our specification is correct, whereas the one in KSVK is wrong. The global picture coming from the comparisons with Tables 6 and7 in KSVK is however very similar to the preceding cases: same order of magnitude for the RM SE for the two cases, which confirms that KSVK is rather robust to these two deviations from normality, and that our approach is not quite well adapted to situations where we have no noise.

Estimation of frontiers

We now extend the marginal analysis studied in the previous sections, to the conditional analysis of the boundary of X given a number of explanatory variables. The observable variables are now (Y, W ), where the vector W = (W t 1 , W t 2 ) t contains variables describing the level of outputs in a cost analysis (denoted by W 1 ) and/or environmental variables (denoted by W 2 ). We are interested in the distribution of X, and in particular in the frontier of X, conditional on W 1 ≥ w 1 and/or W 2 = w 2 . In this paper our approach is based on the estimation of τ for given values of w 1 and w 2 and the estimated τα (w 1 , w 2 ) is the estimated frontier function. This means that we estimate τ (w 1 , w 2 ) without any regularity constraint (smoothness, monotonicity, concavity,...) on this frontier function. If the estimator τ (w 1 , w 2 )

does not satisfy the required constraints, τ (w 1 , w 2 ) may in a second step be approximated by a function that verifies the constraints. Another approach not presented in this paper would be to estimate the parameters λ as a function of w 1 and w 2 by imposing constraints.

For instance, λ may be estimated under smoothness constraints with respect to w 1 and w 2 in order to estimate the smooth frontier τ (w 1 , w 2 ).

Let us now give the precise definition of τ (w 1 , w 2 ). Suppose we have a sample (Y i , W i ), i = 1, . . . , n, of i.i.d. observations having the same distribution as (Y, W ). Conditioning on W 1 ≥ w 1 is done in an elementary way : the sample is restricted to the observations verifying W 1j ≥ w 1 (where the inequality is taken componentwise), and the relevant sample size is no longer n but n(w 1 ) = n j=1 I{W 1j ≥ w 1 }, the number of observations for which W 1j ≥ w 1 . In absence of the W 2 variable the construction is hence identical to our previous marginal construction and the asymptotic normality of λα (w 1 ) is derived from the property that

{n(w 1 )} 1/2 ψY (t|W 1 ≥ w 1 ) -ψ Y (t|W 1 ≥ w 1 )
converges to a Gaussian process with covariance function

ψ Y (s + t|W 1 ≥ w 1 ) -ψ Y (s|W 1 ≥ w 1 )ψ Y (t|W 1 ≥ w 1 )
, where ψY (t|W 1 ≥ w 1 ) = n(w 1 ) -1 n j=1 I{W 1j ≥ w 1 }e itY j . Conditioning on W 2 = w 2 requires kernel smoothing or any other method for estimating a conditional expectation. Our marginal approach studied in the previous sections was based on the (penalized) minimization of the distance between η(•)-η(0) and η λ (•)-η λ (0). We now replace η(•) by η(•|W 2 = w 2 ) and the solution λ will be a function of w 2 . Here, η(t|W 2 = w 2 ) equals ∂ ∂t log ψY (t|W 2 = w 2 ), where ψY (t|W 2 = w 2 ) is given by

ψY (t|W 2 = w 2 ) = 1 n n j=1 K h (w 2 -W 2j ) fW 2 (w 2 ) e itY j ,
where

K(u) = k(u 1 ) × • • • × k(u p 2 ) for u = (u 1 , . . . , u p 2 ) with p 2 the dimension of W 2 , k is a univariate kernel function, K h (u) = h -p 2 k(u 1 /h) × • • • × k(u p 2 /h), h is bandwidth sequence, and fW 2 (w 2 ) = n -1 n j=1 K h (w 2 -W 2j
) is a kernel estimator of the density f W 2 of W 2 . Under usual regularity conditions, it may be proved that

(nh p 2 ) 1/2 ψY (•|W 2 = w 2 ) -ψ Y (•|W 2 = w 2 )
converges weakly in E to a zero-mean Gaussian process with covariance operator characterized by the function

k 2 (u)du f W 2 (w 2 ) ψ Y (s + t|W 2 = w 2 ) -ψ Y (s|W 2 = w 2 )ψ Y (t|W 2 = w 2 ) .
Apart from these modifications on the rate and on the covariance, the analysis is identical to the non-conditional case.

Finally, the combination of the conditions W 1 ≥ w 1 and W 2 = w 2 is obvious : it suffices to replace n by n(w 1 ), f W 2 (w 2 ) by f W 2 (w 2 |W 1 ≥ w 1 ) and the conditional characteristic function

is now ψ Y (t|W 1 ≥ w 1 , W 2 = w 2 ).
Based on this result, the asymptotic normality of the estimators λ(w 1 , w 2 ) and τ (w 1 , w 2 )

can now be derived in the same way as it has been done in Theorems 3.2 and 3.3 in the univariate case. In particular, the rates of convergence are the same as for the estimator ψY .

We conclude this section by some Monte-Carlo experiments to illustrate how the methods work and for allowing the comparison with the Kneip et al. (2015) (KSVK) approach.

We thus follow the KSVK scenario and their strategy for the estimation of the frontier function. The only difference is that they use the assumption of normality of the noise in their approach, whereas we only use the information that the noise is symmetric, as developed above. The Monte-Carlo scenario in KSVK (which was already used in [START_REF] Hall | Estimating a changepoint, boundary or frontier in the presence of observation error[END_REF] can be described as follows. We want to estimate a production frontier (instead of a cost frontier), giving the maximum achievable level of the output Y when using the amount of input W . We use the maximum instead of the minimum, in order to be coherent with what was done in KSVK. This has however no influence on the methodology. The stochastic frontier model can thus be written as

Y = τ (W ) exp(-Z) exp(ε), (5.1)
where we only know that Z ≥ 0 and ε is symmetric around zero. As in KSVK, we simulate n data points from the following model: τ (w) = w 1/2 , Z ∼ Exp(3) and ε ∼ N(0, 0.0667 2 ), with W ∼ U(0, 1), all the variables W, Z and ε are independent in this scenario. Note that ρ nts = 0.20. As in KSVK, we analyze the estimation of τ (w) for 3 values of w = (0.25, 0.50, 0.75)

and we select a neighborhood for each w, by using the same bandwidth h w (given by the usual rule of thumb). Note that taking the logarithms we have

log Y = log τ (W ) -Z + ε, (5.2)
which has exactly the additive form of our model. As in KSVK, we do not use directly the observations log Y i for which |W i -w| ≤ h w but we rather use a "local linear" approximation of τ (•) around the point of interest w. To be more precise, the idea is as follows : if τ (•) is sufficiently smooth we have from a Taylor expansion that log τ (W ) is approximately log τ (w) + β t (w)(W -w) and for small h w we then obtain

log Y i ≈ α(w) + β t (w)(W i -w) -Z i (w) + ε i , if |W i -w| ≤ h w , (5.3)
where α(w) = log τ (w) -E(Z|W = w) and Z i (w) = Z i -E(Z|W = w). Since in the last equation the error term has locally zero mean, α(w) and β(w) can be estimated by ordinary least squares. Now we have for all W i with |W i -w| ≤ h w : We used 200 Monte-Carlo replications to evaluate both the bias and the MSE, and compared the results with those obtained in KSVK for the case n = 100 and n = 500. The results are shown in Table 7. The table deserves some comments. First we see that the results are going in the expected direction: the MSE for both τ (w) and σ(w) decreases when n increases. Comparing with the results in KSVK, we see that for n = 100 the results are, as expected, less accurate in our approach, since we only use the symmetry assumption for the distribution of the noise. Note also that we have very few observations for each value of w (on average only 24 data points for log Y i ). As a result, and comparing with Table 8 in KSVK, we see that the MSE for the estimation of τ (w) has an order augmented by a factor 2. For the estimation of σ(w) the situation is worse and we have a factor 10. The situation is much better when the sample size increases to n = 500. Here for the estimation of τ (w) we reach for the MSE the same order as the values in KSVK. For σ(w) we still have disappointing results with the MSE greater by a factor 7. In this case (n = 500), note that for each value of w, we only use on average 88 data points, which is not much for a procedure that does not exploit the normality of the noise, as opposed to KSVK. This is the price to pay for being more general. 4Only for illustration, we show also in Table 8 the results obtained by conditioning not on W = w but rather on W ≤ w as described above (we are in a production setup here so we look for the upper support of Y given that W ≤ w). In this case there is no local linear smoothing step. To summarize, we see that in this scenario the results are similar for the estimation of τ (w), but they are still worse for the estimation of σ(w), where too many Monte-Carlo estimates were set to zero. For the chosen scenario, the first approach seems more appropriate.

log Y i := log Y i -β t (w)(W i -w) ≈ log τ (w) -Z i + ε i ,
To end this section, we illustrate our procedure with a real data set coming The results are displayed in Figure 1 and the 21 values of log w are shown by the blue circles. The final estimate of the cost frontier is obtained by smoothing these pointwise estimates (we used a standard B-splines smoothing technique, cubic splines, 20 knots, penalty on second derivatives with penalty factor 10, see [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF]. We can see how the cost frontier "envelops" the cloud of points by letting as expected a few points outside the frontier. In particular a few points are really extreme (e.g. around the coordinate (11.7,5))

and our model is quite robust to these outliers. In the logarithmic units that we have chosen, the magnitude of σ ε (w) is stable over the 21 pointwise estimates and is of the order of 10 -4 .

values of log(traffic)

Conclusion and further research

This paper presents a new nonparametric analysis of the model Y = τ + Z + ε (with Z ≥ 0 and ε ∈ R), where ε is only constrained to be symmetric, the distribution of Z is supposed to belong to a flexible parametric class, and Z and ε are independent. We also consider the extension to stochastic frontier models. We provide an estimation approach based on the cumulant functions of odd order and we develop asymptotic properties of our estimator. A tractable family for the distribution of the inefficiency Z is proposed, and the power of our method is shown via simulations.

In this paper we estimate the true frontier, i.e. the minimal value of the cost. It is well known that this estimator is not robust, in particular to outliers. We may argue that the introduction of a stochastic error improves the robustness of the estimation of the frontier, but it is possible to adapt the technology of this paper to a robust estimator. From the information contained in the imaginary part of the characteristic function we may estimate the survival function

S(x|w 1 , w 2 ) = P (X ≥ x|W 1 ≥ w 1 , W 2 = w 2 ) of X given W 1 ≥ w 1 and
W 2 = w 2 , which leads to an estimator of the quantile [START_REF] Daouia | Nonparametric efficiency analysis: a multivariate conditional quantile approach[END_REF] and [START_REF] Cazals | Nonparametric frontier estimation: a robust approach[END_REF]. A deeper analysis of the relation between stochastic frontiers and robustness will be the topic of future research.

τ α (w 1 , w 2 ) = S -1 (1 -α|w 1 , w 2 ), or the m-frontier τ m (w 1 , w 2 ) = ∞ 0 S m (u|w 1 , w 2 ) du, see e.g.
Several extensions of deterministic frontiers may be considered in the case of stochastic frontiers with symmetric errors, like endogeneity, separability, among others.

Finally, note that our model is an example of a more general situation. We want to estimate a function or a set of functions which are non identified (different functions are compatible with the data generating process). However, there exists a sieve family of parametric approximations which satisfies the two following properties. First for any finite dimensional parametric approximation the parameters are identified and, second, there exists a sequence of parametric approximations which converges (for a suitable topology) to any function. Intuitively, if the dimension of the parametric model increases the approximation is more accurate but the model is less identified and the choice of the parametric dimension follows from a trade-off between these two arguments. In our example the lack of identification requires the introduction of the penalty term which introduces an estimation bias.

Two tuning parameters should be determined: the size of the approximation and the penalty parameter. In our paper we assume that the true function belongs to a finite parametric family and we do not develop the general theory.

[31] Van der Vaart, A.W. (1998). Asymptotic Statistics. Cambridge University Press, Cambridge.

[32] Van der Vaart, A.W. and J.A. Wellner (1996). Weak Convergence and Empirical Processes. Springer.

A Appendix A: Proofs of the main results

We start by stating the regularity conditions under which our main asymptotic results are valid.

(C1) T λ = η λ -η λ (0) is an operator from Λ to the space E, T λ is injective and Lipschitz continuous, and

sup λ∈Λ |T λ | 2 w < ∞.
(C2) T λ is Fréchet differentiable and the Fréchet derivative T λ is a linear, injective and bounded operator from

R m to E for λ ∈ Λ. The adjoint T * λ of T λ satisfies T * λ -T * λ 0 ≤ γ λ -λ 0 , and T * λ g (g ∈ E)
is continuously differentiable for λ in the interior of Λ.

(C3) There exists a λ α in the interior of Λ that is the unique solution of min λ∈Λ T λ 0 -

T λ 2 + α λ α 2
, for all α > 0.

(C4) The matrix ∆ (defined in Theorem 3.2) is positive definite.

(C5) The function

u : Λ → R 2 : λ → a(λ) b(λ) is continuously differentiable. Proof of Theorem 3.1. It is easily seen that in the space E the process n 1/2 ( ψY (•) -ψ Y (•))
converges weakly to a zero-mean Gaussian process with covariance function for weak convergence in Hilbert spaces). Next, using the Delta-method it follows that the process

(β(s, t)) s,t given by β(s, t) = ψ Y (s + t) -ψ Y (s)ψ Y (t) (see e.
n 1/2 log ψY (•) -log ψ Y (•)
converges weakly to a zero-mean Gaussian process, whose covariance function is given by

1 4 ψ Y (t -s) ψ Y (t)ψ Y (s) + ψ Y (t -s) ψ Y (t)ψ Y (s) - ψ Y (t + s) ψ Y (t)ψ Y (s) - ψ Y (t + s) ψ Y (t) ψ Y (s) ,
since log ψY (t) = (2i) -1 (log ψY (t)-log ψY (t)). Since the covariance function of this process is twice continuously differentiable, it follows from Theorem 2.2.2 in [START_REF] Adler | The Geometry of Random Fields[END_REF] that

n 1/2 T -T λ 0 = n 1/2 ∂ ∂t log ψY (•) -log ψY (0) -log ψ Y (•) + log ψ Y (0)
converges also weakly to a zero-mean Gaussian process whose covariance operator is given in the statement of the theorem. 2

Proof of Theorem 3.2. For the consistency note that we need to show that (see e.g. Van der Vaart (1998), Theorem 5.7, p. 45)

sup λ∈Λ ( T -T λ ) 2 w + α λ 2 -(T λ 0 -T λ ) 2 w + α λ 2 P → 0.
The integral in the above expression is bounded in absolute value by

T 2 w -T 2 λ 0 w + 2 sup λ ( T -T λ 0 )T λ w ≤ T 2 w -T 2 λ 0 w + 2 ( T -T λ 0 ) 2 w 1/2 sup λ |T λ | 2 w 1/2
, and this tends to zero in probability.

Next we prove the asymptotic normality of λα . It follows from (

( λα -λ α ) + T * λα (T λα -T λα ) -(T * λα -T * λα )( T -T λα ) = T * λα ( T -T λ 0 ). 3.2) and (3.3) that α
Note that

(T * λα -T * λα )( T -T λα ) = (T * λα -T * λα ) ( T -T λ 0 ) + (T λ 0 -T λα ) + (T λα -T λα ) = -∇(T * λ S)| S=T λα -T λ 0 ,λ=λα ( λα -λ α ) + o P ( λα -λ α ), since T * λα -T * λα ≤ γ λα -λ α for
some γ < ∞, and similarly for T λα . Hence,

λα -λ α = αI + T * λα T λα + ∇(T * λ S)| S=T λα -T λ 0 ,λ=λα -1 T * λα ( T -T λ 0 ) + o P ( λα -λ α ),
from which the asymptotic normality follows. 2

Proof of Theorem 3.3. By following similar arguments as in the proof of Theorem 3.1, we can show that the process

n 1/2     Y -E(Y ) n -1 n j=1 Y 2 j -E(Y 2 ) ηY (•) -η Y (•)     = n 1/2     -i ∂ ∂t ψY (t)| t=0 -E(Y ) -∂ 2 ∂t 2 ψY (t)| t=0 -E(Y 2 ) ∂ ∂t log ψY (•) -η Y (•)    
converges weakly to a zero-mean Gaussian process with variance operator Ξ characterized for fixed s and t by

     -∂ 2 ∂s∂t β(s, t)| s,t=0 i ∂ 3 ∂s∂t 2 β(s, t)| s,t=0 -1 2 ∂ 2 ∂s∂t β(s,t) ψ Y (t) -β(s,-t) ψ Y (-t) s=0 i ∂ 3 ∂s 2 ∂t β(s, t)| s,t=0 ∂ 4 ∂s 2 ∂t 2 β(s, t)| s,t=0 -1 2i ∂ 3 ∂s 2 ∂t β(s,t) ψ Y (t) -β(s,-t) ψ Y (-t) s=0 -1 2 ∂ 2 ∂s∂t β(s,t) ψ Y (s) -β(-s,t) ψ Y (-s) t=0 -1 2i ∂ 3 ∂s∂t 2 β(s,t) ψ Y (s) -β(-s,t) ψ Y (-s) t=0 σ(s, t)      , where β(s, t) = ψ Y (s + t) -ψ Y (s)ψ Y (t)
c 1 (s) c 2 (s) σ(s, t)     .
Use the notation

A = a 11 a 12 a 21 a 22 ,
let Σ be the operator characterized by σ(s, t), and let C be the operator from E to R 2 given by

Cg = c 1 , g E c 2 , g E .
Then, Ξ may be partitioned into

A C C * Σ ,
where C * is the adjoint operator of C. Define the 2 × 2 matrix

M = 1 -2E(Y ) 0 1 ,
and let Υ be the linear operator from E to R m that is such that the asymptotic variance of λα for fixed α is given by ΥΣΥ * (this is possible thanks to Theorem 3.2). Also, define the

function u : R m → R 2 such that u(λ) = a(λ) b(λ) ,
and let ∂u ∂λ be the 2 × m matrix of partial derivatives. Then, It now follows that

n 1/2     Y -E(Y ) σ2 Y -[E(Y 2 ) -E(Y ) 2 ] u( λα ) -u(λ)     d → N 0, M AM t M CΥ *
n 1/2 τα -τ α σ2 α -σ 2 α d → N 0, B 11 + B 22 -B 12 -B 21 . 2
Proof of Theorem 3.4. We follow the method of proof given in [START_REF] Engl | Regularization of Inverse Problems[END_REF]. First note that by definition of λα ,

T -T λα 2 ≤ T -T λα 2 + α λα 2 ≤ T -T λ 0 2 + α λ 0 2 , (A.1)
and this converges to zero in probability if

α → 0, since T -T λ 0 2 = O P (n -1 ) thanks to Theorem 3.1. Hence, T λα P → T λ 0 , since T λα -T λ 0 ≤ T λα -T + T -T λ 0 P → 0.
Moreover, λα = O P (1), since by (A.1), λα = O P ((nα) -1 ) + λ 0 2 = O P ((nα) -1 ). It now follows that there exists a subsequence λα,ss of λα that converges in probability to a certain limit λ. Hence, T λα,ss P → T λ . We also know that T λα,ss P → T λ 0 . Hence, by the injectivity of T , λ = λ 0 . Since τα and σ2 α are functions of λα , Y and σ2 Y , the same holds true for τα and σ2 α . 2

B Appendix B: Some properties of Laguerre polynomials

The Laguerre polynomials are easy to build. For k = 0, 1, . . ., they can be written as

v k (z) = k j=0 k j (-1) k-j j! z j , (B.1)
and it is easy to prove that

∞ 0 v k 1 (z)v k 2 (z)e -z dz = 0 if k 1 = k 2 1 if k 1 = k 2 . (B.2)
In Figure 2 we illustrate how flexible the family is with only m = 3. We display 5 densities randomly selected from this family. Here, the parameters θ 1 , θ 2 and θ 3 are randomly generated from independent N (0, 1) variables. We see how even with only m = 3 we cover a wide variety of shapes of the resulting densities. The coefficients a 1 , . . . , a 7 are given in Table 9 for each of the 5 densities, and are derived from formulas (4.3)-(4.5).

In Figures 3 and4 we show how usual parametric densities used for positive random variables are easily approximated by a density in our family (4.1). We consider the halfnormal and the truncated-normal. We see that even with low order for m we have already a good 'visual' approximation.

To conclude we give the proof of Lemma 4.1. Proof of Lemma 4.1. Write

f m (z|θ) = e -z θ 2 m k=0 θ k v k (z) 2 = e -z θ 2 m k=0 θ k k j=0 k j (-1) k-j j! z j 2 = e -z θ 2 m j=0 m k=j θ k k j (-1) k-j j! z j 2 .
The expression between braces can be written as m j=0 c j z j with, for j = 0, . . . , m,

c j = m k=j θ k k j
(-1) k-j j! . ,

where a j = b j Γ(j + 1)/ θ 2 . This provides (4.3). Since f m (z|θ) is a density on R + , it integrates to one, yielding 2m j=0 a j = 1. 2

C Appendix C: Some details about the practical implementation

The method has been implemented in MATLAB, but should be easy to implement in R too.

The basic calculations are summarized in Section 3.1. First, the values appearing in the optimization problem (3.1) are calculated for all t ∈ R as

T (t) = n j=1 iY j exp(itY j ) n j=1 exp(itY j ) -Y , T λ (t) = ψ Z (t|θ)/ψ Z (t|θ) -E(Z|θ),
where explicit expressions for ψ Z (t|θ), ψ Z (t|θ) and E(Z|θ) are given in (4. . This is only a rule of thumb based on the fact that n -1/2 is the order of magnitude of the error when estimating ψ Y (t). 5 For the weight function w(t) we used the Epanechnikov kernel rescaled to the interval [-t * , t * ]. 6 For the optimization procedure we used the unconstrained optimization program "fminunc" from MATLAB.

For the estimation of the variance, we slightly modify the formula given in

-(1/2)t 2 σ 2 to ψ ε (t) in a neighborhood of t = 0. Here ψ ε (t) = ψ Y (t) ψ Z (t) exp(it τ α (t)) ,
and then we can define, as suggested in Delaigle and Hall (2016), the estimator

σ 2 α = arg min σ≥0 t 0 -t 0 ψ ε (t) -1 + (1/2)t 2 σ 2 2 dt,
where t 0 is small and fixed as sup |t|≤t 0 | ψ ε (t)-1| ≤ 0.05. In fact the two estimation approaches gave very similar results, but we report in our tables the values obtained by the Delaigle and Hall approach.

  et al. (2013) in a different context (frailty models) for one-parameter gamma densities, and by Bertrand et al. (2017) for approximating the density of a variable with compact support based on Bernstein polynomials. The extended-exponential model is fixed by the chosen order m of the polynomial expansions. The densities in this family are given by

  et al. (2015) (hereafter denoted by KSVK),where an explicit density (a normal with unknown variance) is used to build the penalized likelihood function (see KSVK for details). For the distribution of Z, KSVK use an histogram based density with fixed number of bins. We will use the same scenarios as in their Examples 1 to 4 and we will compare the results in their Tables 2-7 with ours, where we do not use the information that the noise is normal and we use our approach for approximating the density of Z.The scenarios can be summarized as follows. We have the model Y = τ + Z + ε, where τ = 0 and ε ∼ N(0, σ 2 ) with σ = ρ nts σ Z . Various scenarios will be considered for the density of Z, and we take samples of size n = 50, 100, 500 and noise-to-signal ratios equal to ρ nts = 0, 0.10, 0.25, 0.50 and 0.75. In Examples 1a and 1b in KSVK, Z ∼ Exp(β) with µ Z = σ Z = β with β = 2 and β = 1. In their Examples 2a and 2b, Z ∼ N + (β 1 , β 2 2 ) is a truncated normal, where first β 1 = 0 and β 2 = 0.8 giving an half-normal distribution with µ Z = 0.6383 and σ Z = 0.4822, whereas in their Example 2b, β 1 = 0.6 and β 2 = 0.6 leading to a truncated-normal distribution with µ Z = 0.7726 and σ Z = 0.4761, having its mode far from the frontier point. By analogy with KSVK to select the tuning penalty parameter, we select the values of m (the order of the polynomials) and of α (the penalty parameter) by selecting the pair (m, α) that minimizes the Monte-Carlo sum of the two Root Mean Squared Errors (RM SE): RM SE( τ α ) + RM SE( σ α ).

  which suggests to estimate the frontier function at W = w by the upper boundary of the local linear least squares estimator of log Y i for |W i -w| ≤ h w rather than the upper boundary of the local observations log Y i themselves (see the practical details in Kneip et al., 2015).

  from 2326 post offices (delivery offices) of a French operator, observed in 2013. For each post office we have the level of the activity, the output W (the volume of delivered mail), and the labor cost measured by the quantity of labor (number of hours) X. This data set has been used in a different context (deterministic frontier) in Cazals et al. (2016). We are here interested in the cost frontier, conditional on the value W = w, from the sample of n = 2326 observations (W i , Y i ), where i = 1, . . . , n and Y i are the noisy versions of X i . Due to the shape of the cloud of points (quite asymmetric with heavy right tails for both dimensions) we decided to work on a logarithmic scale for both input and output. The localization in the log W is done by selecting a bandwidth h w = 0.1944 which was obtained by least-squares cross validation for estimation of the conditional distribution of log Y given log W (see Li et al, 2013). As explained above for estimating the cost frontier around log w, we use a local linear approximation, as suggested in Kneip et al. (2015). In particular, we do not use directly the observations log Y i for which | log W i -log w| ≤ h w , but their local linear approximations log Y i = log Y i -β t (w)(log W i -log w) where β(w) is the local slope obtained by local ordinary least squares approximation of log Y i over (log W i -log w) (see details above). The estimates are computed for a fixed grid of 21 values of log w. The values for m are selected over a grid of values m = 0, . . . , 5 and for the penalty term over a grid of 21 equidistant values of α ∈ [0, 3]. The optimal values of the pairs (m, α) are obtained, at each selected value of log w, by minimizing the bootstrap estimate of the RMSE of the frontier estimates (500 bootstrap replications at each point).

Figure 1 :

 1 Figure 1: Example with data from post offices, the curve is a smoothed version of the pointwise estimates (blue circles). The n = 2326 data points are also displayed.

  g. Van der Vaart and Wellner, 1996, Chapter 1.8

Figure 2 : 3 .

 23 Figure 2: The 5 densities d1 to d5 are randomly selected from the parametric family with m = 3. The parameters θ 1 , θ 2 and θ 3 are drawn from independent N (0, 1) variables.

Figure 3 :

 3 Figure 3: Plots of the true half-normal density, f Z (z) (dashed-red), together with the best parametric approximation f m (z|θ) (solid blue) for m = 1, 3, 5 (respectively).

Figure 4 :c k c j-k z j = 2m j=0 b j z j , where b j = j k=0 c k c j-k if j = 0, 1 ,

 41 Figure 4: Plots of the true truncated-normal density, f Z (z) (dashed-red), together with the best parametric approximation f m (z|θ) (solid blue) for m = 1, 3, 5, 10 (respectively).

  (3.4) to obtain positive values for σ 2 α = max{0, σ 2 Y -b( λα )}. As suggested in the supplement of Delaigle and Hall (2016) we obtained slightly better results by fitting the polynomial 1

2

  Appendix B explains how to construct Laguerre polynomials and illustrates by means of some simulated examples how flexible the family is and how well continuous functions can be approximated by functions of this family.

	Another nice feature of the Laguerre family is that we can easily express the density in
	(4.1) as a linear combination of gamma densities. We have indeed the following result. The
	proof is given in Appendix B.
	Lemma 4.1 The density given in (4.1) can be written as

  , and where σ(s, t) is defined in Theorem 3.1. Let us

	rewrite this matrix as		a 11	a 12	c 1 (t)
		  	a 21	a 22	c 2 (t)

  ], where t * is selected following the rule motivated in Delaigle and Hall (2016), i.e., t * is the smallest t > 0 such that | ψ Y (t)| ≤ n -1/4

	7), (4.8) and (4.6),
	respectively. The operator is obtained through the build-in function "imag" from MAT-
	LAB.
	The integral in (3.1) is computed by trapezoidal methods over a grid of 2000 equidistant
	points over the interval [-t

* , t *

Table 1 :

 1 Scenario 1: Z ∼ Exp(2) with µ Z = σ Z = 0.5, and Gaussian noise N(0, σ 2 ), see

Table 2 in

 2 [START_REF] Kneip | Frontier estimation in the presence of measurement error with unknown variance[END_REF]. The results are based on 200 Monte-Carlo trials.

		n = 50				n = 100				n = 500
		τ	σ	m	α	τ	σ	m	α	τ	σ	m	α
	Noise to signal ratio ρ nts = 0.00, σ = 0.0000
	RM SE	0.0652 0.0288	4 0.150 0.0523 0.0172	4 0.150 0.0249	0.0101	2 0.050
	BIAS	0.0443 0.0090			0.0405 0.0039			0.0220	0.0024
	ST D	0.0480 0.0274			0.0333 0.0168			0.0117	0.0099
	Noise to signal ratio ρ nts = 0.10, σ = 0.0500
	RM SE	0.0598 0.0517	4 0.100 0.0387 0.0510	4 0.050 0.0160	0.0406	2 0.000
	BIAS	0.0282 -0.0260			0.0076 -0.0114			-0.0000 -0.0118
	ST D	0.0528 0.0448			0.0380 0.0499			0.0160	0.0390
	Noise to signal ratio ρ nts = 0.25, σ = 0.1250
	RM SE	0.0695 0.0749	4 0.040 0.0472 0.0593	4 0.030 0.0262	0.0289	4 0.020
	BIAS	0.0094 -0.0205			0.0020 -0.0085			0.0049 -0.0038
	ST D	0.0690 0.0723			0.0473 0.0588			0.0258	0.0288
	Noise to signal ratio ρ nts = 0.50, σ = 0.2500
	RM SE	0.0833 0.0776	2 0.000 0.0653 0.0607	2 0.000 0.0449	0.0271	4 0.010
	BIAS	0.0117 0.0003			0.0071 0.0199			0.0040 -0.0036
	ST D	0.0827 0.0778			0.0651 0.0575			0.0448	0.0269
	Noise to signal ratio ρ nts = 0.75, σ = 0.3750
	RM SE	0.1116 0.0851	2 0.000 0.0858 0.0570	2 0.000 0.0459	0.0369	2 0.010
	BIAS	0.0268 -0.0270			0.0156 -0.0039			0.0158	0.0021
	ST D	0.1086 0.0810			0.0846 0.0570			0.0432	0.0369

Table 2 :

 2 Scenario 2: Z ∼ Exp(1) with µ Z = σ Z = 1, and Gaussian noise N(0, σ 2 ), see Table3in[START_REF] Kneip | Frontier estimation in the presence of measurement error with unknown variance[END_REF]. The results are based on 200 Monte-Carlo trials.

		n = 50				n = 100				n = 500
		τ	σ	m	α	τ	σ	m	α	τ	σ	m	α
	Noise to signal ratio ρ nts = 0.00, σ = 0.0000
	RM SE	0.0896	0.1687 10 0.800 0.0633	0.1200 10 0.950 0.0247	0.0605 10 0.950
	BIAS	-0.0200 0.1077			-0.0172 0.0789			-0.0085 0.0381
	ST D	0.0875	0.1302			0.0611	0.0906			0.0232	0.0472
	Noise to signal ratio ρ nts = 0.10, σ = 0.1000
	RM SE	0.0903	0.1360 10 0.590 0.0654	0.1025 10 0.980 0.0279	0.0513 10 0.590
	BIAS	-0.0211 0.0341			-0.0156 0.0090			-0.0092 -0.0094
	ST D	0.0881	0.1320			0.0637	0.1024			0.0264	0.0505
	Noise to signal ratio ρ nts = 0.25, σ = 0.2500
	RM SE	0.1005	0.1230 10 0.250 0.0768	0.0887	8 0.400 0.0362	0.0394 10 0.600
	BIAS	-0.0207 0.0132			-0.0205 0.0153			-0.0064 0.0002
	ST D	0.0986	0.1226			0.0742	0.0876			0.0357	0.0395
	Noise to signal ratio ρ nts = 0.50, σ = 0.5000
	RM SE	0.1318	0.1239 10 0.200 0.0983	0.0807 10 0.350 0.0525	0.0434 10 0.600
	BIAS	-0.0038 -0.0177			-0.0029 0.0016			-0.0033 -0.0019
	ST D	0.1321	0.1229			0.0985	0.0808			0.0525	0.0435
	Noise to signal ratio ρ nts = 0.75, σ = 0.7500
	RM SE	0.1604	0.1592	8 0.600 0.1207	0.1013	8 0.600 0.0630	0.0504 10 0.600
	BIAS	0.0088 -0.0432			0.0039 -0.0123			-0.0015 -0.0066
	ST D	0.1605	0.1537			0.1209	0.1008			0.0631	0.0501

Table 3 :

 3 Scenario 3: Z ∼ N + (0, 0.8 2 ) with µ Z = 0.6383 and σ Z = 0.4822, and Gaussian noise N(0, σ 2 ), see Table4in[START_REF] Kneip | Frontier estimation in the presence of measurement error with unknown variance[END_REF]. The results are based on 200 Monte-Carlo trials.

		n = 50				n = 100				n = 500
		τ	σ	m	α	τ	σ	m	α	τ	σ	m	α
	Noise to signal ratio ρ nts = 0.00, σ = 0.0000
	RM SE	0.0689	0.0007	4 0.500 0.0461	0.0001	2 0.550 0.0250	0.0001	2 0.450
	BIAS	0.0221				0.0013	0.0001			0.0028	0.0001
	ST D	0.0655	0.0007			0.0462	0.0000			0.0249	0.0000
	Noise to signal ratio ρ nts = 0.10, σ = 0.0482
	RM SE	0.0662	0.0474	2 0.550 0.0463	0.0480	2 0.550 0.0254	0.0481	2 0.450
	BIAS	-0.0031 -0.0467			0.0023 -0.0479			0.0015 -0.0480
	ST D	0.0663	0.0086			0.0463	0.0027			0.0254	0.0012
	Noise to signal ratio ρ nts = 0.25, σ = 0.1206
	RM SE	0.0837	0.0820	8 0.050 0.0836	0.0501	4 0.125 0.0656	0.0409	4 0.075
	BIAS	0.0338 -0.0259			-0.0556 -0.0181			-0.0491 -0.0161
	ST D	0.0767	0.0780			0.0625	0.0468			0.0436	0.0377
	Noise to signal ratio ρ nts = 0.50, σ = 0.2411
	RM SE	0.1192	0.0832	2 0.100 0.0933	0.0634	4 0.100 0.0730	0.0496 10 0.025
	BIAS	-0.0876 -0.0412			-0.0605 -0.0362			0.0559 -0.0346
	ST D	0.0811	0.0724			0.0712	0.0522			0.0471	0.0356
	Noise to signal ratio ρ nts = 0.75, σ = 0.3617
	RM SE	0.1290	0.0910	2 0.075 0.1015	0.0778	2 0.100 0.0790	0.0573	6 0.050
	BIAS	-0.0884 -0.0441			-0.0772 -0.0464			0.0473 -0.0421
	ST D	0.0942	0.0798			0.0660	0.0626			0.0634	0.0390

Table 4 :

 4 Scenario 4: Z ∼ N + (0.6, 0.6 2 ) with µ Z = 0.7726 and σ Z = 0.4761, and Gaussian noise N(0, σ 2 ), see Table5in[START_REF] Kneip | Frontier estimation in the presence of measurement error with unknown variance[END_REF]. The results are based on 200 Monte-Carlo trials.

		n = 50				n = 100				n = 500
		τ	σ	m	α	τ	σ	m	α	τ	σ	m	α
	Noise to signal ratio ρ nts = 0.00, σ = 0.0000
	RM SE	0.0673	0.0001	6 1.000 0.0494	0.0001	6 1.200 0.0199	0.0002	8 1.400
	BIAS	-0.0089 0.0001			-0.0054 0.0001			0.0010	0.0002
	ST D	0.0669	0.0000			0.0492	0.0000			0.0199	0.0000
	Noise to signal ratio ρ nts = 0.10, σ = 0.0476
	RM SE	0.0681	0.0475	6 1.000 0.0495	0.0475	6 1.200 0.0204	0.0475	6 1.400
	BIAS	-0.0086 -0.0475			-0.0039 -0.0475			-0.0015 -0.0475
	ST D	0.0678	0.0000			0.0494	0.0000			0.0204	0.0000
	Noise to signal ratio ρ nts = 0.25, σ = 0.1190
	RM SE	0.1038	0.0673 10 0.200 0.0877	0.0515 10 0.250 0.0738	0.0222	8 0.250
	BIAS	-0.0712 -0.0228			-0.0695 -0.0193			-0.0700 -0.0068
	ST D	0.0757	0.0635			0.0537	0.0478			0.0233	0.0212
	Noise to signal ratio ρ nts = 0.50, σ = 0.2381
	RM SE	0.1097	0.0949 10 0.150 0.0994	0.0590	8 0.125 0.0891	0.0280	8 0.125
	BIAS	-0.0684 -0.0589			-0.0810 -0.0259			-0.0846 -0.0129
	ST D	0.0859	0.0746			0.0578	0.0531			0.0281	0.0249
	Noise to signal ratio ρ nts = 0.75, σ = 0.3571
	RM SE	0.1250	0.1027	8 0.075 0.0935	0.0815	6 0.100 0.0884	0.0385	6 0.100
	BIAS	-0.0766 -0.0510			-0.0638 -0.0557			-0.0802 -0.0280
	ST D	0.0990	0.0894			0.0685	0.0597			0.0374	0.0266

Table 5 :

 5 Scenario 5: Z ∼ N + (0, 0.8 2 ) with µ Z = 0.6383 and σ Z = 0.4822, and Student noise C 1 × Stud(4), where C 1 = σ/ √ 2, see Table 6 in Kneip et al. (2015). The results are based on 200 Monte-Carlo trials.Noise to signal ratio ρ nts = 0.00, σ = 0.0000

		n = 50				n = 100				n = 500
		τ	σ	m	α	τ	σ	m	α	τ	σ	m	α
	RM SE	0.0689	0.0007	4 0.500 0.0461	0.0001	2 0.550 0.0250	0.0001	2 0.450
	BIAS	0.0221	0.0001			0.0013	0.0001			0.0028	0.0001
	ST D	0.0655	0.0007			0.0462	0.0000			0.0249	0.0000
	Noise to signal ratio ρ nts = 0.10, σ = 0.0482
	RM SE	0.0663	0.0475	2 0.550 0.0468	0.0480	2 0.550 0.0253	0.0481	2 0.450
	BIAS	-0.0032 -0.0467			0.0023 -0.0479			0.0016 -0.0480
	ST D	0.0664	0.0090			0.0469	0.0028			0.0253	0.0012
	Noise to signal ratio ρ nts = 0.25, σ = 0.1206
	RM SE	0.0814	0.0848	8 0.075 0.0790	0.0569	4 0.150 0.0649	0.0414	4 0.075
	BIAS	0.0394 -0.0371			-0.0478 -0.0335			-0.0488 -0.0173
	ST D	0.0714	0.0764			0.0631	0.0460			0.0428	0.0377
	Noise to signal ratio ρ nts = 0.50, σ = 0.2411
	RM SE	0.0967	0.1155	8 0.050 0.0998	0.0643	2 0.100 0.0731	0.0582 10 0.025
	BIAS	0.0412 -0.0564			-0.0837 -0.0374			0.0540 -0.0420
	ST D	0.0876	0.1010			0.0545	0.0525			0.0494	0.0404
	Noise to signal ratio ρ nts = 0.75, σ = 0.3617
	RM SE	0.1282	0.1137	2 0.075 0.1098	0.0807	2 0.075 0.0768	0.0677	6 0.050
	BIAS	-0.0912 -0.0542			-0.0890 -0.0412			0.0376 -0.0511
	ST D	0.0903	0.1002			0.0645	0.0696			0.0671	0.0445
						36

Table 6 :

 6 Scenario 6: Z ∼ N + (0, 0.8 2 ) with µ Z = 0.6383 and σ Z = 0.4822, and Laplace noise C 2 × Laplace, where C 2 = σ/ √ 2, see Table 7 in Kneip et al. (2015). The results are based on 200 Monte-Carlo trials.Noise to signal ratio ρ nts = 0.00, σ = 0.0000

		n = 50				n = 100				n = 500
		τ	σ	m	α	τ	σ	m	α	τ	σ	m	α
	RM SE	0.0689	0.0007	4 0.500 0.0461	0.0001	2 0.550 0.0250	0.0001	2 0.450
	BIAS	0.0221	0.0001			0.0013	0.0001			0.0028	0.0001
	ST D	0.0655	0.0007			0.0462	0.0000			0.0249	0.0000
	Noise to signal ratio ρ nts = 0.10, σ = 0.0482
	RM SE	0.0673	0.0474	2 0.550 0.0460	0.0480	2 0.550 0.0254	0.0481	2 0.450
	BIAS	-0.0031 -0.0468			0.0012 -0.0480			0.0015 -0.0480
	ST D	0.0674	0.0073			0.0461	0.0014			0.0254	0.0011
	Noise to signal ratio ρ nts = 0.25, σ = 0.1206
	RM SE	0.0888	0.0845	8 0.075 0.0805	0.0570	4 0.150 0.0643	0.0374	4 0.075
	BIAS	0.0416 -0.0394			-0.0521 -0.0315			-0.0498 -0.0143
	ST D	0.0787	0.0749			0.0615	0.0477			0.0407	0.0346
	Noise to signal ratio ρ nts = 0.50, σ = 0.2411
	RM SE	0.1285	0.0777	2 0.075 0.0805	0.0831	8 0.050 0.0703	0.0533 10 0.025
	BIAS	-0.0999 -0.0226			0.0521 -0.0438			0.0547 -0.0356
	ST D	0.0810	0.0745			0.0616	0.0708			0.0444	0.0398
	Noise to signal ratio ρ nts = 0.75, σ = 0.3617
	RM SE	0.1319	0.1033	2 0.075 0.0923	0.0944	8 0.025 0.0761	0.0588 10 0.025
	BIAS	-0.0947 -0.0450			0.0527 -0.0414			0.0586 -0.0380
	ST D	0.0920	0.0932			0.0760	0.0851			0.0487	0.0450

Table 7 :

 7 Bivariate example. Same scenario as for Table 8 in Kneip et al. (2015) and conditioning on W = w. The results are based on 200 Monte-Carlo trials.

	n	w = 0.25 w = 0.50 w = 0.75
	100 τ (w) BIAS	0.0027	0.0071	0.0092
	ST D	0.0272	0.0426	0.0489
	M SE	0.0007	0.0019	0.0025
	σ(w) BIAS	-0.0265	-0.0085	-0.0075
	ST D	0.0573	0.0644	0.0708
	M SE	0.0040	0.0042	0.0051
	500 τ (w) BIAS	0.0003	0.0036	0.0054
	ST D	0.0141		0.0246
	M SE	0.0002	0.0004	0.0006
	σ(w) BIAS	-0.0187	0.0019	0.0028
	ST D		0.0492	0.0500
	M SE	0.0022	0.0024	0.0025

Table 8 :

 8 Bivariate example. Same scenario as for Table 8 in Kneip et al. (2015) and conditioning on W ≤ w. The results are based on 200 Monte-Carlo trials.

	n	w = 0.25 w = 0.50 w = 0.75
	100 τ (w) BIAS	-0.0090	-0.0032	-0.0034
	ST D	0.0437	0.0474	0.0478
	M SE	0.0020	0.0023	0.0023
	σ(w) BIAS	-0.0580	-0.0649	-0.0661
	ST D	0.0441	0.0138	0.0061
	M SE	0.0053	0.0044	0.0044
	500 τ (w) BIAS	0.0002	-0.0027	-0.0001
	ST D	0.0220	0.0215	0.0195
	M SE	0.0005	0.0005	0.0004
	σ(w) BIAS	-0.0661	-0.0666	-0.0665
	ST D	0.0067	0.0000	0.0000
	M SE	0.0044	0.0044	0.0044

Table 9 :

 9 Values of the coefficients a 1 , . . . , a 7 for the densities in Figure2.

	density	a 1	a 2	a 3	a 4	a 5	a 6	a 7
	d 1	0.0086 0.0058 0.1519	0.0564 3.8010 -8.3671 5.3434
	d 2	1.4374 -1.3021 -0.1233	1.1649 0.1756 -0.4860 0.1336
	d 3	0.1928 0.3393 0.4643	0.2048 -0.6054 -1.0000 1.4042
	d 4	0.2183 1.0974 1.0840 -11.9649 25.8999 -25.3179 9.9832
	d 5	0.3496 -0.7125 1.3596 -1.9028 1.5836	0.3062 0.0163

Note that the opposite is not true : if X is symmetric indecomposable, it does not necessarily imply that X is identified.

dz = 0,where A is any set in R + for which there exists a δ > 0 such that f Z (z) > δ for all z ∈ A, which indicates that the densities in (4.1) parametrized by the appropriate θ = (θ 0 , θ 1 , . . . , θ m ) t can

Note the unfortunate typo in[START_REF] Kneip | Frontier estimation in the presence of measurement error with unknown variance[END_REF], where C 1 is wrongly defined as σ/2, but this is only a typo, the calculations have been done in KSVK with the correct value of C 1 .

We might expect better performance of our approach compared to KSVK when the noise is not normal.Unfortunately we know from the results in the univariate case, that for the Student-t and the Laplace cases we need a lot of observations to improve the KSVK approach, see the comments above on Tables

and

Here for each value of w, we only have roughly 24 observations when n = 100 and 88 when n = 500.

Note that we tried also with a grid of only 1000 values and a threshold value for | Y (t)| such as n -1/2 , without significant changes in the results.

We tried also other kernels in some pilot experiments, like uniform, triangular, quartic, with similar results. Symmetric kernels giving more weight near zero gave slightly better results. So the results in our tables correspond to the Epanechnikov case.

Financial support from the European Research Council (2016-2021, Horizon 2020 / ERC grant agreement No. 694409), and from IAP Research Network P7/06 of the Belgian State, are gratefully acknowledged.