Multivariate and Semiparametric Kernel Regression

Multivariate and Semiparametric Kernel Regression

Wolfgang H

ARDLE

Marlene M

ULLER

Institut fur Statistik und Okonometrie, Wirtschaftswissenschaftliche Fakultat Humboldt-Universitat zu Berlin, Germany

March 11, 1997

The paper gives an introduction to theory and application of multivariate and semipara- metric kernel smoothing. Multivariate nonparametric density estimation is an often used pilot tool for examining the structure of data. Regression smoothing helps in investigating the association between covariates and responses. We concentrate on kernel smoothing using local polynomial tting which includes the Nadaraya{Watson estimator. Some the- ory on the asymptotic behavior and bandwidth selection is provided. In the applications of the kernel technique, we focus on the semiparametric paradigm. In more detail we describe the single index model (SIM) and the generalized partial linear model (GPLM).

To appear in: M.G. Schimek (Ed.), Smoothing and Regression. Approaches, Computation and Ap- plication, 1996.

The research for this paper was supported by Sonderforschungsbereich 373 at the Humboldt-University Berlin. The work of M. Muller was supported in part by CentER, Tilburg University (The Netherlands).

The paper is printed using funds made available by the Deutsche Forschungsgemeinschaft.

1

Contents

1 Multidimensional Smoothing with Kernels 3

1.1 Multivariate Kernel Density Estimation . . . 3

1.1.1 Bias, Variance and Asymptotics . . . 5

1.1.2 Bandwidth selection and Graphical Representation . . . 8

1.2 Multivariate Kernel Regression . . . 14

1.2.1 Bias, Variance and Asymptotics . . . 16

1.2.2 Bandwidth Selection and Practical Aspects . . . 18

2 Semiparametric Generalized Regression Models 22

2.1 Generalizing the link function: Single Index Models . . . 24

2.1.1 Average Derivative Estimation . . . 24

2.1.2 Including Discrete Explanatory Variables . . . 26

2.2 Generalizing the index: Generalized Partial Linear Models . . . 28

2.2.1 Semiparametric Maximum Likelihood . . . 28

2.2.2 Practical Application . . . 31

2

Nonparametric smoothing methods serve three essential needs in statistical data analysis.

First they provide a exible analysis tool, often based on interactive graphical data rep- resentation (Scott, 1992). Second they help in constructing a model from observations, for example by comparison with concurrent models (Muller, 1988). Third they provide pilot estimators in adaptation problems, see Newey and Stoker (1993). Here we present the multivariate kernel smoother, examine the asymptotic properties of both density and regression estimators, and review applications of this technique in semiparametric statis- tics.

1 Multidimensional Smoothing with Kernels

In this section we review kernel smoothing methods for density and regression function estimation. Many ideas, in particular for asymptotics, bandwidth choice and graphical representation, are similar for both purposes. We can however only introduce a small part on the available material. In particular, for the regression case we restrict the presentation on the random design case. For a more detailed presentation of the subject we refer to the monographs by Hardle (1990 1991), Scott (1992), Wand and Jones (1995) and Fan and Gijbels (1995).

1.1 Multivariate Kernel Density Estimation

The goal of multivariate nonparametric density estimation is to approximate the probabil- ity density function (pdf) f(t) =f(t¹:::tq) of the random variables T = (T¹:::Tq)^T. The multivariate kernel density estimator in the q-dimensional case is dened as

fbh(t) = 1n

n

X

i⁼¹

h¹:::h1 q^K

Ti^1;t¹

h¹ ::: T^iq;tq

hq

!

(1)

K denoting a multivariate kernel function^K: IR^{q !}IR. Note, that (1) assumes that the bandwidth h is a vector of bandwidthsh = (h¹:::hq)^T.

What form shall the multidimensional kernel function ^K(u) = ^K(u¹:::uq) take on?

The easiest solution is to use a multiplicative kernel

K(u) =K(u¹):::K(uq) (2) withKdenoting an univariate kernel function. For univariate kernels with support ^;11]

(as the Epanechnikov kernelK(u) = 0:75(1^;u²)I(^ju^j1)) observations in a cube around t are used to estimate the density at the point t. An alternative is to use a genuine multivariate kernel function ^K(u), as e.g. the multivariate Epanechnikov

K(u)^/(1^;u^Tu) I(u^Tu1): 3

This type of multivariate kernels can be obtained from univariate by dening

K(u)^/K(^ku^k) (3) where ^ku^k =^pu^Tu denotes the Euclidean norm of the vector u. Note that we use ^/ to indicate that the appropriate constant has to be multiplied. Kernels of the form (3) use observations from a ball aroundt to estimate the pdf at t. This type of kernels is usually called spherical or radialsymmetric since ^K(u) has the same value for all u on a sphere around zero. Figure 1 shows the contour lines from a bivariate product and a bivariate radialsymmetric kernel on the left and right hand side, respectively.

-0.9 -0.6 -0.3 0.0 0.3 0.6 0.9

-0.9-0.6-0.30.00.30.60.9

-0.9 -0.6 -0.3 0.0 0.3 0.6 0.9

-0.9-0.6-0.30.00.30.60.9

Figure 1: Contours from bivariate product (left) and bivariate ra- dialsymmetric (right) Epanechnikov kernel.

Note that the kernel weights in Figure 1 correspond to equal bandwidth in each direc- tion, i.e. h = (h¹h²)^T = (11)^T. When we use dierent bandwidths, the observations aroundt in the density estimatef^bh(x) will be used with dierent weights in both dimen- sions.

Another approach is to use a nonsingular, symmetric bandwidth matrix

H

. The general form for the multivariate density estimator is then

fb^H(t) = 1n

n

X

i⁼¹

det(1

H

^)K

n

H

^{;1(Ti ;t)o= 1}_n^Xn

i⁼¹

K

H(Ti^;t) (4) see Silverman (1986) and Scott (1992). Here we introduce the short notation

K

H() = 1det(

H

)^K(

H

^;1)

analogously toKh in the one{dimensional case. A bandwidth matrix includes all simpler cases as special cases. An equal bandwidth h in all dimensions as in (1) corresponds to

4

H

⁼h

I

q where

I

q denotes the qq identity matrix. Dierent bandwidths as in (1) are equivalent to

H

^{= diag(h1:::h}q), the diagonal matrix with elements h¹:::hq.

What eect has the inclusion of o{diagonal elements? We will see that a good rule of thumb is to use a bandwidth matrix proportional to

^{b;1=2 where}

^b is the covariance matrix of the data. Hence, using such a bandwidth corresponds to a transformation of the data, so that they have an identity covariance matrix. As a consequence we can use bandwidth matrices to correct for correlation between the components of T. We have plotted the contour curves of product and radialsymmetric Epanechnikov weights with bandwidth matrix

H

⁼

0

@ 1 0:5 0:5 1

1

A 1=²

i.e. ^KH(u) =^K(

H

^;1u)=det(

H

), in Figure 2.

-0.9 -0.6 -0.3 0.0 0.3 0.6 0.9

-0.9-0.6-0.30.00.30.60.9

-0.9 -0.6 -0.3 0.0 0.3 0.6 0.9

-0.9-0.6-0.30.00.30.60.9

Figure 2: Contours from bivariate product (left) and bivariate ra- dialsymmetric (right) Epanechnikov kernel. Bandwidth matrix.

In the following we will consider statistical properties as bias, variance, the issue of bandwidth selection and applications for this estimator. We formulate all results for estimators with bandwidth matrices and multivariate kernel function^K.

1.1.1 Bias, Variance and Asymptotics

A consequence of the standard assumption on the non{negative kernel ^K

Z

K(u)du= 1 (5)

is that the estimatef^bH is a density function, i.e. ^R f^bH(t)dt= 1. The estimate is consistent in any point t of continuity of f:

fb^H(t) = 1n

n

X

i^=1KH(Ti^;t) = f(t) +op(1) (6) 5

if n ^!1,

H

^{!0 and} n

H

^!1, see e.g. Ruppert and Wand (1994). The derivation of the mean squared error MSE and the mean integrated squared error MISE is analogous to the one{dimensional case. We will sketch the asymptotic expansions and concentrate on the asymptotic mean integrated squared error AMISE.

As usual, AMISE has a bias part AIB and a variance part AIV . The bias of f^bH(t) is Ef^bH(t)^;f(t) and the integrated squared bias is

IB(

H

) =^{Z f}Ef^bH(t)^;f(t)^g2dt:

The asymptotic integrated squared bias AIB(

H

) is the rst order term of IB(

H

^{), i.e.}

IB(

H

)^;AIB(

H

)

AIB(

H

^{) =o(1)}

as

H

^!0, n^!1 and n

H

^!1. Dene now the integrated variance IV(

H

^{) =Z} E^ff^bH(t)^;Ef^bH(t)^g2dt

and the asymptotic integrated variance IV accordingly to IB. Then the asymptotic mean integrated squared error AMISE can be calculated as

AMISE(

H

) = AIB(

H

) + AIV (

H

): (7) A detailed derivation of the components of AMISE can be found in Scott (1992) or Wand and Jones (1995) and the references therein. As in the univariate case we use a second order Taylor expansion. Here and in the following we denote with ^rf the gradient and with ^Hf the Hessian matrix of second order partial derivatives of a function (here f).

Then the Taylor expansion of f() around t is

f(t+u) = f(t) +u^Trf(t_{) + 12}u^THf(t)t+o(u^Tu) see Wand and Jones (1995, p. 94). This leads to the expression

Ef^bH(t) = ^{Z KH}(u^;t)f(u)du=^{Z K}(s)f(t+

H

s)ds

Z

K(s) f(t) +s^T

H

^Trf(t_{) + 12}s^T

H

^THf(t)

H

s ds: (8) If we assume additionally to (5)

Z u^K(u)du= 0q (9)

Z uu^TK(u)du=²(^K)

I

q (10) then (8) yields Ef^bH(t)^;f(t) ¹²²(^K)tr^f

H

^THf(t)

H

^ghence

AIB(

H

_{) = 14}^{22(K)Z htrf}

H

^THf(t)

H

^{gi2 dt: (11)}

6

As in univariate density estimation, the leading term of the variance part is the second moment of the estimate, i.e.

Varⁿf^bH(t)^o = 1n

Z

fK

H(u^;t)^g2 du^; n¹

nEf^bH(t)^o2

Z 1

ndet(

H

^)K2(s)f(t+

H

s)ds

Z 1

ndet(

H

)^K2(s) ⁿf(t) +s^T

H

^Trf(t)^o ds

ndet(1

H

^{)kK k22f(t) (12)}

with ^{kK k2} denoting the q-dimensionalL²-norm of ^K. Hence

AIV(

H

^{) =} _ndet(¹

H

^{)kK k22 (13)}

and in summary we get the following AMISE formula for the multivariate kernel density estimator

AMISE(

H

_{) = 14}²²(^K)^{Z h}tr^f

H

^THf(t)

H

^gi2 dt+ 1

ndet(

H

^{)kK k22: (14)}

Let us now turn to the problem how to choose the AMISE optimal bandwidth. Again this is the bandwidth which balances bias{variance tradeo in AMISE. Denoteha scalar, such that

H

⁼h

H

^{0 and det(}

H

⁰) = 1. Then AMISE can be written as

AMISE(

H

_{) = 14}^{h422(K)Z htrf}

H

^T0Hf(t)

H

^{0gi2 dt+ 1}_nhq^{kK k22:}

If we only allow changes in h the optimal orders for the smoothing parameter h and AMISE are

h⁰ =O(n^;1=(4+q)) AMISE(h⁰

H

⁰) =O(n^;4=(4+q)):

Hence, this density estimator has a rather slow rate of convergence, especially ifq is large.

If we consider

H

=h

I

q (the same bandwidth in allq dimensions) and we x the sample size n, then the AMISE optimal bandwidth has to be considerably larger than in the one{dimensional case to make sure that the estimate has reasonably small variability.

Some ideas of comparable sample sizes to reach the same quality of the density estimates over dierent dimensions can be found in Silverman (1986, p. 94) and Scott and Wand (1991). Moreover, the computational eort of this technique increases with the number of dimensionsq. Therefore, multidimensional density estimation is usually not practically applied if q5.

7

1.1.2 Bandwidth selection and Graphical Representation

The problem of an automatic, data-driven choice of the bandwidth

H

is of great impor- tance in the multivariate case. In one or two dimensions we may choose an "appropriate"

bandwidth interactively by looking at the sequence of density estimates for dierent band- widths. But how can this be done in three, four or more dimensions? The problem of graphical representation arises, which we address next.

Theoretically the bandwidth selection problem can be handled as in the one{dimensional case. Typically, one searches for a global bandwidth

H

or a local bandwidth

H

^{(t). Two}

approaches are frequently used in both cases

plug{in bandwidths, in particular "rule{of{thumb" bandwidths,

resampling methods, in particular cross{validation and bootstrap.

We will introduce generalizations for Silverman's rule{of{thumb and least squares cross{

validation to stress the analogy with the one{dimensional bandwidth selectors.

Rule{of{thumb Bandwidth

Rule{of{thumb bandwidth selection provides a formula arising from a reference distribution. Obviously, the pdf of a multivariate normal distri- bution Nq(

) is a good candidate for a reference distribution in the multivariate case.

Suppose that the kernel ^K is Gaussian, i.e. the pdf of Nq(0q

I

^q). Note that ²(^K) = 1 and ^{kK k22} = 2^;q;q=2 in this case. Hence, from (14) and the fact that

Z tr^f

H

^THf(t)

H

^g]²dt= 1

2^q+2q=2det(

)¹⁼²

h2tr(

H

^T

^;1

H

)²+^ftr(

H

^T

^;1

H

)^g2i we can easily derive rule{of{thumb formulae for dierent assumptions on

H

^and

^.

In the simplest case, i.e. that we consider

H

^and

to be diagonal matrices

H

⁼

diag(h¹:::hq) and

^{= diag(12:::}q²), this leads to

ehj =

4 q+ 2

!

1=⁽q⁺⁴⁾

n^;1=(q+4)j: (15) Note that this formula coincides with Silverman's rule{of{thumb in the case q = 1, see Silverman (1986, p. 45). Replacing the j's by estimates and noting the rst factor is always between 0.924 and 1.059, we arrive at Scott's rule

hbj =n^;1=(q+4)bj (16) see Scott (1992, p. 152).

8

It is dicult to derive the rule{of{thumb for general

H

^and

. However, (15) shows that it might be a good idea to choose the bandwidth matrix

H

proportional to

^{1=2. In}

this case we get as generalization of Scott's rule

c

H

^=n;1=(q+4)

^{b1=2: (17)}

We remark that this rule is equivalent to apply a Mahalanobis transformation on the data (to transform the estimated covariance matrix to identity), then to compute the kernel estimate with equal bandwidthsh=n^1=(q+4) and nally to retransform the estimated pdf back to the original scale.

But before we go on with applications, let us consider what we can do, if we want to use a kernel dierent from the Gaussian. The idea of canonical kernels by Marron and Nolan (1988) can be easily extended to the multivariate case. Consider a kernel ^K and all equivalent kernel functions ^K = ^;1K(=) with 0. Although that is a scalar, it is working on q{variates arguments of ^K. Now we have ^kKk22 =^{;qkK k22} and ²(^K) =²²(^K). As in the one{dimensional case we choosesuch that the bias-variance tradeo in AMISE(

H

^K) is independent of ^K. This yields

²²(^K0) =^{kK0k22 () 0} =

(

kK k 2

²²(^K2)

)

1=⁽q⁺⁴⁾

:

⁰ again is called canonical bandwidth of the kernel ^K. Denote now ^KA a kernel function with canonical bandwidth ^0A and ^KB a kernel function with canonical bandwidth ^0B. Suppose we have used

H

^A with kernel ^KA and we want to recompute the kernel density estimate with kernel ^KB. Then it holds

AMISE(

H

A^KA) AMISE(

H

B^KB)

if

H

B = ^0B

^0A

H

A (18)

which allows to adjust bandwidths for dierent kernel as in the one{dimensional case.

Let us consider an example. Suppose we want to use the product Quartic kernel ^KQ instead of the q-dimensional Gaussian ^KG which is faster in direct computation because of its compact support on ^;11]. Which is the equivalent rule{of{thumb to (17) in this case? Here we have ^0G = ^f1=(2^p)^gq=(q+4) and ^0Q = (495^q=7^q)^1=(q+4) which gives the canonical bandwidths in Table 1 for dimensions q= 1:::5.

The fourth column of Table 1 gives the factor which the rule{of{thumb bandwidth matrix in (17) needs to be multiplied with to obtain the rule{of{thumb bandwidth for the multiplicative Quartic kernel. Of course all rule{of{thumb bandwidths for other kernel functions can be calculated in a similar way.

9

q ^{0G 0Q 0Q}=^0G 1 0.7764 2.0362 2.6226 2 0.6558 1.7100 2.6073 3 0.5814 1.5095 2.5964 4 0.5311 1.3747 2.5883 5 0.4951 1.2783 2.5820

Table 1: Bandwidth adjusting factors for Gaussian and multiplica- tive Quartic Kernel for dierent dimensionsq.

For a product kernel ^K holds ²(^K) = ²(K) and ^{kK k2} =^kK^kq2 when K denotes the corresponding univariate kernel. A table of values ²(K), ^kK^k22 can be found in Hardle (1991, p.239) for example.

Principally, all plug{in methods for the one{dimensional kernel density estimation can be extended to the multivariate case. See Wand and Jones (1994) for details on multivariate plug{in bandwidth selection.

Cross{validation

As we mentioned before, the cross{validation method is fairly inde- pendent of the special structure of the parameter or function estimate. Considering the bandwidth choice problem, cross{validation techniques allow to adapt to a wider class of density functionsfthan the rule{of{thumb approach. (Remember that the rule{of{thumb bandwidth is optimal for the reference pdf, hence it may fail for multimodal densities for instance.)

Recall, that in contrast to the rule{of{thumb approach, least squares cross{validation for density estimation aims to estimate the ISE optimal bandwidth. Here we approximate the integrated squared error

ISE(

H

^{) = Z ffbH(t);f(t)g2dt}

= ^Z f^bH2(t)dt^;2^Z f^bH(t)f(t)dt+^Z f²(t)dt: (19) Apparently, this is the same formula as in the the one{dimensional case and with the same arguments the last term of (19) can be ignored. The rst term again can be easily calculated from the data. Hence, only the second term of (19) is unknown and has to be estimated. However, observe that ^R f^bH(t)f(t)dt = Ef^bH(T), where the only new aspect now is thatT is q{dimensional. As in the one{dimensional case we estimate this term by a leave{one{out estimator

Efd^bH(T) = 1n

n

X

i⁼¹

fb^H;i(Ti) 10

where

fb^H;i(t) = 1n^;1

n

X

i⁶⁼jj⁼¹

K

H(Tj ^;t):

This yields the multivariate cross{validation criterion as a straightforward generalization of CV in the one{dimensional case:

CV(

H

) = 1

n²det(

H

^{) n}

X

i⁼¹ n

X

j^=1K?^Kn

H

^;1(Tj ^;ti)^o; 2 n(n^;1)

n

X

i⁼¹ n

X

j⁼¹ j⁶⁼i

K

H(Tj ^;Ti): The diculty comes in by the fact that the bandwidth is now a qq matrix

H

. In the most general case, this means, we have to minimize overq(q+1)=2 parameters. Still, if we assume

H

to be a diagonal matrix, this remains a q{dimensional optimization problem.

This holds as well for other cross{validation approaches. Multivariate resampling methods for bandwidth selection are discussed in more detail in Sain, Baggerly and Scott (1994).

Graphical Representation

Consider now the problem to graphically display a multi- variate density estimate. Assume rst q = 2. Here we are still able to show the density estimate in a 3-dimensional plot. This is in particular useful if the estimated function can be rotated on the computer screen interactively. For a two-dimensional presentation a contour plot gives often more insight to the structure of the data.

In the following, we will use the credit data from Fahrmeir and Hamerle (1984), Fahrmeir and Tutz (1994) for illustration. This data set consists of n = 1000 clients, 700 paid a credit back without problems, 300 did not. Among a number of categori- cal variables (running account, previous credits, purpose, personal attributes etc.) three continuous variables are available: duration and amount of credit as well as age.

Figures 3, 4 (upper panels) display a two-dimensional density estimate fbh(t) =f^bh(t¹t²)

for log(duration, log(amount and log(amount, log(age, respectively. We use the subscript h to indicate that we used a diagonal bandwidth matrix

H

^{= diag(}h¹h²).

Additionally, Figures 3, 4 (lower panels) gives contour plots of these density estimates.

It is easily observed, that both distributions are rather symmetric. This is due to the logarithmic transformation. In the duration direction a typical bimodal structure can be recognized. This slightly reproduces in the amount direction. Obviously, both variables are related with positive correlation.

Here, the bandwidth was chosen accordingly to the general \rule{of{thumb" (17), which tends to oversmooth multimodal structures of the data. In fact, the durations of credits are multiples of 6 months in most case. The two clear modes that we observe are those for durations 12 and 24 months. In all applications of this paper we use the

11

X

1.5 2.0 2.5 3.0 3.5 4.0

Y

6.0 7.0 8.0 9.0

Z

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

X: duration Y: amount

Z: density estimate (*10 -1 )

Density: duration & amount

1.5 2.0 2.5 3.0 3.5 4.0

duration

6.06.57.07.58.08.59.09.5

amount

Contours: duration & amount

Figure 3: Two{dimensional density estimate (upper panel) and den- sity contours (lower panel) for duration and amount. h¹ = 0:48, h² = 0:64. Credit data, Fahrmeir and Hamerle (1984).

12

X

6.0 7.0 8.0 9.0

Y

3.0 3.2 3.4 3.83.6 4.0 4.2

Z

0.0 1.0 2.0 3.0 4.0 5.0 6.0

X: amount Y: age

Z: density estimate (*10 -1 )

Density: amount & age

6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5

amount

3.03.23.43.63.84.04.2

age

Contours: amount & age

Figure 4: Two{dimensional density estimate (upper panel) and density contours (lower panel) for amount and age. h¹ = 0:64, h² = 0:25. Credit data, Fahrmeir and Hamerle (1984).

13

X

1.5 2.0

2.5 3.0

3.5 Y 4.0

6.0 7.0 8.0 9.0

Z

3.0 3.2 3.4 3.6 3.8 4.0 4.2

X: duration Y: amount Z: age

0.0 0.5 1.0 1.5 2.0 2.5 3.0

0.02.04.06.08.010.0

Level

Figure 5: Three{dimensional density contours for duration, amount and age. h¹ = 0:56, h² = 0:75, h³ = 0:29. Credit data, Fahrmeir and Hamerle (1984).

Quartic (Biweight) product kernel. Recall that the the univariate Quartic kernel isK(u) = 0:9375(1^;u²)²I(^ju^j1).

For three{dimensional density estimates, it is always possible to hold one variable xed and to plot the density function only in dependence of the other variables. Alternatively, we can again plot contours of the density estimate, which now mean three-dimensional surfaces. Figure 5 shows this for the credit scoring variables. In the original version of this plot, red, green and blue surfaces show the values of the density estimate at the levels (in percent) indicated on the right. Colors and the possibility to rotate the contours on the computer screen eases the exploration of the data structures a lot. Of course, we are restricted to two-dimensional plots here. However, one can clearly recognize the ellipsoidal structure of the contour which indicates a relatively symmetric distribution.

1.2 Multivariate Kernel Regression

Multivariate nonparametric regression aims to estimate the functional relation between a response variable Y and a multivariate explanatory variable T, i.e. the conditional expectation

E(Y^jT) =E(Y^jT¹:::Tq) =m(T) (20) 14

where as before T = (T¹:::Td)^T. The relation E(Y^jT) =^Z yf(y^jt)dy=

R yf(yt)dy f(x)

leads by replacing the multivariate densities f(yt) by the kernel density estimate fbh^H(yt) = 1n

n

X

i⁼¹Kh(Yi^;y)^KH(ti^;t)

and f(t) = fT(t) by (4) to the multivariate generalization of the Nadaraya{Watson esti- mator:

mc^H(t) =

n

P

i^=1KH(Ti^;t)Yi n

P

i^=1KH(Ti^;t) : (21)

Hence, the multivariate kernel regression estimator is just a weighted sum of the observed responses Yi. The denominator ensures that the weights sum up to 1. Depending on the choice of the kernel, ^cm^H(t) is a weighted average of thoseYi whereTi lies in a ball or cube around t.

Note that the multivariate Nadaraya{Watson estimator is a local constant estimator, i.e. the solution of

min

0

n

X

i⁼¹

fYi^;0g2KH(Ti^;t):

Replacing⁰ by a polynomial inTi^;tyields a local polynomial kernel regression estimator.

This denition of local polynomial kernel regression is a straightforward generalization of the univariate case. For details see Ruppert and Wand (1994). Let us illustrate this with the example of a local linear regression estimate. The minimization problem is here

min⁰¹ n

X

i⁼¹

nYi^;0;(Ti^;t)^T1o2KH(Ti^;t): The solution of the problem can hence equivalently be written as

b= (^b0b1T)^T =

T

^T

WT

^;1

T

^T

WY

⁽²²⁾

using the notations

T

=

0

B

B

B

@

1 (T^1;t)^T ... ...

1 (Tn^;t)^T

1

C

C

C

A

Y

=

0

B

B

B

@

Y¹ Y...n

1

C

C

C

A

and

W

^{= diag(KH(T1;t):::KH(Tn;t})). In (22)^b0 estimates the regression function itself, whereas ^b1 estimates the partial derivatives w.r.t. the components T. In the following we denote the multivariate local linear estimator as

cm^1H(t) =^b0(t): (23) 15

1.2.1 Bias, Variance and Asymptotics

The asymptotic conditional variance of the Nadaraya{Watson estimatorm^cH and the local linear ^cm^1H is identical and its derivation can be found in detail in Ruppert and Wand (1994):

Var^fcm^H(t)^jT¹:::Tn^g= 1

ndet(

H

)^{kK k22 2}(t)

f(t) ^f1 +op(1)^g (24) with ²(t) denoting the variance function in Var(Y^jt).

We sketch the derivation of the asymptotic conditional bias since we nd remarkable dierences between both estimators. Denote M the second order Taylor expansion of (m(T¹):::m(Tn))^T, i.e.

M m(t)11n+L(t_{) + 12}Q(t) =

T

0

@ m(t)

rm(t)

1

A+ 12Q(t) (25) with

L(t) =

0

B

B

B

@

(T^1;t)^Trm(t) (Tn^;t...)^Trm(t)

1

C

C

C

A Q(t) =

0

B

B

B

@

(T^1;t)^THm(t)(T^1;t) (Tn^;t)^TH...m(t)(Tn^;t)

1

C

C

C

A: Additionally to (6) it holds

n1

n

X

i^=1KH(Ti ^;t)(Ti^;t) = ²(^K)

HH

^Trf(t) +op(

HH

^T11d) 1n

n

X

i⁼¹

K

H(Ti^;t)(Ti^;t)(Ti^;t)^T = ²(^K)f(t)

HH

^{Trf(t) +op(}

HH

^T)

see Ruppert and Wand (1994). Therefore the denominator of the conditional asymp- totic expectation of the Nadaraya{Watson estimator ^cm^H is approximately f(t). Using E(

Y

^{jT1:::Tn) =M} and the Taylor expansion for M we have

E^fcm^HjT¹:::Tn^g

ff(t) +op(1)^g;11 n

n

X

i⁼¹

K

H(Ti^;t)m(t) +^Xn

i⁼¹

K

H(Ti^;t)(Ti^;t)^Trm(t) +^Xn

i^=1KH(Ti^;t)(Ti^;t)^THm(t)(Ti^;t)

ff(t)^g;1f(t)m(t) +²(^K)^rm

HH

^Trf + 12²(^K)f(t)tr^f

H

^THm(t)

H

^g: This is summarized in the following theorem.

THEOREM 1

The conditional asymptotic bias and variance of the multivariate Nadaraya{Watson kernel 16

regression estimator are

E^fcm^HjT¹:::Tn^g;m(t) ²(^K)^rm(t)^T

HH

^Trf(t)

f(t) + 12²(^K)tr^f

H

^THm(t)

H

^g

Var^fcm^HjT¹:::Tn^g 1

ndet(

H

)^{kK k22 2}(t) f(t) in the interior of the support of fT.

Recall the notation e¹ = (10:::0)^T for the rst unit vector in IR^d. Then we can write the local linear estimator as

cm^1H(t) =e^T1

T

^T

WT

^;1

T

^T

WY

^:

Now we have using (22) and (25)

E^fcm^1HjT¹:::Tn^g;m(t)

= e^T1

T

^T

WT

^;1

T

^T

WT

8

<

: 0

@ m(t)

rm(t)

1

A+ 12Q(t)

9

=

;m(t)

= 12e^T1

T

^T

WT

^;1

T

^T

W

Q(t)

since e^T1m(t)^rm(t)^T] = m(t). Hence, the numerator of the asymptotic conditional bias only depends on the quadratic term. This is one of the key points in asymptotics for local polynomial estimators. If we would use local polynomials of order d and expand M up to order d+ 1, then only the term of order d+ 1 would appear in the numerator of the asymptotic conditional bias. Of course this to be paid with a more complicated structure of the denominator.

THEOREM 2

The conditional asymptotic bias and variance of the multivariate local linear regression estimator are

E^fcm^1HjT¹:::Tn^g;m(t) 1

2²(^K)tr^f

H

^THm(t)

H

^g

Var^fcm^1HjT¹:::Tn^g 1

ndet(

H

)^{kK k22 2}(t) f(t) in the interior of the support of fT.

For all omitted details on the proof of Theorem 2 we refer again to Ruppert and Wand (1994). They also point out that the local linear estimate has same order conditional bias in the interior as well as in the boundary of the support of fT. Fan, Gasser, Gij- bels, Brockmann and Engel (1993) point out that the multivariate local linear t with Epanechnikov kernel is a best linear estimator and has a minimax eciency of at least 89.4% among all estimators.

17

1.2.2 Bandwidth Selection and Practical Aspects

Principally, the methods to choose a smoothing parameter in nonparametric regression are the same as in density estimation. Again, plug{in and resampling ideas are employed for nding a global bandwidth

H

or a local bandwidth

H

(t).

For our presentation, we concentrate on the classical cross{validation bandwidth se- lector. As a motivation, we introduce the residual sum of squares (RSS) as a (naive) way to asses the goodness of t

RSS(

H

) =n^;1Xn

i⁼¹

fYi^;cm^H(Xi)^g2 (26) which is also called resubstitution estimate for the averaged squared error (ASE). Note, that we concentrate on the Nadaraya{Watson estimator in the moment.

There is a problem with the RSS: Yi is used in ^cm^H(Xi) to predict itself. As a conse- quence, ASE(

H

) can be made arbitrarily small by letting

H

^! 0 (in which case m^cH is an interpolation of the Yi's). This leads to the cross{validation function

CV(

H

^{) =n;1Xn}

i⁼¹

fYi^;cm^H;i(Xi)^g2: (27) This function replaces m^cH(Xi) in (26) with the leave{one{out-estimator

cm^H;i(Xi) =

Pj⁶⁼i^KH(Xi^;Xj)Yj

Pj⁶⁼i^KH(Xi^;Xj) : (28) and is equivalent to a dierent approach, which multiplies each term in RSS(

H

) by a penalizing functionthat is correcting for the downward bias of the resubstitution estimate.

For the Nadaraya{Watson estimator CV(

H

) = 1n

n

X

i^=1fYi^;cm^H;i(Xi)^g2

= 1n

n

X

i⁼¹

fYi^;cm^H(Xi)^g2

(Yi^;cm^H;i(Xi) Yi^;cm^H(Xi)

)

2 (29)

and

Yi^;cm^H(Xi) Yi^;cm^H;i(Xi) =

P

j ^KH(Xi^;Xj)Yj^;Yi^Pj ^KH(Xi^;Xj)

P

j⁶⁼i^KH(Xi^;Xj)Yj^;Yij^P6=i^KH(Xi^;Xj)

P

j⁶⁼i^KH(Xi^;Xj)

P

j ^KH(Xi^;Xj)

= 1^{; KH}(0)

P

j ^KH(Xi^;Xj): (30)

Therefore the cross{validation approach is equivalent to the penalizing functions concept and shares the same asymptotic properties. Note that (30) is a function of the i{th

18

diagonal element of the smoother matrix. More precisely, cross{validation is equivalent with generalized cross{validation (Craven and Wahba, 1979) in this case. Hardle, Hall and Marron (1988) show asymptotic optimality of the selected bandwidth, the rate of convergence is slow though. An improved bandwidth selection is discussed in Hardle, Hall and Marron (1992).

We want to remark that (29) and (30) also imply that the computation ofCV(

H

^{) needs}

actually not more computational eort than the computation of m^H(X¹):::m^H(Xn).

However, the optimization over a matrix

H

might be cumbersome, hence diagonal band- width matrices (or even

H

^{= h}

I

q with appropriate standardization of the data) are still preferred in practice.

Before we consider cross{validation bandwidth selection in the local linear case, we want to comment on the practical computation of the estimator. Principally, since multi- variate kernel regression estimators can be expressed as local polynomial estimators, their computation can be done by any statistical package that is able to run weighted least squares regression. However, since we estimate a function, this weighted least squares regression has to be performed in all observation points or on a grid of points in IR^q. Therefore, explicit formulae are useful.

We will give an formula for the multivariate local linear estimator in the following.

Consider for a xed point t the sums

S

0 =^S0(t) = ^Xn

i⁼¹

K

H(Ti^;t)

S

1 =^S1(t) = ^Xn

i^=1KH(Ti^;t)(Ti^;t)

S

2 =^S2(t) = ^Xn

i⁼¹

K

H(Ti^;t)(Ti^;t)(Ti^;t)^T

T

0 =^T0(t) = ^Xn

i^=1KH(Ti^;t)Yi

T

1 =^T1(t) = ^Xn

i^=1KH(Ti^;t)(Ti^;t)Yi:

Note that ^S1 and ^T1 are q{variate vectors and that ^S2 is a qq matrix. Then for the local linear estimate we can write

b=

0

@ S

0 ST

1

S

1 S

2 1

A

;1 0

@ T

0

T

1 1

A: (31)

For the regression function we need only the rst component e^T1b. Applying block{wise matrix inversion we obtain

e^T1

0

@ S

0 ST

1

S

1 S

2 1

A

;1

=^{S0;S1TS2;1S1;1}

1 ^;S1TS2;1

!

19

and hence

cm^1H(t) = ^{T0 ;S1TS2;1T1}

S

0

;ST

1 S

;1

2 S

1

: (32)

The cross{validation criterion here is a weighted RSS as in (29). If we denote the leave{one{out estimator ^cm^1H;i(t) and dene its components accordingly, we observe

S

0^;i = ^S0;KH(0) ^S1;i = ^{S1 S2;}i = ^S2

T

0^;i = ^T0;Yi^KH(0) ^T1;i = ^T1: This means

mc^1H;i(t) = ^T0;Yi^KH(0)^;S1TS2;1T1

S

0

;K

H(0)^;S1TS2;1S1 which yields in analogy to (30)

Yi^;m^cH(Xi)

Yi^;cm^H;i(Xi) = 1^;

K

H(0)

S

0

;ST

1 S

;1

2 S

1

: (33)

As in the Nadaraya{Watson case, (33) is a function of the i{th diagonal element of the smoother matrix. A summary of bandwidth selection methods other than cross{validation can be found in particular in Fan and Gijbels (1995). They also cover rule{of{thumb approaches.

X

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Y

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Z

0.0 1.0 2.0 3.0 4.0

X: x1 Y: x2 Z: m

True Function

X

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Y

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Z

0.0 1.0 2.0 3.0

X: x1 Y: x2 Z: mhat

Nadaraya-Watson

Figure 6: Two{dimensional Nadaraya{Watson Estimate.

20