## ----Init----------------------------------------------------------------
require(SFAt)
## ------------------------------------------------------------------------
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))
## ------------------------------------------------------------------------
# require(GGally)
# ggpairs(dtDT)
plot(data_ex1)
hist(data_ex1$y)
hist(data_ex1$u)
hist(data_ex1$eps)
## ------------------------------------------------------------------------
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)
## ------------------------------------------------------------------------
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))
## ------------------------------------------------------------------------
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))
## ------------------------------------------------------------------------
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))
## ------------------------------------------------------------------------
ex2_ineff_term <- predict(model_ex2,
type = "inefficiency",
estimator = "JLMS")
hist(ex2_ineff_term)
## ------------------------------------------------------------------------
ex2_efficiencies <- predict(model_ex2,
type = "efficiency",
estimator = "JLMS")
hist(ex2_efficiencies)
## ----formula method------------------------------------------------------
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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.