nlexpect: Compute expectation of a function of the coefficients.

Description

Use integration of the joint distribution of the coefficients to compute the expectation of some function of the coefficients. Can be used for non-linear inference tests.

Usage

nlexpect(
  est,
  fun,
  coefs,
  ...,
  tol = getOption("lfe.etol"),
  lhs = NULL,
  cv,
  istats = FALSE,
  flags = list(verbose = 0),
  max.eval = 200000L,
  method = c("hcubature", "pcubature", "cuhre", "suave", "vegas", "divonne"),
  vectorize = FALSE
)

Value

The function nlexpect computes and returns the expectation of the function fun(beta), with beta a vector of coefficients. I.e., if the coefficients beta are bootstrapped a large number of times, nlexpect(est, fun) should be equal to mean(fun(beta)).

Arguments

est: object of class "felm" or "lm", a result of a call to felm() or lm.
fun: function of coefficients to be integrated. Can also be a quoted expression.
coefs: character. Names of coefficients to test. Only needed if fun is a function.
...: other arguments passed to fun or the integration routine.
tol: numeric. Tolerance for the computed integral.
lhs: character. Name of the left hand side, if est has more than one.
cv: Covariance matrix to use in place of vcov(est)
istats: logical. Should convergence information from the integration routine be included as attributes?
flags: list. Additional flags for the underlying integration routine. Not used after the package R2Cuba disappeared.
max.eval: integer. Maximum number of integral evaluations.
method: character. A method specification usable by cubature::cubintegrate. The documentation there says that "pcubature" is good for smooth integrands of low dimensions.
vectorize: logical or numeric. Use vectorized function evaluation from package cubature. This can speed up integration significantly. If method is from the Cuba library (i.e. not pcubature or hcubature), vectorize should be specified as a numeric, a vectorization factor. The default is 128.

Details

The function nlexpect integrates the function fun(x) over the multivariate normal distribution specified by the point estimates and the covariance matrix vcov(est). This is the expectation of fun(beta) if we were to bootstrap the data (e.g. by drawing the residuals anew) and do repeated estimations.

The list of coefficients used by fun must be specified in coefs.

If the function is simple, it can be specified as a quoted expression like quote(a*b+log(abs(d))). In this case, if coefs is not specified, it will be set to the list of all the variables occurring in the expression which are also names of coefficients.

fun may return a vector of values, in which case a vector of expectations is computed, like quote(c(a*b, a^3-b)). However, if the expressions contain different variables, like quote(c(a*b, d*e)), a quite compute intensive 4-dimensional integral will be computed, compared to two cheap 2-dimensional integrals if you do them separately. There is nothing to gain from using vector-valued functions compared to multiple calls to nlexpect().

You may of course also integrate inequalities like quote(abs(x1-0.2) > 0.2) to simulate the probability from t-tests or Wald-tests. See the examples.

The function you provide will get an argument ... if it does not have one already. It will also be passed an argument .z which contains the actual coefficients in normalized coordinates, i.e. if ch is the Cholesky decomposition of the covariance matrix, and pt are the point estimates, the coefficients will be pt + ch \%*\% .z. The first argument is a vector with names corresponding to the coefficients.

If you specify vectorized=TRUE, your function will be passed a list with vectors in its first argument. The function must be able to handle a list, and must return a vector of the same length as the vectors in the list. If you pass an expression like x < y, each variable will be a vector. If your function is vector valued, it must return a matrix where each column is the values.

The tol argument specifies both the relative tolerance and the absolute tolerance. If these should not be the same, specify tol as a vector of length 2. The first value is the relative tolerance, the second is the absolute tolerance. Termination occurs when at least one of the tolerances is met.

The ... can be used for passing other arguments to the integration routine cubature::cubintegrate and the function to be integrated.

Examples

Run this code


N <- 100
x1 <- rnorm(N)
# make some correlation
x2 <- 0.1 * rnorm(N) + 0.1 * x1
y <- 0.1 * x1 + x2 + rnorm(N)
summary(est <- felm(y ~ x1 + x2))
pt1 <- coef(est)["x1"]
pt2 <- coef(est)["x2"]
# expected values of coefficients, should match the summary
# and variance, i.e. square of standard errors in the summary
nlexpect(est, quote(c(x1 = x1, x2 = x2, var = c((x1 - pt1)^2, (x2 - pt2)^2))))
# \donttest{
# the covariance matrix:
nlexpect(est, tcrossprod(as.matrix(c(x1 - pt1, x2 - pt2))))
# }
# Wald test of single variable
waldtest(est, ~x1)["p.F"]
# the same with nlexpect, i.e. probability for observing abs(x1)>abs(pt1) conditional
# on E(x1) = 0.
nlexpect(est, (x1 - pt1)^2 > pt1^2, tol = 1e-7, vectorize = TRUE)
# which is the same as
2 * nlexpect(est, x1 * sign(pt1) < 0)

# Here's a multivalued, vectorized example
nlexpect(est, rbind(a = x1 * x2 < pt1, b = x1 * x2 > 0), vectorize = TRUE, method = "divonne")
# \donttest{

# Non-linear test:

# A simple one, what's the probability that product x1*x2 is between 0 and |E(x1)|?
nlexpect(est, x1 * x2 > 0 & x1 * x2 < abs(pt1), vectorize = TRUE, method = "divonne")
# Then a more complicated one with the expected value of a polynomal in the coefficients
f <- function(x) c(poly = x[["x1"]] * (6 * x[["x1"]] - x[["x2"]]^2))
# This is the linearized test:
waldtest(est, f)["p.F"]
# In general, for a function f, the non-linear Wald test is something like
# the following:
# expected value of function
Ef <- nlexpect(est, f, coefs = c("x1", "x2"))
# point value of function
Pf <- f(c(pt1, pt2))
# similar to a Wald test, but non-linear:
nlexpect(est, function(x) (f(x) - Ef)^2 > Pf^2, c("x1", "x2"), vectorize = TRUE)
# one-sided
nlexpect(est, function(x) f(x) - Ef > abs(Pf), c("x1", "x2"), vectorize = TRUE)
# other sided
nlexpect(est, function(x) f(x) - Ef < -abs(Pf), c("x1", "x2"), vectorize = TRUE)
# }