######################################################################
#Functions to estimate Stochastic Frontier models running JAGS from R.
#Runs on R version 2.5.0
#
#JAGS stores arrays columnwise, like R, not rowwise like BUGS.
#Create data and initial values in R and write to the files using
#the dump function.
#
#Compare combinations of inefficiency term distributions: 
#exponencial, half-normal, truncated normal, gamma, log-normal,
#weibull and generalized gamma and production functions: 
#Cobb-Douglas, Translog, Constant Elasticiy of Substitution (CES) and 
#Generalized Production Function (GPF).
#
#Author: Ricardo Ehlers Jun 2006. Last change Apr 2007.
#####################################################################

start = function() {
 cat('Loading required R packages\n')
 pkgs = c('coda','stats')
 for (i in 1:length(pkgs)) require(pkgs[i],character.only=T,quietly=F)
 options(width=85)
}

plot.mcmc=function(x,saveplot=F,smooth=F,file=NULL,type='trace',ask=T,mfcol=NULL,add=F,select=NULL,...){
#####################################################################
#Traces, acf, histograms and density estimates of MCMC output
#Set saveplot=T and provide file name for redirecting to a Postscript
#file.
#Use select to specify a subset of variables to monitor
#Author: Ricardo Ehlers Nov 2005
#####################################################################
 if (length(select)>0) x=x[,select]     
 nvar=nvar(x)
 names=colnames(x)
 if (!add) {
   if (saveplot) postscript(file=file,horizontal=F)
   if (length(mfcol)==0) mfcol=coda:::set.mfrow(Nchains = nchain(x),Nparms=nvar,nplots=1)
   oldpar = par(mfcol = mfcol)
   oldpar = c(oldpar, par(ask = ask))
 } else {
   dev.set(dev.cur())
 }
 if (type=='trace'){
   traceplot(x,smooth=smooth)
 } else if (type=='autocorr'){
   autocorr.plot(x)
 } else if (type=='hist') {
   for (j in 1:nvar) hist(x[,j],prob=T,xlab='',main=names[j])
 } else if (type=='dens') {
   for (j in 1:nvar) densplot(x[,j])
 }
 if (saveplot) dev.off()
}

diag.mcmc=function(x,select=NULL){
#######################################################
#MCMC Diagnostics.
#Use select to specify a subset of variables to monitor
#Author: Ricardo Ehlers Nov 2005
#######################################################
 if (length(select)>0) x=x[,select]
 cat('\nEffective sample size for estimating the mean\n')
 print(effectiveSize(x))
 cat('\nGeweke convergence diagnostic\n')
 print(geweke.diag(x))
 cat('\nAutocorrelations\n')
 print(autocorr.diag(x))
 cat('\nHeidelberger and Welch convergence diagnostic\n')
 print(heidel.diag(x))
 cat('\nRaftery and Lewis diagnostic\n')
 print(raftery.diag(x))
 if (nchain(x)>1){
   cat('\nGelman and Rubin convergence diagnostic\n')
   print(gelman.diag(x))
 }
}

jags.cmd = function(vars,data,name,mcmc,seed){
###############################################################
#Creates name.cmd file to be called from JAGS.
#Author: Ricardo Ehlers Jan 2006.
# vars: nodes to be monitored.
# data: name for data file
# name: name for files with JAGS commands, model specifications 
# and initial values.
# draw: number of MCMC draws.
# nburn: burn-in size.
# nthin: MCMC thinning.
###############################################################
  file=paste(name,'.cmd',sep='')
  text=NULL
  nvar=length(vars)
  if (!missing(seed)) text=paste(text,"seed ",seed,"\n",sep="")
  text=paste(text,'model in "',name,'.bug"\n',sep='')
  text=paste(text,'data in "',data,'-data.R"\n',sep='')
  text=paste(text,'compile\n',sep='')
  text=paste(text,'inits in "',name,'-inits.R"\ninitialize\n',sep='')
  text=paste(text,'update ',as.integer(mcmc$nburn),'\n',sep='')
  aux = NULL
  if (mcmc$nthin > 1) aux=paste(', ','thin(',mcmc$nthin,')',sep='')
  for(i in 1:nvar) text=paste(text,'monitor set ',vars[i],aux,'\n',sep='')
  text=paste(text,'update ',as.integer(mcmc$niter),"\n",sep='')
  text=paste(text,"coda *\nexit\n",sep="")
  cat(text,file=file)
}

jags.bug = function(name,N,p,y,X,fun,eff,prior,center=F){
#########################################################
#Creates name.bug file to be called from JAGS.
#Author: Ricardo Ehlers Jan 2006.
#########################################################
  file=paste(name,'.bug',sep='')
  a=NULL
  a=paste(a,'## ',fun,' production function ##\n',sep='')
  a=paste(a,'## inefficiency distribution ',eff,' ##\n',sep='')
  a=paste(a,'data {\n',sep='')
  a=paste(a,'  for (i in 1:N) {\n',sep='')
  if (fun=='GPF') a=paste(a,'    zeros[i] <- 0\n',sep='')
  if (center & fun!='CES' & fun!='Translog') {
    a=paste(a,'    X[i,1] <- Xreg[i,1]-xbar[1]\n',sep='')
    a=paste(a,'    X[i,2] <- Xreg[i,2]-xbar[2]\n',sep='')
    a=paste(a,'  }\n',sep='')
    a=paste(a,'  xbar[1] <- mean(Xreg[,1])\n',sep='')
    a=paste(a,'  xbar[2] <- mean(Xreg[,2])\n',sep='')
  } else { 
    a=paste(a,'    X[i,1] <- Xreg[i,1]\n',sep='')
    a=paste(a,'    X[i,2] <- Xreg[i,2]\n',sep='')
    a=paste(a,'  }\n',sep='')
  }
  if (fun=='Translog') {
    a=paste(a,'  for (i in 1:N) {\n',sep='')
    a=paste(a,'    X[i,3] <- X[i,1]*X[i,1]\n',sep='')
    a=paste(a,'    X[i,4] <- X[i,2]*X[i,2]\n',sep='')
    a=paste(a,'    X[i,5] <- X[i,1]*X[i,2]\n',sep='')
    a=paste(a,'  }\n',sep='')
  }
  a=paste(a,'}\n',sep='')
  a=paste(a,'model {\n',sep='')
  a=paste(a,'  for(i in 1:N) {\n',sep='')
  if (fun!='GPF'){
    a=paste(a,'    y[i] ~ dnorm(mu[i],tau)\n',sep='')
  } else {
    a=paste(a,'    zeros[i] ~ dpois(p[i])\n',sep='')
    a=paste(a,'    p[i] <- -0.5*log(tau/(2*3.141593))+0.5*tau*pow(y[i]+gamma*exp(y[i])-mu[i],2)-log(1+gamma*exp(y[i]))+10000\n',sep='')
  }
  if (center & fun!='CES' & fun!='Translog') {
    a=paste(a,'    mu[i] <- alpha0 +',sep='')
  } else {
    a=paste(a,'    mu[i] <- alpha +',sep='')
  }
  if (fun!='CES') {
    a=paste(a,'inprod(beta[],X[i,]) - u[i]\n',sep='')
  } else {
    a=paste(a,'(-v/rho)*log((1-delta)*exp(-rho*X[i,1])+delta*exp(-rho*X[i,2]))',sep='')
    a=paste(a,'-u[i]\n',sep='')
  }
  a=paste(a,'  }\n',sep='')
  if (fun!='CES') {
    a=paste(a,'  for(j in 1:',p,') {\n',sep='')
    a=paste(a,'    beta[j] ~ dnorm(0.0,',prior$tau.beta,')T(0,)\n',sep='')
    a=paste(a,'  }\n',sep='')
  } else {
    a=paste(a,'  v  ~ dnorm(0.0, 1E-3)\n',sep='')
    a=paste(a,'  delta ~ dbeta(1,1)\n',sep='')
    a=paste(a,'  rho  ~ dnorm(0.0, 1E-3)T(-1,)\n',sep='')
  }
  if (center & fun!='CES' & fun!='Translog') {
    a=paste(a,'  alpha0 ~ dnorm(0.0,',prior$tau.beta,')\n',sep='')
  } else {
    a=paste(a,'  alpha ~ dnorm(0.0,',prior$tau.beta,')\n',sep='')
  }
  if (fun=='GPF') a=paste(a,'  gamma ~ dgamma(',prior$gpf[1],',',prior$gpf[2],')\n',sep='')  
  a=paste(a,'  tau ~ dgamma(',prior$tau[1],',',prior$tau[2],')\n',sep='')
  a=paste(a,'  sigma2 <- 1/tau\n',sep='')
  a=paste(a,'  for (i in 1:N){\n',sep='')
  if (eff=='exp'     ) a=paste(a,'    u[i] ~ dexp(lambda)\n',sep='')
  if (eff=='halfnorm') a=paste(a,'    u[i] ~ dnorm(0,invlambda)T(0,)\n',sep='')
  if (eff=='tnorm'   ) a=paste(a,'    u[i] ~ dnorm(zeta,invlambda)T(0,)\n',sep='')
  if (eff=='gamma'   ) a=paste(a,'    u[i] ~ dgamma(phi,lambda)\n',sep='')
  if (eff=='weibull' ) a=paste(a,'    u[i] ~ dweib(c,lambda)\n',sep='')
  if (eff=='lognorm' ) a=paste(a,'    u[i] ~ dlnorm(0.0,lambda)\n',sep='')
# if (eff=='gen.gamma'| eff=='weibull') {
  if (eff=='gen.gamma'                ) {
    a=paste(a,'    uc[i] ~ dgamma(phi,lambda)\n',sep='')
    a=paste(a,'     u[i] <- pow(uc[i],invc)\n',sep='')
  }
  a=paste(a,'    eff[i] <- exp(-u[i])\n',sep='')
  a=paste(a,'  }\n',sep='')
  if (eff=='exp') a=paste(a,'  lambda ~ dexp(-log(',prior$rstar,'))\n',sep='')
  if (eff=='halfnorm') {
    a=paste(a,'  invlambda ~ dgamma(1,1/37.5)\n',sep='')
    a=paste(a,'  lambda <- 1/invlambda\n',sep='')
  }
  if (eff=='tnorm') {
    a=paste(a,'  invlambda ~ dgamma(1,1/37.5)\n',sep='')
    a=paste(a,'  zeta ~ dnorm(0.0, 1E-3)T(0,)\n',sep='')
    a=paste(a,'  lambda <- 1/invlambda\n',sep='')
  }
  if (eff=='gamma') {
    a=paste(a,'  phi <- 1/invphi\n',sep='')
    a=paste(a,'  invphi ~ dgamma(',prior$d[1],',',prior$d[1]+1,')\n',sep='')
    a=paste(a,'  lambda ~ dgamma(phi,-log(',prior$rstar,'))\n',sep='')
  }
  if (eff=='weibull') {
    a=paste(a,'  lambda ~ dexp(pow(-log(',prior$rstar,'), c))\n',sep='')
    a=paste(a,'    invc ~ dgamma(',prior$d[1],',',prior$d[1]+1,')\n',sep='')
    a=paste(a,'       c <- 1/invc\n',sep='')
  }
  if (eff=='lognorm') {
  a=paste(a,'  lambda~ dgamma(',prior$ln.lambda[1],',',prior$ln.lambda[2],')\n',sep='')
  a=paste(a,'  psi2 <- 1/lambda\n',sep='')
  }
  if (eff=='gen.gamma'                ) {
# if (eff=='gen.gamma'| eff=='weibull') {
    a=paste(a,'  lambda ~ dgamma(phi,pow(-log(',prior$rstar,'), 1/invc))\n',sep='')
    a=paste(a,'    invc ~ dgamma(',prior$d[1],',',prior$d[1]+1,')\n',sep='')
    a=paste(a,'       c <- 1/invc\n',sep='')
    if (eff=='gen.gamma'){
      a=paste(a,'  invpsi ~ dgamma(',prior$d[2],',',prior$d[2]+1,')\n',sep='')
      a=paste(a,'     phi <- invc/invpsi\n',sep='')
    } else {
      a=paste(a,'     phi <- 1\n',sep='')
    }
  }
  if (center & fun!='CES' & fun!='Translog') {
    a=paste(a,'# RE-calculating intercept on original scale\n',sep='')
    a=paste(a,'  alpha <- alpha0-inprod(beta[],xbar[])',sep='')
  }
  a=paste(a,'}\n',sep='')
  cat(a,file=file)
}

create.inits = function(name,fun,eff,N,p,inits,center=F) {
####################################################
#Creates name-inits.R file to be called from JAGS.
#Author: Ricardo Ehlers Jan 2006.
####################################################
 file=paste(name,'-inits.R',sep='')
 aux = NULL
 if (center & fun!='CES' & fun!='Translog') {
   alpha0 = inits$beta
   aux = c(aux,'alpha0')
 } else {
   alpha  = inits$beta
   aux = c(aux,'alpha')
 }
 if (fun!='CES') {
   beta  = rep(inits$beta,p)
   aux = c(aux,'beta')
 } else {
   v = 0
   delta = 0.5
   rho = 1
   aux=c(aux,'v','delta','rho')
 }
 if (fun=='GPF') {
   gamma = inits$gpf.gamma
   aux=c(aux,'gamma')
 }
 if (eff=='gen.gamma') {
   uc = inits$u
   aux=c(aux,'uc')
 } else {
   u = inits$u
   aux=c(aux,'u')   
 }
 tau = inits$tau
 aux=c(aux,'tau')
 dump(aux,file='temp')
 if (eff=='lognorm') {
   lambda=inits$ln.lambda
   dump('lambda',file='temp',append=T)
 }
 if (eff=='weibull') {
  lambda = 1
  invc = 1
  dump(c('lambda','invc'),file='temp',append=T)
 }
 if (eff=='halfnorm') {
   invlambda = 0.001
   dump('invlambda',file='temp',append=T)
 }
 expr = paste('sed s/\\`/\\"/g', ' temp >',sep='')
 system(paste(expr, file))
 system('rm temp')
}

create.data = function(name,y,X,N) {
#################################################
#Creates name-data.R file to be called from JAGS.
#Author: Ricardo Ehlers Jan 2006.
#################################################
 file=paste(name,'-data.R',sep='')
 Xreg= X
 dump(c("N","y","Xreg"),file='temp')
 expr = paste('sed s/\\`/\\"/g', ' temp > temp1', sep='')
 system(expr)
 expr = paste('sed s/L//g', ' temp1 > ', sep='')
 system(paste(expr,file))
 system('rm temp temp1')
}

create.rank = function(N,x){
#Ranking efficiency terms from JAGS output
#Author: Ricardo Ehlers Jan 2006.
 a=x[,(ncol(x)-N+1):ncol(x)]
 ranks = as.mcmc(t(apply(a,1,rank,ties='min')))
 b = summary(ranks)
 aux = apply(ranks,2,median)
 names(aux)=NULL
 aux1 = sort(aux,index=T,dec=T)
 b$statistics = b$statistics[aux1$ix,]
 b$quantiles  = b$quantiles [aux1$ix,]
 best = aux1$ix[1]
 worst= aux1$ix[N]
#return(b)
 invisible(list(summary=b, ranks=ranks[,c(best,worst)]))
}

dic = function(sample,y,X,p,fun,eff) {
####################################################
#DIC from JAGS output. Author: Ricardo Ehlers Jun 2006.
#Dbar = posterior mean of -2logL
#Dhat = -2logL at posterior mean of stochastic nodes
####################################################
 N=length(y)
 niter= nrow(sample)
 x = matrix(sample,niter,ncol(sample))
 ind = (ncol(x)-N+1):ncol(x)
 x[,ind] = log(x[,ind])
 mean = apply(x,2,mean)
 if (fun=='Translog') {
   X = cbind(X,X[,1]**2)
   X = cbind(X,X[,2]**2)
   X = cbind(X,X[,1]*X[,2])
 }
 D=matrix(NA,niter,1)

 if (fun=='Cobb-Douglas'| fun=='Translog') {
 for (i in 1:niter) {
   mu  = x[i,1]
   mu  = mu + X %*% (x[i,2:(1+p)]) + x[i,ind]
   D[i]= -2*sum(dnorm(y-mu,0,sqrt(x[i,p+2]),log=T))
 }}
 if (fun=='CES') {
   for (i in 1:niter) {
     mu = x[i,1]
     mu=mu-(x[i,2]/x[i,4])*log((1-x[i,3])*exp(-x[i,4]*X[,1])+x[i,3]*exp(-x[i,4]*X[,2]))
     mu=mu + x[i,ind]
     D[i]= -2*sum(dnorm(y-mu,0,sqrt(x[i,p+2]),log=T))
   }
 }

 if (fun=='GPF'){
   for (i in 1:niter) {
     mu  = x[i,1]
     mu  = mu + X %*% (x[i,2:(1+p)]) + x[i,ind]
     s1 = sum(dnorm(y+x[i,p+3]*exp(y)-mu,0,sqrt(x[i,p+2]),log=T))
     s2 = sum(log(1+x[i,p+3]*exp(y)))
     D[i]= -2*(s1+s2)
   }
 }

 Dbar = mean(D)

 if (fun!='CES') {
   mu = mean[1]
   mu = mu + X %*% mean[2:(1+p)] + mean[ind]
 } else {
   mu = mean[1]   
   mu=mu-(mean[2]/mean[4])*log((1-mean[3])*exp(-mean[4]*X[,1])+mean[3]*exp(-mean[4]*X[,2]))
   mu=mu + mean[ind]
 }
 if (fun!='GPF'){
   Dhat= -2*sum(dnorm(y,mu,sqrt(mean[p+2]),log=T))
 } else {
   s1 =   sum(dnorm(y+mean[p+3]*exp(y),mu,sqrt(mean[p+2]),log=T))
   s2 =   sum(log(1+mean[p+3]*exp(y)))
   Dhat= -2*(s1+s2)
 }
 pD = Dbar - Dhat
 DIC = Dbar + pD
 a = c(Dbar,Dhat,pD,DIC)
 names(a) = c('Dbar','Dhat','pD','DIC')
 write(c(fun,eff,round(a,5)),file='dic.out',append=T,ncolumns=6)
 return(a)
}

spf=function(name,y,X,fun='Cobb-Douglas',eff='exp',mcmc,prior,inits,quiet=T,center=F){
#######################################################################
#sample: sampled values of parameters and efficiencies
#ranks: ranks of the best and worst firms
#Author: Ricardo Ehlers Jun 2006.
#######################################################################
 #Creating required file name.cmd with JAGS commands
 N = length(y)
 p = ifelse(fun=='CES',3,2)
 if (fun=='Translog') p = 5
 if (fun=='Cobb-Douglas'|fun=='Translog') vars=c('beta','sigma2')
 if (fun=='GPF') vars=c('beta','sigma2','gamma')
 if (fun=='CES') vars=c('v','delta','rho','sigma2')
 vars=c('alpha',vars)
 if (eff=='lognorm')   vars=c(vars,'psi2')
 if (eff=='exp'|eff=='halfnorm')vars=c(vars,'lambda')
 if (eff=='tnorm')     vars=c(vars,'zeta','lambda')
 if (eff=='gamma')     vars=c(vars,'phi','lambda')
 if (eff=='gen.gamma') vars=c(vars,'c','phi','lambda')
 if (eff=='weibull')   vars=c(vars,'c','lambda')
 vars=c(vars,'eff')
 jags.cmd(vars=vars,name=name,data=name,mcmc=mcmc)
 jags.bug(name,N,p,y,X,fun=fun,eff=eff,prior,center)
 create.inits(name,fun=fun,eff=eff,N,p,inits,center)
 create.data(name=name,y,X,N)
 M = mcmc$nburn + mcmc$niter
 cat('\nMCMC for Stochastic Production Frontier Models',date(),'\n')
 cat(' Production Function: ',fun,'\n')
 cat(' Inefficiency term distribution:', eff,'\n')
 cat(' Iterations = ',mcmc$nburn+1,':',as.integer(M),sep='','\n')
 cat(' Thinning interval = ',mcmc$nthin,'\n')
 cat(' Sample size per chain = ',ceiling(mcmc$niter/mcmc$nthin),'\n\n')

 # Running JAGS
 start=Sys.time()
 if (quiet) {
   system(paste('jags < ',name,'.cmd > /dev/null',sep=''))
 } else {
   system(paste('jags ',name,'.cmd',sep=''))
 }
 print(Sys.time()-start)
 x=as.mcmc(read.coda('jags.out','jags.ind',quiet=T))
 ind = (ncol(x)-N+1):ncol(x)
 b=summary(x[,-ind])
 print(round(cbind(b$statistics,b$quantiles),4))
 cat('\n Computing ranks...\n')
 b=create.rank(N,x)
 cat('\n 5 best median rank firms\n')
 print(round(cbind(b$summary$statistics[1:5,],b$summary$quantiles[1:5,]),4))
 cat('\n')
 cat('\n 5 worst median rank firms\n')
 print(round(cbind(b$summary$statistics[(N-5):N,],b$summary$quantiles[(N-5):N,]),4))
 cat('\n')
 DIC = dic(x,y,X,p,fun,eff)
 print(DIC)
 invisible(list(sample=x,ranks=b$ranks,fun=fun,eff=eff,dic=DIC))
}

#### End of R functions ####