###################################################
#An example of using the functions in file model.R 
#to estimate stochastic production frontier models.
#Runs on R version 2.5.0
#
#Ricardo Ehlers, Jun 2006
###################################################

#Source the R functions to run JAGS and load required R packages.

source('model.R')
start()

#Reading data file and creating variables
#This data is also available from
#http://econ.queensu.ca/jae/1998-v13.2/zellner-ryu/data.zr

data = matrix(scan(file='data.zr',strip.white=F,skip=10),ncol=4)

#log output, log labour and log capital
y = log(data[,1]/data[,4])
X = cbind(log(data[,2]/data[,4]), log(data[,3]/data[,4]))

## JAGS reads files: name-inits.R, name-data.R, name.bug, name.cmd
name='zellner'

mcmc = list(niter=50000, nburn=50000, nthin=5)

#Prior median efficiency (r*) and hyperparameters d for Gamma, Weibull
#and Generalized Gamma inefficiencies. Parameters a and b for 
#tau~Gamma(a,b) and prior precision tau.beta for the coefficients.

prior = list(rstar=0.8, d=c(3,3),tau=c(0.01,0.01),tau.beta=0.001,ln.lambda=c(0.01,0.01),tau.v=0.001, tau.rho=0.001,delta=c(1,1),gpf=c(0.1,0.01))

#Initial values for the MCMC updates
inits=list(beta=0, tau=1, u=rep(1,length(y)), gpf.gamma=0.5, ln.lambda=1)

eff=c('exp','tnorm','halfnorm','gamma','gen.gamma','weibull','lognorm')

ff=c('Cobb-Douglas','Translog','CES','GPF')

m=spf(name,y,X,fun=ff[1],eff=eff[1],mcmc,prior,inits,quiet=F,center=T)

#The following files should have been created: 
#zellner.bug, zellner.cmd, zellner-data.R, zellner-inits.R
#jags.out, jags.ind

#Histograms for the best firm and worst firm ranks

par(mfrow=c(1,2))
hist(m$ranks[,1],prob=T,nclass=25,xlab='rank',ylab='',main='',col='grey')
hist(m$ranks[,2],prob=T,nclass=25,xlab='rank',ylab='',main='',col='grey')

##Diagnostics and graphical output
##Use option select according to your model choices

diag.mcmc(m,select=1:6)

plot.mcmc(m$sample,type='trace',smooth=T,select=1:6,ask=F)
plot.mcmc(m$sample,type='autocorr',select=1:6,ask=F)
plot.mcmc(m$sample,type='hist',select=1:6,ask=F)


#### End of file ####