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

require(SFAt)

Example dataset

Lets generate some randomly simulated production data...

EX1 <- as.data.frame(sim_data_cs(500))

and explore the production data and simulated inefficiency term.

# require(GGally)
# ggpairs(dtDT)

plot(EX1)

hist(EX1$y)
hist(EX1$eps)
hist(EX1$u)

OLS model

We might be interested in simple OLS model for LR test (and starting values of parameters).

lmfit <- lm(y ~ x1 + x2 + z1 + z2, data = EX1)
summary(lmfit)

SFA model

sfat <- SFAt::SFA(y ~ x1 + x2 ,
                  cm = ~ z1 + z2,
                  data = EX1,
                  dist = "tnorm",
                  grad = "analytic",
                  deb = F,
                  optim_method = "BFGS",
                  optim_control = list(maxit = 3e4, trace = T)
)
summary(sfat)

benchmark with frontier package methods

require(frontier)

frontest <- frontier::sfa(y ~ x1 + x2 | z1 + z2,
                          data = EX1,
                          ineffDecrease = T,
                          maxit = 100)

print(summary(frontest))
print(lrtest(frontest))

fit_sfat <- predict(sfat, type = "efficiency", estimator = "JLMS")
fit_frontier <- frontier::efficiencies(frontest, asInData = FALSE)
fit_true <- exp(-EX1$u)
plot(x = fit_sfat, y = fit_true)
plot(x = fit_frontier, y = fit_true)
plot(x = fit_sfat, y = fit_frontier)
points(x = fit_sfat, y = fit_true, col = "red")
mean(fit_true - fit_sfat)
mean(fit_true - fit_frontier)

Example 2

data("front41Data")

sfa2 <- SFAt::sfa.fit(y = log(front41Data$output), 
                      X = cbind(capital = log(front41Data$capital), 
                                labour = log(front41Data$labour)), 
                      dist = "tnorm",
                      optim_method = "CG",
                      optim_control = list(trace = T),
                      deb = F)

print(summary(sfa2))

sfa2_front <- frontier::sfa( log( output ) ~ log( capital ) + log( labour ),
                             truncNorm = T,
                             data = front41Data )

print(summary(sfa2_front))
print(lrtest(sfa2_front))

vh-d/SFAt documentation built on May 3, 2019, 6:11 p.m.