Time Series Classification Synthetic vs Real Financial Time Series

rm(list = ls())
options(scipen=999)
setwd('C:/Users/Matt/Desktop/Data Science Challenge')
shhh <- suppressPackageStartupMessages

shhh(library(dplyr))
library(readr)
library(TSrepr)
library(ggplot2)
library(data.table)
library(cluster)
library(clusterCrit)
library(fractalrock)
library(cowplot)
library(tidyr)
library(tidyquant)
library(lmtest)
library(aTSA)
library(tsoutliers)
library(tsfeatures)
library(xgboost)
library(caret)
library(purrr)

train_val <- read_csv("train.csv")
test <- read_csv("test.csv")

head(train_val[, 1:5], 1)

## # A tibble: 1 x 5
##   feature1 feature2 feature3 feature4 feature5
##      <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1  0.00629  0.00441  -0.0381   0.0253 -0.00658

colnames(train_val) %>%
  data.frame() %>%
  setNames(c("features")) %>%
  split(as.integer(gl(nrow(.), 20, nrow(.)))) %>%
  kable(caption = "Time series variables") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 12)

######################################################################
################# Clean the data #####################################

# Since the "features" are daily time series, I just choose a random year and rename the feautres into more meaningful names
# Such as "2010-01-01", "2010-01-02", "2010-01-03" instead of "feature1", "feature2", "feature3" etc.
# Theres a "trading dates" package in R to get only the dates which are trading dates.
colnames(train_val) <- getTradingDates('2010-01-01', obs = 260)
colnames(train_val)[ncol(train_val)] <- "class"
colnames(test) <- getTradingDates('2010-01-01', obs = 260)
test$dataset <- "test"
train_val$dataset <- "train"

######################################################################
################# Rearrange the data #################################

# I melt the data for easier analysis, now the data is in a long format.

# "Class" corresponds to whether the asset is Synthetic or Real
# "Dataset" tells me where the data came from
# "row_id" - corresponds to a unique ID assigned to each asset both "(Synthetic & Real)"
# "Variable" is the column names of the original dataset (feature1, feature2, ... , featureN) converted to some date
# "Value" is the daily returns

df <- train_val %>%
  mutate(row_id = row_number()) %>%
  melt(., measure.vars = 1:260) %>%
  arrange(row_id)

head(df)

##   class dataset row_id   variable        value
## 1     0   train      1 2009-01-05  0.006286455
## 2     0   train      1 2009-01-06  0.004407363
## 3     0   train      1 2009-01-07 -0.038088652
## 4     0   train      1 2009-01-08  0.025285012
## 5     0   train      1 2009-01-09 -0.006578773
## 6     0   train      1 2009-01-12  0.005713677

dim(df)

## [1] 3120000       5

# Plot some returns - I only plot a random sample of 20 assets for each Synthetic vs Real.

ret_plot0 <- df %>%
  filter(class == 0) %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>% 
  sample_n(20) %>%
  unnest() %>% 
  ggplot(aes(x = variable, y = value)) +
  geom_line(aes(group = factor(row_id), color = factor(row_id))) +
  ggtitle("Synthetic Financial Time Series") +
  theme_classic() +
  theme(axis.text.x = element_blank(), legend.position = "bottom", legend.title = element_blank())


ret_plot1 <- df %>%
  filter(class == 1) %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>%
  sample_n(20) %>%
  unnest() %>%
  ggplot(aes(x = variable, y = value)) +
  geom_line(aes(group = factor(row_id), color = factor(row_id))) +
  ggtitle("Real Financial Time Series") +
  theme_classic() +
  theme(axis.text.x = element_blank(), legend.position = "bottom", legend.title = element_blank())

plot_grid(ret_plot0, ret_plot1)

ave_box <- df %>%
  group_by(class, row_id) %>%
  summarise(mean = mean(value)) %>%
  ggplot(aes(x = factor(class), y = mean, color = factor(class))) +
  geom_boxplot(show.legend = FALSE) +
  ggtitle("Syn vs Real Average Returns") +
  xlab("Class") +
  ylab("Average Returns") +
  theme_tq()

sd_box <- df %>%
  group_by(class, row_id) %>%
  summarise(sd = sd(value)) %>%
  ggplot(aes(x = factor(class), y = sd, color = factor(class))) +
  geom_boxplot(show.legend = FALSE) +
  ggtitle("Syn vs Real Standard Deviations") +
  xlab("Class") +
  ylab("Standard Deviation") +
  theme_tq()

plot_grid(ave_box, sd_box)

# I calculate the Durbin-Watson statistic and use the "tidy()" function to summarise the key information from the calculation.

dw_test_class_zero <- df %>%
  dplyr::filter(class == 0) %>%
  nest(-row_id) %>%
  mutate(dw_res = map(data, ~ broom::tidy(lmtest::dwtest(value ~ 1, data = .x)))) %>%
  unnest(dw_res) %>%
  mutate(class = "0")

dw_test_class_zero %>% 
  head()

## # A tibble: 6 x 7
##   row_id         data statistic p.value method      alternative            class
##    <int> <list<df[,4>     <dbl>   <dbl> <chr>       <chr>                  <chr>
## 1      1    [260 x 4]      1.98   0.426 Durbin-Wat~ true autocorrelation ~ 0    
## 2      2    [260 x 4]      2.01   0.521 Durbin-Wat~ true autocorrelation ~ 0    
## 3      4    [260 x 4]      2.08   0.747 Durbin-Wat~ true autocorrelation ~ 0    
## 4      5    [260 x 4]      2.49   1.000 Durbin-Wat~ true autocorrelation ~ 0    
## 5      6    [260 x 4]      1.90   0.214 Durbin-Wat~ true autocorrelation ~ 0    
## 6      9    [260 x 4]      1.87   0.138 Durbin-Wat~ true autocorrelation ~ 0

# Here I do the exact same thing as above but this time for the class == 1 data.

dw_test_class_one <- df %>%
  filter(class == 1) %>%
  nest(-row_id) %>%
  mutate(dw_res = map(data, ~ broom::tidy(lmtest::dwtest(value ~ 1, data = .x)))) %>%
  unnest(dw_res) %>%
  mutate(class = "1")

dw_test_class_one %>%
  head()

## # A tibble: 6 x 7
##   row_id         data statistic p.value method      alternative            class
##    <int> <list<df[,4>     <dbl>   <dbl> <chr>       <chr>                  <chr>
## 1      3    [260 x 4]      2.08  0.728  Durbin-Wat~ true autocorrelation ~ 1    
## 2      7    [260 x 4]      1.81  0.0654 Durbin-Wat~ true autocorrelation ~ 1    
## 3      8    [260 x 4]      1.93  0.296  Durbin-Wat~ true autocorrelation ~ 1    
## 4     13    [260 x 4]      2.05  0.644  Durbin-Wat~ true autocorrelation ~ 1    
## 5     15    [260 x 4]      2.07  0.715  Durbin-Wat~ true autocorrelation ~ 1    
## 6     16    [260 x 4]      2.07  0.709  Durbin-Wat~ true autocorrelation ~ 1

# I bind the rows together and plot a box-plot.

bind_rows(dw_test_class_zero, dw_test_class_one) %>%
  group_by(class) %>%
  ggplot(aes(x = factor(class), y = statistic, color = factor(class))) +
  geom_boxplot(show.legend = FALSE) +
  ggtitle("Durbin Watson Box Plot Statistics") +
  xlab("Class") +
  ylab("Durbin Watson") +
  theme_tq()

# Rolling mean and standard deviations
# I only use a random sample of 1 of each class of the grouped observations to save on memory and to make the plot more readable.
# The rollowing window is 10 days
# I use the tq_mutate functionality from the "tidyquant" package to keep things in a "tidy" format as per the "tidyverse" 'rules'.
# In the plot "value" is the returns, "mean_10" is the 10 day rolling mean and "sd_10" is the 10 day rolling standard deviation.

plot0 <- df %>%
  filter(class == 0) %>%
  as_tibble() %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>% 
  sample_n(1) %>%
  unnest() %>%
  mutate(variable = as.Date(variable)) %>%
  tq_mutate(
    select     = value,
    mutate_fun = rollapply,
    width      = 10,
    align      = "right",
    FUN        = mean,
    na.rm      = TRUE,
    col_rename = "mean_10"
    ) %>%
  tq_mutate(
    select     = value,
    mutate_fun = rollapply,
    width      = 10,
    align      = "right",
    FUN        = sd,
    na.rm      = TRUE,
    col_rename = "sd_10"
    ) %>%
  melt(measure.vars = 5:7) %>%
  setNames(c("row_id", "class", "data set", "date", "variable", "value")) %>%
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = variable)) +
  ggtitle("Synthetic Financial Time Series Rolling Mean and Standard Deviation") +
  theme_classic() +
  scale_colour_manual(values = c("#1f77b4", "red", "black")) +
  theme(axis.text.x = element_blank(), legend.position = "bottom", legend.title = element_blank())

plot1 <- df %>%
  filter(class == 1) %>%
  as_tibble() %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>% 
  sample_n(1) %>%
  unnest() %>%
  mutate(variable = as.Date(variable)) %>%
  tq_mutate(
    select     = value,
    mutate_fun = rollapply,
    width      = 10,
    align      = "right",
    FUN        = mean,
    na.rm      = TRUE,
    col_rename = "mean_10"
  ) %>%
  tq_mutate(
    select     = value,
    mutate_fun = rollapply,
    width      = 10,
    align      = "right",
    FUN        = sd,
    na.rm      = TRUE,
    col_rename = "sd_10"
  ) %>%
  melt(measure.vars = 5:7) %>%
  setNames(c("row_id", "class", "data set", "date", "variable", "value")) %>%
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = variable)) +
  ggtitle("Real Financial Time Series Rolling Mean and Standard Deviation") +
  theme_classic() +
  scale_colour_manual(values = c("#1f77b4", "red", "black")) +
  theme(axis.text.x = element_blank(), legend.position = "bottom", legend.title = element_blank())  

plot_grid(plot0, plot1)

#  group_by(row_id) %>%
#  nest() %>%
#  ungroup() %>% 
#  sample_n(1) %>%
#  unnest() %>%

# Dickey Fuller test on the 0 class
# I only randomly sample 1 of the assets for the 0 class to save on output space

df %>%
  filter(class == 0) %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>% 
  sample_n(1) %>%
  unnest() %>%
  nest(-row_id) %>%
  mutate(adf_res = map(data, ~ adf.test(.x$value))) %>%
  unnest(adf_res)

## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag    ADF p.value
## [1,]   0 -17.94    0.01
## [2,]   1 -11.75    0.01
## [3,]   2  -8.66    0.01
## [4,]   3  -7.62    0.01
## [5,]   4  -7.13    0.01
## Type 2: with drift no trend 
##      lag    ADF p.value
## [1,]   0 -17.94    0.01
## [2,]   1 -11.76    0.01
## [3,]   2  -8.67    0.01
## [4,]   3  -7.64    0.01
## [5,]   4  -7.15    0.01
## Type 3: with drift and trend 
##      lag    ADF p.value
## [1,]   0 -18.00    0.01
## [2,]   1 -11.83    0.01
## [3,]   2  -8.77    0.01
## [4,]   3  -7.74    0.01
## [5,]   4  -7.26    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

## # A tibble: 3 x 3
##   row_id           data adf_res          
##    <int> <list<df[,4]>> <named list>     
## 1   7807      [260 x 4] <dbl[,3] [5 x 3]>
## 2   7807      [260 x 4] <dbl[,3] [5 x 3]>
## 3   7807      [260 x 4] <dbl[,3] [5 x 3]>

# Dickey Fuller test on the 1 class
# I only randomly sample 1 of the assets for the 1 class to save on output space

df %>%
  filter(class == 1) %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>% 
  sample_n(1) %>%
  unnest() %>%
  nest(-row_id) %>%
  mutate(adf_res = map(data, ~ adf.test(.x$value))) %>%
  unnest(adf_res)

## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag    ADF p.value
## [1,]   0 -15.99    0.01
## [2,]   1 -10.71    0.01
## [3,]   2  -9.12    0.01
## [4,]   3  -8.74    0.01
## [5,]   4  -7.58    0.01
## Type 2: with drift no trend 
##      lag    ADF p.value
## [1,]   0 -16.10    0.01
## [2,]   1 -10.84    0.01
## [3,]   2  -9.27    0.01
## [4,]   3  -8.93    0.01
## [5,]   4  -7.81    0.01
## Type 3: with drift and trend 
##      lag    ADF p.value
## [1,]   0 -16.27    0.01
## [2,]   1 -10.99    0.01
## [3,]   2  -9.46    0.01
## [4,]   3  -9.18    0.01
## [5,]   4  -8.06    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

## # A tibble: 3 x 3
##   row_id           data adf_res          
##    <int> <list<df[,4]>> <named list>     
## 1  10833      [260 x 4] <dbl[,3] [5 x 3]>
## 2  10833      [260 x 4] <dbl[,3] [5 x 3]>
## 3  10833      [260 x 4] <dbl[,3] [5 x 3]>

# For both classes I take a random sample of 1 observation from each class (Synthetic and Real financial series)

jb_zero <- df %>%
  filter(class == 0) %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>% 
  sample_n(1) %>%
  unnest() %>%
  nest(-row_id) %>%
  mutate(JarqueBeraTest = map(data, ~ JarqueBera.test(.x$value)))

print("Jarque-Bera Test on the 0 - Synthetic class")

## [1] "Jarque-Bera Test on the 0 - Synthetic class"

jb_zero$JarqueBeraTest

## [[1]]
## 
##  Jarque Bera Test
## 
## data:  .x$value
## X-squared = 0.3088, df = 2, p-value = 0.8569
## 
## 
##  Skewness
## 
## data:  .x$value
## statistic = 0.045794, p-value = 0.7631
## 
## 
##  Kurtosis
## 
## data:  .x$value
## statistic = 2.8582, p-value = 0.6406

jb_one <- df %>%
  filter(class == 0) %>%
  group_by(row_id) %>%
  nest() %>%
  ungroup() %>% 
  sample_n(1) %>%
  unnest() %>%
  nest(-row_id) %>%
  mutate(JarqueBeraTest = map(data, ~ JarqueBera.test(.x$value)))

print("Jarque-Bera Test on the 1 - Real class")

## [1] "Jarque-Bera Test on the 1 - Real class"

jb_one$JarqueBeraTest

## [[1]]
## 
##  Jarque Bera Test
## 
## data:  .x$value
## X-squared = 25.14, df = 2, p-value = 0.000003474
## 
## 
##  Skewness
## 
## data:  .x$value
## statistic = 0.084191, p-value = 0.5794
## 
## 
##  Kurtosis
## 
## data:  .x$value
## statistic = 4.514, p-value = 0.0000006251

######################################################################
################# ACF plots ##########################################

# I only use 4 observations for these plots, 2 from the "synthetic" class and 2 from the "real" class.

df %>%
  filter(row_id == 6422 | row_id == 8967 | row_id == 6080 | row_id ==   5734) %>%
  mutate(date = as.Date(variable)) %>%
  ggplot(aes(x = date)) +
  geom_line(aes(y = value), color = "red", alpha = 0.4) +
  geom_hline(yintercept = 0) +
  facet_wrap(~ row_id + class) +
  theme_tq()

acf_data <- df %>%
  dplyr::filter(row_id == 6422 | row_id == 8967 | row_id == 6080 | row_id ==    5734) %>%
  mutate(date = as.Date(variable))

df_acf <- acf_data %>%
  group_by(row_id) %>% 
  summarise(list_acf = list(acf(value, plot=FALSE))) %>%
  mutate(acf_vals = purrr::map(list_acf, ~as.numeric(.x$acf))) %>% 
  select(-list_acf) %>% 
  unnest() %>% 
  group_by(row_id) %>% 
  mutate(lag = row_number() - 1)

df_ci <- acf_data %>% 
  group_by(row_id) %>% 
  summarise(ci = qnorm((1 + 0.95)/2)/sqrt(n()))

ggplot(df_acf, aes(x = lag, y = acf_vals)) +
  geom_bar(stat="identity", width=.05) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color="blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color="blue", linetype="dotted") +
  labs(x="Lag", y="ACF") +
  facet_wrap(~ row_id) +
  theme_tq()

################# Generate Time Series Features ######################

# I create some time series features from the package "tsfeatures". There are 40+ functions in the "tsfeatures" package
# which can generate approximately 106 time series features.
# Due to memory issues I am only able to create a few of the features, therefore I randomly sample 10 features from the
# "tsfeatures" package. We could also add in technical indicators from the "PerformanceAnalytics" or "TTR" packages (I omit these
# here, however creating 'functions2 <- ls("package:TTR")' and adding it to the 'summarise' command will work.)

functions <- ls("package:tsfeatures")[1:42]
# functions <- sample(functions, 20)

Stats <- df %>%
  group_by(row_id, class) %>%
  nest() %>%
  ungroup() %>%
  sample_n(5) %>%
  unnest() %>%
  nest(-row_id, -class) %>%
  group_by(row_id, class) %T>%
  {options(warn = -1)} %>%
  summarise(Statistics = map(data, ~ data.frame(
    bind_cols(
      tsfeatures(.x$value, functions))))) %>%
  unnest(Statistics)

# I saved to whole dataset as "Stats" next I split it between training and test.
Stats <- read.csv("C:/Users/Matt/Desktop/Data Science Challenge/TSfeatures_train_val.csv")

## [1] 12000   109

######################################################################
################# Train and XGBoost model on the TS Features #########

#Stats <- Stats %>%
#  select_if(~sum(!is.na(.)) > 0)

# Split the training set up between train and a small validation set
smp_size <- floor(0.75 * nrow(Stats))
#set.seed(123)
train_ind <- sample(seq_len(nrow(Stats)), size = smp_size)

train <- Stats[train_ind, ]
val <- Stats[-train_ind, ]

# We have 106 time series features for the model to learn from.

x_train <- train %>%
  ungroup() %>%
  select(-class, -row_id, -X) %>%
  as.matrix()

x_val <- val %>%
  ungroup() %>%
  select(-class, -row_id, -X) %>%
  as.matrix()

y_train <- train %>%
  ungroup() %>%
  pull(class)

y_val <- val %>%
  ungroup() %>%
  pull(class)

######################################################################
################# XGBoost Grid Search to locate Optimal Parameters ###

##############################################################################################################################
# NOTE: This section was taken from the first chapter of my PhD where I needed to search over a parameter space to locate the
# most optimal parameters - I have just adapted it for this problem of Time Series Classification.
# Its simple enough to add parameters and different values - I just optimise a few important parameters from domain knowledge
# of the XGBoost model for this task, i.e depth and eta are quite important in gradient boosting.

# 1) I create a "grid" with different parameter values or combinations of parameter values
# 2) I apply cross validation over the parameter space to fine the most optimal values for the XGBoost model.
# 3) I print the model parameters which give the best train / (in-sample test) results in a data table.
##############################################################################################################################

# Grid Search Parameters:
# 1)
searchGridSubCol <- expand.grid(subsample = c(1), #Range (0,1], default = 1, set to 0.5 will prevent overfitting
                                colsample_bytree = c(1), #Range (0,1], default = 1
                                max_depth = c(5, 8, 14), #Range (0, inf], default = 6
                                min_child = c(1), #Range (0, inf], default = 1
                                eta = c(0.1, 0.05, 0.3), #Range (0,1], default = 0.3
                                gamma = c(0), #Range (0, inf], default = 0
                                lambda = c(1), #Default = 1, L2 regularisation on weights, higher the more conservative the model
                                alpha = c(0), #Default = 0, L1 regularisation on weights, higher the more conservative the model
                                max_delta_step = c(0), #Range (0, inf], default = 0 (Helpful for logisitc regression when class is extremely imbalanced, set to value 1-10 may help control the update)
                                colsample_bylevel = c(1) #Range (0,1], default = 1
                                )

ntrees = 200
nfold <- 10                             # I use nfold = 10 which is probably too many folds, 5 should be sufficient.
watchlist <- list(train = dtrain, test = dval)

# 2)
system.time(
  AUCHyperparameters <- apply(searchGridSubCol, 1, function(parameterList){
    #Extract Parameters to test
    currentSubsampleRate <- parameterList[["sub_sample"]]
    currentColsampleRate <- parameterList[["colsample_bytree"]]
    currentDepth <- parameterList[["max_depth"]]
    currentEta <- parameterList[["eta"]]
    currentMinChild <- parameterList[["min_child"]]
    gamma <- parameterList[["gamma"]]
    lambda <- parameterList[["lambda"]]
    alpha <- parameterList[["alpha"]]
    max_delta_step <- parameterList[["max_delta_step"]]
    colsample_bylevel <- parameterList[["colsample_bylevel"]]
    xgboostModelCV <- xgb.cv(data =  dtrain,
                             nrounds = ntrees,
                             nfold = nfold,
                             showsd = TRUE,
                             metrics = c("auc", "logloss", "error"),
                             verbose = TRUE,
                             "eval_metric" = c("auc", "logloss", "error"),
                             "objective" = "binary:logistic", #Outputs a probability "binary:logitraw" - outputs score before logistic transformation
                             "max.depth" = currentDepth,
                             "eta" = currentEta,
                             "gamma" = gamma,
                             "lambda" = lambda,
                             "alpha" = alpha,
                             "subsample" = currentSubsampleRate,
                             "colsample_bytree" = currentColsampleRate,
                             print_every_n = 50, # print ever 50 trees to reduce the outputs printed.
                             "min_child_weight" = currentMinChild,
                             booster = "gbtree", #booster = "dart"  #using dart can help improve accuracy.
                             early_stopping_rounds = 10,
                             watchlist = watchlist,
                             seed = 1234)
    xvalidationScores <<- as.data.frame(xgboostModelCV$evaluation_log)
    train_auc_mean <- tail(xvalidationScores$train_auc_mean, 1)
    test_auc_mean <- tail(xvalidationScores$test_auc_mean, 1)
    train_logloss_mean <- tail(xvalidationScores$train_logloss_mean, 1)
    test_logloss_mean <- tail(xvalidationScores$test_logloss_mean, 1)
    train_error_mean <- tail(xvalidationScores$train_error_mean, 1)
    test_error_mean <- tail(xvalidationScores$test_error_mean, 1)
    output <- return(c(train_auc_mean, test_auc_mean, train_logloss_mean, test_logloss_mean, train_error_mean, test_error_mean, xvalidationScores, currentSubsampleRate, currentColsampleRate, currentDepth, currentEta, gamma, lambda, alpha, max_delta_step, colsample_bylevel, currentMinChild))
    hypemeans <- which.max(AUCHyperparameters[[1]]$test_auc_mean)
    output2 <- return(hypemeans)
    }))

# 3)
output <- as.data.frame(t(sapply(AUCHyperparameters, '[', c(1:6, 20:29))))
varnames <- c("TrainAUC", "TestAUC", "TrainLogloss", "TestLogloss", "TrainError", "TestError", "SubSampRate", "ColSampRate", "Depth", "eta", "gamma", "lambda", "alpha", "max_delta_step", "col_sample_bylevel", "currentMinChild")
colnames(output) <- varnames
data.table(output)

#################################################################################
################# XGBoost Optimal Parameters from Cross Validation ##############

# This is the final training model where I use the most optimal parameters found over the grid space and plug them in here.

watchlist <- list("train" = dtrain)

params <- list("eta" = 0.1, "max_depth" = 5, "colsample_bytree" = 1, "min_child_weight" = 1, "subsample"= 1,
               "objective"="binary:logistic", "gamma" = 1, "lambda" = 1, "alpha" = 0, "max_delta_step" = 0,
               "colsample_bylevel" = 1, "eval_metric"= "auc",
               "set.seed" = 176)

nround <- 95

# Train the XGBoost model

xgb.model <- xgb.train(params, dtrain, nround, watchlist)

## [1]  train-auc:0.700790 
## [2]  train-auc:0.720114 
## [3]  train-auc:0.735281 
## [4]  train-auc:0.741159 
## [5]  train-auc:0.748016 
## [6]  train-auc:0.752070 
## [7]  train-auc:0.754637 
## [8]  train-auc:0.759151 
## [9]  train-auc:0.762538 
## [10] train-auc:0.769652 
## [11] train-auc:0.776582 
## [12] train-auc:0.780015 
## [13] train-auc:0.782065 
## [14] train-auc:0.782815 
## [15] train-auc:0.788966 
## [16] train-auc:0.791026 
## [17] train-auc:0.793545 
## [18] train-auc:0.797363 
## [19] train-auc:0.799069 
## [20] train-auc:0.802015 
## [21] train-auc:0.802583 
## [22] train-auc:0.806938 
## [23] train-auc:0.808239 
## [24] train-auc:0.811255 
## [25] train-auc:0.813142 
## [26] train-auc:0.816767 
## [27] train-auc:0.817697 
## [28] train-auc:0.820239 
## [29] train-auc:0.821589 
## [30] train-auc:0.823343 
## [31] train-auc:0.823939 
## [32] train-auc:0.825701 
## [33] train-auc:0.827316 
## [34] train-auc:0.829365 
## [35] train-auc:0.832646 
## [36] train-auc:0.833297 
## [37] train-auc:0.837006 
## [38] train-auc:0.838857 
## [39] train-auc:0.839923 
## [40] train-auc:0.842968 
## [41] train-auc:0.844877 
## [42] train-auc:0.845940 
## [43] train-auc:0.846583 
## [44] train-auc:0.847330 
## [45] train-auc:0.848292 
## [46] train-auc:0.850215 
## [47] train-auc:0.851641 
## [48] train-auc:0.852670 
## [49] train-auc:0.854706 
## [50] train-auc:0.855752 
## [51] train-auc:0.856772 
## [52] train-auc:0.857806 
## [53] train-auc:0.860245 
## [54] train-auc:0.861337 
## [55] train-auc:0.864178 
## [56] train-auc:0.865290 
## [57] train-auc:0.865808 
## [58] train-auc:0.866386 
## [59] train-auc:0.867751 
## [60] train-auc:0.870032 
## [61] train-auc:0.870500 
## [62] train-auc:0.872442 
## [63] train-auc:0.873391 
## [64] train-auc:0.875188 
## [65] train-auc:0.877767 
## [66] train-auc:0.879196 
## [67] train-auc:0.880079 
## [68] train-auc:0.879969 
## [69] train-auc:0.880638 
## [70] train-auc:0.881389 
## [71] train-auc:0.882066 
## [72] train-auc:0.882515 
## [73] train-auc:0.883854 
## [74] train-auc:0.884654 
## [75] train-auc:0.885104 
## [76] train-auc:0.885922 
## [77] train-auc:0.887100 
## [78] train-auc:0.888646 
## [79] train-auc:0.889833 
## [80] train-auc:0.890387 
## [81] train-auc:0.891815 
## [82] train-auc:0.892281 
## [83] train-auc:0.894417 
## [84] train-auc:0.895006 
## [85] train-auc:0.897079 
## [86] train-auc:0.899254 
## [87] train-auc:0.901114 
## [88] train-auc:0.902460 
## [89] train-auc:0.902939 
## [90] train-auc:0.903763 
## [91] train-auc:0.903792 
## [92] train-auc:0.904433 
## [93] train-auc:0.904986 
## [94] train-auc:0.907339 
## [95] train-auc:0.907761

# Note: Plot AUC on for the in-sample train / validation scores -  this was a note for me at the time of writing this R file - I never did get around to plotting the AUC for the in-sample train / validation scores...

# We can obtain "feature" importance results from the model.
xgb.imp <- xgb.importance(model = xgb.model)
xgb.plot.importance(xgb.imp, top_n = 10)

# I next make the predictions on the 'in-sample' held out test set, that is, originally I took the 12,000 training samples
# and split them between 75% training and 25% 'in-sample' testing (9000 training vs 3000 in-sample testing)

# I plot the probabilities from the model - the "dashed" line is the average predicted probability.
xgb.pred <- predict(xgb.model, dval, type = 'prob')

results <- cbind(y_val, xgb.pred)

results %>%
 as.tibble() %>%
 ggplot(aes(x = xgb.pred)) + 
 geom_density(color = "darkblue", fill = "lightblue") +
 geom_vline(aes(xintercept = mean(xgb.pred)),
            color = "blue", linetype = "dashed", size = 1) +
 geom_histogram(aes(y = ..density..), colour = "black", fill = "white", alpha = 0.1, position = "identity") +
 ggtitle("Predicted probability density plot") +
 theme_tq()

# The average predicted probability sits around 0.48 / 0.49, for simplicity I will just select 0.50 as the cut off threshold.
# That is, all observations <= 0.50 are assigned a "0" class or "synthetic" data and all observations >= are assigned a "1" or
# "real" data.
# Finally I output the confusion matrix on the 'in-sample' testing data.

results <- results %>%
  as_tibble() %>%
  mutate(pred = case_when(
    xgb.pred > 0.5 ~ 1,
    xgb.pred <= 0.5 ~ 0
  ))

confusionMatrix(factor(results$pred), factor(results$y_val))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1041  537
##          1  465  957
##                                              
##                Accuracy : 0.666              
##                  95% CI : (0.6488, 0.6829)   
##     No Information Rate : 0.502              
##     P-Value [Acc > NIR] : <0.0000000000000002
##                                              
##                   Kappa : 0.3319             
##                                              
##  Mcnemar's Test P-Value : 0.0249             
##                                              
##             Sensitivity : 0.6912             
##             Specificity : 0.6406             
##          Pos Pred Value : 0.6597             
##          Neg Pred Value : 0.6730             
##              Prevalence : 0.5020             
##          Detection Rate : 0.3470             
##    Detection Prevalence : 0.5260             
##       Balanced Accuracy : 0.6659             
##                                              
##        'Positive' Class : 0                  
## 

 test_final <- read.csv("C:/Users/Matt/Desktop/Data Science Challenge/test.csv") %>%
   mutate(row_id = row_number()) %>%
   melt(., measure.vars = 1:260) %>%
   arrange(row_id)

functions <- sample(functions, 20)

test_final <- test_final %>%
  group_by(row_id) %>%
#  nest() %>%
#  sample_n(5) %>%
#  ungroup() %>%
#  unnest() %>%
  nest(-row_id) %>%
  group_by(row_id) %T>%
  {options(warn = -1)} %>%
  summarise(Statistics = map(data, ~ data.frame(
    bind_cols(
      tsfeatures(.x$value, functions))))) %>%
  unnest(Statistics)

#print("Generated 106 Time Series features")
#write.csv(test_final, "TSfeatures_test.csv")

# I have already created the features for the training dataset so I can just load them right back in as 
train_final <- read.csv("C:/Users/Matt/Desktop/Data Science Challenge/TSfeatures_train_val.csv")
test_final <- read.csv("C:/Users/Matt/Desktop/Data Science Challenge/TSfeatures_test.csv")

train_final %>%
  head() %>%
  kable(caption = "Final training data set") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 12)

test_final %>%
  head() %>%
  kable(caption = "Final testing data set") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 12)

# previously and run the final training model (to make predictions on the out-of-sample test data)

x_train_final <- train_final %>%
  ungroup() %>%
  select(-class, -row_id, -X) %>%
  as.matrix()

x_test_final <- test_final %>%
  ungroup() %>%
  select(-row_id, -X) %>%
  as.matrix()

y_train_final <- train_final %>%
  ungroup() %>%
  pull(class)

dtrain_final <- xgb.DMatrix(data = as.matrix(x_train_final), label = y_train_final, missing = "NaN")
dtest_final <- xgb.DMatrix(data = as.matrix(x_test_final), missing = "NaN")

watchlist <- list("train" = dtrain_final)

params <- list("eta" = 0.1, "max_depth" = 5, "colsample_bytree" = 1, "min_child_weight" = 1, "subsample"= 1,
               "objective"="binary:logistic", "gamma" = 1, "lambda" = 1, "alpha" = 0, "max_delta_step" = 0,
               "colsample_bylevel" = 1, "eval_metric"= "auc",
               "set.seed" = 176)

nround <- 95

xgb.model_final <- xgb.train(params, dtrain_final, nround, watchlist)

## [1]  train-auc:0.708604 
## [2]  train-auc:0.721700 
## [3]  train-auc:0.723230 
## [4]  train-auc:0.729888 
## [5]  train-auc:0.735542 
## [6]  train-auc:0.738081 
## [7]  train-auc:0.740926 
## [8]  train-auc:0.744105 
## [9]  train-auc:0.746320 
## [10] train-auc:0.748644 
## [11] train-auc:0.754211 
## [12] train-auc:0.756892 
## [13] train-auc:0.761524 
## [14] train-auc:0.763882 
## [15] train-auc:0.767216 
## [16] train-auc:0.772009 
## [17] train-auc:0.772943 
## [18] train-auc:0.774261 
## [19] train-auc:0.775471 
## [20] train-auc:0.777801 
## [21] train-auc:0.780629 
## [22] train-auc:0.784384 
## [23] train-auc:0.787112 
## [24] train-auc:0.788946 
## [25] train-auc:0.791835 
## [26] train-auc:0.793142 
## [27] train-auc:0.795289 
## [28] train-auc:0.798502 
## [29] train-auc:0.799893 
## [30] train-auc:0.802186 
## [31] train-auc:0.804981 
## [32] train-auc:0.805649 
## [33] train-auc:0.807120 
## [34] train-auc:0.809020 
## [35] train-auc:0.810318 
## [36] train-auc:0.812637 
## [37] train-auc:0.814760 
## [38] train-auc:0.816024 
## [39] train-auc:0.817956 
## [40] train-auc:0.819350 
## [41] train-auc:0.821653 
## [42] train-auc:0.822729 
## [43] train-auc:0.824029 
## [44] train-auc:0.824765 
## [45] train-auc:0.826924 
## [46] train-auc:0.827804 
## [47] train-auc:0.828475 
## [48] train-auc:0.831018 
## [49] train-auc:0.832247 
## [50] train-auc:0.833265 
## [51] train-auc:0.834168 
## [52] train-auc:0.835535 
## [53] train-auc:0.836093 
## [54] train-auc:0.837008 
## [55] train-auc:0.837715 
## [56] train-auc:0.839537 
## [57] train-auc:0.840310 
## [58] train-auc:0.841701 
## [59] train-auc:0.842480 
## [60] train-auc:0.843106 
## [61] train-auc:0.844495 
## [62] train-auc:0.845348 
## [63] train-auc:0.845932 
## [64] train-auc:0.847843 
## [65] train-auc:0.849445 
## [66] train-auc:0.850345 
## [67] train-auc:0.851337 
## [68] train-auc:0.852121 
## [69] train-auc:0.852663 
## [70] train-auc:0.854132 
## [71] train-auc:0.855949 
## [72] train-auc:0.856758 
## [73] train-auc:0.857115 
## [74] train-auc:0.857954 
## [75] train-auc:0.858849 
## [76] train-auc:0.859527 
## [77] train-auc:0.859917 
## [78] train-auc:0.860590 
## [79] train-auc:0.861264 
## [80] train-auc:0.862359 
## [81] train-auc:0.863101 
## [82] train-auc:0.863794 
## [83] train-auc:0.864911 
## [84] train-auc:0.866293 
## [85] train-auc:0.866976 
## [86] train-auc:0.867436 
## [87] train-auc:0.869036 
## [88] train-auc:0.869469 
## [89] train-auc:0.869931 
## [90] train-auc:0.870681 
## [91] train-auc:0.872326 
## [92] train-auc:0.873706 
## [93] train-auc:0.875704 
## [94] train-auc:0.876178 
## [95] train-auc:0.876789

# Make the final predictions on the 'test.csv' data and plot the probability density function.

xgb.pred_final <- predict(xgb.model_final, dtest_final, type = 'prob')

xgb.pred_final %>%
 as_tibble() %>%
 setNames(c("Prediction")) %>%
 ggplot(aes(x = Prediction)) + 
 geom_density(color = "darkblue", fill = "lightblue") +
 geom_vline(aes(xintercept = mean(Prediction)),
            color = "blue", linetype = "dashed", size = 1) +
 geom_histogram(aes(y = ..density..), colour = "black", fill = "white", alpha = 0.1, position = "identity") +
 ggtitle("(Out of sample) - Predicted probability density plot") +
 theme_tq()

# Convert the probabilities into a binary class of 0 or 1 by a decision threshold of 0.465.
# Write the predictions to "submission.csv"

xgb.pred_final %>%
  as_tibble() %>%
  setNames(c("Prediction")) %>%
  summarise(mean = mean(Prediction))

## # A tibble: 1 x 1
##    mean
##   <dbl>
## 1 0.465

xgb.pred_final %>%
  as_tibble() %>%
  setNames(c("Prediction")) %>%
  mutate(pred = case_when(
    Prediction > 0.465 ~ 1,
    Prediction <= 0.465 ~ 0
  )) %>%
  write.csv("submission.csv")