sfa

Stochastic frontier analysis in R. A single, consistent interface to a wide range of cross-sectional, panel, zero-inefficiency, two-tier and nonparametric stochastic frontier models, with a common formula syntax for modelling the variance of each error component and a common "sfareg" result object that works with the standard R modelling generics.

Written by David H. Bernstein, Christopher F. Parmeter and Alexander D. Stead.

Installation

The released version from CRAN:

install.packages("sfa")

The development version from GitHub:

# install.packages("remotes")
remotes::install_github("davidhbernstein/sfa")

The development version here is ahead of CRAN; see NEWS.md for what has changed, including one deliberate breaking change to psfm(model_name = "TFE").

Quick start

library(sfa)

## Simulate a cross-section with known true parameters
cs <- data_gen_cs(N = 1000, rand = 1, sig_u = 0.3, sig_v = 0.3,
                  cons = 0.5, beta1 = 0.5, beta2 = 0.5, a = 4, mu = 1)

## Normal-half normal frontier
fit <- sfm(y_pcs ~ x1 + x2, model_name = "NHN", data = cs)

summary(fit)
coef(fit)              # lambda, sigma, (Intercept), x1, x2 -- see note below
logLik(fit)
head(fit$exp_u_hat)    # technical efficiency, E[exp(-u) | e]
head(fit$med_u_hat)    # median efficiency prediction (NHN only)

Which efficiency predictors come back depends on the model: exp_u_hat (Battese and Coelli 1988) is returned by most of sfm()’s models, and the Jondrow et al. (1982) point predictor u_hat = E[u | e] by NE, NTN, NU, NGE, NLN, NW, THT and tHN. See ?sfm for the full list.

A panel model, with a determinant of the inefficiency variance:

pd <- data_gen_p(t = 10, N = 100, rand = 100, sig_u = 1, sig_v = 0.3,
                 sig_r = 0.2, sig_h = 0.4, cons = 0.5, beta1 = 0.5, beta2 = 0.5)

fit_p <- psfm(y_tre_z ~ x1 + x2 | z_gtre, model_name = "TRE_Z",
              data = pd, individual = "name")

The five entry points

The first four return an object of class "sfareg". npsfm() returns "npsfareg" instead — a kernel-estimated frontier has no parameter vector with standard errors, so coef(), vcov() and logLik() would have nothing to return.

sfm() — cross-sectional

sfm() also offers estimator = "cols" — corrected OLS (Olson, Schmidt and Waldman 1980), closed-form and deterministic, for NHN, NE and NG — and robust divergence-based alternatives to MLE via robust = "mlqe" | "psi" | "mdpd" for NHN.

psfm() — panel

GTRE_SEQ1, GTRE_SEQ2 and SSFE are not maximum likelihood, so logLik() (and hence AIC()/BIC()) returns NA for them.

psfm_bootstrap() provides a parametric bootstrap for GTRE-family fits, parallelised over cores.

zsfm() — zero inefficiency

ZISF and ZISF_Z: a mixture of a fully efficient regime and an inefficient frontier regime, with the regime probability optionally parameterised by covariates (ZISF_Z).

ttsfm() — two tier

TTNE (normal–exponential–exponential), TTHN (normal–half normal–half normal), and TTNLS (nonlinear least squares, no distributional assumption beyond the means of the two one-sided components).

npsfm() — nonparametric

Estimates the frontier by kernel regression instead of assuming it linear.

f <- npsfm(y ~ x1 + x2, data = d, method = "FLW", dist = "hn")
head(fitted(f))       # the estimated frontier
head(f$exp_u_hat)     # technical efficiency

PSZ and MY run 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. npsfm() takes a single-part formula and rejects a | z segment: its heteroskedasticity is nonparametric in the covariates themselves.

Kernel regression comes from np, and SZ additionally needs Benchmarking. Both are in Suggests, not Imports, so they are only required if you actually call npsfm():

install.packages(c("np", "Benchmarking"))

Formula syntax

Variance determinants are supplied in extra pipe-delimited segments:

y ~ x1 + x2 | z | zp

• the first segment is the frontier,
• the second parameterises the variance of the first one-sided component (sigma_u),
• the third parameterises a second one-sided component where the model has one (sigma_w in two-tier models, sigma_h in GTRE_Z).

Omitted segments default to 1, i.e. homoskedastic.

The link function differs by model family. sfm()’s NHN_Z/NE_Z and ttsfm()’s TTNE/TTHN use sigma = exp(z'delta), while psfm()’s GTRE_Z/TRE_Z use sigma = sqrt(exp(z'delta)) — that is, delta parameterises the variance rather than the standard deviation. Check which convention applies before interpreting a coefficient on z.

Working with results

"sfareg" objects support the usual generics:

coef(fit); vcov(fit); logLik(fit); nobs(fit); AIC(fit); BIC(fit)
fitted(fit); residuals(fit); predict(fit, newdata = ...)
print(fit); summary(fit)

fit$out is the source of truth — a 3 x p matrix of estimates, standard errors and t-values. Its column names vary by model: several models report the lambda = sigma_u/sigma_v, sigma = sqrt(sigma_u^2 + sigma_v^2) reparameterisation rather than the raw scale parameters.

npsfm() fits are the exception. They carry no out matrix and no standard errors, so only fitted(), residuals(), nobs(), print() and summary() apply; read the frontier, its gradients and the scale estimates off the returned object ($frontier, $frontier.grad, $sigma.u, $sigma.v).

Simulating data

data_gen_cs() and data_gen_p() generate cross-sectional and panel data with known true parameters. Each returns a data frame with one response column per model family (y_pcs, y_pcs_z, y_pcs_r, y_tre_z, …), so a given model_name is matched to the column generated under its own assumptions. These generators are how the package’s estimators are checked against known truth, including npsfm()’s — NPSFM_FLW and NPSFM_SVKZ are registered in the root-n convergence framework and both pass.

Included data

Citation

citation("sfa")

Bernstein, D. H., Parmeter, C. F., and Stead, A. D. (2026). Stochastic Frontier Analysis: The sfa Package. Working Paper.

License

GPL (>= 2). See LICENSE.md.