Stochastic Frontier Analysis
Overview
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.
A feature release. In brief:
- **A Wald test for whether the endogeneity correction was needed**,
endogeneity_test(), from Hou, Ramalho and Roseta-Palma (2025). On anivsfm()fit it testsH0: rho = 0-- that the noise is uncorrelated with the reduced-form errors, so a plainsfm()fit would have been consistent. The null is INTERIOR, since rho lives in the open unit ball, so this is an ordinary chi-square and not the chi-bar-square mixtureinefficiency_test()needs. Size 0.046 at a nominal 5% over 1000 replications; power 1.000 at rho = 0.4 and 0.8. - **
ivsfm()now reports standard errors forrho.** They wereNA, on the ground thatrho = t/sqrt(1 + t't)has a non-diagonal Jacobian and a delta-method value computed as though it were diagonal would be wrong. The first half was right and the conclusion was not: the Jacobian isJ = (I - rho rho')/s, soVar(rho) = J Var(t) J'exactly. Checked against a numerical Jacobian, and validated end-to-end by the size of the test above -- a wrong Jacobian would show up there as a mis-sized test. The full covariance travels with the fit as$vcov_rho, which the joint test needs.
Show all changes in 1.2.0Hide the rest of the changes in 1.2.0
- **
ttsfm(model_name = "TTHN")no longer fabricates a log-density it cannot compute.** The likelihood needsD = Phi2(x1, 0; rho1) - Phi2(x2, 0; rho2), a difference of two bivariate normal CDFs of the same order of magnitude, and in parts of the parameter space it cancels completely: atsigma_v = 0.3, sigma_u = 1, sigma_w = 0.2on a 400-observation draw one observation'sDcomes back as exactly 0 while both CDFs are O(0.1).pmax(D, .Machine$double.xmin)did not repair that -- it inventedlog(D) = -708.4for the observation, and that single observation moved the summed objective by **715 log-units** in a surface whose real curvature is a few units per 0.1 step in log sigma. A cliff that size is what L-BFGS-B reports asABNORMAL_TERMINATION_IN_LNSRCH.Dis now rejected when it falls below the rounding error of its own subtraction, which leaves genuine small probabilities alone: aDof 6e-94 where both CDFs are also ~1e-94 has full significance and is kept. Three smaller repairs in the same branch: the exponentials forsigma_v/sigma_u/sigma_ware clipped atEXP_CLIP_UPPERas theTTNEbranch already clipped its own;rhois held strictly inside(-1, 1), whichsigma_vsmall enough does reach; and the substitution of-sqrt(.Machine$double.xmax / n)-- about-9.5e152at n = 200, 140 orders of magnitude past the finite penalty it sat beside -- is gone, replaced by that penalty. It had also mapped a+Inflog-density to a large NEGATIVE one, so an unusable draw could be rewarded. Measured on eight matched fits: one that previously failed outright now converges, and one that previously returnedsigma_w = 0.044against a truth of 1 now returns 1.047 **with a better likelihood**. The other six are unchanged.sigma_vstill collapses toward zero on some draws; that is a separate problem and is not fixed. - **
ttsfm()no longer throws away a fit because of the optimizer's exit code.** It was the only entry point callingstop()on ANY non-zero convergence code, including L-BFGS-B's code 52,ABNORMAL_TERMINATION_IN_LNSRCH, which it returns whenever its line search meets a discontinuity and usually *having improved* the objective. The guard was written when nothing checkedopt$convergenceat all;opt.optim()now rebuilds at the stage-2 point whenever the value or the Hessian is non-finite, so what was left of it only destroyed usable fits. Two of eight TTHN test fits returned code 52 with sound estimates where thestop()had aborted them. - **What skewness and kurtosis a pair of distributions can actually produce**,
moment_range(), from Papadopoulos and Parmeter (2021).spec_test()already shipped their formal statistic; this is the range check that comes before it. Because the variance shareRlies in(0,1), the composed error can never be more skewed thanuitself -- Normal-Half-Normal is confined togamma1 in (-0.995, 0), Normal-Exponential to(-2, 0)-- so residual moments outside that set refute the pair with no test statistic involved. That is their Type II failure, alongside the Type I failure of wrong skew. Five noise and five inefficiency distributions, including the two whose shape is a free parameter: the truncated normal sweepsgamma1(u) in (0, 2)andgamma2(u) in (-0.243, 6)and is the onlyuhere whose composed error can be platykurtic under normal noise, while the gamma is unbounded in both and so can **never** suffer a Type II failure. The kurtosis range is not the interval between the two components: with leptokurtic noise it dips below both, to0.674for Laplace noise with a half-normaluagainst the half-normal's own0.869. **It is a diagnostic, not a test, and the print method says so.** PP2021's Tables 1-5 replicate -- their population rows exactly, and Tables 2, 3 and 5 cell for cell -- and they show the range check over-rejects: with correctly specified data atn = 200and SNR 3, the sample excess kurtosis leaves its theoretical range 26% of the time. Usespec_test()for the p-value. - **A frontier that tolerates the wrong skewness**,
esfm(), from Hafner, Manner and Simar (2018). When the sample skewness of the residuals comes out positive -- a small-sample accident, not evidence against the model -- the classical MLE collapses:sigma_uis exactly zero, every firm is fully efficient and the fit says nothing. This model keeps the frontier and widens the distribution ofu, with one parametergammawhose SIZE is the scale of inefficiency and whose SIGN chooses the direction of skewness. The classical model is nested atgamma > 0. On a 50-observation wrongly skewed draw,sfm(model_name = "NHN")returnssigma_u = 0.0003and mean efficiency 0.99976;esfm()returnsgamma = -0.57with the slopes still near truth. Becausegamma = 0is an **interior** point rather than a boundary, the LR test of "no inefficiency" is an ordinarychi2(1)-- not the chi-bar-square mixtureinefficiency_test()needs for the classical model. Size over 2000 replications on the paper's design, at a nominal 5%: 0.178 / 0.082 / 0.054 / 0.061 / 0.058 for n = 50 / 100 / 200 / 500 / 1000, against the paper's 0.152 / 0.076 / 0.061 / 0.054 / 0.053. Badly oversized below n = 100, usable from n = 200 -- their own conclusion. **It is not a free improvement everywhere.** The paper reports smaller bias in mean efficiency than the classical model even when the population skewness has the correct sign, and that reproduces when inefficiency is not small: atgamma = 0.5the bias is -0.0001 at n = 100 against the classical +0.021. But atgamma = 0.3, where 22% of samples are wrongly skewed, this model overshoots the other way (-0.040 against +0.036) and is the worse of the two at n = 50 (-0.080 against +0.044). A small truegammalets it fit sizeable inefficiency of either sign to noise. Use it when the wrong skewness is the problem. - **Kim and Schmidt (2008) two-step test**,
uhet_test(). The procedure applied work actually runs -- fit the frontier, regress the JLMS predictor on firm characteristics, test the slopes -- is invalid as usually performed, becauseu-hatis a *generated* dependent variable and the first-step estimation error enters the second-step variance. The correction is not distribution-free, so"NHN"and"NE"are implemented, the two cases the paper works through. What it buys: withzindependent ofxthe corrected and naive tests agree, as the theory requires. Withzstrongly correlated withx(rho = 0.9) the naive standard error is **52% too large** at n = 800 and the naive test's size collapses to **0.000** at a nominal 5%, while the corrected standard error matches the true sampling spread to within half a percent (ratio 0.995) and keeps size 0.032. **Caveat, stated plainly:** the corrected test is conservative in these simulations. On the paper's own base case it rejects 0.024/0.036/0.028 at N = 200/500/1000 against a nominal 0.05, where the paper reports 0.0495 at N = 200. The *naive* branch does reproduce the paper's figure for the uncorrected test (0.023-0.029 against their 0.0237). The implementation was checked against equations (5)-(13), the first-step scores agree with an analytic half-normal score to 1e-11,gamma-hatis unbiased and Gaussian, and an oracle test using the true sampling standard deviation has correct size (0.045) -- so the gap is in the variance estimate, not ingamma-hator the reference distribution. Two candidate explanations were tested and **both refuted**: the homoskedastic second-step variance (the paper's "OLS" against its "OLS-H") changes nothing, and neither does using the Hessian information rather than the outer product when formingr_i. The gap is unexplained rather than merely unverified. Read a rejection as meaningful and a non-rejection as weak evidence. Fits on the wrong-skew boundary are refused rather than answered: everyE[u|eps]is then identical, so there is no estimated inefficiency forzto explain, and the information matrix is singular there (Waldman 1982). - **Chen and Wang (2012) centered-residuals moment test**,
cw_test(). Unlikegof_test(), it never evaluates the composed-error density or distribution function: it compares the empirical characteristic function of the *centred* residuals with its theoretical value, which follows from the characteristic functions of the two components separately -- in closed form for many pairs whose convolution is not, the normal-gamma model being the standard example. Centring is what makes it work. A frontier's intercept is not identified separately fromE[u], so the least-squares intercept estimatesalpha - E[u]; centring cancels it, which is why the test is valid off an ordinary regression rather than requiring maximum likelihood. The variance corrects for both nuisances -- estimating the scales, and the centring itself. Omitting the second understates it and oversizes the test. Size at a nominal 5% over 500 replications, half-normal data: | test | n | tau=1 | tau=1.5 | tau=2 | |---|---|---|---|---| | cosine | 500 | 0.070 | 0.054 | 0.050 | | cosine | 2000 | 0.046 | 0.044 | 0.050 | | sine | 500 | 0.092 | 0.064 | 0.042 | | sine | 2000 | 0.068 | 0.062 | 0.064 | Power against exponential inefficiency attau = 1is 0.231 (n = 500) and 0.952 (n = 2000) for the cosine test against 0.108 and 0.706 for the sine test. Both of the authors' conclusions reproduce: the cosine test is the less sensitive totauand the more powerful, so it is the default. **taudefaults to a single frequency**: passing several at once gives rejection rates of 0.09-0.15 at a nominal 5%, because the cosine moments at nearby frequencies are near-collinear and the covariance matrix is then near-singular. The authors report the same for their combinedtau. A wrongly signed third moment leaves the moment equations with no admissible solution, and that is an error rather than a number. - **The Gauss-Legendre rule is now cached.** It depends on nothing but the node count and was rebuilt by an eigendecomposition on every call -- hundreds of times per
cw_test(), sincenumDerivdifferences a quadrature-based moment function. Caching also keepseigen()out of forked workers: macOS's threaded BLAS is not fork-safe andmclapply()segfaulted inside it before this was added. - **Goodness-of-fit tests for the assumed inefficiency distribution**,
gof_test(), from Wang, Amsler and Schmidt (2011). The distribution ofuis the assumption in this model least often defended and least often tested, and it is testable: hold the normality ofvfixed and the assumeduimplies a distribution for the composed error, so rejecting that is rejecting the assumedu. Kolmogorov-Smirnov and Pearson chi-square, both on the composed error rather than onu-hat-- which is the same test, sinceu-hat = E[u|eps]is monotonic ineps, but far easier, and it avoids the trap the paper warns about: comparing the spread ofu-hatwith the assumed density ofuis a mistake, not a diagnostic, because those are different distributions. Parameter estimation is not ignored. The default parametric bootstrap copies the estimation step exactly -- each replication draws from the fitted model, rebuilds the response, **refits**, forms its own residuals and recomputes the statistic. The"asymptotic"option gives a chi-square p-value that the authors note is conservative at the MLE, and **returnsNArather than a KS p-value**, because with parameters estimated the Kolmogorov distribution does not apply and Bai's (2003) transformation is not implemented. Size and power over 200 replications (half-normal data for size, exponential data fitted as half-normal for power), at nominal 0.10 / 0.05: | | n | KS | chi2 | |---|---|---|---| | size | 200 | 0.065 / 0.025 | 0.085 / 0.065 | | size | 800 | 0.070 / 0.040 | 0.065 / 0.010 | | power | 200 | 0.505 / 0.395 | 0.235 / 0.145 | | power | 800 | 0.985 / 0.965 | 0.810 / 0.655 | | power | 2000 | 1.000 / 1.000 | 0.990 / 0.985 | Both hold their size, mildly conservatively, and **KS dominates chi-square at every sample size** -- the paper's own conclusion, and the reason it is the better default. - **Coelli (1995)'s hypothesis tests, and his COLS standard errors.** New
inefficiency_test()tests H0: no technical inefficiency, reporting the one-sided LR test, the naive LR test, the Wald ratio and the third-moment test side by side. The point of the paper is that **the two tests most often reported are the two with the wrong size**: H0 puts gamma on the boundary, so the LR statistic follows the Gourieroux-Holly-Monfort mixture0.5 chi2(0) + 0.5 chi2(1)(5% critical value 2.71, not 3.84) and the Wald ratio is not asymptotically normal at all. Measured here over 400 replications per cell at a nominal 5%, with sigma_u = 0: | n | LR one-sided | LR naive | Wald | M3T | |---|---|---|---|---| | 100 | 0.060 | 0.030 | 0.230 | 0.055 | | 400 | 0.055 | 0.033 | 0.217 | 0.045 | | 800 | 0.060 | 0.035 | 0.182 | 0.062 | The one-sided LR and third-moment tests 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 the one-sided LR test has the better power of the two valid tests at every sample size, which is Coelli's recommendation.NTNfits get the two-restriction mixture0.25/0.5/0.25overchi2(0,1,2). - **
sfm(estimator = "cols")now reports Coelli's analytic standard errors.** It previously returnedNAforsigvandsiguunless a bootstrap was requested, because the OLS errors are not appropriate for them: those parameters are non-linear functions of the residual momentsm2andm3, not regression coefficients. Coelli's Appendix 1 (A15-A16) gives the correct delta-method variance and it is now filled in forNHN; the intercept additionally carries the error in itsE[u]shift.NEandNGkeepNAand the bootstrap, because the moment inversion -- and therefore the variance -- is distribution-specific. Checked against the nonparametric bootstrap that was already there, on the same data: the analytic/bootstrap ratio forsiguis 1.047 at n = 200, 1.011 at n = 800 and 0.984 at n = 3200, which is an asymptotic approximation behaving as one should. Coelli's A13 inversion turned out to be the moment inversion the package already carried as.gtre_two_step()(equal to ten decimal places, not merely similar), and his A15 variance is the delta method whose GTRE_SEQ2 implementation is corrected elsewhere in this release -- **including the factor of 2 on the covariance term, which A15 independently confirms**, and A16's `mu6 - mu3^2 - 6 mu2 mu4 + 9 mu2^3
The NEWS section is longer than shown here.
From NEWS.md in the release
Install
Health
- OK13 OK2026-08-25
- WARNING12 OK, 1 WARNING2026-08-20
- OK13 OK2026-08-04
- NOTE12 OK, 1 NOTE2026-08-01
- OK13 OK2026-06-09
Show 3 earlier changesHide 3 earlier changes
- ERROR12 OK, 1 ERROR2026-06-08
- WARNING12 OK, 1 WARNING2026-06-07
- OK14 OK2026-03-10
Documentation
- Help pages
- 47includes a package overview page
- Examples that run
- 43%20 run in a plain check · 25 only when donttest examples are enabled (R CMD check --as-cran) · 1 never run
- Exported functions with examples
- 100%
- Documented parameters
- 84%
- Return-value docs
- 100%
- Help pages citing references
- 68%
- Introduction to psfm HTML
Downloads
Repository
- PRs merged (all time)
- 96%
- Median days to close, 50 most recently updated closed issues
- 0
- Median days to close, 50 most recently updated closed PRs
- 0
Stars over time
Forks over time
Issues over time
PRs over time
Commits over time
Commit identities over time
3 development-tooling and community-health practices detected across 2 families in the upstream repository
Checks run against github.com/davidhbernstein/sfa on 2026-10-18.
The R CMD check workflow runs on Linux, macOS and Windows, and on R-devel
This repository shows signs of AI coding tools.
- Claude Codefirst seen 2026-08-20Pull request it wrote: opened by a project memberCommits: at least 211 crediting it in the message
- Claude Opus 5101 commits · Aug to Sep 2026
- Claude Opus 5.56 commits · Sep 2026
Some of these counts are a minimum, because the search stopped before the oldest commits.
Dependencies
Installing sfa adds 27 packages (16 with compiled code, 14 MB of source tarballs). Heaviest: Rcpp 1.8 MB, collapse 1.4 MB, randtoolbox 1.4 MB.
Code & Tests
Datasets
People & History
Previously maintained by Ariane Straub (1.0-0 to 1.0-1).
Pick two releases to compare their code metrics.
- 1.2.0Latest2026-09-07 · current release · diff ↗ ·
added 33 exports
lcsfm, selsfm, ivsfm, copsfm, efficiency_ci, skewness_test, marginal_effects, efficiency, meanefficiency, simulation_se, TIC, vuong, lcsfm_homogeneity, spec_test, spec_test_all, sfma, influence_sfa, hscore, hscore_select, calibrate_c, density_weights, pcomposed, dcomposed, inefficiency_test, pcomposed_model, composed_cdf, gof_test, cw_test, uhet_test, esfm and 3 more - 1.1.5
- RR 4.6.0 released · 2026-04-24
- unarchivedReturned to CRAN2026-01-21
- 1.0.42026-01-21 · diff ↗ ·
removed 5 exports
sfa, dgp, eff, u.sfa, te.eff.sfaadded 7 exports
zsfm, psfm, sfm, data_gen_p, data_gen_cs, summary.sfareg, print.sfareg - RR 4.5.0 released · 2025-04-11
- RR 4.4.0 released · 2024-04-24
- RR 4.3.0 released · 2023-04-21
- archivedRemoved from CRAN2022-06-15check problems were not corrected despite reminders
- RR 4.2.0 released · 2022-04-22
- RR 4.1.0 released · 2021-05-18
- RR 4.0.0 released · 2020-04-24
- RR 3.6.0 released · 2019-04-26
- RR 3.5.0 released · 2018-04-23
- RR 3.4.0 released · 2017-04-21
- RR 3.3.0 released · 2016-05-03
Show 9 earlier eventsHide 9 earlier events
- RR 3.2.0 released · 2015-04-16
- RR 3.1.0 released · 2014-04-10
- 1.0-12014-01-06 · diff ↗
- RR 3.0.0 released · 2013-04-03
- RR 2.15.0 released · 2012-03-30
- RR 2.14.0 released · 2011-10-31
- RR 2.13.0 released · 2011-04-13
- 1.0-02011-01-24
- RR 2.12.0 released · 2010-10-15
Package metadata
- First published
- 2011-01-24
- Total releases
- 5 / 15 yrs
- License
- GPL (>= 2) OSI
- Minimum R
- ≥ 4.0.0
- Bundled data
- 391 KB / 4 files
- Download size
- 1.0 MB
Cite
Cite this package
Run in R for the authors' entry:
citation("sfa")The package provides a citation file, so citation() returns the authors' own entry.
Cite the R Observatory
For a number measured here: a download total, a coverage figure, an archival date.
From data release v2026-10-22, which the citation names so these numbers can be found later. More on citing and the projects behind them.