Semiparametric Single Index Model

npindex computes a semiparametric single index model for a dependent variable and \(p\)-variate explanatory data using the model \(Y = G(X\beta) + \epsilon\), given a set of evaluation points, training points (consisting of explanatory data and dependent data), and a npindexbw bandwidth specification. Note that for this semiparametric estimator, the bandwidth object contains parameters for the single index model and the (scalar) bandwidth for the index function.

Usage

npindex(bws, ..., B = 399)

# S3 method for class 'formula'
npindex(bws,
        data = NULL,
        newdata = NULL,
        y.eval = NULL,
        se = TRUE,
        ..., se.type = c("asymptotic", "bootstrap"))


# Default S3 method
npindex(bws,
        txdat,
        tydat,
        nomad = FALSE,
        se = TRUE,
        ..., se.type = c("asymptotic", "bootstrap"))

# S3 method for class 'sibandwidth'
npindex(bws,
        txdat = stop("training data 'txdat' missing"),
        tydat = stop("training data 'tydat' missing"),
        exdat,
        eydat,
        B = 399,
        se = TRUE,
        gradients = FALSE,
        residuals = FALSE,
        ..., se.type = c("asymptotic", "bootstrap"),
        y.eval = NULL)

Arguments

Data, Bandwidth Inputs And Formula Interface

These arguments identify the bandwidth specification, formula/data interface, and training data.

bws

a bandwidth specification. This can be set as a sibandwidth object returned from an invocation of npindexbw, or as a vector of parameters (beta) with each element \(i\) corresponding to the coefficient for column \(i\) in txdat where the first element is normalized to 1, and a scalar bandwidth (h).

data

This argument belongs to the formula interface (including supported formula-bearing bandwidth-object refits). Native calls must supply their vectors, matrices or data frames directly, or use with(data, ...). Supplying data, including explicit data = NULL, to a native route that does not consume it is an error; see np-data-interfaces.

an optional data frame, list or environment (or object coercible to a data frame by as.data.frame) containing the variables in the model. If not found in data, the variables are taken from environment(bws), typically the environment from which npindexbw was called.

txdat

a \(p\)-variate data frame of explanatory data (training data) used to calculate the regression estimators. Defaults to the corresponding training data recoverable from the bandwidth object in the current session.

tydat

a one (1) dimensional numeric or integer vector of dependent data, each element \(i\) corresponding to each observation (row) \(i\) of txdat. Defaults to the corresponding training data recoverable from the bandwidth object in the current session. For method="kleinspady", response values must be numeric 0 or 1. Factor responses are rejected rather than converted to internal category codes; explanatory variables may still be factors. The same response restriction applies to supplied evaluation outcomes.

Local-Polynomial Degree And Bandwidth Search

This argument controls the recommended automatic local-polynomial NOMAD route, which jointly selects continuous polynomial degree and bandwidths when these are computed inside npindex.

nomad

logical shortcut passed through to npindexbw when bandwidths are computed inside npindex. When TRUE, the single-index bandwidth route fills any missing values among regtype, search.engine, degree.select, bernstein.basis, degree.min, degree.max, degree.verify, and bwtype with the recommended automatic local-polynomial degree-and-bandwidth NOMAD preset documented in npindexbw. Additional NOMAD tuning arguments such as nomad.nmulti may also be supplied through ...; nmulti remains the outer restart count while nomad.nmulti controls inner native crs NOMAD API multistarts within each outer restart. After fitting, inspect fit$bws$nomad.shortcut on the returned object fit to see the normalized shortcut metadata.

Evaluation Data And Returned Quantities

These arguments control where the single-index fit is evaluated and which evaluation quantities are returned.

exdat

a \(p\)-variate data frame of points on which the regression will be estimated (evaluation data). By default, evaluation takes place on the data provided by txdat. If all computed kernel weights are zero at an external evaluation row, the undefined fitted quantities are returned as NA, with one notice per call. This policy is independent of the se and gradients switches; required training computations remain strict. If a bootstrap resample has no support at a requested row, its bootstrap standard error is NA; a supported original point estimate is retained. Replications are not dropped or redrawn.

eydat

a one (1) dimensional numeric or integer vector of the true values of the dependent variable. Optional, and used only to calculate the true errors.

newdata

An optional data frame in which to look for evaluation data. If omitted, the training data are used.

y.eval

Selects the outcomes used for goodness-of-fit diagnostics. The default NULL uses evaluation outcomes when supplied in newdata or eydat. With evaluation predictors only, predictions are returned normally and diagnostics are NA. Explicit TRUE requires evaluation outcomes when evaluating new predictors. Explicit FALSE ignores outcomes in newdata and preserves historical training diagnostics, identified as such in the summary. Explicit eydat takes precedence. With no evaluation inputs, diagnostics describe training data.

Fitted Quantities And Inference

These arguments control residuals, gradients, and asymptotic or bootstrap inference.

B

an integer specifying the number of bootstrap replications to use when se=TRUE, se.type="bootstrap". Defaults to 399.

se

a logical value requesting asymptotic coefficient covariance and standard errors for the conditional mean and, when gradients=TRUE, gradients. Defaults to TRUE for this estimator. Set FALSE to omit all covariance and standard-error calculations, independently of whether gradient estimates are requested.

se.type

the standard-error method for the conditional mean and gradients: "asymptotic" (the default) or "bootstrap". The latter also supplies bootstrap standard errors for average gradients and uses B replications. Coefficient covariance is asymptotic in both modes: the fitted-value/gradient bootstrap holds beta and the selected bandwidth fixed while resampling observations. No bootstrap is run when se=FALSE or se.type="asymptotic".

gradients

a logical value indicating that you want gradients computed and returned. Defaults to FALSE. Their standard errors additionally require se=TRUE. Coefficient covariance does not require this flag.

residuals

a logical value indicating that you want residuals computed and returned in the resulting singleindex object. Defaults to FALSE.

Additional Arguments

Further arguments are passed to the bandwidth-selection counterpart when bandwidths are not supplied.

...

additional arguments supplied to specify the parameters to the sibandwidth S3 method, which is called during estimation.

Details

Factor Predictors and Coefficient Inference

New bandwidth objects expand factors into the retained R contrast design described in npindexbw. Training refits, prediction and plots reuse the same columns, levels and contrasts. Numeric scores and older saved bandwidth objects retain their existing interpretation. Supply raw predictor data to the public fitting and prediction interfaces; expansion is automatic.

Coefficients, covariance and gradients use expanded design-column names. Gradients for factor contrast columns are derivatives of the numerical design extension, \(g'(X^\prime\beta)\beta_j\); they are not the effect of switching between observed categories. Such effects require comparing predictions at the desired categories. No coefficient is estimated for the intercept, which is absorbed into the unknown link.

Klein–Spady coefficient covariance uses the link derivative multiplied by covariates centered on their conditional means given the index. This corrects a legacy raw-covariate information calculation, including for numeric-only inputs. Ichimura retains its conditional-centered sandwich calculation. Coefficient inference requires additional conditional-moment computation; se=FALSE skips this work. This correction does not change the search objective or the conditional fitted-value uncertainty calculation.

Fixed kernel endpoints retained in the bandwidth object are also used by fitting, prediction, inference and bootstrap resamples. They are not recomputed from the evaluation grid or from a resampled training set.

With a formula-derived bandwidth object, supplying only tydat replaces the evaluated response in retained-training-sample order: one row per observation retained after the original subset and missing-value selection. Original excluded rows are not restored; row names or equal lengths do not imply original-data order. The effective na.action is then applied jointly to the retained predictors and replacement response; any new omission indices refer to this retained sample. The selected bandwidths remain fixed, and the returned object retains the replacement sample for subsequent refits. Use data for full formula-data replacement, not together with a partial response replacement. Explicit native training arguments retain their separate positional contract.

For native (non-formula) objects, named newdata columns are matched to the retained variable names and reordered as needed. Extra columns are ignored; missing, duplicate or ambiguous names are errors. An unnamed matrix or vector remains positional for a single evaluation role. Explicit native evaluation arguments retain their positional contract and take precedence over newdata.

For a native bandwidth object, newdata supplies exdat and, in automatic mode, a named response column supplies evaluation outcomes. Explicit native evaluation arguments take precedence. This evaluates at the requested rows without repeating bandwidth selection or re-estimating index coefficients. For held-out assessment, estimate the coefficients and any data-selected bandwidths using the training sample.

For Ichimura fits trained on an unordered factor response (or an ordered response without numeric levels), evaluation outcomes must also be factors compatible with the training levels: raw numeric labels do not identify the factor-code scale used for fitting. In automatic mode such numeric outcomes leave diagnostics unavailable, without blocking predictions; explicit evaluation-outcome requests report an error. For ordered training responses with numeric levels, numeric evaluation outcomes on that same scale are also accepted.

Evaluation outcomes affect diagnostics only. In automatic mode, outcomes that cannot be evaluated, aligned, or validated are treated as unavailable for diagnostics; predictions remain available. Explicit y.eval=TRUE and explicitly supplied eydat retain strict validation. Missing outcomes do not remove predictor-complete prediction rows. Time-series outcomes are matched to the predictor time coordinates; predictions outside the supplied outcome time support remain available but are not scored. Scores use observed outcome/prediction pairs, and are NA when none are available. Formula responses, including transformations, are evaluated only when their variables are present in newdata; an ambient response is never substituted. Predictor omission and row-restoration conventions are unchanged.

Formula variables of class ts, including lagged variables, are aligned on their common time support before fitting. Other variables must already have the length of that intersection. Supply explicitly aligned rows when using native data arguments instead of a formula.

Documentation guide: see np.kernels for kernels, np.options for global options, and plot, plot.np for plotting options.

For S3 plotting help, see plot.np. You can list available plot methods with methods("plot").

A matrix of gradients along with average derivatives are computed and returned if gradients=TRUE is used.

Interpreting \(\beta\) Coefficients

Single-index coefficients are identified only up to scale, so npindexbw normalizes the first coefficient to one. The remaining coefficients describe relative movement in the index \(X^\prime\beta\), not standalone linear-model slopes. Holding other variables fixed, a one-unit increase in \(X_j\) changes the index by the same amount as a \(\beta_j\)-unit increase in \(X_1\). Where the link derivative exists and is nonzero, this is also the local marginal-effect ratio relative to \(X_1\). Choose the first continuous regressor deliberately when coefficient interpretation is important.

For \(E[Y\mid X]=g(X^\prime\beta)\), the component plotted for predictor \(X_j\) is $$\frac{\partial E[Y\mid X]}{\partial X_j} = \beta_j g'(X^\prime\beta)$$ Accordingly, plot(..., gradients=TRUE) displays every component against the common scalar index \(X^\prime\beta\); the y-axis identifies the predictor with respect to which the derivative is taken. These panels are not partial-dependence plots against the individual \(X_j\) values.

Optimizer Choice When Bandwidths Are Computed

When npindex computes its bandwidth object internally, optimizer controls in ... are passed to npindexbw. The default optim.method="Nelder-Mead" is a robust derivative-free starting point for low-dimensional searches, such as one bandwidth plus only a few index coefficients. For larger \(p\), optim.method="BFGS" can be useful as a performance and sensitivity comparison, but users should compare objective values, convergence behavior, and fitted objects because the single-index criterion can be nonconvex and route-sensitive.

For practitioners who want the recommended automatic local-polynomial degree-and-bandwidth NOMAD route without spelling out all LP tuning arguments, npindex(..., nomad=TRUE) and npindexbw(..., nomad=TRUE) expand missing settings to the same documented preset. Explicit incompatible settings fail fast rather than being silently rewritten.

Value

npindex returns a singleindex object. The generic functions fitted, residuals, coef, vcov, se, predict, and gradients, extract (or generate) estimated values, residuals, coefficients, variance-covariance matrix, asymptotic or bootstrap standard errors on estimates, predictions, and gradients, respectively, from the returned object. Furthermore, the functions summary and plot support objects of this type. The returned object has the following components:

eval

evaluation points

mean

estimates of the regression function (conditional mean) at the evaluation points

beta

the model coefficients

betavcov

the asymptotic covariance matrix for the model coefficients

merr

standard errors of the regression function estimates

grad

estimates of the gradients at each evaluation point

gerr

standard errors of the gradient estimates

mean.grad

mean (average) gradient over the evaluation points

mean.gerr

bootstrapped standard error of the mean gradient estimates

R2

if method="ichimura", coefficient of determination (Doksum and Samarov (1995))

MSE

if method="ichimura", mean squared error

MAE

if method="ichimura", mean absolute error

MAPE

if method="ichimura", mean absolute percentage error

CORR

if method="ichimura", absolute value of Pearson's correlation coefficient

SIGN

if method="ichimura", fraction of observations where fitted and observed values agree in sign

diagnostics.sample

"training", "evaluation", or "unavailable", identifying the sample used for goodness-of-fit statistics.

diagnostics.nobs

Number of observed outcome/prediction pairs scored.

training.MSE

For Ichimura fits, the training mean squared error when training fitted values were computed, otherwise NA. Prediction's residual.scale retains this training quantity in its historical MSE units, independently of evaluation outcomes. No additional training fit is performed solely to populate this field.

confusion.matrix

if method="kleinspady", the confusion matrix or NA if outcomes are not available

CCR.overall

if method="kleinspady", the overall correct classification ratio, or NA if outcomes are not available

CCR.byoutcome

if method="kleinspady", a numeric vector containing the correct classification ratio by outcome, or NA if outcomes are not available

fit.mcfadden

if method="kleinspady", the McFadden-Puig-Kerschner performance measure or NA if outcomes are not available

Book And Method Pointers

The single-index model reduces a multivariate predictor to an index, typically written \(E[Y\mid X]=g(X^\prime\beta)\) for continuous outcomes, with a normalization on \(\beta\) for identification. For binary outcomes, the Klein-Spady route estimates the corresponding binary-choice probability through the index. The fitted index regression, gradients, and average derivatives are computed from the selected index direction and kernel regression fit.

For book-length derivations, see Li and Racine (2007), Chapter 8 Semiparametric Single Index Models, especially Sections 8.1, 8.2.1, 8.4, 8.5, and 8.10. The later workflow treatment is Racine (2019), Chapter 8 Semiparametric Conditional Mean Function Estimation, especially the single-index material.

References

Aitchison, J. and C.G.G. Aitken (1976), “Multivariate binary discrimination by the kernel method,” Biometrika, 63, 413-420.

Doksum, K. and A. Samarov (1995), “Nonparametric estimation of global functionals and a measure of the explanatory power of covariates regression,” The Annals of Statistics, 23 1443-1473.

Ichimura, H., (1993), “Semiparametric least squares (SLS) and weighted SLS estimation of single-index models,” Journal of Econometrics, 58, 71-120.

Klein, R. W. and R. H. Spady (1993), “An efficient semiparametric estimator for binary response models,” Econometrica, 61, 387-421.

Li, Q. and J.S. Racine (2007), Nonparametric Econometrics: Theory and Practice, Princeton University Press.

McFadden, D. and C. Puig and D. Kerschner (1977), “Determinants of the long-run demand for electricity,” Proceedings of the American Statistical Association (Business and Economics Section), 109-117.

Wang, M.C. and J. van Ryzin (1981), “A class of smooth estimators for discrete distributions,” Biometrika, 68, 301-309.

Author

Tristen Hayfield tristen.hayfield@gmail.com, Jeffrey S. Racine racinej@mcmaster.ca

Usage Issues

If you are using data of mixed types, then it is advisable to use the data.frame function to construct your input data and not cbind, since cbind will typically not work as intended on mixed data types and will coerce the data to the same type.

For formula fits, character-valued missing-data policies are resolved through the stats model-frame owner; pass a caller-local policy as a function object rather than its name. With replacement training data, an explicit na.action, including NULL, takes precedence over the saved policy. For formula evaluation with newdata, an explicit na.action also applies to the evaluation predictors: na.omit drops incomplete predictor rows, na.exclude restores their positions with NA in predictions, standard errors and requested gradients, and na.fail rejects them. When omitted, evaluation uses the current model-frame missing-data policy. Evaluation outcomes are scored separately and do not remove prediction rows. With both replacement training data and evaluation data, the policy is applied separately to each frame. An evaluation-only request does not refilter retained training data. Explicit native evaluation arguments retain precedence and their existing complete-row handling. Native fitting still requires complete rows, so a permissive policy does not guarantee the same retained sample as stats::model.frame.

When reusing a bandwidth object without explicit training or evaluation data, newly constructed bandwidth objects use retained training values. A fit with explicit replacement training data retains that sample in its returned bandwidth object; evaluation grids do not replace training data. Older objects without retained samples use their stored call environment when available and may require explicit training inputs after saveRDS() and readRDS(). Custom formula transformations used on new evaluation data must still be available.

vcov extracts the asymptotic coefficient covariance computed when se=TRUE, independently of gradients and se.type. Its first row and column are zero because the first coefficient is normalized to one. A singular free-coefficient information matrix raises an informative error; the calculation does not substitute a generalized inverse or bootstrap. In particular, the local-constant uniform kernel has zero first derivative away from its support boundaries and cannot supply this covariance for free coefficients. Use se=FALSE for point estimates in that case. Selecting se.type="bootstrap" does not replace the coefficient covariance calculation, which remains asymptotic. The default inference requires additional training derivatives and moments; it can cost appreciably more than a point-estimate-only fit, especially when bandwidths have already been selected. Use se=FALSE to omit that work. Asymptotic fitted-value and gradient standard errors use the same scalar-index regression calculation as the asymptotic plot intervals, treating the fitted index coefficients as fixed. They inherit the residual normalization of Ruppert et al. (1997) described under Standard Errors in npreg. This is not the separate coefficient covariance returned by vcov(). Missing training residual information can therefore make individual fitted-value or gradient errors unavailable without removing their point estimates. Local-polynomial index fits also inherit the numerical-ridging policy described in npreg. Asymptotic standard errors of average gradients are not supplied; mean.gerr is available only in bootstrap mode.

Missing quantities can be computed without repeating bandwidth search. After model <- npindex(...), use model <- npindex(bws=model$bws, gradients=TRUE, se=TRUE) to add gradients and asymptotic inference, or se=FALSE for gradients alone. This reuses both the fitted beta and bandwidth, just as fitting an object returned by npindexbw() does. The training data must remain recoverable from the retained state/call, or be supplied explicitly. predict(..., se.fit=TRUE) uses the fitted object's stored se.type for new evaluation data unless an explicit se.type is supplied.

See also

np.kernels, np.options, plot, plot.np, npindexbw

Examples

if (FALSE) { # \dontrun{
# EXAMPLE 1 (INTERFACE=FORMULA): Generate a simple linear model then
# estimate it using a semiparametric single index specification and
# Ichimura's nonlinear least squares coefficients and bandwidth
# (default). Also compute the matrix of gradients and average derivative
# estimates.

set.seed(12345)

n <- 100

x1 <- runif(n, min=-1, max=1)
x2 <- runif(n, min=-1, max=1)

y <- x1 - x2 + rnorm(n)

# Note - this may take a minute or two depending on the speed of your
# computer. Note also that the first element of the vector beta is
# normalized to one for identification purposes, and that X must contain
# at least one continuous variable.

bw <- npindexbw(formula=y~x1+x2)

summary(bw)

model <- npindex(bws=bw, gradients=TRUE)

summary(model)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# Or you can visualize the input with plot.

if (interactive()) plot(bw)

if (interactive()) Sys.sleep(5)

# EXAMPLE 1 (INTERFACE=DATA FRAME): Generate a simple linear model then
# estimate it using a semiparametric single index specification and
# Ichimura's nonlinear least squares coefficients and bandwidth
# (default). Also compute the matrix of gradients and average derivative
# estimates.

set.seed(12345)

n <- 100

x1 <- runif(n, min=-1, max=1)
x2 <- runif(n, min=-1, max=1)

y <- x1 - x2 + rnorm(n)

X <- cbind(x1, x2)

# Note - this may take a minute or two depending on the speed of your
# computer. Note also that the first element of the vector beta is
# normalized to one for identification purposes, and that X must contain
# at least one continuous variable.

bw <- npindexbw(xdat=X, ydat=y)

summary(bw)

model <- npindex(bws=bw, gradients=TRUE)

summary(model)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# Or you can visualize the input with plot.

if (interactive()) plot(bw)

if (interactive()) Sys.sleep(5)

# EXAMPLE 2 (INTERFACE=FORMULA): Generate a simple binary outcome linear
# model then estimate it using a semiparametric single index
# specification and Klein and Spady's likelihood-based coefficients and
# bandwidth (default). Also compute the matrix of gradients and average
# derivative estimates.

n <- 100

x1 <- runif(n, min=-1, max=1)
x2 <- runif(n, min=-1, max=1)

y <- ifelse(x1 + x2 + rnorm(n) > 0, 1, 0)

# Note that the first element of the vector beta is normalized to one
# for identification purposes, and that X must contain at least one
# continuous variable.

bw <- npindexbw(formula=y~x1+x2, method="kleinspady")

summary(bw)

model <- npindex(bws=bw, gradients=TRUE)

# Note that, since the outcome is binary, we can assess model
# performance using methods appropriate for binary outcomes. We look at
# the confusion matrix, various classification ratios, and McFadden et
# al's measure of predictive performance.

summary(model)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# EXAMPLE 2 (INTERFACE=DATA FRAME): Generate a simple binary outcome
# linear model then estimate it using a semiparametric single index
# specification and Klein and Spady's likelihood-based coefficients and
# bandwidth (default). Also compute the matrix of gradients and average
# derivative estimates.

n <- 100

x1 <- runif(n, min=-1, max=1)
x2 <- runif(n, min=-1, max=1)

y <- ifelse(x1 + x2 + rnorm(n) > 0, 1, 0)

X <- cbind(x1, x2)

# Note that the first element of the vector beta is normalized to one
# for identification purposes, and that X must contain at least one
# continuous variable.

bw <- npindexbw(xdat=X, ydat=y, method="kleinspady")

summary(bw)

model <- npindex(bws=bw, gradients=TRUE)

# Note that, since the outcome is binary, we can assess model
# performance using methods appropriate for binary outcomes. We look at
# the confusion matrix, various classification ratios, and McFadden et
# al's measure of predictive performance.

summary(model)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# EXAMPLE 3 (INTERFACE=FORMULA): Replicate the DGP of Klein & Spady
# (1993) (see their description on page 405, pay careful attention to
# footnote 6 on page 405).

set.seed(123)

n <- 1000

# x1 is chi-squared having 3 df truncated at 6 standardized by
# subtracting 2.348 and dividing by 1.511

x <- rchisq(n, df=3)
x1 <- (ifelse(x < 6, x, 6) - 2.348)/1.511

# x2 is normal (0, 1) truncated at +- 2 divided by 0.8796

x <- rnorm(n)
x2 <- ifelse(abs(x) < 2 , x, 2) / 0.8796

# y is 1 if y* > 0, 0 otherwise.

y <- ifelse(x1 + x2 + rnorm(n) > 0, 1, 0)

# Compute the parameter vector and bandwidth. Note that the first
# element of the vector beta is normalized to one for identification
# purposes, and that X must contain at least one continuous variable.


bw <- npindexbw(formula=y~x1+x2, method="kleinspady")

# Next, create the evaluation data in order to generate a perspective
# plot

# Create an evaluation data matrix

x1.seq <- seq(min(x1), max(x1), length=50)
x2.seq <- seq(min(x2), max(x2), length=50)
X.eval <- expand.grid(x1=x1.seq, x2=x2.seq)

# Now evaluate the single index model on the evaluation data

fit <- fitted(npindex(exdat=X.eval,
               eydat=rep(1, nrow(X.eval)),
               bws=bw))

# Finally, coerce the fitted model into a matrix suitable for 3D
# plotting via persp()

fit.mat <- matrix(fit, 50, 50)

# Generate a perspective plot similar to Figure 2 b of Klein and Spady
# (1993)

persp(x1.seq,
      x2.seq,
      fit.mat,
      col="white",
      ticktype="detailed",
      expand=0.5,
      axes=FALSE,
      box=FALSE,
      main="Estimated Semiparametric Probability Perspective",
      theta=310,
      phi=25)

# EXAMPLE 3 (INTERFACE=DATA FRAME): Replicate the DGP of Klein & Spady
# (1993) (see their description on page 405, pay careful attention to
# footnote 6 on page 405).

set.seed(123)

n <- 1000

# x1 is chi-squared having 3 df truncated at 6 standardized by
# subtracting 2.348 and dividing by 1.511

x <- rchisq(n, df=3)
x1 <- (ifelse(x < 6, x, 6) - 2.348)/1.511

# x2 is normal (0, 1) truncated at +- 2 divided by 0.8796

x <- rnorm(n)
x2 <- ifelse(abs(x) < 2 , x, 2) / 0.8796

# y is 1 if y* > 0, 0 otherwise.

y <- ifelse(x1 + x2 + rnorm(n) > 0, 1, 0)

# Create the X matrix

X <- cbind(x1, x2)

# Compute the parameter vector and bandwidth. Note that the first
# element of the vector beta is normalized to one for identification
# purposes, and that X must contain at least one continuous variable.


bw <- npindexbw(xdat=X, ydat=y, method="kleinspady")

# Next, create the evaluation data in order to generate a perspective
# plot

# Create an evaluation data matrix

x1.seq <- seq(min(x1), max(x1), length=50)
x2.seq <- seq(min(x2), max(x2), length=50)
X.eval <- expand.grid(x1=x1.seq, x2=x2.seq)

# Now evaluate the single index model on the evaluation data

fit <- fitted(npindex(exdat=X.eval,
               eydat=rep(1, nrow(X.eval)),
               bws=bw))

# Finally, coerce the fitted model into a matrix suitable for 3D
# plotting via persp()

fit.mat <- matrix(fit, 50, 50)

# Generate a perspective plot similar to Figure 2 b of Klein and Spady
# (1993)

persp(x1.seq,
      x2.seq,
      fit.mat,
      col="white",
      ticktype="detailed",
      expand=0.5,
      axes=FALSE,
      box=FALSE,
      main="Estimated Semiparametric Probability Perspective",
      theta=310,
      phi=25)
} # }