# Clear the decks and define the roundup and gmean functions ----
  
  rm(list=ls())
  roundup = function(x,n){sign(x)*trunc(abs(x)*10^n + 0.5)/10^n}
  gmean = function(x){exp(mean(log(x)))}
  
# Set the working directory to the source file location ----
  
  library(rstudioapi)
  setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
  
  setwd ("XXXXXX")
  
  # Install and load packages ----
  #install.packages('Rtools')
  #install.packages('readxl')
  #install.packages('plm')
  #install.packages('frontier')
  #install.packages('TeachingDemos')

  library(readxl)
  library(frontier)
  library(plm)
  library(TeachingDemos)
  library(haven)
  
#txtStart("XXXX.txt") # Uncomment this if want to send output to the log file "XXXX.txt" ---- 

# Data in excel format ----
  
  data <- read_excel("XXX.xls")
  data <- na.omit(data)
  colnames(data)[1] = "Firm"            # rename the firm identifier
  colnames(data)[2] = "Year"            # rename the period identifier
  attach(data)
  summary(data)
  
  q = real_production_value          
  qname = "Output"
  
  x1 = Labour                          # First input should be labour (to calculate LPI later)
  x2 = real_capital
  x3 = real_intermed_input
  xdata = cbind(x1,x2,x3)            
  xnames = c("Labour", "Land", "Machinery", "EMS")
  
  a1 = m1_age
  a2 = m2_gender
  a3 = m3_immigrant
  a4 = m4_edu
  a5 = m4_edu_D1
  a6 = m4_edu_D2
  a7 = m5_kids
  a8 = m6_inc
  #adata = cbind(a1,a2,a3,a4,a5,a6,a7,a8)
  adata = cbind(a1,a2,a3,a5,a6,a7,a8)
  #anames = c("Age", "Gender", "Immigrant", "Edu", "Edu_D1", "Edu_D2", "Kids", "Income")
  anames = c("Age", "Gender", "Immigrant", "Edu_D1", "Edu_D2", "Kids", "Income")
 
  z1 = m10_norrland
  z2 = m11_storstad                
  zdata = cbind(z1,z2)
  znames = c("Norrland", "Storstad")

# Descriptive Statistics ----
  
  N = 1                        # no. of outputs
  M = ncol(xdata); M           # no. of inputs
  J = ncol(zdata); J           # no. of environmental variables
  K = ncol(adata); K           # no. of attributes
  I = length(unique(Firm)); I  # no. of firms
  T = length(unique(Year)); T  # no. of time periods
  nobs = nrow(data); nobs      # no. of observations
  I*T - nobs                   # = 0 if a balanced panel
  
  df = data.frame(Year,Firm,cbind(q,xdata,zdata,adata))
  colnames(df) = c("Year","Firm",qname,xnames,znames,anames)
  summary(df)

# SFA ----
  
  Period = Year - min(Year) + 1
  
  y = log(q)
  t = Period
  lxdata = log(xdata)

  mypanel <- pdata.frame(data, c( "Firm", "Year" ) )
  sfm = sfa(y ~ t + zdata + lxdata | adata, data = mypanel)
  OTE = efficiencies(sfm, asInData = TRUE)
  out = summary(sfm)
  out
  
  bml = coefficients(out)[,1:3]
  alpha = bml[1,1]
  lambda = bml[2,1]
  k1 = 3;        k2 = k1 + J - 1;       delta = bml[k1:k2,1]
  k1 = k2 + 1;   k2 = k1 + M - 1;       beta = bml[k1:k2,1]
  k1 = k2 + 1;   k2 = k1 + K;           phi = bml[k1:k2,1]
  roundup(c(alpha, lambda, delta, beta, phi),3)
  
  eta = sum(beta)
  eta 
  e = y - (alpha + lambda*t + zdata%*%delta + lxdata%*%beta)

# Write coefficient estimates to "XXXX.csv" ----

  write.csv(as.data.frame(coefficients(out)), file = "XXXX_coefficients.csv", row.names = TRUE)

# LPI and TFPI ----

  # LPI
  LP = q/x1
  LPI = LP/LP[1]

  # Multiplicative TFPI
  nbeta = beta/eta
  X = exp(lxdata%*%nbeta)
  TFP=q/X
  TFPI=TFP/TFP[1]
  
  results = data.frame(Firm,Year,roundup(cbind(LPI,TFPI),3))
  colnames(results) = c("Firm","Year","LPI","TFPI")
  head(results)
  tail(results)
  cor(LPI,TFPI)

# TFPI = OTI x OEI x OSMEI x OTEI x SNI ----

  OTI = exp(lambda*(t-1))
  zz = zdata%*%delta
  OEI = exp(zz - zz[1])
  xx = exp(lxdata%*%beta) / q
  OSMEI = TFPI*xx/xx[1]
  OTEI = OTE/OTE[1]
  SNI = TFPI/(OTI*OEI*OSMEI*OTEI)

  # check
  u = -log(OTE)
  v = e + u
  max(abs(SNI - exp(v-v[1])))
  
# OTEI = AI[1] x ... x AI[K] x A_noise ----  
  
  A = sweep(adata, 2, phi[2:(K+1)], `*`)
  AI = exp(sweep(A, 2, A[1,], '-'))
  A_noise = OTEI/apply(AI, 1, prod)
  results = data.frame(Firm,Year,AI,A_noise)
  colnames(results) = c("Firm", "Year", anames,"A_noise")
  head(results)
  k3 = ncol(results)
  max(abs(OTEI  - apply(results[,3:k3], 1, prod)))

# Round index numbers 3 decimal places and write to "XXXX_indexes.csv" ----
  
  results = data.frame(Firm,Year,roundup(cbind(LPI,TFPI,OTI,OEI,OSMEI,OTEI,SNI,AI,A_noise,OTE),3))
  colnames(results) = c("Firm","Year","LPI","TFPI","OTI","OEI","OSMEI","OTEI","SNI",anames,"A_noise","OTE")
  head(results)
  write.csv(results, file = "sni14_indexes.csv", row.names = FALSE)
 
txtStop() # Uncomment this if output has been going to the log file ----