<jupyter_start><jupyter_text>RNN for RegressionNow that you have seen how to use of RNN for classification, we will consider the problem of time-series regression.We will look at real data corresponding to the price of a cryptocurrency and will try to predict it.**Note**: we will be training a lot of networks and the output might clog your jupyter notebook. Remember that you can click in the left margin of an output in order to "hide" it.The next cell does the usual setup.<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.model_selection import train_test_split, GridSearchCV

from sklearn.metrics import accuracy_score, mean_absolute_error,mean_squared_error
    
from sklearn.pipeline import Pipeline
from sklearn.base import TransformerMixin, BaseEstimator

import keras
from keras import regularizers
from keras.layers import LSTM
from keras.layers import Input, Dense, LSTM, SimpleRNN, TimeDistributed
from keras.models import Model

# this is the same as copy pasting a bunch of "def" here.
%run utils/helpers<jupyter_output><empty_output><jupyter_text>Load the dataLoad the data `./data/ethereum_dataset` and have a look.<jupyter_code>eth = pd.read_csv('./data/ethereum_dataset.csv',)
print(eth.shape)
eth.head()<jupyter_output><empty_output><jupyter_text>Data Preprocessing NaNsHow many features contain NaNs and where do the NaNs come from?<jupyter_code>eth.apply(pd.isnull).sum()<jupyter_output><empty_output><jupyter_text>So the `eth_ens_register` has quite a few NaNs as compared to the number of rows in the dataset (~600 vs ~800). Should you just drop this feature?To investigate this further, we will first set the `Date(UTC)` as index then investigate this feature in more depth.  Date/TimeTake the `Date(UTC)` feature, make it into a datetime object thanks to `pd.to_datetime` and use it as an index.<jupyter_code>eth["Date(UTC)"] = pd.to_datetime(eth["Date(UTC)"])
eth = eth.set_index("Date(UTC)")<jupyter_output><empty_output><jupyter_text>Back to NaNsWe saw that `eth_ens_register` (the domain name service for an object on the enthereum blockchain) is `NaN`.As usual we have to deal with missing values and see if we should just drop that feature or maybe impute a sensible value. Start by checking the values for `eth_ens_register`<jupyter_code>eth["eth_ens_register"].unique()<jupyter_output><empty_output><jupyter_text>It turns out there are entries, which are exactly 0, have a look at the corresponding slice of the data, can you spot anything?<jupyter_code>eth[eth["eth_ens_register"]==0]<jupyter_output><empty_output><jupyter_text>All of these records are on or after the 3rd of May 2017. Most probably the ENS did not exist before that, let's confirm by checking the dates corresponding to the `NaN` values.<jupyter_code>eth[eth["eth_ens_register"].isnull()].tail(10)<jupyter_output><empty_output><jupyter_text>It looks like our hypothesis is confirmed. Since we will try to predict future values of ethereum, when the ENS exist, we can safely set them to 0.<jupyter_code>eth.loc[eth["eth_ens_register"].isnull(), "eth_ens_register"] = 0

# you can check it worked as expected
eth[eth["eth_ens_register"]==0].head()<jupyter_output><empty_output><jupyter_text>Exploratory Data AnalysisHave a look at some of the features, in particular have a look at the following features and discuss* `eth_etherprice`* `eth_supply`* `eth_tx`* `eth_difficulty`<jupyter_code>plt.figure(figsize=(8, 6))
eth["eth_etherprice"].plot()
plt.title("ETHERPRICE")

plt.figure(figsize=(8, 6))
eth["eth_supply"].plot()
plt.title("SUPPLY")

plt.figure(figsize=(8, 6))
eth["eth_tx"].plot()
plt.title("TX")

plt.figure(figsize=(8, 6))
eth["eth_difficulty"].plot()
plt.title("DIFFICULTY")<jupyter_output><empty_output><jupyter_text>According to the "Bloom Baster", on October 16, 2017, Ethereum was "renewed" which means that the reward for mining a block fell from 5 `eth` to 3 but the block time was also reduced. Because of this, `difficulty` dropped. If you're not into this cryptocurrency mumbo jumbo, note the discontinuity in "difficulty" and think about what that means for any predictive method using this data.Show that period specifically a drop in `difficulty`.<jupyter_code>eth['2017-10':'2017-11']["eth_difficulty"].plot(marker="o")
plt.ylabel("Difficulty")<jupyter_output><empty_output><jupyter_text>You can clearly see the drop.Apply a standard scaling to the data<jupyter_code>pipeline = Pipeline([
    ('scaling', StandardScaler()),
])

X = pipeline.fit_transform(eth.values)
print(X.shape)<jupyter_output><empty_output><jupyter_text>Build the modelWe want to be able to predict the price tomorrow (prediction horizon = 1 day) given everything but the price time series. Let's start by preprocessing all the values: scale the data. Now the $y$-s become tomorrow's price values, but we also have to drop our last data point.We'll give you that one, make sure it makes sense (especially the indices). Note that we'll exclude the price feature from now on to make it "more interesting".<jupyter_code>id_price = 1 # index of the appropriate column
y = np.expand_dims(X[1:, id_price], -1)

#initially we are going to exclude the price feature from the data set entirely
X_ = X[0:-1, np.arange(X.shape[1]) != id_price]

print(X_.shape)
print(y.shape)<jupyter_output><empty_output><jupyter_text>One of the main limitations of neural networks is they typically need a ton of data for training. 831 transaction is by no means enough. A (sub-optimal) way around this is to decrease the amount of testing data. An obvious caveat is that we could be testing on a very particular test set, not seen in the training data, which could lead to a disappointingly poor generalisation error.The cell below defines a training-testing split for arrays corresponding to time series, it should all look familiar.<jupyter_code>def train_test_split_time_series_regres(X, y, test_ratio=0.15):
    total_samples = X.shape[0]
    train_idx = int(total_samples * (1-test_ratio))
    XTrain = X[:train_idx]
    yTrain = y[:train_idx]
    XTest = X[train_idx:]
    yTest = y[train_idx:]
    return XTrain, yTrain, XTest, yTest
XTrain, yTrain, XTest, yTest = train_test_split_time_series_regres(X_, y)
print(XTrain.shape)
print(yTrain.shape)
print(XTest.shape)
print(yTest.shape)

BATCH_SIZE = 30
XTrain_batch = reshape_to_batches(XTrain, BATCH_SIZE)
yTrain_batch = reshape_to_batches(yTrain, BATCH_SIZE)
print(XTrain_batch.shape)
print(yTrain_batch.shape)

BATCH_SIZE = 30
XTest_batch = reshape_to_batches(XTest, BATCH_SIZE)
yTest_batch = reshape_to_batches(yTest, BATCH_SIZE)
print(XTest_batch.shape)
print(yTest_batch.shape)<jupyter_output><empty_output><jupyter_text>Create the model<jupyter_code>#import all dependencies
from keras.layers import Input, Dense, LSTM, TimeDistributed
from keras.models import Model<jupyter_output><empty_output><jupyter_text>This time we are going to make one little change. Because we would to be able to predict the price of ethereum given only yesterday's information, we will allow our network to accept batches of **any size**. This is accomplished by supplying `None` as a shape parameter.<jupyter_code>inputs = Input(shape=(None, 16))<jupyter_output><empty_output><jupyter_text>Now add an LSTM layer with 64 neurons and no regularisation.Note that* the output layer must have a `linear` activation function (it's a regression, not a classification that we're doing)* the loss is not the cross entropy loss but the `mean_squared_error` lossfit the model with 100 epochs (it will be very quick)<jupyter_code>lstm = LSTM(64, activation='tanh', recurrent_activation='hard_sigmoid', use_bias=True, 
            kernel_initializer='glorot_uniform', recurrent_initializer='orthogonal', 
            bias_initializer='zeros', unit_forget_bias=True, kernel_regularizer=None, 
            recurrent_regularizer=None, bias_regularizer=None, activity_regularizer=None, 
            kernel_constraint=None, recurrent_constraint=None, bias_constraint=None, dropout=0.0, 
            recurrent_dropout=0.0, implementation=1, return_sequences=True, return_state=False, 
            go_backwards=False, stateful=False, unroll=False)(inputs)

# the main difference with classification is that the activation of the last layer is linear
predictions = TimeDistributed(Dense(1, activation='linear'))(lstm)

lstm_model = Model(inputs=inputs, outputs=predictions)

# the second difference is the loss function (now MSE)
lstm_model.compile(optimizer='rmsprop',
                   loss='mean_squared_error',
                   metrics=['accuracy'])
lstm_model.fit(XTrain_batch, yTrain_batch, epochs=100)<jupyter_output><empty_output><jupyter_text>Generate predictions and check the mean squared error!<jupyter_code>y_pred_lstm = lstm_model.predict(XTest_batch)

print(yTest.shape)
print(y_pred_lstm.shape)

print("MAE: {0:.2f}".format(mean_absolute_error(yTest, _3d_to_2d(y_pred_lstm)[:125])))
print("MSE: {0:.2f}".format(mean_squared_error(yTest, _3d_to_2d(y_pred_lstm)[:125])))<jupyter_output><empty_output><jupyter_text>These values don't really tell you much at this point. Remember that you have normalised the values with 0-mean and variance 1 so that the small numbers should not be particularly impressive.We can get a better feeling of the quality of our predictions via some visualisations.Plot the predicted values vs the true values and discuss.<jupyter_code>plt.figure(figsize=(8, 6))
plt.scatter(yTest, _3d_to_2d(y_pred_lstm)[:125])
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=2)
plt.xlabel('Measured', fontsize=12)
plt.ylabel('Predicted', fontsize=12)<jupyter_output><empty_output><jupyter_text>Right... so that's not that great, it looks like our model is generally predicting higher values (bias).On a positive note, it looks like it is capturing the trend.  Adjusting the architectureSince training is very fast here, you can play around with the number of neurons and number of layers to try to get better results. One way to go is to increase the number of neurons (say to 256); feel free to tweak.<jupyter_code>inputs = Input(shape=(None, 16))
lstm = LSTM(256, activation='tanh', recurrent_activation='hard_sigmoid', use_bias=True, 
            kernel_initializer='glorot_uniform', recurrent_initializer='orthogonal', 
            bias_initializer='zeros', unit_forget_bias=True, kernel_regularizer=None, 
            recurrent_regularizer=None, bias_regularizer=None, activity_regularizer=None, 
            kernel_constraint=None, recurrent_constraint=None, bias_constraint=None, dropout=0.0, 
            recurrent_dropout=0.0, implementation=1, return_sequences=True, return_state=False, 
            go_backwards=False, stateful=False, unroll=False)(inputs)

predictions = TimeDistributed(Dense(1, activation='linear'))(lstm)
lstm_model256 = Model(inputs=inputs, outputs=predictions)
lstm_model256.compile(optimizer='rmsprop',
                      loss='mean_squared_error',
                      metrics=['accuracy'])
lstm_model256.fit(XTrain_batch, yTrain_batch, epochs=100)<jupyter_output><empty_output><jupyter_text>Evaluating your modelCopy paste the cell to evaluate your model adjusting for the new model.<jupyter_code>y_pred_lstm256 = lstm_model256.predict(XTest_batch)

print("MAE LSTM 64:  {0:.2f}".format(mean_absolute_error(yTest, _3d_to_2d(y_pred_lstm)[:125])))
print("MSE LSTM 64:  {0:.2f}".format(mean_squared_error(yTest, _3d_to_2d(y_pred_lstm)[:125])))
print("MAE LSTM 256: {0:.2f}".format(mean_absolute_error(yTest, _3d_to_2d(y_pred_lstm256)[:125])))
print("MSE LSTM 256: {0:.2f}".format(mean_squared_error(yTest, _3d_to_2d(y_pred_lstm256)[:125])))

plt.figure(figsize=(8, 6))
plt.plot(yTest, _3d_to_2d(y_pred_lstm)[:125], color="r", 
         ls="none", marker="o", markersize=4, alpha=0.2, label="lstm64")
plt.scatter(yTest, _3d_to_2d(y_pred_lstm256)[:125], label="lstm256")
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=2)
plt.xlabel('Measured', fontsize=12)
plt.ylabel('Predicted', fontsize=12)
plt.legend(fontsize=12)<jupyter_output><empty_output><jupyter_text>Ok so that's better but not hugely so.Try something more complicated, for example:* two layers* some dropout* some regularisation* more epochs<jupyter_code>inputs = Input(shape=(None, 16))

lstm = LSTM(64, activation='tanh', recurrent_activation='hard_sigmoid', use_bias=True, 
            kernel_initializer='glorot_uniform', recurrent_initializer='orthogonal', 
            bias_initializer='zeros', unit_forget_bias=True, 
            # -- reg 
            kernel_regularizer=regularizers.l2(0.001), 
            recurrent_regularizer=regularizers.l2(0.001), 
            bias_regularizer=None, activity_regularizer=None, 
            kernel_constraint=None, recurrent_constraint=None, bias_constraint=None, 
            # -- dropout
            dropout=0.3, 
            recurrent_dropout=0.04, 
            implementation=1, return_sequences=True, return_state=False, 
            go_backwards=False, stateful=False, unroll=False)(inputs)

lstm1 = LSTM(64, activation='tanh', recurrent_activation='hard_sigmoid', use_bias=True, 
            kernel_initializer='glorot_uniform', recurrent_initializer='orthogonal', 
            bias_initializer='zeros', unit_forget_bias=True, 
            # -- reg
            kernel_regularizer=regularizers.l2(0.001), 
            recurrent_regularizer=regularizers.l2(0.001), 
            bias_regularizer=None, activity_regularizer=None, kernel_constraint=None, 
            recurrent_constraint=None, bias_constraint=None, 
            # -- dropout
            dropout=0.3, 
            recurrent_dropout=0.04, 
            implementation=1, return_sequences=True, return_state=False, 
            go_backwards=False, stateful=False, unroll=False)(lstm)

dense =  TimeDistributed(Dense(64, activation='sigmoid'))(lstm1)
predictions = TimeDistributed(Dense(1, activation='linear'))(dense)
model_lstm256x256r = Model(inputs=inputs, outputs=predictions)

model_lstm256x256r.compile(optimizer='rmsprop',
                           loss='mean_squared_error',
                           metrics=['accuracy'])
model_lstm256x256r.fit(XTrain_batch, yTrain_batch, epochs=500)<jupyter_output><empty_output><jupyter_text>Evaluating the model<jupyter_code>y_pred_lstm256x256r = model_lstm256x256r.predict(XTest_batch)

print("MAE LSTM 64:       {0:.2f}".format(mean_absolute_error(yTest, _3d_to_2d(y_pred_lstm)[:125])))
print("MSE LSTM 64:       {0:.2f}".format(mean_squared_error(yTest, _3d_to_2d(y_pred_lstm)[:125])))
print("MAE LSTM 256:      {0:.2f}".format(mean_absolute_error(yTest, _3d_to_2d(y_pred_lstm256)[:125])))
print("MSE LSTM 256:      {0:.2f}".format(mean_squared_error(yTest, _3d_to_2d(y_pred_lstm256)[:125])))
print("MAE LSTM 256x256r: {0:.2f}".format(mean_absolute_error(yTest, _3d_to_2d(y_pred_lstm256x256r)[:125])))
print("MSE LSTM 256x256r: {0:.2f}".format(mean_squared_error(yTest, _3d_to_2d(y_pred_lstm256x256r)[:125])))

plt.figure(figsize=(8, 6))
plt.plot(yTest, _3d_to_2d(y_pred_lstm)[:125], color="r", 
         ls="none", marker="o", markersize=4, alpha=0.2, label="lstm64")
plt.plot(yTest, _3d_to_2d(y_pred_lstm256)[:125], color="c", 
         ls="none", marker="o", markersize=4, alpha=0.4, label="lstm256")
plt.scatter(yTest, _3d_to_2d(y_pred_lstm256x256r)[:125], label="lstm256x256r")
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=2)
plt.xlabel('Measured', fontsize=12)
plt.ylabel('Predicted', fontsize=12)
plt.legend(fontsize=12)<jupyter_output><empty_output><jupyter_text>Let's try some other methods to validate our results. Look at the actual values and compare with the target ones. Remember that we scaled all of the features. To rescale just the price, we are going to pad the rest of the values with zeros.<jupyter_code>plt.plot(_3d_to_2d(y_pred_lstm256x256r)[:125])
plt.plot(yTest)<jupyter_output><empty_output><jupyter_text>The network has learned to capture the general trend, and it has learned that at a certain condition there is a spike but it overestimates the spikes.  GeneratorNow that we have an (more or less) acceptable machinery we can apply the same approach to not only predict the price but also all other features (i.e. multi-dimensional output). To do that instead of a shifted price, we are going to use as target values all features shifted by one (one day horizon).<jupyter_code>y_shift = X[1:, :] # the "future"
X_shift = X[0:-1, :] # the "past"
print(X_shift.shape)
print(y_shift.shape)<jupyter_output><empty_output><jupyter_text>We are going to put aside the testing here (bad!) and focus on trying to have a good model for the problem above and interpret the results. By using all the data we should have the best fit though we will not be able to give any guarantees about generalisability.<jupyter_code># reshape to expected output
BATCH_SIZE = 30
X_batch = reshape_to_batches(X_shift, BATCH_SIZE)
y_batch = reshape_to_batches(y_shift, BATCH_SIZE)
print(X_batch.shape)
print(y_batch.shape)<jupyter_output><empty_output><jupyter_text>Building the netLet's use as a first layer a 512 LSTM without regularisation<jupyter_code>inputs = Input(shape=(None, 17)) # This returns a tensor
lstm = LSTM(512, activation='tanh', recurrent_activation='hard_sigmoid', use_bias=True, 
            kernel_initializer='glorot_uniform', recurrent_initializer='orthogonal', 
            bias_initializer='zeros', unit_forget_bias=True, kernel_regularizer=None, 
            recurrent_regularizer=None, bias_regularizer=None, activity_regularizer=None, 
            kernel_constraint=None, recurrent_constraint=None, bias_constraint=None, dropout=0.0, 
            recurrent_dropout=0.0, implementation=1, return_sequences=True, return_state=False, 
            go_backwards=False, stateful=False, unroll=False)(inputs)<jupyter_output><empty_output><jupyter_text>The only thing which we need to change is the output dimensions of the network from 1 to 17.<jupyter_code>predictions = TimeDistributed(Dense(17, activation='linear'))(lstm)
model_lstm512 = Model(inputs=inputs, outputs=predictions)

model_lstm512.compile(optimizer='rmsprop',
              loss='mean_squared_error',
              metrics=['accuracy'])<jupyter_output><empty_output><jupyter_text>Let's turn the magic on, train for 500 epochs, more may be better but may crash your computer.<jupyter_code>model_lstm512.fit(X_batch, y_batch, epochs=500)<jupyter_output><empty_output><jupyter_text>Trying to predict the futureIn order to make our network tell the future, we need to make it generate the values for tomorrow and then iteratively feed this back into the network.The last data point is from November 7, let's try to predict until the end of January (very dubious from a statistical perspective but we can try).<jupyter_code>#remember that we left out the last X sample, we can start from there
days = 3 * 30

X_last = X[-1, :]
X_batch = np.swapaxes(np.expand_dims(np.expand_dims(X_last, -1), -1), 0, 2)

for day in range(days):
    print("Day #{} - {} data points".format(day, X_batch.shape[1]))
    y_pred = model_lstm512.predict(X_batch)
    # we are only going to use the most recent prediction
    # otherwise the prediction power could quickly deteriorate
    y_pred = np.swapaxes(np.expand_dims(np.expand_dims(y_pred[0, X_batch.shape[1]-1, :], -1), -1), 0, 2)
    X_batch = np.concatenate([X_batch, y_pred], axis=1)
df = pd.DataFrame(
        pipeline.inverse_transform(
            X_batch.reshape(X_batch.shape[1], X_batch.shape[2])))
df.iloc[:, 0] = pd.to_datetime(df.iloc[:, 0], unit='s')

df.head()<jupyter_output><empty_output><jupyter_text>Very interesting, it seems like our model has not even learned that the day increments at every time step. Do you thing this is a bug or a feature? Let's try to give it the entire history to check.<jupyter_code>#remember that we left out the last X sample, we can start from there
days = 3*30

X_last = X[:, :]
X_batch = np.swapaxes(np.swapaxes(np.expand_dims(X_last, -1), 0, 2), 1, 2)
print(X_batch.shape)

for i in range(days):
    print("Day #{} - {} data points".format(i,X_batch.shape[1]))
    y_pred = model_lstm512.predict(X_batch)
    y_pred = np.swapaxes(np.expand_dims(np.expand_dims(y_pred[0,X_batch.shape[1]-1,:], -1), -1), 0, 2)
    X_batch = np.concatenate([X_batch,y_pred], axis=1)

df = pd.DataFrame(
        pipeline.inverse_transform(
                X_batch.reshape(X_batch.shape[1], X_batch.shape[2])))
df.iloc[:,0] = pd.to_datetime(df.iloc[:, 0], unit='s')

df.tail()<jupyter_output><empty_output><jupyter_text>Yes, it is confirmed, it is trying to learn the days as a function of the other feature. This is an illustration that the NN is trying to learn the **joint probability** of the features. In this case, the network seems to have learned that a price of between `$`285-300, along with other features are usually seen around October, so it predicts October as the date. Notice that it has learned that the time should also increase, and in fact the October 04th dates increase with around 30-40 minutesLet's have a look of how our model envisions the price change compared to the general trend.<jupyter_code>df.iloc[:,1].plot()
df.iloc[-90:,1].plot()
#plot just the prediction time series
df.iloc[-90:,1].plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10687201572886434, 'fraction_numerical': 0.025372156165153075, 'mean_word_length': 5.497262773722627, 'pattern_counts': {'":': 0, '<': 125, '<?xml version=': 0, '>': 125, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '242033', 'license': 'MIT', 'n_tokens_mistral': 7452, 'n_tokens_neox': 6893, 'n_words': 1948, 'path': 'ads8-neuralnets_content_students/d2-RNN-for-regression_solution.ipynb', 'repo_name': 'JamesOwers/camsparknn', 'size': '21362'}