A brief overview of the notebook:

Part 1 of the notebook:

• Cleans the data and puts it into a better format for analysis. The data I recieved removed all dates, assest names etc. for anonymity.
• Simple plot of some returns for the Synthetic and Real financial time series.
• Box-plots of average returns and standard deviations.
• Computes the Durbin-Watson test statistics for both Synthetic and Real time series and box-plots.
• Plot the 10 day rolling mean and standard deviations for a random time series for Synthetic and real data.
• Dickey Fuller test on both the Synthetic and real time series.
• Jarque-Bera Test For Normality on the Synthetic and real time series.
• ACF Plots for both the Synthetic and real time series.

Part 2 of the notebook:

• Creates the time series features.
• Splits the train.csv into “train” and “validation” data sets.
• Puts the data into the correct format for XGBoost.
• Sets up and searches over a parameter space to find the most optimal parameters for this data set (on the train data).
• Outputs these parameters into a data frame.
• Train the model using the optimal parameters found from the grid-search.
• Plot the feature importance scores - i.e. the most “important” variables that the model found when making its predictions.
• Assign a cut-off on the probability scores (> 0.5 then assign a 1 - real time series, <= 0.5 then assign a 0 for Synthetic).
• Compute the Confusion Matrix and analyse the ‘in-sample’ validation results.

Part 3 of the notebook:

• Create the “test.csv” features just as before and save as “TSfeatures_test.csv”.
• Load in the “TSfeatures_train_val.csv” and “TSfeatures_test.csv” which were created from “train.csv” and “test.csv”.
• Set up and run the XGBoost model using the optimal parameters found from the cross-validation grid search in “Part 2”.
• Plot the predicted probability density plot as before in “Part 2”.
• Set the cut-off threshold as the mean prediction score (0.465) which is close to the (0.500) score from “Part 2”.
• Save the results as “submission.csv”.

Lets get started…

I often remove all other data in my environment before hand and turn scientific notation off which is what the first 2 lines does. The shhh command is useful for Jupyter Notebooks which outputs all the warning messages, adding shhh suppresses these warning messaged when loading in the packages. (In R markdown I can set warning = FALSE but there is no option on Notebooks. - that I know of - )

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")

NOTE:

I have 2 data sets, the train_Val.csv for training and validation data set and the test.csv data set. I do not touch the test.csv data set until the very end in part 3. All the analysis and optimisation is performed only on the train_val.csv data set. The train_val.csv contains 12,000 observations and the test.csv contains 12,000 observations.

Part 1

The data was given to me in this format:

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

The names of the columns are as follows:

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)

There are 260 “features” in the train data along with a class variable which is excluded from the testing data. With ~253 trading days in a year the feature1, feature2, … featureN were daily time series. From my initial observation (and plots) I believed this data to be “returns” data. I firstly clean a little the data since time series does not do so well with feature1, feature2, … featureN as its input. I chose a year at random and renamed the columns with the function getTradingDates (there is always an R package for everything…).

######################################################################
################# 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"

Here (if I were to do things differently) I would keep to tidy data principles and use test %>% add_column(dataset = "test) and train %>% add_colum(dataset = "train") instead of test$dataset <- "test and train_val$dataset <- "train". But that doesn’t matter much.

How the training data looks after cleaning:

How the testing data looks after cleaning:

The goal: Was to classify which financial time series were real vs which were synthetically created (by some algorithm I have no knowledge of how it generated the synthetic time series)

I re-arranged the data using the melt function in R, however I suggest anybody reading this to use the pivol_longer function from the tidyverse packages. The pivot_longer package was released a few weeks after writing the code for this problem.

######################################################################
################# 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

Note: I call the training data df which in hindsight is probably bad practice and it should be called something related to the train_Val named data set. Just keep in mind that df refers to the train_Val data set. (and does not include the test.csv data set data)

As we can see the data has 3,120,000 rows which is 12,000 assets * 260 trading days. Next I plot the returns series using ggplot.

# 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)

Next I plot boxplots for the Average returns and secondly the standard deviations.

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 next calculate the Durbin-Watson statistic. I mostly code using R’s tidy data principles and therefore use the tidy function from the broom package to tidy the output of the DW statistic up a little. I do this for both the synthetic time series and real time series.

# 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

Next I plot the boxplot statistics for each of the Durbin Watson tests.

# 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()

I compute the 10 day rolling mean and standard deviation using the tq_mutate function from the tidyquant package. value corresponds to the returns of the financial time series and is plotted in blue with the 10 day rolling average and standard deviation plotted over the returns. (I use melt again here but look into the pivot_longer function for a more intuitive application)

# 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)

An important note in the code here is that I randomly sample by group, that is, I do not take a random sample from all observations across all groups. Instead I group_by each time series (each of the 6,000 observations after I filtered by the class == 0, likewise when I filter by the class == 1), I then nest() the data to collapse the daily time series for each asset into a list. From here I will have 6,000 observations, each of which has their time series nested inside a list. Thus, I can sample 1 of the 6,000 observations and then unnest() and obtain a full time series set of one of the random assets selected, - instead of sampling randomly over all assets time series data (which would be completely wrong).

For example the following commented out code group_by() the ID variable and nest() the data, takes a random sample_n() of the grouped data and then unnest() the data to its original form, this time with a random sample of the IDs.

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

Next I compute the Dickey Fuller test on both series for a single random observation, hence the sample_n(1) argument (it’s too computationally expensive to compute it on all 12,000 observations).

For the synthetically created series.

# 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]>

The same but on the real financial series.

# 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]>

Next the Jarque-Bera tests for normality. Firstly on the synthetically created series.

# 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

Also on the real financial series.

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

Autocorrelation plots:

I plot the Autocorrelation Function for a “random” sample of observations time series. I selected 4 observations and filtered the data by them.

######################################################################
################# 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()

Thats enough data analysis I could probably fit the PACF plots also along with a few more exploratory data analysis but I move on to generating the financial time series features using the tsfeatures package.

What I do in the below code is to take a random sample of 5 groups (Using the whole data set takes too long to calculate the time series features) and then apply all the functions in the tsfeatures package to each of the time series assets data which is does by mapping over each assets data and computing the time series features.

This section takes some time to process and compute (especially on the whole sample) and I already saved the results as a csv which I will just work from and load in the pre-computed time series features.

################# 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")

Note: Again, bad practice by me. I just called the df data Stats which consists of only the time series features. This still only refers to the train_val.csv data and not the test.csv data.

The training data looks like: (after computing the time series features). Now each asset has been collapsed from ~260 days down to 1 signal time series feature observation.

Recall the goal here was to classify synthetic time series vs real time series and not what the next days price is going to be. For each asset I have a signal observation and based on this I can train a classifying algorithm to distinguish between real vs synthetic time series.

How the training data looks:

## [1] 12000   109

The dimensions of the data as still 12,000 with 109 features (created from the tsfeatures package). That is we have 6,000 synthetic and 6,000 real financial time series (12,000 * ~260 = 3,120,000 but we applied tsfeatures to collapse the ~260 down to 1 single observation for each asset)

I collapsed this problem down from a time series prediction problem to a pure classification problem. I split the data between training and validation set next… I also split the data into X_train, Y_train… etc.

I split the df/Stats data set into a train set of 75% of the observations and an in-sample test data set of 25% of the observations.

######################################################################
################# 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)

How the training X (input variables) data looks:

How the training Y (predictor variable) data looks:

I set the data up for an XGBoost model:

I create a grid search in order search over a parameter space to locate the optimal parameters for the data set. It needs a little more work but it’s a pretty good starting point. I can just add code to the expand.grid function. That is, say I want to increase the depth of the tree I can add to max_depth = c(5, 8, 14) more parameters such as max_depth = c(5, 8, 14, 1, 2, 3, 4, 6, 7). Note Adding parameters to the grid search increases computational time exponentially. Every parameter you add a value to, the model has to search all possible combinations associated with that parameter. That is, adding an eta = c(0.1) and max_depth = c(5) would give me the optimal parameter for one iteration/loop through the training model, i.e. an eta = c(0.1) mapped onto a max_depth = c(5). Adding an additional value to the eta = c(0.1, 0.3) and max_depth = c(5) would map eta = 0.1 onto max_depth = 5 and eta = 0.3 on to max_depth = 5. If I add another value such that eta = c(0.1, 0.3, 0.4) then all 3 of these values will be mapped to max_depth = c(5). Adding values to the max_depth = c(5) parameter would add an extra layer of complexity to the grid search. This added into the fact that there are many parameters to optimize in an XGBoost model can drastically increase computational complexity. Thus, understanding the statistics behind the models in Machine Learning is important when trying to avoid getting stuck in a local minimum (which any greedy algorithm using gradient descent optimisation can do: greedy algorithm).

######################################################################
################# 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)
    }))

The output of the grid search can be set into a nice data frame using the following code. However I did not save this output to file and therefore cannot read it in. You can view the output on the original Jupyter Notebook In [49] here

# 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)

According to the results at the time the optimal parameters were:

• ntrees = 95,
• eta = 0.1,
• max_depth = 5,

With the other parameters left to default settings for simplicity.

Plug the optimal parameters into the model.

#################################################################################
################# 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

Now that I have the optimal parameters from the cross validation grid search I can train the final XGBoost model on the whole train_val.csv data set. (Whereas before the optimal parameters were obtained from different folds in the model. More info on k-fold cross validation here)

# 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...

What is nice about tree based models is that we can obtain importance scores from the model and find which variables contributed most to the gain in the model. The original paper explains more about the gain in Algorithm 1 and Algorithm 3 here.

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

That is, the XGBoost model found that the spike was the most important variable. The spike comes from the stl_features function of the tsfeatures package in R. It computes various measures of trend and seasonality based on Seasonal and Trend Decomposition (STL) and measures the spikiness of a time series based on the variance of the leave-one-out variances of component e_t.

The second variable is interesting also and comes from the compengine feature set from the CompEngine database. It groups variables as autocorrelation, prediction, stationarity, distribution and scaling.

The ARCH.LM comes from the arch_stat function of the tsfeatures package and is based on the Lagrange Multiplier for Autoregressive Conditional Heteroscedasticity (ARCH) Engle 1982.

These are just a few of the variables the XGBoost model found to be the most important. A full overview and more information of the variables used in the model can be found here.

# 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")