Package {sfa}

Version:	1.1.5
Date:	2026-08-23
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, Benchmarking
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-08-23 12:05:32 UTC; davidbernstein
Author:	David Bernstein [aut, cre], Christopher Parmeter [aut], Alexander Stead [aut]
Repository:	CRAN
Date/Publication:	2026-08-23 18:30: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
USUtilities             USUtilities
data_gen_cs             Generate Cross-Sectional Data for Stochastic
                        Frontier Analysis
data_gen_p              Generate Panel Data for Stochastic Frontier
                        Analysis
npsfm                   Nonparametric Stochastic Frontier Models
panel89                 Panel89
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
sfa-package             Stochastic Frontier Analysis
sfa_diagnostics         Optimizer Diagnostics for Stochastic Frontier
                        Fits
sfareg-methods          Extractor methods for "sfareg" objects
sfm                     sfm
summary.sfareg          Summary method for "sfareg" objects
ttsfm                   Two-Tier Stochastic Frontier Model
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

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)

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)

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)

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)

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"), 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)

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}. A fit that falls into it can have a lower likelihood than the true parameter vector, so it is an optimizer failure rather than a limit of the data. 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.

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:

For SSFE 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 model is likewise not maximum-likelihood based and returns 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)     

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)

max_tre_z   <-  psfm(formula    = y_tre_z ~ x1 +x2| z_gtre,
                     model_name = "TRE_Z",
                     data       = data_trial,
                     individual = "name",
                     PSopt      = TRUE)

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", "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)

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, ...)
## S3 method for class 'sfareg'
logLik(object, ...)
## S3 method for class 'sfareg'
nobs(object, ...)

Arguments

Details

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"),
    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.20, eta = 0.01, alpha = 0.2, verbose = FALSE,
    Nsim = "auto", rand.psoptim = NULL, keep_objective = FALSE,
    estimator = c("mle", "cols"), 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) and the Student's t-half normal (tHN). 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.

Corrected ordinary least squares (estimator = "cols"). The moment estimator of Olson, Schmidt and Waldman (1980, Journal of Econometrics). 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.

"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: f(e) = E_u[\phi((e+u)/\sigma_v)/\sigma_v] is averaged over Nsim Halton draws per observation. The draws are generated once and reused at every optimizer iteration, so the objective is deterministic and reproducible; increase Nsim to reduce simulation error.

Outputs include E[exp(-u)|e] given by exp_u_hat, following Battese and Coelli (1988, JoE), where appropriate. As of this version exp_u_hat is also returned for NE, NTN, NU, NGE, NLN and NW, together with the Jondrow et al. (1982) point predictor u_hat = E[u|e]; previously the exponential and truncated-normal models returned no efficiency prediction at all.

Note on default iteration counts: sfm()'s optimization defaults (thousands of iterations) are much higher than psfm()'s (100/10/10). This is intentional: sfm()'s likelihoods are closed-form (no simulation), so each objective-function evaluation is cheap and many default iterations cost little; several of psfm()'s models require simulated maximum likelihood over Halton draws, where each evaluation is far more expensive. See psfm's Details for more.

Robust divergence-based estimation (robust argument). For model_name = "NHN", sfm() can maximize a robust divergence objective instead of the ordinary log-likelihood, downweighting the influence of outlying observations:

• "mlqe": minimum Lq estimation (Lq, Ferrari and Yang 2010), power parameter c_mlqe (q = 1-c). Objective \sum (f(e_i)^c - 1)/c.

• "psi": Psi-likelihood (Basu et al.-type divergence family), power parameter eta. Objective \sum f(e_i)^c/c - n \, I_{1+c}/(1+c), where I_{1+c} = \int f(e;\theta)^{1+c}\,de is evaluated numerically at the current parameter values on every optimizer iteration.

• "mdpd": minimum density power divergence (Basu et al. 1998), power parameter alpha. Objective ((1+c)/c) \sum f(e_i)^c - n \, I_{1+c}.

All three reduce to "mle" as their power parameter approaches 0. Psi and MDPD produce identical parameter point estimates for the same power parameter (c = \eta = \alpha) – the two objectives are a positive rescaling of one another, \text{MDPD}(c) = (1+c)\cdot \text{Psi}(c), which does not change the argmax; both are exposed for compatibility with the corresponding paper's three-method terminology, not because they differ numerically. MLqE is a genuinely distinct estimator.

Because these are M-estimators rather than ordinary MLE, the naive Hessian-inverse standard error (optHessian = TRUE's usual output) is not statistically valid for them – it assumes the information-matrix equality, which these objectives do not satisfy in general. When robust != "mle", sfm() instead reports a sandwich-form standard error, A^{-1} B A^{-1}, where A is the Hessian of the summed robust objective (as already computed for the point estimate) and B = \sum_i g_i g_i' is built from the per-observation gradient of that objective (obtained via numerical differentiation, jacobian). This requires optHessian = TRUE; std.errors/t.values are NA otherwise, same as for ordinary MLE.

Formulas verified against an independently audited reference implementation (cross-checked by direct numerical integration and finite differences), including the sandwich SE's own verification (the c -> 0 limit against the ordinary MLE Hessian SE, and a Monte Carlo check against the empirical spread of repeated fits). Currently only implemented for model_name = "NHN" – other models will error clearly if robust != "mle" is requested.

Value

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

Note

Standard errors require optHessian set to TRUE. When robust != "mle" (model_name = "NHN" only), a sandwich-form SE is reported instead of the ordinary Hessian-inverse one – see Details.

Author(s)

David H. Bernstein and Alexander Stead

See Also

psfm for panel models, zsfm for zero-inefficiency models, ttsfm for two-tier models, data_gen_cs for simulating data with known true parameters, and sfareg-methods for extracting results.

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)

## Robust divergence-based estimation (NHN only): downweights the influence
## of outlying observations relative to ordinary MLE. std.errors/t.values
## are NA for these methods -- see Details.
cs.nhn.mlqe <-  sfm(formula = y_pcs ~ x1 + x2, model_name = "NHN",
                    data = cs_data_trial, robust = "mlqe", c_mlqe = 0.20)
cs.nhn.psi  <-  sfm(formula = y_pcs ~ x1 + x2, model_name = "NHN",
                    data = cs_data_trial, robust = "psi", eta = 0.20)

Summary method for "sfareg" objects

Description

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

Usage

## S3 method for class 'sfareg'
summary(object, ...)

Arguments

Details

Prints the call, the model name, and the estimated parameters with standard errors and t-values, together with the log-likelihood and the number of observations where the estimator provides them. Which parameters appear depends on the model: several report the \lambda = \sigma_u/\sigma_v, \sigma = \sqrt{\sigma_u^2 + \sigma_v^2} reparameterization rather than the raw scale parameters.

Value

Invisibly returns its object argument unchanged, having printed the summary as a side effect.

Author(s)

David H. Bernstein

See Also

print.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)
summary(cs.nhnz)                    

Two-Tier Stochastic Frontier Model

Description

Code to estimate the two-tier stochastic frontier model (2TSF), which decomposes a symmetric error term into two one-sided components (e.g. buyer/seller informational rents) in addition to the usual two-sided noise term.

Usage

ttsfm(formula, model_name = c("TTNE", "TTHN", "TTNLS"),
data, maxit.bobyqa = 80000, maxit.psoptim = 1000, maxit.optim = 1000,
REPORT = 1, trace = 0, pgtol = 0, start_val = FALSE, PSopt = FALSE,
optHessian = TRUE, inefdec = TRUE, upper = NA,
Method = "L-BFGS-B", logit = TRUE, verbose = FALSE, rand.psoptim = NULL)

Arguments

Details

The two-tier stochastic frontier model (Polachek and Yoon, 1987) models the composed error as v + u - w, where v is a symmetric normal noise term and u, w are independent one-sided components. "TTNE" assumes u and w are exponential; "TTHN" assumes they are half-normal; "TTNLS" makes no distributional assumption and instead estimates \sigma_u/\sigma_w by nonlinear least squares via the scaling property. Determinants of u and w can be included via the second and third parts of a two-pipe formula, which parameterize \sigma_u and \sigma_w through a log link.

Like the other composed-error models in this package, all three two-tier likelihoods can have a genuine (non-error) boundary mode where a variance component is estimated near its numerical floor in some samples; spot-check fits before serious use. "TTHN"'s post-estimation information-deficiency metrics (the metrics component of the returned object) involve a large number of bivariate-normal-CDF evaluations and are more experimental than the parameter estimates themselves – NULL if the underlying calculation fails. For "TTNLS", the reported standard error for the \sigma_v-position parameter (present only for layout consistency with "TTNE"/"TTHN", unused in the NLS objective) is always NA by construction.

Note on default iteration counts: like sfm()/zsfm(), ttsfm()'s likelihoods are closed-form, so its optimization defaults are set much higher than psfm()'s (which uses simulated maximum likelihood for several models and defaults to far fewer iterations for speed). See psfm's Details for more.

Value

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

Note

Standard errors require optHessian set to TRUE

Author(s)

David H. Bernstein and Christopher F. Parmeter

References

Polachek, S. and Yoon, B.J. (1987). A Two-Tiered Earnings Frontier Estimation of Employer and Employee Information in the Labor Market. The Review of Economics and Statistics.

Examples

library(sfa)

data <- data_gen_cs(N = 200, rand = 1, sig_u = 0.3, sig_v = 0.2,
                     cons = 2, beta1 = 0.5, beta2 = 0.3, a = 5, mu = 0)

eq <- y_pcs ~ x1 + x2

ttsfm(formula    = eq,
      model_name = "TTNE",
      data       = data)

Zero-Inefficiency Stochastic Frontier Models

Description

Fits the zero-inefficiency stochastic frontier model of Kumbhakar, Parmeter and Tsionas (2013), in which the sample is a mixture of fully efficient firms and firms drawn from an ordinary stochastic frontier.

Usage

zsfm(formula, model_name = c("ZISF", "ZISF_Z"), 
data, maxit.bobyqa = 10000,maxit.psoptim = 1000, maxit.optim = 1000, 
REPORT = 1, trace = 0, pgtol = 0,start_val = FALSE,PSopt = FALSE, 
optHessian = TRUE, inefdec = TRUE, upper = NA, 
Method = "L-BFGS-B",logit = TRUE,verbose=FALSE,rand.psoptim = NULL)

Arguments

Details

The zero-inefficiency stochastic frontier model treats the sample as a mixture of two regimes: with probability p a firm sits exactly on the frontier (u = 0, noise only), and with probability 1 - p it is drawn from an ordinary normal-half normal frontier. It is the natural specification when some fraction of the sample is thought to be fully efficient, which an ordinary frontier cannot represent because it forces u > 0 almost surely.

Under "ZISF" the mixing probability is a single constant, parameterized as p = \exp(-|\gamma|). That parameterization is exactly symmetric in \gamma: +\gamma and -\gamma give identical fits and the optimizer may return either, so the sign of the reported gamma carries no meaning. Under "ZISF_Z" the probability varies with covariates z through the link chosen by logit.

post.prob reports the posterior probability that each observation belongs to the fully efficient regime, given its residual; jlms reports the Jondrow et al. (1982) inefficiency prediction.

Note on default iteration counts: like sfm(), zsfm()'s likelihood is closed-form, so its optimization defaults (thousands of iterations) are set much higher than psfm()'s (which uses simulated maximum likelihood for several models and defaults to far fewer iterations for speed). See psfm's Details for more.

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

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.

See Also

sfm, psfm, ttsfm, sfareg-methods, panel89

Examples

library(sfa)  

eqz     <- y ~ q1 + q2 + q3 + q4 + q5 + w1 + w2 + w3 + w4 | z

data(panel89)

zsfm(formula    = eqz,
     model_name = "ZISF_Z",
     data       = panel89,
     logit      = TRUE)