<jupyter_start><jupyter_text>Predictive Modelling: XGBoost  Imports<jupyter_code>%load_ext autoreload
%autoreload 2

# Pandas and numpy
import pandas as pd
import numpy as np

#
from IPython.display import display, clear_output
import sys
import time

# Libraries for Visualization
import matplotlib.pyplot as plt
import seaborn as sns
from src.visualization.visualize import plot_corr_matrix, plot_multi, plot_norm_dist, plot_feature_importances

# Some custom tools
from src.data.tools import check_for_missing_vals

# Alpaca API
import alpaca_trade_api as tradeapi

# Pickle
import pickle
import os
from pathlib import Path

# To load variables from .env file into system environment
from dotenv import find_dotenv, load_dotenv

from atomm.Indicators import MomentumIndicators
from atomm.DataManager.main import MSDataManager
from atomm.Tools import calc_open_position
from src.visualization.visualize import plot_confusion_matrix
from atomm.Methods import BlockingTimeSeriesSplit, PurgedKFold           


import time, os

# scikit-learn
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix
from sklearn.model_selection import train_test_split, cross_val_score, TimeSeriesSplit
from xgboost import XGBClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression

# For BayesianHyperparameter Optimization
from src.models.hyperparameter_optimization import search_space, BayesianSearch

# Feature building tool
from src.features.build_features import calcIndicators


# Visualization libraries
import seaborn as sns
import matplotlib.pyplot as plt
from pandas.plotting import scatter_matrix
import matplotlib.gridspec as gridspec
#import matplotlib.style as style
from scipy import stats
# Load environment variables
load_dotenv(find_dotenv())<jupyter_output><empty_output><jupyter_text>Defining functions<jupyter_code>def run_combinations(
        symbol,
        forecast_horizon,
        input_window_size,
        X_train,
        X_test,
        y_train,
        y_test,
        prices,
        model,
        hyper_optimize=False,
        n_eval=10,
        silent = False,
    
    ):
    start = time.time()
    cnt = 1
    df_score = pd.DataFrame(index=forecast_horizon, columns=input_window_size, dtype=float)
    df_score.rename_axis('Forecast horizon (days)').rename_axis('Input window (days)', axis=1)
    df_returns = df_score.copy()
    total = len(forecast_horizon)*len(input_window_size)
    # Loop over forecast horizons
    try:
        for i, fh in enumerate(forecast_horizon):
            # Loop over input / lookback windows
            for j, iw in enumerate(input_window_size):

                params = {}
                if hyper_optimize:
                    _, _, params = BayesianSearch(
                        search_space(model),
                        model,
                        X_train[symbol][[f'{x}_{iw}' for x in ti_list]], 
                        y_train[symbol][f'signal_{fh}'], 
                        X_test[symbol][[f'{x}_{iw}' for x in ti_list]],
                        y_test[symbol][f'signal_{fh}'],
                        num_eval=n_eval,
                        silent=True,
                    )
                clf = model(**params)
                clf.fit(X_train[symbol][[f'{x}_{iw}' for x in ti_list]], y_train[symbol][f'signal_{fh}'])
                y_pred = clf.predict(X_test[symbol][[f'{x}_{iw}' for x in ti_list]])
                # Calculate accuracy score ()
                a_score = accuracy_score(y_test[symbol][f'signal_{fh}'], y_pred)
                # Calculate returns generated by stratedy
                returns = calc_returns(y_pred, prices[symbol])
                cum_returns = returns['Cum_Returns_Strat'][-1]
                df_score.iloc[i, j] = a_score
                df_returns.iloc[i, j] = cum_returns
                if silent == False:
                    clear_output(wait=True)
                    display(f'Leg [{cnt} / {total}]; Test Score {a_score}; Test Returns {cum_returns}')
                cnt += 1
    except (KeyboardInterrupt, SystemExit):
        sys.exit(0)
    print(f'Elapsed time {round(time.time()-start, 0)}s.')
    return df_score, df_returns
def avg_model(  
    symbol_list,
    forecast_horizon,                                 
    input_window_size,                                  
    X_train,    
    X_test,    
    y_train,    
    y_test,    
    prices_test,
    model,
    silent = False
):
    results_dict, returns_dict = {}, {}
    cnt = 1
    start = time.time()
    for symb in symbol_list:
        results, returns = run_combinations(
            symb,                             
            forecast_horizon,                                 
            input_window_size,                                  
            X_train,    
            X_test,    
            y_train,    
            y_test,    
            prices_test,   
            model=model,
            silent=True
        )
        results_dict[symb] = results
        returns_dict[symb] = returns
        avg_results = pd.DataFrame(
            np.array([results_dict[x].values for x in results_dict.keys()]).mean(axis=0),
            columns=results.columns, index=results.index).style.background_gradient(cmap='BuGn')
        
        if silent == False:
            clear_output(wait=True)
            display(f'Leg [{cnt}/{len(symbol_list)}]; Elapsed Time {round(time.time()-start , 0)}s\n', avg_results)
        cnt += 1

    return avg_results, results_dict, returns_dict
def create_two_class_signal(df, symbol_list):
    df_signals = pd.DataFrame(pd.MultiIndex.from_product([[], []]),)# index=df.index)
    sign_dic = {}
    for symbol in symbol_list:
        data = df[symbol]
        dicti = {}
        for n in forecast_horizon:
            diff = (data['Close'].shift(-n) - data['Close'])/data['Close']
            signal = np.where(diff > 0, 1, 0)
            dicti[f'signal_{n}'] = signal
        sign_df = pd.DataFrame.from_dict(dicti)
        sign_dic[symbol] = sign_df
    df_signals = pd.concat(sign_dic.values(), keys=sign_dic.keys(), axis=1)
    return df_signals
def create_three_class_signal(df, symbol_list, threshold=None, multi=1):
    df_signals = pd.DataFrame(pd.MultiIndex.from_product([[], []]),)# index=df.index)
    df_op = df_signals.copy()
    sign_dic, op_dic = {}, {}
    for symbol in symbol_list:
        data = df[symbol]
        dicti, dict_op = {}, {}
        for n, thresh in zip(forecast_horizon, threshold):
            diff = (data['Close'].shift(-n) - data['Close'])/data['Close']
            signal1p = np.where(diff > multi*thresh/100, 1, 0)
            signal1n = np.where(diff < -multi*thresh/100, -1, 0)
            signal = signal1p + signal1n
            open_position = calc_open_position(signal)
            dicti[f'signal_{n}'] = signal
            dict_op[f'open_position_{n}'] = open_position
        sign_dic[symbol] = pd.DataFrame.from_dict(dicti)
        op_dic[symbol] = pd.DataFrame.from_dict(dict_op)
    df_signals = pd.concat(sign_dic.values(), keys=sign_dic.keys(), axis=1)
    df_op = pd.concat(op_dic.values(), keys=op_dic.keys(), axis=1)
    return df_signals, df_op
def calc_returns(signals, prices):
    returns = prices[['Close']].pct_change()
    returns['Signal_Strat'] = calc_open_position(signals)
    returns['Cum_Returns_Strat'] = (returns['Signal_Strat'] * returns['Close']).cumsum()
    returns['Cum_Returns_BH'] = returns['Close'].cumsum()
    return returns<jupyter_output><empty_output><jupyter_text>Loading the data<jupyter_code>from src.data.tools import check_for_missing_vals
forecast_horizon = [1, 3, 5, 7, 10, 15, 20, 25, 30]
input_window_size = [3, 5, 7, 10, 15, 20, 25, 30]
threshold_list = [0.63, 1.15, 1.49, 1.79, 2.14, 2.65, 3.08, 3.48, 3.94]
ti_list = ['macd', 'rsi', 'stoc', 'roc', 'bbu', 'bbl', 'ema', 'atr', 'adx', 'cci', 'williamsr', 'stocd']
ti_list = ['sma', 'rsi', 'stoc', 'roc', 'ema', 'atr', 'adx', 'cci', 'williamsr', 'stocd']

symbol_list = df_moments.columns.get_level_values(0).unique()
symbol_list = ['GPS']
startd = datetime.datetime(2002, 1, 29)
endd = datetime.datetime(2012, 7, 30)
start = time.time()
df = pd.DataFrame(pd.MultiIndex.from_product([[], []]))
dic = {} 
#symb_list = dm.ReturnIndexConstituents('DJIA')+['DJIA']
for symb in symbol_list:
    #print(symb)
    symb_df = dm.ReturnData(
        symb,
        start_date=startd,
        end_date=endd,
        limit=None)
    dic[symb] = symb_df
    
prices = pd.concat(dic.values(), keys=dic.keys(), axis=1)
print(f'Elapsed time: {int(time.time()-start)}s')
check_for_missing_vals(prices)

def calcIndicators2(data, symbol, lookback_windows, ti_list):
    mi = MomentumIndicators(data[symbol])
    df = pd.DataFrame(index=data.index)
    ti_dict = {
        'sma': MomentumIndicators(data[symbol]).calcSMA,
        'macd': MomentumIndicators(data[symbol]).calcMACD,
           'rsi': MomentumIndicators(data[symbol]).calcRSI,
           'stoc': MomentumIndicators(data[symbol]).calcSTOC,
           'roc': MomentumIndicators(data[symbol]).calcROC,
           'bbu': MomentumIndicators(data[symbol]).calcBB,
           'bbl': MomentumIndicators(data[symbol]).calcBB,
           'ema': MomentumIndicators(data[symbol]).calcEMA,
           'atr': MomentumIndicators(data[symbol]).calcATR,
           'adx': MomentumIndicators(data[symbol]).calcADX,
           'cci': MomentumIndicators(data[symbol]).calcCCI,
           'williamsr': MomentumIndicators(data[symbol]).calcWR,
            }
    for n in input_window_size:
        for ti in ti_list:
            if ti == 'stocd':
                df[f'{ti}_{n}'] = mi.calcEMA(3, df[f'stoc_{n}'])
            elif ti == 'log_ret':
                df[f'{ti}_{n}'] = np.log(data[symbol]['Close']).diff(periods=n)
            elif ti in ['autocorr_1', 'autocorr_3', 'autocorr_5']:
                df[f'{ti}_{n}'] = df[f'log_ret_{n}'].rolling(
                    window=n,
                    min_periods=n,
                    center=False
                ).apply(lambda x: x.autocorr(lag=int(ti[-1])), raw=False)
            elif ti == 'vol':
                df[f'vol_{n}'] = df[f'log_ret_{n}'].rolling(
                    window=n,
                    min_periods=n,
                    center=False
                ).std()
            else:
                df[f'{ti}_{n}'] = ti_dict.get(ti)(n)

            # Shynkevich et al 2017 also has SMA, but not Bollinger Bands, MACD
    return df
cnt = 1
start = time.time()
total = len(symbol_list)
mom_dict = {}
ti_list = ['sma', 'rsi', 'stoc', 'roc', 'ema', 'atr', 'adx',
           'cci', 'williamsr', 'stocd', 'log_ret', 'autocorr_1', 'autocorr_3', 'autocorr_5', 'vol']

for symbol in symbol_list:
    res = calcIndicators(prices, symbol, input_window_size, ti_list)
    mom_dict[symbol] = res
    clear_output(wait=True)
    display(f'Leg [{cnt} / {total}]; Last symbol: {symbol}; Elapsed time: {round(time.time()-start, 0)}s')
    cnt += 1
df_moments = pd.concat(mom_dict.values(), keys=mom_dict.keys(), axis=1)
df_moments = df_moments.iloc[max(input_window_size):-max(forecast_horizon)]
df_moments['GPS']['autocorr_1_15'].plot()
check_for_missing_vals(df_moments)<jupyter_output>GPS
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]<jupyter_text>Two class labels<jupyter_code>y_2c = create_two_class_signal(prices, symbol_list)
y_2c = y_2c[max(input_window_size):-max(forecast_horizon)]
y_2c.shape
l1 = []
for y in y_2c.columns.get_level_values(0).unique():
    l = []
    for x in y_2c[y].columns:
        c = y_2c[y][x].value_counts()
        s = c[0] + c[1]
        l.append([c[0]/s, c[1]/s])
    l1.append(l)
pd.DataFrame(np.mean(np.array(l1), axis=0).T, index=pd.Index(['0', '1']), columns=y_2c.columns.get_level_values(1).unique())<jupyter_output><empty_output><jupyter_text>Three class labels<jupyter_code>#%%timeit

sym_list = symbol_list[:]

#dicti = {}
y_3c, y_3op = create_three_class_signal(prices, symbol_list, threshold=threshold_list)
y_3c_test, y_3op_test = create_three_class_signal(prices, symbol_list, threshold=threshold_list)

# Remove first and last 30 entries
y_3c = y_3c[max(input_window_size):-max(forecast_horizon)]
y_3op = y_3op[max(input_window_size):-max(forecast_horizon)]
y_3c.shape
l1 = []
for y in y_3c.columns.get_level_values(0).unique():
    l = []
    for x in y_3c[y].columns:
        c = y_3c[y][x].value_counts()
        s = c[0] + c[1] + c[-1]
        l.append([c[-1]/s, c[0]/s, c[1]/s])
    l1.append(l)
pd.DataFrame(np.mean(np.array(l1), axis=0).T, index=pd.Index(['-1', '0', '1']), columns=y_3c.columns.get_level_values(1).unique())<jupyter_output><empty_output><jupyter_text>Scaling the features<jupyter_code>from sklearn.preprocessing import MinMaxScaler, StandardScaler
#scale = MinMaxScaler()
scale = StandardScaler()
scaled = scale.fit_transform(df_moments)
scaled.shape
X_scaled = pd.DataFrame(data=scaled, columns=df_moments.columns, index=df_moments.index)
#X_scaled = X<jupyter_output><empty_output><jupyter_text>Triple-barrier labelling methods<jupyter_code>from atomm.Filters import cusum_filter
from atomm.Tools import get_daily_vol
data = pd.read_csv('../data/raw/dollar_bars.txt')
data.index = pd.to_datetime(data['date_time'])
data = data.drop('date_time', axis=1)
data.index
dd2 = get_daily_vol(prices['GPS'].Close)
dd2.head(50)
(pd.Series(prices.index[dd - 1], index=prices.index[prices.shape[0] - dd.shape[0]:])).index
close = prices['GPS'][['Close']]
                      
dd.shape<jupyter_output><empty_output><jupyter_text>Train-Test Split<jupyter_code>X_scaled = df_moments.copy()
# Use 70/30 train/test splits
test_p = .15
# Scaled, three-class
test_size = int((1 - test_p) * X_scaled.shape[0])
X_train, X_test, y_train, y_test = X_scaled[:test_size], X_scaled[test_size:], y_3c[:test_size], y_3c[test_size:]
prices_train, prices_test = prices[:test_size], prices[test_size:]
# Unscaled, two-class
test_size = int((1 - test_p) * X_scaled.shape[0])
X_train, X_test, y_train, y_test = X_scaled[:test_size], X_scaled[test_size:], y_2c[:test_size], y_2c[test_size:]
prices_train, prices_test = djia[:test_size], djia[test_size:]
X_test.index.min()
# Scaled, two-class
test_size = int((1 - test_p) * X.shape[0])
X_train, X_test, y_train, y_test = X_scaled[:test_size], X_scaled[test_size:], y_2c[:test_size], y_2c[test_size:]
prices_train, prices_test = djia[:test_size], djia[test_size:]
#test_size = test_p
#X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_3c, test_size=test_size, random_state=101)<jupyter_output><empty_output><jupyter_text>Model<jupyter_code>symbol = 'GPS'
n = 15
# set up cross validation splits
tscv = TimeSeriesSplit(n_splits=5)
btscv = BlockingTimeSeriesSplit(n_splits=5)
#ppcv = PurgedKFold(n_splits=5)<jupyter_output><empty_output><jupyter_text>Single lookback/lookahead combination  RF<jupyter_code>n_estimator = 1000
depth = 2

rf = RandomForestClassifier(max_depth=depth, n_estimators=n_estimator,
                            criterion='entropy',
                            class_weight='balanced_subsample',
                            random_state=101)
rf.fit(X_train[symbol][[f'{x}_{n}' for x in ti_list]], y_train[symbol][f'signal_{n}'])

y_pred_rf = rf.predict(X_test[symbol][[f'{x}_{n}' for x in ti_list]])
print(classification_report(y_test[symbol][f'signal_{n}'], y_pred_rf))
print(confusion_matrix(y_test[symbol][f'signal_{n}'], y_pred_rf))
plot_feature_importances(
    rf.feature_importances_,                         
    columns=X_test[symbol][[f'{x}_{n}' for x in ti_list]].columns,
    model='RF',
    top_count=30,
)<jupyter_output><empty_output><jupyter_text>SVC<jupyter_code>clf = SVC(C=1000, gamma=1E-5)
clf.fit(X_train[symbol][[f'{x}_{n}' for x in ti_list]], y_train[symbol][f'signal_{n}'])

y_pred_svc = clf.predict(X_test[symbol][[f'{x}_{n}' for x in ti_list]])
print(classification_report(y_test[symbol][f'signal_{n}'], y_pred_svc))
print(confusion_matrix(y_test[symbol][f'signal_{n}'], y_pred_svc))
#plot_feature_importances(svc.coef_,
#                         X_scaled[symbol][[f'{x}_{n}' for x in ti_list]].columns, 
#                         model='SVC1', top_count=100)
clf = XGBClassifier(n_jobs=-1)
clf.fit(X_train[symbol][[f'{x}_{n}' for x in ti_list]], y_train[symbol][f'signal_{n}'])
y_pred_xgb = clf.predict(X_test[symbol][[f'{x}_{n}' for x in ti_list]])
print(classification_report(y_test[symbol][f'signal_{n}'], y_pred_xgb))
print(confusion_matrix(y_test[symbol][f'signal_{n}'], y_pred_xgb))
#plot_feature_importances(svc.coef_,
#                         X_scaled[symbol][[f'{x}_{n}' for x in ti_list]].columns, 
#                         model='SVC1', top_count=100)<jupyter_output>precision    recall  f1-score   support

           0       0.41      0.48      0.44       159
           1       0.59      0.52      0.56       229

    accuracy                           0.51       388
   macro avg       0.50      0.50      0.50       388
weighted avg       0.52      0.51      0.51       388

[[ 76  83]
 [109 120]]<jupyter_text>All combinations<jupyter_code>clf = XGBClassifier
score_xgb, returns_xgb = run_combinations(
    symbol=symbol,
    forecast_horizon=forecast_horizon,
    input_window_size=input_window_size,
    X_train=X_train,
    X_test=X_test,
    y_train=y_train,
    y_test=y_test,
    prices=prices_test,
    model=clf,
    hyper_optimize=True,
    n_eval=100
)

plot_corr_matrix(score_xgb, mask_upper=False, show_annot=True, figsize=(8, 8))<jupyter_output>100%|██████████| 100/100 [02:13<00:00,  1.33s/it, best loss: -0.5470588235294118]<jupyter_text>Hyperparamter Optimization: Bayesian Optimization, SVC<jupyter_code>n1=15
n2=15

model = SVC
score, clf_best, parameters_best = BayesianSearch(
    search_space(model),
    model,
    X_train[symbol][[f'{x}_{n1}' for x in ti_list]], 
    y_train[symbol][f'signal_{n2}'], 
    X_test[symbol][[f'{x}_{n1}' for x in ti_list]],
    y_test[symbol][f'signal_{n2}'],
    num_eval=500
)
clf_best<jupyter_output><empty_output><jupyter_text>Hyperparamter Optimization: Bayesian Optimization, XGB<jupyter_code>n1=15
n2=15

model = XGBClassifier
bsearch_xgb = BayesianSearch(
    search_space(XGBClassifier),
    model,
    X_train[symbol][[f'{x}_{n1}' for x in ti_list]], 
    y_train[symbol][f'signal_{n2}'], 
    X_test[symbol][[f'{x}_{n1}' for x in ti_list]],
    y_test[symbol][f'signal_{n2}'],
    num_eval=10
)<jupyter_output>100%|██████████| 10/10 [01:11<00:00,  7.19s/it, best loss: -0.5795454545454546]
##### Results #####
Score best parameters:  -0.5795454545454546
Best parameters:  {'booster': 'gbtree', 'colsample_bytree': 0.6571346286826318, 'cv': 'btscv', 'gamma': 5, 'learning_rate': 43.51596522276907, 'max_depth': 179, 'model': <class 'xgboost.sklearn.XGBClassifier'>, 'n_estimators': 116, 'n_jobs': -1, 'subsample': 1}
Test Score:  0.5902061855670103
Parameter combinations evaluated:  10
Time elapsed:  72.14264488220215<jupyter_text>Hyperparamter Optimization: Bayesian Optimization, RF<jupyter_code>n1=15
n2=15
num_eval=15

model = RandomForestClassifier
score, clf_best, parameters_best = BayesianSearch(
    search_space(model),
    model,
    X_train[symbol][[f'{x}_{n1}' for x in ti_list]], 
    y_train[symbol][f'signal_{n2}'], 
    X_test[symbol][[f'{x}_{n1}' for x in ti_list]],
    y_test[symbol][f'signal_{n2}'],
    num_eval=num_eval
)
plot_feature_importances(
    clf_best.feature_importances_,                         
    columns=X_test[symbol][[f'{x}_{n}' for x in ti_list]].columns,
    model='RF w/ BayesianOptimization',
    top_count=30,
)
## 
n1=15
n2=15
num_eval=100

model = KNeighborsClassifier
bsearch_svm = BayesianSearch(
    search_space(model),
    model,
    X_train[symbol][[f'{x}_{n1}' for x in ti_list]], 
    y_train[symbol][f'signal_{n2}'], 
    X_test[symbol][[f'{x}_{n1}' for x in ti_list]],
    y_test[symbol][f'signal_{n2}'],
    num_eval=num_eval
)<jupyter_output>100%|██████████| 100/100 [01:01<00:00,  1.62it/s, best loss: -0.5802318801952324]
##### Results #####
Score best parameters:  -0.5802318801952324
Best parameters:  {'cv': 'cv', 'model': <class 'sklearn.neighbors._classification.KNeighborsClassifier'>, 'n_jobs': -1, 'n_neighbors': 61, 'normalize': 0, 'p': 8, 'scale': 0}
Test Score:  0.4882903981264637
Parameter combinations evaluated:  100
Time elapsed:  62.092283964157104<jupyter_text>Running on all 50 stocks<jupyter_code>best_params = {'bootstrap': False, 'criterion': 'gini', 'max_depth': 218, 'max_features': 1, 'min_samples_leaf': 19, 'n_estimators': 423}
model_2a = RandomForestClassifier(n_jobs=-1, **best_params)
avg, _, _ = avg_model(  
    symbol_list,
    forecast_horizon,                                 
    input_window_size,                                  
    X_train,    
    X_test,    
    y_train,    
    y_test,    
    prices_test,
    model=model_2a,
    silent = False
)
avg<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1401011868173436, 'fraction_numerical': 0.033807745047047144, 'mean_word_length': 3.359027205276175, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '341743', 'license': 'MIT', 'n_tokens_mistral': 8288, 'n_tokens_neox': 7732, 'n_words': 1419, 'path': 'notebooks/ressources/Insanity_Check.ipynb', 'repo_name': 'robindoering86/capstone_nf', 'size': '21149'}