Keras LSTM Neutal Networks for Univariate Time-Series in R

If any issues, questions or suggestions feel free to reach me out via e-mail wieczynskipawel@gmail.com or Linkedin. You can also visit my Github.

Based on this tutorial I wrote functions which facilitate data preprocessing, fitting LSTM neural network for univariate time-series and interpretation of the results. Although this script is written in R, it runs in the background Keras and Tensorflow environment. Be cautious with the version compatibility. In my case it works with Python 3.10.4 and Tensorflow 2.8

if(!require('pacman')) install.packages('pacman')
pacman::p_load(quantmod, keras, tidyverse, timetk, lubridate, Metrics)

Data normalization

LSTM needs data to be normalized, e.g. centered by the mean and scaled by the standard deviation. Function get_scaling_factors() stores scaling factors as we will use them for data preprocessing as well as to “denormalize” data after making predictions. Function normalize_data() works in both directions: to normalize (reverse=FALSE) and denormalize data (reverse=TRUE).

get_scaling_factors <- function(data){
  out <- c(mean = mean(data), sd = sd(data))
  return(out)
}

normalize_data <- function(data, scaling_factors, reverse = FALSE) {
  
  if (reverse) temp <- (data * scaling_factors[2]) + scaling_factors[1]
  else temp <- (data - scaling_factors[1]) / scaling_factors[2]
  
  out <- temp %>% as.matrix()
  return(out)
}

Data preprocessing

Keras LSTM as input takes data in specific 3D format [samples, timesteps, features]. Function which I named kerasize_data() takes the following parameters:
- data - vector of normalized data
- x - if x=TRUE (default) then regressors are being “kerasized”. If x=FALSE targets are being “kerasized”. In fact, for univariate time series targets are just shifted regressors.
- lag - number of previous observations to be used for time-serie prediction (default=12)
- pred - number of future observiation to be predicted (default=12).

Additionaly kerasize_pred_input() does the same for last n observations which will be used for prediction moving forward.

Those functions can be further generalized to be compatible with multivariate time-series.

kerasize_data <- function(data, x = TRUE, lag = 12, pred = 12) {
  
  if (x) {
    
    temp <- sapply(
      1:(length(data) - lag - pred + 1)
      ,function(x) data[x:(x + lag - 1), 1]
    ) %>% t()
    
    out <- array(
      temp %>% unlist() %>% as.numeric()
      ,dim = c(nrow(temp), lag, 1)
    )
    
  }  else {
    
    temp <- sapply(
      (1 + lag):(length(data) - pred + 1)
      ,function(x) data[x:(x + lag - 1), 1]
    ) %>% t()
    
    out <- array(
      temp %>% unlist() %>% as.numeric()
      ,dim = c(nrow(temp), pred, 1)
    )
    
  }
  
  return(out)
  
}

kerasize_pred_input <- function(data, lag = 12, pred = 12){
  temp <- data[(length(data) - pred + 1):length(data)]
  temp <- normalize_data(temp, get_scaling_factors(data))
  out <- array(temp, c(1, lag, 1))
  return(out)
}

Fitting LSTM neural network

Function lstm_build_model() takes the following parameters:
- x - 3D array of kerasized regressor data, i.e. output of the kerasize_data() function with parameter x=TRUE
- y - 3D array of kerasized target data, i.e. output of the kerasize_data() function with parameter x=FALSE
- units - size of the layer (default=50)
- batch - the number of samples to be propagated thru the network (default=1)
- epochs - how many times algorithm goes thru the dataset (default = 20)
- rate - fraction of the input units to drop (default=0.5)
- seed - sets random seed for R, Numpy and Tensorflow (default=2137).

If you are familiar with the LSTM architecture you can play with function’s definition or just play with some parameters, i.e. units, batch, epochs, rate. There is also extensive documentation of the Keras library available. Depending on your computing power or compatibility of your GPU with Tensorflow it takes more or less time to build the model.

lstm_build_model <- function(x, y, units = 50, batch = 1, epochs = 20, rate = 0.5, seed = 2137){
  
  lag = dim(x)[2]
  
  lstm_model <- keras_model_sequential()

  lstm_model %>%
    layer_lstm(units = units
               ,batch_input_shape = c(batch, lag, 1)
               ,return_sequences = TRUE
               ,stateful = TRUE) %>%
    layer_dropout(rate = rate) %>%
    layer_lstm(units = units
               ,return_sequences = TRUE
               ,stateful = TRUE) %>%
    layer_dropout(rate = rate) %>%
    time_distributed(layer_dense(units = 1))

  lstm_model %>%
    compile(loss = 'mae'
            ,optimizer = 'adam'
            ,metrics = 'accuracy')

  tensorflow::set_random_seed(seed)
  lstm_model %>% fit(
    x = x
    ,y = y
    ,batch_size = batch
    ,epochs = epochs
    ,verbose = 0
    ,shuffle = FALSE)
  
  out <- list(
    model = lstm_model
    ,x = x
    ,batch = batch
    ,lag = lag
    ,pred = dim(y)[2]
  )
  return(out)

}

LSTM prediction

Function lstm_forecast() takes the following parameters:
- x_test - output of the kerasize_pred_input() function
- model - output of the lstm_build_model() function
- scaling_factors - output of the get_scaling_factors() function.

lstm_forecast <- function(x_test, model, scaling_factors){
  
  batch <- model$batch
  
  temp <- model$model %>%
    predict(x_test, batch_size = batch) %>% 
    .[, , 1] %>%
    normalize_data(scaling_factors = scaling_factors, reverse = TRUE)
  
  out <- list(
    forecast = temp
    ,scaling_factors = scaling_factors
  )
  
  return(out)
  
}

Transform results to the forecast object

Function forecast_transform() takes the following parameters:
- data - initial data in xts format
- model - output of the lstm_build_model() function
- forecast - output of the lstm_forecast() function

It will facilitate further analysis and plotting of the results.

forecast_transform <- function(data, model, forecast){
  
  lag <- model$lag
  
  freq <- periodicity(data)$scale
  frequency <- case_when(
    freq == 'weekly' ~ 52
    ,freq == 'monthly' ~ 12
   , freq == 'quarterly' ~ 4
    ,freq == 'yearly' ~ 1
  )
  
  dates <- index(data)
  date_start <- c(year(min(dates)), month(min(dates)), day(min(dates)))
  date_end <- c(year(max(dates)), month(max(dates)), day(max(dates)))
  
  date_start_shifted <- case_when(
    freq == 'weekly' ~ min(dates) %m+% weeks(lag)
    ,freq == 'monthly' ~ min(dates) %m+% months(lag)
    ,freq == 'quarterly' ~ min(dates) %m+% months(3*lag)
    ,freq == 'yearly' ~ min(dates) %m+% years(lag)
    )
  date_start_shifted <- c(year(date_start_shifted), month(date_start_shifted), day(date_start_shifted))
  
  date_end_shifted <- case_when(
    freq == 'weekly' ~ max(dates) %m+% weeks(1)
    ,freq == 'monthly' ~ max(dates) %m+% months(1)
    ,freq == 'quarterly' ~ max(dates) %m+% months(3)
    ,freq == 'yearly' ~ max(dates) %m+% years(1)
    )
  date_end_shifted <- c(year(date_end_shifted), month(date_end_shifted), day(date_end_shifted))

  fitted <- predict(model$model, model$x, batch_size = model$batch) %>% .[, , 1]
  
  if (dim(fitted)[2] > 1) {
    fitted <- c(fitted[, 1], fitted[dim(fitted)[1], 2:dim(fitted)[2]])
  } else {
    fitted <- fitted[, 1]
  }
  
  fitted <- normalize_data(fitted, forecast$scaling_factors, reverse = TRUE)
  fitted <- ts(fitted, start = date_start_shifted, deltat = 1/frequency)
  
  lstm_forecast <- ts(forecast$forecast, start = date_end_shifted, deltat = 1/frequency)
  
  data_trimmed <- data[(model$lag+1):nrow(data)] %>% ts(start = date_start_shifted
                                                        ,deltat = 1/frequency)
  
  
  out <- list(
    model = NULL
    ,method = 'LSTM'
    ,mean = lstm_forecast
    ,x = data_trimmed
    ,fitted = fitted
    ,residuals = as.numeric(data_trimmed) - as.numeric(fitted)
  )
  
  class(out) <- 'forecast'
  
  return(out)
}

Example of use

We download data from FRED database using quantmod package. Below univariate time-serie consists of monthly data on Retail Sales: Beer, Wine, and Liquor Stores from 1992 to 2022 (362 observations).

data <- getSymbols('MRTSSM4453USS', src = 'FRED', auto.assign = FALSE)
dim(data); head(data); tail(data)

[1] 362   1
           MRTSSM4453USS
1992-01-01          1713
1992-02-01          1763
1992-03-01          1753
1992-04-01          1784
1992-05-01          1783
1992-06-01          1782
           MRTSSM4453USS
2021-09-01          5839
2021-10-01          5846
2021-11-01          5934
2021-12-01          5789
2022-01-01          5929
2022-02-01          6088

Here we run functions defined and explained above. We use the default lag=12 and pred=12.

scaling_factors <- get_scaling_factors(data$MRTSSM4453USS)
data_normalized <- normalize_data(data$MRTSSM4453USS, scaling_factors)

x_data <- kerasize_data(data_normalized,  x = TRUE)
y_data <- kerasize_data(data_normalized, x = FALSE)
x_test <- kerasize_pred_input(data_normalized)

model <- lstm_build_model(x_data, y_data)
prediction <- lstm_forecast(x_test, model, scaling_factors)

final_results <- forecast_transform(data, model, prediction)

What returned function forecast_transform()?

class(final_results)

[1] "forecast"

Prediction 12 months ahead.

final_results$mean

          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
2022                   6061.903 6077.350 6070.080 6071.238 6070.613 6070.217 6070.649 6070.571 6070.761 6070.717
2023 6070.605 6071.059                                                                                          

Fitted values from the LSTM model. It consists only 350 observations due to fact that first 12 observations were used to estimate the 13th one.

length(final_results$fitted); head(final_results$fitted, 20)

[1] 350
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
1993 6065.197 2291.137 2455.853 2458.048 2455.469 2452.407 2452.266 2458.632 2458.522 2459.598 2457.225 2456.776
1994 2456.234 2455.839 2455.624 2455.425 2457.487 2459.290 2462.305 2460.422                                    

Empirical values. First 12 values were omitted here to obtain the same vector lenth as theoretical values.

length(final_results$x); head(final_results$x, 20)

[1] 350
      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1993 1840 1836 1814 1792 1791 1796 1799 1786 1767 1768 1780 1773
1994 1796 1809 1825 1832 1858 1844 1862 1840                    

Residuals, i.e. differences beetwen fitted and empirical values.

head(final_results$residuals, 20)

 [1] -4225.1973  -455.1367  -641.8534  -666.0480  -664.4688  -656.4067  -653.2656  -672.6322  -691.5219  -691.5978
[11]  -677.2250  -683.7761  -660.2342  -646.8386  -630.6242  -623.4253  -599.4868  -615.2900  -600.3050  -620.4218

On the above objects we can calculate different error metrics, especially using Metrics package.

cat(paste0(
  'Mean Absolute Error: ', mae(final_results$x, final_results$fitted) %>% round (2), '\n'
  ,'Mean Absolute Percent Error: ', mape(final_results$x, final_results$fitted) %>% round (4), '\n'
  ,'Root Mean Squared Error: ', rmse(final_results$x, final_results$fitted) %>% round (2), '\n'
  ,'R squared: ', (cor(final_results$x, final_results$fitted)[1])^2 %>% round (4)
)) 

Mean Absolute Error: 524.59
Mean Absolute Percent Error: 0.171
Root Mean Squared Error: 672.08
R squared: 0.9037

We can use also function forecast::autoplot on the output of the forecast_transform() function.

forecast::autoplot(final_results)

Example with quarterly granularity

Above set of functions should work well with weekly, quarterly and anually univariate time-series as well.

As a time-serie with quarterly frequency let’s take Median Sales Price of Houses Sold for the United States. It contains 236 observations from 1963 to 2021.

data_quarterly <- getSymbols('MSPUS', src = 'FRED', auto.assign = FALSE)
dim(data_quarterly); head(data_quarterly); tail(data_quarterly)

[1] 236   1
           MSPUS
1963-01-01 17800
1963-04-01 18000
1963-07-01 17900
1963-10-01 18500
1964-01-01 18500
1964-04-01 18900
            MSPUS
2020-07-01 337500
2020-10-01 358700
2021-01-01 369800
2021-04-01 382600
2021-07-01 411200
2021-10-01 408100

Here we use lag=4 and pred=4.

scaling_factors <- get_scaling_factors(data_quarterly$MSPUS)
data_normalized <- normalize_data(data_quarterly$MSPUS, scaling_factors)

x_data <- kerasize_data(data_normalized, x = TRUE, lag = 4, pred = 4)
y_data <- kerasize_data(data_normalized, x = FALSE, lag = 4, pred = 4)
x_test <- kerasize_pred_input(data_normalized, lag = 4, pred = 4)

model <- lstm_build_model(x_data, y_data)
prediction <- lstm_forecast(x_test, model, scaling_factors)

final_results <- forecast_transform(data_quarterly, model, prediction)

final_results$mean

         Qtr1     Qtr2     Qtr3     Qtr4
2022 345461.4 346901.4 348643.6 349260.7

head(final_results$fitted, 20)

          Qtr1      Qtr2      Qtr3      Qtr4
1964 310151.41  88043.11  61899.75  69332.46
1965  73579.86  73847.98  73653.61  73673.55
1966  73789.27  73805.36  73899.39  74061.21
1967  74173.41  74251.08  74263.75  74384.09
1968  74497.44  74599.57  74633.95  74818.42

head(final_results$x, 20)

      Qtr1  Qtr2  Qtr3  Qtr4
1964 18500 18900 18900 19400
1965 20200 19800 20200 20300
1966 21000 22100 21500 21400
1967 22300 23300 22500 22900
1968 23900 24900 24800 25600

cat(paste0(
  'Mean Absolute Error: ', mae(final_results$x, final_results$fitted) %>% round (2), '\n'
  ,'Mean Absolute Percent Error: ', mape(final_results$x, final_results$fitted) %>% round (4), '\n'
  ,'Root Mean Squared Error: ', rmse(final_results$x, final_results$fitted) %>% round (2), '\n'
  ,'R squared: ', (cor(final_results$x, final_results$fitted)[1])^2 %>% round (4)
)) 

Mean Absolute Error: 34914.53
Mean Absolute Percent Error: 0.6166
Root Mean Squared Error: 46684.88
R squared: 0.8996

forecast::autoplot(final_results)