Prediction interval around LSTM time series forecast
Is there a method to calculate the prediction interval (probability distribution) around a time series forecast from an LSTM (or other recurrent) neural network?
Say, for example, I am predicting 10 samples into the future (t+1 to t+10), based on the last 10 observed samples (t-9 to t), I would expect the prediction at t+1 to be more accurate than the prediction at t+10. Typically, one might draw error bars around the prediction to show the interval. With an ARIMA model (under the assumption of normally distributed errors), I can calculate a prediction interval (e.g. 650
95%) around each predicted value. Can I calculate the same, (or something that relates to the prediction interval) from an LSTM model?
I'm been working with LSTMs in Keras/Python, following lots of examples from machinelearningmastery.com, from which my example code (below) is based on. I'm considering reframing the problem as classification into discrete bins, as that produces a confidence per class, but that seems a poor solution.
There are a couple of similar topics (such as the below), but nothing seems to directly address the issue of prediction intervals from LSTM (or indeed other) neural networks:
https://stats.stackexchange.com/questions/25055/how-to-calculate-the-confidence-interval-for-time-series-prediction
https://datascience.stackexchange.com/questions/12721/time-series-prediction-using-arima-vs-lstm
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from math import sin
from matplotlib import pyplot
import numpy as np
# Build an LSTM network and train
def fit_lstm(X, y, batch_size, nb_epoch, neurons):
X = X.reshape(X.shape[0], 1, X.shape[1]) # add in another dimension to the X data
y = y.reshape(y.shape[0], y.shape[1]) # but don't add it to the y, as Dense has to be 1d?
model = Sequential()
model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
model.add(Dense(y.shape[1]))
model.compile(loss='mean_squared_error', optimizer='adam')
for i in range(nb_epoch):
model.fit(X, y, epochs=1, batch_size=batch_size, verbose=1, shuffle=False)
model.reset_states()
return model
# Configuration
n = 5000 # total size of dataset
SLIDING_WINDOW_LENGTH = 30
SLIDING_WINDOW_STEP_SIZE = 1
batch_size = 10
test_size = 0.1 # fraction of dataset to hold back for testing
nb_epochs = 100 # for training
neurons = 8 # LSTM layer complexity
# create dataset
#raw_values = [sin(i/2) for i in range(n)] # simple sine wave
raw_values = [sin(i/2)+sin(i/6)+sin(i/36)+np.random.uniform(-1,1) for i in range(n)] # double sine with noise
#raw_values = [(i%4) for i in range(n)] # saw tooth
all_data = np.array(raw_values).reshape(-1,1) # make into array, add anothe dimension for sci-kit compatibility
# data is segmented using a sliding window mechanism
all_data_windowed = [np.transpose(all_data[idx:idx+SLIDING_WINDOW_LENGTH]) for idx in np.arange(0,len(all_data)-SLIDING_WINDOW_LENGTH, SLIDING_WINDOW_STEP_SIZE)]
all_data_windowed = np.concatenate(all_data_windowed, axis=0).astype(np. 856