### First example with simulated data
#Set seed, sample size and type of function
set.seed(1337)
N=500 #Sample size
s=-1 #Set to production function
#Generate covariates
x1<-runif(N,-1,1); x2<-runif(N,-1,1); x3<-runif(N,-1,1)
x4<-runif(N,-1,1); x5<-runif(N,-1,1)
#Set parameters of the distribution
mu=2+0.75*x1+0.4*x2+0.6*x2^2+6*log(x3+2)^(1/4) #production function parameter
sigma_v=exp(-1.5+0.75*x4) #noise parameter
sigma_u=exp(-1+sin(2*pi*x5)) #inefficiency parameter
#Simulate responses and create dataset
y<-rcomper(n=N, mu=mu, sigma_v=sigma_v, sigma_u=sigma_u, s=s, distr="normhnorm")
dat<-data.frame(y, x1, x2, x3, x4, x5)
#Write formulae for parameters
mu_formula<-y~x1+x2+I(x2^2)+s(x3, bs="ps")
sigma_v_formula<-~1+x4
sigma_u_formula<-~1+s(x5, bs="ps")
#Fit model
model<-dsfa(formula=list(mu_formula, sigma_v_formula, sigma_u_formula),
data=dat, family=comper(s=s, distr="normhnorm"), optimizer = c("efs"))
#Model summary
summary(model)
#Smooth effects
#Effect of x3 on the predictor of the production function
plot(model, select=1) #Estimated function
lines(x3[order(x3)], 6*log(x3[order(x3)]+2)^(1/4)-
mean(6*log(x3[order(x3)]+2)^(1/4)), col=2) #True effect
#Effect of x5 on the predictor of the inefficiency
plot(model, select=2) #Estimated function
lines(x5[order(x5)], -1+sin(2*pi*x5)[order(x5)]-
mean(-1+sin(2*pi*x5)),col=2) #True effect
### Second example with real data
# \donttest{
data("RiceFarms", package = "plm") #load data
RiceFarms[,c("goutput","size","seed", "totlabor", "urea")]<-
log(RiceFarms[,c("goutput","size","seed", "totlabor", "urea")]) #log outputs and inputs
RiceFarms$id<-factor(RiceFarms$id) #id as factor
#Set to production function
s=-1
#Write formulae for parameters
mu_formula<-goutput ~ s(size, bs="ps") + s(seed, bs="ps") + #non-linear effects
s(totlabor, bs="ps") + s(urea, bs="ps") + #non-linear effects
varieties + #factor
s(id, bs="re") #random effect
sigma_v_formula<-~1
sigma_u_formula<-~bimas
#Fit model with normhnorm dstribution
model<-dsfa(formula=list(mu_formula, sigma_v_formula, sigma_u_formula),
data=RiceFarms, family=comper(s=-1, distr="normhnorm"), optimizer = "efs")
#Summary of model
summary(model)
#Plot smooths
plot(model)
# }
### Third example with real data of cost function
# \donttest{
data("electricity", package = "sfaR") #load data
#Log inputs and outputs as in Greene 1990 eq. 46
electricity$lcof<-log(electricity$cost/electricity$fprice)
electricity$lo<-log(electricity$output)
electricity$llf<-log(electricity$lprice/electricity$fprice)
electricity$lcf<-log(electricity$cprice/electricity$fprice)
#Set to cost function
s=1
#Write formulae for parameters
mu_formula<-lcof ~ s(lo, bs="ps") + s(llf, bs="ps") + s(lcf, bs="ps") #non-linear effects
sigma_v_formula<-~1
sigma_u_formula<-~s(lo, bs="ps") + s(lshare, bs="ps") + s(cshare, bs="ps")
#Fit model with normhnorm dstribution
model<-dsfa(formula=list(mu_formula, sigma_v_formula, sigma_u_formula),
data=electricity, family=comper(s=s, distr="normhnorm"),
optimizer = "efs")
#Summary of model
summary(model)
#Plot smooths
plot(model)
# }
Run the code above in your browser using DataLab