Robust functional estimation in the multivariate partial linear model
Michael Levine
math.ST stat.TH






   






Robust functional estimation in the multivariate partial linear model
    Michael Levine
    
=====================================================================




We consider the problem of adaptive estimation of the functional component in a multivariate partial linear model where the argument of the function is defined on a q-dimensional grid. Obtaining an adaptive estimator of this functional component is an important practical problem in econometrics where exact distributions of random errors and the parametric component are mostly unknown and cannot safely assumed to be normal. An estimator of the functional component that is adaptive in the mean squared sense over the wide range of multivariate Besov classes and robust to a wide choice of distributions of the linear component and random errors is constructed. It is also shown that the same estimator is locally adaptive over the same range of Besov classes and robust over large collections of distributions of the linear component and random errors as well. At any fixed point, this estimator also attains a local adaptive minimax rate. The procedure needed to obtain such an estimator turns out to depend on the choice of the right shrinkage approach in the wavelet domain. We show that one possible approach is to use the multivariate version of the classical BlockJS method. The multivariate version of BlockJS is developed in the manuscript and is shown to represent an independent interest. Finally, the Besov space scale over which the proposed estimator is locally adaptive is shown to depend on the dimensionality of the domain of the functional component; the higher the dimension, the larger the smoothness indicator of Besov spaces must be.  





§ INTRODUCTION


In this manuscript, we consider a partial linear multivariate model defined as  

    Y_i⃗=a+X_i⃗^'β⃗+f(U_i⃗)+ξ_i⃗

where X_i⃗∈^p and U_i⃗∈^q, β is an unknown p× 1 vector of parameters, a an unknown intercept term, f(·) is an unknown function, and ξ_i⃗ are independent and identically distributed random variables. No moment conditions are imposed on ξ_i⃗ but we do assume for convenience that the median of ξ_i⃗ is equal to zero. X_i⃗ is defined as a p-dimensional continuous random variable. In this manuscript, we consider a special case where each U_i⃗ is viewed as a q-dimensional vector with each coordinate defined on an equispaced grid on [0,1].  The sequence {X_i⃗} is assumed to be independent of {ξ_i⃗}. We use bold font for indices since  the most convenient notation for this model involves multivariate indices; the detailed description of these indices is postponed until Section (<ref>). In our manuscript, we only consider the case where q>1 that has been relatively little explored in statistical literature. More specifically, the only papers in the statistical literature that we are aware of discussing the multivariate case are <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. Econometric literature discusses the multivariate case to some extent. In particular, <cit.> contains a review of some possible applications, clearly showing practical utility of considering the case of q>1. 

Partial linear models are, in many cases, preferable to purely nonparametric regression model because of the well-known “curse of dimensionality". For the most part, the parametric part can be estimated at the √(n) rate where n is the sample size. At the same time, the estimation precision of the nonparametric component usually decays as the dimensionality of its argument grows. Partial linear models have a long history of application in both econometrics and statistics. From the practical viewpoint, they are often of considerable interest  because relationships between the response and predictors in the same model may be of a very different nature. Some of these relationships can often presumed to be linear while others are harder to parameterize. In the most common case, a small subset of the variables are presumed to have an unknown nonlinear relationship with the response while the rest are assumed to have a linear relationship with it. A good example can be found in <cit.> that considered a partial linear model to analyze gasoline household consumption in the US. In this model the demand for gasoline, measured in log number of gallons, is assumed to depend linearly on the number of drivers in the family, the household size, residence, region, and lifecycle. At the same time, the response is assumed to depend nonlinearly on the two remaining covariates, log of the household income and log of the age of the head of household.  

Most often, the estimation of the nonparametric component is conducted in order to suggest a possible parametric form for this component where no prior rationale for choosing such a form is available. This, in turn, allows a researcher to describe the data more parsimoniously. For example, the above mentioned <cit.> fits the function g that describes the dependence of the log demand for gasoline Y on the log of the household income Z^1 and the log of the age of the head of household Z^2. When the function is fit, <cit.> plots it first as a function of Z^1 for several fixed values of Z^2 and as a function of Z^2 for several fixed values of Z^1. The plots suggest that a possible parsimonious representation of the nonparametric component can be a piecewise linear function in both Z^1 and Z^2; later diagnostic testing confirms this conclusion.  An in-depth discussion of this issue can be found in <cit.>. 

Our main goal in this manuscript is to construct an estimator of the nonparametric component f that is adaptive over a range of functional classes for f and robust with respect to the wide choice of distributions of X and ξ. In other words, we want to develop the estimator of the function f that achieves the optimal, or nearly optimal, rate of convergence over a wide range of functional classes for f and that is reasonably robust to choices of a distribution of X and that of random errors ξ. To the best of our knowledge, this question has not been considered before in statistical literature. The farthest step in this direction seems to have been made in <cit.> who constructed an asymptotically efficient estimator of f in the model (<ref>). Also, a somewhat related result in the existing literature can be found in <cit.>. They obtained some optimal convergence rates for a link function in the partial linear single index problem.  

There are a number of reasons why the estimator robust to the wide choice of possible distributions of X and ξ is of interest. First, typical theory for models of the type (<ref>) assumes that errors ξ_i are independent, identically distributed (i.i.d.) normal random variables. From the practical viewpoint, normality may not always be satisfactory; see, for example, <cit.> and <cit.>. Moreover, the use of Gaussian formulation often implies that the specifically mean regression is of interest; in other words, the value of f(U_i) is viewed as a conditional expectation [Y_i⃗|X_i⃗]-a-X_i⃗^'β. In general, however, it is desirable to be able to handle more general formulations, such as median regression and other quantile regressions, as well. This has been noticed as early as <cit.> who suggested using an M-type objective function in order to treat mean regression, median regression, and other quantile regressions in one setting. Until the present time, if normality has not been required, it has been the usual approach in most statistical and econometric research to impose some moment assumptions that are necessary to obtain asymptotic results for estimators of the parametric component β⃗; see, e.g. <cit.> and <cit.>. In this manuscript, we argue that this is also unnecessary and that the estimator robust to the wide choice of possible distributions for ξ and achieving minimax or nearly minimax rate of convergence can be constructed without any moment assumptions being imposed on the distribution of ξ.  

Another issue lies in the fact that the distribution of X is often not known in practice; in particular, it need not be a multivariate Gaussian. In econometric practice, in particular, it is exceedingly common to have to deal with a vector X that includes at least some discrete components. An earlier mentioned model <cit.>  has the vector X that consists of discrete components only. In general, <cit.>, p. 53 notes when introducing partial linear models in econometric context that “...X may be discrete or continuous". Thus, from the practical viewpoint it seems inadvisable to limit the distribution of X to the multivariate normal only.  In our manuscript, we develop an estimation method that is adaptive over a wide range of functional classes for f and robust over a large collection of error distributions for ξ_i⃗ and design vector distributions for X_i⃗.  This means, in particular, that an exact optimal rate of convergence for both mean squared error risk and the squared error at a point risk is achieved over a wide range of multivariate Besov classes for f without  the prior knowledge of the distribution for the design vector X or the random error distribution. 

The method that we propose is based on the idea of treating the sum of the parametric part X_i⃗β⃗ and the random error ξ_i⃗ as the “new" random error ρ_i⃗ and viewing the model (<ref>) as a nonparametric regression with unknown random errors. The usefulness of this idea lies in the fact that the resulting nonparametric regression model is somewhat similar to the nonparametric regression with an unknown error distribution considered in <cit.>. The main difference between it and the model of <cit.> is that the argument of the function f now is multivariate while it was univariate in <cit.>. Thus, although the origin of our problem is rather different from that of <cit.>, solutions of the two problems turn out to be connected. To obtain an adaptive estimator of the function f, we divide the q-dimensional cube [0,1]^q into a number of equal volume bins, take the median of observations in each of these bins, and then apply a wavelet based procedure to these local medians together with a bias correction. Out of several possible procedures, we choose a multivariate generalization of the Block JS procedure developed in <cit.>. To the best of our knowledge, the multivariate BlockJS procedure has not been properly described before and therefore may be of independent interest. Note that, in general, the multivariate extension of the adaptive estimation procedure described in <cit.> is highly non-trivial. This is due to the necessity of selecting a correct blockwise shrinkage procedure that can be applied to the empirical wavelet coefficients. Our contribution consists of, first, designing a multivariate version of the by now classical BlockJS procedure of <cit.>, and, second, of showing that it is the right choice that produces an adaptive and robust estimator of the functional component f. 

Our manuscript is organized as follows. In Section 2, we define our proposed procedure exactly and establish several auxiliary results. The asymptotic properties of our procedure are established in Section 3. Section 4 contains some further discussion of our work, while the formal proofs are relegated to the Appendix. 





§ GENERAL APPROACH TO ESTIMATION OF THE FUNCTIONAL COMPONENT




 §.§ Methodology for adaptive estimation of the functional component


As mentioned in the introduction, we begin with defining a random variable ρ_i⃗=X_i⃗^'β⃗+ξ_i⃗ and rewriting the model (<ref>) as 

    Y_i⃗=a+f(U_i⃗)+ρ_i⃗.
 
In this form, (<ref>) is simply a multivariate nonparametric regression with an unknown error distribution. Note that in this model the intercept a cannot be absorbed in the design matrix X due to identifiability issues; in order to ensure that the model is identifiable, we have to require that an identifiability condition ∫_[0,1]^q f(u)du=0 is satisfied. Otherwise, one can add and subtract ∫_[0,1]^q f(u)du to the right hand side of the model with the new constant becoming a^'=a+∫_[0,1]^q f(u)du. Without loss of generality, we will adopt the following simpifying assumptions. First, we assume without loss of generality that a is known and equal to 0. Second, we assume that the vector median of X_i⃗ is also equal to zero. Indeed, if this is not the case, <cit.> has shown that there exists an estimator â such that â=O_p(n^-1/2) and that the parameter β can be estimated at the same rate as well. Therefore, in a general model, we will have an additional term A:=a+(med X_i⃗)^'β that can be estimated at the same fast parametric rate of convergence and its presence will not influence the optimal rate of convergence of the adaptive estimator of f that we construct. Appropriate remarks will be made in the proofs in mathematical Appendix as well.  Therefore, from now on we will work with the model 

    Y_i⃗=f(U_i⃗)+ρ_i⃗

where it is assumed that the median of ρ_i⃗ is equal to zero. Most of the classical nonparametric regression has been developed under the assumption of independent and identically distributed (i.i.d.) errors. In particular, a variety of smoothing techniques, such as wavelet thresholding techniques, were developed and shown to be highly adaptive in the Gaussian case. When errors are heavy-tailed, these techniques are typically not applicable. <cit.> worked with the model (<ref>) when q=1 and noticed that, for example, if ρ_i⃗ is Cauchy distributed, the maximum observation (out of n) will be of the order n, instead of log n, which is the case when the errors are normally distributed. This invalidates classical denoising approches, such as the wavelet thresholding, and suggests the need for a different take on this problem. 

Similarly to <cit.>,  to estimate the function f adaptively, we bin Y_i⃗ according to the values of coordinates of U_i⃗. The sample median is then computed within each bin. Now, bin centers can be treated as independent variables in a multivariate nonparametric regression, bin medians being dependent variables. The number of bins has to be chosen in a suitable range; in our case, it turns out that the number of bins V≍ n^3/4 where n is the original sample size, is a suitable choice. The resulting model can be viewed as a Gaussian multivariate nonparametric regression and a multivariate version of the BlockJS method of <cit.> can be used to obtain an adaptive estimator of the function f. To the best of our knowledge, such a generalization of BlockJS method has not been implemented before. The implementation of the proposed procedure is not difficult, since the number of bins can be chosen as a power of 2. We will show that the resulting estimator enjoys excellent adaptivity properties over a wide range of multivariate Besov balls and is robust over wide ranging sets of distributions of X and ξ.  

Before the detailed description of our procedure, it is necessary to specify the exact notation that will be used. For simplicity, we start with values of U_i⃗ defined on an equispaced grid. More specifically, we define a point U_i⃗=(i_1/m, …, i_q/m)^'∈ℝ^q where each i_k∈{0,1,…,m}, k=1,…,q  for some positive m. Note that the total sample size is n=(m+1)^q. This assumption ensures that m=o(n) as n→∞.   In the model (<ref>), the multivariate index is i⃗=(i_1,…,i_q)^'. Throughout this article, we will use bold font for all multivariate indices and a regular font for scalar ones.  We will say that two multivariate indices i⃗^1=(i^1_1,…, i^1_q) ≤i⃗_2=(i^2_1,…,i^2_q) if i^1_k≤ i^2_k for any k=1, …,q; the relationship between i⃗^1 and i⃗^2 is that of partial ordering.  For convenience, we denote a q-dimensional vector n⃗=(m+1,…,m+1)^' and a q-dimensional vector 0⃗=(0,…,0)^'.  The l^2 norm of a vector will be denoted ||·||_2. To define the bins we use for aggregating observations, we start with defining J=⌊1/qlog_2n^3/4⌋ and T=2^J. We split the jth edge of the q-dimensional cube [0,1]^q into T intervals of the type (l_j-1/T,l_j/T] where l_j=1,2,…,T for any j=1,2,…,q. All of the n observations are split into bins in the following way: all of Y_i⃗ such that the jth coordinate of the corresponding U⃗_i⃗ belongs in an interval (l_j-1/T,l_j/T], l_j=1,2,…,T,  are assigned to the bin with the multivariate index l⃗=(l_1,…,l_q)^'. We will denote the l⃗th bin D_l⃗. For convenience of notation, we denote a q-dimensional vector all of whose coordinates are equal to 1/T 1⃗/⃗T⃗=(1/T,…,1/T)^', a q-dimensional vector consisting of 1's as 1⃗=(1,…,1)^', and T⃗=(T,…,T)^'. We also denote l⃗/T=(l_1/T,…,l_q/T)^'. Clearly, the total number of such bins is V=T^q; note that, due to selection of J and T, the total number of bins V≍ n^3/4. Also, we define an approximate number of observations in each bin κ=n/V; clearly, κ≍ n^1/4.  Let us denote η_l⃗ the median of all ρ_i⃗ such that the corresponding observation Y_i⃗ belongs in the l⃗th bin D_l⃗; also, we denote the expectation of the sample median η_l⃗ b_l⃗:= η_l⃗. 



In order to approximate the median of observations in each multivariate bin with a normal random variable, we will have to develop a multivariate median coupling inequality. What follows is a very brief general discussion of median coupling; for more details, see <cit.>. Note that our first step will have to be estimation of the expectation of the median η_l⃗ for l⃗th bin D_l⃗. Let h(x) be the density function of ρ_i⃗, 0⃗≤i⃗≤n⃗. Let 𝒳_1,…,𝒳_n be independent random variables with the density h(x). We will need the following assumption on the density h; this assumption comes from <cit.> but is given here in full for convenience. 

∫_-∞^0h(x)=1/2, h(0)>0, and h(x) is locally Lipschitz at x=0. The Lipschitz condition at zero for h(x) means that there exists a constant C>0 such that | h(x)-h(0)|≤ C| x | in an open neighborhood of 0. 

Note that the assumption that h(0)>0 guarantees uniqueness of the median of the distribution and the asymptotic normality of the sample median; see e.g. <cit.> p. 483.  



In order to approximate a median of the l⃗th bin, we will use a coupling inequality of <cit.>. We only give its statement here without a proof.  

Let 𝒳_1,…,𝒳_n be iid random variables with the density function h(x) that satisfies the Assumption (<ref>) while Z is a standard normal random variable. We assume that n=2k+1 for some integer k≥ 1. Then, there exists a mapping X̃(Z):→ such that the distribution law of X̃_med  L(X̃_med(Z))= L(𝒳_med) and 

    |√(4n)h(0)X̃_med-Z|≤C/√(n)+C/√(n)|√(4n)h(0)X̃_med|^2

when |X̃_med|≤ where C,>0 depend on the density h but not on n. 
 
The bound given in the Theorem (<ref>) can also be expressed in terms of Z as follows. 

Under the assumptions of Theorem (<ref>), the mapping  X̃_med(Z) in Theorem (<ref>) also satisfies 

    |√(4n)h(0)X̃_med(Z)-Z|≤C/√(n)(1+| Z|^2)

when | Z|≤√(n) where C, >0 do not depend on n.


Note that the assumption that the number of observations n is odd has been made for convenience. The Remark 1 on p. 2059 of <cit.> shows that it can be dispensed by redefining the median as 𝒳_med=𝒳_(k)+𝒳_(k+1)/2 when n=2k and following a similar argument. In the future, we will always assume that the number of observations whose median is considered is odd for simplicity. Also, we can say that the Theorem (<ref>) lets us approximate each bin median with a normal random variable that has the mean f(l⃗/T)+b_l⃗ and the variance 1/4κ h^2(0).

A more general situation that we described briefly earlier would be if the intercept a and/or the median of X_i⃗ was not equal to zero. If that was the case, we would be looking at approximating the median of all observations Y_i⃗ from l⃗th bin with a normal random variable that has the mean f(l⃗/T)+b_l⃗+A and the variance 1/4κ h^2(0). 

However, an interesting question remains open. What kind of distributions of X_i⃗ and ξ_i⃗ will result in the density h(x) of ρ_i⃗ that is going to satisfy the Assumption (<ref>)? The following simple proposition identifies a reasonably wide range of distributions of X_i⃗ and ξ_i that will satisfy this Assumption. This proposition takes the form of a convenient sufficient condition. However, before stating it, we need to define a family of elliptical distributions.

A random vector X=(X_1,…,X_p)^' is said to have an elliptical distribution if its characteristic function ϕ(t) can be expressed as ϕ(t)=e^it^'μψ(t^'Σ t) where μ is the p× 1 vector and p× p matrix Σ=AA^' for some matrix A. The function ψ is called characteristic generator of X. 

Some examples of elliptical distributions are normal, t-distribution, Laplace distribution, Cauchy distribution, and many others when q=1 and their multivariate analogs when q>1. In general, elliptical distributions may not have a density; however, when they do have one, it has the form 

    f(X)=c_n/√(|Σ|)g_n((X-μ)^'Σ^-1(X-μ))
 where the function g_n is called a density generator, c_n is a normalizing constant, and μ is the median vector (which is equal to the mean if the latter exists). It is a common practice to refer to an elliptical distribution with a density function as E_n(μ,Σ, g_n). The following simple property of elliptical distributions with a density function is easy to establish and so it is given here without a proof for brevity. For detailed discussion, see e.g. <cit.>.  

Let X∼ E_n(μ,Σ, g_n). Let B be an r× q matrix while b∈^r. Then, any affine transformation of X is an elliptical distribution as well:

    b+BX∼ E_r(b+Bμ,BΣ B^',g_r).


With the above in mind, we can now state our Proposition. 

Assume that the density of random errors ξ_i⃗ h_1(x) is Lipschitz continuous at zero, symmetric around zero and strictly positive in an open neighborhood of zero. We also assume that X_i⃗ has an elliptic distribution with a density function that is continuous in an open neighborhood of zero. Moreover, without loss of generality (as discussed earlier) we assume that the median of X_i⃗ is zero. Then, the Assumption (<ref>) is satisfied.
 

Due to Lemma (<ref>), X_i⃗^'β⃗ has a univariate elliptical distribution with the median equal to zero. Clearly, the distribution of X_i⃗^'β also has a density function. The form of the density function of any  elliptical distribution given in (<ref>) clearly implies that, if the median of X_i⃗^'β⃗ is zero, the distribution of X_i⃗^'β⃗ is symmetric  around zero. Let us denote h_3(x) the distribution of X_i⃗^'β. Since X_i⃗ and ξ_i are independent, the distribution of X_i⃗^'β⃗+ξ_i⃗ is a convolution of the two densities: h(x)=∫_-∞^∞h_1(x-u)h_3(u) du.  It follows directly from the definition of Lipschitz continuity at zero that, since h_1(x) is Lipschitz continuous at zero then so is h(x); moreover, it is not hard to check that the symmetry of h_1(u) and h_3(u) around zero implies the symmetry of their convolution h(x) around zero as well. Since h(x) is symmetric around zero, we have immediately that ∫_-∞^0h(x)=1/2. Finally, since the density h_3(x) is continuous in an open neighborhood of zero and the corresponding median is unique, there exists an open neighborhood of zero where h_3(x)>0. This, together with strict positivity of h_1(x) in some open neighborhood of zero, guarantees that h(0)=∫_-∞^∞h_1(x)h_3(x) dx>0.  


Note that Proposition (<ref>) covers a wide variety of distributions of X_i⃗ and ξ_i⃗. Concerning X_i⃗, distributions such as multivariate normal, multivariate t-distribution, multivariate Laplace, multivariate Cauchy, and, in general, any symmetric stable distribution, are elliptical distributions with a density. It is also of interest that, when the index of stability less than 1, which is the case for Cauchy and Levy distributions, there is no finite mean yet the resulting distribution of X_i⃗ still satisfies requirements of the Proposition (<ref>). The requirements for the distribution of ξ_i⃗ also do not contain any moment requirements; note that, for example, a (univariate) Cauchy distribution for ξ_i⃗ satisfies assumptions of the Proposition (<ref>).






 §.§ Wavelet procedure for a binned semiparametric model




We will start with a brief introduction into the multivariate wavelet bases. First, let {ϕ,ψ} be a pair of univariate father and mother wavelets. We need to assume that both of them are compactly supported on [0,1] and ∫ϕ=1. As is known, translation and dilation of ϕ and ψ generates an orthonormal wavelet basis in L_2[0,1]; thus, in order to obtain a convenient periodized wavelet basis, we have ϕ_j,k^p(t)=∑_l=-∞^∞ϕ_j,k(t-l) and ψ_j,k^p(t)=∑_l=-∞^∞ψ_j,k(t-l) where ϕ_j,k(t)=2^j/2ϕ(2^jt-k), and ψ_j,k(t)=∑_l=-∞^∞2^j/2ψ(2^jt-k). The primary resolution level j_0 should be selected large enough to ensure that the support of the scaling functions (father wavelets) and mother wavelets at level j_0 does not cover the whole [0,1]. We also assume r regularity of our wavelet system for some positive r. This means that the first ⌊ r⌋ moments of the wavelet function ψ are equal to zero. The periodized wavelets generate a curtailed multiresolution ladder:

    V_0^p⊂ V_1^p⊂ V_2^p⊂⋯

where spaces V_j^p are spanned by ϕ_j,k^p(t). As in any regular multiresolution analysis, we have W_j^p⊕ V_j^p=V_j+1^p where the space W_j^p is spanned by ψ_j,k^p. In the future, we will suppress the superscript p from the notation for convenience. An orthonormal wavelet basis has an associated orthogonal Discrete Wavelet Transform (DWT) whose function is to transform sampled data into the wavelet coefficients. A square integrable function f on [0,1] can be expanded into a wavelet series

    f(t)=∑_k=1^2^j_0θ_j_0,kϕ_j_0,k(t)+∑_j=j_0^∞∑_k=1^2^jθ_j,kψ_j,k(t)

where θ_j_0,k=⟨ f,ϕ_j_0,k⟩ and θ_j,k=⟨ f,ψ_j,k⟩ are the wavelet coefficients of f. 
We will use a one-dimensional orthogonal wavelet basis to define a q-dimensional one. There are a number of ways to construct such a basis based on the one-dimensional one; we will use the one that is based on a tensor product construction and preserves a multiresolution analysis (MRA) in a q-dimensional space. By using q univariate orthogonal MRA's

    V_0,(i)⊂ V_1,(i)⊂ V_2,(i)⊂⋯⊂ L_2[0,1],

i=1,2,…,q, we can define a q-dimensional multiresolution analysis

    _0⊂_1⊂_2⊂⋯⋯ L_2[0,1]^q

in which _j=⊗_i=1^qV_j,(i)⊂ L_2[0,1]^q. The resulting q-dimensional multiresolution analysis corresponds to, first, one q-variate scaling function ϕ(u⃗)≡ϕ(u_1,…,u_q)=∏_i=1^qϕ_(i)(u_i) where ϕ_(i) is an ith copy of the father wavelet, and, second,  2^q-1 q-variate wavelets

    ψ^i(u⃗)≡ψ^i(u_1,…,u_q)=∏_i=1^qξ_(i)(u_i),
 
i=1,2,…,2^q-1 where ξ_(i) is either ϕ or ψ but not all of ξ_(i) are equal to the father wavelet ϕ. Now, let k⃗=(k_1,…,k_q)∈^q be the k-dimensional lattice. To complete the description of our notation, we also introduce rescaled and translated versions ϕ_j_0,k⃗(u⃗)=2^j_0q/2∏_m=1^qϕ_(m)(2^j_0u_m-k_m) and ψ_j,k⃗^i(u⃗)=2^jq/2∏_m=1^qξ_(m)(2^ju_m-k_m) where ξ=ϕ or ψ but not all ξ=ϕ. Finally, any function f∈ L_2[0,1]^q can be represented as 

    f(u⃗)=∑_1⃗≤k⃗≤2⃗^j_0θ_j_0,k⃗ϕ_j_0,k⃗(u⃗)
       +∑_j≥ j_0∑_1⃗≤k⃗≤2⃗^j∑_i=1^2^q-1θ_j,k⃗^iψ_j,k⃗^i(u⃗)

For a more detailed discussion of multivariate wavelet bases see, for example,  <cit.> and <cit.>. 

At this point, we can give a detailed description of our estimator of the functional component f. The first step is to bin observations Y_i⃗ according to values of coordinates of U_i⃗ as described in Chapter (<ref>). Then, the sample median must be computed within each bin. We will use the notation g(l⃗/T)=f(l⃗/T)+b_l⃗ where b_l⃗ is the median of errors in the l⃗th bin as described earlier. We will denote medians of observations in l⃗th bin Q_l⃗, 1⃗≤l⃗≤T⃗. Our next step consists of applying a discrete wavelet transform to all of the medians Q_l⃗. In the multivariate case that we consider, a discrete wavelet transform is a tensor. After the application of the discrete wavelet transform, we can describe the transformed data as 

    U=(y_j_0,1,…,y_j_0,2^j_0,y_j_0,1^1,…,y_j_0,2^j_0^1,…,y_J-1,1^2^q-1,…,y_J-1,2^J-1^2^q-1).

Here, y_j_0,r, r=1,…,2^j_0 are the gross structure terms at the lowest resolution level, while y_j,k^i, i=1,…,2^q-1, j=j_0,…,J-1, k=1,…,2^j are empirical wavelet coefficients at level j corresponding to the wavelet i that represent fine structure at scale 2^j. The empirical wavelet coefficients can be written as 

    y_j,k⃗^i=θ̆_j,k⃗^i+_j,k⃗^i+1/2h(0)√(n)z_j,k⃗^i+ξ_j,k⃗^i

where θ̆_j,k⃗^i are the discrete wavelet coefficients of g(l⃗/T), _j, k⃗^i are deterministic errors, z_j,k⃗^i are iid N(0,1) and ξ_j,k⃗^i are stochastic errors. Later, we will show that both deterministic errors  _j,k⃗^i and stochastic errors ξ_j,k⃗^i are negligible in a certain sense. Ignoring them, we end up with 

    y_j,k⃗^i≈θ̆_j,k⃗^i+1/2h(0)√(n)z_j,k⃗^i

which is essentially an idealized sequence model with the noise level σ=1/2h(0)√(n). As a next step, we propose a multivariate generalization of the BlockJS procedure of <cit.> and apply it to the empirical coefficients y_j,k⃗^i as if they were generated by (<ref>). At each resolution level j the empirical wavelet coefficients y_j,k⃗^i  are grouped into nonoverlapping blocks of length L. We define each block as B_j,u^i  consisting of observations y_j,k⃗^i  such that, for any choice of ith wavelet and jth resolution level, only (u-1)L+1≤ k_s≤ uL are included where s=1,2,…,q and j=j_0,…,J-1.  We also define S_j,u,i^2≡∑_k⃗∈ B_j,u^i(y_j,k⃗^i)^2 the sum of squared empirical wavelet coefficients included in uth block. Let ĥ^2(0) be an estimator of the squared value of the density function h at the point zero. The following shrinkage rule is then applied to each block B_j,u^i:

    θ̂_j,k⃗^i=(1-λ_*L/4ĥ^2(0)nS_j,u,i^2)_+y_j,k⃗^i

where λ_* is the solution of the equation λ_*-logλ_*=3 and 4nĥ^2(0) is present due to the fact that the noise level is equal to σ=1/2h(0)√(n). For the gross structure terms, we define an estimator θ̂_j_0,k⃗=y_j_0,k⃗. Now, we reconstruct the estimate of the function g at the points l⃗/T by applying the inverse discrete wavelet transform to the shrunk empirical wavelet coefficients. In other words, the entire function g(u⃗) can be estimated for any u⃗=l⃗/T, 1⃗≤l⃗≤T⃗ as 

    ĝ(u⃗)=∑_1⃗≤k⃗≤2⃗^j_0θ̂_j_0, k⃗ϕ_j_0,k⃗(u⃗)+∑_j=j_0^J-1∑_1⃗≤k⃗≤2⃗^j∑_i=1^2^q-1θ̂_j,k⃗^iψ^i_j,k⃗(u⃗).

The last remaining step is to estimate the median b_l⃗ in order to obtain an estimate of the function f. The following procedure is employed for that purpose. Recall that the jth edge of the q-dimensional cube [0,1]^q has been split into T intervals; each of these intervals contained ⌊m+1/T⌋ observations. We begin with splitting each of these T intervals for a jth edge of a q-dimensional cube in two in such a way that the smaller half contains ⌊m+1/2T⌋ observations. Now, we define a set of q-dimensional bins in the following way. Each bin in this set consists of smaller halves of all one-dimensional intervals for each of the q dimensions. Note that the cardinality of that set of bins will still be V=T^q. Next, we define Q_l⃗^* to be the median of the l⃗th “halfbin" with Q_l⃗ being the median of all observations from the corresponding l⃗th “full" bin. Then, the median expectation is defined to be 

    b̂_l⃗=1/V∑_1⃗≤l⃗≤T⃗[Q_l⃗^*-Q_l⃗].

This lets us define 

    f̂_n(u⃗)=ĝ(u⃗)-b̂_l⃗=∑_1⃗≤k⃗≤2⃗^j_0θ̂_j_0,k⃗ϕ_j_0,k⃗(u⃗)+∑_j=j_0^J-1∑_1⃗≤k⃗≤2⃗^j∑_i=1^2^q-1θ̂_j,k⃗^iψ^i_j,k⃗(u⃗)-b̂_l⃗.






§ ADAPTIVITY OF THE PROCEDURE


We study the theoretical properties of our procedure over the Besov spaces in ^q and over a suitable class of distributions of X_i⃗ and ξ_i⃗ as defined in (<ref>). Besov spaces in ^q make up a very rich class of spaces that incorporates functions of very significant spatial inhomogeneity and includes Hölder and Sobolev spaces as special cases. Our discussion of Besov spaces is necessarily brief; for details, see for example <cit.>. Since multivariate Besov spaces are not used in the statistical literature very often, we will briefly describe them first. Let e_i=(δ_i1,…,δ_iq)^' where δ_ij=I(i=j) is the ith unit vector. We define the first difference of the function f in the ith direction as Δ_i,hf(U)=f(U+he_i)-f(U) and the second difference as Δ_i,h^2f(U)=Δ_i,h(Δ_i,hf(U))=f(U+2he_i)-2f(U+he_i)+f(U). Let α>0 be a positive integer and 𝔞=⌊α⌋. We denote the increment in each direction i h where | h|<1. We also need to define g_i,h=(0,1)^i-1×(0,0∨(1-2h))× (0,1)^q-i. The so-called Besov norm in the direction i is then defined as 

    f_b^α_i,s,t=(∫_0^1| h|^(α-𝔞)t-1Δ^2_i,h(∂^𝔞/∂ u_i^𝔞f)^t_L_s(g_i,h) dh)^1/t

for t<∞  and 

    f_b^α_i,s,∞=sup_0≤ h ≤ 1{| h|^α-𝔞Δ^2_i,h(∂^𝔞/∂ u_i^𝔞f)_L_s(g_i,h)}

otherwise. The Besov norm is now defined as the sum of the L_s norm of the function f and Besov norms in all possible directions:

    f_B^α_s,t=f_L_s[0,1]^q+∑_i=1^qf_b^α_i,s,t.
 A Besov class, sometimes also called a Besov ball,  can be defined as B^α_s,t(M)≐{f|f_B^α_s,t≤ K} for some positive K. 

Besov norm of a function can also be characterized in terms of its wavelet coefficients. For a fixed primary resolution level j_0, the Besov sequence norm of the wavelet coefficients of a function f can be defined the following way. First,  let θ_j_0 be a vector of the father wavelet coefficients at the primary resolution level j_0, θ_j be the vector of wavelet coefficients at level j, and w=α+q(1/2-1/s)>0. Then, we define first

    ||θ_j_0||_s=(∑_k⃗|θ_j_0,k⃗|^s)^1/s

and

    ||θ_j||_s=(∑_k⃗∑_i=1^2^q-1|θ_j,k⃗^i|^s)^1/s.


With the above in mind, a sequence norm for a function f∈ B^α_s,t can be defined as as 

    ||f||_B^α_s,t=||θ_j_0||_s+(∑_j=j_0^∞(2^w||θ_j||_s)^t)^1/t.

We know that the Besov sequence norm is equivalent to the Besov function norm defined in (<ref>) and, therefore, the Besov class can also be defined as B^α_s,t(M)={f;f_B^α_s,t≤ M}; for details, see e.g. <cit.>. It is also known that, in case of Gaussian noise, the minimax risk of estimating f over the Besov body B^α_s,t(M), 

    R^*(B^α_s,t(M))=inf_f̂sup_f∈ B^α_s,t f̂-f^2_2,

converges to zero at the rate of n^-2α/2α+q as n→∞; see e.g. <cit.>. The following theorem shows that our estimator achieves optimal global adaptation for a wide range of multivariate Besov classes B^α_s,t(M). Moreover, this adaptation is also uniform over a range of distributions of the design vector X_i⃗ and the error term _i⃗. To state this theorem properly, we need to define appropriate classes of distributions of X_i⃗ and ξ_i⃗. We begin with the following Assumption on the density of ρ_i⃗ h(x). This Assumption guarantees the existence of some low order moment of ρ_i⃗. 

∫ |x|^𝒜h(x) dx<∞ for some 𝒜>0. 

For any 0<_1<1, _2>0, we define the class of densities H__1,_2 by 

    H__1,_2={h:∫_-∞^0h(x) dx=1/2,_1≤ h(x)≤1/_1,| h(x)-h(0)|≤| x|/_1| x|≤_2}.

We will also need another, more narrow class of densities that is defined as H≡ H__1,_2,_3,_4 for some 0<_1<1, _i>0, i=2,3,4 where 

    H={h:h∈ H__1,_2, | h^(3)(x)|≤_4| x|≤_3∫| x|^_3h(x) dx<_4}.

One simple way to ensure that h(x)∈ H is to require, first, that  densities h_1(u) and h_2(u) satisfy assumptions of the Proposition (<ref>). In addition, we also have to require the existence of a finite moment of some order λ>0 for the density h_1(u) and for the distribution of X_i⃗ plus the existence of the finite third derivative for h_1(u) in a small open neighborhood of zero. Then, the moment assumption is satisfied because, for any 0<π≤ 1, and random variables X_1,…,X_n, we have |X_1+⋯+X_n|^π≤ |X_1|^π+⋯+ |X_n|^π.  Moreover, differentiability of h_1(u) ensures the same property for the convolution h(u). With the above in mind, we can define the class G_1={h_1(u): h_1(u) , h_1(u)  u=0, h_1(u)>_1 |u|<_2, |h^(3)_1(u)|≤_3 |u|≤_4, ∫ |u|^_3h_1(u) du <_4}. As for X_i⃗, we define G_2 as a class of all p-dimensional elliptical distributions with finite moments of a small order _3>0. With these definitions in mind, we can now formulate the global adaptivity result. 


 Suppose the wavelet ψ is r-regular. Define d=min(α-q/s,1). Then the estimator f̂_n defined in (<ref>) satisfies, for s≥ 2, α≤ r, and 3d/2q>2α/2α+q,
 
    sup_h_1∈ G_1,h_2∈ G_2sup_f∈ B^α_s,t(M) f̂_n-f^2_2≤ Cn^-2α/2α+q,

 and for 1≤ s <2, α≤ r and 3d/2q>2α/2α+q, 
 
    sup_h_1∈ G_1,h_2∈ G_2sup_f∈ B^α_s,t(M) f̂_n-f^2_2≤ Cn^-2α/2α+q(log n)^2-s/[s(2α+q)+2(1-q)].

 
Effectively, we show that the estimator attains the optimal rate of convergence over a wide range of Besov classes for f and a large collection of the unknown error distributions for _i, as well as for a large collection of distributions of the design vector X_i⃗. Note also that when q=1, these rates are reduced to those in the one dimensional case, n^-2α/2α+1 and n^-2α/2α+1(log n)^2-s/[s(2α+1), respectively; see, e.g. <cit.>. 

Since functions belonging to Besov classes B^α_s,t(M) are highly spatially inhomogeneous, local adaptivity at an arbitrary point u⃗_0∈ [0,1]^q should also be investigated. To measure such a spatial adaptivity, the local mean squared risk 

    R(f̂_n(u⃗_0), f(u⃗_0))≡(f̂_n(u⃗_0)-f(u⃗_0))^2

is used. The precise way of measuring the local smoothness of the function f(u⃗) at a given point u⃗=u⃗_0 is by the use of its local Hölder smoothness index, that is a characteristic of the point u⃗_0. From the technical viewpoint, it is more straightforward to assume that the function f on the entire cube [0,1]^q belongs to a Lipschitz class Λ^α(M) with the same smoothness index α for every point. To define such a class, we define first for a multivariate index i⃗ |i⃗|=i_1+…+i_q. Next, let us select a constant M>0, and, for a q-dimensional index i=(i_1,…, i_q), define i(l)={ i:| i|=i_1+…+i_q=l}. Then, for any function f:ℝ^q→ℝ, the mixed partial derivative of order l, D^ i(l)f/∂ u_1^i_1…∂ u_q^i_q is defined for all i such that | i|=l. With partial derivative thus defined, the Lipschitz class Λ^α(M) consists of all functions f(u⃗):[0,1]^q→ℝ such that |D^ i(l)f(u⃗)|≤ M for l=0,1,…,⌊α⌋ and |D^ i(⌊α⌋)f(v)-D^ i(⌊α⌋)f(w)|≤ M||v-w||^α^' with α^'=α-⌊α⌋. The Hölder smoothness index α will be used to measure local smoothness of the function f. Then, the following theorem shows that our estimator achieves optimal local adaptation uniformly over the same families of distributions of X⃗_i and _i as before. 

Assume that f∈Λ^α(M) on [0,1]^q. Suppose the wavelet ψ is r-regular, r≥α>q/6, and u⃗_0∈ (0,1)^q is a fixed point. Then, the estimator f̂_n defined in (<ref>), satisfies

    sup_h_1∈ G_1,h_2∈ G_2sup_f∈Λ^α(M) f̂_n(u⃗_0)-f(u⃗_0)^2_2≤ C(log n/n)^2α/2α+q.



Note that the smoothness index of the Besov space scale that we used must remain above the ratio q/6. This can be explained by the additional difficulty of uniform estimation over a wide ranging scale of Besov spaces in higher dimensions. In other words, the higher the dimension is, the smoother must the functional scale remain to enable a uniform in mean squared error risk estimation over it.




§ DISCUSSION AND A FUTURE RESEARCH


We constructed an estimator of the functional component f in the multivariate partial linear model (<ref>) under relatively simple conditions to stress the basic ideas. For example, the assumption of equispaced grid on each edge of the q-dimensional cube [0,1]^q (standard fixed design) reflects the situation where the data points U_i⃗ are obtained experimentally.  A more realistic assumption would be to allow for non-equispaced grid; our approach, however, can be modified to account for such a possibility. 

A simple option would be to ignore the fact that the grid is non-equispaced, and keep using the same binning approach as we defined earlier. In this case, bins will have different numbers of observations and the resulting medians will have different variances. This implies that a multivariate wavelet smoothing procedure that accounts for heteroskedasticity will have to be used. To the best of our knowledge, such a procedure has been proposed so far only in the univariate case in <cit.>. However, a generalization of their procedure to the multivariate case is fairly straightforward. The other possible approach to this problem would be to bin observations in such a way that each bin contains the same number of observations; this would imply that, for each edge of the cube [0,1]^q, it will be necessary to adjust length of intervals according to the density of the design. This will result in irregularly spaced but homoskedastic data. There are a number of approaches that have been proposed for handling of irregularly spaced data in the wavelet shrinkage context; some of the better known ones are <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. Almost all of them treat the univariate case only; however, a multivariate wavelet thresholding procedure that can adapt to local changes in the smoothness of the regression function and to the distribution of the design has been proposed in <cit.>. That procedure was proposed for a multivariate random design. However, since the true density of the design is not usually known, in either case one can simply construct an orthonormal wavelet basis of the space L_2[0,1]^q(μ_n) where μ_n is the measure corresponding to the empirical distribution function of U_i⃗. 

Another issue of importance is the fact that the multivariate curve estimation itself is sometimes regarded as problematic since, even for relatively small number of dimensions q, the “curse of dimensionality" begins to manifest itself and minimax rates of convergence becomes unacceptably slow. In practice, however, the true complexity of the multivariate curve may be considerably lower simply due to the fact that the real multivariate data are not often truly isotropic. As an example, <cit.> remarks in Chapter 2 that “Multivariate data in ^d are almost never d-dimensional. That is, the underlying structure of data in ^d is almost always of dimensions lower than d".  This would imply that it makes sense to model the function f in our model (<ref>) as a member of some anisotropic functional class, e.g. the anisotropic Besov classes B^α⃗_s⃗,t(M) where parameters α⃗ and s⃗ are now both q-dimensional vectors. Such classes were rigorously defined in <cit.> and the multivariate wavelet thresholding procedure for a multivariate anisotropic Gaussian white noise model was proposed in <cit.>; for a more recent approach to this topic see also <cit.>. The next sensible step is to allow for a function f to belong to an anisotropic functional class and to obtain an estimator of f that is, again, adaptive to a range of such functional classes and robust with respect to a range of distributions of the design vector X_i⃗ and of the distributions of random errors ξ_i⃗. The results of this ongoing work will be reported elsewhere. 



§ APPENDIX


In order to prove our adaptation results, we need to start with the following Proposition that provides us with an expansion of the median of binned observations. That expansion shows that, up to small stochastic and deterministic errors, that median has an approximately normal distribution. An important fact that we need to use in order to prove this Proposition is that any Besov ball B^α_s,t(M) can be embedded into a Hölder ball with the smoothness index d=min(α-q/s,1); for details see e.g. <cit.>. We also remark here that, in this section, we use the notation C for a generic positive constant that can be different from one line to another. 

Let the function f∈ B^α_s,t(M) and d≐min(α-q/s,1). The median Q_l⃗ of observations that belong to the l⃗th bin can be written as 

    √(κ)Q_l⃗=√(κ)f(l⃗/T)+√(κ)b_l⃗+1/2Z_l⃗+_l⃗+ζ_l⃗

where 

  * Z_l⃗∼ N(0,1/h^2(0));

  * _l⃗ are constants such that |_l⃗|≤ C√(κ)q^d/2T^-d;

  * ζ_l⃗ are independent random variables such that for any r>0 

    |ζ_l⃗|^r≤ C_rκ^-r/2+C_rκ^r/2q^dr/2T^-dr

where C_r is a positive constant that depends on r only; moreover, for any a>0

    P(|ζ_l⃗|>a)≤ C_r(a^2κ)^-r/2+C_r(a^2T^2d/κ q^d)^-r/2.



 

In this proof, we denote the cdf of the standard normal distribution Φ(·). Define Z_l⃗=1/h(0)Φ^-1(G(η_l⃗)) where G is the distribution of the median η_l⃗. Due to Theorem (<ref>), we know that the rescaled median of errors √(4κ)η_l⃗ can be well approximated by a mean zero normal random variable with the variance equal to 1/h^2(0). Next, we  define 

    _l⃗=√(κ) Q_l⃗-√(κ)f(l⃗/T)-√(κ)b_l⃗
       = {√(κ)Q_l⃗-√(κ)f(l⃗/T)-√(κ)η_l⃗}.

What we have in the above is the deterministic component of the approximation error due to binning. Clearly, for the l⃗th bin D_l⃗, we have 

    min_u_i⃗∈ D_l⃗[f(u_i⃗)-f(l⃗/T)]
       ≤ Q_l⃗-η_l⃗-f(l⃗/T)≤max_u_i⃗∈ D_l⃗[f(u_i⃗)-f(l⃗/T)]

Since the function f is in a Hölder ball with the smoothness index d=min(α-q/s,1),  we have 

    |_l⃗|≤√(κ) | Q_l⃗-f(l⃗/T)-η_l⃗|
       ≤√(κ)max_u_i⃗∈ D_l⃗| f(u_i⃗)-f(l⃗/T)|≤ C√(κ)q^d/2T^-d


Now, it becomes necessary to characterize the random error of our approximation. First, we define ζ_l⃗=√(κ)Q_l⃗-√(κ)f(l⃗/T)-√(κ)b_l⃗-_l⃗-1/2Z_l⃗. Note that ζ_l⃗=0 and this random error can be represented as the sum of two components, ζ_1l⃗=√(κ)Q_l⃗-√(κ)f(l⃗/T)-√(κ)η_l⃗-_l⃗ and ζ_2l⃗=√(κ)η_l⃗-√(κ)b_l⃗-1/2Z_l⃗. The first component ζ_1l⃗ represents the random error resulting from the binning of observations while ζ_2l⃗ is the error resulting from approximation of the median of random errors with a normal random variable. First, the error ζ_1l⃗ is bounded due to (<ref>) as |ζ_1l⃗|≤ C√(κ)q^d/2T^-d. Next, using the Corollary (<ref>), we can bound the absolute value of the second random error term as |ζ_2l⃗|≤C/κ^1/2(1+| Z_l⃗|^2) when | Z_l⃗|≤√(κ) for some >0.  Thus, for any fixed r≥ 0, 

    |ζ_2l⃗|^r= |ζ_2l⃗|^rI(| Z_l⃗|≤√(κ))+ |ζ_2l⃗|^rI(| Z_l⃗| > √(κ))
       ≤ Cκ^-r/2+ |ζ_2l⃗|^rI(| Z_l⃗| > √(κ))
       ≤ Cκ^-r/2+{ |ζ_2l⃗|^2r}^r/2exp(-^2/κ)

since, due to the Mill's ratio inequality, we can quickly verify that P(| Z_l⃗| > √(κ))≤exp(-^2/κ).  

To continue, we will need to use the Assumption (A2). Using the same argument as in <cit.>, we can show that the density of the centered median η_l⃗-b_l⃗ g(x) is such that the sample centered median has any finite moments and, therefore, |√(κ)(η_l⃗-b_l⃗)|^2r≤κ^r |η_l⃗|^2r≤ D_rκ^r for some positive constant D_r that does not depend on n. This allows us to conclude that  |ζ_2l⃗|^r≤ C_rκ^-r/2 since the normal random variable Z_j has finite moments of any order. Finally,  

    |ζ_l⃗|^r≤ 2^r-1( |ζ_1l⃗|^r+ |ζ_2l⃗|^r)
       ≤ C_rκ^-r/2+C_r⃗κ^r/2q^dr/2T^-dr.
 
The inequality (<ref>) can now be obtained immediately using the Markov's inequality. 


In the proof of our adaptation results, we will assume everywhere that h(0) is known and equal to 1. This can be done because h^-2(0) can always be estimated in such a way that the difference between h^-2(0) and ĥ^-2(0) is bounded in probability as O_p(n^-δ)  for some δ>0 and, moreover, P(|ĥ^-2(0)-h^-2(0)|≥ n^-δ)=c_ln^-l for any l ≥ 1. Such an estimator can be constructed as a properly normalized sum of ordered squared difference of medians of observations Q_l⃗, 1⃗≤l⃗≤ V; to order medians, one can use, for example, a lexicographical order of their indices l⃗.. To make the notation easier, we denote two successive medians Q_2k-1 and Q_2k, using a scalar index to avoid confusion here. Then, the needed estimator of h^-2(0) will be proportional to ∑_k(Q_2k-1-Q_2k)^2. When such an estimator is constructed, one can immediately check that everywhere in proofs asymptotic properties do not change if λ^*(1+O(n^-δ)) is used instead of λ^*. The details are similar to the argument of <cit.> and are omitted for conciseness.  

The next proposition is needed to obtain a uniform bound of the mean squared error risk of estimating the expected error median. Its proof is very similar to that of Lemma 5 in <cit.> and is omitted for brevity.  

Let the expectation of the error median for l⃗th bin and its estimate be b_l⃗ and b̂_l⃗, as defined earlier in (<ref>). Then, 

    sup_h∈ H| b_l⃗+h^'(0)/8h^3(0)κ|≤ Cκ^-2,

and

    sup_h∈ Hsup_f∈ B^α_s,t(M) (b̂_l⃗-b_l⃗)^2≤ Cmax{q^dT^-2d,κ^-4}

for any index l⃗. 
 


Without loss of generality, assume that κ=2ν+1. Then, the expectation of the l⃗th median is

    η_l⃗=∫ x2ν+1/(ν!)^2H^ν(x)[1-H(x)]^ν dH(x)

where H(x) is the distribution function corresponding to h(x). For any δ>0, define a set A_δ={x:| H(x)-1/2|≤δ}. It follows from the definition of the class H  that there exists a constant δ>0 such that for some >0 we have h^(3)(x)≤1/ and ≤ h(x)≤1/ for any x∈ A_δ uniformly for all h∈ H. This property implies that H^-1(x) is well defined and is differentiable up to the fourth order for any x∈ A_δ. Now we can expand the expectation of the median into two parts:

    η_l⃗=(∫_A_δ+∫_A_δ^c)x2ν+1/(ν!)^2H^ν(x)[1-H(x)]^ν dH(x)

Since we earlier established that all of the moments of the median are finite, Q_2 goes to zero exponentially as μ→∞. Next, we find that 

    Q_1=∫_1/2-δ^1/2+δ(H^-1(x)-H^-1(1/2))(2ν+1)!/(ν!)^2x^ν(1-x)^ν dx
       =∫_1/2-δ^1/2+δ[1/2(H^-1)^”(1/2)(x-1/2)^2+(H^-1)^(4)(τ)/24(x-1/2)^4]
       ×2ν+1/(ν!)^2x^ν(1-x)^ν dx

since x^ν(1-x)^ν is symmetric around 1/2. The expression 2ν+1/(ν!)^2x^ν(1-x)^ν is the density function of the Beta(ν+1,ν+1) with the mean equal to 1/2; this, and the fact that (H^-1)^(4)(τ) is bounded uniformly for all h∈ H , implies that, in the same way as in <cit.>, Q_1=-h^'(0)/8h^3(0)κ+O(1/κ^2). 
Recall that we earlier established that (see Proposition (<ref>)) that 

    Q_l⃗=f(l⃗/T)+b_l⃗+1/2√(κ)Z_l⃗+1/√(κ)_l⃗+1/√(κ)ζ_l⃗.

In a similar way, we can write for the median of the “half" l⃗ th bin that 

    Q_l⃗^*=f(l⃗-1⃗/⃗2⃗/T)+b_l⃗^*+1/2√(ν)Z_l⃗^*+1/√(ν)_l⃗^*+1/√(ν)ζ_l⃗^*

where 1⃗/⃗2⃗ is a q-dimensional vector (1/2,…,1/2)^', b_l⃗^* is the  expected median of the errors of all observations in the “half" l⃗th bin, and where Z_l⃗^*, _l⃗^*, and ζ_l⃗^* satisfy the Proposition (<ref>). Then, the error from the median estimation b̂_l⃗-b_l⃗=1/V∑_l⃗(Q_l⃗^*-Q_l⃗)-b_l⃗ can be written as follows:

    b̂_l⃗-b_l⃗=1/V∑_l⃗(f(l⃗-1⃗/⃗2⃗/T)-f(l⃗/T))+(b_l⃗^*-2b_l⃗)
       +[1/√(ν)1/V∑_l⃗_l⃗^*-1/√(κ)1/V∑_l⃗_l⃗]
       +[1/2√(ν)1/V∑_l⃗Z_l⃗^*-1/2√(κ)1/V∑_l⃗Z_l⃗]
       +[1/√(ν)1/V∑_l⃗ζ_l⃗^*-1/√(κ)1/V∑_l⃗ζ_l⃗]
       ≡ R_1+R_2+R_3+R_4+R_5.

Due to the embedding of the Besov ball B_s,t^α(M) into the Hölder ball with the smoothness index d=min(α-q/s,1), the first term is uniformly bounded: sup_f∈ B^α_s,t(M)R_1^2≤ CT^-2d. The second term is bounded as sup_h∈ HR_2^2≤ Cκ^-4. By Proposition (<ref>), the third term is also bounded as sup_h∈ H, f∈ B^α_s,t(M)R_3^2≤ CT^-d. Since Z_l⃗^*-Z_l⃗ are always independent, we have R_4^2≤1/h^2(0)(1/κ+1/ν)1/T≤ Cn^-1.  Finally, by the Proposition (<ref>), we have R_5^2=o(n^-1). Thus, the overall bound is 

    sup_h∈ H,f∈ B^α_s,t(M) (b̂_l⃗-b_l⃗)^2≤max{T^-2d,κ^-4}.
 






 Due to the Proposition (<ref>), we can write

    1/√(V)Q_l⃗=g(l⃗/T)/√(V)+_l⃗/√(n)+Z_l⃗/2√(n)+ζ_l⃗/√(n).


Let Q be the vector of all bin medians Q_l⃗; such a vector will have the length V. Applying the discrete wavelet transform to both sides of (<ref>), we can expand the empirical wavelet coefficients y^i_j,k⃗ as

    y^i_j,k⃗=θ̆_j,k⃗^i+_j,k⃗^i+1/2h(0)√(n)z_j,k⃗^i+ξ_j,k⃗^i

where θ̆_j,k⃗^i are the discrete wavelet coefficients of g(l⃗/T)_1⃗≤l⃗≤T⃗ that are approximately equal to the true wavelet coefficients of g θ_j,k⃗^i, _j,k⃗^i are “small" deterministic approximation errors, z_j,k⃗^i are i.i.d N(0,1), and ξ_j,k⃗^i are “small" stochastic errors. If it can be assumed that _j,k⃗^i and ξ_j,k⃗^i are both negligible in some sense, we may be able to treat the model in the wavelet domain as an idealized sequence model

    y_j,k⃗^i≈θ̆_j,k⃗^i+1/2h(0)√(n)z_j,k⃗^i
 
where 1/2h(0)√(n) plays the role of the noise level.  At this point, we can define a simple estimation procedure for the function f. Some auxiliary results are necessary before stating the main result. The first of these results is needed in order to bound the difference between the true wavelet coefficients θ^i_j,k⃗ of the function f(l⃗/T)_1⃗≤l⃗≤T⃗ and discrete wavelet coefficients θ̆^i_j,k. Its proof is straightforward and is, therefore, omitted. 

Let T=2^J and define f_J(u)=1/√(V)∑_l⃗≤T⃗∑_i=1^2^q-1f(l⃗/T)ψ_J,l⃗^i(u). Then, 

    sup_f∈ B^α_s,t(M)||f_J-f||^2_2≤ C· q^dT^-2d

where d=min(α-q/s,1). Moreover, |θ̆_j_0,k⃗^i-θ_j,k⃗^i|≤ CT^-d2^-j/2 and so ∑_j=j_0^J-1∑_1⃗≤ k ≤2⃗^j(θ̆_j,k⃗^i-θ_j,k⃗^i)^2≤ CT^-2d. 
 
 Another result that we need is the following proposition that studies the risk of our proposed procedure. 
 
 Let the empirical wavelet coefficients y_j,k⃗^i=θ̆_j,k⃗^i+_j,k⃗^i+1/2h(0)√(n)z_j,k⃗^i+ξ_j,k⃗^i be as given in (<ref>) and let estimated block thresholding coefficients θ̂_j,k⃗^i be as defined in  (<ref>). Then, for some constant C>0, 
 
    ∑_k⃗∈ B_j,u^i(θ̂_j,k⃗^i-θ̆_j,k⃗^i)^2≤
       ≤min{4∑_k⃗∈ B^i_j,u([θ̆_j,k⃗^i]^2,8λ L n^-1}+6∑_k⃗∈ B^i_j,u[_j,k⃗^i]^2+C L n^-2;

 also, for any 0<τ<1, there exists a constant C_τ>0 that depends on τ only such that for all k⃗∈ B^i_j,u
 
    (θ̂_j,k⃗^i-θ̆_j,k⃗^i)^2≤ C_τmin{max_k⃗∈ B^i_j,u{(θ̆_j,k⃗^i+_j,k⃗^i)^2},Ln^-1}+n^-2+τ.

  
 
 This proof is similar to the proof of Proposition 2 in <cit.>; therefore, we only give a brief proof for the inequality (<ref>).  First, we recall that |_l⃗|≤ C√(κ)q^d/2T^-d. The discrete wavelet transform of _l⃗/√(n) in our case is equal to ^i_j,k⃗=∑_l⃗∈^q_l⃗/√(n)∫ϕ_j,l⃗ψ_j,k⃗^i. The Proposition (<ref>) suggests that ∑_j∑_k⃗∑_i[_j,k⃗^i]^2=1/n∑_l⃗∈^q_l⃗^2≤ Cq^dT^-2d for some positive constant C due to orthogonality of the discrete wavelet transform. Thus, we have for any r>0 
 
    |ξ_j,k⃗^i|^r≤ C_r(κ n)^-r/2+C_r(κ n )^-r/2q^dr/2T^-dr

and, for any a>0

    P(|ξ_j,k⃗^i|>a)≤ C_r^'(a^2κ n)^-r/2+C_r^'(a^2n T^2d/κ q^d)^-r/2.

where C_r and C_r^' are constants that do not depend on n. 
At this point, we need to use Lemma 2 of <cit.> with the number of observations being the size of each block L. For an ith wavelet function, i=1,…,2^q-1, we have the expectation of the risk over each block bounded as 

    ∑_k⃗∈ B^i_j,u(θ̂_j,k⃗^i-θ̆_j,k^i)^2≤
       ≤min{4∑_k⃗∈ B^i_j,u(θ̆_j,k⃗^i)^2,8λ Ln^-1}+6∑_k⃗∈ B^i_j,u(_j,k⃗^i)^2+
       + 2n^-1∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2I(∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2>λ L).

Denote by A the event that all |ξ_j,k⃗^i| are bounded by 1/2√(n)L, that is

    A={|2√(n)ξ_j,k⃗^i|≤ L^-1k⃗∈ B^i_j,u}.
 
Then it follows from (<ref>) that for any r≥ 1, the probability of a complement of A is

    P(A^')≤∑_k⃗∈ B^i_j,uP(| 2√(n)ξ_j,k⃗^i| >L^-1)
       ≤ C_r^'(L^-2κ)^-r/2+C_r^'(L^-2T^d/κ q^d)^-r/2.


Thus, we have 

    D= ∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2I(∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2>λ L)
       =∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2I(A∩∑_k⃗∈ B^i_j,u(z_j,k^i+2√(n)ξ_j,k^i)^2>λ L)
       +∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2I(A^c∩∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2>λ L)≡ D_1+D_2.

Recall that (x+y)^2≤ 2x^2+2y^2 for any x and y. At the next step, we have to use the inequality from Lemma 3 of <cit.> (with the value λ̃=λ L-λ-1/L), and Hölder's inequality to obtain 

    D_1= ∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2I(A∩∑_k⃗∈ B^i_j,u(z_j,k⃗^i+2√(n)ξ_j,k⃗^i)^2>λ L)
       ≤ 2 ∑_k⃗∈ B^i_j,u[z_j,k⃗^i]^2I(∑_k⃗∈ B^i_j,u[z_j,k⃗^i]^2>λ L-λ-1)
       +8n ∑_k⃗∈ B^i_j,u[ξ_j,k⃗^i]^2I(∑_k⃗∈ B^i_j,u[z_j,k⃗^i]^2>λ L-λ-1)
       ≤ 2(λ L-λ-1)e^-L/2(λ-(λ+1)L^-1-log (λ-(λ+1)L^-1)-1)
       +8n∑_k⃗∈ B^i_j,u( [ξ_j,k⃗^i]^2v)^1/v(P(∑_k⃗∈ B^i_j,u[z_j,k⃗^i]^2>λ L-λ-1))^1/ω

where v,ω>1 and 1/v+1/ω=1. Recall that κ=n^1/4 and choose 1/ω=1-1/4=3/4. This lets us conclude that D_1≤ CLn^-1. Arguing similarly, we conclude that D_2≤ n^-1 and so the final inequality is obtained.
 

Note that the tail probability P(|ξ_j,k⃗^i|>a) must decay faster than any polynomial in n to ensure that the contribution of ξ_j,k⃗^i to the squared risk of the proposed procedure is negligible compared to that of z_j,k⃗^i. Recall that κ=n^1/4. Then, T^2d/κ=n^6d-q/4q and we have to require that 6d-q>0, or d=min(α-q/s,1)>q/6. Since d characterizes smoothness of the Hölder ball the Besov ball is embedded into, it may be a good idea to describe this requirement in terms of the original smoothness indicator α. Note that (see Remark (<ref>)) that, due to approximation error over multivariate Besov spaces, we must have 3d/2q> 2α/2α+q. To guarantee that d>q/6, we may require that 4α q/3(2α+q)> q/6 which is equivalent to α> q/6. This is the origin of the lower bound on α in the statement of Theorems (<ref>) and (<ref>). 
   

Proof of the Theorem (<ref>)
First, note that 

    f̂_n-f^2_2≤ 2 ĝ_n-g^2_2+2 (b̂_l⃗^2-b_l⃗)^2.

By selecting κ≍ n^1/4, we ensure that (b̂_l⃗^2-b_l⃗)^2=o(n^-2α/2α+q) and so we need only to focus on bounding ĝ_n-g^2_2. Note that, if the intercept a and/or the median of the vector X_i⃗ in the model (<ref>) are non-zero, an appropriate term in the model can be estimated at the rate of n^-1=o(n^-2α/2α+q). Using the notation of Section (<ref>), let A=a+(med X_i⃗)^'β and A_n an asymptotically normal √(n) convergent estimator of A. Such an estimator can be easily obtained as in <cit.>. In that case we will have f̂_n-f^2_2≤ 2 ĝ_n-g^2_2+2 (b̂_l⃗^2-b_l⃗)^2+2E (Â_n-A)^2 where E (Â_n-A)^2=o(n^-2α/2α+q). Since functions f and g only differ by a constant b_l⃗, their wavelet coefficients coincide, that is, θ_j,k⃗^i=∫_[0,1]^q fψ⃗_j,k⃗^i=∫_[0,1]^q]gψ_j,k⃗^i. To make our analysis straightforward, we expand ĝ_n-g^2_2 as follows:

    ĝ_n^2-g^2_2=∑_1⃗≤k⃗≤2⃗^⃗j⃗_⃗0⃗∑_i=1^2^q-1 [θ̂_j_0,k⃗-θ_j_0,k⃗]^2+∑_j=j_0^J-1∑_1⃗≤k⃗≤2⃗^⃗j⃗∑_i=1^2^q-1 (θ̂_j,k⃗^i-θ_j,k⃗^i)^2
       +∑_j=J^∞∑_1⃗≤k⃗≤2⃗^⃗j⃗∑_i=1^2^q-1(θ_j,k⃗^i )^2≡ S_1+S_2+S_3.

First, we note that the term S_1 is asymptotically small. Indeed, by definition, θ̂_j_0,k⃗=y_j_0,k⃗. Since  θ̂_j_0,k⃗-θ_j_0,k⃗=(y_j_0,k⃗-θ̆_j_0,k⃗)+(θ̆_j_0,k⃗-θ_j_0,k⃗), we have

    S_1≤ C· 2^j_0qn^-1^2+CT^-2d=o(n^-2α/2α+q).

The term S_3 is also small asymptotically. To show that this is true, we note first that 2^j(α+q(1/2-1/s))(∑_1⃗≤k⃗≤2⃗^⃗j⃗∑_i=1^2^q-1|θ_j,k^i|^s)^1/s≤ M for any function f∈ B^α_s,t(M). Then, using the inequality ||x||_p_2≤||x||_p_1≤ q^1/p_1-1/p_2||x||_p_2 for any 0<p_1≤ p_2≤∞ and x∈^q, we obtain that

    S_3≤ 2^-2Jmin(α,α+q(1/2-1/s))=o(n^-2α/2α+q)

due to assumptions on J and α. In the next step, we will use the Proposition (<ref>) to analyze the term S_2. Next, we find out that 

    S_2≤ 2∑_j=j_0^J-1∑_1⃗≤k⃗≤2⃗^⃗j⃗∑_i=1^2^q-1 (θ̂^i_j,k⃗-θ̆_j,k⃗^i)^2+2∑_j=j_0^J-1∑_i=1^2^q-1(θ̆_j,k⃗^i-θ^i_j,k⃗)^2
       ≤∑_j=1^J-1∑_u=1^2^qj/L∑_i=1^2^q-1min{8∑_k⃗∈ B^i_j,u[θ̆_j,k⃗^i]^2,8λ^*Ln^-1}+6∑_j=j_0^J-1∑_1⃗≤k⃗≤2⃗^⃗j⃗∑_i=1^2^q-1[_j,k⃗^i]^2
       +Cn^-1+2∑_j=j_0^J-1∑_1⃗≤k⃗≤2⃗^⃗j⃗∑_i=1^2^q-1[θ̆_j,k⃗^i-θ_j,k⃗^i]^2
       ≤∑_j=1^J-1∑_u=1^2^qj/L∑_i=1^2^q-1min{8∑_k⃗∈ B^i_j,u[θ̆_j,k⃗^i]^2,8λ^*Ln^-1}+Cn^-1+CT^-2d.

At this point, we consider two different cases. First, start with p≥ 2. Select J_1=⌊q/2α+qlog_2n⌋ which implies that 2^J_1≈ n^q/2α+q. Thus, using the result of Lemma (<ref>). we obtain

    S_2≤ 8λ^*∑_j=j_0^J_1-1∑_u=1^2^qj/L∑_i=1^2^q-1Ln^-1+8∑_j=J_1^J-1∑_1⃗≤k⃗≤2⃗^⃗j⃗∑_i=1^2^q-1[θ̆_j,k⃗^i]^2+Cn^-1+CT^-2d≤ Cn^-2α/2α+q.

Next, consider the case of p<2. First, note that 

    ∑_u=1^2^jq/L∑_i=1^2^q-1(∑_k⃗∈ B^i_j,u[θ_j,k⃗^i]^2)^s/2
       ≤∑_1⃗≤k⃗≤ 2^j⃗([θ_j,k⃗^i]^2)^s/2≤ M2^-jws.

Select J_2 such that 2^J_2≍ n^1/2α+q(log n)^2-s/s(2α+q)+2(1-q). Using Lemma 6 from <cit.>, one obtains 

    ∑_j=J_2^J-1∑_u=1^2^j/Lmin{8∑_k⃗∈ B^i_j,u[θ̆_j,k⃗^i]^2,8λ^*Ln^-1}
       ≤ Cn^-2α/2α+q(log n)^2-s/s(2α+q)+2(1-q).

On the other hand, we also have 

    ∑_j=j_0^J_2-1∑_u=1^2^qj/Lmin{8∑_k⃗∈ B^i_j,u[θ̆_j,k⃗^i]^2,8λ^*Ln^-1}
       ∑_j=j_0^J_2-1∑_u=1^2^qj/L8λ_*Ln^-1≤ Cn^-2α/2α+q(log n)^2-s/[s(2α+q)+2(1-q)].

Thus, we can now confirm that the L_2 risk in the case p<2 is bounded from above uniformly as

    f̂_n-f^2_2≤ Cn^-2α/2α+q(log n)^2-s/[s(2α+q)+2(1-q)]


In order to ensure that the risk of b̂_l⃗ is negligible, we need to have κ^-4=o(n^-2α/2α+q); note that κ=n^1/4 satisfies this assumption. Also, to make the approximation error f_J-f^2_2 negligible, we need to have T^-2d=O(n^-2α/2α+q). It is easy to see that this is guaranteed by the inequality 3d/2q>2α/2α+q. Note that the latter, rather ponderous, assumption, is needed due to approximation over the q-dimensional Besov spaces.



Proof of Theorem (<ref>)
As in the proof of the Theorem (<ref>), and without loss of generality, we can assume that the med(X_i) is identically equal to zero; if this is not the case, an additional term can be estimated at the rate of n^-1=o(log n/n)^α/2α+q. Next, note that for all f∈Λ^α(M), the absolute values of its wavelet coefficients are |θ_j,k⃗^i|=|⟨ f, ψ_j,k⃗^i⟩ |≤ C2^-j(q/2+α) for some constant C>0 that does not depend on f. Also, note that for any random variables X_i, i=1,…,n, (∑_i=1^nX_i)^2≤(∑_i=1^n( X_i^2)^1/2)^2. Then, we have 

    (f̂_n(u⃗_0)-f(u⃗_0))^2
       = [∑_1⃗≤k⃗≤2⃗^j_0(θ̂_j_0,k⃗-θ_j_0,k⃗)ϕ_j_0,k⃗^i(u⃗_0)
    +∑_j=j_0^∞∑_1⃗≤k⃗≤2⃗^j∑_i=1^2^q-1(θ̂_j,k⃗^i-θ_j,k⃗^i)ψ_j,k⃗^(i))(u⃗_0)
       -(b̂_l⃗-b_l⃗)]^2
       ≤[( (b̂_l⃗-b_l⃗)^2)^1/2+∑_1⃗k⃗≤2⃗^j_0( (θ̂_j_0,k⃗-θ_j_0,k⃗)^2ϕ_j_0,k⃗^2(u⃗_0))^1/2
       +∑_j=j_0^J-1∑_1⃗≤k⃗≤2⃗^j∑_i=1^2^q-1( (θ̂_j,k⃗^i-θ_j,k⃗^i)^2[ψ_j,k⃗^i]^2(u⃗_0))^1/2+∑_j=J^∞∑_1⃗≤k⃗≤2⃗^j∑_i=1^2^q-1|θ_j,k⃗^iψ_j,k⃗^i(u⃗_0)|]^2
       ≡ (Q_1+Q_2+Q_3+Q_4)^2.

First of all, we note that, due to Proposition (<ref>), we have Q_1 =( (b̂_l⃗-b_l⃗)^2)^1/2=o(n^-α/2α+q). Next, clearly we have Q_2=∑_1⃗≤k⃗≤2⃗^j_0( (θ̂_j_0,k⃗-θ_j_0,k⃗)^2|ϕ_j_0,k⃗(u⃗_0)|=O(n^-1).  Recall that for any sequence of translated and rescaled ith wavelet ψ_j,k⃗^i there are at most N that are nonvanishing at the point u⃗_0; here, N is the length of support of ψ^i. Mathematically, we have K(t_0,j)={k⃗:ψ_j,k⃗^(i)(u⃗_0) 0} such that | K_t_0,j|≤ N.  
Thus, we have 

    Q_4=∑_j=J^∞∑_1⃗≤k⃗≤2⃗^j∑_i=1^2^q-1|θ_j,k⃗^i||ψ_j,k⃗^i(u⃗_0)|≤∑_j=J^∞N 2^qψ_∞2^jq/2C2^-j(q/2+α)≤ C T^-α.

Finally, if we select a sufficiently small τ, and use the second inequality from the Proposition (<ref>), we have 

    Q_3≤∑_j=j_0^J-1∑_K(t_0,j)∑_i=1^2^q-12^jq/2ψ_∞( (θ̂_j,k⃗^i-θ_j,k⃗^i)^2)^1/2
       ≤ C∑_j=j_0^J-12^jq/2[min(2^-j(q+2α)+T^-2α∧ 12^-jq,Ln^-1)+n^-2+τ]^1/2
       ≤  C(log n/n)^α/2α+q.

The final statement of this theorem is then obtained combining all inequalities for Q_1, Q_2, Q_3, and Q_4.




























apalike 

