Package {sfa}

Version:	1.2.0
Date:	2026-09-07
Title:	Stochastic Frontier Analysis
Type:	Package
Maintainer:	David Bernstein <davebernstein1@gmail.com>
Description:	Provides a user-friendly framework for estimating a wide variety of cross-sectional and panel stochastic frontier models. Suitable for a broad range of applications, the implementation offers extensive flexibility in specification and estimation techniques.
Suggests:	knitr, rmarkdown, testthat, sn, np, lpSolve, sandwich, lmtest, texreg
Imports:	pso, numDeriv, gsl, plm, minqa, randtoolbox, mnormt, truncnorm, tmvtnorm, Formula, methods, pbapply, pracma, MASS, parallel, graphics, grDevices
Depends:	R (≥ 4.0.0)
License:	GPL-2 | GPL-3 [expanded from: 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
Packaged:	2026-09-07 19:32:03 UTC; davidbernstein
Author:	David Bernstein [aut, cre], Christopher Parmeter [aut], Alexander Stead [aut]
Repository:	CRAN
Date/Publication:	2026-09-07 23:40:02 UTC

Stochastic Frontier Analysis

Description

Provides a user-friendly framework for estimating a wide variety of cross-sectional and panel stochastic frontier models. Suitable for a broad range of applications, the implementation offers extensive flexibility in specification and estimation techniques.

Details

The DESCRIPTION file:

Index of help topics:

FinnishElec             FinnishElec
Indian                  Indian
PL80_MVTN               Pitt and Lee (1981) Model III: the multivariate
                        truncated normal panel likelihood
TIC                     Choosing between non-nested stochastic frontier
                        specifications
USUtilities             USUtilities
calibrate_c             Fixed robustness parameter by weight matching
copsfm                  Stochastic Frontier with Dependence Between
                        Noise and Inefficiency
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
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
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                 Panel89
pcomposed               The distribution of the composed error
pcomposed_model         Composed-error distribution function for any
                        cross-sectional model
predict.sfareg          Predictions, fitted values and residuals for
                        stochastic frontier fits
print.sfareg            Print method for "sfareg" objects
psfm                    psfm
psfm_bootstrap          Parametric Bootstrap for 'psfm()' Panel
                        Stochastic Frontier Models
selsfm                  Stochastic Frontier Models with Sample
                        Selection or Technology Choice
sfa-package             Stochastic Frontier Analysis
sfa_diagnostics         Optimizer Diagnostics for Stochastic Frontier
                        Fits
sfareg-methods          Extractor methods for "sfareg" objects
sfm                     sfm
sfma                    Average over inefficiency distributions instead
                        of choosing one
simulation_se           Simulation Standard Errors for Simulated-ML
                        Fits
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
ttsfm                   Two-Tier Stochastic Frontier Model
uhet_test               Test whether technical inefficiency depends on
                        firm characteristics
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)
               

FinnishElec

Description

Cross-sectional data on Finnish electricity distribution firms, including annual averages of expenditure and output measures over a four-year regulatory period.

Usage

data("FinnishElec")

Format

A data frame with 89 observations on the following 6 variables.

id

a character vector containing a unique identifier for each distribution firm

x

a numeric vector containing total expenditure (TOTEX*) (1000 Euros)

y1

a numeric vector containing weighted energy transmitted through the network (GWh of 0.4 kV equivalents)

y2

a numeric vector containing total length of the network (km)

y3

a numeric vector containing total number of customers connected to the network

z

a numeric vector containing the proportion of underground cables in the total network length.

Details

*TOTEX includes capital expenditure (CAPEX), controllable operational expenditure (OPEX), and estimated external cost of interruptions.

Source

Kuosmanen, T. (2012). 'Stochastic semi-nonparametric frontier estimation of electricity distribution networks: Application of the StoNED method in the Finnish regulatory model.' Energy Economics, 34(6), pp. 2189-2199. doi:10.1016/j.eneco.2012.03.005

Examples

data(FinnishElec)
plot(FinnishElec)

Indian

Description

Panel data on 34 paddy farmers from Aurepalle, India, collected over ten years (1975-76 to 1984-85). The panel is unbalanced: farmers are observed between 2 and 10 times. Includes farmer characteristics (age, schooling) and production variables (output, land, labor, bullocks, input costs).

Usage

data("Indian")

Format

A data frame with 273 observations (an unbalanced panel of 34 farmers over 10 years) on the following 10 variables.

id

a numeric vector containing a unique identifier for each farmer

yr

a numeric vector containing the year of the observation

age

a numeric vector containing the age of the primary decision maker

school

a numeric vector containing the number of years of schooling of the primary decision maker

yvar

a numeric vector containing the natural logarithm of the total value of output (rupees)

Lland

a numeric vector containing the natural logarithm of the total area of land operated (ha)

PIland

a numeric vector containing the proportion of land that is irrigated

Llabor

a numeric vector containing the natural logarithm of the total number of hours of hired and family labour used

Lbull

a numeric vector containing the natural logarithm of the number of hours of bullock labour used

Lcost

a numeric vector containing the natural logarithm of the value of inputs including fertilizer, manure, pesticides, machinery, etc.

Source

Battese, G.E. and Coelli, T.J. (1995) 'A model for technical inefficiency effects in a stochastic frontier production function for panel data', Empirical Economics, 20(2), pp. 325-332. doi:10.1007/BF01205442.

References

Battese, G.E. and Coelli, T.J. (1992) 'Frontier production functions, technical efficiency and panel data: With application to paddy farmers in India', Journal of Productivity Analysis, 3(1-2), pp. 153-169. doi:10.1007/BF00158774.

Examples

data(Indian)

Pitt and Lee (1981) Model III: the multivariate truncated normal panel likelihood

Description

Notes on psfm(model_name = "PL80_MVTN"), which fits the likelihood derived in Appendix 2 of Pitt and Lee (1981) but never used there.

Details

The model is

y_{it} = x_{it}'\beta + u_{it} + v_{it}, \qquad u_{it} \le 0,

with the firm's inefficiency vector u_i = (u_{i1},\dots,u_{iT})' drawn from a T-variate normal N(0,\Sigma) truncated to the negative orthant, and v_{it} iid N(0,\sigma_v^2) independent of u. Unlike model_name = "PL80", which holds inefficiency fixed over time, here it varies across periods and is correlated within a firm – \Sigma carries that dependence.

Why it was never used. Pitt and Lee derived this likelihood and then set it aside, writing that it “is difficult to evaluate since the quantities P_0 and P(y_i - x_i\beta) involve T-dimensional numerical integrals”, and estimating Model III by Zellner seemingly-unrelated regression instead. Those quantities are orthant probabilities of a multivariate normal. mnormt::sadmvn() evaluates one in about 3 milliseconds at T = 6, so a likelihood evaluation costs roughly N+1 of them – about 0.3 s at N = 100, and about 12 s for a whole fit at N = 80, T = 4. What was intractable in 1981 is merely slow now.

Writing Q^{-1} = \Sigma^{-1} + I/\sigma_v^2 and \mu_i = Q\varepsilon_i/\sigma_v^2, the per-firm log density is

\log f(\varepsilon_i) = -\tfrac{T}{2}\log 2\pi - T\log\sigma_v - \tfrac{1}{2}\log|\Sigma| + \tfrac{1}{2}\log|Q| - \log P_0 - \tfrac{1}{2}\left(\tfrac{\varepsilon_i'\varepsilon_i}{\sigma_v^2} - \tfrac{\varepsilon_i'Q\varepsilon_i}{\sigma_v^4}\right) + \log P_i,

where P_0 = \Pr(w \le 0) for w \sim N(0,\Sigma) is the truncation constant and P_i = \Pr(w \le 0) for w \sim N(\mu_i, Q).

\Sigma is equicorrelated, not unrestricted. It is parameterized as \Sigma = \sigma_u^2[(1-\rho)I + \rho \mathbf{1}\mathbf{1}'], costing two parameters. An unrestricted \Sigma costs T(T+1)/2 – 21 at T=6, 55 at T=10 – on top of \beta and \sigma_v, every one identified only through orthant probabilities. The equicorrelated form captures what the general \Sigma was introduced for: dependence of a firm's inefficiency across periods. \rho = 0 gives inefficiency independent over time; \rho \to 1 approaches the time-invariant "PL80" case. Because the form is equicorrelated every matrix quantity above is closed form (Sherman-Morrison), so only the orthant probabilities are numerical.

Requirements and limits. A balanced panel with T \ge 2: \Sigma is a single T \times T matrix shared by every firm, and with T = 1 there is no cross-period dependence for it to describe. An unbalanced panel is an error pointing at "PL80". Reported parameters are sigv, sigu, rho and the frontier coefficients.

u_hat is the posterior mean of u from the Gaussian part, floored at zero, not E[u_{it}\mid\varepsilon_i]: the exact conditional mean of a truncated multivariate normal is another T-dimensional integral. It is exact where the truncation does not bind and an approximation where it does.

Author(s)

David H. Bernstein

References

Pitt, M.M. and Lee, L.-F. (1981) 'The measurement and sources of technical inefficiency in the Indonesian weaving industry', Journal of Development Economics, 9(1), pp. 43-64. See Appendix 2 for the likelihood.

See Also

psfm, data_gen_p for the y_pl_mvtn column that tests it.

Examples

library(sfa)
d <- data_gen_p(t = 4, N = 60, rand = 5, sig_u = 1, sig_v = 0.3, sig_r = 0.2,
                sig_h = 0.4, cons = 0.5, beta1 = 0.5, beta2 = 0.5)
f <- psfm(y_pl_mvtn ~ x1 + x2, model_name = "PL80_MVTN",
          data = as.data.frame(d), individual = "name")
f$out

Choosing between non-nested stochastic frontier specifications

Description

Model selection for stochastic frontier models when the candidates are not nested and none of them need be correctly specified. TIC is Takeuchi's information criterion, a version of AIC whose penalty is estimated from the data rather than assumed equal to the number of parameters. vuong is Vuong's pairwise test, which returns a p-value and is allowed to conclude that two specifications are indistinguishable.

Usage

TIC(object, detail = FALSE)

vuong(object1, object2, correction = c("none", "aic", "tic"), level = 0.05)

Arguments

Details

Why not AIC. sfa offers many model_names for the same data, and AIC assumes the fitted model is correctly specified – exactly the assumption in doubt when choosing between "NHN", "NE", "NG" and "NTN". Most such pairs are not nested either, so the ordinary likelihood ratio test has no chi-square limit.

What TIC does. Takeuchi's criterion replaces AIC's penalty d with \mathrm{tr}[H(\hat\theta) I(\hat\theta)^{-1}], where I is the sample Fisher information and H the outer product of the per-observation scores. Under correct specification the information matrix equality gives H = I, the trace collapses to d, and TIC equals AIC. The gap between penalty and df – reported as ratio under detail = TRUE – is therefore a readable diagnostic in its own right: a ratio near one says the distributional assumption is not doing visible damage.

What vuong does. It tests H_0 that the two models are equally close to the truth in the Kullback-Leibler sense. Neither model has to be correct. With m_i the difference of per-observation log-likelihoods, the statistic is n^{-1/2}\sum_i m_i / \hat\sigma, asymptotically standard normal. Because it is a test rather than a criterion it can return "neither", which is frequently the honest answer and is not available from AIC.

Requirement. Both need per-observation log-likelihoods, so the fits must be made with keep_objective = TRUE. This is supported by sfm (for "NHN", "NE", "NR", "NG", "NNAK", "THT", "NTN", "NHN_Z", "NE_Z", "NU", "NGE", "NLN", "NW", "tHN" and "TSL") and by psfm. vuong additionally requires both fits to use the same observations, which it checks by sample size and warns about when the data arguments differ.

Estimators that maximise nothing. psfm's "GTRE_SEQ1", "GTRE_SEQ2" and "SSFE", and ivsfm's "C2SLS", are moment-based or FE and have no log-likelihood, so no likelihood-based criterion applies to them and both functions refuse rather than returning a number built on NA.

Value

TIC returns a single numeric value, or with detail = TRUE a list with components TIC, AIC, logLik, df, penalty and ratio.

vuong returns an object of class "sfa_vuong" with components statistic, p.value, n, lr, lr_adjusted, omega, correction, penalty, critical, level, models, favoured and exact_null. Lower TIC is better; favoured is the name of the winning model_name, or "neither".

References

Lai, H.-p. and Huang, C.J. (2010). Likelihood ratio tests for model selection of stochastic frontier models. Journal of Productivity Analysis, 34(1), 3–13.

Takeuchi, K. (1976). Distribution of informational statistics and a criterion of model fitting. Suri-Kagaku (Mathematical Sciences), 153, 12–18.

Vuong, Q.H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica, 57(2), 307–333.

See Also

sfm, skewness_test, sfa_diagnostics, AIC

Examples

set.seed(42)
n  <- 400
x1 <- rnorm(n); x2 <- rnorm(n)
d  <- data.frame(y = 1 + 0.5 * x1 + 0.5 * x2 +
                     rnorm(n, 0, 0.4) - abs(rnorm(n)),
                 x1 = x1, x2 = x2)

a <- sfm(y ~ x1 + x2, model_name = "NHN", data = d, keep_objective = TRUE)
b <- sfm(y ~ x1 + x2, model_name = "NE",  data = d, keep_objective = TRUE)

## Lower is better, and comparable across non-nested models.
TIC(a)
TIC(b)

## Is the penalty near the parameter count? If so, TIC has little to add.
TIC(a, detail = TRUE)$ratio

## And a test, which may decline to pick a winner.
vuong(a, b)

USUtilities

Description

Panel data on U.S. investor-owned fossil fuel-fired steam electric utilities for the period 1986-1997. These data include measures of output, capital, labour and maintenance, and fuel.

Usage

data("USUtilities")

Format

A data frame with 972 observations (a balanced panel of observations on 81 utilities over 12 years) on the following 7 variables.

firmID

a numeric vector containing a unique firm identifier

year

a numeric vector containing the year of the observation

q

a numeric vector containing net steam electric power generation (MWh)

K

a numeric vector containing capital stock, calculated using a method described by Christensen and Jorgenson (1970)

L

a numeric vector containing quantity of labor and maintenance, calculated as cost divided by price index

F

a numeric vector containing quantity of fuel used, calculated as fuel costs divided by fuel price index

trend

a numeric vector containing an annual time trend (1992=100)

Details

The dataset covers 81 investor-owned utilities after aggregating subsidiaries and excluding plants in states with partial deregulation plans. Data sources include the Energy Information Administration (EIA), Federal Energy Regulatory Commission (FERC), and Bureau of Labor Statistics (BLS). Output is net steam electric generation from fossil fuel-fired boilers.

Source

Rungsuriyawiboon, S. and Stefanou, S.E. (2007). 'Dynamic Efficiency Estimation: An Application to U.S. Electric Utilities.' Journal of Business & Economic Statistics, 25(2), pp. 226-238. doi:10.1198/073500106000000288

References

Christensen, L.R. and Jorgenson, D.W. (1970). 'U.S. Real Product and Real Factor Input, 1928-1967.' Review of Income and Wealth, 16(1), pp. 19-50. doi: 10.1111/j.1475-4991.1970.tb00695.x

Examples

data(USUtilities)

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

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 Dependence Between Noise and Inefficiency

Description

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

Usage

copsfm(formula, data,
       copula = c("gaussian", "fgm", "frank",
                  "clayton", "clayton90", "clayton180", "clayton270",
                  "gumbel", "gumbel90", "gumbel180", "gumbel270",
                  "joe", "joe90", "joe180", "joe270"), inefdec = TRUE,
       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,
       rand.psoptim = NULL)

Arguments

Details

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 exactly that and nothing else: the marginals remain normal and half-normal. 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.

Why n_nodes defaults to 128. 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.

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:

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.

Choosing a family, and a warning about doing so. 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 paragraphs establish: 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.

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

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:

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.

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.

Families. Only densities that could be verified are offered. Each is checked in the tests two ways: it integrates to 1 over the unit square, and it returns exactly 1 at the independence parameter. Frank, Clayton and Gumbel are deliberately absent rather than transcribed without a source to check against.

Value

An object of class "sfareg", with copula and copula_par alongside the usual components. jlms is E[u \mid \varepsilon] computed by the same quadrature the likelihood used.

References

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

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

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:

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

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:

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)

Arguments

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} is generated using a heteroskedastic truncated-normal structure, reflecting an alternative model type.

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.

Usage

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

Arguments

Value

A numeric vector of weights in (0, 1], one per observation.

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

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

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). 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

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

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.

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)

Stochastic frontier tolerating the wrong skewness

Description

Fits the extended stochastic frontier model of Hafner, Manner and Simar (2018), in which one parameter \gamma carries both the scale of inefficiency (its size) and the direction of skewness (its sign), so a wrongly skewed sample gives a well-defined non-degenerate fit instead of collapsing.

Usage

esfm(formula, data, dist = c("halfnormal", "exponential"), inefdec = TRUE,
     start = NULL)

symmetry_test(object)

Arguments

Details

The classical frontier implies a negatively skewed composed error. When the sample skewness comes out positive – common at small n or large \sigma_v/\sigma_u, and a small-sample accident rather than evidence against the model – the MLE collapses: \sigma_u is exactly zero, every firm is fully efficient, and the fit is uninformative. That is the Type I failure skewness_test reports.

This model keeps the frontier and widens the distribution of u. For \gamma > 0, u has the classical density with scale |\gamma|; for \gamma < 0 that density is mirrored about zero, shifted right by B = a_0|\gamma| and truncated to [0, B], giving negative skewness for u and so positive skewness for the composed error. The truncation constant a_0 is fixed, not estimated: it is chosen so E[u] = k_1|\gamma| whichever sign \gamma takes, which keeps the model a one-parameter scale family and avoids the identification difficulties of bounded-inefficiency models. a_0 = 1.3892032925 for the half-normal and 1.5936242600 for the exponential.

The classical model is nested at \gamma > 0, and \gamma = 0 is an interior point rather than a boundary. So the likelihood ratio test of H_0: \gamma = 0 – no inefficiency, symmetric errors – is an ordinary \chi^2(1), not the chi-bar-square mixture that inefficiency_test needs for the classical model.

What it does on a wrongly skewed sample. On a 50-observation draw with positive residual skewness, sfm(model_name = "NHN") returns \sigma_u = 0.0003 and mean efficiency 0.99976 – every firm on the frontier. esfm() returns \hat\gamma = -0.57 with slopes still near their true values.

Size of the LR test, 2000 replications against the paper's 100,000, on its design (\sigma_v = 1, \gamma = 0):

The test is badly oversized below n = 100 and usable from n = 200, which is the paper's own conclusion.

Efficiency bias, and a caveat. The paper reports smaller bias in mean efficiency than the classical model even when the population skewness has the correct sign. That reproduces when inefficiency is not small: at \gamma = 0.5, \sigma_v = 0.25 the bias in E[\exp(-U)] is -0.0001 at n = 100 against the classical model's +0.021. But at \gamma = 0.3, where 22% of samples are wrongly skewed, this model overshoots the other way: bias -0.040 against the classical +0.036, and it is the larger of the two at n = 50 (-0.080 against +0.044). A small true \gamma lets the extended model fit sizeable inefficiency of either sign to noise. It is the better choice when the wrong skewness is the problem, not a free improvement everywhere.

Value

An object of class "esfm" (and "sfareg"), whose out matrix reports gamma, sigv and the frontier coefficients, with Efficiency holding E[\exp(-U) \mid W] for each firm.

References

Hafner, C. M., Manner, H. and Simar, L. (2018). The “wrong skewness” problem in stochastic frontier models: A new approach. Econometric Reviews 37(4), 380–400.

Almanidis, P., Qian, J. and Sickles, R. C. (2014). Stochastic frontier models with bounded inefficiency. In Festschrift in Honor of Peter Schmidt.

See Also

skewness_test, inefficiency_test, sfm

Examples

set.seed(3)
n <- 200
lx1 <- rnorm(n, 1.5, 0.3); lx2 <- rnorm(n, 1.8, 0.3)
y <- 0.9 + 0.6 * lx1 + 0.5 * lx2 + rnorm(n, 0, 0.25) - abs(rnorm(n, 0, 0.5))
d <- data.frame(y, lx1, lx2)

f <- esfm(y ~ lx1 + lx2, data = d)
f
symmetry_test(f)

Goodness-of-fit test for the assumed inefficiency distribution

Description

Tests the distributional assumption on u, holding the normality of v fixed, following Wang, Amsler and Schmidt (2011).

Usage

gof_test(object, data = NULL, test = c("ks", "chisq"),
         null = c("bootstrap", "asymptotic"), B = 199, cells = 10,
         seed = NULL)

Arguments

Details

The distribution of u is the assumption in this model least often defended and least often tested, and it is testable: with the noise distribution maintained, the assumed u implies a distribution for the composed error, so rejecting that is rejecting the assumed u.

The test is on \varepsilon, not on \hat u. Wang, Amsler and Schmidt point out that \hat u = E[u \mid \varepsilon] is a monotonic function of \varepsilon, so the KS test is identical either way and the chi-square test is identical when the cells are defined conformably – but \varepsilon is much easier to work with. It is a mistake, not a diagnostic, to compare the observed spread of \hat u with the assumed density of u: those are different distributions.

Parameter estimation cannot be ignored. Both statistics are evaluated at \hat\theta, which changes their null distributions. The default "bootstrap" handles this by copying the estimation step exactly: each replication draws a composed error from the fitted model, rebuilds the response, refits, forms its own residuals, and recomputes the statistic. "asymptotic" returns a \chi^2(k-1-m) p-value for the chi-square statistic, which the authors note is conservative when evaluated at the MLE – that reference distribution belongs to the minimum-chi-square estimator. There is no asymptotic KS p-value here: with \theta estimated the Kolmogorov distribution does not apply, and Bai's (2003) martingale transformation is not implemented, so NA is returned rather than a number that would look valid.

Size and power. Measured in this package over 200 replications with B = 49, half-normal data fitted as half-normal for size and exponential data fitted as half-normal for power:

Both tests hold their size, mildly conservatively. The KS test dominates the chi-square test at every sample size, which is the paper's own conclusion and the reason it is the better default: the chi-square statistic discards the within-cell information that KS uses. Power against the half-normal/exponential pair – one of the harder distinctions, since both are one-parameter and similarly shaped – is about 0.4 at n = 200 and essentially 1 by n = 800, so a failure to reject in a small sample is weak evidence.

Value

A data frame with one row per test: test, statistic, null, p.value and B (replications actually used).

References

Wang, W. S., Amsler, C. and Schmidt, P. (2011). Goodness of fit tests in stochastic frontier models. Journal of Productivity Analysis 35, 95–118.

Bai, J. (2003). Testing parametric conditional distributions of dynamic models. Review of Economics and Statistics 85, 531–549.

See Also

spec_test, pcomposed_model, sfm

Examples

d <- as.data.frame(data_gen_cs(N = 200, rand = 5, sig_u = 1, sig_v = 0.5,
                               cons = 0.5, beta1 = 0.5, beta2 = 0.5,
                               a = 5, mu = 0.1))
f <- sfm(y_pcs ~ x1 + x2, model_name = "NHN", data = d)

## Every bootstrap replication REFITS, so keep B small in an example; use
## B = 999 for a real test.
gof_test(f, data = d, B = 39, seed = 1)

gof_test(f, data = d, test = "chisq", null = "asymptotic")

Hyvarinen score for a Normal–Half-Normal frontier fit

Description

Evaluates the Hyvarinen score used to choose the robustness tuning parameter of the maximum L_q-likelihood and minimum density-power divergence estimators, following Sugasawa and Yonekura (2021). Lower is better.

Usage

hscore(e, sigma_v, sigma_u, c, stable = TRUE)

Arguments

Details

The score is

H(c) = n^{-1}\sum_i \{2 D''(e_i) + D'(e_i)^2\},\qquad D = (f^c-1)/c,

with derivatives taken with respect to the composed residual.

The score is evaluated in logarithms. The natural-scale expression contains the fitted density raised to the power c-2; since c<1 that exponent is close to -2, so a single observation whose density underflows makes the score non-finite. The failure is silent and directional: the tuning candidates that remain evaluable are those closest to maximum likelihood, so a naive implementation reports too little robustness. Set stable = FALSE to reproduce that behaviour for comparison.

Value

A single numeric value, the mean observationwise score.

References

Hyvarinen, A. (2005). Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research 6, 695–709.

Sugasawa, S. and Yonekura, S. (2021). On selection criteria for the tuning parameter in robust divergence. Entropy 23(9), 1147.

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, density_weights

Examples

set.seed(1)
e <- rnorm(500, 0, 0.3) - abs(rnorm(500, 0, 0.6))
hscore(e, sigma_v = 0.3, sigma_u = 0.6, c = 0.20)

## the two evaluations agree wherever both are computable
hscore(e, 0.3, 0.6, 0.20) - hscore(e, 0.3, 0.6, 0.20, stable = FALSE)

## but only one of them survives an observation in the far tail
e2 <- c(e, -60)
hscore(e2, 0.3, 0.6, 0.15)                  # finite
hscore(e2, 0.3, 0.6, 0.15, stable = FALSE)  # not finite

Select the robustness tuning parameter by Hyvarinen score

Description

Refits a Normal–Half-Normal stochastic frontier at each candidate value of the robustness tuning parameter and returns the value minimising the Hyvarinen score of hscore. The search is coarse-to-fine and includes the maximum-likelihood endpoint c=0, so the criterion is free to decline robustification altogether.

Usage

hscore_select(object, method = c("mlqe", "psi", "mdpd"),
              range = c(0, 0.60), coarse = 0.01, fine = 0.001,
              half_width = 0.02, multistart = TRUE, verbose = FALSE)

Arguments

Details

Two failure modes are guarded here, both of which return a converged fit and a finite criterion value if they are not.

Silent grid loss. The criterion is evaluated in log space (see hscore). Evaluated naively it is non-finite wherever any fitted density underflows, which discards candidates non-randomly – those nearest maximum likelihood survive – and so understates the robustness the data call for.

Warm-start capture. Walking the grid upward and warm-starting each fit from the previous one is efficient but unguarded: if one fit falls into the degenerate \sigma_u \to 0 basin, every later fit inherits it. With multistart = TRUE each candidate is fitted from both the warm start and a fixed reference start, keeping the better objective.

Value

An object of class "sfa_hscore": the selected tuning parameter c, its criterion value hscore, the whole search path as a data frame, the refitted model fit, and n_finite, how many candidates were scorable out of how many tried.

Reading the result

A flat region around the minimum is not non-identification in the usual sense; it means the criterion discriminates weakly among nearby tuning values. The drop from c=0 to the minimum is a descriptive score difference, not a test of maximum likelihood against a robust alternative. Where the selected value matters, compare it with the fixed calibration from calibrate_c and check its stability by resampling: the selection can be bimodal across resamples even when the full-sample criterion clearly prefers a robust solution.

References

Sugasawa, S. and Yonekura, S. (2021). On selection criteria for the tuning parameter in robust divergence. Entropy 23(9), 1147.

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

See Also

hscore, calibrate_c, density_weights, sfm

Examples

set.seed(3)
n <- 200
x <- runif(n, 1, 10)
y <- 1 + 0.5 * log(x) + rnorm(n, 0, 0.3) - abs(rnorm(n, 0, 0.6))
dat <- data.frame(y = y, x = log(x))
fit <- sfm(y ~ x, data = dat, model_name = "NHN")

sel <- hscore_select(fit, method = "mlqe")
sel
plot(sel$path$c, sel$path$hscore, type = "l",
     xlab = "c", ylab = "Hyvarinen score")
abline(v = sel$c, lty = 3)

Test whether there is any technical inefficiency

Description

Tests H_0: no technical inefficiency, i.e. \gamma = 0 in the Battese-Corra parameterization \gamma = \sigma_u^2/(\sigma_u^2 + \sigma_v^2), following Section 3 of Coelli (1995).

Usage

inefficiency_test(object, test = c("lr_1sided", "lr", "wald", "m3t"),
                  level = 0.05)

Arguments

Details

The two tests most often reported are the two with the wrong size. H_0 puts \gamma on the boundary of the parameter space, so the usual asymptotics fail. The likelihood ratio statistic is not \chi^2(1) but the Gourieroux, Holly and Monfort (1982) mixture \tfrac12\chi^2(0) + \tfrac12\chi^2(1), whose 5% critical value is 2.71 rather than 3.84; using the \chi^2(1) value makes the naive LR test conservative. The Wald ratio \hat\gamma/s_{\hat\gamma} is worse: it is not asymptotically standard normal at the boundary at all.

Measured in this package on 400 replications per cell, half-normal data with \sigma_u = 0 and a nominal 5% level:

The one-sided LR test and the third-moment test hold their size; the naive LR test rejects about half as often as it should; the Wald test rejects three to four times too often and does not improve with n. Under the alternative (\sigma_u = \sigma_v = 1) the one-sided LR test has the better power of the two valid tests at every sample size, which is Coelli's recommendation. The Wald column's apparently high power is not comparable, because its size is wrong.

For model_name = "NTN" the null restricts both \gamma and \mu, and the mixture becomes \tfrac14\chi^2(0) + \tfrac12\chi^2(1) + \tfrac14\chi^2(2).

The "m3t" statistic is the same one skewness_test reports as test = "coelli"; it is repeated here so the four tests can be compared in a single table.

Value

A data frame with one row per test: test, statistic, null (the reference distribution), p.value and reject. The log-likelihoods under both hypotheses are attached as attributes.

References

Coelli, T. (1995). Estimators and hypothesis tests for a stochastic frontier function: A Monte Carlo analysis. Journal of Productivity Analysis 6, 247–268.

Gourieroux, C., Holly, A. and Monfort, A. (1982). Likelihood ratio test, Wald test, and Kuhn-Tucker test in linear models with inequality constraints on the regression parameters. Econometrica 50, 63–80.

See Also

skewness_test, spec_test, sfm

Examples

d <- data_gen_cs(N = 300, rand = 7, sig_u = 1, sig_v = 0.5, cons = 0.5,
                 beta1 = 0.5, beta2 = 0.5, a = 5, mu = 0.1)
f <- sfm(y_pcs ~ x1 + x2, model_name = "NHN", data = as.data.frame(d))
inefficiency_test(f)

Which observations are driving this frontier fit?

Description

Empirical influence function and case-influence measures for a fitted stochastic frontier. Answers which observations move the estimates, and whether the specification lets any single observation move them without bound.

Usage

influence_sfa(object, scale = TRUE)

Arguments

Details

The standard outlier rules do not transfer to this model. Least-median-of-squares and its relatives ignore the asymmetry of the composed error, so a genuinely large u_i – an inefficient firm, the thing being measured – reads as an outlier and is discarded, which defeats the analysis.

The influence function is the right object instead. Maximum likelihood is an M-estimator, so its influence function is a linear transformation of the score, IF_i \propto I(\theta)^{-1}\psi(y_i,\theta), and the estimator is B-robust exactly when that is bounded. Both pieces are already available: estfun supplies the per-observation scores and vcov the inverse information.

Which sensitivity to read. The raw sup-norm depends on how the parameters happen to be scaled, so it is not comparable between models with different parameter vectors. On a clean sample of 200 it reads 72.6 for "NHN" and 1322.5 for "tHN", which appears to say the robust specification is far worse; all it reflects is tHN's extra \nu on its own scale. sensitivity_std measures the same influence function in the information metric and is invariant to reparameterisation. On that scale, contaminating a single response takes "NHN" from 13.5 to 32.3 and "tHN" from 13.1 to 17.2 – Stead, Wheat and Greene's finding, and invisible in the raw number.

Neither number has an absolute threshold. The comparison to make is between specifications on the same data, not either number against a cut-off.

What to do about a large value. Nothing, by deletion. A large sensitivity is a property of the specification: under an unbounded influence function the next most extreme point simply replaces the one removed. The remedies are a specification with bounded influence – Stead, Wheat and Greene's Student's t noise term, model_name = "tHN" – or a robust divergence criterion, sfm(robust = ), whose tuning is chosen by hscore_select or calibrate_c.

What this cannot see. Like density_weights, it is computed at the fitted surface. An observation mis-recorded in a regressor can bend that surface toward itself and so appear unremarkable here. Pair this with a leverage check before concluding that an observation is not driving a result.

Value

An object of class "sfa_influence":

References

Stead, A.D., Wheat, P. and Greene, W.H. (2023). Robust maximum likelihood estimation of stochastic frontier models. European Journal of Operational Research, 309(1), 188–201.

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

See Also

density_weights, hscore_select, sfm, sfareg-methods

Examples

set.seed(3)
n <- 120
x <- runif(n, 1, 10)
d <- data.frame(y = 1 + 0.5 * log(x) + rnorm(n, 0, 0.3) - abs(rnorm(n, 0, 0.6)),
                x = log(x))

fit <- sfm(y ~ x, data = d, model_name = "NHN", keep_objective = TRUE)
influence_sfa(fit)

## Contaminate one response and compare specifications on the
## self-standardised scale.
d2 <- d; d2$y[1] <- d2$y[1] - 5
a <- sfm(y ~ x, data = d2, model_name = "NHN", keep_objective = TRUE)
b <- sfm(y ~ x, data = d2, model_name = "tHN", keep_objective = TRUE)
c(NHN = influence_sfa(a)$sensitivity_std,
  tHN = influence_sfa(b)$sensitivity_std)

Stochastic Frontier Models with Endogenous Regressors

Description

Fits the stochastic frontier model of Amsler, Prokhorov and Schmidt (2016) in which one or more regressors is endogenous, in the simultaneous-equation sense: correlated with the statistical noise v. Three estimators are provided – full information maximum likelihood, a two-step control function, and corrected 2SLS. With uhet the model becomes that of Amsler, Prokhorov and Schmidt (2017), in which environmental variables entering the inefficiency scale may themselves be endogenous.

Usage

ivsfm(formula, endogenous, instruments, data,
      model_name = c("IVLIML", "IVCF", "C2SLS"), uhet = NULL,
      inefdec = TRUE, maxit.bobyqa = 10000, maxit.psoptim = 1000,
      maxit.optim = 1000, start_val = FALSE, PSopt = FALSE,
      optHessian = TRUE, Method = "L-BFGS-B", verbose = FALSE,
      rand.psoptim = NULL)

Arguments

Details

Unlike the rest of the package neither formula nor endogenous nor instruments takes a | segment. A pipe means “variance determinant” in sfm and psfm, and reusing it here for ivreg-style instrument syntax would give the same character two meanings across entry points, so a pipe is an error rather than being reinterpreted.

The model. Writing p_i for the endogenous block and z_i for the instruments,

y_i = \beta'x_i + v_i - u_i, \qquad p_i = \Pi'z_i + \xi_i

u_i = \sigma_{U,i}|U_i|, \quad U_i \sim N(0,1), \qquad (v_i, \xi_i') \sim N(0, \Omega)

with u_i independent of (v_i, \xi_i). Endogeneity is the correlation between \xi and v: the regressors are correlated with the noise, not with the inefficiency.

Conditioning v on \xi turns this into an ordinary normal–half normal frontier in a shifted residual with a smaller noise variance,

\mu_{c,i} = \Sigma_{v\xi}\Sigma_{\xi\xi}^{-1}\xi_i, \qquad \sigma_c^2 = \sigma_V^2 - \Sigma_{v\xi}\Sigma_{\xi\xi}^{-1}\Sigma_{\xi v}

which is what makes the likelihood closed-form and needs no simulation.

The estimators.

"IVLIML"

APS (2016) Equation (13): maximum likelihood over the frontier and the reduced form jointly, so \Pi and \Sigma_{\xi\xi} are estimated alongside \beta. Efficient, and its standard errors are correct as reported.

"IVCF"

The two-step control function of Kutlu (2010), described in APS (2016) Section 4.4. \Pi and \Sigma_{\xi\xi} are fixed at their reduced-form least-squares values and only the frontier block is maximized. Cheaper, still consistent, but it discards the information about the reduced form carried in the frontier likelihood; its conventional standard errors understate the true ones unless \Sigma_{v\xi} = 0. Measured over 30 replications at n = 1000: with \rho = 0.6 the mean reported standard error on the endogenous coefficient is 0.912 of the Monte Carlo standard deviation, while "IVLIML"'s is 1.005; with \rho = 0 the two are identical, as the theory says they should be. Bootstrap, or use "IVLIML", when the standard errors matter.

"C2SLS"

APS (2016) Section 4.1: 2SLS, then \sigma_U and \sigma_V from the second and third moments of the residuals, then the intercept corrected by \sqrt{2/\pi}\,\hat\sigma_U. The residuals use the actual endogenous regressors, not their fitted values, which the paper flags explicitly. Moment-based rather than likelihood-based, so the fit carries no $opt and logLik returns NA with a warning – the same footing as psfm's "GTRE_SEQ1". It estimates no \rho, since it never models the correlation it corrects for.

Inefficiency. jlms implements APS (2016) Equations (14)–(15), which condition on \xi as well as on \varepsilon. Since \xi is correlated with v it carries information about u even though u is independent of \xi, so this is a strictly better predictor than the plain Jondrow et al. (1982) formula.

What is assumed, and what is not covered. u is assumed independent of (v, \xi). The case where the reduced-form error is correlated with the inefficiency as well requires the copula likelihood of APS (2016) Section 4.5 and is not implemented; the paper notes there that two intuitive approaches to that case do not work.

Value

An object of class "sfareg". out is a p x 3 matrix, one row per parameter, columns par, st_err and t-val. Parameters are the frontier coefficients, sigma_u, sigma_v, any delta_ terms from uhet, and one rho_ per endogenous variable (absent for "C2SLS"). Alongside the usual components:

Standard errors for \rho arrived in 1.2.0. They were previously reported as NA, on the ground that \rho = t/\sqrt{1+t't} has a non-diagonal Jacobian and a delta-method value computed as though it were diagonal would be wrong. That much was true, but the Jacobian itself is short – J = (I - \rho\rho')/s with s = \sqrt{1+t't} – so \mathrm{Var}(\hat\rho) = J\,\mathrm{Var}(\hat t)\,J' is available exactly. Checked against a numerical Jacobian in the tests, and validated end-to-end by the size of endogeneity_test, which is 0.046 at a nominal 5% over 1000 replications.

References

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

Amsler, C., Prokhorov, A. and Schmidt, P. (2017). Endogenous environmental variables in stochastic frontier models. Journal of Econometrics, 199(2), 131–140.

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

See Also

endogeneity_test for the Wald test of whether the correction was needed at all; sfm for the frontier without one.

Examples

set.seed(1)
n   <- 600
w1  <- rnorm(n); w2 <- rnorm(n); x1 <- rnorm(n)
eta <- rnorm(n)
v   <- 0.5 * (0.6 * eta + sqrt(1 - 0.6^2) * rnorm(n))
x2  <- 0.9 * w1 - 0.7 * w2 + 0.5 * x1 + eta
y   <- 0.5 + 0.8 * x1 - 0.6 * x2 + v - abs(rnorm(n))
dat <- data.frame(y = y, x1 = x1, x2 = x2, w1 = w1, w2 = w2)

fit <- ivsfm(y ~ x1 + x2, endogenous = ~x2, instruments = ~ w1 + w2,
             data = dat)
fit$out

## Ignoring the endogeneity biases the coefficient on x2 toward zero.
sfm(y ~ x1 + x2, model_name = "NHN", data = dat)$out

Latent Class Stochastic Frontier Models

Description

Fits the latent class stochastic frontier model of Greene (2005) and Orea and Kumbhakar (2004), in which several unobserved technologies coexist in one sample and each firm contributes to every one of them, weighted by its class probability.

Usage

lcsfm(formula, model_name = c("LCM", "LCM_Z", "LCM_CN"), 
data, n_class = 2, maxit.bobyqa = 10000, maxit.psoptim = 1000, 
maxit.optim = 1000, REPORT = 1, trace = 0, pgtol = 0, 
start_val = FALSE, PSopt = FALSE, optHessian = TRUE, 
inefdec = TRUE, penalty_c = 0, upper = NA, Method = "L-BFGS-B", 
verbose = FALSE, rand.psoptim = NULL)

Arguments

Details

"LCM_CN", the contaminated normal frontier. A different use of the same machinery: instead of separating technologies, it gives the noise heavier tails. The noise density is a scale mixture of normals, f_v(v) = \sum_j p_j \phi(v/\sigma_{vj})/\sigma_{vj}, while the frontier, \sigma_u and the mixing proportions stay common across components. At J = 2 this is the scale-contaminated normal, an alternative to a Student-t noise term.

It has a closed form, which is why it is cheap. Since f_\epsilon(e) = \int f_v(e + Su) f_u(u)\,du is linear in f_v, a mixture noise density passes straight through:

f_\epsilon(e) = \sum_j p_j\, f_{NHN}(e; \sigma_{vj}, \sigma_u)

a mixture of ordinary normal/half-normal densities sharing one \sigma_u. No extra integration is required, and the identity is verified against direct numerical integration in the tests.

The parameter vector is (\sigma_{v1},\dots,\sigma_{vJ}, \sigma_u, \beta, \mathrm{logit}) rather than one block per class.

It is also the specification lcsfm_homogeneity”'s asymptotic null is stated for: one scalar parameter differing between components. On "LCM", where everything varies, the \chi^2_{0:1} null rejected 63.5% of the time at a nominal 5% under a true null; on "LCM_CN" the same null gives 8.0%. Still mildly liberal, and the test says so, but usable.

The latent class stochastic frontier supposes that J distinct technologies are present in one sample and that which technology a firm operates is unobserved. Each class j has its own frontier and its own two scales,

y_i = x_i'\beta_j + v_{ij} - u_{ij}, \quad v_{ij} \sim N(0, \sigma_{vj}^2), \quad u_{ij} \sim N^+(0, \sigma_{uj}^2),

and class membership follows a multinomial logit,

P(j \mid q_i) = \exp(q_i'\delta_j) / \sum_m \exp(q_i'\delta_m), \quad \delta_J = 0,

with class J the reference. Every firm contributes to every class, weighted by its class probability, so the log-likelihood mixes the J composed densities rather than assigning firms to groups:

\log L_i = \log \sum_j P(j \mid q_i) f_j(\varepsilon_{ij}).

Under "LCM" the q_i are a constant, so the class probabilities are the same for every firm; under "LCM_Z" they are the covariates in the second part of the formula.

This is a different object from fitting separate frontiers to known subgroups: the groups are not known, and the uncertainty about them is carried through into the efficiency predictions. It is also more general than zsfm's zero-inefficiency model, which is the restricted two-class case in which one class has no inefficiency at all.

post.prob is an n \times J matrix of posterior class probabilities, jlms_class holds the class-conditional Jondrow et al. (1982) predictions, and jlms is their posterior-weighted average – a firm's inefficiency is only defined relative to a frontier, and which frontier it faces is itself uncertain, so reporting the modal class alone would discard that. class gives the modal class and class_prob the average prior class shares.

Label switching. Class labels in a finite mixture are identified only up to permutation: relabelling the classes and permuting \delta gives the identical likelihood. Starting values are built by splitting the OLS residuals at their J-quantiles and refitting within each group, which makes a run reproducible on the same data – and avoids the saddle point that identical starting components would sit at, since equal posterior class probabilities give a zero score with respect to the class split – but it does not pin the labels to any external ordering. Do not compare class 1 across two fits, or against a simulation's true labels, without matching the classes first.

Choosing J. The likelihood cannot fall when a class is added, and the usual regularity conditions for a likelihood-ratio test of J against J+1 do not hold, since the null puts a class probability on the boundary. Use AIC()/BIC(), the separation of the fitted frontiers, and whether the extra class is economically interpretable, rather than a p-value.

Value

An object of class "sfareg" containing the following components:

Note

Standard errors require optHessian set to TRUE

Author(s)

Christopher F. Parmeter and David H. Bernstein

References

Greene, W. (2005) 'Reconsidering heterogeneity in panel data estimators of the stochastic frontier model', Journal of Econometrics, 126(2), pp. 269-303. doi:10.1016/j.jeconom.2004.05.003.

Orea, L. and Kumbhakar, S.C. (2004) 'Efficiency measurement using a latent class stochastic frontier model', Empirical Economics, 29(1), pp. 169-183. doi:10.1007/s00181-003-0184-2.

Caudill, S.B. (2003) 'Estimating a mixture of stochastic frontier regression models via the em algorithm: a multiproduct cost function application', Empirical Economics, 28(3), pp. 581-598. doi:10.1007/s001810200142.

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. doi:10.1016/0304-4076(82)90004-5.

See Also

lcsfm_homogeneity, zsfm, sfm, psfm, ttsfm, sfareg-methods

Examples

library(sfa)

## Two coexisting technologies, class probabilities constant.
set.seed(42)
n   <- 400
x1  <- rnorm(n); x2 <- rnorm(n)
cls <- rbinom(n, 1, 0.4) + 1L
b   <- rbind(c(1, 0.5, 0.5), c(3, 1.0, 0.2))
sv  <- c(0.2, 0.2); su <- c(0.5, 1.0)
y   <- b[cls, 1] + b[cls, 2] * x1 + b[cls, 3] * x2 +
       rnorm(n, 0, sv[cls]) - abs(rnorm(n, 0, su[cls]))

fit <- lcsfm(y ~ x1 + x2, model_name = "LCM",
             data = data.frame(y, x1, x2), n_class = 2)
fit$class_prob        # average class shares
head(fit$post.prob)   # posterior class membership, one row per firm

Test a latent class frontier against homogeneity

Description

Tests whether a lcsfm fit's latent classes are needed at all, against the null that one class describes the data. This is the first question a latent class fit raises and the ordinary likelihood ratio test cannot answer it, because under the null the second class's parameters are unidentified, the mixing proportion sits on a boundary, and the classes are interchangeable.

Usage

lcsfm_homogeneity(object, null = c("bootstrap", "chisq01"),
                  B = 199, c = 1, level = 0.05,
                  seed = NULL, quiet = FALSE, envir = parent.frame())

Arguments

Details

The statistic. lcsfm is refitted maximising the modified likelihood – the ordinary log-likelihood plus Chen et al.'s penalty c\log(J^J\prod_j p_j), which diverges to -\infty as any class probability approaches 0 or 1 and so keeps every class occupied. The statistic is twice the difference between that maximum and the maximised log-likelihood of the one-class model, which is sfm(model_name = "NHN") on the same formula, since a single latent class is the ordinary frontier. The penalty is exactly zero at equal class probabilities, so the null model is recovered without a handicap.

Why the default null is a bootstrap. Stead, Wheat and Greene report the modified statistic is asymptotically \chi^2_{0:1}, a 50:50 mixture of a point mass at zero and \chi^2_1. That result is stated for the case in which a scalar parameter differs between classes – Zhu and Zhang (2004) require \theta_1 and \theta_2 scalar, and the paper's own application holds every parameter except the noise scale common across classes.

lcsfm's "LCM" does not do that: it lets \sigma_v, \sigma_u and the whole frontier vector vary by class. The limit therefore does not carry over, and this was measured rather than assumed. Over 200 replications of a one-class data generating process at n = 400:

with the sampling distribution tracking \chi^2_5 – exactly the parameter-count difference – rather than \chi^2_{0:1}, whose median is 0 against an observed 4.21.

A parametric bootstrap makes no appeal to an asymptotic distribution: it simulates from the fitted null and re-runs the whole procedure, so it is valid whatever the parameter structure. Over 60 replications of the same one-class process with B = 99:

both within a binomial standard error of nominal, with the p-values indistinguishable from uniform (KS against Uniform(0,1), p = 0.80). On a sample that "chisq01" rejected at p = 0.005, the bootstrap returns p = 0.32, and on genuinely two-class data the statistic is 17 times the largest of 99 null draws.

Cost. The bootstrap is B+1 pairs of refits. Budget accordingly, and use seed for reproducibility.

"LCM_Z" is refused, by both this function and lcsfm(penalty_c = ): it makes the class probabilities depend on covariates, and the penalty is defined on a scalar mixing proportion.

Value

An object of class c("sfa_mlrt", "htest"), so it prints like any other R test. Alongside statistic and p.value it carries null, B, boot (the bootstrap null draws), boot_ok, penalty_c, penalty_at_max, the penalised and unpenalised J-class log-likelihoods, logLik_1, class_prob, and reject.

References

Chen, H., Chen, J. and Kalbfleisch, J.D. (2001). A modified likelihood ratio test for homogeneity in finite mixture models. Journal of the Royal Statistical Society B, 63(1), 19–29.

Stead, A.D., Wheat, P. and Greene, W.H. (2023). On hypothesis testing in latent class and finite mixture stochastic frontier models, with application to a contaminated normal-half normal model. Journal of Productivity Analysis, 60(1), 37–48.

Zhu, H.-T. and Zhang, H. (2004). Hypothesis testing in mixture regression models. Journal of the Royal Statistical Society B, 66(1), 3–16.

See Also

lcsfm, TIC, skewness_test

Examples

set.seed(21)
n   <- 300
x1  <- rnorm(n)
cls <- rbinom(n, 1, 0.5)
d   <- data.frame(y = ifelse(cls == 1, 4, 1) + x1 +
                      rnorm(n, 0, 0.4) - abs(rnorm(n, 0, 0.6)),
                  x1 = x1)

f <- lcsfm(y ~ x1, model_name = "LCM", data = d, n_class = 2)

## Small B for the example only; use the default in real work.
lcsfm_homogeneity(f, B = 19, seed = 1, quiet = TRUE)

Marginal Effects of the Variance Determinants on Inefficiency

Description

For a model whose inefficiency scale depends on covariates (the _Z models), the marginal effects of those covariates on the mean and variance of the inefficiency term, \partial E[u]/\partial z_k and \partial \mathrm{Var}[u]/\partial z_k.

Usage

marginal_effects(object, average = FALSE, component = c("u", "h"))

Arguments

Details

A \delta coefficient is not interpretable on its own: it sits in a log link for a scale parameter, so its units are neither those of u nor those of z, and its magnitude is not comparable across models whose link differs. The marginal effect is on the scale of u, and is what applied work reports.

Writing s = \sigma_u(z), every derivative reduces to \partial s/\partial z_k times a constant that cancels back into the moment itself. For the half-normal and the exponential alike:

\textrm{SD link:} \quad \partial E[u]/\partial z_k = \delta_k E[u], \qquad \partial \mathrm{Var}[u]/\partial z_k = 2\delta_k \mathrm{Var}[u]

\textrm{variance link:} \quad \partial E[u]/\partial z_k = (\delta_k/2) E[u], \qquad \partial \mathrm{Var}[u]/\partial z_k = \delta_k \mathrm{Var}[u]

The factor of two matters. sfm()'s "NHN_Z" and "NE_Z" place the linear predictor on the standard deviation, \sigma_u = \exp(z'\delta), while psfm()'s "TRE_Z" and "GTRE_Z" place it on the variance, \sigma_u = \sqrt{\exp(z'\delta)}. Reading a \delta from one family with the other's convention in mind is wrong by exactly this factor. Reporting the marginal effect instead of the coefficient removes the trap, since the effect is on the scale of u either way.

The link is read from the fit rather than assumed, so a sfm model fitted with z_link = "var" and a psfm model are handled the same way and their effects are directly comparable.

psfm(model_name = "GTRE_Z") separates persistent inefficiency h_i from transient u_{it} and gives each its own determinant block, so both are available: component = "u" for the transient effects and component = "h" for the persistent ones. Each block is located by NAME rather than by position, because the \sigma_h coefficients are the trailing ones and a positional rule would report them under a \sigma_u label.

Effects are per observation, because \sigma_u varies with z. Constant columns of the z design are dropped: the derivative with respect to an intercept is not a marginal effect.

Standard errors are deliberately not reported. The effect is a nonlinear function of both \delta and z, so its sampling distribution requires the delta method through the full covariance matrix or a bootstrap; a number printed without one would be read as inference.

Value

With average = FALSE, a data frame with one row per observation: sigma_u, E_u, Var_u, then one dE_u.d<z> and one dVar_u.d<z> column per non-constant variance determinant. Attributes "average" (the average marginal effects), "link" and "family" record the convention used.

With average = TRUE, just the named vector of average marginal effects.

See Also

sfm, efficiency_ci

Examples

dat <- data_gen_cs(N = 300, rand = 1, sig_u = 1, sig_v = 0.3,
                   cons = 0.5, beta1 = 0.5, beta2 = 0.5, a = 5, mu = 0.5)
fit <- sfm(y_pcs_z ~ x1 + x2 | z, data = dat, model_name = "NHN_Z")

me <- marginal_effects(fit)
head(me)

## The conventional summary.
marginal_effects(fit, average = TRUE)

Model-Implied Mean Efficiency

Description

The mean efficiency implied by the fitted distribution of inefficiency, E[\exp(-U)], and the mean efficiency among the most efficient p of that distribution.

Usage

meanefficiency(object,
               p = c(0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99),
               use_closed_form = TRUE)

Arguments

Details

A different quantity from the reported average. print and summary report the average of the per-observation conditional predictions E[\exp(-u_i) \mid \varepsilon_i]. That is a property of the sample in hand. E[\exp(-U)] is a property of the fitted model, and it is the number that can be compared across studies.

Closed forms. Available for the half-normal (2e^{\sigma_U^2/2}\Phi(-\sigma_U)), exponential (1/(1+\sigma_U)), truncated normal, Rayleigh, gamma ((1+\theta)^{-k}), uniform ((1-e^{-\theta})/\theta) and generalized exponential (2/((\sigma_U+1)(\sigma_U+2))). The lognormal, Weibull and Nakagami have none in elementary terms and are integrated numerically; method in the returned list says which route was taken.

Every closed form is checked against numerical integration of the same density in the package's tests, and the density is separately checked against its own quantile function by simulation. The first catches an algebraic slip; the second catches a parameter being read off the fit incorrectly, which would otherwise make both routes agree and both be wrong.

Supra-percentile means are always integrated. Efficiency is decreasing in u, so the most efficient p are the smallest p of u: the reported value is p^{-1}\int_0^{Q(p)} e^{-u}f(u)\,du. It necessarily exceeds the unconditional mean and falls as p rises.

Value

A list with model, distribution, the parameters used, mean_efficiency, the method used to obtain it, and supra, a data frame of p against mean efficiency.

See Also

efficiency for the per-observation predictors.

Examples

d <- data_gen_cs(N = 400, 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, model_name = "NHN", data = d)

r <- meanefficiency(fit)
r$mean_efficiency
r$supra

What skewness and excess kurtosis can a pair of distributions produce?

Description

Reports, for each pair of noise and inefficiency distributions, the set of skewness and excess kurtosis values the composed error \varepsilon = v \pm u can attain at all, and says which pairs the residual moments rule out. This is the Type II failure check of Papadopoulos and Parmeter (2021): a pair whose attainable set does not contain the sample moments is refuted without any test statistic. It needs no frontier fit.

Usage

moment_range(object = NULL,
             noise = c("normal", "logistic", "laplace", "uniform", "t"),
             inefficiency = c("halfnormal", "exponential", "genexponential",
                              "truncatednormal", "gamma"),
             inefdec = NULL, df = 8)

Arguments

Details

The two equations. With R = \mathrm{SNR}^2/(1+\mathrm{SNR}^2) the share of composed-error variance contributed by u, Papadopoulos and Parmeter (2021, eqs. 5 and 10) give

|\gamma_1(\varepsilon)| = \gamma_1(u)\,R^{3/2}, \qquad \gamma_2(\varepsilon) = \gamma_2(v)(1-R)^2 + \gamma_2(u)R^2.

Because R lies in (0,1), the composed skewness can never exceed the skewness of u itself, whatever the noise: the Normal-Half-Normal pair is confined to \gamma_1(\varepsilon) \in (-0.995, 0) and the Normal-Exponential pair to (-2, 0). The kurtosis equation is quadratic in R, so its range is not simply the interval between the two components' values – with leptokurtic noise it dips below both before rising, reaching \gamma_2(v)\gamma_2(u)/(\gamma_2(v)+\gamma_2(u)) at R = \gamma_2(v)/(\gamma_2(v)+\gamma_2(u)).

Which u can fail. Half-normal, exponential and generalized exponential have fixed skewness and kurtosis, so their attainable sets are the narrowest. The truncated normal is indexed by \mu/\sigma_u and sweeps \gamma_1(u) \in (0,2) and \gamma_2(u) \in (-0.243, 6) – it is the only one here whose composed error can be platykurtic under normal noise. The gamma has \gamma_1(u) = 2/\sqrt{k} and \gamma_2(u) = 6/k, unbounded as k \to 0, so a gamma inefficiency term can never suffer a Type II failure; the price, as PP2021 note, is that k < 1 makes the density log-convex.

This is a diagnostic, not a test. PP2021 section 2.3 shows the range check over-rejects, and the package's own replication of their Tables 1-5 confirms it: with correctly specified Normal-Half-Normal data at n = 200 and \mathrm{SNR} = 3, the sample excess kurtosis falls outside its theoretical range 26% of the time, and the sample skewness 16% of the time. Use spec_test for the formal statistic, which is built on the same two equations but studentised and bootstrapped.

Relation to the wrong-skew problem. A Type I failure – residual skewness of the wrong sign, Waldman (1982) – is a property of the data and is reported for every row at once. Type II failures are pair-specific: they are what tells you that the residuals are too skewed, or too heavy-tailed, for the particular u you assumed. Wrong skew is handled by skewness_test and, as an estimator, by esfm.

Value

A data frame of class "sfa_moment_range", one row per pair, with the attainable interval for the composed skewness (g1_lo, g1_hi) and excess kurtosis (g2_lo, g2_hi), logical skew_ok and kurt_ok, and a verdict of "ok", "type I (wrong skew)" or "type II (...)". The sample moments, sample size and orientation are attributes.

References

Papadopoulos, A. and Parmeter, C. F. (2021). Type II failure and specification testing in the Stochastic Frontier Model. European Journal of Operational Research 293(3), 990-1001.

Waldman, D. M. (1982). A stationary point for the stochastic frontier likelihood. Journal of Econometrics 18(2), 275-279.

See Also

spec_test for the formal test built on the same two equations, skewness_test for the Type I failure, esfm for an estimator that tolerates it.

Examples

## The reference table on its own: no data needed.
moment_range(noise = "normal")

set.seed(4)
d <- data_gen_cs(N = 400, rand = 4, sig_u = 1, sig_v = 0.5, cons = 1,
                 beta1 = 0.5, beta2 = 0.5, a = 1, mu = 0.5)

## Exponential inefficiency, read through a normal-half-normal lens: the
## residuals are more skewed than a half-normal can make them.
fit <- sfm(y_pcs_e ~ x1 + x2, data = d, model_name = "NHN")
moment_range(fit, noise = "normal")

Nonparametric Stochastic Frontier Models

Description

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.

Usage

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)

Arguments

Details

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.

Value

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:

Note

Requires the np package, which is in Suggests. Bandwidth cross-validation dominates the run time and scales quadratically in the sample size.

Author(s)

Christopher F. Parmeter and David H. Bernstein

References

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.

See Also

sfm for parametric cross-sectional frontiers, psfm for panel models, and data_gen_cs for simulating data with known true parameters.

Examples

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)
}

Panel89

Description

The dataset is a cross-section of U.S. commercial banks for 1989, extracted from the panel dataset used by Kumbhakar, Parmeter and Tsionas (2013) and based on the Federal Reserve Bank of Chicago's Reports of Condition and Income. It contains detailed cost data with inputs and outputs defined under the intermediation approach, and input prices constructed as expense-quantity ratios.

Usage

data("panel89")

Format

A data frame with 4,985 observations on the following 11 variables.

y

a numeric vector containing the natural logarithm of total cost*

q1

a numeric vector containing the natural logarithm of installment loans

q2

a numeric vector containing the natural logarithm of real estate loans

q3

a numeric vector containing the natural logarithm of business loans

q4

a numeric vector containing the natural logarithm of federal funds sold and securities purchased

q5

a numeric vector containing the natural logarithm of other assets

w1

a numeric vector containing the natural logarithm of the price of labour*

w2

a numeric vector containing the natural logarithm of the price of capital*

w3

a numeric vector containing the natural logarithm of the price of purchased funds*

w4

a numeric vector containing the natural logarithm of the price of interest-bearing deposits in total transaction accounts*

z

a numeric vector containing the natural logarithm of total assets

Details

*The cost and input price variables are normalised by that of a fifth input: the price of interest-bearing deposits in total non-transaction accounts. Total cost is defined as the sum of total expenses for each input. Input prices are derived by dividing the total expense for each input by the corresponding input quantity.

Source

Kumbhakar, S.C., Parmeter, C.F. and Tsionas, E.G. (2013) 'A zero inefficiency stochastic frontier model', Journal of Econometrics, 172(1), pp. 66-76. doi:10.1016/j.jeconom.2012.08.021.

References

Kumbhakar, S.C. and Tsionas, E.G. (2005) 'Measuring technical and allocative inefficiency in the translog cost system: a Bayesian approach', Journal of Econometrics, 126(2), pp. 355-384. doi:10.1016/j.jeconom.2004.05.006.

Examples

data(panel89)
  plot(panel89)

The distribution of the composed error

Description

Density and distribution function of the stochastic frontier composed error \varepsilon = v - u (production) or v + u (cost), with v \sim N(0, \sigma_v^2) and u \sim N^+(0, \sigma_u^2) independent. The density is a skew-normal and closed form; the distribution function is not, and is obtained here by quadrature that stays accurate far into the tails.

Usage

dcomposed(x, sigma_u, sigma_v, inefdec = TRUE, log = FALSE)

pcomposed(q, sigma_u, sigma_v, inefdec = TRUE, lower.tail = TRUE,
          log.p = FALSE, method = c("quadrature", "simulate"),
          n_nodes = 128, R = 1e6, seed = NULL)

Arguments

Details

Why this exists. Every likelihood in sfa evaluates the composed error's density and none needs its distribution function, which is why there was not one. Copula models need it: a Gaussian copula on the composed error contains \Phi^{-1}(F(\varepsilon)), so an F that saturates at 0 or 1 does not merely lose precision, it returns an infinity and takes the copula density with it.

How it is computed. Conditioning on u leaves a normal distribution function in closed form, so

F(q) = E_u\left[\Phi\left((q + su)/\sigma_v\right)\right],

with s = +1 for a production frontier and -1 for a cost frontier. The noise is therefore never drawn or integrated, and the tails inherit the accuracy of pnorm, which is reliable twenty standard deviations out. The expectation over u is taken by Gauss-Legendre quadrature on u = \sigma_u t/(1-t), and accumulated by log-sum-exp so that terms which would underflow individually still contribute.

The upper tail. Requesting lower.tail = FALSE negates the argument of \Phi rather than forming 1 - F. At q = 30 with \sigma_u = \sigma_v = 1 the complement is exactly 1 in double precision while the direct computation returns a finite log-probability below -10^{2}.

Accuracy. There is no closed form except at \lambda = \sigma_u/\sigma_v equal to 0 (the normal) or 1, where Amsler, Schmidt and Tsay show P(Q) = \Phi(Q/(\sqrt{2}\sigma_u))^2 for the cost case. Against that standard the default quadrature is accurate to about 10^{-13} relative at probabilities as small as 10^{-115}. It is deliberately used everywhere rather than special-cased at \lambda = 1: the production form of that identity is a complement, which would reintroduce the saturation the function exists to avoid.

Value

A numeric vector the length of x or q.

References

Amsler, C., Schmidt, P. and Tsay, W.-J. (2019). Evaluating the cdf of the distribution of the stochastic frontier composed error. Journal of Productivity Analysis, 52(1), 29–35.

See Also

copsfm, sfm

Examples

## Density integrates to one.
integrate(function(x) dcomposed(x, 1, 0.5), -Inf, Inf)$value

## Distribution function, and a tail that does not saturate.
pcomposed(c(-2, 0, 2), sigma_u = 1, sigma_v = 0.5)
pcomposed(30, 1, 1, lower.tail = FALSE, log.p = TRUE)

## Checked against the one exact standard there is (cost frontier).
su <- sv <- sqrt(0.5)
c(exact = pnorm(-6 / (sqrt(2) * su))^2,
  computed = exp(pcomposed(-6, su, sv, inefdec = FALSE, log.p = TRUE)))

Composed-error distribution function for any cross-sectional model

Description

P(\varepsilon \le q) for the composed error of any cross-sectional model, where pcomposed covers only the half-normal case.

Usage

pcomposed_model(q, model_name, par, inefdec = TRUE, lower.tail = TRUE,
                log.p = FALSE, n_nodes = 128)

composed_cdf(object, q = NULL, data = NULL, ...)

Arguments

Details

The probability is computed as an expectation over the inefficiency term,

F(q) = E_u\left[F_v\left((q + su)/\sigma_v\right)\right],

with s = +1 for production and -1 for cost. Conditioning on u leaves the noise CDF in closed form, so the noise is never integrated and the tails inherit pnorm()'s (or pt()'s) own accuracy.

Three parameterizations do not read the way their names suggest, and are taken from each likelihood in sfm() rather than from the label: "THT" lists sigu before sigv, the only model that does; "NG"'s sigu is the gamma scale and its mu the shape; "NNAK"'s sigu is the Nakagami spread, so \Omega = \code{sigu}^2; and "NLN"'s mu is a meanlog.

"THT" is not an independent convolution. Tancredi's (2002) composed error is skew-t, a scale mixture in which v and u are divided by the same \sqrt{V/a}. Treating it as t noise plus an independent half-normal – which is what "tHN" actually is – gets the log density wrong by up to 1.74, so the two models take different paths here despite looking alike.

Verified two ways: against pcomposed for the half-normal case, agreeing to 3.5\times10^{-15} relative including a lower tail at \log F = -209; and for all thirteen models against a four-million-draw simulation, with a maximum absolute error of 4\times10^{-4} against a Monte Carlo standard error of 7.5\times10^{-4}. Numerically differentiating this CDF reproduces each model's own stored log-density to about 10^{-8}, except for the simulated-ML models "NLN" and "NW", where the difference is the simulation error of their likelihood rather than a disagreement about the model.

Value

A numeric vector of probabilities (or their logs).

See Also

pcomposed, dcomposed, sfm

Examples

p <- c(sigv = 0.5, sigu = 1.1)
pcomposed_model(c(-2, -1, 0, 1), "NE", p)

d <- as.data.frame(data_gen_cs(N = 200, rand = 5, sig_u = 1, sig_v = 0.5,
                               cons = 0.5, beta1 = 0.5, beta2 = 0.5,
                               a = 5, mu = 0.1))
f <- sfm(y_pcs ~ x1 + x2, model_name = "NHN", data = d)
head(composed_cdf(f, data = d))

Predictions, fitted values and residuals for stochastic frontier fits

Description

Standard modelling generics for objects of class "sfareg", as returned by sfm, psfm, zsfm and ttsfm.

Usage

## S3 method for class 'sfareg'
predict(object, newdata = NULL,
        type = c("frontier", "response", "efficiency"), ...)
## S3 method for class 'sfareg'
fitted(object, ...)
## S3 method for class 'sfareg'
residuals(object, ...)

Arguments

Details

Only the first part of a pipe formula (y ~ x | z | zp) describes the frontier; later parts parameterize variances and take no part in prediction.

Frontier coefficients are identified by matching the fitted coefficient names against the frontier design matrix, so these methods work for every model in the package without depending on where a given model happens to place its auxiliary parameters (lambda/sigma for NHN, sigv/sigu for NE, sigmaSq/gamma for PL80, and so on).

type = "efficiency" and type = "response" are available only for the estimation sample, since both require the composed residual and therefore the response; requesting them with newdata raises an error rather than returning something misleading. "efficiency" additionally requires a model that returns an efficiency prediction – see sfm for which ones do.

Value

A numeric vector, one element per observation.

Author(s)

David H. Bernstein

See Also

sfm, psfm, sfareg-methods

Examples

library(sfa)
d <- data_gen_cs(N = 200, rand = 1, sig_u = 0.8, sig_v = 0.3,
                 cons = 1, beta1 = 0.5, beta2 = 0.5, a = 5, mu = 0.5)
f <- sfm(y_pcs ~ x1 + x2, model_name = "NHN", data = d)

head(fitted(f))
head(residuals(f))
head(predict(f, type = "efficiency"))
head(predict(f, newdata = d[1:5, ]))

Print method for "sfareg" objects

Description

Compact display of a stochastic frontier fit, as returned by sfm, psfm, zsfm or ttsfm.

Usage

## S3 method for class 'sfareg'
print(x, ...)

Arguments

Details

Prints the model name, the estimated coefficients with their standard errors and t-values (the out matrix), and the log-likelihood where the estimator has one. For a fuller display, including the frontier and variance-parameter blocks, use summary.sfareg.

Value

No return value, called for side effects

Author(s)

David H. Bernstein

See Also

summary.sfareg, sfareg-methods

Examples

library(sfa)     

cs_data_trial   <- 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)

cs.nhnz     <-  sfm(formula    = y_pcs_z ~ x1 +x2| z,    model_name = "NHN_Z",
                    data       = cs_data_trial,          PSopt      = TRUE)
print(cs.nhnz)

psfm

Description

Function to implement various panel data stochastic frontier estimators

Usage

psfm(formula, model_name = c("TRE_Z", "GTRE_Z", "TRE",
                    "GTRE", "GTRE_FML", "TFE", "TFE_WMLE", "FD",
                    "GTRE_SEQ1", "GTRE_SEQ2", "SSFE",
                    "PL80", "PL80_MVTN", "BC92", "K1990", "K1990modified",
                    "CSS", "LS", "SSRE", "SSCRE", "KSS"), data,
                    maxit.bobyqa = 5000, maxit.nlminb = 500,
                    maxit.psoptim = 100, maxit.optim = 1000,
                    REPORT = 1, trace = 3, pgtol = 0, individual,
                    halton_num = NULL, start_val = FALSE, gamma = FALSE,
                    PSopt = FALSE, optHessian = TRUE, inefdec = TRUE, Method = "L-BFGS-B",
                    verbose = FALSE, rand.gtre = NULL, rand.psoptim = NULL,
                    OPG_calc = FALSE,
                    estimator = c("fiml", "sml", "seq1", "seq2"),
                    collinear_action = c("start_only", "error", "warn_drop"),
                    time = NULL, tfe_lambda_max = 100,
                    keep_objective = FALSE, kss_L = "auto",
                    kss_smooth = "auto", kss_L_max = 7L, mundlak = NULL)

Arguments

Details

The generalized true random effects model (GTRE, 4-component model) and true random effects models (TRE) are both estimated by simulated maximum likelihood based on the paper by the Fillipini and Greene (2016, JPA). The TRE_Z and GTRE_Z allow for modeling the u-component of the GTRE and TRE with determinants of inefficiency. The first-difference estimator (FD) of Wang and Ho (2010, JoE) is also available.

Choosing how GTRE is estimated. Since version 1.1.4 the four routes to the four-component model are selected with estimator rather than by picking among four model_name values that looked like four different models. model_name = "GTRE" now defaults to estimator = "fiml"; through 1.1.3 it fitted the simulated-ML estimator, which is estimator = "sml". Scripts that pass "GTRE" therefore get a different estimator than they did and are warned once per call. The explicit names "GTRE_FML", "GTRE_SEQ1" and "GTRE_SEQ2" are unaffected.

Two ways of estimating the same four-component model. "GTRE" and "GTRE_FML" both fit

y_{it} = \beta_0 + x_{it}'\beta + r_i - h_i + v_{it} - u_{it},

with r_i, v_{it} two-sided and h_i, u_{it} one-sided (half-normal), separating persistent inefficiency h_i from transient inefficiency u_{it}. They differ in how the likelihood is obtained.

"GTRE" uses simulated maximum likelihood over Halton draws. "GTRE_FML" uses full information maximum likelihood through the model's closed-skew-normal representation: with A = -[1_T, I_T], V = \mathrm{diag}(\sigma_h^2, \sigma_u^2 I_T) and \Sigma = \sigma_v^2 I_T + \sigma_r^2 1_T 1_T', the firm-level density is available in closed form and nothing is integrated by Monte Carlo. The likelihood is therefore deterministic, and repeated fits to the same data give identical answers.

The (T+1)-dimensional normal CDF this requires would normally be the bottleneck, but its covariance is diagonal-plus-rank-one for this model, so it collapses to a one-dimensional Gauss-Hermite integral – the same reduction used for "TFE".

"GTRE_FML" reports raw standard deviations sigr, sigv, sigh, sigu (individual effect, noise, persistent inefficiency, transient inefficiency), not the \lambda/\sigma reparameterization "GTRE" uses, and it requires a balanced panel – the closed-skew-normal pieces are built at a single T. It errors rather than silently mis-estimating if the panel is unbalanced; use "GTRE" in that case.

Starting values for "GTRE_FML". The closed-skew-normal likelihood carries a boundary optimum at \sigma_h = 0, where the model collapses to "TRE" and the intercept absorbs the missing E[h] = \sigma_h\sqrt{2/\pi}. Two different things can put a fit there, and the likelihood tells them apart. A fit whose boundary solution has a lower likelihood than the true parameter vector reached it by getting stuck, and is an optimizer failure; one whose boundary solution has a higher likelihood reached it because that is where the maximum is, and is the correct answer for that sample (see ‘When a persistent scale collapses to zero’ below). The start search addresses only the first kind. Following Colombi (2010) and Colombi, Martini and Vittadini (2011), who recommend seeding this likelihood from the two-step moment estimates, psfm() builds two candidate starts – the random-effects decomposition shared with "GTRE", and the two-step estimator reported by "GTRE_SEQ2" – evaluates the likelihood at each, and begins from the better one. Which was chosen, and the log-likelihood at both, is returned in start_search. Supplying start_val overrides the search.

Note that \sigma_r and \sigma_h enter the density only through \sigma_r^2 + \sigma_h^2; their separation rests entirely on the skewness term. That pair is therefore the least precisely estimated part of the model even when the fit is at the global optimum, which is consistent with the confounding between latent heterogeneity and persistent inefficiency reported in the source papers.

"GTRE_SEQ1" and "GTRE_SEQ2" are large-T estimators. Both take the random-effects decomposition from plm and hand its two output series to a second stage – "GTRE_SEQ1" fits an intercept-only normal–half-normal likelihood to each, "GTRE_SEQ2" inverts their second and third moments. The second stage treats those series as if they were draws from v_{it}-u_{it} and r_i-h_i, and at finite T they are not: alpha_hat is the shrunken predictor c\,(\alpha_i+\bar e_i) with c = T\sigma_\alpha^2/(T\sigma_\alpha^2+\sigma_\epsilon^2), and eps_hat is the quasi-demeaned residual (1-\theta)\alpha_i + e_{it} - \theta\bar e_i with \theta = 1-\sigma_\epsilon/\sqrt{T\sigma_\alpha^2+\sigma_\epsilon^2}.

Relative to the latent draws this attenuates the third central moment on the \epsilon side by (1-\theta/T)^3+(T-1)(-\theta/T)^3 (about 0.81 at T=10) and scales the \alpha side by c^2 and c^3. None of these factors tends to 1 as N grows with T fixed, so both estimators converge to the wrong constants and are inconsistent for fixed T; each call warns accordingly. On the package's own convergence design (T=10, \sigma_h=0.8, \sigma_r=0.4) the estimates settle at \gamma_{uv}\approx0.87 against 0.917 and \sigma^2_{hr}\approx0.68 against 0.80, and the log MSE-on-log N slopes are near zero rather than -1. The same routines applied to the latent draws themselves recover the truth at the \sqrt{N} rate, which locates the problem in the first-stage decomposition rather than in either second stage. Raising T removes the bias; raising N only sharpens it. For a consistent alternative use model_name = "GTRE", which estimates all four variance components jointly.

The standard errors these two report are delta-method errors built from the iid moment variances, so they describe sampling variability around each estimator's probability limit, not around the truth. They are also optimistic on the \epsilon side, whose n residuals are correlated within firm.

The two true-fixed-effects estimators, and a change of meaning in "TFE". Both "TFE" and "TFE_WMLE" fit the same model,

y_{it} = \alpha_i + x_{it}'\beta + v_{it} - u_{it},\quad v \sim N(0,\sigma_v^2),\ u \sim N^+(0,\sigma_u^2),

with \alpha_i an unrestricted individual effect; they differ in how they estimate it.

"TFE" is Greene's (2005, JoE/JPA) true fixed effects estimator: the ordinary composed-error likelihood with the N individual intercepts estimated jointly with (\lambda, \sigma, \beta). This package concentrates the \alpha_i out at every likelihood evaluation (each is the unique root of a strictly decreasing score, solved by safeguarded Newton to machine precision), so the reported parameter vector, its standard errors and the out layout are the same 2 + K as "TFE_WMLE". It is the estimator subject to the incidental parameters problem: with T fixed, \sigma_u is biased upward and efficiency scores correspondingly pessimistic.

"TFE_WMLE" is Chen, Schmidt and Wang's (2014, JoE) within maximum-likelihood estimator, which applies maximum likelihood to the deviations from individual means. The \alpha_i drop out of the likelihood entirely, so the estimator is free of the incidental parameters problem; the firm effects are recovered afterwards from the mean residual with the usual \sqrt{2/\pi}\,\sigma_u adjustment.

The two are not nested and their log-likelihoods are not on the same scale – "TFE" evaluates the likelihood of the data, "TFE_WMLE" that of the within-transformed deviations – so logLik(), AIC() and BIC() must not be used to choose between them.

Through sfa 1.1.2, model_name = "TFE" selected the Chen-Schmidt-Wang estimator. It now selects Greene's, matching standard usage; the Chen-Schmidt-Wang estimator is "TFE_WMLE". Existing scripts that pass "TFE" therefore get a different estimator than they did, and are warned once per call.

Degeneracy of the Greene likelihood. Because \alpha_i is unrestricted, Greene's likelihood always has a supremum on the \sigma_v \to 0 boundary: setting \alpha_i = \max_t (y_{it} - x_{it}'\beta) makes every composed error weakly negative, \Phi(-\lambda e/\sigma) \to 1, and what remains is the deterministic-frontier likelihood, which is finite and attained. On simulated panels this supremum is reached along a path that first dips below the interior local maximum and then rises past it, so an unconstrained search that escapes the interior basin returns \sigma_v \approx 0 together with meaningless standard errors; when \lambda is large and T small, no interior maximum need exist at all. tfe_lambda_max bounds the search accordingly, and a fit that pins at the bound warns. This is a property of the estimator, not of the implementation, and is one of the motivations for "TFE_WMLE".

Both problems are small-T phenomena. Simulations in which T grows with N fixed at 100 (truth \lambda = 3.333) recover \hat\lambda = 88.5, 4.8, 3.6, 3.6, 3.5 at T = 10, 20, 30, 40, 50: the fit is pinned at the bound at T = 10 and has converged by T = 50. The degeneracy and the incidental-parameters bias dissolve together as T grows, so tfe_lambda_max matters on short panels and is irrelevant on long ones. "TFE_WMLE", having no incidental parameters, is unaffected by T either way.

"SSFE" is the classical Schmidt and Sickles (1984, JBES) fixed-effects (LSDV) estimator: a standard within/fixed-effects panel regression, with technical inefficiency read directly off the estimated firm effects (each firm's distance from the best-performing firm in the sample) rather than estimated by maximum likelihood. It requires no numerical optimization and returns no log-likelihood – logLik()/AIC()/BIC() are not defined for it (see sfareg-methods). Note this is a different model from "TFE"/"TFE_WMLE" above, despite the similar names: both of those are maximum-likelihood estimators that allow a separate stochastic noise term on top of the fixed effects, while "SSFE" attributes all within-firm deviation to inefficiency directly.

"K1990" and "K1990modified" (Kumbhakar, 1990) extend the same error-components frontier with richer time paths for inefficiency: "K1990" uses B_{it} = (1 + \exp(b t + c t^2))^{-1} and "K1990modified" uses B_{it} = 1 + d(t - T_i) + e(t - T_i)^2, against "BC92"'s \exp(-\eta(t - T_{ref})) and "PL80"'s time-invariant B_{it} = 1. All four share one likelihood, differing only in B_{it}. Note that "K1990"'s b and c are weakly identified – quite different pairs produce near-identical fitted time paths – so the fitted B_{it} path is more interpretable than either coefficient on its own.

"PL80" (time-invariant panel inefficiency, Pitt and Lee, 1980, JoE) and "BC92" (time-varying decay, Battese and Coelli, 1992, JPA) are both natively estimated maximum-likelihood implementations of the same "error components frontier": y_it = x_it'beta + v_it - B_it*u_i (production; sign flips for cost), with v_it iid normal, u_i iid half-normal and constant across time within a firm, and B_it = 1 for "PL80" or B_it = exp(-eta*(t - Tref)) for "BC92" (Tref being the last time period in the whole panel). The closed-form log-likelihood is obtained by integrating u_i out of the joint density of (epsilon_i, u_i); it nests the ordinary cross-sectional normal-half-normal likelihood when every firm has a single observation. Verified against frontier::sfa() (matching coefficients, log-likelihood, and frontier::efficiencies()-style predicted efficiencies to numerical precision across balanced/unbalanced panels and production/cost specifications) before that dependency was removed. The time argument above supplies each observation's period for "BC92"'s decay term; if left NULL, a sequential within-individual time index is constructed from row order, which assumes data is already sorted chronologically within each individual.

Note on default iteration counts: psfm()'s defaults were raised in version 1.1.3, from 100/10/10 for maxit.bobyqa/maxit.psoptim/maxit.optim to 5000/100/1000, and maxit.nlminb (previously hard-coded) is now an argument. The old values were binding rather than merely economical: the K1990/K1990modified branch already carried a note that 100 bobyqa evaluations left its seven-parameter fits several log-likelihood units short of the optimum purely on the iteration cap, and GTRE_FML at N = 500, T = 10 roughly halves its root-mean-square error against known true values when the caps are lifted, at about 1.5 times the run time. They remain below sfm()'s defaults because several psfm() models (GTRE, GTRE_Z, TRE, TRE_Z) are estimated by simulated maximum likelihood over Halton draws, so each objective evaluation is far more expensive than the closed-form likelihoods used elsewhere. Lower them explicitly if a fast approximate call is wanted.

Value

An object of class "sfareg" containing components that vary by model. All models return:

Additional model-specific components:

For GTRE_FML model:

For GTRE and GTRE_Z models:

For GTRE_Z model:

For GTRE, GTRE_Z, TRE and TRE_Z models:

For TFE and TFE_WMLE models:

For TFE model only:

For FD model:

For GTRE_SEQ1 and GTRE_SEQ2 models:

When a persistent scale collapses to zero. GTRE carries two persistent components – \sigma_h, persistent inefficiency, and \sigma_r, the firm random effect – and the likelihood cannot always separate them in a given sample. When it cannot, it merges them: one scale goes to zero and the other absorbs its variation, coming back inflated.

This is common rather than exceptional. On simulated data at T = 6 with both components genuinely present (\sigma_r = 0.20, \sigma_h = 0.40), one of the two collapsed in 37% of 87 replications, and the commoner direction was \sigma_r \to 0 rather than \sigma_h \to 0. The share with both interior does improve with N (15, 19 and 21 of about 29 at N = 50, 100, 200) but slowly.

A collapse is usually the right answer, not a numerical failure. On one such replication – whose realized \sigma_h scale was 0.3695, so persistent inefficiency was certainly present – the boundary solution had a log-likelihood higher than the true parameter vector by 3.66, a likelihood ratio of about 39. This is the panel counterpart of the cross-sectional wrong-skewness result (Waldman 1982), where \sigma_u = 0 is the correct maximum likelihood estimate rather than an optimizer failure; see skewness_test. This is not in tension with the start search described above under "GTRE_FML": a boundary fit that beats the truth on likelihood is the estimate, while one that loses to it is a trapped optimizer, and the search exists to remove the second so that what remains is the first.

What a collapse does to the rest of the fit, exactly. One rule covers it: a collapse conserves the persistent variance and relocates the persistent mean into the intercept. Those are two different destinations. h is half-normal, so it carries both a mean E[h] = \sigma_h\sqrt{2/\pi} and a variance \mathrm{Var}(h) = \sigma_h^2(1 - 2/\pi), and the model has a separate place for each: the frontier intercept takes the mean, since E[y] = \beta_0 - E[h] - E[u], and the surviving persistent scale takes the variance. The intuitive shortcut – that the survivor comes back at \sqrt{\sigma_h^2 + \sigma_r^2} – ignores that split and is wrong.

At \sigma_h = 0.40, \sigma_r = 0.20 that predicts 0.3133 and -0.3192 in the first row, 0.5197 and +0.0955 in the second, with nothing fitted. Over 1000 replications at each of five sample sizes the measurements are 0.2973, -0.3035, 0.5058 and +0.0931 – within 2.5% to 5.1%, and all four slightly below prediction, which is what selection produces rather than a wrong constant.

The practical consequence is that a collapsed fit does not merely mislabel the persistent split. It displaces the intercept by a known and substantial amount, so any level read off coef() is displaced with it.

What a collapse is telling you about your data. It is wrong skewness, measured between firms rather than within them. Take the OLS residuals, average them within each firm – that average is the persistent part of the error – and look at their skew. Over 1000 replications, the samples whose between-firm skew came out positive (the wrong way for a production frontier) lost \sigma_h to the boundary 69% of the time, against 3.6% for the rest; the area under the ROC curve for that single statistic is 0.955. Not one wrong-skew sample collapsed \sigma_r instead. So a collapsed \sigma_h is not a numerical accident to be retried with different starting values: this sample carries no between-firm evidence of persistent inefficiency, and the estimate is reporting that.

Note that this is a reading of the sign, not a test. The firm mean of the residuals contains -\bar u_i as well as -h_i, and the mean of T half-normals is itself skewed the same way, so on short panels the between-firm skew leans negative even when \sigma_h = 0 – 82% of the time at T = 4, 63% at T = 10, 53% at T = 30. skewness_test is therefore not applicable to the between-firm residuals as it stands.

Whether firms or periods help depends on the estimator, and the two answers are opposite. Crossed (N, T) designs on the same data-generating process, fitting the collapse rate on \log N and \log T:

Under the default, firms buy about 1.8 times what periods do; under simulated ML it is the other way round, periods buying about 1.6 times what firms do. Restricting the "fiml" design to N \le 200, cell for cell the same grid the "sml" one used, leaves it at -1.64 (0.22), so this is the estimator and not the design; the difference in the N coefficient is significant (z = 2.8) and the difference in the T coefficient is not.

The practical reading, for the estimator you get by default: add firms. At T = 4 the collapse rate falls 68%, 36%, 18%, 6% across N = 50, 100, 200, 400, and reaches 0% at N = 400, T = 16. Under "sml" the same row is flat – 60%, 39%, 41% – which is why a short panel really is a dead end there however many firms it has. If you must use estimator = "sml" (an unbalanced panel forces it), then length is what buys the split, and a four-period panel will not give it to you.

How to tell in advance whether your panel can support the split. One quantity governs it, and it is not \sigma_h/\sigma_r. What the likelihood has to do is read the persistent part out of firm means that the transient block contaminates, so write the firm-mean composite as r_i - h_i - \bar u_i + \bar v_i and compare its persistent variance with the transient variance left in it:

S/N = \frac{\sigma_r^2 + \sigma_h^2(1 - 2/\pi)}{[\sigma_u^2(1 - 2/\pi) + \sigma_v^2]/T}

Across six designs varying all four scales in different combinations (50 replications each, N = 100, T = 6), this orders the collapse rate perfectly – Spearman -1.000 under both estimators, and a logistic coefficient on \log(S/N) of -1.15 (0.20):

As a rule of thumb: below S/N \approx 1 about half your fits will merge the two components, and it takes roughly S/N > 5 to get under 10%.

Two things follow that are easy to get backwards. First, raising \sigma_h/\sigma_r does not help – it only decides which scale dies. Going from a ratio of 2 to 4 by shrinking \sigma_r takes \sigma_h \to 0 down (0.15 to 0.10) and \sigma_r \to 0 up (0.17 to 0.40), and the net rate rises from 31% to 50%, because shrinking \sigma_r removes persistent variance and lowers S/N from 1.30 to 0.90. Holding the ratio at 2 and scaling both components instead takes the rate from 31% to 8%. Second, more transient inefficiency hurts: \sigma_u enters the denominator, so a noisier transient block contaminates the firm means and makes the persistent split harder.

Note what this quantity does and does not explain. T enters it directly, so it covers the contamination channel by which longer panels help. It says nothing about N, which helps through a different route – sampling error in what is effectively a deconvolution from N firm-level draws – which is why the two dimensions have to be measured separately, as above.

One limit on all of that, stated because it is easy to over-read. Those elasticities were measured at \sigma_r = 0.20 against \sigma_h = 0.40, where persistent inefficiency is the larger of the two time-invariant components. On a design where the random effect dominates instead (\sigma_r = 0.50, \sigma_h = 0.40, N = 100, T = 6) the level ordering reverses: "fiml" collapses on 33% of replications and "sml" on 20%. Which estimator handles the split better is not a fact about the estimators alone – it depends on how much persistent inefficiency there is to find.

psfm therefore reports the condition rather than preventing it, in $sigh_at_bound and $sigr_at_bound and in a warning naming the surviving scale. Bounding either parameter away from zero would corrupt precisely the samples where the boundary is the answer. What the warning asks of you is to read the two persistent scales together and to treat the split between them as unidentified in that sample – their sum remains informative even where the division does not.

Profiling a fitted panel likelihood. With keep_objective = TRUE the fitted object carries $objective, the same negative-summed-log-likelihood the optimizer minimised. Evaluating it away from the optimum is the only way to ask whether a parameter is weakly identified, as distinct from badly estimated:

fit <- psfm(y ~ x1 + x2, "GTRE", d, individual = "id",
            estimator = "sml", keep_objective = TRUE)
ih   <- which(names(coef(fit)) == "sigh")
grid <- seq(0.05, 1, by = 0.05)
ll   <- vapply(grid, function(v) {
          p <- fit$opt$par; p[ih] <- v; -fit$objective(p) }, numeric(1))
plot(grid, ll - max(ll), type = "l")

A flat curve means no estimator could recover the parameter from these data; a sharp one means the difficulty lies elsewhere. For GTRE's sigh the curve spans over 200 log units at N = 200 and sharpens with N, which is what ruled out weak identification as the explanation for its convergence behaviour.

The Schmidt–Sickles family. "SSFE", "SSRE" and "SSCRE" are the within, GLS and correlated-random-effects members of one family: all three read inefficiency off a firm effect as distance from the best firm, and none assumes a distribution for u. The choice among them is the standard panel trade-off, with a frontier-specific twist. Random effects is the more efficient estimator and is the only one of the three that identifies coefficients on time-invariant regressors, which the within transformation sweeps away – but it requires the effects to be uncorrelated with the regressors, and in a frontier model that assumption says a firm's inefficiency is unrelated to its input choices, which is close to saying managers do not respond to their own productivity.

"SSCRE" is the answer to that. Adding the within-firm means of the time-varying regressors (Mundlak 1978) models the correlation instead of assuming it away; the slopes on the original regressors then equal the within estimates exactly – not approximately – while time-invariant regressors remain identified. The coefficients on the added means are reported with a .mean_ prefix, and testing them jointly against zero is the Mundlak form of the Hausman test: if they are indistinguishable from zero, "SSRE" was safe. A regressor already constant within firm has a group mean identical to itself, so it contributes no Mundlak term; those are exactly the regressors random effects exists to identify, and they are left alone.

Kneip, Sickles and Song (2012). "KSS" writes the firm effect as \alpha_{it} = \sum_{r=1}^{L}\theta_{ir} g_r(t), where both the basis functions g_r and the dimension L are estimated from the data rather than assumed. That makes explicit what the other time-varying estimators here are special cases of:

Estimation follows the three steps of the original: cross-sectional centering by period (which imposes \sum_i u_i(t) = 0 and removes any common function of time), then smoothing of each firm's residual trajectory followed by an eigendecomposition of their empirical covariance to give \hat g_r, then loadings by least squares and the effects rebuilt from them – iterated with \beta until the residual sum of squares settles. Smoothing before the eigendecomposition is what makes the estimated basis a set of smooth functions rather than T arbitrary period effects, and is the whole difference between this and an unrestricted factor model.

Two tuning choices are therefore unavoidable, and both are exposed. kss_smooth defaults to generalized cross-validation. kss_L defaults to the Bai and Ng (2002) IC_{p2} criterion – not the threshold rule of the original paper, which is stated here because it matters: IC_{p2} is the standard criterion for this problem and the literature that motivates the estimator treats Bai and Ng (2002) and Bai (2009) as the comparable methods, but a different rule can select a different L. Set kss_L explicitly to take the choice out of the criterion's hands, and inspect $kss$eigenvalues to see whether the choice was close.

Mundlak adjustment terms (mundlak). The true random effects model assumes the firm effect is independent of the regressors. In production data it usually is not – better managers choose different input levels – and the effect is then split between the slopes and inefficiency, biasing the technology parameters and understating inefficiency. Mundlak's (1978) device models the correlation instead of assuming it away, \alpha_i = \pi'\bar{x}_i + \delta_i. Substituting into the frontier leaves the same likelihood with a wider design matrix, which is why this is an argument rather than a new model_name.

Karagiannis and Kellermann (2019) widen it to \alpha_i = \pi'\bar{x}^*_i + \gamma'z_i + \delta_i, adding environmental variables – conditions outside the producer's control, such as agro-ecological conditions or the geography an airline serves. Those enter the formula as ordinary time-invariant regressors; mundlak supplies the piece a user cannot easily build by hand, the within-firm means computed against the same panel index the fit uses. Naming a subset of the regressors gives their Hausman-Taylor variant, in which only some regressors are taken to correlate with the effect.

Measured over 40 replications with N = 80, T = 6, on a design where the firm effect is correlated with \bar{x} by construction, and on a control where it is not:

So it removes the bias where there is one and nearly halves RMSE – and where there is no correlation it costs about 23% in RMSE for nothing. It is not free, and it is not a default.

A variable that is already constant within firm gets no mean column: its own mean is itself, so the two would be exactly collinear. Every mundlak fit also emits a collinearity warning from the starting-value regression, which is expected and harmless – in a between-individual regression x and \bar{x} are the same column by construction. The starting values fall back to pooled OLS and the likelihood is estimated with the full requested formula.

The criterion needs a long enough panel, and says so. IC_{p2} works because the residual variance flattens out once the real factors are in and the penalty then dominates. On a short panel that never happens: if L is allowed up to T-1 the factors can span nearly the whole time dimension, the variance collapses toward zero, and no penalty competes. Measured at n = 100 over five seeds and two true ranks, an uncapped criterion returned its maximum on all ten designs at T = 6 and again at T = 8. The automatic search is therefore capped at \lfloor T/2 \rfloor, which recovers the true rank exactly from T = 8 upward; at T = 6 it is right for a one-factor structure and occasionally one too high for a two-factor one. Whenever the automatic choice lands on that cap psfm warns, because a selection at the boundary is not a free choice – read it as "at least this many", inspect the eigenvalues, and prefer a longer panel to a larger kss_L_max.

One consequence of the period centering is worth knowing before comparing models: "KSS"'s alpha_hat carries no period mean, while "CSS"'s and "LS"'s do. The two therefore differ by a per-period constant by construction and should not be correlated against each other directly. u_hat is a within-period contrast and is invariant to that shift, so it is the quantity to compare across these models – and is what the estimators are for.

Time-varying fixed effects without a distribution for u. "CSS" and "LS" sit between "SSFE", which holds inefficiency fixed over time, and the maximum-likelihood panel models, which let it move but only by assuming a distribution for u and a decay path for it. Both instead let the firm effect itself vary with time and read inefficiency off it, so neither assumes any distribution and neither returns a log-likelihood.

They differ in what they are willing to spend on flexibility. Cornwell, Schmidt and Sickles (1990) give every firm its own quadratic in time, \alpha_{it} = \theta_{i1} + \theta_{i2}t + \theta_{i3}t^2: 3N parameters, firms may overtake one another, and T_i \ge 4 is required before a firm contributes to \beta at all (a firm observed three times is fitted exactly by its own quadratic and its residuals are identically zero – sfa warns rather than letting such firms silently dilute the estimate). Lee and Schmidt (1993) impose one common temporal pattern scaled per firm, \alpha_{it} = \delta_t\alpha_i: N + T - 1 parameters, and the ranking of firms cannot change, only the spread. Cross-over versus no cross-over is the substantive choice between them, not a technical one.

"CSS" is fitted by residualizing y and X on each firm's own (1, t, t^2) rather than by building 3N dummy columns; a regressor spanned by those quadratics – a pure time trend, for instance – is a rank failure and is refused rather than reported as a small coefficient. "LS" is a rank-one factor model, so it is fitted by alternating least squares, which handles an unbalanced panel with no special-casing; \delta_1 = 1 fixes the scale that \delta_t\alpha_i leaves free, with a fallback to a norm constraint when the first period's effect is near zero.

For both, inefficiency is identified only relative to the best firm in each period – u_{it} = \max_j \alpha_{jt} - \alpha_{it} – which is what makes them time-varying at all; comparing against a single overall best would collapse them back to a time-invariant ranking. Neither absorbs an intercept separately: it goes into the firm effect, as in "SSFE", so coef() returns slopes only.

For SSFE model:

For the GTRE model:

For CSS and LS models:

For SSRE and SSCRE models:

For the KSS model:

For PL80 and BC92 models:

Note

Standard errors require optHessian set to TRUE

Note

The GTRE_SEQ1 and GTRE_SEQ2 models use sequential estimation methods and do not return optimization objects or starting values. The SSFE, CSS, LS, SSRE, SSCRE and KSS models are likewise not maximum-likelihood based and return no optimization object, starting values, or log-likelihood. All panel models require the individual argument to identify panel units.

Author(s)

David H. Bernstein

References

Fillipini and Greene (2016, JPA); Colombi, Kumbhakar, Martini and Vittadini (2014, JPA) on the closed-skew-normal four-component likelihood, with Colombi (2010) and Colombi, Martini and Vittadini (2011) on the two-step estimator used to seed it; Greene (2005, JoE and JPA); Wang and Ho (2010, JoE); Chen, Schmidt and Wang (2014, JoE); Kumbhakar (1990, JoE); Pitt and Lee (1980, JoE); Battese and Coelli (1992, JPA); Battese and Coelli (1988, JoE); Jondrow, Lovell, Materov and Schmidt (1982, JoE); Schmidt and Sickles (1984, JBES)

See Also

sfm for cross-sectional models, psfm_bootstrap for bootstrap standard errors, data_gen_p for simulating panel data with known true parameters, and sfareg-methods for extracting results.

Examples

library(sfa)     

## A deliberately small panel with few simulation draws, so the examples
## build quickly. Do not read the estimates as a serious fit. TRE_Z becomes
## unstable below about 60 firms, so this is close to the floor.
data_trial <- data_gen_p(t=6, N=60, rand = 100,
                         sig_u = 1,  sig_v = 0.3,
                         sig_r = .2, sig_h = .4,
                         cons = 0.5, beta1 = 0.5,
                         beta2 = 0.5)

## rand.psoptim is not optional with PSopt = TRUE: the particle swarm draws
## from the session RNG, so results otherwise change between builds.
max_tre_z   <-  psfm(formula    = y_tre_z ~ x1 +x2| z_gtre,
                     model_name = "TRE_Z",
                     data       = data_trial,
                     individual = "name",
                     halton_num = 50,
                     rand.psoptim = 7L,
                     PSopt      = TRUE)

## Time-varying fixed effects, no distribution assumed for u. CSS gives each
## firm its own quadratic in time; LS gives all firms one common pattern.
css_fit     <-  psfm(formula = y_tre ~ x1 + x2, model_name = "CSS",
                     data = data_trial, individual = "name", time = "year")
ls_fit      <-  psfm(formula = y_tre ~ x1 + x2, model_name = "LS",
                     data = data_trial, individual = "name", time = "year")
ls_fit$delta      # the common temporal pattern, normalized to delta[1] = 1

## The Schmidt-Sickles family: within, GLS, and correlated random effects.
## SSCRE's slopes equal SSFE's within slopes exactly (Mundlak 1978); the
## .mean_ coefficients are the Mundlak form of the Hausman test.
ssre_fit    <-  psfm(formula = y_tre ~ x1 + x2, model_name = "SSRE",
                     data = data_trial, individual = "name", time = "year")
sscre_fit   <-  psfm(formula = y_tre ~ x1 + x2, model_name = "SSCRE",
                     data = data_trial, individual = "name", time = "year")

## Kneip, Sickles and Song (2012): the temporal basis and its dimension are
## both estimated. Inspect what the criterion chose before trusting it.
kss_fit     <-  psfm(formula = y_tre ~ x1 + x2, model_name = "KSS",
                     data = data_trial, individual = "name", time = "year")
kss_fit$kss$L                 # selected dimension
kss_fit$kss$eigenvalues       # was the choice close?

Parametric Bootstrap for psfm() Panel Stochastic Frontier Models

Description

Performs a parametric bootstrap for panel stochastic frontier models fit with psfm. Supports model_name %in% c("GTRE_Z", "TRE_Z", "GTRE", "GTRE_FML", "TRE", "TFE", "FD") – every psfm() model except the moment-based/LSDV estimators ("GTRE_SEQ1", "GTRE_SEQ2", "SSFE", which are not maximum-likelihood and have no $opt to bootstrap around in the same sense) and the sfa-wrapped models ("PL80", "BC92", whose $out/efficiency-score fields don't follow this package's own conventions). Each supported model's frontier/variance-component parameter blocks are read directly from psfm_object$out and psfm_object$formula; for each of BOOT replications, a new response is simulated from the fitted parameters using a data-generating process appropriate to that specific model (see Details), the model is re-estimated in parallel, and bootstrap standard errors / t-values are returned for every parameter in psfm_object$out.

Usage

psfm_bootstrap(psfm_object,
               numCores,
               BOOT,
               individual,
               h_type        = c("auto", "none", "scalar", "parametric"),
               maxit.psoptim = 1000,
               seed_offset   = 0,
               write_back    = TRUE,
               pkgs          = c("sfa", "Formula", "pbapply", "truncnorm"),
               inefdec,
               rand.gtre     = NULL,
               rand.psoptim  = NULL,
               maxit.bobyqa  = 1,
               maxit.optim   = 1)

Arguments

Details

Supported models fall into three families, each with its own data-generating process and $out row layout (see R/psfm_bootstrap.R's header comment for the full derivation, including how each was cross-checked against that model's own likelihood in psfm.R):

"randeff" family (GTRE_Z, TRE_Z, GTRE, TRE)

y* = X %*% beta_hat + v + u + r + h (sign per inefdec), with v ~ N(0, sigma_v^2), u half-normal (covariate-driven sigma for GTRE_Z/TRE_Z, a single homoskedastic scalar for bare GTRE/TRE, which have no z pipe at all), r ~ N(0, sigma_r^2) (one draw per individual, repeated over time), and h per h_type (GTRE/GTRE_Z only). $out row layout: sigv, sigr, [x-block], [z-block], [h-block] for GTRE_Z/TRE_Z; lambda, sigma, sigr, [sigh], [x-block] for bare GTRE/TRE (note the different reparameterization – see psfm's own model-family sigma-link inconsistency note).

"tfe" family (TFE)

Fixed-effects model: the individual effects are held FIXED at psfm_object$r_hat_m's original point estimates across every replication (not redrawn – redrawing them would contradict treating them as fixed rather than random), and only v ~ N(0, sigma_v^2) and u ~ |N(0, sigma_u^2)| (both derived from the fitted lambda/sig, or gamma/sig if the original fit used psfm(..., gamma = TRUE)) are redrawn. $out row layout: lambda (or gamma), sig, [x-block] (no intercept row – the within-transformation removes it).

"fd" family (FD)

Wang and Ho (2010) first-difference estimator: a single time-invariant u_i ~ TruncatedNormal(mean = mu, sd = sqrt(sig_u2), lower = 0) draw per individual (via rtruncnorm), scaled by a deterministic time-varying factor h_it = exp(z_it %*% delta), plus v_it ~ N(0, sig_v2). $out row layout: sig_u2, sig_v2, mu, [x-block], [z-block] (no intercept row in either the x- or z-block).

psfm_bootstrap validates that nrow(psfm_object$out) matches the expected layout for the resolved model/h_type before running any bootstrap replications, and stops with an informative error if it does not. For GTRE/GTRE_Z it similarly validates up front that length(psfm_object$H) equals the number of unique individuals implied by the individual column ($H is one value per individual, not one per observation) – TRE/TRE_Z/TFE/FD do not return $H at all, so boot_eff_h is simply NULL for those four models (not an error).

The model is re-estimated on each simulated dataset via sfa::psfm(), executed in parallel across replications using parLapply (or pblapply for a live progress bar, when the pbapply package is installed). Whether the psoptim (particle-swarm) optimizer stage is used on each refit is model-dependent: it is enabled for the "randeff" family (matching this function's original behavior, tested clean), but disabled for "TFE" – testing found psoptim's wide exploration made TFE repeatedly collapse onto a degenerate variance-boundary solution (reproduced even refitting the SAME real data across different psoptim seeds), a problem that disappeared entirely once that stage was skipped (matching psfm()'s own default of PSopt = FALSE). As a defense-in-depth safety net across every family, any replication whose refit still lands on a degenerate variance-boundary mode (a scale parameter pinned near its numerical floor – the same well-understood boundary-mode phenomenon documented for ttsfm()'s TTNE/TTHN models) is treated as a per-replication FAILURE rather than silently corrupting the bootstrap standard error with a meaningless draw.

Re-estimation failures (optimizer errors, or the degenerate-mode check just described) are caught per replication, recorded with their message, and set to NA in the output matrices rather than stopping the bootstrap; a single warning at the end lists every failed replication index together with its message. Bootstrap standard errors are computed with na.rm = TRUE, so they remain well-defined even when some replications failed.

Requires the parallel and Formula packages (and truncnorm when bootstrapping an "FD" fit), and requires sfa (and any other packages named in pkgs) to be installed and loadable on each cluster worker.

Value

A list with components:

Author(s)

David H. Bernstein and Christopher F. Parmeter

See Also

psfm, parLapply, pblapply, Formula

Examples

## Not run: 
## Not run: a meaningful bootstrap needs BOOT in the hundreds, each
## replication refitting the model on a PSOCK cluster, so this takes
## minutes rather than seconds.

data_trial <- data_gen_p(t = 5, N = 50, 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)

GTRE_ZZ <- psfm(y_gtre_zz ~ x1 + x2 | z_gtre | zp_gtre,
                model_name = "GTRE_Z",
                data       = data_trial,
                individual = "name")

boot_results <- psfm_bootstrap(
  psfm_object = GTRE_ZZ,
  numCores    = 2,
  BOOT        = 199,
  individual  = "name",
  h_type      = "auto",
  inefdec     = FALSE
)

## Bootstrap standard errors / t-values for every model parameter
boot_results$se
boot_results$tval

## Bootstrap draws of time-varying (U) and time-invariant (H)
## efficiency/inefficiency scores, one row per replication
boot_results$boot_eff
boot_results$boot_eff_h

## Fitted object with bootstrap SEs/t-values written into $out
boot_results$model$out

## End(Not run)

Stochastic Frontier Models with Sample Selection or Technology Choice

Description

Two models, chosen by model_name.

"greene" (the default) fits Greene's (2010) sample-selection stochastic frontier: a normal–half normal frontier observed only for the subset of units that select into the sample, where the unobservables driving selection are correlated with the noise in the frontier. Estimation follows the two-step scheme of the paper – a probit for the selection equation on the full sample, then maximum simulated likelihood for the frontier on the selected subsample.

"kts" fits Kumbhakar, Tsionas and Sipilainen (2009), which is a different model rather than a variant. Both technologies are observed, each with its own frontier and its own two scales, and the technology choice depends on inefficiency itself: u lowers output through the frontier and shifts the adoption decision through the choice equation. Fitted by single-step maximum likelihood over the whole sample.

Usage

selsfm(selection, frontier, data, model_name = c("greene", "kts"),
       n_nodes = 64, Nsim = "auto",
       sim_type = c("halton", "sobol", "torus", "uniform"),
       antithetics = FALSE, seed = NULL, inefdec = TRUE,
       maxit.bobyqa = 10000, maxit.psoptim = 1000, maxit.optim = 1000,
       start_val = FALSE, PSopt = FALSE, optHessian = TRUE,
       Method = "L-BFGS-B", verbose = FALSE, rand.psoptim = NULL)

Arguments

Details

The model is Greene's Equation (10):

d_i = 1[\alpha'z_i + w_i > 0], \quad w_i \sim N(0,1)

y_i = \beta'x_i + v_i - u_i, \quad \mathrm{observed\ only\ when\ } d_i = 1

u_i = \sigma_U|U_i|, \quad U_i \sim N(0,1), \qquad v_i = \sigma_V V_i, \quad V_i \sim N(0,1)

(w_i, v_i) \sim N_2[(0,0), (1, \rho\sigma_V; \rho\sigma_V, \sigma_V^2)]

Selection operates through the noise v, not through the inefficiency u. That is the substantive assumption separating this model from Kumbhakar, Tsionas and Sipilainen (2009), where the selection disturbance is correlated with u instead; the two answer different questions and are not nested.

Because u cannot be integrated out in closed form, the log-likelihood is simulated (Greene's Equation (15)): for each selected observation R draws of |U_{ir}| are taken and

\log L = \sum_{d_i=1} \log \frac{1}{R} \sum_{r=1}^{R} \phi\!\left(\frac{v_{ir}}{\sigma_V}\right)\frac{1}{\sigma_V} \Phi\!\left(\frac{\rho v_{ir}/\sigma_V + \alpha'z_i}{\sqrt{1-\rho^2}}\right)

with v_{ir} = y_i - \beta'x_i + \sigma_U|U_{ir}|. Everything is evaluated in logs and accumulated by log-sum-exp, since the summand underflows in the tails well before the total does.

A note on the published formula. Equations (11) and (13)–(15) of Greene (2010) write \sigma_\varepsilon in the denominator of the \Phi argument. It must be \sigma_V: the term being standardized is v_{ir} = y_i - \beta'x_i + \sigma_U|U_{ir}|, which is the noise, and the conditional distribution of w given v for a bivariate normal with \mathrm{Var}(w) = 1 and correlation \rho is N(\rho v/\sigma_V,\, 1-\rho^2). This implementation uses \sigma_V.

Standard errors are conditional on the first stage. Following the paper's two-step scheme, \alpha is taken as given in the frontier likelihood. The reported standard errors come from the second-stage Hessian and therefore do not include the sampling variability of \hat\alpha; Greene applies the Murphy–Topel correction at this point, which is not implemented here. Treat them as a lower bound, and bootstrap over both stages when the precision matters.

Value

An object of class "sfareg". As elsewhere in the package out is a p x 3 matrix, one row per parameter, with columns par, st_err and t-val; the parameters are ordered sigma_u, sigma_v, rho, then the frontier coefficients. Alongside the usual components:

The KTS model

For technology I_i \in \{0,1\},

y_i \mid I_i, u_i \sim N(x_i'\beta_{I_i} - u_i,\ \sigma_{v,I_i}^2), \qquad P(I_i = 1 \mid u_i) = \Phi(z_i'\gamma + \delta u_i), \qquad u_i \mid I_i \sim N^+(0, \sigma_{u,I_i}^2).

\delta is the parameter the model exists for: it is the dual role of inefficiency, and \delta > 0 means the less efficient are more likely to choose technology 1. Neither two-step order works – the choice equation cannot be a probit because u is unobserved, and fitting the frontiers first ignores the endogeneity of the choice – which is why the paper estimates in one step.

Because u's distribution depends on I while I's probability depends on u, the marginal density of u is a mixture of the two regimes' half-normals, with weight \phi_i solving the paper's equation (9). The observation's density is then the same mixture of two one-dimensional integrals, its equation (14).

Recovery. On data simulated from this exact model with \delta = 0.9, the frontier coefficients, the four scales and \delta all recover, and improve with the sample:

\delta is the parameter that needs the most data, as one would expect of a coefficient identified only through the coupling between the two equations, but it is properly identified: holding the rest at the truth, the mean log-likelihood falls monotonically away from \delta = 0.9 in both directions, by about 0.069 per observation across [0.1, 1.8].

References

Kumbhakar, S.C., Tsionas, E.G. and Sipilainen, T. (2009). Joint estimation of technology choice and technical efficiency: an application to organic and conventional dairy farming. Journal of Productivity Analysis, 31(3), 151–161.

Greene, W. (2010). A stochastic frontier model with correction for sample selection. Journal of Productivity Analysis, 34(1), 15–24.

Heckman, J. (1979). Sample selection bias as a specification error. Econometrica, 47(1), 153–161.

See Also

sfm for the frontier without a selection correction.

Examples

set.seed(1)
n  <- 800
z1 <- rnorm(n); x1 <- rnorm(n)
w  <- rnorm(n)
d  <- as.numeric(0.2 + 0.8 * z1 + w > 0)
v  <- 0.5 * (0.6 * w + sqrt(1 - 0.6^2) * rnorm(n))
y  <- 0.5 + 0.8 * x1 + v - abs(rnorm(n))
y[d == 0] <- NA
dat <- data.frame(y = y, x1 = x1, z1 = z1, d = d)

fit <- selsfm(selection = d ~ z1, frontier = y ~ x1, data = dat, Nsim = 100)
fit$out

Optimizer Diagnostics for Stochastic Frontier Fits

Description

Reports how the optimizer actually finished and how well the likelihood identifies the parameters: the convergence code and message, the eigenvalue spectrum and conditioning of the Hessian, the implied parameter correlations, and – where the objective was retained – the gradient at the reported optimum and a likelihood slice through each parameter.

Usage

sfa_diagnostics(object, ...)

## S3 method for class 'sfareg'
plot(x, which = 1:4, n_grid = 41, span = 0.25, ...)

## S3 method for class 'sfadiag'
print(x, ...)

Arguments

Details

Every stochastic frontier fit in this package runs a staged minimizer and stores optim()'s convergence code, message, evaluation counts and numerically differentiated Hessian. None of that was reported: a fit that stopped on the iteration cap printed exactly like a converged one. sfa_diagnostics() is the reporting side of that.

What the Hessian tells you. optim() returns the Hessian of the negative log-likelihood, so at a genuine interior optimum it is positive definite. A non-positive eigenvalue means the reported point is not a minimum in that direction. A very small one means the likelihood is nearly flat along it, which is what weak identification looks like numerically – $hessian$flat_direction names the parameters that load on that flattest axis, and they are usually the pair the data cannot separate.

Slices, not profiles. Panel 3 varies one parameter and holds the others at their estimates. A true profile would re-optimize the remaining parameters at every grid point, at the cost of one optimization per point. The slice is the right object for reading curvature at the optimum and for the simplest check there is: the estimate should sit at the peak of its own slice. If it does not, the optimizer stopped early.

Retaining the objective. The gradient and the slices need the likelihood itself, which an ordinary fit does not keep. Pass keep_objective = TRUE to sfm to store it. This is off by default because a closure carries its enclosing environment, so a fit saved with one serializes the estimation data alongside the results – on a 200-observation example the fitted object grows from roughly 38 KB to 1.7 MB.

The convergence code alone is not diagnostic, and the flags do not treat it as though it were. Across NHN, NE and NTN at n = 150, 500 and 1500, code 52 (“ABNORMAL_TERMINATION_IN_LNSRCH”) turns up routinely alongside a maximum relative gradient of about 10^{-6} and a positive definite Hessian: the staged minimizer had already found the optimum and the final L-BFGS-B stage simply could not take a step away from it. The same code on NTN at n = 150 came with a relative gradient of 5 \times 10^{7} and an indefinite Hessian, which is a genuine failure. What separates the two is the gradient and the Hessian, so the verdict is built from those and the code is reported as context:

benign

a line-search code (51/52) together with a small relative gradient and a positive definite Hessian. Reported as noise.

unverified

a line-search code with no retained objective, so there is no gradient to settle it. Reported as an open question rather than resolved in the fit's favour – refit with keep_objective = TRUE.

failure

anything else, including code 1 (the iteration limit), which is never excused: it means the search was still moving when it ran out of budget, whatever the curvature at the stopping point.

Flags are raised for a non-zero convergence code, a Hessian that is not positive definite or is worse than 10^{8} conditioned, any parameter pair correlated beyond |0.95|, and a large relative gradient. They are phrased as statements about what the fit means, not as generic warnings.

Value

sfa_diagnostics() returns an object of class "sfadiag": a list with model_name, call, pnames, estimates, convergence (code, meaning, message, counts, logLik), hessian (eigenvalues, condition, pos_def, flat_direction), correlation (matrix, worst_pair, worst_value); convergence additionally carries benign_nonzero and unverified_nonzero, gradient (gradient, max_abs, relative, max_rel; NULL without a retained objective), flags, and has_objective.

plot.sfareg() is called for its side effect and returns the "sfadiag" object invisibly. print.sfadiag() returns its argument invisibly.

Author(s)

David H. Bernstein

See Also

sfm, sfareg-methods, summary.sfareg

Examples

library(sfa)
d <- data_gen_cs(N = 300, rand = 3, sig_u = 1, sig_v = 0.3,
                 cons = 0.5, beta1 = 0.5, beta2 = 0.5, a = 5, mu = 0.5)

## Retaining the objective enables the gradient and the likelihood slices
f <- sfm(y_pcs ~ x1 + x2, model_name = "NHN", data = d, keep_objective = TRUE)
sfa_diagnostics(f)

## Four-panel diagnostic display
plot(f)

## Works without the objective too, on the Hessian alone
g <- sfm(y_pcs ~ x1 + x2, model_name = "NHN", data = d)
sfa_diagnostics(g)$hessian$condition

Extractor methods for "sfareg" objects

Description

Standard extractor methods for objects of class "sfareg", as returned by psfm, sfm, zsfm, and ttsfm. These let "sfareg" fits be used with the usual R modeling idioms (e.g. coef(fit), logLik(fit), AIC(fit), BIC(fit)) instead of reaching into the returned list by hand.

AIC() and BIC() are not defined explicitly for "sfareg" objects; they fall through to their default methods in the stats package, which only require a working logLik method (and, for BIC(), a working nobs method) to work correctly.

Usage

## S3 method for class 'sfareg'
coef(object, ...)
## S3 method for class 'sfareg'
vcov(object, type = c("hessian", "bhhh"), ...)
## S3 method for class 'sfareg'
logLik(object, ...)
## S3 method for class 'sfareg'
nobs(object, ...)

Arguments

Details

BHHH standard errors. vcov(fit, type = "bhhh") returns the inverse of the outer product of the per-observation scores instead of the inverse Hessian. Its value is that it needs no Hessian at all: it is defined whenever the scores are, which is exactly the case where the default path fails. When the Hessian is singular vcov otherwise falls back to a diagonal approximation built from the reported standard errors, discarding every covariance; BHHH keeps them. It requires keep_objective = TRUE for the same reason estfun does. It is not the default, because it is only as good as the information-matrix equality and so disagrees with the Hessian under misspecification – which is useful when the two are compared on purpose and misleading when they are not.

Robust and clustered standard errors. vcov returns the inverse Hessian, which is valid only if the likelihood is correctly specified and the observations are independent – neither being safe to assume in applied frontier work. Registering bread and estfun makes the sandwich package operate on these fits, so sandwich::sandwich() gives Huber–White errors and sandwich::vcovCL() gives clustered ones, and either can be handed to lmtest::coeftest():

fit <- sfm(y ~ x1 + x2, model_name = "NHN", data = d,
           keep_objective = TRUE)
lmtest::coeftest(fit, vcov. = sandwich::vcovCL(fit, cluster = d$firm))

Two requirements. The fit must have been made with keep_objective = TRUE, since the scores are differenced from the retained likelihood; and the robust divergence estimators are excluded, because they do not maximise a log-likelihood and so have no score. sandwich::vcovHC() does not work and is not expected to: its finite-sample corrections are built from hat values, which are not defined for a nonlinear likelihood whose parameter vector includes variance components.

The scores are taken by central differences with a per-parameter step, there being no analytic gradient for most of these models. Where one does exist (NHN) the two agree to about 3 \times 10^{-7}, and they agree with numDeriv to the same order.

coef.sfareg returns the named vector of estimated parameters (object$coefficients).

vcov.sfareg returns the estimated variance-covariance matrix. If the fit was estimated with optHessian = TRUE and the Hessian is invertible, this is solve(object$opt$hessian). Otherwise it falls back to a diagonal matrix built from the reported standard errors (object$std.errors), with a warning; if neither is available it returns a matrix of NAs with a warning.

logLik.sfareg returns the maximized log-likelihood as an object of class "logLik", with "df" and "nobs" attributes set so that AIC()/BIC() work automatically. Every estimator in this package minimizes the negative summed log-likelihood internally, so logLik = -object$opt$value. Exception: psfm()'s "GTRE_SEQ1" and "GTRE_SEQ2" models are estimated by a moment-based procedure and "SSFE" by fixed-effects (LSDV) regression, not by maximum likelihood; none of the three carries an $opt component, so calling logLik() (or, transitively, AIC()/BIC()) on those fits returns NA with a warning rather than erroring.

nobs.sfareg returns the number of observations used in estimation, recovered from object$data when the fitting function stored it, or by re-evaluating the data argument of the original call otherwise. Returns NA if neither is available.

Value

coef.sfareg returns a named numeric vector.

vcov.sfareg returns a numeric matrix.

logLik.sfareg returns an object of class "logLik".

nobs.sfareg returns a single integer (or NA).

See Also

psfm, sfm, zsfm, ttsfm, print.sfareg, summary.sfareg

sfm

Description

Implementation of the cross-sectional stochastic frontier model across an array of distributional assumptions for both v and u (user specified). For panel models, see the psfm() call.

Usage

sfm(formula, model_name = c("NHN", "NHN_Z", "NE", "NE_Z", "NR",
    "THT", "NTN", "NG", "NNAK", "NU", "NGE", "NLN", "NW", "tHN",
    "TSL"), data, maxit.bobyqa = 10000, maxit.nlminb = 500,
    maxit.psoptim = 1000, maxit.optim = 1000, REPORT = 1, trace = 2,
    pgtol = 0, start_val = FALSE, PSopt = FALSE, use.nlminb = "auto",
    use.bobyqa = "auto", optHessian = TRUE, inefdec = TRUE, upper =
    NA, Method = "L-BFGS-B", robust = c("mle", "mlqe", "psi",
    "mdpd"), c_mlqe = 0.2, eta = 0.01, alpha = 0.2, verbose = FALSE,
    start_from = NULL, Nsim = "auto", z_link = c("sd", "var"), vhet =
    NULL, uhet = NULL, muhet = NULL, scaling = NULL, shapehet = NULL,
    weights = NULL, wscale = TRUE, sim_type = c("halton", "sobol",
    "torus", "uniform"), antithetics = FALSE, sim_burn = NULL,
    sim_scrambling = 0L, sim_prime = NULL, sim_seed = NULL,
    rand.psoptim = NULL, keep_objective = FALSE, estimator = c("mle",
    "cols", "mols"), cols_boot = 0, rand.cols = NULL)

Arguments

Details

The options include the Normal-Half Normal (NHN), Normal-exponential (NE), Student's t-Half t (THT), the Normal-Truncated Normal (NTN), the Normal-Uniform (NU), the Normal-Generalized Exponential (NGE), the Normal-Lognormal (NLN), the Normal-Weibull (NW), the Student's t-half normal (tHN) and the Normal-Truncated Skew-Laplace (TSL). NHN_Z and NE_Z are extensions for the NHN and NE models that allow for modeling the u-component of those models with determinants of inefficiency.

Modified ordinary least squares (estimator = "cols" or "mols"). The moment estimator of Olson, Schmidt and Waldman (1980, Journal of Econometrics), standardly called MOLS. Note the name: this is not Winsten's corrected OLS, which shifts the intercept by the largest residual and estimates no variance parameters at all. The package has shipped this estimator under the spelling "cols" since it was introduced, so that spelling continues to work and is not deprecated; "mols" is accepted as an exact synonym so the method can be found under its usual name. OLS is consistent for the slopes of a composed-error frontier whatever the one-sided distribution; only the intercept is biased, by E[u], because the composed error does not have zero mean. COLS therefore keeps the OLS slopes, inverts the central moments of the OLS residuals for the scale parameters, and shifts the intercept up by the implied E[u]. Writing m_2, m_3 for the second and third central moments of e = v - u: for "NHN", m_3 = \sigma_u^3\sqrt{2/\pi}(1-4/\pi) and m_2 = \sigma_v^2 + \sigma_u^2(1-2/\pi) with E[u] = \sigma_u\sqrt{2/\pi}; for "NE", m_3 = -2\sigma_u^3 and m_2 = \sigma_v^2+\sigma_u^2 with E[u] = \sigma_u; for "NG", from the third and fourth cumulants, \sigma_u = -k_4/(3k_3) and \mu = -k_3/(2\sigma_u^3) with E[u] = \mu\sigma_u.

No optimizer runs and the answer is deterministic, which makes COLS a natural robustness check against a maximum-likelihood fit that may have settled at a local optimum. It is less efficient than ML when the distributional assumption holds, and the "NG" inversion in particular leans on a fourth cumulant and is noticeably noisier than the other two at moderate sample sizes.

Wrong skew under COLS. A production frontier implies m_3 < 0. When a sample comes out with m_3 \ge 0 the moment equations have no admissible solution: sfm() warns, reports sigu as 0 with the whole residual variance assigned to sigv, and returns no efficiency predictions. Read that as no evidence of inefficiency in these data, not as an estimate of zero. This is the Type I failure of Olson, Schmidt and Waldman and is common in small samples; $wrong_skew and $residual_moments record it.

NR is not NHN. "NR" places a Rayleigh distribution on u, not a half-normal, and the two are separate families rather than reparameterizations of one another: the Rayleigh density carries an extra factor of u, and no rescaling removes it. The practical consequence is that the standardized skewness each contributes to the composed error is a different constant, -0.631 against the half-normal's -0.995, so a sample generated under one cannot be fitted by the other. \sigma_u is reported on the second-raw-moment convention E[u^2] = \sigma_u^2, matching "NHN" and making the two directly comparable; the Rayleigh scale is then \sigma_u/\sqrt{2}, with E[u] = \sigma_u\sqrt{\pi}/2 and Var(u) = (1-\pi/4)\sigma_u^2. Starting values come from inverting these moments, since the Rayleigh skewness is a constant and so m_3 identifies Var(u) outright. Use data_gen_cs()'s y_pcs_r column to test it, not y_pcs.

tHN is not THT. In THT (Tancredi 2002) a single scale mixture is shared by both error components, so the noise and the inefficiency are both Student's t with the same degrees of freedom, and that shared mixture is what makes the composed error a closed-form skew-t. In tHN the noise is Student's t and the inefficiency is half-normal, drawn independently, so the two components have different tail behaviour, there is no closed form, and the density is the convolution f(e) = \int_0^\infty f_v(e+u) f_u(u) du evaluated by Gauss-Legendre quadrature. tHN is the heavy-tailed-noise model, which is what makes it the natural parametric comparison for the density-power robust estimators (robust = "mlqe", "psi", "mdpd"); THT cannot play that role because its inefficiency term is heavy-tailed too. Parameters are reported as (sigv, sigu, nu), the conventional order, not THT's inverted (sigu, sigv, a).

Two properties of tHN are worth knowing before using it. First, the degrees of freedom are weakly identified: on data simulated from the model at n = 1000 with a true \nu = 5, the profile log-likelihood varies by only about 0.24 across \nu from 10 to 100 and peaks near 20, so a single reported \nu should not be taken at face value. Profile over a grid of fixed \nu and report the path. Because of that flat ridge sfm() always fits tHN from several widely separated starting values, keeps the best, and records the outcome in thn_starts; it warns when they reach different optima. Second, on real data \sigma_u can collapse onto zero, the heavy noise tail absorbing the entire one-sided component and leaving mean predicted efficiency near one. That is a property of the model rather than a numerical failure, so it is reported rather than suppressed: sfm() warns and sets thn_sigma_u_at_bound. Fits are also substantially slower than the closed-form models, since every likelihood evaluation is a quadrature whose node count grows with \sigma_u/\sigma_v.

"NU" assumes u \sim U(0,\theta) (Li 1996; Nguyen 2010), giving the closed-form composed density f(e) = \theta^{-1}[\Phi((e+\theta)/\sigma_v) - \Phi(e/\sigma_v)]. The reported theta is the upper bound of the inefficiency support, so predicted inefficiency is bounded by construction. "NGE" assumes a generalized exponential F(u) = (1-e^{-\lambda u})^2, whose composed density is a difference of two exponentially-tilted normal tail terms; it nests the exponential model and permits a non-monotonic inefficiency density. NGE reports sigu on the same scale as NE (the exponential mean, \lambda = 1/\sigma_u) so the two are directly comparable.

The variance-determinant scale. sfm's _Z models place z'\delta on the standard deviation, \sigma_u = \exp(z'\delta), while psfm's TRE_Z/GTRE_Z place it on the variance, \sigma_u = \sqrt{\exp(z'\delta)}, as do the competing packages. Since \exp(\eta) = \exp(\eta/2)^2, fitting the same data under the two conventions returns the same \sigma_u, the same log-likelihood and the same marginal effects, but \delta under the SD link is exactly half \delta under the variance link. Reading a \delta from one family as though it came from the other therefore doubles or halves every reported effect. z_link exists so the two can be aligned; the default preserves sfm's own convention so existing results do not move.

Heteroskedasticity in more than one component. vhet, uhet and muhet let \sigma_v, \sigma_u and (for "NTN") the pre-truncation mean \mu each depend on their own covariates, in any combination. They are named formulas rather than further pipe segments deliberately: pipe position already carries different meanings in different model families – the second segment is \sigma_w in ttsfm but \sigma_h in psfm's GTRE_Z – and a fourth position would not be readable. The existing | z segment keeps its meaning, so nothing that worked before changes.

Two reasons to reach for vhet in particular. Ignored heteroskedasticity in v does not merely cost efficiency: Caudill, Ford and Gropper (1995) show that it biases the estimated frontier, because the composed error's mean depends on the scales. And a \sigma_v that varies with observed characteristics is one of the few specifications that can absorb apparent wrong skewness without forcing \sigma_u onto the zero boundary. muhet with model_name = "NTN" is the Battese and Coelli (1995) inefficiency-effects model, in which covariates shift the location of inefficiency rather than its dispersion; that is usually the model wanted when the question is “what makes firms inefficient”, since \delta_\mu then has a sign that reads directly.

These fits carry every positive quantity on the log scale, so the parameter vector is unconstrained and no variance can be driven onto a boundary by the optimizer; that is also why this path uses nlminb followed by an unbounded BFGS rather than the bounded four-stage scaffold the other models use. Coefficients are reported prefixed by block – Zv., Zu., Zmu. – since the same covariate may appear in more than one. marginal_effects reads the \sigma_u block of these fits with no extra arguments. estimator = "cols" and the robust divergences are moment-based and homoskedastic respectively, so neither combines with these arguments; both refuse rather than silently ignoring them.

Simulated-ML draws. "NLN" and "NW" have no closed-form composed density and are fitted by maximum simulated likelihood, averaging over draws built by sim_type. Two properties are enforced whatever the sequence. Each observation receives its own contiguous block of the sequence rather than a strided subsequence: a column-major fill hands observation i the elements h_i, h_{i+n}, \ldots of a van der Corput sequence, which is not equidistributed – at n = 500 and Nsim = 100 the first observation's block spanned only [0.50, 0.75], integrating over a quarter of the inefficiency distribution. And the draws are clamped away from 0 and 1, since qnorm() of an exact endpoint is infinite.

Blocking by observation is also what makes the draws worth using: Train (2002, ch. 9) attributes the advantage of a low-discrepancy sequence to its coverage and to the negative correlation it induces across observations, and the second only arises when observations receive different blocks. On draw counts, Bhat (2001) found 100 Halton draws more precise than 1000 pseudorandom ones, with 125 Halton draws giving half the simulation error of 1000 random draws; Nsim = "auto" scales the count with n and warns if a supplied count is too small.

"TSL" assumes a truncated skew-Laplace inefficiency (Wang 2012, eq. 3.7; the distribution is due to Aryal and Rao 2005, and the model was published as Nguyen, Gupta and Wang 2014), with density f(u) = \frac{1+\lambda}{\sigma_u(2\lambda+1)}\left[2e^{-u/\sigma_u} - e^{-(1+\lambda)u/\sigma_u}\right] for u \ge 0. This is a signed mixture of two exponentials – the second enters negatively – so it nests the exponential model as \lambda \to 0 while allowing a non-monotonic inefficiency density, and the composed density is the corresponding signed combination of two normal-exponential terms. It is evaluated as a difference in log space rather than by subtracting the raw terms: both carry \sigma_v^2/(2\sigma_u^2), which overflows once \sigma_u becomes small relative to \sigma_v, a region the optimizer routinely visits. lambda is a shape parameter and is the least sharply identified of the three. This is a property of the family, not of the fit, and it is worth understanding before reporting lambda at all. The shape parameter changes the standardised distribution remarkably little: E[U]/\sigma_u equals 1 in both limits \lambda \to 0 and \lambda \to \infty, and its maximum over all \lambda is only 4-2\sqrt{2} \approx \bold{1.1716}, attained at \lambda = 1/\sqrt{2}; the standardised skewness moves only from about 1.95 down to 1.79 and back. So the data identify the product \sigma_u \times m(\lambda) far better than they identify \sigma_u and \lambda separately, and the two trade off against each other.

Measured on the package's own convergence design (\sigma_u = 1, \lambda = 1.5, 60 replications at each of n = 500, 1000, 2000, 4000), estimating one parameter with the other held at its true value:

Each parameter attains the root-n rate when the other is known and fails it when the other is free, which is the signature of weak scale–shape separation rather than of an optimiser problem. \lambda's dispersion falls fortyfold once \sigma_u is fixed.

What to do about it. Report the functionals rather than the primitives where you can: E[U] converges at the root-n rate even with both parameters free (slope -1.03, sd 0.031 at n = 4000), and E[e^{-U}] nearly so (-0.89). Both are what an efficiency analysis actually uses. If \lambda itself is of interest, prefer profile-based inference over its reported t-value, which is a Wald statistic on a parameter whose sampling distribution is neither symmetric nor tightly concentrated.

"NLN" (u lognormal, reported as sigu = sdlog and mu = meanlog) and "NW" (u Weibull, reported as sigu = scale and k = shape) have no closed-form composed density and are estimated by simulated maximum likelihood. The integral f(e) = \int_0^\infty f_u(u)\,\phi((e+u)/\sigma_v)/\sigma_v\,du is a product of two densities, and simulating from either one alone fails wherever the other is the narrow factor: draws from f_u miss the normal kernel for a very inefficient unit, while draws from the kernel miss a sharply peaked f_u. Half the draws are therefore taken from each, and combined by the balance heuristic of Veach and Guibas (1995), which weights every draw by the proposal likelier to have produced it. This needs no rule for choosing between the two and remains accurate in both regimes. The draws are generated once and reused at every optimizer iteration, so the objective is deterministic and reproducible; increase Nsim to reduce simulation error.