The aim of this tutorial is to show the use of TensorFlow with KERAS for classification and prediction in Time Series Analysis. The latter just implement a Long Short Term Memory (LSTM) model (an instance of a Recurrent Neural Network which avoids the vanishing gradient problem).

Introduction

The code below has the aim to quick introduce Deep Learning analysis with TensorFlow using the Keras back-end in R environment. Keras is a high-level neural networks API, developed with a focus on enabling fast experimentation and not for final products. Keras and in particular the keras R package allows to perform computations using also the GPU if the installation environment allows for it.

Installing KERAS and TensorFlow in Windows … otherwise it will be more simple

1. Install Anaconda: https://www.anaconda.com/download/
2. Install Rtools(34): https://cran.r-project.org/bin/windows/Rtools/

GPU-TensorFlow

To use the option GPU-TensorFlow, you need CUDA Toolkit that matches the version of your GCC compiler: https://developer.nvidia.com/cuda-toolkit

If you have Python (i.e. Anaconda) just

install.packages("keras")
library(keras)
install_keras()

and this will install the Google Tensorflow module in Python.

If you want it working on GPU and you have a suitable CUDA version, you can install it with tensorflow = "gpu" option

install_keras(tensorflow = "gpu")

Simple check

library(keras)
to_categorical(0:3)

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1

Background on Neural Networks

Example old faithful IRIS data

Consider the well-known IRIS data set

rm(list=ls())
data(iris)
plot(iris$Petal.Length,
 iris$Petal.Width, col = iris$Species)

head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

We want to build an iris specie classifier based on the observed four iris dimensions. This is the usual classification (prediction) problem so we have to consider a training sample and evaluate the classifier on a test sample.

Data in TensorFlow

Data are

1. Matrices “`matrix´´´ of doubles.
2. Categorical variables need to be codified in dummies: one hot encoding.

onehot.species = to_categorical(as.numeric(iris$Species) - 1)
iris = as.matrix(iris[, 1:4])
iris = cbind(iris, onehot.species)

Training and Test Data Sets

Define training and test

set.seed(17)
ind <- sample(2, nrow(iris),
replace = TRUE, prob = c(0.7, 0.3))
iris.training <- iris[ind == 1, 1:4]
iris.test <- iris[ind == 2, 1:4]
iris.trainingtarget <- iris[ind == 1, -seq(4)]
iris.testtarget <- iris[ind == 2, -(1:4)]

Model building

Initialize the model

model <- keras_model_sequential()

and suppose to use a very simple one

model %>%
 layer_dense(units = ncol(iris.trainingtarget), activation = 'softmax',
 input_shape = ncol(iris.training))
summary(model)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #
## ===========================================================================
## dense_1 (Dense)                  (None, 3)                     15
## ===========================================================================
## Total params: 15
## Trainable params: 15
## Non-trainable params: 0
## ___________________________________________________________________________

This is the model structure:

model$inputs

## [[1]]
## Tensor("dense_1_input:0", shape=(?, 4), dtype=float32)

model$outputs

## [[1]]
## Tensor("dense_1/Softmax:0", shape=(?, 3), dtype=float32)

and a graphical representation:

library(igraph)

##
## Attaching package: 'igraph'

## The following objects are masked from 'package:keras':
##
##     %<-%, normalize

## The following objects are masked from 'package:stats':
##
##     decompose, spectrum

## The following object is masked from 'package:base':
##
##     union

g = graph_from_literal(Sepal.Length:Sepal.Width:Petal.Length:Petal.Width---Species,simplify = TRUE)
layout <- layout_in_circle(g, order = order(degree(g)))
plot(g,layout = layout,vertex.color = c(2,2,2,2,3))

In the plot, blue colors stand for input and green ones for output.

Its analytic representation is the following one:

$$Species_j = act.func(\mathbf{w}_j,\mathbf{x} = (PW,PL,SW,SL)),$$

where the activation function is the softmax (the all life logistic!):

$$act.func(\mathbf{w}_j,\mathbf{x}) = \frac{e^{\mathbf{x}^T\mathbf{w}}}{\sum e^{\mathbf{x}^T\mathbf{w}}}$$

which estimates $Pr(Specie = j|\mathbf{x} = (PW,PL,SW,SL))$.

Model fitting: fit() and the optimizer

Estimation consists in finding the weights $\mathbf{w}$ that minimizes a loss function. For instance, if the response $Y$ were quantitative, then

$$w = \arg\min \sum_{i = 1}^m(y_i-wx_i)^2,$$

whose solution is given by the usual equations of derivatives $w$:

$$\frac{\partial \sum_{i = 1}^n(y_i-wx_i)^2}{\partial w} = 0,$$

Note however, that

$$\partial \sum (y_i-wx_i)^2 = \sum \partial (y_i-wx_i)^2,$$

(Is parallelizable in batches of samples (of length batch_size), that is

$$\sum \partial (y_i-wx_i)^2 = \sum{\partial\sum (y_i-wx_i)^2}$$

where $n_l$ is the batch_size.

Suppose in general a non-analytical loss function (the usual case in more complicated networks) $Q(w) = \sum_{i = 1}^m(y_i-wx_i)^2,$ and suppose that $\frac{\partial Q(w)}{\partial w} = 0$ is not available analytically. Then we would have to use “Newton-Raphson” optimizer family (or gradient optimizers) whose best known member in Deep Learning (DL) is the Stochastic Gradient Descent (SGD):

Starting form an initial weight $w^{(0)}$ at step $m$:

$$w^{(m)} = w^{(m-1)}-\eta\Delta Q_i(w),$$ where $\eta>0$ is the Learning Rate: the lower (bigger) $\eta$ is, the more (less) steps are needed to achieve the optimum with a greater (worse) precision.

It is stochastic in the sense that the index $i$ of the sample is random (avoids overfitting): $\Delta Q(w) : = \Delta Q_i(w)$. This also induces complications when (if) dealing with time series.

There are other variations to the SGD: Momentum, Averaging, AdaGrad, Adam, …

Using SGD with $\eta = 0.01$ we have to set:

sgd <- optimizer_sgd(lr = 0.01)

and then this is plugged in into the model and used afterwards in compilation. Once it is established, the loss function $Q$ (here we use the categorical_crossentropy because the response is a non-binary categorical variable):

model %>% compile(
 loss = 'categorical_crossentropy',
 optimizer = sgd,
 metrics = 'accuracy'
)

we have to train it in epochs (i.e. the $m$ steps above) using a portion of the training sample, validation_split, to verify eventual overfitting (i.e. the model is fitted and the loss evaluated in that random part of the sample which is finally not used for training):

history <- model %>% fit(
 x = iris.training,
 y = iris.trainingtarget,
 epochs = 100,
 batch_size = 5,
 validation_split = 0.2,
 verbose = 0
)

The result of the trained model is:

plot(history)

Validation on the test sample:

classes <- model %>% predict_classes(iris.test)
table(iris.testtarget%*%0:2, classes)

##    classes
##      0  1  2
##   0 20  0  0
##   1  0 12  0
##   2  0  7  8

with a validation score

(score <- model %>% evaluate(iris.test, iris.testtarget))

## $loss
## [1] 0.3292574
##
## $acc
## [1] 0.8510638

Another example: Classification of breast cancer

We have 10 variables (all factors) and a binary response: benign versus malign.

library(mlbench)
data(BreastCancer)
dim(BreastCancer)

## [1] 699  11

levels(BreastCancer$Class)

## [1] "benign"    "malignant"

head(BreastCancer)

##        Id Cl.thickness Cell.size Cell.shape Marg.adhesion Epith.c.size
## 1 1000025            5         1          1             1            2
## 2 1002945            5         4          4             5            7
## 3 1015425            3         1          1             1            2
## 4 1016277            6         8          8             1            3
## 5 1017023            4         1          1             3            2
## 6 1017122            8        10         10             8            7
##   Bare.nuclei Bl.cromatin Normal.nucleoli Mitoses     Class
## 1           1           3               1       1    benign
## 2          10           3               2       1    benign
## 3           2           3               1       1    benign
## 4           4           3               7       1    benign
## 5           1           3               1       1    benign
## 6          10           9               7       1 malignant

str(BreastCancer)

## 'data.frame':	699 obs. of  11 variables:
##  $ Id             : chr  "1000025" "1002945" "1015425" "1016277" ...
##  $ Cl.thickness   : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 5 5 3 6 4 8 1 2 2 4 ...
##  $ Cell.size      : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 1 4 1 8 1 10 1 1 1 2 ...
##  $ Cell.shape     : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 1 4 1 8 1 10 1 2 1 1 ...
##  $ Marg.adhesion  : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 1 5 1 1 3 8 1 1 1 1 ...
##  $ Epith.c.size   : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 2 7 2 3 2 7 2 2 2 2 ...
##  $ Bare.nuclei    : Factor w/ 10 levels "1","2","3","4",..: 1 10 2 4 1 10 10 1 1 1 ...
##  $ Bl.cromatin    : Factor w/ 10 levels "1","2","3","4",..: 3 3 3 3 3 9 3 3 1 2 ...
##  $ Normal.nucleoli: Factor w/ 10 levels "1","2","3","4",..: 1 2 1 7 1 7 1 1 1 1 ...
##  $ Mitoses        : Factor w/ 9 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 5 1 ...
##  $ Class          : Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...

Data in matrices

tt = BreastCancer[complete.cases(BreastCancer),2:11]
x = NULL
for(i in seq(9)) x = cbind(x,to_categorical(as.numeric(tt[,i])-1))
y = to_categorical(as.numeric(tt[,10])-1)
head(y)

##      [,1] [,2]
## [1,]    1    0
## [2,]    1    0
## [3,]    1    0
## [4,]    1    0
## [5,]    1    0
## [6,]    0    1

Set training and test

set.seed(17)
ind <- sample(2, nrow(x),
 replace = TRUE, prob = c(0.7, 0.3))

x.train = x[ind == 1,]
y.train = y[ind == 1,]
x.test = x[ind == 2,]
y.test = y[ind == 2,]

Let’s build the DL model with tree layers of neurons:

# Initialize a sequential model
model <- keras_model_sequential()

# Add layers to model
model %>%
 layer_dense(units = 8, activation = 'relu', input_shape = ncol(x.train)) %>%
 layer_dense(units = 5, activation = 'relu') %>%
 layer_dense(units = ncol(y.train), activation = 'softmax')

summary(model)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #
## ===========================================================================
## dense_2 (Dense)                  (None, 8)                     720
## ___________________________________________________________________________
## dense_3 (Dense)                  (None, 5)                     45
## ___________________________________________________________________________
## dense_4 (Dense)                  (None, 2)                     12
## ===========================================================================
## Total params: 777
## Trainable params: 777
## Non-trainable params: 0
## ___________________________________________________________________________

As activation function (being the response binary) we use a user defined relu ($f(x) = x^+$):

Let’s use the adam optimizer

# Compile the model
model %>% compile(
 loss = 'categorical_crossentropy',
 optimizer = "adam",
 metrics = 'accuracy'
)

Train the model

history <- model %>% fit(
 x = x.train,
 y = y.train,
 epochs = 50,
 batch_size = 50,
 validation_split = 0.2,
 verbose = 2
)
plot(history)

Validate it on the test set

classes <- model %>% predict_classes(x.test)
table(y.test%*%0:1, classes)

##    classes
##       0   1
##   0 121   6
##   1   5  70

also with a score

(score <- model %>% evaluate(x.test, y.test))

## $loss
## [1] 0.1734528
##
## $acc
## [1] 0.9455446

LSTM model

Here we apply the DL to time series analysis: it is not possible to draw train and test randomly and they must be random sequences of train and test of length batch_size.

Data

From Yahoo Finance let’s download the IBEX 35 time series on the last 15 years and consider the last 3000 days of trading:

library(BatchGetSymbols)
tickers <- c('%5EIBEX')
first.date <- Sys.Date() - 360*15
last.date <- Sys.Date()

YAHOO database query and the ACF of the considered IBEX 35 series is here:

myts <- BatchGetSymbols(tickers = tickers,
 first.date = first.date,
 last.date = last.date,
 cache.folder = file.path(tempdir(),
 'BGS_Cache') ) # cache in tempdir()

##
## Running BatchGetSymbols for:
##    tickers = %5EIBEX
##    Downloading data for benchmark ticker | Not Cached
## %5EIBEX | yahoo (1|1) | Not Cached - Feliz que nem lambari de sanga!

print(myts$df.control)

##    ticker   src download.status total.obs perc.benchmark.dates
## 1 %5EIBEX yahoo              OK      3787            0.9903304
##   threshold.decision
## 1               KEEP

y = myts$df.tickers$price.close
myts = data.frame(index = myts$df.tickers$ref.date, price = y, vol = myts$df.tickers$volume)
myts = myts[complete.cases(myts), ]
myts = myts[-seq(nrow(myts) - 3000), ]
myts$index = seq(nrow(myts))

library(plotly)

plot_ly(myts, x = ~index, y = ~price, type = "scatter", mode = "markers", color = ~vol)

acf(myts$price, lag.max = 3000)

Training and Testing samples

Data must be standardized

msd.price = c(mean(myts$price), sd(myts$price))
msd.vol = c(mean(myts$vol), sd(myts$vol))
myts$price = (myts$price - msd.price[1])/msd.price[2]
myts$vol = (myts$vol - msd.vol[1])/msd.vol[2]
summary(myts)

##      index            price               vol
##  Min.   :   1.0   Min.   :-2.20595   Min.   :-1.2713
##  1st Qu.: 750.8   1st Qu.:-0.73810   1st Qu.:-1.2689
##  Median :1500.5   Median :-0.06936   Median : 0.1166
##  Mean   :1500.5   Mean   : 0.00000   Mean   : 0.0000
##  3rd Qu.:2250.2   3rd Qu.: 0.36329   3rd Qu.: 0.6992
##  Max.   :3000.0   Max.   : 3.00692   Max.   : 4.7057

Let’s use the first 2000 days for training and the last 1000 for test. Remember that the ratio between the number of train samples and test samples must be an integer number as also the ratio between these two lengths with batch_size. This is why ²⁰⁰⁰⁄₁₀₀₀, ²⁰⁰⁰⁄₅₀ and ¹⁰⁰⁰⁄₅₀:

datalags = 10
train = myts[seq(2000 + datalags), ]
test = myts[2000 + datalags + seq(1000 + datalags), ]
batch.size = 50

Data for LSTM

Predictor $X$ is a 3D matrix: 1. first dimension is the length of the time series 2. second is the lag; 3. third is the number of variables used for prediction $X$ (at least 1 for the series at a given lag).

Response $Y$ is a 2D matrix: 1. first dimension is the length of the time series 2. second is the lag;

x.train = array(data = lag(cbind(train$price, train$vol), datalags)[-(1:datalags), ], dim = c(nrow(train) - datalags, datalags, 2))
y.train = array(data = train$price[-(1:datalags)], dim = c(nrow(train)-datalags, 1))

x.test = array(data = lag(cbind(test$vol, test$price), datalags)[-(1:datalags), ], dim = c(nrow(test) - datalags, datalags, 2))
y.test = array(data = test$price[-(1:datalags)], dim = c(nrow(test) - datalags, 1))

The LSTM model codified with Keras

model <- keras_model_sequential()

model %>%
 layer_lstm(units = 100,
 input_shape = c(datalags, 2),
 batch_size = batch.size,
 return_sequences = TRUE,
 stateful = TRUE) %>%
 layer_dropout(rate = 0.5) %>%
 layer_lstm(units = 50,
 return_sequences = FALSE,
 stateful = TRUE) %>%
 layer_dropout(rate = 0.5) %>%
 layer_dense(units = 1)

model %>%
 compile(loss = 'mae', optimizer = 'adam')

model

## Model
## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #
## ===========================================================================
## lstm_1 (LSTM)                    (50, 10, 100)                 41200
## ___________________________________________________________________________
## dropout_1 (Dropout)              (50, 10, 100)                 0
## ___________________________________________________________________________
## lstm_2 (LSTM)                    (50, 50)                      30200
## ___________________________________________________________________________
## dropout_2 (Dropout)              (50, 50)                      0
## ___________________________________________________________________________
## dense_5 (Dense)                  (50, 1)                       51
## ===========================================================================
## Total params: 71,451
## Trainable params: 71,451
## Non-trainable params: 0
## ___________________________________________________________________________

Let’s train in 2000 steps. Remember: for being the model stateful (stateful = TRUE), which means that the signal state (the latent part of the model) is trained on the batch of the time series, you need to manually reset the states (batches are supposed to be independent sequences (!) ):

for(i in 1:2000){
 model %>% fit(x = x.train,
 y = y.train,
 batch_size = batch.size,
 epochs = 1,
 verbose = 0,
 shuffle = FALSE)
 model %>% reset_states()
}

The prediction

pred_out <- model %>% predict(x.test, batch_size = batch.size) %>% .[,1]

plot_ly(myts, x = ~index, y = ~price, type = "scatter", mode = "markers", color = ~vol) %>%
 add_trace(y = c(rep(NA, 2000), pred_out), x = myts$index, name = "LSTM prediction", mode = "lines")

## Warning: line.color doesn't (yet) support data arrays

## Warning: line.color doesn't (yet) support data arrays

more on validation:

plot(x = y.test, y = pred_out)

Some notes on Deep Learning

A deep learning (DL) model is a neural network with many layers of neurons (Schmidhuber 2015), it is an algorithmic approach rather than probabilistic in its nature, see (Breiman and others 2001) for the merits of both approaches. Each neuron is a deterministic function such that a neuron of a neuron is a function of a function along with an associated weight $w$. Essentially for a response variable $Y_i$ for the unit $i$ and a predictor $X_i$ we have to estimate $Y_i = w_1f_1(w_2f_2(…(w_kf_k(X_i))))$, and the larger $k$ is, the “deeper” is the network. With many stacked layers of neurons all connected (a.k.a. dense layers) it is possible to capture high non-linearities and all interactions among variables. The approach to model estimation underpinned by a DL model is that of composition function against that od additive function underpinned by the usual regression techniques including the most modern one (i.e. $Y_i = w_1f_1(X_i)+w_2f_2(X_i)+…+w_kf_k(X_i)$). A thorough review of DL can be found at (Schmidhuber 2015).

Likely the DL model can be also interpreted as a maximum a posteriori estimation of $Pr(Y|X,Data)$ (Polson, Sokolov, and others 2017) for Gaussian process priors. Despite this and because of its complexity it cannot be evaluated the whole distribution $Pr(Y|X,Data)$, but only its mode.

Estimating a DL consists in just estimating the vectors $w_1,\ldots,w_k$. The estimation requires to evaluate a multidimensional gradient which is not possible to be evaluated jointly for all observations, because of its dimensionality and complexity. Recalling that the derivative of a composite function is defined as the product of the derivative of inner functions (i.e. the chain rule $(f\circ g)’ = (f’\circ g)\cdot g’$) which is implemented for purposes of computational feasibility as a tensor product. Such tensor product is evaluated for batches of observations and it is implemented in the open source software known as Google Tensor Flow (Abadi et al. 2015).

Fundamentals of LSTM can be found here (Sherstinsky 2018) (it needs some translation to the statistical formalism).

References

• Abadi, Martín, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, et al. 2015. “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems.” URL: https://www.tensorflow.org/.

• Breiman, Leo, and others. 2001. “Statistical Modeling: The Two Cultures (with Comments and a Rejoinder by the Author).” Statistical Science 16 (3). Institute of Mathematical Statistics: 199-231.

• Polson, Nicholas G, Vadim Sokolov, and others. 2017. “Deep Learning: A Bayesian Perspective.” Bayesian Analysis 12 (4). International Society for Bayesian Analysis: 1275-1304.

• Schmidhuber, Jürgen. 2015. “Deep Learning in Neural Networks: An Overview.” Neural Networks 61. Elsevier: 85-117.

• Sherstinsky, Alex. 2018. “Fundamentals of Recurrent Neural Network (Rnn) and Long Short-Term Memory (Lstm) Network,” August. URL: http://arxiv.org/abs/http://arxiv.org/abs/1808.03314v4.