Nonparametric Instrumental Regression

npregiv computes nonparametric estimation of an instrumental regression function $\varphi$ defined by conditional moment restrictions stemming from a structural econometric model: $E [Y - \varphi (Z,X) | W ] = 0$, and involving endogenous variables $Y$ and $Z$ and exogenous variables $X$ and instruments $W$. The function $\varphi$ is the solution of an ill-posed inverse problem.

When method="Tikhonov", npregiv uses the approach of Darolles, Fan, Florens and Renault (2011) modified for local polynomial kernel regression of any order (Darolles et al use local constant kernel weighting which corresponds to setting p=0; see below for details). When method="Landweber-Fridman", npregiv uses the approach of Horowitz (2011) again using local polynomial kernel regression (Horowitz uses B-spline weighting).

Usage

npregiv(y, ...)

# Default S3 method
npregiv(y,
        z,
        w,
        x = NULL,
        zeval = NULL,
        xeval = NULL,
        alpha = NULL,
        alpha.iter = NULL,
        alpha.max = 1e-01,
        alpha.min = 1e-10,
        alpha.tol = .Machine$double.eps^0.25,
        bw = NULL,
        constant = 0.5,
        iterate.diff.tol = 1.0e-08,
        iterate.max = 1000,
        iterate.Tikhonov = TRUE,
        iterate.Tikhonov.num = 1,
        method = c("Landweber-Fridman","Tikhonov"),
        nmulti = NULL,
        optim.abstol = .Machine$double.eps,
        optim.maxattempts = 10,
        optim.maxit = 500,
        optim.method = c("Nelder-Mead", "BFGS", "CG"),
        optim.reltol = sqrt(.Machine$double.eps),
        p = 1,
        penalize.iteration = TRUE,
        random.seed = 42,
        return.weights.phi = FALSE,
        return.weights.phi.deriv.1 = FALSE,
        return.weights.phi.deriv.2 = FALSE,
        smooth.residuals = TRUE,
        start.from = c("Eyz","EEywz"),
        starting.values  = NULL,
        stop.on.increase = TRUE,
        regtype = NULL,
        degree = NULL,
        nomad = FALSE,
        ...)

# S3 method for class 'formula'
npregiv(y,
        data = NULL,
        subset,
        na.action,
        newdata = NULL,
        zeval = NULL,
        xeval = NULL,
        ...)

# S3 method for class 'npregiv'
plot(x, gradients = FALSE,
        data.overlay = TRUE, data.rug = FALSE, ...)

Arguments

Data, Bandwidth Inputs And Formula Interface

These arguments identify the response, endogenous variables, instruments, exogenous covariates, and evaluation data.

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 used with the formula interface. Formula roles are written explicitly as y ~ z | w or y ~ z | w | x.

subset

an optional expression selecting training rows in the formula interface.

na.action

a function specifying missing-value handling for the formula model frame.

newdata

an optional data frame supplying fit-time evaluation values for the formula z and, when present, x partitions. Explicit zeval and xeval take precedence. This is fit-time evaluation, not a post-fit predict() method.

w

a $q$-variate data frame of instruments. The data types may be continuous, discrete (unordered and ordered factors), or some combination thereof.

x

for npregiv.default, an $r$-variate data frame of exogenous regressors. The data types may be continuous, discrete (unordered and ordered factors), or some combination thereof. For plot.npregiv, an npregiv object.

xeval

an $r$-variate data frame of exogenous regressors on which the regression will be estimated (evaluation data). By default, evaluation takes place on the data provided by x.

y

a one (1) dimensional numeric or integer vector of dependent data, or a formula using the explicit IV grammar described under data. Each response element $i$ corresponds to row $i$ of z.

z

a $p$-variate data frame of endogenous regressors. The data types may be continuous, discrete (unordered and ordered factors), or some combination thereof.

zeval

a $p$-variate data frame of endogenous regressors on which the regression will be estimated (evaluation data). By default, evaluation takes place on the data provided by z.

gradients

logical; display the stored first partial derivative instead of the structural function. The current plot method supports one endogenous coordinate. The default is FALSE.

data.overlay

logical; overlay the training response on a structural-function display. The default is TRUE, as for plot.npregression. It is ignored for derivative displays because response values and derivatives do not share a common vertical scale. When limits are not supplied, the structural curve and response are both included in the displayed range.

data.rug

logical; add a rug showing the training z support. The default is FALSE. A rug may accompany either structural-function or derivative displays; categorical endogenous coordinates omit the numerical rug.

Landweber-Fridman Iteration Controls

These arguments control the Landweber-Fridman iteration path.

constant: the constant to use when using method="Landweber-Fridman".
iterate.diff.tol: the search tolerance for the difference in the stopping rule from iteration to iteration when using method="Landweber-Fridman" (disable by setting to zero).
iterate.max: an integer indicating the maximum number of iterations permitted before termination occurs when using method="Landweber-Fridman".
iterate.Tikhonov: a logical value controlling whether alpha is reselected after the first response-adaptive bandwidth/operator refinement. It does not change the order of the fixed order-two criterion used to select alpha.
iterate.Tikhonov.num: a positive integer giving the number of outer response-adaptive bandwidth/operator refinement cycles. It is not the order of iterated Tikhonov regularization.
method: the regularization method employed (defaults to "Landweber-Fridman", see Horowitz (2011); see Darolles, Fan, Florens and Renault (2011) for details for "Tikhonov").
nmulti: integer number of times to restart the process of finding extrema of the cross-validation function from different (random) initial points.

Optimization Controls

These arguments control numerical optimization for the inverse problem.

optim.abstol

the absolute convergence tolerance used by optim. Only useful for non-negative functions, as a tolerance for reaching zero. Defaults to .Machine$double.eps.

optim.maxattempts

maximum number of attempts taken trying to achieve successful convergence in optim. Defaults to 10.

optim.maxit

maximum number of iterations used by optim. Defaults to 500.

optim.method

method used by optim for minimization of the objective function. See ?optim for references. Defaults to "Nelder-Mead".

the default method is an implementation of that of Nelder and Mead (1965), that uses only function values and is robust but relatively slow. It will work reasonably well for non-differentiable functions.

method "BFGS" is quasi-Newton method (also known as a variable metric algorithm), specifically that published simultaneously in 1970 by Broyden, Fletcher, Goldfarb and Shanno. This uses function values and gradients to build up a picture of the surface to be optimized.

method "CG" is a conjugate gradients method based on that by Fletcher and Reeves (1964) (but with the option of Polak-Ribiere or Beale-Sorenson updates). Conjugate gradient methods will generally be more fragile than the BFGS method, but as they do not store a matrix they may be successful in much larger optimization problems.

optim.reltol

relative convergence tolerance used by optim. The algorithm stops if it is unable to reduce the value by a factor of 'reltol * (abs(val) + reltol)' at a step. Defaults to sqrt(.Machine$double.eps), typically about 1e-8.

p

the order of the local polynomial regression (defaults to p=1, i.e. local linear). This legacy argument is retained. If it is supplied together with modern smoothing controls, the specifications must be equivalent.

regtype

an optional modern fixed-degree spelling: "lc", "ll", or "lp". Omitting both regtype and degree preserves the historical p=1 default exactly.

degree

one nonnegative integer supplied with regtype="lp". The scalar controls the existing local-polynomial p engine; vectors are not recycled across IV stages.

nomad

logical or "auto". Only FALSE is currently supported by the bespoke npregiv objective. TRUE and "auto" fail explicitly rather than being silently ignored.

Returned Weights And Smooth Residuals

These arguments control returned kernel weights, starting values, residual smoothing, and iteration stopping behavior.

penalize.iteration: a logical value indicating whether to penalize the norm by the number of iterations or not (default TRUE)
random.seed: an integer used to seed R's random number generator. This ensures replicability of the numerical search. Defaults to 42.
return.weights.phi: a logical value (defaults to FALSE) indicating whether to return the weight matrix which when postmultiplied by the response $y$ delivers the instrumental regression
return.weights.phi.deriv.1: a logical value (defaults to FALSE) indicating whether to return the weight matrix which when postmultiplied by the response $y$ delivers the first partial derivative of the instrumental regression with respect to $z$. For more than one continuous coordinate, the returned component is a named list containing one weight matrix per pure coordinate partial.
return.weights.phi.deriv.2: a logical value (defaults to FALSE) indicating whether to return the weight matrix which when postmultiplied by the response $y$ delivers the second partial derivative of the instrumental regression with respect to $z$. For more than one continuous coordinate, the returned component is a named list containing one weight matrix per pure coordinate partial.
smooth.residuals: a logical value indicating whether to optimize bandwidths for the regression of $(y-\varphi(z))$ on $w$ (defaults to TRUE) or for the regression of $\varphi(z)$ on $w$ during iteration
start.from: a character string indicating whether to start from $E(Y|z)$ (default, "Eyz") or from $E(E(Y|w)|z)$ (this can be overridden by providing starting.values below)
starting.values: an optional finite numeric vector supplying the training-row structural initialization. Its derivative state is estimated from the corresponding local-polynomial interpolant using the selected smoothing contract; it is not set to zero. When weights are requested, this supplied initialization is treated as fixed with respect to $y$.
stop.on.increase: a logical value (defaults to TRUE) indicating whether to halt iteration if the stopping criterion (see below) increases over the course of one iteration (i.e. it may be above the iteration tolerance but increased)

Tikhonov Regularization Controls

These arguments control Tikhonov regularization and its bandwidth.

alpha: a numeric scalar that, if supplied, is used rather than numerically solving for alpha, when using method="Tikhonov".
alpha.iter: a positive numeric scalar that, if supplied, is used for the ordinary Tikhonov solve after the outer bandwidth/operator refinement rather than selecting a refined regularization parameter numerically.
alpha.max: maximum of search range for $\alpha$, the Tikhonov regularization parameter, when using method="Tikhonov".
alpha.min: minimum of search range for $\alpha$, the Tikhonov regularization parameter, when using method="Tikhonov".
alpha.tol: the search tolerance for optimize when solving for $\alpha$, the Tikhonov regularization parameter, when using method="Tikhonov".
bw: an object which, if provided, contains bandwidths and parameters (obtained from a previous invocation of npregiv) required to re-compute the estimator without having to re-run cross-validation and/or numerical optimization which is particularly costly in this setting (see details below for an illustration of its use). A Landweber-Fridman replay returns the recomputed stopping diagnostics and uses convergence="BANDWIDTH_REPLAY", including when norm.index=1.

Additional Arguments

Further arguments are passed to lower-level kernel-sum and estimation routines.

...: Only fixed bandwidths are currently supported. Supplying bwtype="generalized_nn" or bwtype="adaptive_nn" raises an informative error; use bwtype="fixed" (the default). Additional kernel and bandwidth-scaling arguments supplied to the private IV regression and kernel-sum stages. In particular, bandwidth.divide, ukertype, and okertype retain their ordinary lower-level meanings and are owned once per call, avoiding duplicate-argument collisions. Categorical bandwidth searches use the valid upper bound for the selected categorical kernel. For the plot method, graphical arguments control the fitted curve and base-graphics frame as applicable. Explicit xlim and ylim are honored; automatic limits include every active curve, overlay, and rug layer.

Details

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.

Local-polynomial nuisance fits may use numerical diagonal ridging. For its Seifert–Gasser motivation, see npreg. This is distinct from regularization of the instrumental-variable inverse problem and does not apply the ordinary-regression residual uncertainty correction to IV inference.

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

The formula interface uses explicit role order: y ~ z | w and y ~ z | w | x. Additive transformations and factors are retained through model.frame. Dot expansion and formula interaction operators are rejected because they do not identify a unique mixed-data kernel role; create an explicit transformed variable (for example with I(a * b)) instead. Native vector/data-frame calls remain supported unchanged.

For Landweber-Fridman regularization, state $N$ means the complete state after exactly $N$ residual-correction updates. Column $N$ of phi.mat, norm.stop[N], norm.index, the returned structural curve, derivatives, and optional weights therefore all refer to that same state. With penalize.iteration=TRUE, the stopping value is $N$ times the normalized squared conditional-residual norm evaluated at state $N$.

Multivariate derivative components are pure coordinate partials in the order of the continuous columns in the augmented structural block data.frame(z, x). Thus the first-derivative matrix has one named column per continuous endogenous or appended exogenous coordinate; it does not contain a mixed derivative.

The Tikhonov route solves the Darolles et al. operator equation with $r=E[Y|W]$ and right-hand side $T^*r$. Its default data-driven $\alpha$ selector is the fixed order-two criterion of Feve and Florens (2010), Equation (3.10); the returned curve is the ordinary Tikhonov solution. iterate.Tikhonov.num counts outer bandwidth and operator refinements. The current implementation selects alpha.iter during the first outer cycle and reuses it in later cycles; it does not implement a per-cycle alpha policy or the alternative Darolles extended- residual selector.

Tikhonov regularization requires computation of weight matrices of dimension $n\times n$ which can be computationally costly in terms of memory requirements and may be unsuitable for large datasets. Landweber-Fridman will be preferred in such settings as it does not require construction and storage of these weight matrices while it also avoids the need for numerical optimization methods to determine $\alpha$.

method="Landweber-Fridman" uses an optimal stopping rule based upon $||E(y|w)-E(\varphi_k(z,x)|w)||^2 $. However, if local rather than global optima are encountered the resulting estimates can be overly noisy. To best guard against this eventuality set nmulti to a larger number than the default nmulti=min(2,p) for the first iteration, where p is the dimension of the current smoothing problem.

Note that for subsequent Landweber-Fridman iterations, a “warm start” strategy is employed. The optimal bandwidths from the previous iteration are used as starting values for the current iteration. The user-supplied nmulti is respected for all iterations. For iterations after the first successful one, these optimal bandwidths serve as the first of the multiple initial points (a warm start), while any remaining restarts are cold starts. If nmulti is not explicitly supplied by the user, it defaults to min(2,p) for the first iteration and to 1 for all subsequent iterations. This strategy provides a balance between computational efficiency and robustness, allowing the numerical optimizer to refine the structural bandwidths as the residuals evolve incrementally while still guarding against local optima.

When using method="Landweber-Fridman", iteration will terminate when either the change in the value of $||(E(y|w)-E(\varphi_k(z,x)|w))/E(y|w)||^2 $ from iteration to iteration is less than iterate.diff.tol or we hit iterate.max or $||(E(y|w)-E(\varphi_k(z,x)|w))/E(y|w)||^2 $ stops falling in value and starts rising.

The option bw= would be useful, say, when bootstrapping is necessary. Note that when passing bw, it must be obtained from a previous invocation of npregiv. For instance, if model.iv was obtained from an invocation of npregiv with method="Landweber-Fridman", then the following needs to be fed to the subsequent invocation of npregiv:



    model.iv <- npregiv(\dots)

    bw <- NULL
    bw$bw.E.y.w <- model.iv$bw.E.y.w
    bw$bw.E.y.z <- model.iv$bw.E.y.z
    bw$bw.resid.w <- model.iv$bw.resid.w
    bw$bw.resid.fitted.w.z <- model.iv$bw.resid.fitted.w.z
    bw$norm.index <- model.iv$norm.index

    foo <- npregiv(\dots,bw=bw)

If, on the other hand model.iv was obtained from an invocation of npregiv with method="Tikhonov", then the following needs to be fed to the subsequent invocation of npregiv:



    model.iv <- npregiv(\dots)

    bw <- NULL
    bw$alpha <- model.iv$alpha
    bw$alpha.iter <- model.iv$alpha.iter
    bw$bw.E.y.w <- model.iv$bw.E.y.w
    bw$bw.E.E.y.w.z <- model.iv$bw.E.E.y.w.z
    bw$bw.E.phi.w <- model.iv$bw.E.phi.w
    bw$bw.E.E.phi.w.z <- model.iv$bw.E.E.phi.w.z

    foo <- npregiv(\dots,bw=bw)

Or, if model.iv was obtained from an invocation of npregiv with either method="Landweber-Fridman" or method="Tikhonov", then the following would also work:



    model.iv <- npregiv(\dots)

    foo <- npregiv(\dots,bw=model.iv)

When exogenous predictors x (xeval) are passed, they are appended to both the endogenous predictors z and the instruments w as additional columns. If this is not desired, one can manually append the exogenous variables to z (or w) prior to passing z (or w), and then they will only appear among the z or w as desired.

Value

npregiv returns a npregiv object. The generic functions print, summary, plot, fitted, residuals, and gradients support objects of this type. fitted returns the selected structural function and gradients returns its stored first derivatives. Both are training-row accessors even when fit-time evaluation data were supplied; evaluation results remain in the phi.eval and phi.deriv.eval.* fields. Residuals use the aligned training structural curve. For formula fits using an exclusion-style na.action, these accessors restore omitted rows. summary returns a structured summary.npregiv object with a dedicated print method. The printed smoothing description reports the effective continuous-stage policy using the same Kernel Regression Estimator: labels as npreg: Local-Constant, Local-Linear, or, for regtype="lp", the local-polynomial degree and basis details. The retained legacy p field is not printed separately.

The plot method displays the structural level by default, uses a separate fit-time evaluation grid when phi.eval is available, and overlays the training response by default without replacing that evaluation curve. Set gradients=TRUE for the first derivative and data.rug=TRUE for the training-support rug.

npregiv returns a list with components phi, phi.mat and either alpha when method="Tikhonov" or norm.index, norm.stop and convergence when method="Landweber-Fridman", among others. The compact bws, smoothing.spec, and stage.specs components record stage bandwidth payloads and requested/effective fixed-degree smoothing state without replacing the legacy bandwidth fields.

In addition, if any of return.weights.* are invoked (*=1,2), then phi.weights and phi.deriv.*.weights return weight matrices for computing the instrumental regression and its partial derivatives. Note that these weights, post multiplied by the response vector $y$, will deliver the estimates returned in phi, phi.deriv.1, and phi.deriv.2 (the latter only being produced when p is 2 or greater). When invoked with evaluation data, similar matrices are returned but named phi.eval.weights and phi.deriv.eval.*.weights. These weights can be used for constrained estimation, among others.

When method="Landweber-Fridman" is invoked, bandwidth objects are returned in bw.E.y.w (scalar/vector), bw.E.y.z (scalar/vector), and bw.resid.w (matrix) and bw.resid.fitted.w.z, the latter matrices containing bandwidths for each iteration stored as rows. When method="Tikhonov" is invoked, bandwidth objects are returned in bw.E.y.w, bw.E.E.y.w.z, and bw.E.phi.w and bw.E.E.phi.w.z.

References

Carrasco, M. and J.P. Florens and E. Renault (2007), “Linear Inverse Problems in Structural Econometrics Estimation Based on Spectral Decomposition and Regularization,” In: James J. Heckman and Edward E. Leamer, Editor(s), Handbook of Econometrics, Elsevier, 2007, Volume 6, Part 2, Chapter 77, Pages 5633-5751

Darolles, S. and Y. Fan and J.P. Florens and E. Renault (2011), “Nonparametric instrumental regression,” Econometrica, 79, 1541-1565.

Feve, F. and J.P. Florens (2010), “The practice of non-parametric estimation by solving inverse problems: the example of transformation models,” Econometrics Journal, 13, S1-S27.

Florens, J.P. and J.S. Racine and S. Centorrino (2018), “Nonparametric instrumental derivatives,” Journal of Nonparametric Statistics, 30 (2), 368-391.

Fridman, V. M. (1956), “A method of successive approximations for Fredholm integral equations of the first kind,” Uspeskhi, Math. Nauk., 11, 233-334, in Russian.

Horowitz, J.L. (2011), “Applied nonparametric instrumental variables estimation,” Econometrica, 79, 347-394.

Landweber, L. (1951), “An iterative formula for Fredholm integral equations of the first kind,” American Journal of Mathematics, 73, 615-24.

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

Li, Q. and J.S. Racine (2004), “Cross-validated Local Linear Nonparametric Regression,” Statistica Sinica, 14, 485-512.

Author

Jeffrey S. Racine racinej@mcmaster.ca, Samuele Centorrino samuele.centorrino@univ-tlse1.fr

Note

This function should be considered to be in ‘beta test’ status until further notice.

Examples

if (FALSE) { # \dontrun{
## This illustration was made possible by Samuele Centorrino
## <samuele.centorrino@univ-tlse1.fr>

set.seed(42)
n <- 500
trim <- 0.005

## Endogenous-regressor DGP:
##
## y = phi(z) + u
## E(u | z) != 0
## z = 0.2 w + v
## E(u | w) = 0

v <- rnorm(n, mean = 0, sd = 0.27)
eps <- rnorm(n, mean = 0, sd = 0.05)
u <- -0.5 * v + eps
w <- rnorm(n)

phi <- function(z) z^2

z <- 0.2 * w + v
y <- phi(z) + u

ivdata <- data.frame(y, z, w)
ivdata <- ivdata[order(ivdata$z), ]

model.iv <- npregiv(y ~ z | w, data = ivdata, regtype = "ll")
phi.iv <- fitted(model.iv)

## Compare with the non-IV local-linear estimator of E(y | z).

model.noniv <- npreg(y ~ z, data = ivdata, regtype = "ll")
mean.noniv <- fitted(model.noniv)

## Restrict the display to the central 99 percent of the z observations.

zlim <- quantile(ivdata$z, c(trim, 1 - trim))
keep <- ivdata$z >= zlim[1L] & ivdata$z <= zlim[2L]
ylim <- range(ivdata$y[keep], phi.iv[keep], mean.noniv[keep])

plot(ivdata$z[keep], phi(ivdata$z[keep]),
     ylim = ylim,
     ylab = "Y",
     xlab = "Z",
     main = "Nonparametric Instrumental Kernel Regression",
     type = "l",
     lwd = 2)

points(ivdata$z[keep], ivdata$y[keep], cex = 0.25, col = "grey")
lines(ivdata$z[keep], phi.iv[keep], col = "blue", lwd = 2, lty = 2)
lines(ivdata$z[keep], mean.noniv[keep], col = "red", lwd = 2, lty = 4)

legend("topright",
       c(expression(varphi(z)),
         expression(paste("Nonparametric ", hat(varphi)(z))),
         "Nonparametric E(y | z)"),
       lty = c(1, 2, 4),
       col = c("black", "blue", "red"),
       lwd = c(2, 2, 2),
       bty = "n")
} # }