Fits a stochastic frontier whose frontier function is estimated by kernel regression rather than assumed linear. Two estimators are available: the two-step estimator of Fan, Li and Weersink (1996), and the local method-of-moments estimator of Simar, Van Keilegom and Zelenyuk (2017), which additionally lets both variance components vary with the covariates.
npsfm(formula, data, method = c("FLW", "SVKZ", "PSZ", "KPST", "MY", "SZ"),
dist = c("hn", "exp", "gamma", "unif"),
regtype = c("lc", "ll"), bw.sel = c("cv.ls", "cv.aic"),
bw = NULL, cost = FALSE, eff = TRUE,
maxit = 5000, tol = 1e-3, iter = 25,
rts = c("vrs", "crs", "drs", "irs"),
prior.fit = NULL, log.form = TRUE, verbose = FALSE)An object of class "npsfareg". This is deliberately not an "sfareg" object: there is no parameter vector with standard errors, so coef(), vcov() and logLik() would have nothing meaningful to return. fitted(), residuals(), nobs(), print() and summary() are provided. Components:
The estimated frontier \(\hat{m}(x)\), i.e. the kernel fit shifted up by the estimated \(E[u]\).
Matrix of estimated frontier gradients, one row per observation and one column per covariate.
The uncorrected kernel fit of \(E[y|x]\), before the mean shift.
Composed residuals measured against the corrected frontier, \(y - \hat{m}(x)\). These are negative up to noise, unlike the centered residuals of the underlying kernel regression.
The estimated \(E[u]\): a scalar for "FLW", a vector for "SVKZ".
Estimated scale parameters. Scalars under "FLW"; vectors of \(\sigma_u(x_i)\), \(\sigma_v(x_i)\) under "SVKZ".
Returned by "FLW" with dist = "hn" only: \(\lambda = \sigma_u/\sigma_v\) and \(\sigma = \sqrt{\sigma_u^2+\sigma_v^2}\).
Returned by "FLW" with dist = "exp" or "gamma": the rate parameter of the one-sided term.
Returned by "FLW" with dist = "unif": the estimated upper bound of the uniform.
Returned by "SVKZ": the gradient of \(\sigma_u(x)\), and a logical vector flagging observations whose local third moment had the wrong sign. "PSZ" also returns sigma.u.grad and sigma.v.grad, the local-linear slopes of the two log variance functions.
Returned by "PSZ": the minqa::bobyqa status code from each observation's local optimization, 0 for success. A large share of non-zero codes means maxit is too low and the local fits should not be trusted.
Returned by "MY": how many outer iterations ran, whether the tolerance was met, and the final squared change in \((\lambda,\sigma)\).
Returned by "SZ": the smooth frontier that was monotonized, the DEA efficiency scores, and the returns-to-scale assumption used.
Jondrow et al. (1982) inefficiency predictions \(E[u|\varepsilon]\) and Battese-Coelli (1988) efficiency predictions \(E[\exp(-u)|\varepsilon]\). Returned when eff = TRUE and dist is "hn" or "exp".
The bandwidth object(s) used: one for "FLW", a list of three (r1, r2, r3) for "SVKZ".
Settings and bookkeeping.
A single-part symbolic description of the frontier, y ~ x1 + x2. Unlike the package's parametric entry points npsfm() takes no | z segment: heteroskedasticity is handled nonparametrically through the covariates themselves under method = "SVKZ", and a pipe is an error rather than something silently ignored.
A data frame containing the variables named in formula.
Which estimator to use.
"FLW"Fan, Li and Weersink (1996). Estimates \(E[y|x]\) by kernel regression, then recovers the scale parameters from the residuals -- by maximizing their concentrated likelihood in \(\lambda\) when dist = "hn", and by inverting central moments otherwise. \(\sigma_u\) and \(\sigma_v\) are constants; only the frontier is smooth.
"SVKZ"Simar, Van Keilegom and Zelenyuk (2017). Runs three local-linear regressions -- \(y\) on \(x\), then the squared and cubed residuals on \(x\) -- and inverts the local moments, so \(\sigma_u(x)\) and \(\sigma_v(x)\) both vary with the covariates. No optimizer runs. Normal-half normal only.
"PSZ" (or "KPST")Park, Simar and Zelenyuk. Local maximum likelihood: at every evaluation point the frontier and both log variance components are given local-linear expansions and the kernel-weighted normal-half normal likelihood is maximized in those \(3(k+1)\) parameters. One optimization per observation, so it is much slower than the two above. \(\sigma_u(x)\) and \(\sigma_v(x)\) vary with the covariates.
"MY"Martins-Filho and Yao. Iterative local likelihood: alternates a local-linear fit of the frontier at every evaluation point, holding \((\lambda,\sigma)\) fixed, with a global update of \((\lambda,\sigma)\) from the resulting composed residuals, until the scale parameters stop moving. Cost is (observations \(\times\) iterations) optimizations. \(\sigma_u\) and \(\sigma_v\) are constants; only the frontier is local.
"SZ"Simar and Zelenyuk (2011). Not an estimator in its own right: it takes an already-estimated smooth frontier and passes its fitted values through an output-oriented DEA, imposing the monotonicity and (under rts = "vrs"/"crs") convexity that a kernel fit does not guarantee. Supply the prior fit through prior.fit, or leave it NULL to fit "SVKZ" first. The DEA is solved internally as one linear program per unit and needs the lpSolve package; it is implemented for production frontiers only.
Matching ignores case, so "nhn", "NHN" and "Nhn" are the same choice.
Distribution of the one-sided inefficiency term, for method = "FLW" only: "hn" (half normal, the default), "exp" (exponential), "gamma", or "unif" (uniform on \([0,b]\)). Every other method is derived for the half-normal case and errors for anything else.
Kernel regression type passed to np::npregbw: "lc" (local constant, the default) or "ll" (local linear). Applies to "FLW"; "SVKZ" is local linear throughout by construction.
Bandwidth selection method: "cv.ls" (least-squares cross-validation, the default) or "cv.aic". Ignored when bw is supplied.
Optional numeric vector of bandwidths, one per covariate. When supplied, cross-validation is skipped and these are used directly -- useful for undersmoothing, for sensitivity analysis, or simply to avoid paying for bandwidth selection repeatedly in a simulation.
Logical. FALSE (the default) fits a production frontier, in which inefficiency is subtracted; TRUE fits a cost frontier.
Logical. Compute observation-level inefficiency predictions (u_hat, exp_u_hat). Available for dist = "hn" and dist = "exp"; the gamma and uniform cases return moment estimates only. Defaults to TRUE.
Maximum function evaluations for each local optimization under "PSZ" (passed to minqa::bobyqa) and maximum iterations for each local fit under "MY" (passed to optim). Ignored by the other methods. Defaults to 5000. Raise it if convergence reports many non-zero codes.
Convergence tolerance for "MY": the iteration stops when the squared change in \((\lambda,\sigma)\) falls below this. Defaults to 1e-3.
Maximum number of outer iterations for "MY". Defaults to 25.
Returns-to-scale assumption for the DEA step under "SZ": "vrs" (the default), "crs", "drs" or "irs". These restrict \(\sum_j \lambda_j\) to be unrestricted, \(= 1\), \(\le 1\) and \(\ge 1\) respectively.
For "SZ": a numeric vector of already-estimated frontier values, one per observation, to be monotonized. If NULL (the default) an "SVKZ" fit is computed first and used.
For "SZ": whether the data are in logs, as is conventional in this literature. When TRUE (the default) the frontier and covariates are exponentiated before the DEA step and the result is returned to the log scale.
Logical. Report progress through the per-observation loops of "PSZ" and the outer iterations of "MY". Defaults to FALSE.
Christopher F. Parmeter and David H. Bernstein
Both estimators relax the parametric frontier of sfm while keeping the composed-error structure \(y = m(x) + v - u\). Because a kernel regression of \(y\) on \(x\) estimates \(E[y|x] = m(x) - E[u]\) rather than \(m(x)\), both proceed by fitting that conditional mean and then shifting it back up by an estimate of \(E[u]\).
Where they differ. "FLW" treats \(\sigma_u\) and \(\sigma_v\) as constants, so its correction is a single number and the fitted gradients are those of the conditional mean. "SVKZ" estimates \(\sigma_u(x)\) from the local third moment, so the correction varies across observations and the frontier gradient picks up an extra term through the chain rule.
Least squares versus local likelihood. "FLW" and "SVKZ" both begin from a least-squares kernel regression, which estimates \(E[y|x] = m(x) - E[u]\) and therefore needs the mean shift described above. "PSZ" and "MY" instead maximize the composed-error likelihood locally, in which the local intercept is \(m(x)\) directly and no shift is applied. They pay for that with one numerical optimization per observation -- for "MY", per observation per iteration -- so expect them to be one to two orders of magnitude slower than "FLW". Both are seeded from an "FLW" fit.
Which to use. On a simulated nonlinear frontier with \(\sigma_u = 0.6\), \(\sigma_v = 0.25\) (3 replications), mean absolute frontier error at \(n = 300\) was 0.070 for "MY", 0.079 for "FLW", and 0.116 for both "SVKZ" and "PSZ". "FLW" is the steadiest for its cost; "MY" is the most accurate if the run time is acceptable; "SVKZ" and "PSZ" earn their keep when \(\sigma_u\) genuinely varies with \(x\), which this design does not test.
Wrong skew. The identification of \(\sigma_u\) rests on the residuals being negatively skewed. "FLW" inverts a single sample moment and either succeeds or warns. "SVKZ" inverts a local third moment, which is far noisier, and at any point where the estimated skew has the wrong sign the implied \(\sigma_u(x)^3\) is negative; following the paper those points are floored at \(\sigma_u(x) = 0\) and their contribution to the frontier gradient is set to zero. wrong.skew records which observations these were. A large share of them means the local third moment is too noisy to be informative and the "SVKZ" fit should not be trusted -- "FLW" is much steadier at moderate sample sizes.
The np dependency. Both estimators need kernel regression and bandwidth selection from the np package, which is listed under Suggests rather than Imports because nothing else in sfa requires it. npsfm() checks for it and stops with an install instruction if it is missing. Bandwidth selection by cross-validation is the dominant cost and scales quadratically in the sample size; supply bw to skip it.
Fan, Y., Li, Q. and Weersink, A. (1996) 'Semiparametric estimation of stochastic production frontier models', Journal of Business & Economic Statistics, 14(4), pp. 460-468.
Simar, L., Van Keilegom, I. and Zelenyuk, V. (2017) 'Nonparametric least squares methods for stochastic frontier models', Journal of Productivity Analysis, 47(3), pp. 189-204.
Jondrow, J., Lovell, C.A.K., Materov, I.S. and Schmidt, P. (1982) 'On the estimation of technical inefficiency in the stochastic frontier production function model', Journal of Econometrics, 19(2-3), pp. 233-238.
Battese, G.E. and Coelli, T.J. (1988) 'Prediction of firm-level technical efficiencies with a generalized frontier production function and panel data', Journal of Econometrics, 38(3), pp. 387-399.
sfm for parametric cross-sectional frontiers, psfm for panel models, and data_gen_cs for simulating data with known true parameters.
# \donttest{
if (requireNamespace("np", quietly = TRUE)) {
set.seed(42)
n <- 150
x1 <- runif(n, 1, 4)
x2 <- runif(n, 1, 4)
m <- 1 + 0.6 * log(x1) + 0.4 * sqrt(x2) # nonlinear frontier
d <- data.frame(y = m + rnorm(n, 0, 0.25) - abs(rnorm(n, 0, 0.6)),
x1 = x1, x2 = x2)
## Fan, Li and Weersink, normal-half normal
f <- npsfm(y ~ x1 + x2, data = d, method = "FLW", dist = "hn")
f
head(fitted(f))
head(f$exp_u_hat)
## Simar, Van Keilegom and Zelenyuk: sigma_u and sigma_v vary with x
g <- npsfm(y ~ x1 + x2, data = d, method = "SVKZ")
summary(g$sigma.u)
}
# }
Run the code above in your browser using DataLab