Monte Carlo Dropout for Predicting Prices with Deep Learning and Tensorflow Applied to the IBEX35

    	
            
# !pip install tensorflow-gpu==1.15
 
# # # You might need to run the below:
# # !pip install --upgrade pip
# # !pip install --upgrade tensorflow-gpu

        
        Copy
    

    	
            
import eikon as ek
import pandas as pd
import sys
import plotly.graph_objects as go
import plotly as pyo
import numpy as np
import warnings
import tensorflow as tf
from tensorflow.keras import layers
import tensorflow.experimental.numpy as tnp

        
        Copy
    

    	
            
print(sys.version) # 3.8.8 (default, Feb 24 2021, 15:54:32) [MSC v.1928 64 bit (AMD64)]
for i, j in zip(["eikon", "pandas", "numpy", "plotly", "tensorflow"],
                [ek, pd, np, pyo, tf]):
    print(f"The Python library '{i}' imported in this script is version '{j.__version__}'")

        
        Copy
    

    	
            
def get_seq_model(neurons = 16, dr = 0.5, act = "tanh", cell = layers.LSTM):
    inputs = layers.Input(shape = (sequence, X_train.shape[1]))
    x = cell(neurons, activation = act, return_sequences = True)(inputs)
    x = layers.BatchNormalization()(x)
    x = layers.Dropout(dr)(x, training = True)
    x = cell(neurons, activation = act)(x)
    x = layers.Flatten()(x)
    x = layers.BatchNormalization()(x)
    x = layers.Dropout(dr)(x, training = True)
    x = layers.Dense(neurons // 2, activation = act)(x)
    outputs = layers.Dense(1, activation = "elu")(x)
    model = tf.keras.models.Model(inputs = inputs, outputs = outputs)
    return model
 
def get_conv_model(neurons = 16, k = 5, dr = 0.5, act = "tanh"):
    inputs = layers.Input(shape = (sequence, X_train.shape[1]))
    x = layers.Conv1D(neurons, k, activation = act)(inputs)
    x = layers.BatchNormalization()(x)
    x = layers.Dropout(dr)(x, training = True)
    x = layers.Conv1D(neurons, k, activation = act)(x)
    x = layers.Flatten()(x)
    x = layers.BatchNormalization()(x)
    x = layers.Dropout(dr)(x, training = True)
    x = layers.Dense(neurons // 2, activation = act)(x)
    outputs = layers.Dense(1, activation = "elu")(x)
    model = tf.keras.models.Model(inputs = inputs, outputs = outputs)
    return model

        
        Copy
    

    	
            
# The key is placed in a text file so that it may be used in this code without showing it itself:
eikon_key = open("eikon.txt","r")
ek.set_app_key(str(eikon_key.read()))
# It is best to close the files we opened in order to make sure that we don't stop any other services/programs from accessing them if they need to:
eikon_key.close()

        
        Copy
    

    	
            dates = pd.date_range('1995', '2020', freq = 'YS')
        
        Copy
    

    	
            
df_full = pd.DataFrame()
chain = "0#.IBEX" # Constituents of the Ibex 35
fields = ["TR.GICSSector"]

        
        Copy
    

    	
            
for date in dates:
    year = int(str(date.to_period('Y')))
    day = f"{year}-01-01"
    params = {'SDate': day, "EDate": day, "Frq": "FY"} 
 
    # We obtain the constituents of the Ibex 35 at the start 
    # of each year and store that date as a column
 
    df = ek.get_data(f'{chain}({day})', fields, params)[0]
    df["Date"] = f'{date.to_period("D")}'
 
    # Then we drop the missing values and store them in
    # the previously created dataframe
 
    df = df.replace('', np.nan).dropna()
    df_full = df_full.append(df)

        
        Copy
    

    	
            df_full.T
        
        Copy
    

    	
            
# We store the names of the companies and the year they were listed in the index.
stocks = df_full.groupby("Instrument")["Date"].min().index
starts = df_full.groupby("Instrument")["Date"].min()
df_ts = pd.DataFrame() # We will fill this dataframe with all the data

        
        Copy
    

    	
            
# This lambda function takes a year and returns that date for the next step years
step = 10
to_date = lambda x: (f"{x}-01-01", f"{x+step}-12-31")
 
# e.g.:
[to_date(x) for x in range(2000, 2020, step)]

        
        Copy
    

    	
            
cols = ['OPEN', 'CLOSE', 'HIGH', 'LOW', 'VOLUME']
 
for stock, start in zip(stocks, starts):
    year = int(start[:4])
    date_range = [to_date(x) for x in range(year, 2020, step)]
    for st, end in date_range:
        try:
            df = ek.get_timeseries(stock, cols, start_date=st, end_date = end, interval = "daily")
        except ek.EikonError:
            try:
                df = ek.get_timeseries(stock, cols, start_date=st, end_date = None, interval = "daily")
            except ek.EikonError:
                continue
        df["stock"] = stock
        df[cols[1:]] = df[cols[1:]].pct_change().shift(1)
        df["target"] = df["OPEN"].pct_change().shift(-1)
        df["OPEN"] = df["OPEN"].pct_change()
 
        # We create a multi-index appending the identifier of the stock to the 'dates' index.
        df = df.set_index("stock", append = True)
        df_ts = df_ts.append(df.dropna())
 
# Finally we sort the data by its index
df_ts = df_ts.sort_index()

        
        Copy
    

    	
            df_ts
        
        Copy
    

    	
            
for col in df_ts.columns:
    df_ts[col] = df_ts[col].astype(np.float32)
 
cols = df_ts.drop(columns = "target").columns
 
X_train = df_ts.loc[:"2011", slice(None), :].drop(columns = "target")
df_ts[cols] = (df_ts[cols] - X_train.min()) / (X_train.max() - X_train.min())

        
        Copy
    

    	
            df_ts
        
        Copy
    

    	
            
X_train = df_ts.loc[:"2011", slice(None), :].drop(columns = "target")
y_train = df_ts.loc[:"2011", slice(None), :].pop("target")
X_val = df_ts.loc["2012":"2015", slice(None), :].drop(columns = "target")
y_val = df_ts.loc["2012":"2015", slice(None), :].pop("target")
X_test = df_ts.loc["2016":, slice(None), :].drop(columns = "target")
y_test = df_ts.loc["2016":, slice(None), :].pop("target")
X_full = df_ts.drop(columns = "target")

        
        Copy
    

    	
            
sequence = 22
batch = 256
stocks = df_ts.index.get_level_values(1).unique()

        
        Copy
    

    	
            
loc = lambda x, s: x.loc[slice(None), s, :]
to_ds = lambda x, y: tf.keras.preprocessing.timeseries_dataset_from_array(
    x, y, sequence_length=sequence, batch_size=batch)
 
ds_train = to_ds(loc(X_train, stocks[0]), loc(y_train, stocks[0]))
ds_val = to_ds(loc(X_val, stocks[0]), loc(y_val, stocks[0]))
ds_test = to_ds(loc(X_test, stocks[0]), loc(y_test, stocks[0]))
ds_err = {'train': [], 'val': [], 'test': []}
 
for stock in stocks[1:]:
    try:
        ds_train = ds_train.concatenate(to_ds(loc(X_train, stock), loc(y_train, stock)))
    except KeyError as err:
        ds_err['train'].append(err)
    try:
        ds_val = ds_val.concatenate(to_ds(loc(X_val, stock), loc(y_val, stock)))
    except KeyError as err:
        ds_err['val'].append(err)
    try:
        ds_test = ds_test.concatenate(to_ds(loc(X_test, stock), loc(y_test, stock)))
    except KeyError as err:
        ds_err['test'].append(err)
        continue

        
        Copy
    

    	
            
def loss(y_true, y_pred):
    dif =  - tf.reduce_mean(y_true * y_pred - 0.0025 * tf.math.abs(y_pred)) /  (tf.math.reduce_std(y_true * y_pred) + 1e-10) 
    return dif

        
        Copy
    

    	
            
early_stop_val = tf.keras.callbacks.EarlyStopping(
    monitor='val_loss', 
    patience=15, 
    verbose=0, 
    restore_best_weights=True,
    min_delta = 0.01
)
 
early_stop_train = tf.keras.callbacks.EarlyStopping(
    monitor='loss', 
    patience=7, 
    verbose=0, 
    restore_best_weights=True,
    min_delta = 0.01
)

        
        Copy
    

    	
            
def lr_sc(ep, lr):
    if ep <= 10 : return 1.1*lr
    return lr * tf.math.exp(-(ep-10) / 100)
 
mod_lr = tf.keras.callbacks.LearningRateScheduler(lr_sc)

        
        Copy
    

    	
            
history = []
models = [get_seq_model(cell = layers.SimpleRNN), get_seq_model(cell = layers.GRU), 
          get_seq_model(), get_conv_model()]
 
tf.random.set_seed(0)
 
for model in models:
    tf.keras.backend.clear_session()
    optimizer=tf.optimizers.RMSprop(3e-4)
    model.compile(optimizer = optimizer, loss = loss)
    history.append(
        model.fit(
            ds_train,
            epochs = 100, 
            validation_data=ds_val, 
            verbose = 0, # This value can be increased to visualize the progress.
            callbacks = [early_stop_val, mod_lr, early_stop_train],
            use_multiprocessing = True,
            )
        )

        
        Copy
    

    	
            
def make_predictions(samples = 100):
    stocks = X_full.index.get_level_values(1).unique()
    for i, stock in enumerate(stocks):
        for m, model in enumerate(models):
            data = to_ds(X_full.loc[(slice(None), stock), :], None)
            predictions = [model.predict(data, use_multiprocessing = True) for _ in range(samples)]
            mean = np.append([np.nan]*(sequence - 1), tf.reduce_mean(predictions, axis = 0).numpy())
            std = np.append([np.nan]*(sequence - 1), tf.math.reduce_std(predictions, axis = 0).numpy())
            first = np.append([np.nan]*(sequence - 1), predictions[0])
            df_ts.loc[(slice(None), stock), f"pred_mean_{m}"] = mean
            df_ts.loc[(slice(None), stock), f"pred_{m}"] = first
            df_ts.loc[(slice(None), stock), f"pred_std_{m}"] = std
        if (i % 5 == 0): print(f"{100*i/len(stocks):.2f} %")

        
        Copy
    

    	
            
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    make_predictions()

        
        Copy
    

    	
            
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    for m, model in enumerate(models):
        df_ts[f"pred_weighted_{m}"] = df_ts.groupby("Date")[f"pred_{m}"].apply(lambda x: x / (x.abs().sum()))

        
        Copy
    

    	
            
cps = 0.0001
names = ["SimpleRNN", "GRU", "LSTM", "Convolutional"]
 
mask = lambda x, y : (df_ts.dropna()[f"pred_std_{x}"] * y < df_ts.dropna()[f"pred_mean_{x}"].abs()) 
r_mean = lambda x, y: (df_ts.dropna()[mask(x, y)].groupby("Date")[f"pred_mean_{x}"].apply(lambda z: z / (np.sum(np.abs(z)))) * df_ts.dropna()[mask(x, y)]["target"]).groupby("Date").sum()
costs_mean = lambda x, y: (df_ts.dropna()[mask(x, y)].groupby("stock")[f"pred_mean_{x}"].apply(lambda z: z / (np.sum(np.abs(z)))).diff(1).abs() * cps).groupby("Date").sum()
r_mean_total = lambda x, y: r_mean(x, y) - costs_mean(x, y)
r_single = lambda x: (df_ts.dropna()[f"pred_weighted_{x}"] * df_ts.dropna()["target"]).groupby("Date").sum()
costs_single = lambda x: (df_ts.dropna().groupby("stock")[f"pred_weighted_{x}"].diff(1).abs() * cps).groupby("Date").sum()
r_single_total = lambda x: r_single(x) - costs_single(x)

        
        Copy
    

    	
            
# # Setting up sigma range up to a max of 10:
 
max_sigma_l = []
for x in range(len(names)):
    max_sigma = 0
    while True in (mask(x, max_sigma + 1)).to_list() and max_sigma < 10:
        max_sigma += 1
    max_sigma_l.append(max_sigma)
max_sigma = min(max_sigma_l)
 
if max_sigma >= 10: sigmas = [0, 2, 3, 5, 10]
elif max_sigma > 5: sigmas = [0, 2, 3, 5, max_sigma]
elif max_sigma <= 5:
    sigmas = [0, 2, 3, 5]
    print(f"Have too few sigmas, only up to {max_sigma}")

        
        Copy
    

    	
            rets_mean = [r_mean_total(x, sigma) + 1 for x, _ in enumerate(names) for sigma in sigmas]
        
        Copy
    

    	
            
rets_single = [r_single_total(x) + 1 for x, _ in enumerate(names)]
 
names_mean = ["Average prediction " + name + f" {y} sigmas filter" for name in names for y in sigmas]
names_single = ["Single prediction " + name for name in names]

        
        Copy
    

    	
            
df_ibex = ek.get_timeseries(chain[2:], "OPEN", start_date="2016-01-01", end_date = "2020-12-31", interval = "daily")
df_ibex["OPEN"] = df_ibex["OPEN"].pct_change()
 
rets_b = df_ibex.dropna()["OPEN"] + 1

        
        Copy
    

    	
            
def plot_returns(title = "", date = "2016", single = True, mean = True, std = max_sigma):
    
    rets_mean_ = [rets_mean[x] for x in range(sigmas.index(std), len(rets_mean), len(sigmas))]
    names_mean_ = [names_mean[x] for x in range(sigmas.index(std), len(names_mean), len(sigmas))]
    
    fig = go.Figure()
 
    if mean:
        [fig.add_trace(go.Scatter(
            x = r.loc[date:"2020"].index,
            y = (r.loc[date:"2020"]).cumprod(),
            name = n,
            mode = "lines"))
         for r, n in zip(rets_mean_, names_mean_)]
    if single:
        [fig.add_trace(go.Scatter(
            x = r.loc[date:"2020"].index,
            y = (r.loc[date:"2020"]).cumprod(),
            name = n))
         for r, n in zip(rets_single, names_single)]
    fig.add_trace(go.Scatter(
        x=rets_b.loc[date:"2020"].index,
        y=(rets_b.loc[date:"2020"]).cumprod(),
        name = "BuyHold"))
 
    fig.update_annotations(font = dict(size = 20))
    
    fig.update_layout(
        title = f"{title}",
        xaxis_title = "Date",
        yaxis_title = "Returns (%)",
        yaxis_tickformat = '%',
        template = "plotly_dark",
        font = dict(size = 16),
        
    )    
    fig.show()

        
        Copy
    

    	
            
plot_returns(title = "First prediction of the models", mean = False)
for std in sigmas:
    plot_returns(title = f"Average predictions using {std} sigmas filter", single = False, std = std)

        
        Copy
    

    	
            
def sharpe_ratio(returns):
    return np.sqrt(returns.groupby(returns.index.year).count().mean()) * (returns - 1).mean() / (returns - 1).std()
 
def information_rate(returns, benchmark):
    return (returns - benchmark).mean() / (returns - benchmark).std()
 
def sortino_ratio(returns):
    return np.sqrt(returns.groupby(returns.index.year).count().mean()) * (returns - 1).mean() / (returns[(returns -1) < 0] - 1).std()

        
        Copy
    

    	
            
date = "2016"
df_rets = pd.DataFrame(columns = ["Sharpe ratio", "IR", "Sortino ratio"], dtype = np.float64)
for i, ret in enumerate(rets_mean):
    ret = ret.loc[date:]
    rets_b = rets_b.loc[date:]
    df_rets.loc[names_mean[i]] = sharpe_ratio(ret), information_rate(ret, rets_b), sortino_ratio(ret)
    
for i, ret in enumerate(rets_single):
    ret = ret.loc[date:]
    rets_b = rets_b.loc[date:]
    df_rets.loc[names_single[i]] = sharpe_ratio(ret), information_rate(ret, rets_b), sortino_ratio(ret)
df_rets.loc["buyhold"] = sharpe_ratio(rets_b.loc[date:]), 0, sortino_ratio(rets_b.loc[date:])

        
        Copy
    

    	
            df_rets.sort_values(by = "Sharpe ratio", ascending = False)
        
        Copy
    

    	
            
import pickle # need to ' pip install pickle-mixin '


import inspect # Needed to check what object is or isn't a module object; pickle cannot 'pickle' module objects.

        
        Copy
    

    	
            
alll, alll_err, alll_names = [[], []], [], []
pickle_out = open("MCDropout_005_module_test.pickle","wb")
for i in dir():
    try:
        exec(f"pickle.dump({i}, pickle_out)")
        alll_names.append(i)
    except:
        alll_err.append(i)
pickle_out.close()
 
for i in alll_names:
    alll[0].append(i)
    exec(f"alll[1].append({i})")

        
        Copy
    

    	
            
# To save data out:
pickle_out = open("MCDropout_005.pickle","wb") # This creates the '.pickle' file where our data of choice will be saved. ' wb ' stand for 'write bytes'.
pickle.dump(alll, pickle_out) # ' pickle_out ' specifies the '.pickle' file in which we want to write (actually: overwrite - everything previously in that file will be overwritten)
pickle_out.close() # We need to close this '.pickle' file; leaving it open could corrupt it.
 
# # To load data in:
# pickle_in = open("MCDropout_005.pickle","rb") # ' rb ' stand for 'read bytes'.
# allll = pickle.load(pickle_in)
# pickle_in.close() # We ought to close the file we opened to allow any other programs access if they need it.

        
        Copy