In [1]:
dpath = 'data/basic/'
#path_to_dataset = dpath + 'household_power_consumption.txt'
In [2]:
%mkdir -p dpath
!wget -P $dpath https://raw.githubusercontent.com/jaungiers/LSTM-Neural-Network-for-Time-Series-Prediction/master/sinwave.csv
!wget -P $dpath https://raw.githubusercontent.com/jaungiers/LSTM-Neural-Network-for-Time-Series-Prediction/master/sp500.csv
In [3]:
import matplotlib.pyplot as plt
%matplotlib inline
import time
import warnings
import numpy as np
from numpy import newaxis
from keras.layers.core import Dense, Activation, Dropout
from keras.layers.recurrent import LSTM
from keras.models import Sequential
warnings.filterwarnings("ignore")
In [83]:
import pandas as pd
def load_ts(filename):
df = pd.read_csv(filename, header=None)
data = df[0].tolist()
return data
#filename = 'sinwave.csv'
filename = 'sp500.csv'
print('> Loading data...: ', dpath+filename)
#X_train, y_train, X_test, y_test = load_data(dpath+'sp500.csv', seq_len, True)
ts = load_ts(dpath + filename)
ts[:10]
Out[83]:
In [84]:
plt.plot(ts)
Out[84]:
In [85]:
seq_len = 50
def load_data(ts, seq_len, normalise_window):
sequence_length = seq_len + 1
result = []
# create gliding window
for index in range(len(ts) - sequence_length):
result.append(ts[index: index + sequence_length])
if normalise_window:
result = normalise_windows(result)
result = np.array(result)
print("Data shape: ", result.shape)
print(result[:4, :])
row = round(0.9 * result.shape[0])
train = result[:row, :]
test = result[row:, :]
print("Test shape: ", test.shape)
#np.random.shuffle(train)
x_train = train[:, :-1]
y_train = train[:, -1]
x_test = test[:, :-1]
print("xtest shape: ", x_test.shape)
y_test = test[:, -1]
x_train = np.reshape(x_train, (x_train.shape[0], x_train.shape[1], 1))
x_test = np.reshape(x_test, (x_test.shape[0], x_test.shape[1], 1))
return [x_train, y_train, x_test, y_test]
def normalise_windows(window_data):
normalised_data = []
for window in window_data:
normalised_window = [((float(p) / float(window[0])) - 1) for p in window]
normalised_data.append(normalised_window)
return normalised_data
print('> Loading data...')
#X_train, y_train, X_test, y_test = load_data(ts, seq_len, False)
X_train, y_train, X_test, y_test = load_data(ts, seq_len, True)
X_train.shape, y_train.shape, X_test.shape, y_test.shape
X_train[0,:seq_len, 0]
X_train[1,:seq_len, 0]
y_train[:5]
Out[85]:
Out[85]:
Out[85]:
Out[85]:
In [86]:
#plt.plot(X_train[0,:,0])
#plt.plot(y_train)
fig = plt.figure(facecolor='white')
ax = fig.add_subplot(211)
ax.plot(X_train[0, :, 0], label='X')
ax.legend()
ax = fig.add_subplot(212)
ax.plot(y_train, label='y')
ax.legend()
Out[86]:
Out[86]:
Out[86]:
Out[86]:
In [66]:
def build_model(layers):
model = Sequential()
model.add(LSTM(
input_dim=layers[0],
output_dim=layers[1],
return_sequences=True))
model.add(Dropout(0.2))
model.add(LSTM(
layers[2],
return_sequences=False))
model.add(Dropout(0.2))
model.add(Dense(
output_dim=layers[3]))
model.add(Activation("linear"))
start = time.time()
model.compile(loss="mse", optimizer="rmsprop")
print("Compilation Time : ", time.time() - start)
return model
print('> Data Loaded. Compiling...')
model = build_model([1, 50, 100, 1])
model.summary()
In [87]:
start_time = time.time()
epochs = 50
model.fit(
X_train,
y_train,
batch_size=512,
nb_epoch=epochs,
validation_split=0.05)
print("Training time: ", time.time() - start_time)
Out[87]:
In [ ]:
def plot_results(predicted_data, true_data, figsize=(12,6)):
fig = plt.figure(facecolor='white', figsize=figsize)
ax = fig.add_subplot(111)
ax.plot(true_data, label='True Data')
plt.plot(predicted_data, label='Prediction')
plt.legend()
plt.show()
If you’re observant you’ll have noticed in our load_data() function above we split the data in to train/test sets as is standard practice for machine learning problems. However what we need to watch out for here is what we actually want to achieve in the prediction of the time series.
If we were to use the test set as it is, we would be running each window full of the true data to predict the next time step. This is fine if we are only looking to predict one time step ahead, however if we’re looking to predict more than one time step ahead, maybe looking to predict any emergent trends or functions (e.g. the sin function in this case) using the full test set would mean we would be predicting the next time step but then disregarding that prediction when it comes to subsequent time steps and using only the true data for each time step. s You can see below the graph of using this approach to predict only one time step ahead at each step in time:
In [73]:
def predict_point_by_point(model, data):
#Predict each timestep given the last sequence of true data, in effect only predicting 1 step ahead each time
predicted = model.predict(data)
predicted = np.reshape(predicted, (predicted.size,))
return predicted
start_time = time.time()
predicted = predict_point_by_point(model, X_test)
predicted[:5]
print("Prediction time: ", time.time() - start_time)
plot_results(predicted, y_test)
Out[73]:
If however we want to do real magic and predict many time steps ahead we only use the first window from the testing data as an initiation window. At each time step we then pop the oldest entry out of the rear of the window and append the prediction for the next time step to the front of the window, in essence shifting the window along so it slowly builds itself with predictions, until the window is full of only predicted values (in our case, as our window is of size 50 this would occur after 50 time steps). We then keep this up indefinitely, predicting the next time step on the predictions of the previous future time steps, to hopefully see an emerging trend.
In [71]:
def predict_sequence_full(model, data, window_size):
#Shift the window by 1 new prediction each time, re-run predictions on new window
curr_frame = data[0]
predicted = []
# loop over entire testdata
for i in range(len(data)):
predicted.append(model.predict(curr_frame[newaxis,:,:])[0,0]) #get element from shape(1,1)
curr_frame = curr_frame[1:] #move window
curr_frame = np.insert(curr_frame, [window_size-1], predicted[-1], axis=0) #fill frame with prediction
return predicted
start_time = time.time()
predicted = predict_sequence_full(model, X_test, seq_len)
predicted[:5]
print("Prediction time: ", time.time() - start_time)
plot_results(predicted, y_test)
Out[71]:
Overlaid with the true data we can see that with just 1 epoch and a reasonably small training set of data the LSTM has already done a pretty damn good job of predicting the sin function. You can see that as we predict more and more into the future the error margin increases as errors in the prior predictions are amplified more and more when they are used for future predictions. As such we see that the LSTM hasn’t got the frequency quite right and it drifts the more we try to predict it. However as the sin function is a very easy oscillating function with zero noise it can predict it to a good degree.
We predicted a several hundred time steps of a sin wave on an accurate point-by-point basis. So we can now just do the same on a stock market time series and make a shit load of money right?
Well, no.
A stock time series is unfortunately not a function that can be mapped. It can best described more as a random walk, which makes the whole prediction thing considerably harder. But what about the LSTM identifying any underlying hidden trends? Well, let’s take a look.
Here is a CSV file where I have taken the adjusted daily closing price of the S&P 500 equity index from January 2000 – August 2016. I’ve stripped out everything to make it in the exact same format as our sin wave data and we will now run it through the same model we used on the sin wave with the same train/test split.
There is one slight change we need to make to our data however, because a sin wave is already a nicely normalized repeating pattern it works well running the raw data points through the network. However running the adjusted returns of a stock index through a network would make the optimization process shit itself and not converge to any sort of optimums for such large numbers.
So to combat this we will take each n-sized window of training/testing data and normalize each one to reflect percentage changes from the start of that window (so the data at point i=0 will always be 0). We’ll use the following equations to normalise and subsequently de-normalise at the end of the prediction process to get a real world number out of the prediction:
n = normalised list [window] of price changes p = raw list [window] of adjusted daily return prices
In [ ]:
start_time = time.time()
predicted = predict_point_by_point(model, X_test)
predicted[:5]
print("Prediction time: ", time.time() - start_time)
In [100]:
plot_results(predicted, y_test, figsize=(20,10))
Running the data on a single point-by-point prediction as mentioned above gives something that matches the returns pretty closely. But this is deceptive! Why? Well if you look more closely, the prediction line is made up of singular prediction points that have had the whole prior true history window behind them. Because of that, the network doesn’t need to know much about the time series itself other than that each next point most likely won’t be too far from the last point. So even if it gets the prediction for the point wrong, the next prediction will then factor in the true history and disregard the incorrect prediction, yet again allowing for an error to be made.
We can’t see what is happening in the brain of the LSTM, but I would make a strong case that for this prediction of what is essentially a random walk (and as a matter of point, I have made a completely random walk of data that mimics the look of a stock index, and the exact same thing holds true there as well!) is “predicting” the next point with essentially a Gaussian distribution, allowing the essentially random prediction to not stray too wildly from the true data.
So what would we look at if we wanted to see whether there truly was some underlying pattern discernable in just the price movements? Well we would do the same as for the sin wave problem and let the network predict a sequence of points rather than just the next one.
Doing that we can now see that unlike the sin wave which carried on as a sin wave sequence that was almost identical to the true data, our stock data predictions converge very quickly into some sort of equilibrium.
In [ ]:
start_time = time.time()
predicted = predict_sequence_full(model, X_test, seq_len)
predicted[:5]
print("Prediction time: ", time.time() - start_time)
In [101]:
plot_results(predicted, y_test, figsize=(20,10))
In [102]:
def plot_results_multiple(predicted_data, true_data, prediction_len, figsize=(12,6)):
fig = plt.figure(facecolor='white', figsize=figsize)
ax = fig.add_subplot(111)
ax.plot(true_data, label='True Data')
#Pad the list of predictions to shift it in the graph to it's correct start
for i, data in enumerate(predicted_data):
padding = [None for p in range(i * prediction_len)]
plt.plot(padding + data, label='Prediction'+str(i))
plt.legend()
plt.show()
In [ ]:
def predict_sequences_multiple(model, data, window_size, prediction_len):
#Predict sequence of 50 steps before shifting prediction run forward by 50 steps
prediction_seqs = []
for i in range(len(data)//prediction_len):
curr_frame = data[i*prediction_len]
predicted = []
for j in range(prediction_len):
predicted.append(model.predict(curr_frame[newaxis,:,:])[0,0])
curr_frame = curr_frame[1:]
curr_frame = np.insert(curr_frame, [window_size-1], predicted[-1], axis=0)
prediction_seqs.append(predicted)
return prediction_seqs
start_time = time.time()
predictions = predict_sequences_multiple(model, X_test, seq_len, 50)
#predicted = predict_sequence_full(model, X_test, seq_len)
#predicted = predict_point_by_point(model, X_test)
print("Prediction time: ", time.time() - start_time)
In [103]:
plot_results_multiple(predictions, y_test, 50, figsize=(20,10))
In fact when we take a look at the graph above of the same run but with the epochs increased to 400 (which should make the model mode accurate) we see that actually it now just tries is predict an upwards momentum for almost every time period!
In [ ]: