Eric Gautier
email: eric.gautier@ensae-paristech.fr

Yuichi Kitamura
email: yuichi.kitamura@yale.edu

NONPARAMETRIC ESTIMATION IN RANDOM COEFFICIENTS BINARY CHOICE MODELS

Keywords:

published or not. The documents may come

Nonparametric estimation in random coefficients binary choice models

Introduction

Consider a binary choice model (1.1) Y = I X ′ β ≥ 0 where I denotes the indicator function and X is a d-vector of covariates. We assume that the first element of X is 1, therefore the vector X is of the form X = (1, X′ ) ′ . The vector β is random. The random element (Y, X, β) is defined on some probability space (Ω, F, P), and (y i , xi , β i ), i = 1, ..., N denote its realizations. The econometrician observes (y i , xi ), i = 1, ..., N , but β i , i = 1, ..., N remain unobserved. The vectors X and β correspond to observed and unobserved heterogeneity across agents, respectively. Note that the first element of β in this formulation absorbs the usual scalar stochastic shock term as well as a constant in a standard binary choice model with non-random coefficients.

This formulation is used in [START_REF] Ichimura | Maximum Likelihood Estimation of a Binary Choice Model with Random Coefficients of Unknown Distribution[END_REF], and is convenient for the subsequent development in this paper. Our basic model maintains exogeneity of the covariates X:

Assumption 1.1. β is independent of X, Section 5.3 considers ways to relax this assumption. Under (1.1) and Assumption 1.1, the choice probability function is given by

r(x) = P(Y = 1|X = x) (1.2) = E β [I {x ′ β > 0}].
Discrete choice models with random coefficients are useful in applied research since it is often crucial to incorporate unobserved heterogeneity in modeling the choice behavior of individuals. There is a vast and active literature on this topic. Recent contributions include Briesch, Chintagunta and Matzkin (1996), [START_REF] Brownstone | Forecasting New Product Penetration with Flexible Substitution Patterns[END_REF], [START_REF] Chesher | Taste Variation in Discrete Choice Models[END_REF], [START_REF] Hess | Random Covariance Heterogeneity in Discrete Choice Models[END_REF], Harding and Hausman (2006), [START_REF] Athey | Discrete Choice Models with Multiple Unobserved Choice Characteristics[END_REF], [START_REF] Bajari | Linear Regression Estimation of Discrete choice Models with Nonparametric Distribution of Random Coefficients[END_REF] and [START_REF] Train | Discrete Choice Methods with Simulation[END_REF]. A common approach in estimating random coefficient discrete choice models is to impose parametric distributional assumptions. A leading example is the mixed Logit model, which is discussed in details by [START_REF] Train | Discrete Choice Methods with Simulation[END_REF]. If one does not impose a parametric distributional assumption, the distribution of β itself is the structural parameter of interest. The goal for the econometrician is then to recover it nonparametrically from the information about r(x) obtained from the data.

Nonparametric treatments for unobserved heterogeneity distributions have been considered in the literature for other models. [START_REF] Heckman | A Method for Minimizing the Impact of Distributional Assumptions in Econometric Models for Duration Data[END_REF] study the issue of unobserved heterogeneity distributions in duration models and propose a treatment by a nonparametric maximum likelihood estimator (NPMLE). [START_REF] Elbers | True and Spurious Duration Dependence: The Identifiability of the Proportional Hazard Models[END_REF] also develop some identification results in such models. [START_REF] Beran | Estimating Coefficient Distributions in Random Coefficient Regression[END_REF] and Hoderlein et al. (2007) discuss nonparametric estimation of random coefficients linear regression models. Despite the tremendous importance of random coefficient discrete choice models, as exemplified in the above references, nonparametrics in these models is relatively underdeveloped. In their important paper, [START_REF] Ichimura | Maximum Likelihood Estimation of a Binary Choice Model with Random Coefficients of Unknown Distribution[END_REF] propose an NPMLE for the CDF of β. They present sufficient conditions for identification and prove the consistency of the NPMLE. The NPMLE requires high dimensional numerical maximization and can be computationally intensive even for a moderate sample size. [START_REF] Berry | Nonparametric Identification of Multinomial Choice Demand Models with Heterogeneous Consumers[END_REF] explore nonparametric identification problems in a random coefficients multinomial choice model that often arises in empirical IO. This paper considers nonparametric estimation of the random coefficients distribution, using a novel approach that shares some similarities with standard deconvolution techniques. This allows us to reconsider the identifiability of the model and obtain a constructive identification result. Moreover, we develop a simple plug-in estimator for the density of β that requires no numerical optimization or integration. It is easy to implement in empirical applications, while being flexible about the treatment of unobserved heterogeneity.

Since the scale of β is not identified in the binary choice model, we normalize it so that β is a vector of Euclidean norm 1 in R d . The vector β then belongs to the d -1 dimensional sphere S d-1 . This is not a restriction as long as the probability that β is equal to 0 is 0. Also, since only the angle between X and β matters in the binary decision I{X ′ β ≥ 0}, we can replace X by X/ X without any loss of information. We therefore assume that X is on the sphere S d-1 as well in the subsequent analysis. Results from the directional data literature are thus relevant to our analysis. We aim to recover the joint probability density function f β of β with respect to the uniform spherical measure σ over S d-1 from the random sample (y 1 , x 1 ), . . . , (y N , x N ) of (Y, X).

The problem considered here is a linear ill-posed inverse problem. We can write (1.3) r(x) = b∈S d-1

I x ′ b ≥ 0 f β (b)dσ(b) = H(x) f β (b)dσ(b) := H (f β ) (x)
where the set H(x) is the hemisphere {b : x ′ b ≥ 0}. The mapping H is called the hemispherical transformation. Inversion of this mapping was first studied by [START_REF] Funk | Uber Eine Geometrische Anwendung der Abelschen Integralgleichung[END_REF] and later by [START_REF] Rubin | Inversion and Characterization of the Hemispherical Transform[END_REF]. [START_REF] Groemer | Geometric Applications of Fourier Series and Spherical Harmonics[END_REF] also discusses some of its properties. H is not injective without further restrictions and conditions need to be imposed to ensure identification of f β from r. Even under a set of assumptions that guarantee identification, however, the inverse of H is not a continuous mapping, making the problem ill-posed. To see this, suppose we restrict f β to be in L 2 (S d-1 ). Since the kernel of H is square integrable by compactness of the sphere, it is Hilbert-Schmidt and thus compact.

Therefore if the inverse of H were continuous, H -1 H would map the closed unit ball in L 2 (S d-1 ) to a compact set. But the Riesz theorem states that the unit ball is relatively compact if and only if the vector space has finite dimension. The fact that L 2 (S d-1 ) is an infinite dimensional space contradicts this. Therefore the inverse of H cannot be continuous. In order to overcome this problem, we use a one parameter family of regularized inverses that are continuous and converge to the inverse when the parameter goes to infinity. This is a common approach to ill-posed inverse problems in statistics (see, e.g. [START_REF] Carrasco | Linear Inverse Problems in Structural Econometrics Estimation Based on Spectral Decomposition and Regularization[END_REF].

Due to the particular form of its kernel that involves the scalar product x ′ b, the operator H is an analogue of convolution in R d , as illustrated in a simple example in Section A.1.1 of Supplemental Appendix. This analogy provides a clear insight into the identification issue. In particular, our problem is closely related to the so-called boxcar deconvolution (see, e.g. [START_REF] Groeneboom | Density Estimation in the Uniform Deconvolution Model[END_REF] and [START_REF] Johnstone | Periodic Boxcar Deconvolution and Diophantine Approximation[END_REF]), where identifiability is often a significant problem. The connection with deconvolution is also useful in deriving an estimator based on a series expansion on the Fourier basis on S 1 or its extension to higher dimensional spheres called Fourier-Laplace series.

These bases are defined via the Laplacian on the sphere, and they diagonalize the operator H on L 2 S d-1 . Such techniques are used in [START_REF] Healy | An Empirical Bayes Approach to Directional Data and Efficient Computation on the Sphere[END_REF] for nonparametric empirical Bayes estimation in the case of the sphere S 2 . The kernel of the integral operator H, however, does not satisfy the assumptions made by Healy and Kim. Unlike [START_REF] Healy | An Empirical Bayes Approach to Directional Data and Efficient Computation on the Sphere[END_REF], we make use of socalled "condensed" harmonic expansions. The approach replaces a full expansion on a Fourier-Laplace basis by an expansion in terms of the projections on the finite dimensional eigenspaces of the Laplacian on the sphere. This is useful since an explicit expression of the kernel of the projector is available.

It enables us to work in any dimension and does not require a parametrization by hyperspherical coordinates nor the actual knowledge of an orthonormal basis. This approach, to the best of our knowledge, appears to be new in the econometrics literature.

The paper is organized as follows. Section 2 provides a practical guide for our procedure, which is easy to implement. Section 3 deals with identification while introducing basic notions used throughout the paper. We derive the convergence rates of the estimators in all the L q spaces for q ∈ [1, ∞] and also prove a pointwise CLT in Section 4. Some extensions, such as estimation of marginals, treatments of models with non-random coefficients, and the case with endogenous regressors are presented in Section 5. Simulation results are reported in Section 6. Section 7 concludes. Supplemental

Appendix presents analysis of a toy model, technical tools used in the main text, estimators for choice probabilities that are used to construct our density estimators, and the proofs of the main results.

A Brief Guide for Practical Implementation

This section presents our basic estimation procedure when a random sample {(y i , xi )} generated from the model (1.1) is available. As noted in Section 1, normalize covariates data and define

x i = (1, x′ i )/ (1, x′ i ) ∈ S d-1 , i = 1, ..., N .
To estimate the joint density of the random vector β, use the following formula:

(2.1) fβ (b) = max   2 |S d-1 | T N -1 p=0 χ(2p + 1, 2T N )h(2p + 1, d) λ(2p + 1, d)C ν(d) 2p+1 (1)   1 N N i=1 (2y i -1)C ν(d) 2p+1 (x ′ i b) max fX (x i ), m N   , 0   .
The factors |S d-1 |, χ, h and λ are constants that do not depend on data and trivial to compute.

The surface area S d-1 of S d-1 is given by S

d-1 = 2π d/2
Γ(d/2) where Γ denotes the Gamma function. The constants h, ν and λ are obtained via the numerical formulas h(n, d)

= (2n+d-2)(n+d-2)! n!(d-2)!(n+d-2) , ν(d) = (d -2)/2 and λ(2p + 1, d) = (-1) p |S d-2 |1•3•••(2p-1) (d-1)(d+1)•••(d+2p-1)
, respectively. The function χ is defined on N × N and used for smoothing. This is to be chosen by the user: see Proposition A.3 as well as the numerical example reported in Section 6 for examples of χ. The truncation parameter T N needs to be chosen so that it grows with the sample size with a sufficiently slow rate. The trimming factor m N is also user-defined, and it is chosen so that it goes to zero as the sample size increases. The notation C ν n (•) signifies the Gegenbauer polynomial 1 ; They, for example, correspond to the Chebychev polynomials of the first kind in the case of one random slope (i.e. the case with d = 2) 2 . The only remaining factor which needs to be calculated in the above formula is the nonparametric density estimator fX for f X on S d-1 . For example, the following nonparametric estimator can be used:

(2.2) fX (x) = max   1 |S d-1 | T ′ N n=0 χ(n, T ′ N )h(n, d) C ν(d) n (1) 1 N N i=1 C ν(d) n (x ′ i x) , 0  
where T ′ N is an another truncation parameter, playing a role similar to T N . Our estimator fβ requires neither numerical integration nor optimization. This is a clear advantage over existing estimators for random coefficient binary choice models, including many parametric estimators. This is our main proposal, on which the rest of the paper focuses. In Section 4 we explain how the formula (2.1) is derived, and investigate its asymptotic properties. 1 The Gegenbauer polynomials are given by

C ν n (t) = [n/2] l=0 (-1) l (ν) n-l l!(n -2l)! (2t) n-2l , ν > -1/2, n ∈ N
where (a)0 = 1 and for n in N \ {0}, (a)n = a(a + 1)

• • • (a + n -1) = Γ(a + n)/Γ(a)
. See Section A.1.2 for further properties of the Gegenbauer polynomials.

2 When d = 2, the following relations can be used in (2.1) and (2.2)

∀p ≥ 0, 1 |S d-1 | h(2p + 1, 2)C 0 2p+1 (x ′ i b) λ(2p + 1, 2)C 0 2p+1 (1) = (-1) p (2p + 1) π cos (2p + 1) arccos(x ′ i b) , ∀n ≥ 0, 1 |S d-1 | h(n, 2)C 0 n (x ′ i b) C 0 n (1) = 1 π cos n arccos(x ′ i b) .

Identification Analysis

In this section we address the following two questions:

(Q1) Under what conditions is f β identified? (Q2) Does the random coefficients model impose restrictions?

To answer these questions it is useful to introduce the notion of the odd and even part of a function defined on the sphere.

Definition 3.1. We denote the odd part and the even part of a function f by

f -(b) = (f (b) -f (-b))/2 and f + (b) = (f (b) + f (-b))/2, respectively, for every b in S d-1 .
Let us start with the question (Q1). As noted in Section A.1.4, operating H reduces the even part of a function to a constant 1 and therefore it is impossible to recover f + β from the knowledge of r, which is what observations offer. Our identification strategy is therefore as follows: (Step 1) Assume conditions that guarantee the identification of f - β ; then (Step 2) Show that f β is uniquely determined from f - β under a reasonable assumption. We first consider Step 1. Define

H + = H(n) = {x ∈ S d-1 : x ′ n ≥ 0}
, where n = (1, , 0, ..., 0) ′ , that is, the northern hemisphere of S d-1 . For later use, also define its southern hemisphere H -= H(-n). Since the model we consider has a constant as the first element of the covariate vector before normalization, the same vector after normalization is necessarily an element of H + . We make the following assumption, which also appears in [START_REF] Ichimura | Maximum Likelihood Estimation of a Binary Choice Model with Random Coefficients of Unknown Distribution[END_REF], and show that it achieves Step 1.

Assumption 3.1. The support of X is H + .
This assumption demands that X, the vector of non-constant covariates in the original scale, is supported on the whole space R d-1 . It rules out discrete or bounded covariates; see Section 5 for a potential approach to deal with regressors with limited support. In what follows we assume that the law of X is absolutely continuous with respect to σ and denote its density by f X . Step 1 of our identification argument is to show that the knowledge of r(x) on H + , which is available under Assumption 3.1, identifies f - β . The problem at hand calls for solving r = Hf β = 1 2 + Hf - β for f - β , and the inversion formula derived in (4.1) is potentially useful for the purpose. A direct application of the formula to r is inappropriate, however, since it requires integration of r on the whole sphere S d-1 , but r is defined only on H + even when X has full support on R d-1 . An appropriate extension of r(x), x ∈ H + to the entire S d-1 is in order. Using the random coefficients model (1.1) and Assumption 1.1, then noting that f β is a probability density function, conclude

(3.1) H(f β )(-x) = H(-x) f β (b)dσ(b) = 1 -H(f β )(x) = 1 -r(x)
for x in H + . This suggests an extension R of r to S d-1 as follows:

(3.2) ∀x ∈ H + , R(x) = r(x), and ∀x ∈ H -, R(x) = 1 -r(-x) = 1 -R(-x).
The function R is well-defined on the whole sphere under Assumption 3.1. Later we derive a formula for f - β in terms of R(x), x ∈ S d-1 , which shows the identifiability of f - β under Assumption 3.1. Note that

R(x) = R + (x) + R -(x) (3.3) = 1 2 [R(x) + R(-x)] + R -(x) = 1 2 [R(x) + (1 -R(x))] + R -(x) by (3.2) = 1 2 + R -(x)
thus R is completely determined by its odd part and therefore,

R(x) = 1 2 + H f - β (x), or (3.4) R -= Hf - β .
We can invert this equation to obtain f - β . Now we turn to Step 2 in our identification argument. Obviously f - β does not uniquely determine f β without further assumptions. This is a fundamental identification problem in our model.

We need to identify f β from the choice probability function r, but we can choose an appropriate even function g so that f β + g is a legitimate density function (see the proof of Proposition 3.1 for such a construction). Then r = H (f β + g), and the knowledge of r identifies f β only up to such a function g.

Ichimura and Thompson (1998, Theorem 1) give a set of conditions that imply the identification of the model (1.1). One of their assumptions postulates that there exists c on S d-1 such that P(c ′ β > 0) = 1. This, in our terminology, means that: Assumption 3.2. The support of β is a subset of some hemisphere.

As noted by [START_REF] Ichimura | Maximum Likelihood Estimation of a Binary Choice Model with Random Coefficients of Unknown Distribution[END_REF], Assumption 3.2 does not seem too stringent in many economic applications. It is often reasonable to assume that an element of the random coefficients vector, such as a price coefficient, has a known sign. If the j-th element of β has a known sign (and positive), then Assumption 3.2 holds with c being a unit vector with its j-th element being 1. This is a case in which the location of the hemisphere in Assumption 3.2 is known a priori, though the knowledge about its location is not necessary for identification. Assumption 3.2 implies the following mapping from f - β to f β developed in (A.24):

(3.5) f β (b) = 2f - β (b)I f - β (b) > 0 .
This is useful because it shows that Assumption 3.2 guarantees identification if f - β is identified. Moreover, it will be used in the next section to develop a key formula that leads to a simple and practical estimator for f β that is guaranteed to be non-negative. Remark 3.1. Assumption 3.2 is testable since it imposes restrictions on f - β , which is identified under weak conditions. For example, for values of b with f - β (b) > 0, f - β (-b) < 0 must hold. Or, it implies that f - β integrates to 1/(2|S d-1 |) on a hemisphere H(x) for some x, and -1/(2|S d-1 |) on the other H(-x).

The subsequent result, Proposition 3.1, answers question (Q2), and a proof is given in Supplemental Appendix.

Notation. We use the notation L 2 (S d-1 ) for the space of square integrable complex valued functions equipped with the hermitian product (f, g) L 2 (S d-1 ) = S d-1 f (x)g(x)dσ(x), and more generally use L p (S d-1 ) for p ∈ [1, ∞] the Banach space of p-integrable functions and • p the corresponding norm.

We also use the notation W s p (S d-1 ) (and H s (S d-1 ) for p = 2) to signify the corresponding Sobolev spaces with norm • p,s defined as The global smoothness assumption that R belongs to H d/2 (S d-1 ) imposes substantial restriction on the property of observables, that is, the behavior of the choice probability function r. Note that the smoothness condition in this proposition is stated in terms of R, and even if the choice probability function r is sufficiently smooth on the support of X, which is H + , it is not necessarily consistent with the random coefficients binary choice model (1.1) unless its extension is smooth globally on S d-1 . In particular, the Sobolev embedding of H s (S d-1 ) into the space of continuous functions for s > (d -1)/2 implies that if the extension R is in H d/2 (S d-1 ), it has to be continuous on S d-1 . This, in turn, means that the corresponding r has to satisfy certain matching conditions at a boundary point x of H + (i.e.

f p,s = f p + -∆ S s
x ′ n = 0) and its opposite point -x.

4. Nonparametric Estimation of f β 4.1. Derivation of the closed form estimation formula. This section discusses how the closed form estimation formula (2.1) is derived. Suppose an odd function f -defined on S d-1 satisfies an integral equation f -= Hg with g square integrable with respect to the spherical measure. In Section A.1.4 we show that the solution to this equation is given by: (4.1)

H -1 (f -)(y) = ∞ p=0 1 λ(2p + 1, d) S d-1 q 2p+1,d (x, y)f -(x)dσ(x)
where expressions for λ and q are provided in Proposition A.4 and Theorem A.1, respectively. If an appropriate estimator Rof R -is available, an application of the inversion formula (4.1) to (3.4) suggests the following estimator for f - β :

f - β = H -1 R- (4.2) = ∞ p=0 1 λ(2p + 1, d) S d-1 q 2p+1,d (•, x) R-(x)dσ(x).
Then use the mapping (3.5) to define

(4.3) fβ (b) = 2 f - β (b)I f - β (b) > 0
as an estimator for f β .

We use the following notation in the rest of the paper:

Notation. For two sequences of positive numbers (a n ) n∈N and (b n ) n∈N , we write a n ≍ b n when there exists a positive M such that M -1 b n ≤ a n ≤ M b n for every positive n.

Proposition A.6 implies that if f - β -f - β ∈ H s (S d-1 ) then R--R -∈ H σ (S d-1 ), σ = s + d 2 and for v ∈ [0, s], (4.4) f - β -f - β 2,v ≍ R--R - 2,v+d/2 .
As discussed earlier, the estimation of f β is related to deconvolution in S d-1 , and the degree of illposedness in our model is d/2, which is indeed the rate at which the absolute values of the eigenvalues of H (c.f. Proposition A.4) λ(n, d), n = 2p + 1, p ∈ N converges to zero as p grows, as shown in (A.27).

Existing results for deconvolution problems (see, for example, [START_REF] Fan | On the Optimal Rates of Convergence for Nonparametric Deconvolution Problems[END_REF]Kim and[START_REF] Kim | Directional Mixture Models and Optimal Esimation of the Mixing Density[END_REF] then suggest that we should be able to estimate f β at the rate

N -s 2s+2d-1 in the L 2 (S d-1 ) provided that f β ∈ H s (S d-1
). The relationship (4.4), evaluated at v = 0, implies that this can be achieved if

we can estimate R -at the rate N -σ-d 2 2σ+d-1 in the • 2,d/2
norm. The latter is the usual nonparametric rate for estimation of densities on d -1 dimensional smooth submanifolds of R d (see, for example, [START_REF] Hendriks | Nonparametric Estimation of a Probability Density on a Riemannian Manifold Using Fourier Expansions[END_REF].

The estimation formula given in (4.2) is natural and reasonable, though it typically requires numerical evaluation of integrals to implement it. Moreover, in practice one needs to evaluate the infinite sum in (4.2), for example, by truncating the series. This results in a general estimator that can be written in the following two equivalent forms

f - β = H -1 P TN R- (4.5) = T N p=0 1 λ(2p + 1, d) S d-1 q 2p+1,d (•, x) R-(x)dσ(x)
for suitably chosen TN that goes to infinity with N and P TN defined in (A.20). The sequence H -1 • P TN , N = 1, 2, ... can be interpreted as regularized inverses of H, with the spectral cut-off method often used in statistical inverse problems.

We now discuss how to obtain Rin the calculation of (4.5). The following choice is particularly convenient:

(4.6) R-(x) = 1 N N i=1 (2y i -1)K - 2T N (x i , x) max fX (x i ), m N
where m N is a trimming factor going to 0 with the sample size, K -(x i , •) denotes the odd part (of the second argument) of the kernel function K(x i , •) defined in (A.23) and fX is a nonparametric density estimator for f X . See Section A.1.5 of Supplemental Appendix for the derivation of the above formula. Various nonparametric estimators for f X can be used in (4.6), since estimation of densities on compact manifolds have been studied by several authors, using histogram [START_REF] Ruymgaart | Strong Uniform Convergence of Density Estimators on Spheres[END_REF]), projection estimators (see, e.g. [START_REF] Devroye | Nonparametric Density Estimation. The L1-View[END_REF] for the circle and [START_REF] Hendriks | Nonparametric Estimation of a Probability Density on a Riemannian Manifold Using Fourier Expansions[END_REF] for general compact Riemannian manifolds) or kernel estimators (see, e.g. [START_REF] Devroye | Nonparametric Density Estimation. The L1-View[END_REF] for the case of the circle, and [START_REF] Hall | Kernel Density Estimation with Spherical Data[END_REF] and [START_REF] Klemelä | Estimation of Densities and Derivatives with Directional Data[END_REF] for higher dimensional spheres).

Note also that [START_REF] Baldi | Adaptive Density Estimation for Directional Data Using Needlets[END_REF] develops an adaptive density estimator on the sphere using needlet thresholding. In the simulation experiment we use

(4.7) fX (x) = max 1 N N i=1 K T ′ N (x i , x), 0
for a suitably chosen T ′ N that depends on the sample size and the smoothness of f X and K T ′ N is a kernel of the form (A.23) satisfying Assumption A.1. Note that its rate of convergence in supnorm can be obtained in the same manner as the proof of Theorem 4.1. This estimator is in the spirit of the projection estimators of [START_REF] Hendriks | Nonparametric Estimation of a Probability Density on a Riemannian Manifold Using Fourier Expansions[END_REF], but here we are able to derive a closed form using the condensed harmonic expansions together with the Addition Formula. Note also that K T N is a smoothed projection kernel (note the factor χ in (A.23)), which is used here in order to have good approximation properties in the L q (S d-1 ) norms with arbitrary q ∈ [1, ∞], in particular in the L ∞ (S d-1 ) norm. Using (4.5) and (4.6) with TN = 2T N , define

f - β = H -1 R-= H -1   1 N N i=1 (2y i -1)K - 2T N (x i , •) max fX (x i ), m N   .
Computing fβ is straightforward. First, note that the estimator (4.6) for R -resides in a finite dimensional space T N p=0 H 2p+1,d , therefore P 2T N R-= Rholds. Consequently, unlike in (4.5) where a general estimator for R -is considered, we do not need to apply any additional series truncation to Rprior to the inversion of H. Second, the estimator requires no numerical integration. To see this, note the formula

H -1 K - 2T N (x i , •) (b) = T N -1 p=0 χ(2p + 1, 2T N ) λ(2p + 1, d) q 2p+1,d (x i , b),
which follows from

S d-1 q 2p+1,d (x, b)K - 2T N (x, x i )dσ(x) = S d-1 q 2p+1 (x, b) T N -1 p ′ =1 χ(2p ′ + 1, 2T N )q 2p ′ +1,d (x, x i )dσ(x) = χ(2p + 1, 2T N )q 2p+1,d (b, x i ).
which, in turn, can be seen by the definition of K T in (A.23), the fact that the integral operators with q as kernels are projections and (A.16). Thus

f - β (b) = 1 N N i=1 2y i -1 max fX (x i ), m N T N -1 p=0 χ(2p + 1, 2T N ) λ(2p + 1, d) q 2p+1,d (x i , b).
Using (4.3) and the Addition formula (Theorem A.1), we arrive at an estimator for f β with the following explicit form:

fβ (b) = 2 f - β (b)I{ f - β (b) > 0}, (4.8) where f - β (b) = 1 |S d-1 | T N -1 p=0 χ(2p + 1, 2T N )h(2p + 1, d) λ(2p + 1, d)C ν(d) 2p+1 (1)   1 N N i=1 (2y i -1)C ν(d) 2p+1 (x ′ i b) max fX (x i ), m N   .
This is equivalent to the formula (2.1) previously presented in Section 2. Likewise, using the definition of the smoothing kernel (A.23) and the Addition Theorem in the above definition (4.7) of fX , we obtain the formula (2.2) as well.

4.2.

Rates of Convergence in L q (S d-1 )-norms. Now we analyze the rate of our estimator fβ . The following assumption is weak and reasonable.

Assumption 4.1. f X ∈ L ∞ .
The proofs of the following theorems and corollaries in the rest of this section are given in Section A.1.6 of Supplemental Appendix. Theorem 4.1 (Upper bounds in L q (S d-1 )). Suppose Assumptions A.1, 3.1 and 4.1 hold, and choose T N that does not grow more than polynomially fast in N . If f - β belongs to W s q (S d-1 ) with q in [1, ∞] and s > 0, and

(4.9) max i=1,...,N f X (x i ) -fX (x i ) = O p (m N ) , then, for any 1 ≤ r ≤ q, fβ -f β q = O p m -1 N N -1/2 T (2d-1)/2 N (log N ) (1/2-1/q)I{q≥2} + T -s N + T d/2 N m -2 N max i=1,...,N f X (x i ) -fX (x i ) T d/2+(d-1)(1-1/r) n σ (f X < m N ) 1/q-1/r+1 . (4.10)
When there exists m > 0 such that f X ≥ m σ a.e. on H + , the following holds for the estimator without the trimming factor (i.e. m N = 0) when the estimator fX which is consistent in sup norm:

fβ -f β q = O p N -1/2 T (2d-1)/2 N (log N ) (1/2-1/q)I{q≥2} +T -s N + T d/2 N max i=1,...,N f X (x i ) -fX (x i ) .
The first term in (4.10) is the stochastic error, the second term is the approximation bias, the third the plug-in error and the fourth the trimming bias. Note that Theorem 4.1 imposes the mild assumption (4.9); otherwise, we need to replace

T d/2 N m -2 N max i=1,...,N f X (x i ) -fX (x i ) in (4.10) with T d/2 N m -2 N max i=1,...,N f X (x i ) -fX (x i ) 1 + (log N ) (1/2-1/q)I{q≥2} N -1/2 T (d-1)/2 N . Since max i=1,...,N f X (x i ) -fX (x i ) ≤ f X -fX ∞ ,
this term can be made of order

O P N log N -v/(2v+d-1) when f X ∈ W v
∞ with a suitably chosen parameter T ′ N if we take (4.7) as an estimator. The proof of the latter statement is classical and can be obtained simplifying the proof of Theorem 4.1 and Corollary 4.1. Equation (4.10) yields that, for proper choices of m N going to zero and T N to infinity, fβ is consistent given that f X has some smoothness in the Sobolev scales.

Though the additional condition of f X being bounded away from 0 in the last statement of Theorem 4.1 is convenient, it is restrictive. To see this, consider the d = 2 case. In polar coordinates,

f X (cos(θ), sin(θ)) = f X (tan(θ))(1 + tan 2 (θ), thus, assuming f X ≥ m on H + , which does not require trimming, yields ∀x ∈ R, f X (x) ≥ m 1 + x 2 .
It implies that X has tails larger than Cauchy tails and all moments are infinite. The introduction of the trimming factor m N allows us to relax the assumption f X ≥ m, though it introduces bias. As is clear from (4.10), the condition for the trimming bias to go to zero with N depends both on T N and m N . The quantity σ(f X < m N ) should decay to zero with N sufficiently fast. We can check, for example, that when X is standard Gaussian then σ(f

X < m N ) = O (-log m N ) -1/2 , when it is Laplace then σ(f X < m N ) = O (-log m N ) -1 and when f X is proportional to (1 + x 2 ) -k with k > 1 we obtain that σ(f X < m N ) = O m 1/(2(k-1)) N
. In all these cases, it is possible to adjust adequately T N and m N and to obtain rates of convergence. The upper bound on the rates become slower as the tail of f X becomes thinner. Let us now return to the general case of d -1 regressors. The assumptions below allow us to obtain rates that differ slightly from the rates that we would obtain in the ideal case where f X ≥ m σ a.e. for positive m on H + . Assumption 4.2. Suppose for q in [1, ∞], there exist positive τ and r X such that

(i ) σ(f X < h) = O(h τ ) and f X ∈ L ∞ ,
and either

(ii ) max i=1,...,N f X (x i ) -fX (x i ) = O p N (log N ) (1-2/q)I{q≥2} -r X or, (iii ) for some constant C, lim N →∞ N (log N ) (1-2/q)I{q≥2} r X max i=1,...,N f X (x i ) -fX (x i ) ≤ C a.s.

holds.

As seen before, Assumption 4.2 (ii) or (iii) are very mild. (i) holds for a reasonable class of distributions for f X . In the above example where f X is proportional to (1 + x 2 ) -k with k > 1, we have the relation τ = ρ/(2(k -1)). This allows for a higher order moment to exist for a large k. Corollary 4.1. Assume that f - β belongs to W s q (S d-1 ) with q in [1, ∞] and s > 0. Let assumptions A.1, 3.1, 4.1 and 4.2 (i) and (ii) hold, and take

m N ≍ N (log N ) (1-2/q)I{q≥2} -ρ , T N ≍ N (log N ) (1-2/q)I{q≥2} γ(ρ)
where ρ yields a maximum γ of

γ(ρ) = min 1 -2ρ 2s + 2d -1 , 2ρτ 2s + d + 2(d -1)(1 -1/q) , 2r X -4ρ 2s + d , 1 d -1 .
We then have

(4.11) fβ -f β q = O p N (log N ) (1-2/q)I{q≥2} -γs .
Moreover, if, instead of Assumption 4.2 (ii), Assumption 4.2 (iii) holds with q = ∞, then there exists a constant C such that

(4.12) lim N →∞ N log N γs fβ -f β ∞ ≤ C a.s.
The rate γs in Corollary 4.1 accounts for the dimension d -1, the degree of smoothing d/2 of the operator and features of the density of the covariates (i.e. its smoothness and tail behavior).

We now make stronger assumptions on f X and its estimate that yield, up to a logarithmic term, the convergence rate N -s 2s+2d-1 . We need to be able to trim the estimate of f X with a term which is logarithmic in N : m N = (log N ) -ρ for some positive ρ.

Assumption 4.3. Suppose for q in [1, ∞], and positive r σ and r X ,

(i ) σ(f X < (log N ) -ρ ) = O N (log N ) 2ρ+(1-2/q)I{q≥2}
-rσ , and either

(ii ) max i=1,...,N f X (x i ) -fX (x i ) = O p (log N ) -2ρ N (log N ) 2ρ+(1-2/q)I{q≥2} -r X or, (iii ) for some constant C, lim N →∞ (log N ) 2ρ N (log N ) 2ρ+(1-2/q)I{q≥2} r X max i=1,...,N f X (x i ) -fX (x i ) ≤ C a.s. Corollary 4.2. Assume that f - β belongs to W s q (S d-1
) with q in [1, ∞] and s > 0. Let assumptions A.1, 3.1, 4.1 and 4.3 (i)-(ii), hold, and take

T N ≍ N (log N ) 2ρ+(1-2/q)I{q≥2} γ where γ = min 1 2s + 2d -1 , 2r σ 2s + d + 2(d -1)(1 -1/q) , 2r X 2s + d then we have (4.13) fβ -f β q = O p N (log N ) 2ρ+(1-2/q)I{q≥2} -γs .
Moreover, if, instead of Assumption 4.3 (ii), Assumption (iii) holds with q = ∞, then there exists a constant C such that

(4.14) lim N →∞ N (log N ) 2ρ+1 γs fβ -f β ∞ ≤ C a.s. When f X ∈ W s+d/2+ǫ ∞ (S d-1
) for any positive ǫ then 2r X 2s+d > 1 2s+2d-1 and γ in Corollary 4.2 is simply min

1 2s+2d-1 , 2rσ 2s+d+2(d-1)(1-1/q) . Recall that the smoothness s + d/2 is related to the smoothness of R. Indeed, we have seen in Section 3 that R ∈ W s+d/2 2 (S d-1 ) if and only if f β ∈ W s 2 (S d-1
). Consider now the most restrictive case where f X ≥ m σ a.e., then the estimator without the trimming factor (i.e. m N = 0) satisfies the following:

Corollary 4.3. Assume that f - β belongs to W s q (S d-1
) with q in [1, ∞] and s > 0. Let assumptions A.1, 3.1 and 4.1 hold, and suppose, for positive r X ,

(4.15) max i=1,...,N f X (x i ) -fX (x i ) = O p N (log N ) (1-2/q)I{q≥2} -r X . Take T N ≍ N (log N ) (1-2/q)I{q≥2} γ where γ = min 1 2s + 2d -1 , 2r X 2s + d then we have (4.16) fβ -f β q = O p N (log N ) (1-2/q)I{q≥2} -γs .
Moreover, if we replace (4.15) by for some positive C (4.17)

N (log N ) (1-2/q)I{q≥2} r X max i=1,...,N f X (x i ) -fX (x i ) ≤ C a.s. then (4.18) lim N →∞ N log N γs fβ -f β ∞ ≤ C a.s. When f X belongs to W s-d/2+ǫ ∞
, for arbitrary positive ǫ, γ = [START_REF] Healy | An Empirical Bayes Approach to Directional Data and Efficient Computation on the Sphere[END_REF] who study deconvolution on S 2 for non-degenerate kernels. [START_REF] Kim | Directional Mixture Models and Optimal Esimation of the Mixing Density[END_REF] prove that the rate in Healy and [START_REF] Healy | An Empirical Bayes Approach to Directional Data and Efficient Computation on the Sphere[END_REF] is optimal in the minimax sense. Their statistical problem, however, involves neither a plug-in method nor trimming. Also, somewhat less importantly, it does not cover the case when the convolution kernel is given by an indicator function, which appears in our operator H. Hoderlein et al. ( 2010) study a linear model of the form W = X ′ β where β is a d-vector of random coefficients. They obtain a nonparametric random coefficients density estimator that has the L 2 -rate N -s 2s+2d-1 when f X ≥ mσ a.e. for positive m3 when f X is assumed to be bounded from below and thus no trimming is required. They also consider trimming but the approach is slightly different and rates of convergence are not given. Unlike the previous results, we cover L q loss for all q ∈ [1, ∞].

4.3.

Pointwise Asymptotic Normality. This section discusses the asymptotic normality property of our estimator.

N 1/2 T -(d-1)/2 N m -2 N max i=1,...,N f X (x i ) -fX (x i ) = o p (1), (4.19) N -1/2 T (d-1)/2 N m -(1+ǫ) N = o(1) for some ǫ > 0, (4.20) N 1/2 T -2s+2d-1 2 N = o(1), (4.21) N 1/2 T (d-1)/2 N σ ({f X < m N }) = o(1) (4.22) then (4.23) N 1 2 s -1 N (b) fβ (b) -f β (b) d → N (0, 1) holds for b such that f β (b) = 0, where s 2 N (b) := var(Z N (b)), Z N (b) = 2 (2Y -1)H -1 K - 2T N (X,•) (b) max(f X (X),m N )
.

The standard error s N (b) is the standard deviation of (4.24)

Z N (b) = 2 |S d-1 | T N -1 p=0 χ(2p + 1, 2T N )h(2p + 1, d) λ(2p + 1, d)C ν(d) 2p+1 (1) (2Y -1)C ν(d) 2p+1 (X ′ b) max(f X (X), m N )
(see equation (4.8)), which can be estimated using an estimate fX of f X .

The next theorem is concerned with the restrictive case where the density of the covariates is bounded from below and hence the trimming factor m N is set at zero.

Theorem 4.3 (Asymptotic normality when the density of the covariates is bounded from below).

Suppose f - β belongs to W s ∞ (S d-1 ) with s > 0, and Assumptions A.1, 3.1 and 4.1 hold. If fX , f X and T N satisfy

N 1/2 T -(d-1)/2 N max i=1,...,N f X (x i ) -fX (x i ) = o p (1), (4.25) N -1/2 T (d-1)/2 N = o(1) (4.26) N 1/2 T -2s+2d-1 2 N = o(1), (4.27) then (4.28) N 1 2 s -1 N (b) fβ (b) -f β (b) d → N (0, 1)
holds for b such that f β (b) = 0, where s 2 N (b) := var(Z N (b)), Z N (b) = 2 (2Y -1)H -1 K - 2T N (X,•) (b) f X (X)
.

A formula for Z N for this case is obtained by replacing max(f X (X), m N ) with f X (X) in (4.24).

Discussion

S d-1 g(b)dσ d-1 (b) = B k S d-1-k g ρ b u, b dσ d-1-k (u) dP * σ d-1 b where B k is the k dimensional ball of radius 1. Setting g = |S d-1 |f β (b)I b ∈ A for A Borel set of B k
shows that the marginal density of β with respect to the dominating measure P * σ d-1 is given by (5.2)

f β b = |S d-1 | S d-1-k f β ρ b u, b dσ d-1-k (u).
One can use deterministic methods to compute the integral (e.g., [START_REF] Narcowich | Decomposition of Besov and TriebelLizorkin spaces on the sphere[END_REF] for quadrature methods on the sphere) or for example one may use a Monte-Carlo method, by forming

(5.3) f M β b = |S d-1 | M M j=1 fβ ρ b u j , b
where u j , j = 1, ..., M are draws from independent uniform random variables on S d-1-k .

5.2.

Treatment of Non-Random Coefficients. It may be useful to develop an extension of the method described in the previous sections to models that have non-random coefficients, at least for two reasons. 4 First, the convergence rate of our estimator of the joint density of β slows down as the dimension d of β grows, which is a manifestation of the curse of dimensionality. Treating some coefficients as fixed parameters alleviates this problem. Second, our identification assumption in Section 3 precludes covariates with discrete or bounded support. This may not be desirable as many random coefficient discrete choice models in economics involve dummy variables as covariates. As we shall see shortly, identification is possible in a model where the coefficients on covariates with limited support are non-random, provided that at least one of the covariates with "large support" has a non-random coefficient as well. More precisely, consider the model:

(5.4)

Y i = I{β 1i + β ′ 2i X 2i + α 1 Z 1i + α ′ 2 Z 2i ≥ 0}
where

β 1 ∈ R and β 2 ∈ R d X -1 are random coefficients, whereas the coefficients α 1 ∈ R and α 2 ∈ R d Z -1 are nonrandom. The covariate vector (Z 1 , Z ′ 2 ) ′ is in R d Z , though the (d Z -1)
-subvector Z 2 might have limited support: for example, it can be a vector of dummies. The covariate vector (X ′ 2 , Z 1 ) ′ is assumed to be, among other things, continuously distributed. Normalizing the coefficients vector and the vector of covariates to be elements of the unit sphere works well for the development of our procedure, as we have seen in the previous sections. The model (5.4), however, is presented "in the original scale" to avoid confusion. 4 Hoderlein et al. (2010) suggest a method to deal with non-random coefficients in their treatment of random coefficient linear regression models.

Define β

* 1 (Z 2 ) := β 1 + α ′ 2 Z 2 .
We also use the notation

τ (Z 2 ) := (β * 1 (Z 2 ), α 1 , β 2 ) ′ (β * 1 (Z 2 ), α 1 , β ′ 2 ) ∈ S d X +1 , W := (1, Z 1 , X ′ 2 ) ′ (1, Z 1 , X ′ 2 ) ′ ∈ S d X +1 . Then (5.4) is equivalent to: Y = I{(β * 1 (Z 2 ), α 1 , β 2 )(1, Z 1 , X ′ 2 ) ′ ≥ 0} = I τ (Z 2 ) ′ W ≥ 0 .
This has the same form as our original model if we condition on Z 2 = z 2 . We can then apply previous results for identification and estimation under the following assumptions. First, suppose (β 1 , β ′ 2 ) ′ and W are independent, instead of Assumption 1.1. Second, we impose some conditions on f W |Z 2 =z 2 , the conditional density of W given Z 2 = z 2 . More specifically, suppose there exists a set Z 2 ⊂ R d Z -1 , such that Assumption 3.1 holds if we replace f X and d with f W |Z 2 =z 2 and d X + 1 for all z 2 ∈ Z 2 . If Z 2 is a vector of dummies, for example, Z 2 would be a discrete set. By (A.30) and (4.1) we obtain

(5.5) f - τ (Z 2 )|Z 2 =z 2 (t) = ∞ p=0 1 λ(2p + 1, d X + 1) E (2Y -1)q 2p+1,d X +1 (W, t) f W |Z 2 =z 2 (W ) Z 2 = z 2
for all z 2 ∈ Z 2 , where the right hand side consists of observables. This determines f τ (Z 2 )|Z 2 =z 2 . That is, the conditional density

f (β * 1 (Z 2 ), α 1 , β 2 ) (β * 1 (Z 2 ), α 1 , β 2 ) ′ Z 2 = z 2 is
identified for all z 2 ∈ Z 2 (Here and henceforth we use the notation f (•|•) to denote conditional densities with appropriate arguments when adding subscripts is too cumbersome). This obviously identifies (5.6)

f (β * 1 (Z 2 ), α 1 , β 2 ) β 2 Z 2 = z 2
for all z 2 ∈ Z 2 as well. If we are only interested in the joint distribution of β 2 under a suitable normalization, we can stop here. The presence of the term α 1 Z 1 in (5.4) is unimportant so far. Some more work is necessary, however, if one is interested in the joint distribution of the coefficients on all the regressors. Notice that the distribution (5.6) gives

f β * 1 (Z 2 ) β 2 Z 2 = z 2 = f β 1 + α ′ 2 Z 2 β 2 Z 2 = z 2 ,
from which we can, for example, get

E β * 1 (Z 2 ) β 2 Z 2 = z 2 = E β 1 β 2 + E 1 β 2 α ′ 2 z 2 for all z 2 ∈ Z 2 .
Define a constant

c := E 1 β 2
then we can identify cα 2 as far as z 2 ∈ Z 2 has enough variation and

E α 1 β 2 = cα 1 is identified as well. Let (5.7) f (β ′ 2 , α 1 , α ′ 2 ) ′ β 2
denote the joint density of all the coefficient (except for β 1 , which corresponds to the conventional disturbance term in the original model (5.4), normalized by the length of β 2 ). Then

f (β ′ 2 , α 1 , α ′ 2 ) ′ β 2 = f           I d X -1 0 0 1 . . . cα 2 cα 1        β 2 β 2 α 1 β 2        .
In the expression on the right hand side, f ((β ′ 2 , α 1 ) ′ / β 2 ) is available from (5.6), and cα 1 and cα 2 are identified already, therefore the desired joint density (5.7) is identified. Obviously (5.7) also determines the joint density of (β ′ 2 , α 1 , α ′ 2 ) ′ under other suitable normalizations as well. The density (5.5) is estimable: when Z 2 is discrete, one can use the estimator of Section 4 to each subsample corresponding to each value of Z 2 . If Z 2 is continuous we can estimate f W |Z 2 =z 2 and the conditional expectation by nonparametric smoothing. An estimator for the density (5.6) can be then obtained numerically. 5.3. Endogenous Regressors. Assumption 1.1 is violated if some of the regressors are endogenous in the sense that the random coefficients and the covariates are not independent. This problem can be solved if an appropriate vector of instruments is available. To be more specific, suppose we observe (Y, X, Z) generated from the following model

(5.8) Y = I{β 1 + β′ X ≥ 0} with (5.9) X = ΓZ + V
where V is a vector of reduced form residuals and Z is independent of (β, V ). Note that Hoderlein et al. (2010) utilize a linear structure of the form (5.9) in estimating a random coefficient linear model.

The equations (5.8) and (5.9) yield

Y = I{ β 1 + V ′ β + Z ′ Γ ′ β}.
Suppose the distribution of ΓZ satisfy Assumption 3.1. It is then possible to estimate the density of τ = τ / τ where τ = β 1 + V ′ β, β ′ by replacing Γ with a consistent estimator, which is easy to obtain under the maintained assumptions. This yields an estimator for the joint density of β/ τ , the random coefficients on the covariates under scale normalization.

Numerical Examples

The purpose of this section is to illustrate the performance of our new estimator in finite samples using simulated data. We consider the model of the form (1.1) with d = 3. The covariates are specified to be X = (1, X 1 , X 2 ) where (X 2 , X 3 ) ′ ∼ N ( 0 0 , 2•I 2 ). The coefficients vector β = (β 1 , β 2 , 1) ′ is set random except for the last element. Fixing the last component constant fulfills Assumption 3.2 for identification. Two specifications for the random elements (β 1 , β 2 ) are considered. In the first specification (Model 1) we let (β 1 , β 2 ) ′ ∼ N ( 0 0 , 0.3 • I 2 ). In the second (Model 2) we consider a two point mixture of normals ) is defined on S 2 in this case, and we performed an appropriate transformation to plot it as a density on R 2 . In the case of model 1, with the reasonable sample size, the location of the peak of the density, as well as its shape, are successfully recovered by our procedure. For model 2, again, our procedure works well: the estimated surface plot nicely captures the locations of the two peaks and their shapes of the true density, thereby exhibiting the underlying mixture structure. While further experimentations are necessary, these results seem to indicate our estimator's good performance in practical settings.

β 1 β 2 ∼ λN   µ -µ ,   σ 2 ρσ 2 ρσ 2 σ 2     + (1 -λ)N   -µ µ ,   σ 2 ρσ 2

Conclusion

In this paper we have considered nonparametric estimation of a random coefficients binary choice model. By exploiting (previously unnoticed) connections between the model and statistical deconvolution problems and applying results of integral transformation on the sphere, we have developed a new estimator that is practical and possesses desirable statistical properties. It requires neither numerical optimization nor numerical integration, and as such its computational cost is trivial and local maxima and other difficulties in optimization need not be of concern. Its rate of convergence in the L q norm for all q ∈ [1, ∞] is derived. Our numerical example suggests that the new procedure works well in finite samples, consistent with its good theoretical properties. It is of great theoretical and practical interest to obtain an adaptive procedure for choosing the smoothing parameters of our estimator, though it is a task we defer to subsequent investigations. 5 With appropriate under-smoothing, the estimator is shown to be asymptotically normal, providing a theoretical basis for nonparametric statistical inference for the random coefficients distribution. 5 Gautier and Le Pennec (2011) consider a needlet-based procedure and discuss its rate optimality in a minimax sense and adaptation. As is often the case when Fourier series techniques are used, we consider spaces of complex valued functions. Let L p (S 1 ) denote the Banach space of Lebesgue p-integrable functions and its norm by

• p .
In the case of L 2 (S 1 ), the norm is derived from the hermitian product 2π 0 f (θ)g(θ)dθ. Let R θ and f φ denote the extension R of r according to (3.2) and f β after the reparameterization. Our task is then to obtain f φ from the knowledge of R θ . Rewrite (1.3) using these definitions, then divide both sides by π, to get:

(A.1) R θ π (θ) = H(f β ) π (θ) = 2π 0 1 π I {|θ -φ| < π/2} f φ (φ)dφ.
If we further define f θ := R θ /π and f η (η) := 1 π I{|η| < π/2}, then using the standard notation for convolution, (A.1) can be written as f θ = f η * f φ . It is now obvious that the estimation of f φ (thus f β ) is linked to the following statistical deconvolution problem: unobservable random variables φ and η with densities f φ and f η are related to an observable random variable θ according to θ = η + φ, and one wishes to recover f φ from f θ , the density of θ, when f η is known (and it is Uniform[-π/2, π/2] in this case). 1 The problem of deconvolution on the unit circle can be conveniently solved using Fourier series. The set of functions exp(-int)/ √ 2π n∈Z is the orthonormal basis of L 2 (S 1 ) used to define Fourier series. This system is also complete in L 1 (S 1 ). Reparameterize a function f ∈ L 1 (S 1 ) it using angles as above, and denote it by f t . Denoting the Fourier coefficients of f ∈ L 1 (S 1 ) by 1 It is also useful to note that the inversion of H is closely related to differentiation. Differentiating the right hand-side of expression (A.1) with respect to θ identifies f φ (θ + π/2) -f φ (θ -π/2) where f φ is defined on the line by periodicity.

If f φ is supported on a semicircle, with an assumption that is elaborated further in Section 3, f φ (which is positive) is identified. Thus if the model is identified the inverse of H is a differential operator and as such unbounded.

c n (f t ) = 2π 0 f t (t) exp(-int)dt/(2π), (A.2) f t (θ) = n∈Z c n (f t ) exp(int)
holds in the L 1 (S 1 ) sense. Recall also that for f and g in L 1 (S 1 ), after the same reparameterization,

(A.3) c n (f t * g t ) = 2πc n (f t )c n (g t ).
Using equation (A.3) we obtain the following proposition.

Proposition A.1. c 0 (R θ ) = πc 0 (f φ ) and for n ∈ Z \ {0}, c n (R θ ) = c n (f φ ) 2 sin (nπ/2) /n.
As in classical deconvolution problems on the real line, our aim is to obtain f t (thus f β ) using equation (A.2) and Proposition A.1. Proposition A.1 shows that c 2p (R θ ) = 0 holds for all non-zero p's, regardless of the values of c 2p (f φ ), p ∈ Z\{0}. Thus from r(x) = R θ (θ) one can only recover the Fourier coefficients c n (f φ ) for n = 0 (which is easily seen to be 1/2π, by integrating both sides of (A.1)

and noting that f β is a probability density function) and n = 2p + 1, p ∈ Z. The same phenomenon occurs in higher dimensions, as explained in Section A.1.3.

Remark A.1. The vector spaces H 2p+1,2 = span exp(i(2p + 1)t)/ √ 2π, exp(-i(2p + 1)t)/ √ 2π , p ∈ N are eigenspaces of the compact self-adjoint operator H on L 2 (S 1 ). These eigenspaces are associated with the eigenvalues 2(-1) p 2p+1 . Also, p∈N\{0} H 2p,2 is the null space ker H.

A.1.2. The Gegenbauer polynomials. We summarize some results on the Gegenbauer polynomials, which are used in various parts of the paper. These can be found in Erdélyi et al. (1953) and [START_REF] Groemer | Geometric Applications of Fourier Series and Spherical Harmonics[END_REF]. When ν = 0 and d = 2, it is related to the Chebychev polynomials of the first kind, as

∀n ∈ N \ {0}, C 0 n (t) = 2 n T n (t)
and

C 0 0 (t) = T 0 (t) = 1 hold for T n (t) = cos (n arccos(t)) , n ∈ N.
When ν = 1 and d = 4, C 1 n (t) coincides with the Chebychev polynomial of the second kind U n (t), which is given by

U n (t) = sin[(n + 1)

arccos(t)] sin[arccos(t)]

, n ∈ N.

sense:

f = ∞ n=0 h(n,d) l=1 (f, Y n,l ) L 2 (S d-1 ) Y n,l .
Thus h(n, d) is the multiplicity of the eigenvalue ζ n,d , and H n,d is the corresponding eigenspace.

Lemma A.1 (i), (ii) and (iv) give the decomposition

L 2 (S d-1 ) = n∈N H n,d .
The space of surface harmonics of degree 0 is the one dimensional space spanned by 1. A series expansion on an orthonormal basis of surface harmonics is called a Fourier series when d = 2, a

Laplace series when d = 3 and in the general case a Fourier-Laplace series.

Orthonormal bases of surface harmonics usually involve parametrization by angles, such as the spherical coordinates when d = 3 as used by [START_REF] Healy | An Empirical Bayes Approach to Directional Data and Efficient Computation on the Sphere[END_REF] or hyperspherical coordinates for d > 3. Instead, here we work with the decomposition of a function on the spaces H n,d as presented in the next definition so that we avoid specific expressions of basis functions.

Definition A.2. The condensed harmonic expansion of a function

f in L 1 (S d-1 ) is the series ∞ n=0 Q n,d f , where Q n,d is the projector from L 2 (S d-1 ) to H n,d .
This leads to a simple method both in terms of theoretical developments and practical implementations. The projector Q n,d can be expressed as an integral operator with kernel (A. [START_REF],2 of the Bateman Manuscript Project[END_REF])

q n,d (x, y) = h(n,d) l=1 Y n,l (x)Y n,l (y), where (Y n,l ) h(n,d) l=1
is any orthonormal basis of H n,d . The kernel has a simple expression given by the addition formula:

Theorem A.1 (Addition Formula). For every x and y ∈ S d-1 , we have

(A.17) q n,d (x, y) = ♭ q n,d (x ′ y), ♭ q n,d (t) := h(n, d)C ν(d) n (t) |S d-1 |C ν(d) n (1) where C ν n are Gegenbauer a and ν(d) = (d -2)/2.
The Sobolev spaces are defined in the Fourier-Laplace domain through the fractional Laplacian defined on a certain subset of L p (S d-1 ) as

(A.18) -∆ S s/2 f := ∞ n=0 ζ s/2 n,d Q n,d f.
For the case where p = 2, in stead of the definition of the norm • p,s given in Section 3 it is also possible to use an equivalent norm, the square of which is equal to

∞ n=0 (1 + ζ n,d ) s Q n,d f 2 2 .
The following integration by parts holds for functions f in H 1 (S d-1 ) (A. [START_REF] Fisher | Statistical Analysis of Spherical Data[END_REF])

- S d-1 f (x)∆ S f (x)dσ(x) = S d-1 ∇ S x f ′ ∇ S x f dσ(x)
and as a consequence for the second definition of the norm of H 1 (S d-1 ) we have

f 2 2,1 = f 2 2 + ∇ S f 2 2 .
In Section A.1.1 we observed the close relationship between the random coefficients binary choice model and convolution for d = 2. This connection remains valid in higher dimensions. Suppose a function f (x, y) defined on S d-1 ⊗ S d-1 depends on x and y only through the spherical distance d(x, y) = arccos(x ′ y) (that is, f is a zonal function). Consider the following integral:

h(x) = S d-1 f (x, y)g(y)dσ(y) := f * g(x),
then the function h is a convolution on the sphere. We now see that the choice probability function

r(x) = H(f β )(x) = S d-1 I{x ′ b ≥ 0}f β (b)dσ(b
) is a special case of h and therefore can also be regarded as convolution. Obtaining f β from r (or, inverting H) is therefore a deconvolution problem.

In what follows we often write f (x, ⋆) when a function f on S d-1 ⊗ S d-1 is regarded as a function of ⋆. Also, the notation f (x, ⋆) p is used for the L p norm of f (x, ⋆), that is, f (x, ⋆) p = S d-1 |f (x, y)| p dσ(y). Note that if f is a zonal function as in the above definition of spherical convolution, its L p norm f (x, ⋆) p does not depend on x. The following Young inequalities for convolution on the sphere (see, for example, [START_REF] Kamzolov | The Best Approximation of the Class of Functions W α p (S n ) by Polynomials in Spherical Harmonics[END_REF]) are useful: Proposition A.2 (Young inequalities). Suppose f (x, ⋆) and g belong to L r (S d-1 ) and L p (S d-1 ), respectively. Then h(x) = f * g(x) is well-defined in L q (S d-1 ) and h q ≤ f r g p , where 1 ≤ p, q, r ≤ ∞ and 1 q = 1 p + 1 r -1.

Let P T denote the projection operator onto T n=0 H n,d , i.e.

(A.20)

P T f (x) = T n=0 Q n,d f (x) = S d-1 D T (x, y)f (y)dσ(y) Assumption A.1. (i ) K T (x, ⋆) 1 is uniformly bounded in T .
(ii ) There exists constants C and α such that for all x, y, z ∈ S d-1 ,

|K T (z, x) -K T (z, y)| ≤ C x -y T α ,
where • denotes the Euclidean norm.

(iii ) For p ∈ [1, ∞] and s > 0, there exists a constant C such that for every f in W s p (S d-1 ),

f (•) - S d-1 K T (•, y)f (y)dσ(y) p ≤ CT -s f p,s .
(iv ) χ(•, T ) takes values in [0, 1] and is such that there exists c > 0 such that for all 0 ≤ n ≤ ⌊T /2⌋,

χ (n, T ) ≥ c.
The smoothed projection kernel K T (x, y) depends on x and y only through d(x, y), thus the value of the norm K T (x, ⋆) 1 in Assumption (i) does not depend on x ∈ S d-1 . Assumption (i) could be relaxed, but imposing this on K T allows us to make relatively weak assumptions on the smoothness of the density of the covariates later in this paper. Assumption (ii) is used to establish the L ∞ -rates of convergence of our estimators. Assumption (iii) provides bounds for approximation errors. Under

this condition, K T * f approximates f ∈ L p (S d-1
) with an error of the same order as that of the best n-th degree spherical harmonic approximation of a function f ∈ L p (S d-1 ) in W s p (S d-1 ) (see e.g. [START_REF] Kamzolov | The Best Approximation of the Class of Functions W α p (S n ) by Polynomials in Spherical Harmonics[END_REF][START_REF] Ditzian | Fractional Derivatives and Best Approximation[END_REF]. This is useful in our treatment of the bias terms in our estimators.

As concrete examples, the following two choices for the weight function χ in (A.23) satisfy Assumption A.1, as shown in the appendix. The first and the second choices of χ correspond to the Riesz kernel and the delayed means kernel, respectively. Proposition A.3. In the definition of the smoothed kernel (A.23), let

χ(n, T ) = 1 - ζ n,d ζ T,d + 1 s/2 l , where l is an integer satisfying l > (d -2)/2, or χ(n, T ) = ψ(n/T ) where ψ : [0, ∞) → [0, ∞) is infinitely differentiable, nonincreasing, such that ψ(x) = 1 if x ∈ [0, 1], 0 ≤ Ψ(x) ≤ 1 if x ∈ [1, 2], ψ(x) = 0 if x ≥ 2. Then K T satisfies Assumption A.1.
The delayed means kernel has the nice property that it does not require prior knowledge of the regularity s in Assumption A.1. The Dirichlet kernel satisfies (ii), (iii) (for p = 2) and (iv) of Remark A.2. If f is a probability density function, the coefficient of degree 0 in the expansion of f on surface harmonics is 1/|S d-1 |. Conversely, any harmonic polynomial or series such that its degree

0 coefficient is 1/|S d-1 | integrates to one.
The next theorem shows that Fourier-Laplace series on the sphere is a natural tool for the study of the operator H.

Theorem A.2 (Funk-Hecke Theorem). If g belongs to H n,d for some n, and a function F on (-1, 1)

satisfies 1 -1 |F (t)| 2 (1 -t 2 ) (d-3)/2 dt < ∞, then (A.25) S d-1 F (x ′ y)g(y)dσ(y) = λ n (F )g(x)
where

λ n (F ) = |S d-2 |C ν(d) n (1) -1 1 -1
F (t)C ν(d) n (t)(1 -t 2 ) d-3 2 dt.
In other words, the kernel operator defined by

f ∈ L 2 (S d-1 ) → x → S d-1 F (x ′ y)f (y)dσ(y) ∈ L 2 (S d-1 )
is, in the subspace H n,d , equivalent to the multiplication by λ n (F ). Thus a basis of surface harmonics diagonalizes an integral operator if its kernel is a function of the scalar product x ′ y.

Remark A.3. Healy and [START_REF] Healy | An Empirical Bayes Approach to Directional Data and Efficient Computation on the Sphere[END_REF] use Fourier-Laplace expansions to analyze a deconvolution problem on S 2 . As we shall see below, the Addition Formula along with condensed harmonic expansions provide a general treatment that works for arbitrary dimensions.

(i ) λ(0, d) = |S d-1 | 2 (ii ) λ(1, d) = |S d-2 | d-1 A-11 (iii ) ∀p ∈ N \ {0}, λ(2p, d) = 0 (iv ) ∀p ∈ N, λ(2p + 1, d) = (-1) p |S d-2 |1•3•••(2p-1) (d-1)(d+1)•••(d+2p-1) .
Proof of Proposition A.

C ν(d) n (t)(1 -t 2 ) (d-3)/2 dt, thus using (A.7), α(n, d) = (-2) -n (d -2) n n! ((d -1)/2) n 1 0 d n dt n (1 -t 2 ) n+(d-3)/2 dt.
Therefore for n ≥ 1 and d ≥ 3,

α(n, d) = - (-2) -n (d -2) n n! ((d -1)/2) n d n-1 dt n-1 (1 -t 2 ) n-1+(d-3)/2 dt t=0
since the term on the right hand-side is equal to 0 for t = 1. To prove that the coefficients α(2p, d)

are equal to zero for p positive it is enough to prove

d 2p+1 dt 2p+1 (1 -t 2 ) 2p+1+m t=0 = 0, ∀m ≥ 1, p ≥ 0.
The Faá di Bruno formula gives that this quantity is equal to

k 1 +2k 2 =2p+1
(-1) 2p+1-k 2 (2p + 1)!(m + 1)

• • • (2p + 1 + m) k 1 !k 2 ! (1 -t 2 ) m+k 2 (2t) k 1 t=0
.

and the result follows since k 1 in the sum cannot be equal to 0.

When n = 2p+1 for p ∈ N we obtain, again using the Faá di Bruno formula, that the derivative at t = 0 is equal to

(-1) p (2p)! p! [(2p + 1 + (d -3)/2)(2p + (d -3)/2) • • • (p + 2 + (d -3)/2)] .
Together with (A.9), the desired result follows. For the case d = 2 we use Proposition A.1.

Define L 2 odd (S d-1 ) and H s odd (S d-1 ) as the restrictions of L 2 (S d-1 ) and H s (S d-1 ) to odd functions and similarly L 2 even (S d-1 ) and H s even (S d-1 ) for even functions. The following corollary is a direct consequence of the Funk-Hecke Theorem and Proposition A. Corollary A.1. The null space of the hemispherical transform H is given by

ker H = ∞ p=1 H 2p,d = f ∈ L 2
even (S d-1 ) :

S d-1 f (x)dσ(x) = 0 ,
when H is viewed as an operator on L 2 (S d-1 ). The spaces H 0,d and H 2p+1,d for p ∈ N are the eigenspaces associated with the non-zero eigenvalues of H.

As a consequence of Proposition A.4, H is not injective and restrictions have to be imposed in order to ensure identification of f β . Section 3 presents sufficient conditions that allows us to reconstruct f β from f - β . The following proposition can be found in [START_REF] Rubin | Inversion and Characterization of the Hemispherical Transform[END_REF]. When d ≥ 4 we have

h(n, d) = 2 (d -2)! (n + (d -2)/2)[(n + 1)(n + 2) • • • (n + d -3)],
and the results follow.

Next we turn to (A.27). When d is even and p ≥ d/2

|λ(2p + 1, d)| = κ d (2p + 1)(2p + 3) • • • (2p + d -1)
where

κ d = |S d-2 |1 • 3 • • • (d -1)
d -1 and (A.27) follows. Sterling's double inequality (see [START_REF] Feller | An Introduction to Probability Theory and Its Applications -Third edition[END_REF] p.50-53), that is,

√ 2πn n+1/2 exp -n + 1 12n + 1 < n! < √ 2πn n+1/2 exp -n + 1 12n , implies that (2 p p!) 2 (2p)! ≍ √ p and therefore 1 • 3 • • • (2p -1) ≍ √ p2 • 4 • • • (2p).
Thus for p ≥ d/2 and d odd we have

|λ(2p + 1, d)| ≍ √ p (2p + 2)(2p + 4) • • • (2p + d -1)
and (A.27) holds for both even and odd d.

We can now easily check that Proposition A.6. For all s > 0, there exists positive constants C l and C u such that for all f in

H s (S d-1 ) C l f - 2,s ≤ H(f -) 2,s+d/2 ≤ C u f - 2,s .
Proof of Proposition A. [START_REF] Billingsley | Probability and Measure -Third Edition[END_REF]. By definition we have

H f -2 2,s+d/2 = ∞ p=0 (1 + ζ 2p+1,d ) s+d/2 Q 2p+1,d H(f -) 2 2
where according to the Funk-Hecke Theorem

Q 2p+1,d H(f -) = Q 2p+1,d H   ∞ q=0 Q 2q+1,d f   = Q 2p+1,d   ∞ q=0 λ(2q + 1, d)Q 2q+1,d f   = λ(2p + 1, d)Q 2p+1,d f.
The result follows since Lemma A.2 implies that (1

+ ζ 2p+1,d ) s+d/2 λ 2 (2p + 1, d) ≍ (1 + ζ 2p+1,d ) s .
The factor d/2 in Proposition A.6 corresponds to the degree of "regularization" due to smoothing by H. Now the inverse of an odd function f -is given by (A.28)

H -1 (f -)(y) = ∞ p=0 1 λ(2p + 1, d) S d-1 q 2p+1,d (x, y)f -(x)dσ(x).
This is straightforward given our results at hand: for example, operate H on the RHS to see:

H   ∞ p=0 1 λ(2p + 1, d) S d-1 q 2p+1,d (x, y)f -(x)dσ(x)   = ∞ p=0 1 λ(2p + 1, d) HQ 2p+1,d f - = ∞ p=0 λ(2p + 1, d) λ(2p + 1, d) Q 2p+1,d f -(by the Funk-Hecke Theorem) = f -. If f -belongs to H d/2 (S d-1 ), then H -1 (f -)(b) is a well-defined L 2 (S d-1
) function. Otherwise it should be understood as a distribution and is only defined in a Sobolev space with negative exponent.

Moreover, if d is a multiple of 4, it is possible to relate the inverse of the operator H with differentiation as in the case of d = 2:

Proposition A.7. If d is a multiple of 4,

H -1 = |S d-2 | d/4 k=1 [-∆ S + 2(k -1)(d -2k)].
Proof of Proposition A.7. If we consider the case where d is even, we know from Proposition A.4,

that 1 λ(2p + 1, d) = (-1) p |S d-2 |(2p + 1)(2p + 3) . . . (d + 2p -1).
Thus if d is a multiple of 4,

1 λ(2p + 1, d) = |S d-2 | d/4 k=1 [-ζ 2p+1,d + 2(k -1)(d -2k)].
Using this and (4.1),

H -1 = ∞ p=0 1 λ(2p + 1, d) Q 2p+1,d = ∞ p=0 |S d-2 |  d/4 k=1 [-ζ 2p+1,d + 2(k -1)(d -2k)]   Q 2p+1,d .
Recall (A.18) and the proposition is proved.

This connection between the inverse of H and differentiation suggests that a Bernstein-type inequality might hold for H -1 . Indeed, even though the above inversion formula is concerned with d's that are multiples of 4, the following Bernstein inequality holds for every dimension.

Theorem A.3 (Bernstein inequality). For every d ≥ 2 and every q ∈ [1, ∞], there exists a positive constant B(d, q) such that for all P in T p=0 H 2p+1,d , (A.29)

H -1 P q ≤ B(d, q)T d/2 P q .
Proof of Theorem A.3. We can write

H -1 = P 1 (D) -P 2 (D)
where P 1 (D) and P 2 (D) are defined for all odd function f -by

P 1 (D)f -= ∞ p=0 1 λ(4p + 3) S d-1 q 4p+3 (x, y)f -(x)dσ(x) P 2 (D)f -= - ∞ p=0 1 λ(4p + 1) S d-1 q 4p+1 (x, y)f -(x)dσ(x).
P 1 (D) and P 2 (D) are two unbounded operators on B = L q odd (S d-1 ) with non-positive eigenvalues. We apply Theorem 3.2. of Ditzian (1998) to -P 1 (D) and -P 2 (D) choosing α = 1. Condition (1.6) of Ditzian (1998) can be verified using Proposition 2.2 with r = 1 and p = q and the fact that for the Cesaro kernels C l h are uniformly bounded in L 1 (S d-1 ) for l > d-2 2 (see, e.g. [START_REF] Bonami | Sommes de Cesàro et Multiplicateurs des Développements en Harmoniques Sphériques[END_REF]. We see, using the triangle inequality, that for all P in T p=0 H 2p+1,d ,

H -1 P q ≤ C 1 λ 2 (2T + 1, d) P q ≤ CT d P q .
The last inequality follows from (A.27).

Rubin (1999) gives other inversion formulas for the Hemispherical transform in terms of differential operators. The fact that the inversion roughly corresponds to differentiation is another manifestation of the ill-posedness of our problem at hand. The inverse operator H -1 is indeed unbounded. We call the factor d/2 in (A.29) the degree of ill-posedness of the inverse problem. For the case q = 2, there exists a lower bound for H -1 P q in (A.29) of order T d/2 as well, implying that the upper bound T d/2 in the order of T obtained in Theorem A.3 is tight.

A.1.5. Estimators for the Choice Probability Function. This section considers estimation of the choice probability function r and its extension R. We propose an estimator for r, which, in turn, yields a computationally simple estimator for f β . Also the asymptotic results presented here are useful for the next section where we study the limiting properties of our estimator for the random coefficients density f β .

Since R is square integrable on S d-1 , it has a condensed harmonic expansion which enables us to obtain the expressions in the next theorem.

Theorem A.4. For x in S d-1 , we have

(A.30) R(x) = 1 2 + ∞ p=0 E (2Y -1) f X (X) q 2p+1,d (X, x) .
This suggests an estimator of the form R1

(x) = 1 2 + R- 1 with R- 1 (x) = 1 N N i=1 (2y i -1) fX (x i ) T N p=0 q 2p+1,d (x i , x)
where fX is an estimator of f X and T N is a suitably chosen sequence diverging to infinity with N .

Note that the second summation corresponds to the Dirichlet kernel. We can generalize this, by introducing a class of estimators of the form

(A.31) R- 2 (x) = 1 N N i=1 (2y i -1) fX (x i ) K - 2T N (x i , x)
where K - 2T N is the odd part of a kernel of the form (A.23) satisfying Assumption A.1, such as the two kernels in Proposition A.3.

The estimator (A.31) is convenient, though the plug-in term fX has to be treated with care.

We avoid restrictive assumptions on the distributions of covariates and allow f X (x) to decay to zero as x approaches the boundary of its support H + . To deal with the latter problem, we modify (A.31) by

(A.32) R-(x) = 1 N N i=1 (2y i -1)K - 2T N (x i , x) max fX (x i ), m N
where m N is a trimming factor going to 0 with the sample size. Our estimator for R is then

(A.33) R = 1 2 + R-.
Remark A.4. Alternative estimators of R -are available. For example, one may use kernel regression on the sphere to estimate r in order to obtain an estimator for R -. As noted before, however, we then need to use numerical integration to evaluate (4.5) to calculate fβ .

Proof of Theorem A.4. R has the following condensed harmonic expansion

R(x) = 1 2 + ∞ p=1 (Q 2p+1,d R)(x).
We then write using (3.2), changing variables and using (A.11),

(Q 2p+1,d R)(x) = S d-1 q 2p+1,d (x, z)R(z)dσ(z) = H + q 2p+1,d (x, z)r(z)dσ(z) + H - q 2p+1,d (x, z)(1 -r(-z))dσ(z) = H + q 2p+1,d (x, z)r(z)dσ(z) - H + q 2p+1,d (x, z)(1 -r(z))dσ(z) = H + q 2p+1,d (x, z)(2r(z) -1)dσ(z) = H + q 2p+1,d (x, z)E 2Y -1 f X (z) X = z f X (z)dσ(z) = E (2Y -1)q 2p+1,d (x, X) f X (X) .
A.1.6. Proofs of Main Results.

Proof of Proposition 3.1. It is straightforward that the model (1.1) and Assumption 1.1 imply that the choice probability function r given by (1.2) is homogeneous of degree 0. Proposition A.5

along with the fact that R =

1 2 + H f - β with f - β ∈ L 2 odd (S d-1
) implies that R belongs to H d/2 (S d-1 ). We now turn to the proof of sufficiency. If the extension R given by (3.2) belongs to H d/2 (S d-1 ) then so does R -and Proposition A.5 shows that there exists a unique odd function

f -in L 2 (S d-1 ) such that R = 1 2 + H f -= H 1 |S d-1 | + f -.
Moreover, since 0 ≤ R(x) ≤ 1 holds for every x ∈ S d-1 , the above relationship implies that

1 2 ≥ Hf -(x), ∀x ∈ S d-1 . But Hf -(x) ≥ {f -(b)≥0} f -(b)dσ(b) holds for some x. Therefore we conclude that 1 2 ≥ {f -(b)≥0} f -(b)dσ(x) = -{f -(b)≤0} f -(b)dσ(b), thus S d-1 |f -(b)|dσ(b) ≤ 1. Also, following the discussion in Section A.1.3, 1 |S d-1 | + f -integrates to 1.
We have seen in Corollary A.1 that for even function g that has 0 as the coefficient of degree 0 in its expansion on the surface harmonics (i.e.

an even function that integrates to zero over the sphere),

R = H g + 1 |S d-1 | + f - holds. Now consider g = |f -| - 1 |S d-1 | S d-1 |f -(b)|dσ(b),
then this certainly is even and integrates to zero. Using this, define

f * β := g + 1 |S d-1 | + f -= 2f -I{f -> 0} + 1 |S d-1 | 1 - S d-1 |f -(b)|dσ(b) ≥ 0.
Obviously f * β -= f -. This function f * β is non-negative and integrates to one, and thus it is a proper probability density function (pdf). It is indeed bounded from below by

1 |S d-1 | 1 -S d-1 |f -(b)|dσ(b) . As a consequence, there exists a pdf f * β such that R = H f * β = 1 2 + H f * β - and for all x in H + , r(x) = H f * β (x).
Proof of Theorem 4.1. We use the shorthand notation

I(b) := I{f - β (b) > 0} and Î(b) := I{ f - β (b) > 0}. Then f β = 2f - β I and fβ = 2 f - β Î. We write f - β,T (b) = 1 N N i=1 (2y i -1)H -1 K - 2T N (x i , •) (b) max (f X (x i ), m N )
A similar proof can be carried out replacing L q (S d-1 ) by L ∞ (S d-1 ). Thus it is enough to consider the behavior of fβ -f - β instead of fβ -f β . As noted above, the former can be decomposed into four terms, S p , S e , B t and B a .

We start with the analysis of S p . Note that for q ∈ [1, ∞]

S p q = H -1   1 N N i=1 (2y i -1)K - 2T N (x i , •) max(f X (x i ), m N )   max (f X (x i ), m N ) max fX (x i ), m N -1     q ≤ B(d, q)T d/2 N 1 N N i=1 (2y i -1)K - 2T N (x i , •) max(f X (x i ), m N )   max (f X (x i ), m N ) max fX (x i ), m N -1   q
(by Theorem A.3)

≤ B(d, q)T d/2 N m -1 N 1 N N i=1 |K 2T N (x i , •)| q max i=1,...,N max (f X (x i ), m N ) max fX (x i ), m N -1 ≤ B(d, q)T d/2 N m -2 N 1 N N i=1 |K 2T N (x i , •)| q max i=1,...,N f X (x i ) -fX (x i )
holds, where we have used the triangle inequality. The L q -norm on the right hand side is bounded from above by q + E |K 2T N (X, •)| q := T 1 q + T 2 q .

First consider the term T 1 q . We begin with the case of q ∈ [START_REF] Athey | Discrete Choice Models with Multiple Unobserved Choice Characteristics[END_REF][START_REF] Bajari | Linear Regression Estimation of Discrete choice Models with Nonparametric Distribution of Random Coefficients[END_REF]. By the Hölder inequality,

E T 1 q q = S d-1 E [T 1 (x) q ] dσ(x) ≤ S d-1
E T 1 (x) 2 q/2 dσ(x)

where N m -2 N T 1 q = O p (log N ) 1/2-1/q m -2 N N -1/2 T (2d-1)/2 N .

E T 1 (x) 2 ≤ 1 N E (K 2T N (X, x)) 2 (A.36) ≤ C N K 2T N (⋆ 2 , x)
Concerning T 2 q , q ∈ [1, ∞], since f X is bounded by assumption, there exists a positive C such that

T 2 q ≤ C K 2T N (⋆ 1 , ⋆ q ) 1 q
where integration in • 1 is with respect to argument ⋆ 1 and integration in • q is with respect to ⋆ q . But K 2T N (⋆ 1 , ⋆ q ) 1 is a constant and does not depend on ⋆ q , as previously noted. Thus

K 2T N (⋆ 1 , ⋆ q ) 1 q = |S d-1 | 1/q K 2T N (⋆ 1 , ⋆ q ) 1
and we conclude that this term is O(1) using Assumption A.1 (i) on the kernel, thus

T d/2 N m -2 N T 2 q = O m -2 N T d/2 N .
Analogously to our treatment of T 1 q , we can prove that when q ∈ [1, 2],

S e q = O p m -1 N N -1/2 T (2d-1)/2 N , while for q ∈ (2, ∞] S e q = O p m -1 N (log N ) 1/2-1/q N -1/2 T (2d-1)/2 N .

Let us now turn to the bias term induced by trimming

B t (b) = E (2Y -1)H -1 K - 2TN (X, •) (b) f X (X) f X (X) max(f X (X), m N ) -1 = {z∈S d-1 : fX (z)<mN } E[2Y -1|X = z]H -1 K - 2TN (z, •) (b) f X (z)m -1 N -1 dσ(z).
2 1 f X (z) I{f X ≥ m N } + f X (z)m -2 N I{f X < m N } dσ(z) ≥ 4 1 f X ∞ H + TN -1 p=0 χ(2p + 1, 2T N ) q 2p+1,d (z, b) λ(2p + 1, d) 2 dσ(z) -4 1 f X ∞ {fX <mN } TN -1 p=0 χ(2p + 1, 2T N ) q 2p+1,d (z, b) λ(2p + 1, d) 2 dσ(z)
With similar computations as (A.36), using as well (A.27), we know that there exists a constant C such that ) and L 1 (S d-1 ) norms of the kernels using the Hölder inequality: We now need to prove that the remaining terms S p , B t and B a , multiplied by N 1/2 s -1 N (b), are o p (1). The term S p is treated in a similar manner as in the proof of Theorem 4.1.

E |Z N | 2+δ ≤ f X ∞ m -(2+δ) N H -1 K - 2T N (z, •) 2+δ 2+δ ≤ f X ∞ m - (
|S p (b)| ≤ 2   1 N N i=1 H -1 K - 2T N (x i , •) (b) max(f X (x i ), m N )   max i=1,...,N max (f X (x i ), m N ) max f N X (x i ), m N -1 .
Using the Markov inequality, the empirical average in the parenthesis is of the stochastic order of

m -1 N H -1 K - 2T N (⋆, •) 1 . But m -1 N H -1 K - 2T N (⋆, •) 1 ≤ B(d, 1)T d/2 N m -1 N K - 2T N (⋆, •) 1 ≤ B(d, 1)T d/2 N m -1 N K 2T N (⋆, •) 1
where the first inequality follows from Theorem A.3 and the second is obtained using the definition of the odd part and the triangle inequality. Note that the term K 2T N (⋆, •) 1 in the last line

  /2 f p where ∆ S denotes the Laplacian on the Sphere S d-1 : See Section A.1.3 for further discussions. Proposition 3.1. A [0, 1]-valued function r is compatible with the random coefficients model (1.1) with f β in L 2 (S d-1 ) and Assumption 1.1 if and only if r is homogeneous of degree 0 and its extension R according to (3.2) belongs to H d/2 (S d-1 ).

  Nonparametric estimation of the regression function with random degenerate design, in the sense that the density of regressors can be low on its support, is a difficult issue. It has been studied for the pointwise risk in Hall et al. 1997, Gaïffas 2005, Gaïffas 2009 and Guerre 1999. Extension to inverse problems setting is a widely open problem. We tackle this problem for our specific inverse problem. Future research includes the study of lower bounds from the minimax point of view that account for the degeneracy of the design.

1 2s+2d- 1

 11 in Corollary 4.3, and we recover the L 2 convergence rate of N s 2s+2d-1 , the rate mentioned in Section 4.1. It is in accordance with the L 2 rate in Healy and

Theorem 4 . 2 (

 42 Asymptotic normality). Suppose f - β belongs to W s ∞ (S d-1) with s > 0, and Assumptions A.1, 3.1 and 4.1 hold. If fX , f X , m N and T N satisfy

5. 1 . 2 .

 12 Estimation of Marginals. In Section 3 we have provided an expression for the estimator of the full joint density of β, from which an estimator for a marginal density can be obtained. Let σ k denote the surface measure and σ k = σ k /|S k | the uniform probability measure on S k . We writeβ = β ′ , β ′ ′and wish to obtain the density of the marginal of β which is a vector of dimension k. Also define P and P the projectors such that β = P β and β = P β and denote by P * σ d-1 and P * σ d-1 the direct image probability measures. One possibility is to define the marginal law of β as the measure P * P β , where dP β = f β dσ. This may not be convenient, however, since the uniform distribution over S d-1 would have U-shaped marginals. The U-shape becomes more pronounced as the dimension of β increases.In order to obtain a flat density for the marginals of the uniform joint distribution on the sphere it is enough to consider densities with respect to the dominating measure P * σ d-1 . Notice that sampling U uniformly on S d-1 is equivalent to sampling U according to P * σ d-1 and then given U forming ρ U V where V is a draw from the uniform distribution σ d-1-k on S d-1-k and ρ U = 1 -U Indeed given U , U /ρ U is uniformly distributed on S d-1-k . Thus, when g is an element of L 1 (S d-1 ) we can write for k in {1, . . . , d -1},(5.1)

ρσ 2 σ 2 

 2 where µ = 0.7, σ 2 = 0.3, ρ = 0.5 and λ = 0.5. Random samples of size 500 from each of the two specifications are generated, then the new estimator (4.8) is computed. It is implemented using the Riesz kernel with s = 2 and l = 3 (see Proposition A.3). The truncation parameter T N is set at 3, and the trimming parameter r is 2. It also requires a nonparametric estimator for f X , and we use the projection estimator (4.7) based on the same Riesz kernel (i.e. s = 2, l = 3) and T N = 10.

Figure 1 .

 1 Figure 1. Nonparametric estimator of f β for Model 1 (left) and Model 2 (right)

A- 1 SUPPLEMENTAL

 1 APPENDIX FOR "NONPARAMETRIC ESTIMATION IN RANDOM COEFFICIENTS BINARY CHOICE MODELS" ERIC GAUTIER AND YUICHI KITAMURA A.1.1. A Toy Model. As noted in the main text, the key insight for our estimation procedure lies in the fact the estimation of f β in (1.3) is mathematically equivalent to a statistical deconvolution problem. To see this, it is useful to first consider the case with d = 2. We parameterize the vectors b = (b 1 , b 2 ) ′ and x = (x 1 , x 2 ) ′ on S 1 by their angles φ = arccos (b 1 ) and θ = arccos (x 1 ) in [0, 2π).

A. 1 . 4 .

 14 The Hemispherical Transform. The hemispherical transform H, defined by Hf (x) = S d-1 I{x ′ y ≥ 0}f (y)dσ(y), plays a central role in our analysis. It is a special case of the operator considered in the Funk-Hecke theorem above, with F (t) = I{t ∈ [0, 1]}, therefore the next proposition follows. Notation. We define λ(n, d) = λ n (I {t ∈ [0, 1]}) for d ≥ 3 and λ(n, 2) = 2 sin(nπ/2) n . Proposition A.4. When d ≥ 2, the coefficients λ(n, d) have the following expressions

4 .

 4 Define α(n, d) := C ν(d) n (1)|S d-2 | -1 λ n (I {t ∈ [0, 1]}). By the Funk-Hecke theorem α(n, d) = 1 0

  4, and corresponds to an observation made in Remark A.1 for the d = 2 case.

Proposition A. 5 .

 5 H is a bijection from L 2 odd (S d-1 ) to H d/2 odd (S d-1 ). Lemma A.2. h(n, d) ≍ n d-2 , (A.26) |λ(2p + 1, d)| ≍ p -d/2 . (A.27)Proof. Estimate (A.26) is clearly satisfied when d = 2 and 3 since h(n, 2) = 2 and h(n, 3) = 2n + 1.

|K

  2T N (x i , •)| -E |K 2T N (X, •)|

q 2 ( 2 N m - 2 N sup x∈S d- 1 |T 1 / 2 N m - 2 N T 1 ∞

 22211221 2T N )q n,d (⋆ 2 , x) n,d (⋆ 2 , x) 2 by Assumption A.1(iv)) positive constants C and C 2 , (A.44) P v N T d/(x)| ≥ M ≤ C exp (log N )(C 1 -C 2 M 2 )holds. For a large enough M , C 1 -C 2 M 2 < 0 and the right hand side of (A.44) converges to zero, so (A.39) follows. In summary, we have just shown thatT d= O p (log N ) 1/2 m -2 N N -1/2 T (2d-1)/2 Nand with (A.37) and (A.38) we also conclude thatT d/2

1 HS d- 1 H 2 N m - 2 N sup x∈S d- 1 |T 1

 112211 This yields|B t (b)| ≤ S d--1 K - 2TN (z, •) (b) I z ∈ S d-1 : f X (z) < m N dσ(z) = -1 K - 2TN (b, •) (z) I z ∈ S d-1 : f X (z) < m N dσ(z) (using the condensed Harmonic expansion), in Corollary 4.1, noticing that the bias terms B t and B a are not stochastic and bounded after proper scaling, we just have to focus on S p and S e appearing in the proof of Theorem 4.1. Concerning S p ,proceed as before and note that taking M large enough so that C 1 -C 2 M 2 < -1 implies summability of the left hand side in (A.44). We conclude from the first Borel-Cantelli lemma that the probability that the events occur infinitely often is zero thus with probability onelim N →∞ v -1 N B(d, ∞)T d/(x)| < M.The term T 2 is non-stochastic and its treatment in our previous analysis remains valid, therefore we can use the same non-stochastic upper bound. We then use Assumption 4.2 (iii) instead of Assumption 4.2 (ii) to show almost sure uniform boundedness of S p after proper rescaling. The treatment of S e is analogous to that of T 1 . The proof is the same in Corollaries 4.2 and 4.3.Proof of Theorem 4.2. We first prove that the Lyapounov condition holds: there exists δ > 0 such that for N going to infinity,(A.47) E |Z N (b) -E [Z N (b)]| 2+δN δ/2 (var (Z N (b))) 1+δ/2 → 0 (see, e.g. Billingsley, 1995). We start from deriving a lower bound on var (Z N (b)). Since E[Z N (b)] converges to f - β (b), it is enough to obtain a lower bound on + 1, 2T N ) q 2p+1,d (z, b) max (f X (z), m N ) λ(2p + 1, d) + 1, 2T N ) q 2p+1,d (z, b) λ(2p + 1, d)

T N - 1 p=0χ(2p + 1 , 2 ≤ CT 2d- 1 N,

 1121 2T N ) q 2p+1,d (z, ⋆) λ(2p + 1, d) therefore using Proposition A.2 with p = q = r = 1 we obtain E[Z 2 N ](b) ≥ 4 f X ∞ TN -1 p=0 χ(2p + 1, 2T N ) 2 H + q 2p+1,d (z, b) 2 λ(2p + 1, d) 2 dσ(z) -CT 2d-1 N σ (f X < m N ) .

A- 25 Using 2 i 1 N. 1 N.

 25211 Assumption A.1 (iv), the first term on the right hand side can be bounded from below byC ⌊(T N -1)/2⌋ p=0 q 2p+1,d (z, b) λ(2p + 1, d)2 .e. by CT 2d-Thus as m N decays to zero, σ (f X < m N ) decays to zero and (A.48) E[Z 2 N ](b) ≥ CT 2d-We now derive an upper bound of E |Z N (b)| 2+δ using Theorem A.3 and interpolation between L ∞ (S d-1

.- 2 . 1 N

 21 2+δ) N B(d, 2 + δ) 2+δ T d(2+δ)/2 N K - 2T N (z, •)By this and (A.48) an upper bound for the ratio appearing in (A.47) is given by m Therefore the Lyapounov condition is satisfied if (4.20) holds, and it follows that N 1/2 s -

Note that the dimension of their estimator is d, whereas that of ours is d -1. On the other hand, in their problem W is observable, and it is obviously more informative than our binary outcome Y , which causes difficulties both in identification and estimation.

for helpful comments. Yuhan Fang and Xiaoxia Xi provided excellent research assistance. Kitamura acknowledges financial support from the National Science Foundation via grants SES-0241770, SES-0551271 and SES-0851759. Gautier is grateful for support from the Cowles Foundation as this research was initiated during his visit as a postdoctoral associate.

The Gegenbauer polynomials are orthogonal with respect to the weight function (1 -t 2 ) ν-1/2 dt on [-1, 1]. Note that C ν 0 (t) = 1 and C ν 1 (t) = 2νt for ν = 0 while C 0 1 (t) = 2t. Moreover, the following recursion relation holds (A.4) (n + 2)C ν n+2 (t) = 2(ν + n + 1)tC ν n+1 (t) -(2ν + n)C ν n (t).

Implementation of our estimator requires evaluation of the Gegenbauer polynomials for a series of successive values of n. The recursion relation (A. [START_REF] Beran | Estimating Coefficient Distributions in Random Coefficient Regression[END_REF]) is therefore a powerful tool. The Gegenbauer polynomials are related to each other through differentiation, that is, they satisfy

for ν > 0 and

For ν = 0 the Rodrigues formula states that (A.7)

The following results are also used in the paper:

These orthogonal polynomials are normalized such that (A.12)

A.1.3. Tools for Higher Dimensional Spheres. Let us introduce some concepts used for the treatment of the general case d ≥ 2. We consider functions defined on the sphere S d-1 , which is

signifies the surface area of the unit sphere.

Recall that the basis functions exp(±int)/ √ 2π are eigenfunctions of -d dt 2 associated with eigenvalue n 2 . In a similar way, the Laplacian on the sphere S d-1 , d ≥ 2, denoted by ∆ S , can be used to obtain an orthonormal basis for higher dimensional spheres. It can be defined by the formula

here ∆ is the Laplacian in R d , fˇthe radial extension of f , that is fˇ(x) = f (x/ x ), and fˆthe restriction of f to S d-1 . Likewise the gradient on the sphere is given by:

Definition A.1. A surface harmonic of degree n is the restriction of a homogeneous harmonic polynomial (a homogeneous polynomial p whose Laplacian ∆p is zero) of degree n in R d to S d-1 .

The reader is referred to [START_REF] Müller | Spherical Harmonics[END_REF] and [START_REF] Groemer | Geometric Applications of Fourier Series and Spherical Harmonics[END_REF] for clear and detailed expositions on these concepts and important results concerning spherical harmonics used in this paper. Erdélyi (iii ) The dimension of the vector space H n,d of surface harmonics of degree n is

(iv ) A system formed of orthonormal bases (Y n,l )

of H n,d for each degree n = 0, . . . , ∞ is complete in L 1 (S d-1 ), that is, for every f ∈ L 1 (S d-1 ) the following equality holds in the L 1 (S d-1 ) where

q n,d (x, y).

The kernel D T extends the classical Dirichlet kernel on the circle to the sphere S d-1 . The sum over T in the definition of D T also has the simple closed form in terms of derivatives of Gegenbauer polynomials; see Equation (52) in [START_REF] Müller | Spherical Harmonics[END_REF]. The linear form f → S d-1 D T (x, y)f (y)dσ(y) converges to S d-1 f (y)dδ x (y) = f (x) as T goes to infinity, where δ x denotes the Dirac measure. The Dirichlet kernel yields the best approximation

) by polynomials that belong to T n=0 H n,d , but is known to have flaws. For example, D T does not satisfy

that is, the sequence D T , T = 0, 1, ... is not an approximate identity (see, e.g., [START_REF] Devroye | Nonparametric Density Estimation. The L1-View[END_REF] in L 1 (S d-1 ). Indeed, the L 1 (S d-1 ) norm of the kernel is not uniformly bounded; more precisely, we have

when d ≥ 3 and

when d = 2 (as noted above, these norms do not depend on the value of x ∈ S d-1 ). These bounds can be found in [START_REF] Gronwall | On the Degree of Convergence of Laplace's Series[END_REF] for d = 3 and Ragozin (1972) and Colzani and Traveglini (1991) for higher dimensions. Also, D T does not have good approximation properties in L ∞ (S d-1 ); in particular, we do not have

Near the points of discontinuity of f , D T * f has oscillations which do not decay to zero as T grows to infinity, known as the Gibbs oscillations. This phenomenon deteriorates as the dimension increases.

These problems can be addressed by using kernels that involves extra smoothing instead of the Dirichlet kernel D T . To this end, define a general class of kernel

for some sequence χ(n, T ). These are called smoothed projection kernels. Typically the function χ is chosen so that it puts more weight on lower frequencies. In particular we impose the following conditions:

Assumption A.1. Like the delayed means kernel, it achieves the optimal rate of approximation without the prior knowledge of s.

Proof of Proposition A.

Ditzian's (1998) Theorem 6.1, which gives an upper bound on the best spherical harmonic approximation in L p (S d-1 ) to functions in W s p (S d-1 ) (see also [START_REF] Kamzolov | The Best Approximation of the Class of Functions W α p (S n ) by Polynomials in Spherical Harmonics[END_REF]), then apply equation (4.1) in Ditzian (1998) again to obtain the desired result.

If the function f is in L 2 (S d-1 ) then Equations (A.17) and (A.11)

Consequently, the odd order terms in the condensed harmonic expansions of f , f + and f -satisfy

Likewise, for the even order terms in the condensed harmonic expansions of these functions Q 2p f + = Q 2p f and Q 2p f -= 0 hold. We conclude that the sum of the odd order terms in the condensed harmonic expansion corresponds to f -and that of the even order terms to f + . As anticipated from the analysis of the d = 2 case, the operator H reduces the even part of f β to a constant 1 2 , therefore Fourier-Laplace series expansions for f β derived later involve only odd order terms.

We now provide a formula that is used to obtain our estimator for f β . If a non-negative function f has its support included in some hemisphere of S d-1 then (A. [START_REF] Guerre | Efficient Random Rates for Nonparametric Regression Under Arbitrary Designs[END_REF])

Denote the support of f by suppf and let -suppf = {x| -x ∈ suppf }, then this formula follows from the fact that f -(x) = f + (x) ≥ 0 on suppf while f -(x) = -f + (x) ≤ 0 on -suppf and both f -and

and use the decomposition

and denote the terms on the right hand side by S p (stochastic component due to plug-in), S e (stochastic component of the infeasible estimator f β,T ), B t (trimming bias) and B a (approximation bias). Take q ∈ [1, ∞),

Obviously

and A 4 = 0. Also,

Similarly,

and given

By the Markov inequality, (A.37)

, providing a convergence rate for T 1 q , q ∈ [START_REF] Athey | Discrete Choice Models with Multiple Unobserved Choice Characteristics[END_REF][START_REF] Bajari | Linear Regression Estimation of Discrete choice Models with Nonparametric Distribution of Random Coefficients[END_REF]. So if we can establish a similar rate for T 1 ∞ , all L q (S d-1 ) convergence rates of T 1 for q ∈ (2, ∞] can be interpolated between the L 2 (S d-1 ) and L ∞ (S d-1 ) convergence rates using the following inequality:

To see this, note

We can thus focus on T 1 ∞ . We cover the sphere S d-1 by N(N, r, d) geodesic balls (caps) (B i )

and radius R(N, r, d), that is, We now prove that for every ǫ > 0 positive, there exists a positive M such that (A.39)

holds for an appropriately chosen sequence v N ↑ ∞. Write

where the last inequality is obtained using Assumption A.1 (ii) on the kernel and letting R(N, r, d) ≍

M (where α is given in Assumption A.1 (ii)). Notice

where

N ≤ L (using (A.17) and (A.8)).

The bound L in the last line is obtained by noting that

N , which follows from (A.17 . ω is the leading term in the denominator of the exponent in the last inequality.

If we take v

Also, use this v N in our choice of R(N, r, d) made above to get:

for some constant C 1 that might be greater than 1 2 (d -1), depending on the value of α. Indeed, T N does not grow more than polynomially fast in N . (A.40), (A.41), (A.42) and (A.43) imply that, for a thus, for every 1 ≤ r ≤ q,

(from Proposition A.2)

where in the last inequality we use Theorem A.3 and calculate an upper bound on the L r -norm of the kernel by interpolation, using Hölder's inequality, between the uniformly bounded L 1 -norm and the upper bound on the sup norm of the order of T d-1 N seen previously, C is a constant. We finally treat B a using Assumption A.1 (iii) with the condition that f - β ∈ W s q (S d-1 ):

In the case where f X ≥ m σ a.e., we use the decomposition

and we can treat the terms Sp and Se on this event.

Proof of the corollaries 4.

which is what we get from (4.9) and (4.10). Indeed, it is enough to find γ(ρ) as the minimum of

The first is an increasing function of γ while the second and fourth are decreasing Let us now consider the bias term induced by the trimming procedure. In the proof of Theorem