!pip install pmdarima

import os
import io
import math
import random
import requests
from tqdm import tqdm
import numpy as np
import pandas as pd
import sklearn
import matplotlib.dates as mdates
import datetime as dt
from sklearn.model_selection import train_test_split
#from sklearn.metrics import r2_score

import matplotlib.pyplot as plt
%matplotlib inline

plt.style.use('seaborn-colorblind')

import warnings
warnings.filterwarnings("ignore")

"""
Explanation: QUESTION: Should we buy or sell in 3 days?
The idea behind this mini project:

Create data features from the price of the day
Check the performance with different models
ARIMA
Random Forest
LSTM
Obtain a model that can predict in some manner how the price is going to move in the following three days.

Import Libraries
End of explanation
"""


def download_file(url, filename):
    r = requests.get(url, stream=True)

    total_size = int(r.headers.get('content-length', 0)); 
    block_size = 1024
    total_kb_size = math.ceil(total_size//block_size)
    
    wrote = 0 
    with open(filename, 'wb') as f:
        for data in tqdm(r.iter_content(block_size), total=total_kb_size , unit='KB', unit_scale=True):
            wrote = wrote  + len(data)
            f.write(data)

datafile = "eth-eur.csv"

#import from server
if not os.path.exists(datafile):
    download_file("https://www.coingecko.com/price_charts/export/279/eur.csv", datafile)

data = pd.read_csv(datafile)

# Only we take care of the last year to date
data = data[-365:]
data.info()

data.snapped_at = data.snapped_at.apply(lambda x: x.split()[0])
data.set_index('snapped_at', inplace=True)
data.index = pd.to_datetime(data.index)

"""
Explanation: Getting the Data
End of explanation
"""


import statsmodels.api as sm
seas_d = sm.tsa.seasonal_decompose(data.price, model='additive')

fig = seas_d.plot()
fig.set_figheight(10)
fig.set_figwidth(15)

seas_d = sm.tsa.seasonal_decompose(data.price, model='multiplicative')

fig = seas_d.plot()
fig.set_figheight(10)
fig.set_figwidth(15)

"""
Explanation: Analysis
Decomposing Time Series into several components- Trend, Seasonality and Random noise
The additive model is Y[t] = T[t] + S[t] + e[t]
The multiplicative model is Y[t] = T[t] * S[t] * e[t]
End of explanation
"""


# Multiplicative fits better, so changing to log to convert mult by sum
data['log_price'] = np.log(data.price)

"""
Explanation: Insights:
* Multiplicative fits better
* There seems to be a week seasonality
End of explanation
"""


data['future'] = data.price.shift(-3)
data['log_future'] = data.log_price.shift(-3)
data = data[:-3]
data.tail()

"""
Explanation: Data Engineering
End of explanation
"""


#Checking if we have NaN in another place
data[data.isnull().any(axis=1)].count()

#Fill the blanks with the mean between the previous and the day after
filled = False

for idx in data.index:
    dayloc = data.index.get_loc(idx)
    day = data.loc[idx]
    if day.hasnans:
        #updating
        rg = slice(dayloc-1, dayloc+2)
        data.loc[idx] = data.iloc[rg].mean()
        filled = True
        print("Day <{}> filled".format(idx))

if not filled:
    print("Nothing to fill")

"""
Explanation: Filling de NaNs
End of explanation
"""


for i in [1,2,3,5,7,10]:
    column = "diff{}".format(i)
    data[column] = data.log_price.diff(i)

data.dropna(inplace=True)

"""
Explanation: Adding Difference Features
End of explanation
"""


# 1 -> up, 0 -> down
data['move'] = data.diff1.apply(lambda x: 1 if  x >= 0 else 0)

"""
Explanation: Adding Movement Feature
End of explanation
"""


data.tail()

"""
Explanation: Final Dataset
End of explanation
"""


# 10 days for validation, 10 days for for testing and the rest for training
idx1 = -20
idx2 = -10
data_train, data_val, data_test = data[:idx1], data[idx1:idx2], data[idx2:]

print("Size data_train: {}".format(data_train.shape[0]))
print("Size data_val: {}".format(data_val.shape[0]))
print("Size data_test: {}".format(data_test.shape[0]))

plt.figure(figsize=(10,7))
plt.plot(data_train.log_future, label='Training Data')
plt.plot(data_val.log_future, 'orange', label='Validation Data')
plt.plot(data_test.log_future, 'red', label='Testing Data')
plt.legend()

data_train.corr()

plt.plot(data_train.diff1)

"""
Explanation: Training, Validation and Testing Datasets
End of explanation
"""


X_train = data_train.log_price
X_val = data_val.log_price
y_val = data_val.log_future

import pmdarima as pm

arima = pm.auto_arima(X_train,
                      test='adf',
                      #d=1,
                      #m = 1, # week seasonality
                      seasonal=True,
                      trace=True,
                      error_action='ignore',  
                      suppress_warnings=True, 
                      stepwise=True) #False)

arima.summary()

start = len(X_train)
end = start + len(X_val)

# Rolling step with re-estimation
pred = pd.Series()
for i in range(start, end):
    yhat = arima.predict(n_periods=3)[-1]
    pred[data.index[i]] = yhat

    # updating the model
    arima = arima.update(data.log_price[i])

def mean_absolute_percentage_error(y_true, y_pred): 
    y_true, y_pred = np.array(y_true), np.array(y_pred)
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100


y_val = np.exp(y_val)
pred = np.exp(pred)
pred1 = pred.copy()

error = mean_absolute_percentage_error(y_val, pred)
print("{} Validation MAPE: {:.3f}".format('ARIMA', error))

# plot
plt.figure(figsize=(13,7))
plt.plot(y_val, label='Actual')
plt.plot(pred, color='red', label='Prediction')
plt.legend()
plt.show()

"""
Explanation: It could be interesing to use the first order of differencing to remove the trend. There're also some outliers that should be removed.
ARIMA
Just in case, we're trying to find any other best order using 'auto ARIMA' from pmdarima library.
End of explanation
"""


X_train = data_train[['log_price', 'diff1','diff2','diff3','diff5','diff7','diff10','move']]
y_train = data_train.log_future

from sklearn.ensemble import RandomForestRegressor

rf = RandomForestRegressor(random_state=0, criterion='mae')
rf.fit(X_train, y_train)


importances = rf.feature_importances_
indices = np.argsort(importances)

plt.title('Feature Importances')
plt.barh(range(len(indices)), importances[indices], color='g', align='center')
plt.yticks(range(len(indices)), X_train.columns[indices])
plt.xlabel('Relative Importance')


# Define RF features to be used
RF_FEATURES = X_train.columns[indices[-5:]]

# Refit the model
rf.fit(X_train[RF_FEATURES], y_train)


X_val = data_val[RF_FEATURES]
y_val = data_val.log_future
pred = rf.predict(X_val)

y_val = np.exp(y_val)
pred = np.exp(pred)
pred2 = pred.copy()

error = mean_absolute_percentage_error(y_val, pred)
print("Random Forest Validation MAPE: {:.3f}".format(error))

# plot
plt.figure(figsize=(13,7))
plt.plot(y_val, label='Actual')
plt.plot(data_val.index, pred, color='red', label='Prediction')
plt.legend()
plt.show()

"""
Explanation: Random Forest
End of explanation
"""


from distutils.version import LooseVersion
import warnings
import tensorflow as tf

# Check TensorFlow Version
assert LooseVersion(tf.__version__) >= LooseVersion('1.0'), 'Please use TensorFlow version 1.0 or newer'
print('TensorFlow Version: {}'.format(tf.__version__))

# Check for a GPU
if not tf.test.gpu_device_name():
    warnings.warn('No GPU found. Please use a GPU to train your neural network.')
else:
    print('Default GPU Device: {}'.format(tf.test.gpu_device_name()))

"""
Explanation: LSTM
End of explanation
"""


FEATURES_LEN = len(RF_FEATURES)

from sklearn.preprocessing import StandardScaler

X_train = data_train[RF_FEATURES]
y_train = data_train.log_future

X_val = data_val[RF_FEATURES]
y_val = data_val.log_future

# Scale the features to keep all of them in the same range
scaler = StandardScaler()
X_train[X_train.columns] = scaler.fit_transform(X_train[X_train.columns])
X_val[X_val.columns] = scaler.transform(X_val[X_val.columns])

# Adapt to LSTM sequence shape
X_train_reshaped = X_train.values.reshape( (X_train.shape[0], 1 , FEATURES_LEN) )
X_val_reshaped = X_val.values.reshape( (X_val.shape[0], 1 , FEATURES_LEN) )

X_train.describe()

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import Dropout
from keras.callbacks import EarlyStopping

lstm = Sequential()
lstm.add(LSTM(64, activation="tanh", input_shape=(1, FEATURES_LEN) ))
lstm.add(Dropout(0.2))
lstm.add(Dense(1))
lstm.compile(loss='mean_absolute_percentage_error', optimizer='adam')

es_callback = EarlyStopping(monitor="loss", min_delta=0, patience=5, restore_best_weights=True)
history = lstm.fit(X_train_reshaped, y_train, validation_data=(X_val_reshaped, y_val), shuffle=False, epochs=100, batch_size=32, verbose=1, callbacks=[es_callback])

def visualize_loss(history):
    loss = history.history["loss"]
    val_loss = history.history["val_loss"]
    epochs = range(len(loss))
    plt.figure(figsize=(8,4.5))
    plt.plot(epochs, loss, label="Training loss")
    plt.plot(epochs, val_loss, label="Validation loss")
    plt.title("Training and Validation Loss")
    plt.xlabel("Epochs")
    plt.ylabel("Loss")
    plt.legend()

def visualize_plot(model, X, y, title):
    index = y.index
    pred = model.predict(X)
    y = np.exp(y)
    pred = np.exp(pred)
    plt.figure(figsize=(8,4.5))
    plt.plot(index, y, label='Actual')
    plt.plot(index, pred, "red", label='Prediction')
    plt.title(title)
    plt.xlabel("Dates")
    plt.ylabel("Price (€)")
    plt.legend()
    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))
    plt.gca().xaxis.set_major_locator(mdates.DayLocator(interval=15))
    plt.gcf().autofmt_xdate()

def print_mape(model, X, y, title):
    pred = model.predict(X)
    y = np.exp(y)
    pred = np.exp(pred)
    error = mean_absolute_percentage_error(y, pred)
    print("{} MAPE in LSTM: {:.2f}".format(title, error))

visualize_loss(history)
visualize_plot(lstm, X_val_reshaped, y_val, "Validation Data")
print_mape(lstm, X_val_reshaped, y_val, 'Validation')

"""
Explanation: The same dataset with differences is going to be used as before (Random Forest model) but now it's going to split the training data into a validation set as well.
End of explanation
"""


start = len(data_train) + len(data_val)
end = start + len(data_test)

## ARIMA
# Rolling step with re-estimation
pred1 = pd.Series()
for i in range(start, end):
    yhat = arima.predict(n_periods=3)[-1]
    pred1[data.index[i]] = yhat
    arima = arima.update(data.log_price[i])

## Random Forest
pred2 = rf.predict( data_test[RF_FEATURES] )

## LSTM
X_test = data_test[RF_FEATURES]
X_test[X_test.columns] = scaler.transform(X_test[X_test.columns])
X_test = X_test.values.reshape((X_test.shape[0], 1, FEATURES_LEN))
pred3 = lstm.predict(X_test)

#
pred1 = np.exp(pred1)
pred2 = np.exp(pred2)
pred3 = np.exp(pred3)

#
mape1 = mean_absolute_percentage_error(data_test.future, pred1)
mape2 = mean_absolute_percentage_error(data_test.future, pred2)
mape3 = mean_absolute_percentage_error(data_test.future, pred3)

#
print("Arima MAPE: {:.2f}".format(mape1))
print("Random Forest MAPE: {:.2f}".format(mape2))
print("LSTM MAPE: {:.2f}".format(mape3))

##
plt.figure(figsize=(16,9))
plt.plot(data_test.index, data_test.future, label='Actual')
plt.plot(data_test.index, pred1, "green", linestyle='dashed', label='Arima')
plt.plot(data_test.index, pred2, "orange", linestyle='dashed', label='Random Forest')
plt.plot(data_test.index, pred3, "red", linestyle='dashed', label='LSTM')
plt.xlabel("Dates")
plt.ylabel("Price (€)")
plt.legend()
plt.show()

"""
Explanation: Testing the Models
End of explanation
"""


