rm(list = ls())
#######################################################################################################
###ENTER LOCATIONS OF NECESSARY FILES
#######################################################################################################
FuncPath = ""
dataPath = ""
outPath = "/output.txt"

#######################################################################################################
###ENTER LOCATIONS OF NECESSARY FILES
#######################################################################################################

#Load in functions that carry out SUR estimation
source(FuncPath)
#read benchmarking data into memory
bmdata <-
  read.csv(datapath,
           stringsAsFactors = F)


#create price variables and make adjustments as needed. We divide by one of the prices
#in order to make cost function homogenous in input prices.
bmdata$wm <- bmdata$womhl
bmdata$wk <- bmdata$wkhl / bmdata$wm
bmdata$c <- bmdata$ch / bmdata$wm

#create share variables
bmdata$sk <- bmdata$skh
bmdata$sm <- bmdata$soh



#create lists of variables to be logged and variables to be demeaned
dvarlist <-
  c("wk",
    "ych",
    "yvh"
  )
lvarlist <- c(dvarlist, 
              "c"
)

#demean variables
dbmdataout <- meanscale(bmdata, dvarlist)
dbmdata <- dbmdataout[[1]]
dbmmeans <- dbmdataout[[2]]

#log variables - NB: Do not need to add 1 to 0 values. logvars function does that.
ldbmdataout <- logvars(dbmdata, lvarlist, T)
ldbmdata <- ldbmdataout[[1]]

#generate constant and trend variable
ldbmdata$constant <- 1
ldbmdata$trend <- bmdata$year - 1994


#create parameters for hetsur
peglhs <- c("c", "sk")
pegrhs1 <-
  c(
    "constant",
    "wk",
    "ych",
    "yvh",
    "wk*wk/2",
    "ych*ych/2",
    "yvh*yvh/2",
    "wk*ych",
    "wk*yvh",
    "ych*yvh",
    "trend"
  )


#second share equation, should correspond to second variables in peglhs
pegrhs2 <- c("constant",
             "wk",
             "ych",
             "yvh")



pegrhs <- c(pegrhs1, pegrhs2)
pegnumvar <- c(length(pegrhs1), length(pegrhs2))

pegfactors <- c("k")


#create the pegR and pegz matrices
pegrestrict <-
  surcostRmat(
    lhs_vars = peglhs,
    rhs_vars = pegrhs,
    num_var = pegnumvar,
    factors = pegfactors,
    constantvar = "constant",
    verbose = F
  )

pegR <- pegrestrict[[1]]
pegz <- pegrestrict[[2]]

pegRnames <- as.data.frame(pegR)
names(pegRnames) <- pegrhs


peghmodel_1 <-
  hetsur(
    lhs_vars = peglhs,
    rhs_vars = pegrhs,
    num_var = pegnumvar,
    data = ldbmdata,
    idvar = "pegid",
    tvar = "year",
    restrict.R = pegR,
    restrict.z = pegz,
    het = T,
    ac = "ar1",
    lrtol = 0.00000001
  )

peghmodel_2 <- est_share_new(hetest = peghmodel_1, newsharename = "sm")

peghmodel_3 <- prehetsur(peghmodel_2)

peggof <- systemR2(peghmodel_3)


#######################################################################################################################################################
###Format Output
#######################################################################################################################################################


###create output
OutputToTxt = list()

#Coefficient estimate results
OutputToTxt[["Models"]] = peghmodel_1[[2]][[1]]

#Goodness of fit measures
OutputToTxt[["gof"]] = data.frame("Goodness of Fit Stat" = sapply(names(peggof), function(n){n}),
                                  "Value" = sapply(names(peggof), function(n){peggof[[n]]}))



#Write output to file
sink(file = outPath,
     type = "output")

cat("\n********************************************************************************\n\n")
cat("Benchmarking Estimation Parameters and Results\n\n")
cat(paste0("Date: ", Sys.time(), "\n\n"))
cat("********************************************************************************\n\n")

cat("  Years Used for Parameter Estimation\n")
cat("********************************************************************************\n")
cat("1995-2014")
cat("\n********************************************************************************\n\n")

cat("  REGRESSION WITH GROUPWISE HETEROSKEDASTICITY & AUTOCORRELATION CORRECTION\n")
cat("********************************************************************************\n\n")
print(format(OutputToTxt[["Models"]], justify = "right", digits = 4))
cat("\n")
print.data.frame(format(OutputToTxt[["gof"]], justify = "left"), row.names = FALSE, right = FALSE)
cat("\n********************************************************************************\n\n")



sink()