In vh-d/SFAt: SFAt: Stochastic frontier analysis on cross-section and panel data

require(SFAt)

Basic SFA models

Example dataset

Lets generate some randomly simulated production data,

data_ex1 <- as.data.frame(sim_data_cs(1000, 
                                      ineff = -1, # production function
                                      x_coeff = c(50, 4, 2), 
                                      z_coeff = c(20), 
                                      sigma_u = 10, sigma_v = 3))

and explore the production data and simulated inefficiency term.

# require(GGally)
# ggpairs(dtDT)

plot(data_ex1)
hist(data_ex1$y)
hist(data_ex1$u)
hist(data_ex1$eps)

SFA -- normal/t-normal model

model_ex1 <- sfa.fit(y = data_ex1$y,
                     X = cbind(a = data_ex1$x1, b = data_ex1$x2),
                     CM = NULL,   
                     dist = "tnorm",
                     optim_control = list(temp = 1, tmax = 20, maxit = 6e4)
                     # deb = T, debll = T
)

print(summary(model_ex1))
print(lrtest(model_ex1))

ineff_fit <- inefficiencyTerm(model_ex1, estimator = "JLMS")
hist(ineff_fit)

Check for separation of frontier and actual production data

yhat <- predict(model_ex1, type = "frontier")
hist(yhat,       xlim = c(0, max(data_ex1$y)))
hist(data_ex1$y, xlim = c(0, max(data_ex1$y)))

hist(inefficiencyTerm(model_ex1, estimator = "JLMS"))

hist(efficiency(model_ex1))

SFA model with endogenous inefficiency mean

Generate random cross-section dataset

data_ex2_cm <- as.data.frame(sim_data_cs(1000, 
                                         x_coeff = c(50, 4, 2), 
                                         z_coeff = c(20, -2, 2), 
                                         sigma_u = 10, sigma_v = 5))

Fit BC95 model

model_ex2 <- sfa.fit(y = data_ex2_cm$y,
                     X = cbind(a = data_ex2_cm$x1, b = data_ex2_cm$x2),
                     CM = cbind(c = data_ex2_cm$z1, d = data_ex2_cm$z2),
                     dist = "tnorm",
                     # deb = T, debll = T,
                     optim_method = "BFGS",
                     optim_control = list(maxit = 50e3))

# print(model_ex2$coeff_frontier)
# print(model_ex2$coeff_cm)
# print(model_ex2$coeff_cv_u)
# print(model_ex2$coeff_cv_v)
print(summary(model_ex2))
print(lrtest(model_ex2))

Fitted inefficiency term ($E(u_i|\epsilon_i)$)

ex2_ineff_term <- predict(model_ex2, 
                          type = "inefficiency", 
                          estimator = "JLMS")

hist(ex2_ineff_term)

Fitted efficiency (%)

ex2_efficiencies <- predict(model_ex2, 
                            type = "efficiency", 
                            estimator = "JLMS")

hist(ex2_efficiencies)

sfa1 <- SFA(y ~ x1 + x2,
            data = data_ex1,
            dist = "tnorm",
            form = "production",
            # cm = ~ 0,  
            # grad = "fd",
            grad = "analytic",
            opt_strategy = 1,
            nlopt_opts = list("algorithm" = "NLOPT_LN_PRAXIS", print_level = 1, maxeval = 1e4, xtol_rel = 1e-10,
                              "local_opts" = list("algorithm" = "NLOPT_LN_PRAXIS", print_level = 1, maxeval = 1e3, xtol_rel = 1e-10)),
            nlopt_bounds = list(lb = -5, ub = 60),
            optim_method = "BFGS",
            optim_control = list(trace = 1, maxeval = 10000),
            maxLik_method = "NR",
            maxLik_control = list(printLevel = 3)
            #,deb = T#, debll = T
)

print(summary(sfa1))
print(lrtest(sfa1))