Benchmarking (version 0.27)

sfa: Stochastic frontier estimation

Description

Estimate a stochastic frontier production function using a maximum likelihood method.

Usage

sfa(x, y, beta0 = NULL, lambda0 = 1, resfun = ebeta, 
    TRANSPOSE = FALSE, DEBUG=FALSE,
    control=list(), hessian=2)


te.sfa(object)
teBC.sfa(object)
teMode.sfa(object)
teJ.sfa(object)
te.add.sfa(object)
sigma2u.sfa(object)
sigma2v.sfa(object)
sigma2.sfa(object)
lambda.sfa(object)

Arguments

x

Input as a K x m matrix of observations on m inputs from K firms; (firm x input); MUST be a matrix. No constant for the intercept should be included in x as it is added by default.

y

Output; K times 1 matrix (one output)

beta0

Optional initial parameter values

lambda0

Optional initial ratio of variances

resfun

Function to calculate the residuals, default is a linear model with an intercept. Must be called as resfun(x,y,parm) where parm=c(beta,lambda) or parm=c(beta), and return the residuals as an array of length corresponding to the length of output y.

TRANSPOSE

If TRUE, data is transposed, i.e. input is now m x K matrix

DEBUG

Set to TRUE to get various debugging information written on the console

control

List of control parameters to ucminf

hessian

How the Hessian is delivered, see the ucminf documentation

object

Object of class ‘sfa’ as output from the function sfa

Value

The values returned from sfa is the same as for ucminf, i.e. a list with components plus some especially relevant for sfa:

par

The best set of parameters found c(beta,lambda).

value

The value of minus log-likelihood function corresponding to 'par'.

beta

The parameters for the function

sigma2

The estimate of the total variance

lambda

The estimate of lambda

N

The number of observations

df

The degrees of freedom for the model

residuals

The residuals as a K times 1 matrix/vector, can also be obtained by residuals(sfa-object)

fitted.values

Fitted values

vcov

The variance-covarians matrix for all estimated parameters incl. lambda

convergence

An integer code. '0' indicates successful convergence. Some of the error codes taken from ucminf are

'1' Stopped by small gradient (grtol).

'2' Stopped by small step (xtol).

'3' Stopped by function evaluation limit (maxeval).

'4' Stopped by zero step from line search

More codes are found in ucminf

message

A character string giving any additional information returned by the optimizer, or 'NULL'.

o

The object returned by ucminf, for further information on this see ucminf

Details

The optimization is done by the R method ucminf from the package with the same name. The efficiency terms are assumed to be half--normal distributed.

Changing the maximum steplength, the trust rgion, might be important, and this can be done by the option 'control= list(stepmax=0.1)'. The default value is 0.1 and that value is suitable for parameters around 1; for smaller parameters a lower value should be used. Notice that the step length is updated by the optimizing program and thus must be set for every call of the function sfa if it is to be set.

The generic functions print.sfa, summary.sfa, fitted.sfa, residuals.sfa, logLik.sfa, and coef.sfa all work as expected.

The methods te.sfa, teMode.sfa etc. calculates the efficiency corresponding to different methods

References

Bogetoft and Otto; Benchmarking with DEA, SFA, and R, Springer 2011; chapters 7 and 8.

See Also

See the method ucminf for the possible optimization methods and further options to use in the option control.

The method sfa in the package frontier gives another way to estimate stochastic production functions.

Examples

# NOT RUN {
# Example from the book by Coelli et al.
# d <- read.csv("c:/0work/rpack/front41Data.csv", header = TRUE, sep = ",")
# x <- cbind(log(d$capital), log(d$labour))
# y <- matrix(log(d$output))

n <- 50
x1 <- 1:50 + rnorm(50,0,10)
x2 <- 100 + rnorm(50,0,10)
x <- cbind(x1,x2)
y <- 0.5 + 1.5*x1 + 2*x2 + rnorm(n,0,1) - pmax(0,rnorm(n,0,1))
sfa(x,y)
summary(sfa(x,y))


# Estimate efficiency for each unit
o <- sfa(x,y)
eff(o)

te <- te.sfa(o)
teM <- teMode.sfa(o)
teJ <- teJ.sfa(o)
cbind(eff(o),te,Mode=eff(o, type="Mode"),teM,teJ)[1:10,]


sigma2.sfa(o)       # Estimated varians
lambda.sfa(o)       # Estimated lambda
# }