require(SFAt)
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)
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))
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)
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.