# Transformation functions with tag and timestamp

## a. Rolling standard deviation (of arbitrary window)
### input: dataset, variable name (it can be a vector of names of several vars) and window size
roll_sd <- function(data, varnames, window){
  fun <- "roll_sd"              ## function tag
  x <- data.frame(1:nrow(data)) ## create an empty df
  for(i in 1:length(varnames)){ ## calculate for one var at once
    a <- data[,varnames[i]]     ## collect data for the var
    calc <- rollapply(a, window, sd) ## apply function
    calc <- as.data.frame(c((rep(NA, window-1)), ## NA fill - keep vec size
                            calc)) 
    colnames(calc) <- paste0(varnames[i], "_rollsd") ## create new var transformed name
    x <- cbind(x, calc)               ## attach it
    timestamp <- format(Sys.time(),   ## create a timestamp
                        "%m %d, %Y %H:%M") 
    track <- t(c(varnames[i], fun, window, timestamp)) ## create a track vec with metadata                                                           info
    write.table(track, file = "trackmetadata.csv", sep = "\t",
                append = TRUE, col.names = FALSE, row.names = FALSE) ## write metainfo in                                                                         a file
  }
  x[,1] <- NULL ## take off x index
  return(x)
}

## b. Rolling mean (of arbitrary window)
### input: dataset, variable name (it can be a vector of names of several vars) and window size
roll_mean <- function(data, varnames, window){
  fun <- "roll_mean" ## tag
  x <- data.frame(1:nrow(data)) ## empty df
  for(i in 1:length(varnames)){ ## for one var at once
    a <- data[,varnames[i]]     ## collect data
    calc <- rollapply(a, window, mean) ## apply transformation
    calc <- as.data.frame(c((rep(NA, window-1)), calc)) ## set vec size
    colnames(calc) <- paste0(varnames[i], "_rollmean") ## create name of new transf var
    x <- cbind(x, calc) ## attach it
    timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## get timestamp
    track <- t(c(varnames[i], fun, window, timestamp)) ## get track info
    write.table(track, file = "trackmetadata.csv", sep = "\t",
                append = TRUE, col.names = FALSE, row.names = FALSE) ## add it in the file
  }
  x[,1] <- NULL
  return(x)
}

## c. Lagging (of arbitrary order)
### input: dataset name (dataframe), varnames (string names) (it can be a vector), lag and dif values (pattern is 1)
lagging <- function(data, varnames, lag=1, dif=1){
  fun <- "lag" ## function tag
  x <- data.frame(1:nrow(data)) ## empty df
  for(i in 1:length(varnames)){ ## for each variable
    a <- data[,varnames[i]] ## collect var data
    calc <- as.data.frame(lag(a, lag=lag, dif=dif)) ## apply transformation
    colnames(calc) <- paste0(varnames[i], "_lag") ## create new var name
    x <- cbind(x, calc) ## attach it
    timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## get timestamp
    track <- t(c(varnames[i], fun, paste0(lag,"_","dif"), timestamp)) ## get track info
    write.table(track, file = "trackmetadata.csv", sep = "\t",
                append = TRUE, col.names = FALSE, row.names = FALSE) ## add it in the file
  }
  x[,1] <- NULL ## remove x indexes
  return(x)
}

## d. Leading (of arbitrary order)

leading <- function(data, varnames, k=1, dif=1){
  fun <- "lead" ## function tag
  x <- data.frame(1:nrow(data)) ## empty df
  for(i in 1:length(varnames)){ ## for var
    a <- data[,varnames[i]] ## collect var data
    calc <- dplyr::lead(a, k=k, dif=dif) ## apply transform
    calc <- as.data.frame(calc)
    colnames(calc) <- paste0(varnames[i], "_lead") ## put new var name
    x <- cbind(x, calc) ## attach it
    timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## get timestamp
    track <- t(c(varnames[i], fun, paste0(k,"_",dif), timestamp)) ## get tracker
    write.table(track, file = "trackmetadata.csv", sep = "\t",
                append = TRUE, col.names = FALSE, row.names = FALSE) ## add it in the file
  }
  x[,1] <- NULL
  return(x)
}

## e. Differencing

diffts <- function(data, varnames, lag=1, dif=1){
  fun <- "diff" ## function tag
  x <- data.frame(1:nrow(data)) ## empty df
  for(i in 1:length(varnames)){ ## for each var
    a <- data[,varnames[i]] ## het var data
    calc <- diff(a, lag=lag, dif=dif) ## apply transformation
    calc <- as.data.frame(c((rep(NA, lag)), calc))
    calc <- apply(calc, 2, as.numeric)
    colnames(calc) <- paste0(varnames[i], "_diff") # create new varname
    x <- cbind(x, calc) ## attach it
    timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## get timestamp
    track <- t(c(varnames[i], fun, window, timestamp)) ## get tracker
    write.table(track, file = "trackmetadata.csv", sep = "\t",
                append = TRUE, col.names = FALSE, row.names = FALSE) ## add it in the file
  }
  x[,1] <- NULL
  return(x)
}

## f. Spread (between two input drivers)

spread <- function(data, var1, var2){
  fun <- "spread" ## function tag
  x <- as.data.frame(data[,var1]-data[,var2]) ## calculate spread
  colnames(x) <- paste0(var1, ".", var2, "_spr") ## new var name
  timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## get timestamp
  track <- t(c(paste0(var1, "_", var2), fun, 0, timestamp)) ## get tracker
  write.table(track, file = "trackmetadata.csv", sep = "\t",
              append = TRUE, col.names = FALSE, row.names = FALSE) ## add tracker in the file
  return(x)
}

## g. Ratio (between two input drivers)

ratio <- function(data, var1, var2){
  fun <- "ratio" ## function tag
  x <- as.data.frame(data[,var1]/data[,var2]) ## get ratio values
  colnames(x) <- paste0(var1, ".", var2, "_rt") ## new var name
  timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## timestamp
  track <- t(c(paste0(var1, "_", var2), fun, 0, timestamp)) ## get tracker
  write.table(track, file = "trackmetadata.csv", sep = "\t",
              append = TRUE, col.names = FALSE, row.names = FALSE) ## add tracker in the file
  return(x)
}


## h. Product (between two input drivers)

prod <- function(data, var1, var2){
  fun <- "prod" ## function tag
  x <- as.data.frame(data[,var1]*data[,var2]) ## get the product values
  colnames(x) <- paste0(var1, ".", var2, "_prod") ## new var name
  timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## timestamp
  track <- t(c(paste0(var1, "_", var2), fun, 0, timestamp)) ## get tracker
  write.table(track, file = "trackmetadata.csv", sep = "\t",
              append = TRUE, col.names = FALSE, row.names = FALSE) ## add tracker in the file
  return(x)
}

## Between diff
betdiff <- function(data, var1, var2){
  fun <- "betdiff" ## fun tag
  x <- as.data.frame(data[,var1]-data[,var2]) ## get diff between 2 vars
  x <- x*1 ## just a transformation to format var type
  colnames(x) <- paste0(var1, ".", var2, "_betdiff") ## new var name
  timestamp <- format(Sys.time(), "%m %d, %Y %H:%M") ## time stamp
  track <- t(c(paste0(var1, "_", var2), fun, 0, timestamp)) ## get tracker
  write.table(track, file = "trackmetadata.csv", sep = "\t",
              append = TRUE, col.names = FALSE, row.names = FALSE) ## add tracker to the file
  return(x)
}

oil[,ncol(oil)+1] <-  betdiff(oil, "OILPRICE", "CL2_log")

head(betdiff(oil, "OILPRICE", "CL2_log") %>%
  roll_sd("OILPRICE.CL2_log_betdiff", 2)) # Note that here it is used the name resulted by                                             previous transformation, following the pattern.

# define function composition
composition <- function(data, list_args){
  x <- data.frame(1:nrow(data))
  for(i in seq(1,length(list_args),2)){
      result <- do.call(list_args[[i]], list_args[[i+1]]) # apply functions chosen
      x <- cbind(x, result)
  }
  x[,1] <- NULL
  data <- cbind(data, x) # attach it to the data
  return(data)
}

lista <- list(betdiff, list(oil, "OILPRICE", "respLAG"),
              lagging, list(oil, "OILPRICE", 2),
              roll_mean, list(oil, c("CL3_log", "CL2_log", "DX1_log"), 2),
              leading, list(oil, "CL15_log", 5))

# apply a composition
head(composition(oil, lista))

# Getting stats values of the tests: Stationarity (KPSS & DF); Normality (Shapiro-Wilk);
# Pearson Correlation with target variable
# input: dataset, name of the variables to be tested, name of the target variable
getStats <- function(data, varnames, target){
  statsTest <- structure(list(variable = character(),
                         test = character(),
                         value = numeric()),
                        class = "data.frame")
  for(i in 1:length(varnames)){
    dat <- replace(x = data[,varnames[i]], list = is.na(data[,varnames[i]]),
                  values = 0)
    # Kwiatkowski-Phillips-Schmidt-Shin - Stacionarity
    kpss <- suppressWarnings(kpss.test(dat)[[3]])
    statsTest <- rbind(statsTest, c(varnames[i], "kpss", kpss))
    # Dickey Fuller - Stacionarity
    dickeyfuller <- suppressWarnings(adf.test(dat))[[4]]
    statsTest <- rbind(statsTest, c(varnames[i], "dickeyfuller",
                                    dickeyfuller))
    # Shapiro-Wilk Test - Normality
    shapiro <- suppressWarnings(shapiro.test(dat))[[2]]
    statsTest <- rbind(statsTest, c(varnames[i], "shapiro", shapiro))
    # Pearson Correlation with target variable
    targetcor <- cor(data[,target], dat)
    statsTest <- rbind(statsTest, c(varnames[i], "targetcor",
                                    targetcor)) ## bind to the                                                                               object
  }
  colnames(statsTest) <- c("var", "test", "value")
  write.table(statsTest, "statsmetadata.txt", sep = "\t",
              col.names = FALSE, append = TRUE)  ## write stats in a file
  return(statsTest)
}

# apply a transformation on data
#oil[,ncol(oil)+1] <-  betdiff(oil, "OILPRICE", "CL2_log")

# get stats for variables
getStats(oil, c("OILPRICE.CL2_log_betdiff", "CL2_log", "CL3_log"), "OILPRICE")

laglist <- list(lagging, list(oil, c("BDIY_log",
                                     "SPX_log", "DX1_log", "GC1_log",
                                     "HO1_log", "USCI_log", "GNR_log",
                                     "SHCOMP_log", "FTSE_log")))

oillag <- composition(oil, laglist)
oillag <- na.omit(oillag) # exclude NA values generated by the lag
oillag[,3:24] <- NULL # excluding CL log columns and indexes not lagged

rm(oil)
invisible(gc())

oilts <- ts(oillag, frequency = 365, start = c(2018, 1, 1)) # random date just to help

autoplot(oilts[,c("OILPRICE", "CL2_log", "BDIY_log_lag", "SPX_log_lag",
                  "DX1_log_lag", "GC1_log_lag", "HO1_log_lag",
                  "USCI_log_lag", "GNR_log_lag", "SHCOMP_log_lag",
                  "FTSE_log_lag", "respLAG")]) +
ggtitle("Plot of the Oil Time-Series") +
theme(plot.title = element_text(hjust = 0.5),
      axis.title.y=element_blank())

# create time slices for the data
oilfinal <- na.omit(oilstationary)

timeSlices <- createTimeSlices(1:nrow(oilfinal), 
                   initialWindow = 700, horizon = 200, fixedWindow = TRUE)

trainSlices <- timeSlices[[1]]
#testSlices <- timeSlices[[2]]

traintransf = as.data.frame(oilfinal[trainSlices[[1]],])
testtransf = as.data.frame(oilfinal[-trainSlices[[1]],])


# original dataset
oilts1 <- oilts[,c("OILPRICE", "CL2_log", "BDIY_log_lag",
                                     "SPX_log_lag", "DX1_log_lag", "GC1_log_lag",
                                     "HO1_log_lag", "USCI_log_lag", "GNR_log_lag",
                                     "SHCOMP_log_lag", "respLAG")]

trainori = as.data.frame(oilts1[trainSlices[[1]],])
testori = as.data.frame(oilts1[-trainSlices[[1]],])

rm(timeSlices, oilts, oilstationary, trainSlices)
invisible(gc())

base_pred <- predict(model_base, newdata = testori, type = "response")
base_cs_pred <- predict(model_base_cs, newdata = testori, type = "response")

sel_pred <- predict(model_sel, newdata = testori, type = "response")
sel_cs_pred <- predict(model_sel_cs, newdata = testori, type = "response")

sel_t_pred <- predict(model_sel_t, newdata = testtransf, type = "response")
sel_t_cs_pred <- predict(model_sel_t_cs, newdata = testtransf, type = "response")

t_pred <- predict(model_t, newdata = testtransf, type = "response")
t_cs_pred <- predict(model_t_cs, newdata = testtransf, type = "response")

accuracy_values <- rbind(
                         accuracy(base_pred, testori[,"OILPRICE"]),
                         accuracy(base_cs_pred, testori[,"OILPRICE"]),
                         accuracy(sel_pred, testori[,"OILPRICE"]),
                         accuracy(sel_cs_pred, testori[,"OILPRICE"]),
                         accuracy(sel_t_pred, testtransf[,"OILPRICE_diff"]),
                         accuracy(sel_t_cs_pred, testtransf[,"OILPRICE_diff"]),
                         accuracy(t_pred, testtransf[,"OILPRICE_diff"]),
                         accuracy(t_cs_pred, testtransf[,"OILPRICE_diff"]))
        
row.names(accuracy_values) <- c("base_pred", "base_cs_pred", "sel_pred",
                                "sel_cs_pred", "sel_t_pred", "sel_t_cs_pred",
                                "t_pred", "t_cs_pred")
accuracy_values

# Limits for 99% of confindence
critval <-  2.80
upr <- sel_pred + (critval * (sel_pred-testori$OILPRICE))
lwr <- sel_pred - (critval * (sel_pred-testori$OILPRICE))

sel_predts <- ts(sel_pred, frequency = 365, start = c(2019,336))
trainorits <- ts(trainori, frequency = 365, start = c(2018,1,1))
testorits <- ts(testori, frequency = 365, start = c(2019,336))
upper <- ts(upr, frequency = 365, start = c(2019,336))
lower <- ts(lwr, frequency = 365, start = c(2019,336))
upr <- c(trainori$OILPRICE, upr)
lwr <- c(trainori$OILPRICE, lwr)
g <- autoplot(oilts1[,"OILPRICE"]) +
  autolayer(sel_predts, color="red") +
#  autolayer(upper) +
#  autolayer(lower) +
  xlab("Year") + ylab("Oil Price") +
  ggtitle("Probabilistic Forecast with Confidence of 99%") +
  guides(colour=guide_legend(title="Probabilistic Forecast for Oil Price"))

g + geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "gray", alpha = 0.6)

# OILPRICE ~ BDIY_log_lag + SPX_log_lag + GC1_log_lag +
#            HO1_log_lag + GNR_log_lag +
#            SHCOMP_log_lag + respLAG,
testing <- data.frame("BDIY_log_lag"=6.579251, 
                        "SPX_log_lag"=7.626482,
                        "GC1_log_lag"=7.166073,
                        "HO1_log_lag"=0.4987733,
                        "GNR_log_lag"=3.771841,
                        "SHCOMP_log_lag"=8.117241,
                        "respLAG"=3.809990)

#calculate the number of model parameters - 1
k <- length(model_sel$coefficients)-1

#calculate sum of squared residuals
SSE <- sum(model_sel$residuals**2)

#calculate total observations in dataset
n <- length(model_sel$residuals)

#calculate residual standard error
stderror <- sqrt(SSE/(n-(1+k)))
conf <- round(1.96*stderror, 2)

print(paste0("The Oil Price tomorrow will be ", round(predict(model_sel, newdata = testing, type = "response"),2), "+/-", conf, " with 95% of confindence. Considering the Mean Absolute Percentage Error, the prediction are different from the real value less than 1%. Then, it will probably be between ", round(predict(model_sel, newdata = testing, type = "response"),2)-round(predict(model_sel, newdata = testing, type = "response"),2)*0.01, " and ",  round(predict(model_sel, newdata = testing, type = "response"),2)+round(predict(model_sel, newdata = testing, type = "response"),2)*0.01))