Package 'sfa'

Stochastic Frontier Analysis

Description

Maximum likelihood estimation of stochastic frontier models, in which observed output falls short of an unobserved production frontier by a one-sided inefficiency term in addition to symmetric noise, following Aigner, Lovell and Schmidt (1977) <doi:10.1016/0304-4076(77)90052-5> and Meeusen and van den Broeck (1977) <doi:10.2307/2525757>. Cross-sectional models cover seventeen noise and inefficiency distribution pairs, with heteroskedasticity permitted in either error component. Panel models include the true fixed and random effects estimators of Greene (2005) <doi:10.1016/j.jeconom.2004.05.003> and the four-component model of Colombi, Kumbhakar, Martini and Vittadini (2014) <doi:10.1007/s11123-014-0386-y>. Further entry points address zero inefficiency, latent classes, two-tier frontiers, sample selection, nonparametric frontiers, dependence between the error components, and endogenous regressors after Amsler, Prokhorov and Schmidt (2016) <doi:10.1016/j.jeconom.2015.06.013>. Technical efficiency prediction, marginal effects, and specification tests for wrong skewness, for the assumed inefficiency distribution and for heteroskedasticity are provided.

Details

The DESCRIPTION file:

Package:	sfa
Version:	1.2.1
Date:	2026-09-07
Title:	Stochastic Frontier Analysis
Type:	Package
Authors@R:	c(person("David", "Bernstein", email = "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2267-5741")), person("Christopher", "Parmeter", role = c("aut")), person("Alexander", "Stead", role = c("aut")))
Maintainer:	David Bernstein <[email protected]>
Description:	Maximum likelihood estimation of stochastic frontier models, in which observed output falls short of an unobserved production frontier by a one-sided inefficiency term in addition to symmetric noise, following Aigner, Lovell and Schmidt (1977) <doi:10.1016/0304-4076(77)90052-5> and Meeusen and van den Broeck (1977) <doi:10.2307/2525757>. Cross-sectional models cover seventeen noise and inefficiency distribution pairs, with heteroskedasticity permitted in either error component. Panel models include the true fixed and random effects estimators of Greene (2005) <doi:10.1016/j.jeconom.2004.05.003> and the four-component model of Colombi, Kumbhakar, Martini and Vittadini (2014) <doi:10.1007/s11123-014-0386-y>. Further entry points address zero inefficiency, latent classes, two-tier frontiers, sample selection, nonparametric frontiers, dependence between the error components, and endogenous regressors after Amsler, Prokhorov and Schmidt (2016) <doi:10.1016/j.jeconom.2015.06.013>. Technical efficiency prediction, marginal effects, and specification tests for wrong skewness, for the assumed inefficiency distribution and for heteroskedasticity are provided.
Suggests:	knitr, rmarkdown, testthat, sn, np, DEA, sandwich, lmtest, texreg, pbapply, gsl
Imports:	pso, numDeriv, plm, minqa, randtoolbox, mnormt, truncnorm, mvtnorm, Formula, methods, parallel, graphics, grDevices
Depends:	R (>= 4.0.0)
License:	GPL (>= 2)
Language:	en-US
URL:	https://www.davidharrybernstein.com/software, https://github.com/davidhbernstein/sfa
BugReports:	https://github.com/davidhbernstein/sfa/issues
LazyLoad:	yes
NeedsCompilation:	no
VignetteBuilder:	knitr
Config/testthat/edition:	3
Config/testthat/parallel:	true
Config/pak/sysreqs:	make
Repository:	https://davidhbernstein.r-universe.dev
Date/Publication:	2026-10-14 22:03:28 UTC
RemoteUrl:	https://github.com/davidhbernstein/sfa
RemoteRef:	HEAD
RemoteSha:	6aaaa4cd388d0c8f78b572696571d6ce28b0e9bc
Author:	David Bernstein [aut, cre] (ORCID: <https://orcid.org/0000-0002-2267-5741>), Christopher Parmeter [aut], Alexander Stead [aut]

Index of help topics:

calibrate_c             Fixed robustness parameter by weight matching
copsfm                  Stochastic Frontier with Dependent, and
                        Optionally Skewed, Error Components
cw_test                 Chen and Wang (2012) centered-residuals moment
                        test
data_gen_cs             Generate Cross-Sectional Data for Stochastic
                        Frontier Analysis
data_gen_p              Generate Panel Data for Stochastic Frontier
                        Analysis
density_weights         Density-power weights from a robust frontier
                        fit
efficiency              Technical Efficiency Scores
efficiency_ci           Confidence Intervals for Individual
                        Inefficiency and Efficiency
endogeneity_test        Were the regressors endogenous after all?
esfm                    Stochastic frontier tolerating the wrong
                        skewness
FinnishElec             Expenditure and output of Finnish electricity
                        distribution firms
gof_test                Goodness-of-fit test for the assumed
                        inefficiency distribution
hscore                  Hyvarinen score for a Normal-Half-Normal
                        frontier fit
hscore_select           Select the robustness tuning parameter by
                        Hyvarinen score
Indian                  Paddy farming in Aurepalle, India, 1975-76 to
                        1984-85
inefficiency_test       Test whether there is any technical
                        inefficiency
influence_sfa           Which observations are driving this frontier
                        fit?
ivsfm                   Stochastic Frontier Models with Endogenous
                        Regressors
lcsfm                   Latent Class Stochastic Frontier Models
lcsfm_homogeneity       Test a latent class frontier against
                        homogeneity
marginal_effects        Marginal Effects of the Variance Determinants
                        on Inefficiency
meanefficiency          Model-Implied Mean Efficiency
moment_range            What skewness and excess kurtosis can a pair of
                        distributions produce?
npsfm                   Nonparametric Stochastic Frontier Models
panel89                 Costs of U.S. commercial banks, 1989
pcomposed               The distribution of the composed error
pcomposed_model         Composed-error distribution function for any
                        cross-sectional model
PL80_MVTN               Pitt and Lee (1981) Model III: the multivariate
                        truncated normal panel likelihood
predict.sfareg          Predictions, fitted values and residuals for
                        stochastic frontier fits
print.sfareg            Print method for "sfareg" objects
psfm                    Panel Stochastic Frontier Models
psfm_bootstrap          Parametric Bootstrap for 'psfm()' Panel
                        Stochastic Frontier Models
selsfm                  Stochastic Frontier Models with Sample
                        Selection or Technology Choice
sfa_diagnostics         Optimizer Diagnostics for Stochastic Frontier
                        Fits
sfa-package             Stochastic Frontier Analysis
sfa-wrongskew           The wrong-skewness suite: what to do when the
                        residuals lean the wrong way
sfareg-methods          Extractor methods for "sfareg" objects
sfm                     Cross-Sectional Stochastic Frontier Models
sfma                    Average over inefficiency distributions instead
                        of choosing one
simulation_se           Simulation Standard Errors for Simulated-ML
                        Fits
skewness_decomp         Where the skewness of the composed error comes
                        from
skewness_test           Formal tests for wrong skewness in the OLS
                        residuals
spec_test               Is this pair of distributions defensible?
summary.sfareg          Summary method for "sfareg" objects
TIC                     Choosing between non-nested stochastic frontier
                        specifications
tt_boundary_report      Is a collapsed two-tier noise scale the
                        maximum, or did the optimizer stop short?
ttsfm                   Two-Tier Stochastic Frontier Model
uhet_test               Test whether technical inefficiency depends on
                        firm characteristics
USUtilities             U.S. fossil-fuel steam electric utilities,
                        1986-1997
zsfm                    Zero-Inefficiency Stochastic Frontier Models

See Also

https://www.davidharrybernstein.com/software

https://github.com/davidhbernstein/sfa

Examples

## Simple application of the generalized true random effects estimator.
library(sfa)

data_trial <- data_gen_p(t=10,N=100,  rand = 100, 
                         sig_u = 1,   sig_v = 0.3, 
                         sig_r = .2,  sig_h = .4, 
                         cons  = 0.5, beta1 = 0.5,
                         beta2 = 0.5)

psfm(formula    = y_gtre ~ x1 + x2,    
     model_name = "GTRE", 
     data       = data_trial,
     individual = "name",
     PSopt      = FALSE)

---

Fixed robustness parameter by weight matching

Description

Chooses the robustness tuning parameter so that a residual a given number of scale units from the reference point retains a common share of its clean-data influence, comparably across the robust criteria. This is the transparent alternative to the data-driven selection of hscore_select.

Usage

calibrate_c(sigma_v, sigma_u, method = c("mlqe", "psi", "mdpd", "all"),
            target = 0.10, k = 3, range = c(1e-4, 0.6))

Arguments

sigma_v, sigma_u	Positive scalars at which to calibrate, typically the estimates from a preliminary fit or the values of a simulation design.
method	Robust criterion, or "all" for a named vector.
target	Share of clean-data influence the reference residual should retain. Default 0.10.
k	Reference residual in scale units. Default 3.
range	Search interval for the tuning parameter.

Details

The calibration matches the norm of the observationwise estimating-equation contribution in the intercept, $\log\sigma_v$ and $\log\sigma_u$ directions. Matching on the intercept direction alone would be vacuous: there the criteria share the density-power factor $f^{c-1}f'$, so it would force the same tuning value for all of them by construction. They differ through the Fisher-consistency correction, which enters the scale directions.

Note that $\sigma=(\sigma_v^2+\sigma_u^2)^{1/2}$ is the scale parameter of the composed density, not the standard deviation of the composed error, which is $\{\sigma_v^2+(1-2/\pi)\sigma_u^2\}^{1/2}$. The argument k is measured in the former.

The ratio is not monotone, and there can be two answers. For the Fisher-consistency-corrected criteria the influence ratio falls below the target and is then pulled back up again by the integral correction, so it crosses the target twice. At $\sigma_v=0.30$, $\sigma_u=0.60$ the $\Psi$/MDPD solutions are $c=0.2387$ and $c=0.5385$. The larger is returned: it is the calibration reported by Bernstein, Parmeter and Wright, and it lies on the far side of the dip, so the target down-weighting holds there rather than being passed through on the way. Both are in the "roots" attribute; the non-monotonicity is reported rather than resolved silently. MLqE has a single solution, $c=0.2166$.

"psi" and "mdpd" always calibrate to the same value: the two objectives differ by the factor $(1+c)$, which is constant in the parameters and cancels from the ratio being matched.

Value

A named numeric vector of calibrated tuning parameters, carrying an attribute "roots" that lists every solution found for each criterion.

References

Bernstein, D.H., Parmeter, C.F. and Wright, I.A. (2026). On Robust Estimation of the Stochastic Frontier Model. Working paper.

See Also

hscore_select, density_weights

Examples

## the calibration used in Bernstein, Parmeter and Wright
calibrate_c(sigma_v = 0.30, sigma_u = 0.60, method = "all")

---

Stochastic Frontier with Dependent, and Optionally Skewed, Error Components

Description

Fits a cross-sectional stochastic frontier in which the noise $v$ and the inefficiency $u$ may be dependent, with the dependence carried by a copula, and in which the noise itself may be skewed. Every other model in the package assumes independence and a symmetric $v$.

Usage

copsfm(formula, data,
       copula = c("gaussian", "fgm", "frank",
                  "clayton", "clayton90", "clayton180", "clayton270",
                  "gumbel", "gumbel90", "gumbel180", "gumbel270",
                  "joe", "joe90", "joe180", "joe270", "independent"),
       inefdec = TRUE,
       udist = c("hnormal", "exponential"),
       vdist = c("normal", "logistic", "glogistic"),
       n_nodes = 128, maxit.bobyqa = 10000, maxit.psoptim = 1000,
       maxit.optim = 1000, start_val = FALSE, PSopt = FALSE,
       optHessian = TRUE, Method = "L-BFGS-B", verbose = FALSE,
       keep_objective = FALSE, rand.psoptim = NULL)

Arguments

formula	Two-sided formula for the frontier, y ~ x1 + x2. A | segment is an error: variance determinants are not yet supported alongside a copula.
data	A data.frame.
copula	The dependence family. See ‘Choosing a family’ below. family parameter independence at dependence it can express "gaussian" $\rho \in (-1,1)$ 0 both signs, no tail dependence "fgm" $\theta \in [-1,1]$ 0 both signs but weak: Spearman $\rho = \theta/3$ "frank" $\theta \in \mathbb{R}$ 0 both signs, full range, no tail dependence "clayton" $\theta > 0$ $\theta \to 0$ positive, lower-tail dependence "gumbel" $\theta \ge 1$ 1 positive, upper-tail dependence "joe" $\theta \ge 1$ 1 positive, heavier upper tail than Gumbel Clayton, Gumbel and Joe carry only positive dependence, which is a real restriction: nothing rules out a negative association between noise and inefficiency, and against negatively dependent data those families can only report the independence boundary. The 90 and 270 suffixes are the rotations that reverse the sign, and 180 is the survival copula, which preserves it. So "clayton270" gives lower-tail dependence with a negative association. "independent" sets $c \equiv 1$ and estimates no dependence parameter. It is not a null option: combined with vdist it fits the skewed-noise frontier on its own, which is a different model from anything sfm offers.
inefdec	TRUE for a production frontier, FALSE for cost.
udist	Marginal distribution of the inefficiency $u$. "hnormal" is half-normal with scale $\sigma_u$; "exponential" has mean and standard deviation $\sigma_u$.
vdist	Marginal distribution of the noise $v$. "normal" is $N(0, \sigma_v^2)$. "logistic" and "glogistic" are the generalized logistic of Bonanno and Domma (2022), centred so that $E[v] = 0$ for every parameter value: $G(v) = \left(1 + e^{-z}\right)^{-\alpha_v}, \qquad z = v/\delta_v + \psi(\alpha_v) - \psi(1).$ "logistic" fixes $\alpha_v = 1$ and "glogistic" estimates it, bounded to $[0.05, 20]$. $\alpha_v < 1$ skews the noise negative and $\alpha_v > 1$ positive. The reported scale is named delta_v, not sigma_v, because it is not the standard deviation: $\mathrm{Var}(v) = \delta_v^2[\psi'(\alpha_v) + \psi'(1)]$, so a logistic with $\delta_v = 1$ has standard deviation $\pi/\sqrt{3} = 1.814$.
n_nodes	Gauss–Legendre nodes for the integral over $u$.
keep_objective	Logical. Retain the likelihood closure on the fitted object as $objective, so the score-based tools can difference it: vcov(type = "bhhh") (OPG/BHHH standard errors), influence_sfa, the sandwich methods, TIC and vuong. New in 1.2.1. FALSE by default, because a closure carries its enclosing environment, so a fit saved with one serializes the estimation data and the quadrature nodes alongside the results.
maxit.bobyqa, maxit.psoptim, maxit.optim, start_val, PSopt, optHessian, Method, verbose, rand.psoptim	As in sfm.

Details

In brief:

• The frontier slopes behave normally. The dependence parameter does not: it needs several thousand observations, and at $n = 2000$ its sign can come out wrong, taking $\sigma_u$ with it (0.57 against a truth of 1 in the design below).

• Below a crossover between $n = 400$ and $n = 1000$, plain sfm estimates $\sigma_u$ with a lower RMSE than the true copula does.

• Prefer copula = "frank" or "clayton". The other families often return the independence bound on data generated from themselves, and copsfm warns when you pick one.

• $\alpha_v$ (vdist = "glogistic") is consistent but noisy, and below about $n = 2000$ it can return its upper bound.

• Keep n_nodes at 64 or more. At 32 the estimates are biased.

Each point is measured in the subsections below.

The model

The package's taxonomy already writes the general joint density as

$f_{V,U}(v,u) = f_V(v)\,f_U(u)\,c(F_V(v), F_U(u); \rho),$

with every existing specification setting $c \equiv 1$. This function relaxes that, and lets the two marginals vary. Composing over $v = \varepsilon + Su$,

$f_\varepsilon(\varepsilon) = \int_0^\infty f_V(\varepsilon + Su)\,f_U(u)\, c\!\left(F_V(\varepsilon + Su), F_U(u)\right) du,$

a one-dimensional integral evaluated by Gauss–Legendre quadrature on $u = t/(1-t)$ rather than by simulation.

Quadrature nodes

Measured against the closed-form normal/half-normal density at independence, the largest absolute error in $\log f$ is 1.3e-2 at 32 nodes, 2.5e-5 at 64, and 8.2e-14 at 128. Sixty-four looks adequate and is not: 1e-5 per observation is 0.02 in a log-likelihood over 2000 points, which is the scale at which competing modes are compared.

Node counts, measured in paired runs. Forty datasets sent to every node count, so the differences are the quadrature and not the draw. Against the 128-node default, the mean difference in $\hat\rho$ is +0.0002 at 256 nodes ($t = 1.37$) and +0.0046 at 64 ($t = 1.47$) – neither distinguishable – but -0.0783 at 32 ($t = -2.79$, $p = 0.008$), with per-dataset differences as large as 0.52. So the quadrature has converged by 64, and 32 must not be used as a fast option: it is biased, not merely noisier. For scale, $\hat\rho$'s own sampling standard deviation is 0.385 at $n = 600$, which dwarfs every node-count effect here.

Skewed noise, and what it is for

Bonanno and Domma (2022) argue that the “wrong skewness” anomaly is a specification failure rather than a small-sample accident. The third central moment of the composed error is

$E[(\varepsilon - E\varepsilon)^3] = J\,E[(u-Eu)^3] + E[(v-Ev)^3] + 3\mathrm{cov}(u^2, v) + 3J\,\mathrm{cov}(u, v^2) - 6[E(u)+J\,E(v)]\mathrm{cov}(u,v),$

so it depends on the skew of $v$ and on the dependence between $v$ and $u$, not only on the skew of $u$. In a model where $v$ is symmetric and independent of $u$ the first term is the only one, and a positive residual skew has nowhere to go but $\sigma_u = 0$. vdist = "glogistic" gives it somewhere else to go.

Their four specifications are four (copula, vdist) pairs, all with udist = "exponential":

their label	call
(I,S)	copula = "independent", vdist = "logistic"	independence, symmetric $v$
(I,A)	copula = "independent", vdist = "glogistic"	independence, skewed $v$
(D,S)	copula = "fgm", vdist = "logistic"	dependence, symmetric $v$
(D,A)	copula = "fgm", vdist = "glogistic"	dependence, skewed $v$

Verification against the closed form, and two corrected moment formulas

The paper derives the composed density in closed form for this case, as a sum of four Gauss hypergeometric terms – Theorem 1 for a production frontier and Theorem 2 for a cost frontier. This function computes it by quadrature instead, and the two agree to about $10^{-14}$ in $\log f$ across $\alpha_v$, $\theta$ and both orientations; the check is in test-copsfm.R. Note that the paper's summary moment formulas are not all consistent with its own appendix: Eq. (9)'s dependence term is half what the appendix's $\mathrm{cov}(U,V)$ implies, and Eq. (10)'s leading term is $-2\delta_u^3$ where an exponential's third central moment is $+2\delta_u^3$. Simulation follows the appendix. The density theorems – which are what is implemented here – are correct.

The same two slips appear in Eq. (10) and Eq. (11) of the earlier production paper (Bonanno, De Giovanni and Domma 2017), and the corrected expressions have since been derived and pinned. Writing $D = \psi(2\alpha_v) - \psi(\alpha_v)$, the covariance is $\mathrm{cov}(u,v) = \theta \delta_u \delta_v D / 2$ – twice the published value – so for a production frontier

$\mathrm{Var}(\epsilon) = \delta_u^2 + \delta_v^2[\psi'(\alpha_v) + \psi'(1)] - \theta \delta_u \delta_v D$

and

$E[(\epsilon - E\epsilon)^3] = -2\delta_u^3 + \delta_v^3[\psi''(\alpha_v) - \psi''(1)] + \tfrac{3}{2}\theta \delta_u \delta_v \{\delta_u D + \delta_v[\psi'(\alpha_v) - \psi'(2\alpha_v) - D^2]\}.$

The published dependence bracket carries an extra $\psi'(\alpha_v) + \psi'(1)$ on top of the factor of two. Both corrections rest on the FGM copula's exact factorization $E[g(v)h(u)] = E[g]E[h] + \theta\, E[g(v)(1-2G_v)]\, E[h(u)(1-2F_u)]$, whose two generalized-logistic weights are $-\delta_v D$ and $\delta_v^2[\psi'(\alpha_v) - \psi'(2\alpha_v) - D^2]$; those are verified to $10^{-9}$ in test-wrongskew-suite.R. skewness_decomp computes the three pieces for any fit.

How well the noise shape is estimated

$\alpha_v$ is consistent but noisy, and below about $n = 2000$ it can return its upper bound. Measured over 12 replications at each size, data generated at $\alpha_v = 2.5$, $\sigma_u = 0.9$, $\delta_v = 0.4$ and (for FGM) $\theta = 0.7$:

copula	$n$	median $\hat\alpha_v$	IQR	on the bound	RMSE
independent	1000	2.54	[1.49, 10.21]	3/12	1.88
independent	4000	3.05	[2.38, 3.31]	0/12	1.20
independent	16000	2.77	[2.21, 2.90]	0/12	0.48
FGM	1000	3.53	[1.92, 4.61]	1/12	4.27
FGM	4000	2.78	[2.63, 3.90]	0/12	1.25
FGM	16000	2.63	[2.38, 2.95]	0/12	0.61

$\sigma_u$ and $\delta_v$ are unaffected throughout – their medians stay within 2% of the truth at every size, as does $\hat\theta$ by $n = 16000$ (0.689 against 0.7). So the cost of a poorly determined $\alpha_v$ falls on $\alpha_v$ alone.

On an individual sample the profile can be monotone all the way to the bound, and the boundary point can genuinely attain a higher likelihood than the truth. On one such sample at $n = 2000$ the profile log-likelihood rose from -2627.5 at $\alpha_v = 1.5$ to -2622.5 at $\alpha_v = 20$ (the truth attaining -2626.3), with $\hat\theta$ falling from 0.73 to 0.38 along the way – $\alpha_v$ trades off against the dependence parameter and the intercept, not against the two scales. Read a boundary $\hat\alpha_v$ as “the noise is skewed, direction known, magnitude not”, and compare against vdist = "logistic" rather than reporting the number.

Choosing a family

The families differ in two things that matter: whether they admit negative dependence, and where they concentrate it. Gaussian and Frank spread association evenly and have no tail dependence; Clayton concentrates it in the lower tail and Gumbel and Joe in the upper; FGM can only express weak association at all.

Against that, weigh what the next subsection establishes: the dependence parameter is estimated very imprecisely. On a Gaussian design at $n = 600$ its sampling standard deviation is 0.385 against a truth of 0.5. Offering fifteen families does not make that better, and it makes one thing worse – with fifteen candidates, some family will attain the highest likelihood on independent data by chance. On a sample generated with independent $v$ and $u$, the fitted log-likelihoods across families span only about 0.8, and the best of them reports $\theta = 1.83$ for Gumbel. Treat a family comparison as descriptive, not as evidence of a dependence structure, and prefer a family chosen for a reason to one chosen by its likelihood.

How well the dependence parameter is estimated

The dependence parameter needs a large sample. This is the thing to know before using it. The frontier slopes behave normally, but $\rho$ does not. Most of these families do not recover their own dependence parameter, and copsfm warns when you pick one. A correct density does not imply an estimable parameter. Measured on 25 samples per family, generated from that family at $n = 400$ and refitted with it:

family	mean estimate (truth)	fits on the independence bound
"frank"	5.43 (5)	0%
"clayton"	2.24 (2)	0%
"gumbel"	1.29 (2)	36%
"joe"	1.52 (2)	40%
"clayton270"	0.23 (2)	56%
"gumbel90"	1.26 (2)	60%

On data generated from a Gumbel copula with Spearman $\rho = 0.685$ at $n = 2000$, every family – including the true one – returns the independence boundary, and their log-likelihoods differ by less than 0.04. The likelihood is close to flat in the dependence parameter. This is a property of the model rather than of the implementation: each density is checked against the second mixed partial of its own CDF, and each sampler against the family\'s theoretical Spearman correlation. The remaining rotations were not measured, and warn that they were not, rather than implying either outcome.

Prefer "frank" or "clayton", and treat an estimate from any other family as exploratory.

On a Gaussian-copula design with true $\rho = 0.6$, $\sigma_U = 1$, $\sigma_V = 0.4$:

$n$	$\hat\rho$	$\hat\sigma_U$	$\hat\sigma_V$
2000	-0.247	0.571	0.275
8000	0.528	0.921	0.379
20000	0.571	0.971	0.389

It is consistent – the estimates converge on the truth – but at $n = 2000$ the sign can come out wrong. The cause is joint identification rather than a bad optimiser: with $\beta$ and the scales held at their true values the likelihood peaks exactly at $\rho = 0.6$, so $\rho$ is well identified conditionally; it is the trade-off against $\sigma_U$, $\sigma_V$ and the intercept that defeats it in moderate samples, and the fitted point genuinely attains a higher likelihood than the truth.

Treat $\hat\rho$ as informative only in large samples, and compare against an independent fit (sfm) before reading anything into it.

And below a few hundred observations, the independent fit is often the better estimate – even when the dependence is real. Fitting the true copula gives a nearly unbiased but very noisy $\hat\sigma_U$; ignoring the copula gives a precise but biased one. RMSE of $\hat\sigma_U$ (truth 1.0) on Gaussian-copula data, 25 replications:

design	$n$	true copula	independence
$\rho = 0.5$	400	0.431	0.240
$\rho = 0.5$	1000	0.151	0.220
$\rho = -0.5$	400	0.251	0.196
$\rho = -0.5$	1000	0.175	0.203

The crossover sits between $n = 400$ and $n = 1000$. Below it the variance cost of estimating $\rho$ jointly with $\sigma_U$ exceeds the bias it removes, so plain sfm wins on RMSE; above it copsfm wins. The independence bias does not shrink with $n$ (about 22% of $\sigma_U$, since it is misspecification rather than noise), while the copula fit's RMSE falls as it should.

Which densities are offered

Only densities that could be verified are offered. Each is checked in test-copsfm.R: it integrates to 1 over the unit square, and it returns exactly 1 at the independence parameter. Frank, Clayton, Gumbel and Joe, whose densities concentrate in a corner where the integral check is loose, are also checked against the second mixed partial of their own CDF, a test a wrong formula cannot pass. The rotations are those same densities reflected.

Value

An object of class "sfareg", with copula, copula_par (NA when copula = "independent"), udist, vdist and alpha_v alongside the usual components. jlms is $E[u \mid \varepsilon]$ and exp_u_hat is $E[e^{-u} \mid \varepsilon]$, both by the same quadrature the likelihood used; efficiency is exp(-jlms), which is what this function has always returned under that name and is not the Battese–Coelli predictor.

References

Smith, M. D. (2008). Stochastic frontier models with dependent error components. The Econometrics Journal, 11(1), 172–192.

Bonanno, G. and Domma, F. (2022). Analytical derivations of new specifications for stochastic frontiers with applications. Mathematics, 10(20), 3876.

See Also

sfm for the independent case.

Examples

set.seed(4)
## n_nodes = 64 rather than the default 128: the quadrature is converged by 64
## (see Details), and this is the slowest entry point in the package.
n   <- 600
x1  <- rnorm(n); x2 <- rnorm(n)
z1  <- rnorm(n); z2 <- 0.6 * z1 + sqrt(1 - 0.6^2) * rnorm(n)
y   <- 0.5 + 0.8 * x1 - 0.4 * x2 + 0.4 * z1 - qnorm((1 + pnorm(z2)) / 2)
dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit <- copsfm(y ~ x1 + x2, data = dat, n_nodes = 64)
fit$out

---

Chen and Wang (2012) centered-residuals moment test

Description

Tests the distributional assumption on the composed error using the characteristic function of the centred residuals, without ever evaluating the composed-error density or distribution function.

Usage

cw_test(x, model_name = NULL, tau = 1, type = c("cosine", "sine"),
        data = NULL, n_nodes = 256L)

Arguments

x	An "sfareg" fit, or a numeric vector of regression residuals – ordinary least squares residuals are enough.
model_name	The distributional pair under the null. Required when x is a numeric vector; taken from the fit otherwise. One of "NHN", "NE", "NR", "NU", "NGE".
tau	Frequency at which the characteristic function is compared. Chen and Wang suggest a value around 1, and a single frequency: combining several oversizes the test (see Details).
type	"cosine" (recommended) or "sine".
data	The data the model was fitted to, when x is a fit.
n_nodes	Quadrature nodes for the theoretical moments.

Details

Why this test rather than gof_test. It never touches the composed-error PDF or CDF. It compares the empirical characteristic function of the centred residuals with its theoretical value, and that follows from the characteristic functions of the two components separately, which are in closed form for many pairs whose convolution is not – the normal-gamma model being the standard example of a pair with no closed-form composed density.

Centring is what makes it work. A frontier's intercept is not identified separately from $E[u]$, so the least-squares intercept estimates $\alpha - E[u]$ rather than $\alpha$. Centring the residuals cancels that inconsistent intercept outright, which is why the test is valid off an ordinary regression and does not require maximum likelihood. The variance then has to correct for two nuisances – estimating the scale parameters, and the centring itself. Omitting the second understates the variance and oversizes the test.

Use one frequency, not several. Measured here over 500 replications at a nominal 5%, half-normal data:

test	$n$	$\tau=1$	$\tau=1.5$	$\tau=2$
cosine	500	0.070	0.054	0.050
cosine	2000	0.046	0.044	0.050
sine	500	0.092	0.064	0.042
sine	2000	0.068	0.062	0.064

Against exponential inefficiency the cosine test's power at $\tau = 1$ is 0.231 at $n = 500$ and 0.952 at $n = 2000$, against 0.108 and 0.706 for the sine test. Both of Chen and Wang's conclusions reproduce: the cosine test is the less sensitive to the choice of $\tau$, and the more powerful. Passing several frequencies at once gives rejection rates of 0.09 to 0.15 at a nominal 5%, because the cosine moments at nearby frequencies are close to collinear and the covariance matrix is then near-singular; the authors report the same for their combined $\tau$.

The moment estimator solves the second- and third-moment equations, which for "NHN" and "NE" is the closed form of the paper's equation (32) and is the same inversion sfm(estimator = "cols") uses. A wrongly signed third moment leaves it with no admissible solution, and that is an error rather than a number: see skewness_test.

Value

An object of class "cw_test": the statistic D, its degrees of freedom and p-value, the moment estimates of $\sigma_v$ and $\sigma_u$, and the frequencies used.

References

Chen, Y.-T. and Wang, H.-J. (2012). Centered-residuals-based moment estimator and test for stochastic frontier models. Econometric Reviews 31(6), 625–653.

See Also

gof_test, spec_test, skewness_test

Examples

set.seed(4)
n <- 500
X <- matrix(rnorm(n * 3), n, 3)
y <- X %*% c(1, 1, 1) + rnorm(n, 0, 1) - abs(rnorm(n, 0, 1))
e <- residuals(lm(y ~ X))

cw_test(e, "NHN")            # true pair
cw_test(e, "NE")             # wrong inefficiency distribution

---

Generate Cross-Sectional Data for Stochastic Frontier Analysis

Description

data_gen_cs generates simulated cross-sectional data based on the stochastic frontier model, allowing for different distributional assumptions for the one-sided technical inefficiency error term ($u$) and the two-sided idiosyncratic error term ($v$). The model has the general form: $Y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + v - u$ where $u \geq 0$ and represents inefficiency. All variants are produced so that the user can select those that they want.

Usage

data_gen_cs(N, rand, sig_u, sig_v, cons, beta1, beta2, a, mu, sig_w = sig_u,
            shape_g = 2, m_nak = 1, mu_ln = -0.5, k_w = 1.5, lam_tsl = 1.5)

Arguments

N	A single integer specifying the number of observations (cross-sectional units).
rand	A single integer to set the seed for the random number generator, ensuring reproducibility.
sig_u	The standard deviation parameter ($\sigma_u$) for the base distribution of the one-sided error term $u$.
sig_v	The standard deviation parameter ($\sigma_v$) for the base distribution of the two-sided error term $v$.
cons	The value of the constant term (intercept) in the model.
beta1	The coefficient for the $x_1$ variable.
beta2	The coefficient for the $x_2$ variable.
a	The degrees of freedom parameter for the t half-t distribution (u_t and v_t, respectively). Requires the rt function.
mu	The mean parameter ($\mu$) for the normal truncated normal distribution (u_tn). Requires the rtruncnorm function.
sig_w	The standard deviation/scale parameter ($\sigma_w$) for the second one-sided error component $w$, used to generate two-tier model columns (see ttsfm). Defaults to sig_u.
shape_g	Shape of the gamma inefficiency draw behind y_pcs_g, which targets sfm(model_name = "NG"). The scale is set to sig_u/shape_g, so $E[u] = \sigma_u$ regardless of the shape chosen. sfm()'s NG likelihood reports the shape as mu and the scale as sigu, so the values to recover are shape_g and sig_u/shape_g.
m_nak	Nakagami shape $m$ behind y_pcs_nak, which targets sfm(model_name = "NNAK"). The spread is $\Omega = \sigma_u^2$, so $\sigma_u$ is the root-mean-square of $u$. The default 1 is the Rayleigh case; 0.5 would collapse the Nakagami onto the half-normal and so duplicate y_pcs.
mu_ln	Log-scale mean of the lognormal inefficiency draw behind y_pcs_ln, which targets sfm(model_name = "NLN"). The log-scale standard deviation is sig_u. The default -0.5 with sig_u = 1 gives $E[u] = 1$.
k_w	Weibull shape behind y_pcs_wb, which targets sfm(model_name = "NW"). The scale is sig_u.
lam_tsl	Skew parameter $\lambda$ behind y_pcs_tsl, which targets sfm(model_name = "TSL"). The scale is sig_u; $\lambda \to 0$ recovers the exponential case.

Details

The function simulates two explanatory variables, $x_1$ and $x_2$, as transformations of uniform random variables.

The function generates several different frontier models by combining various distributions for $u$ and $v$:

• **$u$ Distributions (Inefficiency):** Half-Normal (HN), Truncated Normal (TN), Half-T (HT), Half-Cauchy (HC), Exponential (E), Half-Uniform (HU).

• **$v$ Distributions (Idiosyncratic):** Normal (N), t, Cauchy (C).

**Specific Model Outputs (y_pcs variants):**

• y_pcs: Normal-Half-Normal (N-HN): $v \sim N(0, \sigma_v^2)$, $u \sim |N(0, \sigma_u^2)|$.

• y_pcs_z: N-HN with Heteroskedastic $\sigma_u$: $\sigma_{u,i} = \exp(0.9 + 0.6 Z_i)$, where $Z$ is a uniform variable.

• y_pcs_t: T-Half T (T-HT): $v \sim T(\text{df}=a) \cdot \sigma_v$, $u \sim |T(\text{df}=a)| \cdot \sigma_u$.

• y_pcs_st: the column that actually matches sfm's THT. y_pcs_t above draws its two components as two independent rt() variates, which shares the degrees of freedom but not the mixing variable, so the composed error is not skew-t and THT cannot recover it. y_pcs_st uses the single common $\lambda \sim \Gamma(a/2, a/2)$ of Tancredi (2002), i.e. $y = h + (\epsilon - z)/\sqrt{\lambda}$. y_pcs_t is retained unchanged only for backward compatibility.

• y_pcs_thn: Student's t–half-normal, matching sfm's tHN: $v \sim T(\text{df}=a) \cdot \sigma_v$ and $u \sim |N(0,\sigma_u^2)|$, drawn independently. Heavy-tailed noise with a conventional one-sided term – distinct from both y_pcs_t and y_pcs_st, in which the inefficiency is heavy-tailed too.

• y_pcs_tn: Normal-Truncated Normal (N-TN): $v \sim N(0, \sigma_v^2)$, $u \sim TN(\mu, \sigma_u^2)$ on $[0, \infty)$.

• y_pcs_e: Normal-Exponential (N-E): $v \sim N(0, \sigma_v^2)$, $u \sim Exp(\phi)$, where $\phi = 1/\sigma_u$.

• y_pcs_c: Cauchy-Half Cauchy (C-HC): $v \sim Cauchy(0, \sigma_v)$, $u \sim |Cauchy(0, \sigma_u)|$.

• y_pcs_u: Normal-Half Uniform (N-HU): $v \sim N(0, \sigma_v^2)$, $u \sim U(0, \sigma_u)$.

• y_pcs_w: Normal + Cauchy - Half-Normal: $v \sim N(0, \sigma_v^2) + Cauchy(0, \sigma_v)$, $u \sim |N(0, \sigma_u^2)|$. This introduces a composite $v$ term.

**Note:** The rtruncnorm function is required for y_pcs_tn and loads with the package. In isolation it could be loaded by using library(truncnorm).

Value

A data frame containing $N$ observations with the following columns:

name	Individual identifier (simply $1$ to $N$).
cons	The constant term value.
x1	Simulated explanatory variable $x_1$.
x2	Simulated explanatory variable $x_2$.
u, uz, u_t, u_c, u_e, u_u, u_tn, u_st, u_thn	The simulated one-sided error terms under different distributions.
u_r	Rayleigh one-sided error term, for y_pcs_r. Drawn on the convention sfm's "NR" uses, $E[u^2] = \sigma_u^2$, so the Rayleigh scale is $\sigma_u/\sqrt{2}$, $E[u] = \sigma_u\sqrt{\pi}/2$ and $Var(u) = (1-\pi/4)\sigma_u^2$.
v, v_t, v_c, v_st, v_thn	The simulated two-sided error terms under different distributions.
lam_st	The common Gamma(a/2, a/2) scale-mixing variable behind y_pcs_st; both of that column's components are divided by its square root, which is what makes the composed error skew-t.
y_pcs, y_pcs_t, y_pcs_st, y_pcs_thn, y_pcs_e, y_pcs_ez, y_pcs_c, y_pcs_u, y_pcs_z, y_pcs_w, y_pcs_tn	The dependent variable $Y$ under the corresponding SFA model distributions.
y_pcs_r	Normal-Rayleigh dependent variable, for sfm's "NR". Use this column, not y_pcs: the normal-Rayleigh and normal-half-normal composed errors are different families, not reparameterizations of one another, and "NR" cannot recover y_pcs's parameters. True values are $(\sigma_v, \sigma_u) =$ (sig_v, sig_u).
z	The auxiliary variable used for heteroskedasticity in y_pcs_z, y_pcs_ez, y_tthn_z, and y_zisf_z.
uz_e	Exponential one-sided error term with heteroskedastic scale, for y_pcs_ez.
w_tt, w_tt_hn, wz_hn	The second one-sided error component ($w$) for the two-tier columns below, under exponential, homoskedastic half-normal, and heteroskedastic half-normal distributions respectively.
zp	A second auxiliary variable used for heteroskedasticity of $\sigma_w$ in y_tthn_z.
y_ttne	Homoskedastic two-tier dependent variable ($v + w - u$, both one-sided components exponential), for ttsfm's "TTNE"/"TTNLS".
y_tthn	Homoskedastic two-tier dependent variable (both one-sided components half-normal), for ttsfm's "TTHN" with no pipes.
y_tthn_z	Heteroskedastic two-tier dependent variable ($\sigma_u$ a function of z, $\sigma_w$ a function of zp), for ttsfm's "TTHN" with formula ~x1+x2|z|zp.
eff_ind, eff_ind_z	Indicator (1 = efficient, $u$ forced to 0) used to build the zero-inefficiency columns below.
prob_z_true	The true heteroskedastic "efficient regime" probability used to draw eff_ind_z.
y_zisf	Zero-inefficiency dependent variable with constant efficient-regime probability, for zsfm's "ZISF".
y_zisf_z	Zero-inefficiency dependent variable with z-dependent efficient-regime probability, for zsfm's "ZISF_Z".
u_g, y_pcs_g	Gamma one-sided error term and the corresponding dependent variable, for sfm's "NG". Shape shape_g, scale sig_u/shape_g, so $E[u] = \sigma_u$ whatever shape is chosen.
u_nak, y_pcs_nak	Nakagami one-sided error term and dependent variable, for sfm's "NNAK". Shape m_nak, spread $\Omega = \sigma_u^2$. The default m_nak = 1 is the Rayleigh case, so this column coincides in distribution with y_pcs_r at the default.
u_ge, y_pcs_ge	Generalized-exponential one-sided error term and dependent variable, for sfm's "NGE", drawn as $u = -\sigma_u\log(1 - \sqrt{U})$ with $U \sim \mathrm{Unif}(0,1)$.
u_ln, y_pcs_ln	Lognormal one-sided error term and dependent variable, for sfm's "NLN". Log-scale mean mu_ln and log-scale standard deviation sig_u.
u_w, y_pcs_wb	Weibull one-sided error term and dependent variable, for sfm's "NW". Shape k_w, scale sig_u.
con	A constant column set to 1, potentially for use in estimation.

Author(s)

David H. Bernstein

See Also

data_gen_p for the panel generator, and sfm, zsfm and ttsfm for the estimators these columns are built for.

rnorm, runif, rt, rexp, rcauchy, rtruncnorm (if available).

Examples

# Generate 100 observations of SFA data
data_sfa <- data_gen_cs(
  N     = 100,
  rand  = 123,
  sig_u = 0.5,
  sig_v = 0.2,
  cons  = 5,
  beta1 = 1.5,
  beta2 = 2.0,
  a     = 5,   # degrees of freedom for T/Half-T
  mu    = 0.1  # mean for Truncated Normal
)

# Display the first few rows of the generated data
head(data_sfa)

# Example of a Normal-Half-Normal SFA model data
summary(data_sfa$y_pcs)
plot(density(data_sfa$y_pcs))

---

Generate Panel Data for Stochastic Frontier Analysis

Description

data_gen_p generates simulated panel data for estimating various panel stochastic frontier models, including the Generalized True Random Effects (GTRE), True Random Effects (TRE), Pooled Cross-Section (PCS), and True Fixed Effects (TFE) models. The function returns the data as a pdata.frame. All variants are produced so that the user can select those that they want.

Usage

data_gen_p(t, N, rand, sig_u, sig_v, sig_r, sig_h, cons, tau = 0.5,
                  mu = 0, beta1, beta2, eta = 0.1, b_k90 = 0.05, c_k90 = 0.01,
                  d_k90 = 0.05, e_k90 = -0.005, rho_mvtn = 0.5,
                  delta_fd = mu, fd_draw = c("folded", "truncated"))

Arguments

t	The number of time periods.
N	The number of individuals.
rand	A seed for the random number generator to ensure reproducibility.
sig_u	The standard deviation ($\sigma_u$) for the one-sided error component ($u_{it}$).
sig_v	The standard deviation ($\sigma_v$) for the two-sided error component ($v_{it}$).
sig_r	The standard deviation ($\sigma_r$) for the two-sided individual effect ($r_i$).
sig_h	The standard deviation ($\sigma_h$) for the one-sided individual effect ($h_i$).
cons	The constant term ($\beta_0$) for the frontier models.
tau	The dependence parameter ($\tau$) used for the y_tfe (TFE) model formulation, default is 0.5. See Chen, Schmidt, and Wang (2014, Journal of Econometrics).
mu	The pre-truncation mean ($\mu$) of the firm-level inefficiency $u_i^*$ behind y_fd. Default 0. See Wang and Ho (2010, Journal of Econometrics).
delta_fd	The coefficient $\delta$ in the scaling function $u_{it} = \exp(\delta z_{it}) u_i^*$ behind y_fd. Defaults to mu, which reproduces this function's historical output exactly; set it separately for a design in which $\mu$ and $\delta$ differ.
fd_draw	How $u_i^*$ is drawn for y_fd. "folded" (the default, and the historical behaviour) takes $|N(\mu, \sigma_u^2)|$; "truncated" takes $N^+(\mu, \sigma_u^2)$, which is what Wang and Ho (2010) and psfm(model_name = "FD") assume. The two coincide exactly at mu = 0 and differ otherwise. Because the truncated draw consumes random numbers differently, it also changes the columns generated after it (z_gtre, u_gtre, u_inv and those built from them), so do not mix columns from the two modes.
beta1	The coefficient for the x1 variable ($\beta_1$).
beta2	The coefficient for the x2 variable ($\beta_2$).
eta	The decay-rate parameter ($\eta$) for the time-varying inefficiency component u_bc92, following Battese and Coelli (1992): $u_{it} = u_i \exp(-\eta(T_i - t))$. Default is 0.1. See Battese and Coelli (1992, Journal of Productivity Analysis).
b_k90, c_k90	Coefficients of the Kumbhakar (1990) inefficiency path $B_{it} = (1 + \exp(bt + ct^2))^{-1}$ behind y_k1990, which targets psfm(model_name = "K1990"). The defaults give a path running from 0.485 to 0.182 over ten periods. Note that $b$ and $c$ are weakly identified: quite different pairs produce near-identical fitted paths.
rho_mvtn	Within-firm correlation of inefficiency across periods for the y_pl_mvtn column, which targets psfm(model_name = "PL80_MVTN"). The inefficiency vector is drawn from $N(0, \sigma_u^2[(1-\rho)I + \rho \mathbf{1}\mathbf{1}'])$ truncated to the negative orthant, so 0 gives inefficiency independent over time and values near 1 approach the time-invariant "PL80" case. Default 0.5. Note the realized correlation of the drawn column is ATTENUATED relative to rho by the truncation (about 0.29 at rho = 0.5, $T = 4$) – rho parameterizes the underlying normal, not the truncated draw.
d_k90, e_k90	Coefficients of the modified Kumbhakar path $B_{it} = 1 + d(t - T_i) + e(t - T_i)^2$ behind y_k1990m, which targets psfm(model_name = "K1990modified"). The defaults give a path running from 0.145 to 1 over ten periods.

Details

A pdata.frame object with $N \times t$ observations, containing the following columns:

• name Individual identifier.

• year Time period identifier.

• cons The constant term used in the data generation.

• x1, x2 Explanatory variables generated from a log-uniform distribution.

• x1_w, x2_w Explanatory variables with dependence parameter $\tau$ and linkage with $r_i$, used for the TFE model.

• u, v, r, h The generated error and individual effect components.

• y_gtre, y_tre, y_pcs, y_tfe Output variables for the Production Frontier models, including the constant.

• y_gtre_nc, y_tre_nc, y_pcs_nc Output variables for the Production Frontier models, excluding the constant.

• c_gtre, c_tre, c_pcs, c_tfe Output variables for the Cost Frontier models, including the constant.

• c_gtre_nc, c_tre_nc, c_pcs_nc Output variables for the Cost Frontier models, excluding the constant.

• y_fd Output variable for the first difference model (see Wang and Ho, 2010).

• x_fd Explanatory variable for the y_fd model.

• u_fd_star, z_fd, r_fd, u_fd Components used to generate y_fd.

• u_gtre, z_gtre, zp_gtre, y_gtre_z,y_gtre_zz,y_tre_z Variables for models with heteroskedastic inefficiency ($\sigma_{u,i} = \exp(0.33 + 0.61 Z_{i}))$. Variables for models with heteroskedastic persistent inefficiency ($\sigma_{h,i} = \exp(0.25 + 0.21 Z_{i}))$.

• u_inv, y_ssfe Time-invariant firm-level inefficiency ($u_i \sim |N(0,\sigma_u^2)|$, constant within individual across $t$) and the corresponding output variable, for psfm's "SSFE" and "PL80".

• u_bc92, y_bc92 Time-varying inefficiency following a Battese-Coelli (1992) decay structure ($u_{it} = u_i \exp(-\eta(T_i - t))$) and the corresponding output variable, for psfm's "BC92".

• B_k90, u_k90, y_k1990 The Kumbhakar (1990) inefficiency path $B_{it} = (1 + \exp(b t + c t^2))^{-1}$ (controlled by b_k90, c_k90), the resulting inefficiency $u_{it} = B_{it} u_i$, and the corresponding output variable, for psfm's "K1990".

• B_k90m, u_k90m, y_k1990m The modified Kumbhakar path $B_{it} = 1 + d(t - T_i) + e(t - T_i)^2$ (controlled by d_k90, e_k90), the resulting inefficiency, and the corresponding output variable, for psfm's "K1990modified".

• h_z The heteroskedastic persistent inefficiency $h_i$ expanded to one value per observation (constant within individual), used to build y_gtre_zz.

The data is generated based on standard Stochastic Frontier Analysis (SFA) formulations, primarily for a **Production Frontier** where the one-sided error component $u_{it}$ is subtracted:

• y_gtre: GTRE model: $y_{it} = \beta_0 + \beta_1 x_{1,it} + \beta_2 x_{2,it} + r_i - h_i + v_{it} - u_{it}$

• y_tre: TRE model: $y_{it} = \beta_0 + \beta_1 x_{1,it} + \beta_2 x_{2,it} + r_i + v_{it} - u_{it}$

• y_pcs: PCS model: $y_{it} = \beta_0 + \beta_1 x_{1,it} + \beta_2 x_{2,it} + v_{it} - u_{it}$

• y_tfe: TFE model: $y_{it} = \beta_1 x_{1,it}^w + \beta_2 x_{2,it}^w + r_i + v_{it} - u_{it}$

• y_gtre_z: GTRE with Heteroskedastic $u_{it}$: $\sigma_{u,i} = \exp(0.33 + 0.61 Z_i)$.

• y_gtre_zz: GTRE with Heteroskedastic $u_{it}$: $\sigma_{u,i} = \exp(0.33 + 0.61 Z_i)$ and $\sigma_{h,i} = \exp(0.25 + 0.21 Z_{i})$.

For **Cost Frontier** models, the one-sided error component $u_{it}$ is added (e.g., c_gtre).

The error terms are generated as:

• $r_i \sim N(0, \sigma_r^2)$ (individual two-sided effect)

• $h_i \sim |N(0, \sigma_h^2)|$ (individual one-sided effect)

• $v_{it} \sim N(0, \sigma_v^2)$ (two-sided noise)

• $u_{it} \sim |N(0, \sigma_u^2)|$ (one-sided inefficiency)

The First-Difference estimation model (y_fd) uses a variation where $r_{i,fd} \sim U(0,1)$ and $u_{it,fd} = \exp(\delta z_{it}) u_i^*$, with $u_i^*$ time-invariant: folded normal by default, truncated normal under fd_draw = "truncated" (see fd_draw).

Value

A pdata.frame object containing $N \times t$ observations suitable for Stochastic Frontier Analysis (SFA).

Author(s)

David H. Bernstein

References

Chen, Y., Schmidt, P., & Wang, H. (2014). Consistent estimation of the fixed effects stochastic frontier model. Journal of Econometrics, 181(2), 65-76.

Filippini, M., & Greene, W. H. (2016). Persistent and transient productive inefficiency: a maximum simulated likelihood approach. Journal of Productivity Analysis, 45, 187-196.

Wang, H., & Ho, C. M. (2010). Estimating fixed-effect panel stochastic frontier models by model transformation. Journal of Econometrics, 157(2), 286-296.

See Also

data_gen_cs for the cross-sectional generator, and psfm for the estimators these columns are built for.

Examples

library(sfa) 
# Generate a dataset 
data_trial <- data_gen_p(t=10, N=100, rand = 100, 
                       sig_u = 1,  sig_v = 0.3, 
                       sig_r = .2, sig_h = .4, 
                       cons = 0.5, tau = 0.5,
                       mu= 0.5, beta1 = 0.5,
                       beta2 = 0.5)
 # See the first few rows 
 head(data_trial)

---

Density-power weights from a robust frontier fit

Description

Returns the weight each observation receives in the estimating equation of a maximum $L_q$-likelihood or minimum density-power divergence fit. The weight is the fitted density raised to the power $c$, normalised so the largest is one. It is formed as $\exp(c \log f)$ rather than $(\exp \log f)^c$, so an observation whose density underflows the natural scale still receives its own weight and keeps its place in the ordering.

Usage

density_weights(object, sigma_v = NULL, sigma_u = NULL, c = NULL,
                normalize = TRUE, log = FALSE)

Arguments

object	A fitted frontier model, or a numeric vector of composed residuals.
sigma_v, sigma_u	Scale parameters. Required when object is a numeric vector; taken from the fit otherwise.
c	Robustness tuning parameter. Taken from the fit when available.
normalize	Logical; scale so the maximum weight is one. Default TRUE.
log	Logical; return $\log w$ instead of $w$. Default FALSE. Use this to rank or threshold observations: the weights of the most heavily discounted observations span hundreds of orders of magnitude and are not all distinguishable once exponentiated.

Value

A numeric vector of weights in $(0, 1]$, one per observation, or of log weights in $(-\infty, 0]$ when log = TRUE.

What this diagnostic can and cannot find

The weight is a function of the fitted density at the observation, so it responds to observations lying far from the fitted surface. It has no purchase on an observation that is wrong in a regressor and, because the fitted surface bends toward it, ends up close to that surface. A mis-recorded input can therefore be highly influential and still receive a weight near one. This is not a defect of the implementation: the bounded influence of these estimators is bounded with respect to response contamination, conditional on a fixed or bounded design. Pair this with a leverage diagnostic, and with case deletion, before concluding that an observation is or is not driving a result.

References

Bernstein, D.H., Parmeter, C.F. and Wright, I.A. (2026). On Robust Estimation of the Stochastic Frontier Model. Working paper.

See Also

hscore_select, calibrate_c

Examples

set.seed(1)
e <- rnorm(200, 0, 0.3) - abs(rnorm(200, 0, 0.6))
w <- density_weights(e, sigma_v = 0.3, sigma_u = 0.6, c = 0.217)
summary(w)
head(order(w), 10)   # the observations the estimator discounts most

## Rank on the log scale when the fit is sharply peaked: there the smallest
## weights are far below the smallest positive double and cannot be told
## apart once exponentiated.
lw <- density_weights(e, sigma_v = 0.02, sigma_u = 0.9, c = 0.217, log = TRUE)
head(order(lw), 10)

---

Technical Efficiency Scores

Description

Extracts observation-specific technical efficiency from a fitted frontier, choosing among the three usual point predictors and between the two scales the dependent variable may be on.

Usage

efficiency(object, type = c("bc", "jlms", "mode"), logDepVar = TRUE,
           newdata = NULL)

Arguments

object	An "sfareg" fit.
type	Which predictor of $u_i$ to use. "bc" (the default) is Battese–Coelli, $E[\exp(-u_i) \mid \varepsilon_i]$; "jlms" is Jondrow et al., $\exp(-E[u_i \mid \varepsilon_i])$; "mode" uses the mode of the posterior of $u_i$.
logDepVar	TRUE (the default) if the dependent variable is on a log scale, so that efficiency is $\exp(-u_i)$; FALSE if it is on a level scale, giving $1 - u_i / f(x_i)$ for a production frontier and $f(x_i) / (f(x_i) + u_i)$ for a cost frontier (inefdec = FALSE).
newdata	Data to rebuild the fitted frontier from, needed only when logDepVar = FALSE. Supply it when the data used for the fit can no longer be recovered from the stored call.

Details

Which predictor. The three differ in what they report about the same posterior, not in how well they are estimated. $E[\exp(-u)] \ge \exp(-E[u])$ by Jensen's inequality, so "bc" is never below "jlms"; the gap widens as the posterior spreads out. "mode" is the only one that can equal exactly 1, and it does so for every observation whose posterior mean is negative – typically much of the efficient tail. That is a property of the predictor rather than a defect: where the posterior mean of $u$ is negative the single most likely value really is the boundary. It is available only for the models whose posterior is a truncated normal ("NHN", "NHN_Z", "NE", "NTN").

Which scale. The package elsewhere assumes throughout that the dependent variable is logged, which is the usual case and makes $u$ a proportional shortfall. logDepVar = FALSE instead forms $(f(x_i) - u_i)/f(x_i)$ for a production frontier, and for a cost frontier, where $y = f(x) + u + v$, cost efficiency $f(x_i)/(f(x_i) + u_i)$: minimum over actual cost, the level-scale counterpart of $\exp(-u_i)$. The orientation is read from the fit's inefdec, including when it was passed as a variable. That is a ratio of two estimated quantities and is less well behaved than the log version: it is undefined where the fitted frontier is zero and negative where $u$ exceeds it, both of which happen near zero output. A warning is issued if any score comes back negative or non-finite, since silence there would be the real failure.

Value

A numeric vector of efficiency scores, one per observation used in the fit.

See Also

efficiency_ci for interval estimates of the same quantity.

Examples

set.seed(2)
n  <- 400
x1 <- rnorm(n); x2 <- rnorm(n)
y  <- 6 + 0.5 * x1 + 0.5 * x2 + rnorm(n, 0, 0.3) - abs(rnorm(n, 0, 0.8))
d  <- data.frame(y = y, x1 = x1, x2 = x2)

fit <- sfm(y ~ x1 + x2, model_name = "NHN", data = d)

summary(efficiency(fit))                      # Battese-Coelli, log scale
summary(efficiency(fit, type = "jlms"))       # never above the above
summary(efficiency(fit, type = "mode"))       # reaches exactly 1
summary(efficiency(fit, logDepVar = FALSE, newdata = d))

---

Confidence Intervals for Individual Inefficiency and Efficiency

Description

Horrace and Schmidt (1996) intervals for the inefficiency term $u_i$ of a fitted stochastic frontier model, and the implied intervals for technical efficiency $\exp(-u_i)$.

Usage

efficiency_ci(object, level = 0.95, type = c("both", "u", "te"))

Arguments

object	an object of class "sfareg", as returned by sfm.
level	coverage of the interval, a single number strictly between 0 and 1. Defaults to 0.95.
type	which intervals to return: "both" (the default), "u" for inefficiency only, or "te" for technical efficiency only.

Details

Conditional on the fitted parameters, $u_i \mid e_i$ is normal with mean mu_star and standard deviation sigma_star, truncated below at zero – the same posterior the Jondrow et al. (1982) and Battese and Coelli (1988) point predictors average over. Inverting it at $(1-\alpha)/2$ and $1-(1-\alpha)/2$ gives the bounds in closed form, so the interval costs no estimation beyond the fit itself. The efficiency bounds follow by monotonicity, with the endpoints swapped, since $\exp(-u)$ is decreasing in $u$.

These intervals condition on the estimated parameters. They describe where $u_i$ sits given $e_i$ and a known frontier, and make no allowance for sampling error in the slopes or the variance parameters. Horrace and Schmidt are explicit about this, and it is why the intervals do not narrow as the sample grows: they measure the irreducible difficulty of splitting a single residual into noise and inefficiency, not estimation uncertainty. Read a wide interval as a warning against reading much into that unit's rank.

Available for the models whose posterior really is a truncated normal: "NHN", "NHN_Z", "NE" and "NTN". Every other model_name has a posterior of a different shape, for which these formulas do not hold; the function reports that rather than returning a misleading number.

Value

A data frame with one row per observation, carrying u_lower, u_hat, u_upper and/or te_lower, te_hat, te_upper according to type, and a "level" attribute. u_hat is the Jondrow et al. (1982) posterior mean and te_hat the Battese and Coelli (1988) score – exactly the values the fit already reports, taken from the fitted object rather than recomputed.

References

Horrace, W. C. and Schmidt, P. (1996). Confidence statements for efficiency estimates from stochastic frontier models. Journal of Productivity Analysis, 7, 257–282.

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, 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, 387–399.

See Also

sfm, sfa_diagnostics

Examples

dat <- data_gen_cs(N = 200, rand = 1, sig_u = 1, sig_v = 0.3,
                   cons = 0.5, beta1 = 0.5, beta2 = 0.5, a = 1, mu = 0.5)
fit <- sfm(y_pcs ~ x1 + x2, data = dat, model_name = "NHN")

ci <- efficiency_ci(fit, level = 0.90)
head(ci)

## How much of an efficiency ranking is real? Compare the spread of the point
## predictions against the width of a single unit's interval.
diff(range(ci$te_hat))
median(ci$te_upper - ci$te_lower)

---

Were the regressors endogenous after all?

Description

A Wald test of $H_0\!: \rho = 0$ on an ivsfm fit – the hypothesis that the noise $v$ is uncorrelated with the reduced-form errors of the endogenous regressors, so that the endogeneity correction was not needed and a plain sfm fit would have been consistent. This is the Durbin–Wu–Hausman idea carried into the stochastic frontier, and is Equation (25) of Hou, Ramalho and Roseta-Palma (2025).

Usage

endogeneity_test(object, level = 0.05)

Arguments

object	A fit from ivsfm with model_name = "IVLIML" or "IVCF". "C2SLS" is a moment estimator and carries no correction parameter, so it is refused.
level	Size used for the reported verdict.

Details

The statistic. With $\hat\rho$ the vector of correlations between $v$ and the reduced-form errors,

$W = \hat\rho' \,[\widehat{\mathrm{Var}}(\hat\rho)]^{-1}\, \hat\rho \;\sim\; \chi^2(q),$

$q$ the number of endogenous regressors. With a single endogenous regressor it is exactly the squared z-ratio already printed for rho_ in the fit.

Why this is an ordinary chi-square. $\rho$ lives in the open unit ball, so $\rho = 0$ is an interior point of the parameter space. That is the opposite of the situation inefficiency_test faces, where $\sigma_u = 0$ sits on a boundary and the null distribution is a chi-bar-square mixture. No mixture is needed here.

Standard errors for $\rho$. ivsfm parameterises the correlations as $\rho = t/\sqrt{1 + t't}$, which keeps them inside the unit ball for any real $t$. Before version 1.2.0 the standard errors of $\rho$ were reported as NA on the ground that this map has a non-diagonal Jacobian. It does, but the Jacobian is short –

$\partial \rho_i/\partial t_j = \delta_{ij}/s - t_i t_j / s^3, \qquad s = \sqrt{1 + t't},$

that is $J = (I - \rho\rho')/s$ – so $\mathrm{Var}(\hat\rho) = J\,\mathrm{Var}(\hat t)\,J'$ by the delta method. At the null $J = I$, which is what makes the test well calibrated there. ivsfm now reports these standard errors and carries the full $\mathrm{Var}(\hat\rho)$ as $vcov_rho, which is what the joint test needs.

Size and power. Measured in this package over 1,000 replications per cell on the design of Hou, Ramalho and Roseta-Palma's section 4.1.1, with $T = 20$ and $N$ varying. The $\rho = 0$ block is size; the rest is power. Monte Carlo standard errors at 1,000 replications are 0.0095, 0.0069 and 0.0031 at the three levels.

$\rho$	$n$	10%	5%	1%	$se/sd$
0.0	500	0.101	0.051	0.014	0.999
0.0	1000	0.098	0.051	0.006	1.008
0.0	2000	0.100	0.053	0.009	1.009
0.0	4000	0.092	0.043	0.004	1.045
0.1	500	0.422	0.298	0.116	1.027
0.1	1000	0.629	0.492	0.270	0.991
0.1	2000	0.860	0.787	0.578	0.983
0.1	4000	0.984	0.969	0.897	0.970
0.2	500	0.879	0.807	0.611	0.976
0.2	1000	0.986	0.976	0.902	0.951
0.2	2000	1.000	1.000	0.998	0.971
0.2	4000	1.000	1.000	1.000	0.974

The test holds its size at every sample size tried, and power rises monotonically to one at both alternatives. The last column is the diagnostic behind that: it is the mean reported standard error of $\hat\rho$ divided by the actual sampling standard deviation of $\hat\rho$ across replications. The statistic is $\hat\rho^2/\widehat{\mathrm{Var}}(\hat\rho)$, so a rejection rate away from nominal comes either from $\hat\rho$ itself or from the variance, and only the second would be a defect here. At 0.95 to 1.05 throughout, the variance estimate is sound – which is what separates this test from uhet_test, where the same ratio is 0.93 and the test is conservative in consequence.

At the paper's own design ($N = 100$, $T = 20$) the larger alternatives it reports are saturated: power is 1.000 at every level at both $\rho = 0.4$ and $\rho = 0.8$. The $\rho = 0.1$ and $0.2$ rows are included precisely because a column of ones says nothing about how power builds. Replication scripts: spec_tests/hou2025_table2.R and spec_tests/hou2025_endog_size.R.

Reading a non-rejection. Failing to reject is not evidence for exogeneity. $\rho$ is identified through the instruments, so weak instruments make $\hat\rho$ imprecise and the test powerless; the print method says so rather than letting the p-value speak alone.

Value

An object of class c("sfa_endog_test", "htest") carrying statistic (the Wald $W$), parameter (its degrees of freedom, the number of endogenous regressors), p.value, the estimated rho with its standard errors and z-ratios, and reject.

References

Hou, Z., Ramalho, J. J. S. and Roseta-Palma, C. (2025). Dealing with endogeneity in stochastic frontier models: A comparative assessment of estimators. Energy Economics 151, 108922.

Amsler, C., Prokhorov, A. and Schmidt, P. (2016). Endogeneity in stochastic frontier models. Journal of Econometrics 190(2), 280-288.

Kutlu, L. (2010). Battese-Coelli estimator with endogenous regressors. Economics Letters 109(2), 79-81.

See Also

ivsfm for the fit, inefficiency_test for the boundary-null case, spec_test and moment_range for the distributional assumptions.

Examples

set.seed(9)
n  <- 600
z  <- rnorm(n); x1 <- rnorm(n); eps <- rnorm(n)
x2 <- z + eps
v  <- 0.6 * eps + sqrt(1 - 0.36) * rnorm(n)
d  <- data.frame(y = 0.5 * x1 + 0.5 * x2 + v - abs(rnorm(n)),
                 x1 = x1, x2 = x2, z = z)

fit <- ivsfm(y ~ x1 + x2, endogenous = ~ x2, instruments = ~ z,
             data = d, model_name = "IVCF")
endogeneity_test(fit)

---