Nonparametric econometrics: The np package

We describe the R np package via a series of applications that may be of interest to applied econometricians. The np package implements a variety of nonparametric and semiparametric kernel-based estimators that are popular among econometricians. There are also procedures for nonparametric tests of significance and consistent model specification tests for parametric mean regression models and parametric quantile regression models, among others. The np package focuses on kernel methods appropriate for the mix of continuous, discrete, and categorical data often found in applied settings. Data-driven methods of bandwidth selection are emphasized throughout, though we caution the user that data-driven bandwidth selection methods can be computationally demanding.


Introduction
Devotees of R (R Development Core Team 2008) are likely to be aware of a number of nonparametric kernel1 smoothing methods that exist in R base (e.g., density) and in certain R packages (e.g., locpoly in the KernSmooth package Wand and Ripley 2008).These routines deliver nonparametric smoothing methods to a wide audience, allowing R users to nonparametrically model a density or to conduct nonparametric local polynomial regression, by way of example.
The appeal of nonparametric methods, for applied researchers at least, lies in their ability to reveal structure in data that might be missed by classical parametric methods.Nonparametric kernel smoothing methods are often, however, much more computationally demanding than their parametric counterparts.
In applied settings we often encounter a combination of categorical and continuous datatypes.Those familiar with traditional nonparametric kernel smoothing methods will appreciate that these methods presume that the underlying data is continuous in nature, which is frequently not the case.One approach towards handling the presence of both continuous and categorical data is called a 'frequency' approach, whereby data is broken up into subsets ('cells') corresponding to the values assumed by the categorical variables, and only then do you apply say density or locpoly to the continuous data remaining in each cell.Nonparametric frequency approaches are widely acknowledged to be unsatisfactory as they often lead to substantial efficiency losses arising from the use of sample splitting, particularly when the number of cells is large.
Recent theoretical developments offer practitioners a variety of kernel-based methods for categorical data only (i.e., unordered and ordered factors), or for a mix of continuous and categorical data.These methods have the potential to recapture the efficiency losses associated with nonparametric frequency approaches as they do not rely on sample splitting, rather, they smooth the categorical variables in an appropriate manner; see Li and Racine (2007a) and the references therein for an in-depth treatment of these methods, and see also the references listed in the bibliography.
The np package implements recently developed kernel methods that seamlessly handle the mix of continuous, unordered, and ordered factor datatypes often found in applied settings.The package also allows the user to create their own routines using high-level function calls rather than writing their own C or Fortran code. 2 The design philosophy underlying np aims to provide an intuitive, flexible, and extensible environment for applied kernel estimation.We appreciate that there exists tension among these goals, and have tried to balance these competing ends, with varying degrees of success.np is available from the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project.org/package=np.
In this article, we illustrate the use of the np package via a number of empirical applications.Each application is chosen to highlight a specific econometric method in an applied setting.We begin first with a discussion of some of the features and implementation details of the np package in Section 2. We then proceed to illustrate the functionality of the package via a series of applications, sometimes beginning with a classical parametric method that will likely be familiar to the reader, and then performing the same analysis in a semi-or nonparametric framework.It is hoped that such comparison helps the reader quickly gauge whether or not there is any value added by moving towards a nonparametric framework for the application they have at hand.We commence with the workhorse of applied data analysis (regression)  in Section 3, beginning with a simple univariate regression example and then moving on to a multivariate example.We then proceed to nonparametric methods for binary and multinomial outcome models in Section 4. Section 5 considers nonparametric methods for unconditional probability density function (PDF) and cumulative distribution function (CDF) estimation, while Section 6 considers conditional PDF and CDF estimation, and nonparametric estimators of quantile models are considered in Section 7. A range of semiparametric models are then considered, including partially linear models in Section 8, single-index models in Section 9, and finally varying coefficient models are considered in Section 10.

Important implementation details
In this section we describe some implementation details that may help users navigate the methods that reside in the np package.We shall presume that the user is familiar with the traditional kernel estimation of, say, density functions (e.g., Rosenblatt 1956;Parzen 1962) and regression functions (e.g., Nadaraya 1965;Watson 1964) when the underlying data is continuous in nature.However, we do not presume familiarity with mixed-data kernel methods hence briefly describe modifications to the kernel function that are necessary to handle the mix of categorical and continuous data often encountered in applied settings.These methods, of course, collapse to the familiar estimators when all variables are continuous.

The primacy of the bandwidth
Bandwidth selection is a key aspect of sound nonparametric and semiparametric kernel estimation.It is the direct counterpart of model selection for parametric approaches, and should therefore not be taken lightly.np is designed from the ground up to make bandwidth selection the focus of attention.To this end, one typically begins by creating a 'bandwidth object' which embodies all aspects of the method, including specific kernel functions, data names, datatypes, and the like.One then passes these bandwidth objects to other functions, and those functions can grab the specifics from the bandwidth object thereby removing potential inconsistencies and unnecessary repetition.For convenience these steps can be combined should the user so choose, i.e., if the first step (bandwidth selection) is not performed explicitly then the second step will automatically call the omitted first step bandwidth selection using defaults unless otherwise specified, and the bandwidth object could then be retrieved retroactively if so desired.Note that the combined approach would not be a wise choice for certain applications such as when bootstrapping (as it would involve unnecessary computation since the bandwidths would properly be those for the original sample and not the bootstrap resamples) or when conducting quantile regression (as it would involve unnecessary computation when different quantiles are computed from the same conditional cumulative distribution estimate).
Work flow therefore typically proceeds as follows: 1. compute data-driven bandwidths; 2. using the bandwidth object, proceed to estimate a model and extract fitted or predicted values, standard errors, etc.; 3. optionally, plot the object.
In order to streamline the creation of a set of complicated graphics objects, plot (which calls npplot) is dynamic; i.e., you can specify, say, bootstrapped error bounds and the appropriate routines will be called in real time.Be aware, however, that bootstrap methods can be computationally demanding hence some plots may not appear immediately in the graphics window.

Data-driven bandwidth selection methods
We caution the reader that data-driven bandwidth selection methods can be computationally demanding.We ought to also point out that data-driven (i.e., automatic) bandwidth selection procedures are not guaranteed always to produce good results due to perhaps the presence of outliers or the rounding/discretization of continuous data, among others.For this reason, we advise the reader to interrogate their bandwidth objects with the summary command which produces a table of bandwidths for the continuous variables along with a constant multiple of σ x n α , where σ x is a variable's standard deviation, n the number of observations, and α a known constant that depends on the method, kernel order, and number of continuous variables involved, e.g., α = −1/5 for univariate density estimation with one continuous variable and a second order kernel.Seasoned practitioners can immediately assess whether undersmoothing or oversmoothing may be present by examining these constants, as the appropriate constant (called the 'scaling factor') that is multiplied by σ x n α often ranges from between 0.5 to 1.5 for some though not all methods, and it is this constant that is computed and reported by summary.Also, the admissible range for the bandwidths for the categorical variables is provided when summary is used, which some readers may also find helpful.
We caution users to use multistarting for any serious application (multistarting refers to restarting numerical search methods from different initial values to avoid the presence of local minima -the default is the minimum of the number of variables or 5 and can be changed via the argument nmulti =), and do not recommend overriding default search tolerances (unless increasing nmulti = beyond its default value).
We direct the interested reader to the frequently asked questions document on the author's website (http://socserv.mcmaster.ca/racine/np_faq.pdf) for a range of potentially helpful tips and suggestions surrounding bandwidth selection and the np package.

Interacting with np functions
A few words about the R data.frameconstruct are in order.Data frames are fundamental objects in R, defined as "tightly coupled collections of variables which share many of the properties of matrices and of lists, used as the fundamental data structure by most of R's modeling software."A data frame is "a matrix-like structure whose columns may be of differing types (numeric, logical, factor and character and so on A few words on the formula interface are in order.We use the standard formula interface as it provides capabilities for handling missing observations and so forth.This interface, however, is simply a convenient device for telling a routine which variable is, say, the outcome and which are, say, the covariates.That is, just because one writes x 1 + x 2 in no way means or is meant to imply that the model will be linear and additive (why use fully nonparametric methods to estimate such models in the first place?).It simply means that there are, say, two covariates in the model, the first being x 1 and the second x 2 , we are passing them to a routine with the formula interface, and nothing more is presumed nor implied.This will likely be obvious to most R users, but we point it out simply to avoid any potential confusion for those unfamiliar with kernel smoothing methods.

Writing your own functions
We have tried to make np flexible enough to be of use to a wide range of users.All options can be tweaked by the user (kernel function, kernel order, bandwidth type, estimator type and so forth).One function, npksum, allows you to create your own estimators, tests, etc.The function npksum is simply a call to specialized C code, so you get the benefits of compiled code along with the power and flexibility of the R language.We hope that incorporating the npksum function renders the package suitable for teaching and research alike.

Generalized product kernels
As noted above, traditional nonparametric kernel methods presume that the underlying data is continuous in nature, which is frequently not the case.The basic idea underlying the treatment of kernel methods in the presence of a mix of categorical and continuous data lies in the use of what we call 'generalized product kernels', which we briefly summarize.
Suppose that you are interested in kernel estimation for 'unordered' categorical data, i.e., you are presented with discrete data X d ∈ S d , where S d denotes the support of X d .We use x d s and X d is to denote the sth component of x d and X d i (i = 1, . . ., n), respectively.Following Aitchison and Aitken (1976), for (1) Note that when ) becomes an indicator function, and if Observe that these weights sum to 1.
If, instead, some of the discrete variables are ordered (i.e., are 'ordinal categorical variables'), then we should use a kernel function which is capable of reflecting their ordered status.Assuming that x d s can take on c s different ordered values, {0, 1, . . ., c s − 1}, Aitchison and Aitken (1976, p. 29) suggested using the kernel function given by where Observe that these weights sum to 1.
If instead a variable was continuous, you could use, say, the second order Gaussian kernel, namely Observe that these weights integrate to 1.
The 'generalized product kernel function' for a vector of, say, one unordered, one ordered, and one continuous variable is simply the product of l u (•), l o (•), and w(•), and multivariate versions with differing numbers of each datatype are created in the same fashion.We naturally allow the bandwidths to differ for each variable.For what follows, for a vector of mixed variables, we define K γ (x, X i ) to be this product, where γ is the vector of bandwidths (e.g., the hs and λs) and x the vector of mixed datatypes.For further details see Li and Racine (2003) who proposed the use of these generalized product kernels for unconditional density estimation and developed the underlying theory for a data-driven method of bandwidth selection for this class of estimators.The use of such kernels offers a seamless framework for kernel methods with mixed data.For further details on a range of kernel methods that employ this approach, we direct the interested reader to Li and Racine (2007a, Chapter 4).
If the user wishes to apply kernel functions other than those provided by the default settings, the kernel function can be changed by passing the appropriate arguments to ukertype, okertype, and ckertype (the first letter of each representing unordered, ordered, and continuous, respectively), while the 'order' for the continuous kernel (i.e., the first non-zero moment) can be changed by passing the appropriate (even) integer to ckerorder.We support a variety of unordered, ordered, and continuous kernels along with a variety of high-order kernels, the default being ckerorder = 2. Using these arguments, the user could select an eighth order Epanechnikov kernel, a fourth order Gaussian kernel, or even a uniform kernel for the continuous variables in their application.See ?np for further details.

Nonparametric regression
We shall start with the workhorse of applied data analysis, namely, regression models.In order to introduce nonparametric regression, we shall first consider the simplest possible regression model, one that involves one continuous dependent variable, y, and one continuous explanatory variable, x.We shall begin with a popular parametric model of a wage equation, and then move on to a fully nonparametric regression model.Having compared and contrasted the parametric and nonparametric approach towards univariate regression, we then proceed to multivariate regression.

Univariate regression
We begin with a classic dataset taken from Pagan and Ullah (1999, p. 155)  If we first consider the model's goodness of fit, the model has an unadjusted R 2 of 23.1%.Next, we consider the local linear nonparametric method proposed by Li and Racine (2004) employing cross-validated bandwidth selection using the method of Hurvich et al. (1998) Using the measure of goodness of fit introduced in the next section, we see that this method produces a better in-sample model, at least as measured by the R 2 criterion, having an R 2 of 32.5%. 3o far we have summarized the model's goodness-of-fit.However, econometricians also routinely report the results from tests of significance.There exist nonparametric counterparts to these tests that were proposed by Racine (1997), and extended to admit categorical variables by Racine et al. (2006) Having conducted this test we observe that, as was the case for the linear parametric model, the explanatory variable age is significant at all conventional levels in the local linear nonparametric model.

Assessing goodness of fit for nonparametric models
The reader may have wondered what the difference is between the R 2 measures reported by the linear and nonparametric models summarized above, or perhaps how the R 2 was generated for the nonparametric model.It is desirable to use a unit-free measure of goodness-of-fit for nonparametric regression models which is comparable to that used for parametric regression models, namely R 2 .Note that this will be a within-sample measure of goodness-of-fit.Given the known drawbacks of computing R 2 based on the decomposition of the sum of squares (such as possible negative values), there is an alternative definition and method for computing R 2 which can be used that is directly applicable to any model, linear or nonlinear.Letting y i denote the outcome and ŷi the fitted value for observation i, we may define R 2 as follows: and this measure will always lie in the range [0, 1] with the value 1 denoting a perfect fit to the sample data and 0 denoting no predictive power above that given by the unconditional mean of the outcome.It can be demonstrated that this method of computing R 2 is identical to the standard measure computed as n i=1 (ŷ i − ȳ) 2 / n i=1 (y i − ȳ) 2 when the model is linear, fitted with least squares, and includes an intercept term.This useful measure will permit direct comparison of within-sample goodness-of-fit subject to the obvious qualification that this is by no means a model selection criterion, rather, simply a summary measure that some may wish to report.This measure can, of course, also be computed using out-ofsample predictions and out-of-sample realizations.Were we to consider models estimated on a randomly selected subset of data and evaluated on an independent sample of hold-out data, this measure computed for the hold-out observations might serve to guide model selection, particularly when averaged over a number of independent hold-out datasets.

Graphical comparison of parametric and nonparametric models
Often, a suitable graphical comparison of models allows the user to immediately appreciate features present in both the data and the estimated models. 4The upper left plot in Figure 1 presents the fitted parametric and nonparametric regression models along with the data for the cps71 example.The following code can be used to construct the plots appearing in Figure 1.
R> plot(cps71$age, cps71$logwage, xlab = "age", ylab = "log(wage)", cex=.1)R> lines(cps71$age, fitted(model.np),lty = 1, col = "blue") 4 By suitable, we mean those that display both point estimates and variability estimates.q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 20 30 40 The upper right plot in Figure 1 compares the gradients for the parametric and nonparametric models.Note that the gradient for the parametric model will be given by ∂logwage i /∂age i = ∂ β2 + 2 β3 age i , i = 1, . . ., n, as the model is quadratic in age.
Of course, it might be preferable to also plot error bars for the estimates, either asymptotic or resampled.We have automated this in the function npplot which is automatically called by plot.The lower left and lower right plots in Figure 1 present pointwise error bars using asymptotic standard error formulas for the regression function and gradient, respectively.
Often, however, distribution free (bootstrapped) error bounds may be desired, and we allow the user to readily do so as we have written np to leverage the boot package (Canty and Ripley 2008).By default ('iid'), bootstrap resampling is conducted pairwise on (y, X, Z) (i.e., by resampling from rows of the (y, X) data or (y, X, Z) data where appropriate).Specifying the type of bootstrapping as 'inid' admits general heteroskedasticity of unknown form via the wild bootstrap (Liu 1988), though it does not allow for dependence.'fixed' conducts the block bootstrap (Künsch 1989) for dependent data, while 'geom' conducts the stationary bootstrap (Politis and Romano 1994).

Generating predictions from nonparametric models
Once you have obtained appropriate bandwidths and estimated a nonparametric model, generating predictions is straightforward involving nothing more than creating a set of explanatory variables for which you wish to generate predictions.These can lie inside the support of the original data or outside should the user so choose.We have written our routines to support the predict function in R, so by using the newdata = option one can readily generate predictions.It is important to note that typically you do not have the outcome for the evaluation data, hence you need only provide the explanatory variables.However, if by chance you do have the outcome and you provide it, the routine will compute the out-of-sample summary measures of predictive ability.This would be useful when one splits a dataset into two independent samples, estimates a model on one sample, and wishes to assess its performance on the independent hold-out data.
By way of example we consider the cps71 data, generate a set of values for age, two of which lie outside of the support of the data (10 and 70 years of age), and generate the parametric and nonparametric predictions using the generic predict function.

Multivariate regression with qualitative and quantitative data
Based on the presumption that some readers will be unfamiliar with the kernel smoothing of qualitative data, we next consider a multivariate regression example that highlights the potential benefits arising from the use of kernel smoothing methods that smooth both the qualitative and quantitative variables in a particular manner.For what follows, we consider an application taken from Wooldridge (2003, p. 226) that involves multiple regression analysis with qualitative information.
We consider modeling an hourly wage equation for which the dependent variable is log(wage) (lwage) while the explanatory variables include three continuous variables, namely educ (years of education), exper (the number of years of potential experience), and tenure (the number of years with their current employer) along with two qualitative variables, female ('Female'/'Male') and married ('Married'/'Notmarried').For this example there are n = 526 observations.The classical parametric approach towards estimating such relationships requires that one first specify the functional form of the underlying relationship.We start by first modelling this relationship using a simple parametric linear model.By way of example, Wooldridge (2003, p. 222)  This model is, however, restrictive in a number of ways.First, the analyst must specify the functional form (in this case linear) for the continuous variables (educ, exper, and tenure).Second, the analyst must specify how the qualitative variables (female and married) enter the model (in this case they affect the model's intercepts only).Third, the analyst must specify the nature of any interactions among all variables, quantitative and qualitative (in this case, there are none).Should any of these presumptions be incorrect, then the estimated model will be biased and inconsistent potentially leading to faulty inference.
One might next test the null hypothesis that this parametric linear model is correctly specified using the consistent model specification test described in Hsiao et al. (2007) that admits both categorical and continuous data (we have overridden the default search tolerances for this example as the objective function is well-behaved via tol = 0.1 and ftol = 0.1, which should not be done in general -this example will likely take a few minutes on a desktop computer as it uses bootstrapping and cross-validated bandwidth selection): Note that it might appear at first blush that the user needs to do some redundant typing as the X data in this example is the same as the regressors in the model.However, in general X could differ which is why the user must specify this object separately as we cannot assume that the explanatory variables in the model will be equal to X.
This naïve linear model is rejected by the data (the p value for the null of correct specification is < 0.001), hence one might proceed instead to model this relationship using kernel methods.
As noted, the traditional nonparametric approach towards modeling relationships in the presence of qualitative variables requires that you first split your data into subsets containing only the continuous variables of interest (lwage, exper, and tenure).For instance, we would have four such subsets, a) n = 132 observations for married females, b) n = 120 observations for single females, c) n = 86 observations for single males, and d) n = 188 observations for married males.One would then construct smooth nonparametric regression models for each of these subsets and proceed with analysis in this fashion.However, this may lead to a loss in efficiency due to a reduction in the sample size leading to overly variable estimates of the underlying relationship.
Instead, however, we could construct smooth nonparametric regression models by i) using a kernel function that is appropriate for the qualitative variables as outlined in Section 2.5 and ii) modifying the nonparametric regression model as was done by Li and Racine (2004).One can then conduct sound nonparametric estimation based on the n = 526 observations rather than resorting to sample splitting.The rationale for this lies in the fact that doing so may introduce potential bias, however it will always reduce variability thus leading to potential finite-sample efficiency gains.Our experience has been that the potential benefits arising from this approach more than offset the potential costs in finite-sample settings.
Next, we consider using the local-linear nonparametric method described in Li and Racine (2004).For the reader's convenience we supply precomputed cross-validated bandwidths which are automatically loaded when one loads the wage1 dataset (recall being cautioned about the computational burden associated with multivariate data-driven bandwidth methods).
In the example that follows, we use the precomputed bandwidth object bw.all that contains the data-driven bandwidths for the local linear regression produced below using all observations in the sample.
R> bw.all <-npregbw(lwage ~female + married + educ + exper + tenure, + regtype = "ll", bwmethod = "cv.aic",data = wage1) R> model.np<-npreg(bws = bw.all)R> summary(model.np)Note again that the bandwidth object is the only thing you need to pass to npreg as it encapsulates the kernel types, regression method, and so forth.You could also use npreg and manually specify the bandwidths using a bandwidth vector if you so choose. 6We have tried to make each function as flexible as possible to meet the needs of a variety of users.
The goodness of fit of the nonparametric model (R 2 = 51.5%) is better than that for the parametric model (R 2 = 40.4%).In order to investigate whether this apparent improvement reflects overfitting or simply that the nonparametric model is in fact more faithful to the underlying data generating process, we shuffled the data and created two independent samples, one of size n 1 = 400 and one of size n 2 = 126.We fit the models on the n 1 training observations then evaluate the models on the n 2 independent hold-out observations using the predicted square error criterion, namely n −1 , where the y i s are the lwage values for the hold-out observations and the ŷi s are the predicted values.Finally, we compare the parametric model, the nonparametric model that smooths both the categorical and continuous variables, and the traditional frequency nonparametric model that breaks the data into subsets and smooths the continuous data only.For this example we use the precomputed bandwidth object bw.subset which contains the data-driven bandwidths for a random subset of data.
Next, we display partial regression plots in Figure 2. 7 We also plot bootstrapped variability bounds, where the bootstrapping is done via the boot package thereby facilitating a variety of bootstrap methods.The following code will generate Figure 2.
R> plot(model.np,plot.errors.method= "bootstrap", + plot.errors.boot.num= 25) Figure 2 reveals that (holding other regressors constant at their median/mode), males have higher expected wages than females, there does not appear to be a significant difference between the expected wages of married and nonmarried individuals, wages increase with education and tenure and first rise then fall as experience increases, other things equal.

Nonparametric binary outcome and count data models
For what follows, we adopt the conditional probability estimator proposed in Hall et al. (2004) to estimate a nonparametric model of a binary outcome when there exist a number of categorical covariates.
For this example, we use the birthwt data taken from the MASS package (Venables and Ripley 2002), and compute a parametric logit model and a nonparametric conditional mode model.We then compare their confusion matrices8 and assess their classification ability.Note that many of the variables in this dataset have not been classed so we must do this upon invocation of a function.The outcome is an indicator of low infant birthweight (0/1).

Nonparametric unconditional PDF and CDF estimation
The Old Faithful Geyser is a tourist attraction located in Yellowstone National Park.This famous dataset containing n = 272 observations consists of two variables, eruption duration in minutes (eruptions) and waiting time until the next eruption in minutes (waiting).This dataset is used by the park service to model, among other things, expected duration conditional upon the amount of time that has elapsed since the previous eruption.Modeling the joint distribution is, however, of interest in its own right, and the underlying bimodal nature of the joint PDF and CDF is readily revealed by the kernel estimator.For this example, we load the old faithful geyser data and compute the density and distribution functions.To plot the conditional distribution we use cdf = TRUE in plot.Results are presented in Figure 3.Note that in this example we conduct bandwidth selection and estimation in one step.
R> plot(f.faithful,xtrim = -0.2,view = "fixed", main = "") R> plot(f.faithful,cdf = TRUE, xtrim = -0.2,view = "fixed", main = "") If one were to instead model this density with a parametric model such as the bivariate normal (being symmetric, unimodal, and monotonically decreasing away from the mode), one would of course fail to uncover the underlying structure readily revealed by the kernel estimate.

Nonparametric conditional PDF and CDF estimation
We consider Giovanni Baiocchi's (Baiocchi 2006) Italian GDP growth panel for 21 regions covering the period 1951-1998(millions of Lire, 1990 = base) = base).There are n = 1, 008 observations in total, and two variables, gdp and year.First, we compute the bandwidths.Note that this may take a minute or two depending on the speed of your computer.We override the default tolerances for the search method as the objective function is well-behaved (do not of course do this in general), then we compute the npcdens object.To plot the conditional distribution we use cdf = TRUE in plot.Note that in this example we conduct bandwidth selection and estimation in one step.R> plot(fhat, view = "fixed", main = "", theta = 300, phi = 50) R> plot(fhat, cdf = TRUE, view = "fixed", main = "", theta = 300, phi = 50) Figure 4 reveals that the distribution of income has evolved from a unimodal one in the early 1950s to a markedly bimodal one in the 1990s.This result is robust to bandwidth choice, and is observed whether using simple rules-of-thumb or data-driven methods such as likelihood cross-validation.The kernel method readily reveals this evolution which might easily be missed were one to use parametric models of the income distribution (e.g., the unimodal lognormal distribution is commonly used to model income distributions).

Nonparametric quantile regression
We again consider Giovanni Baiocchi's Italian GDP growth panel.First, we compute the likelihood cross-validation bandwidths (default).We override the default tolerances for the search method as the objective function is well-behaved (do not of course do this in general).Then we compute the resulting conditional quantile estimates using the method of Li and Racine (2008).By way of example, we compute the 25th, 50th, and 75th conditional quantiles.Note that this may take a minute or two depending on the speed of your computer.Note that for this example we first call npcdensbw to avoid unnecessary re-computation of the bandwidth object.R> plot(Italy$year, Italy$gdp, main = "", xlab = "Year", + ylab = "GDP Quantiles") R> lines(Italy$year, model.q0.25$quantile, col = "red", lty = 1, lwd = 2) R> lines(Italy$year, model.q0.50$quantile, col = "blue", lty = 2, lwd = 2) R> lines(Italy$year, model.q0.75$quantile, col = "red", lty = 3, lwd = 2) R> legend (ordered(1951), 32, c("tau = 0.25", "tau = 0.50", "tau = 0.75"), + lty = c(1, 2, 3), col = c("red", "blue", "red")) One nice feature of this application is that the explanatory variable is ordered and there exist multiple observations per year.Using the plot function with ordered data produces a boxplot which readily reveals the non-smooth 25th, 50th, and 75th quantiles.These non-smooth quantile estimates can then be directly compared to those obtained via direct estimation of the smooth CDF which are plotted in Figure 5.

Semiparametric partially linear models
Suppose that we consider the wage1 dataset from Wooldridge (2003, p. 222), but now assume that the researcher is unwilling to presume the nature of the relationship between exper and lwage, hence relegates exper to the nonparametric part of a semiparametric partially linear model.The partially linear model was proposed by Robinson (1988) and extended to handle the presence of categorical covariates by Racine and Liu (2007).
Before proceeding, we ought to clarify a common misunderstanding about partially linear models.Many believe that, as the model is apparently simple, its computation ought to also be simple.However, the apparent simplicity hides the perhaps under-appreciated fact that bandwidth selection for partially linear models can be orders of magnitude more computationally burdensome than that for fully nonparametric models, for one simple reason.

Continuous outcomes (Ichimura with cross-validation)
Next, we consider applying Ichimura (1993)'s single-index method which is appropriate for continuous outcomes, unlike that of Klein and Spady (1993) (we override the default number of multistarts for the user's convenience as the global minimum appears to have been located in the first attempt).We again consider the wage1 dataset found in Wooldridge (2003, p. 222) It is interesting to compare this model with the parametric and nonparametric models presented in Section 3.2 as it provides an in-sample fit (43.1%) that lies in between the misspecified parametric model (40.4%) and fully nonparametric model (51.5%).Whether this model yields improved out-of-sample predictions could also be explored.

Semiparametric varying coefficient models
We revisit the wage1 dataset found in Wooldridge (2003, p. 222), but assume that the researcher believes that the model's parameters may differ depending on one's sex.In this case, one might adopt a varying coefficient approach such as that found in Li and Racine (2007b) and Racine et al. (2007).We compare a simple linear model with the semiparametric varying coefficient model.Note that the X data in the varying coefficient model must be of type numeric, so we create a 0/1 dummy variable from the qualitative variable for X, but of course for the nonparametric component we can simply treat these as unordered factors.Note that we do bandwidth selection and estimation in one step.

Writing your own kernel-based functions
The function npksum exists so that you can create your own kernel objects with or without a variable to be weighted (defaults to 1).With the options available, you could create new nonparametric tests or even new kernel estimators.The convolution kernel option would allow you to create, say, the least squares cross-validation function for kernel density estimation.npksum implements a variety of methods for computing multivariate kernel sums (p-variate) defined over a set of possibly continuous and/or discrete data.
By way of example, we construct a local constant kernel estimator with a bandwidth of, say, two years.Figure 6 plots the resulting estimate.
R> plot(cps71$age, cps71$logwage, xlab = "Age", ylab = "log(wage)") R> lines(cps71$age, fit.lc, col = "blue") npksum is exceedingly flexible, allowing for leave-one-out sums, weighted sums of matrices, raising the kernel function to different powers, the use of convolution kernels, and so forth.See ?npksum for further details.q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 20 30 40 50 60

A parallel implementation
Data-driven bandwidth selection methods are, by their nature, computationally burdensome.However, many bandwidth selection methods lend themselves well to parallel computing approaches.High performance computing resources are becoming widely available, and multiple CPU desktop systems have become the norm and will continue to make inroads.
When users have a large amount of data, serial bandwidth selection procedures can be infeasible.For this reason, we have developed an MPI-aware version of the np package that uses some of the functionality of the Rmpi package which we have tentatively called the npRmpi package and is available from the authors upon request (the Message Passing Interface (MPI) is an open library specification for parallel computation available for a range of computing platforms).The functionality of np and npRmpi is identical, however, using npRmpi you could take advantage of a cluster computing environment or a multi-core/multi-CPU desktop machine thereby alleviating the computational burden associated with the nonparametric analysis of large datasets.Installation of this package, however, requires knowledge that goes beyond that which even seasoned R users may possess.Having access to local expertise would be necessary for many users therefore this package is not available via CRAN.

Summary
The np package offers users of R a variety of nonparametric and semiparametric kernel-based methods that are capable of handling the mix of categorical and continuous data typically encountered by applied researchers.In this article we have described the functionality of the np package via a series of illustrative applications that may be of interest to applied econometricians interested in becoming familiar with these methods.We do not delve into details of the underlying estimators, rather we provide references where appropriate and direct the interested reader to those resources.
The help files accompanying many functions found in the np package contain numerous examples which may be of interest to some readers, and we encourage readers to experiment with these examples in addition to those contained herein.
Finally, we encourage readers who have successfully implemented new kernel-based methods using the npksum function to send such functions to us so that they can be included in future versions of the np package with appropriate acknowledgment of course.

Figure 1 :
Figure 1: The figure on the upper left presents the parametric (quadratic, dashed line) and the nonparametric estimates (solid line) of the regression function for the cps71 data.The figure on the upper right presents the parametric (dashed line) and nonparametric estimates (solid line) of the gradient.The figures on the lower left and lower right present the nonparametric estimates of the regression function and gradient along with their variability bounds, respectively.

Figure 2 :
Figure 2: Partial local linear nonparametric regression plots with bootstrapped pointwise error bounds for the wage1 dataset.

Figure 3 :Figure 4
Figure 3: Nonparametric multivariate PDF and CDF estimates for the Old Faithful data.

Figure 4 :
Figure 4: Nonparametric conditional PDF and CDF estimates for the Italian GDP panel.
Figure 5 plots the resulting quantile estimates.The following code will generate Figure 5.

Figure 5 :
Figure 5: Nonparametric quantile regression on the Italian GDP panel.

Figure 6 :
Figure 6: A local constant kernel estimator generated with npksum for the cps71 dataset.
'***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 who consider Canadian cross-section wage data consisting of a random sample taken from the 1971 Canadian Census Public Use Tapes for male individuals having common education (Grade 13).There are n = 205 observations in total, and 2 variables, the logarithm of the individual's wage (logwage) and their age (age).The traditional wage equation is typically modelled as a quadratic in age.First, we begin with a simple parametric model for this relationship.
. Note that in this example we conduct bandwidth selection and estimation in one step.
. Note that in this example we conduct bandwidth selection and estimation in one step.
It is again interesting to compare this model with the parametric and nonparametric models presented in Section 3.2.It can be seen that the average values of the model's coefficients are in agreement with those from the linear specification, while the in-sample goodness of fit (42.6%) lies in between the misspecified parametric model (40.4%) and fully nonparametric model (51.5%).