technical_indicators_clustering.py

# Import dependencies
import pickle
import datetime as dt
import os
import bs4 as bs
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pandas_datareader as pdr
import requests
import seaborn as sns
from IPython.display import clear_output
from scipy.stats import mstats
from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestClassifier
from sklearn.mixture import GaussianMixture
from sklearn.model_selection import (
    GridSearchCV,
    RandomizedSearchCV,
    TimeSeriesSplit,
    validation_curve,
)

# Obtain list of S&P 100 companies from Wikipedia
resp = requests.get("https://en.wikipedia.org/wiki/S%26P_100")
convert_soup = bs.BeautifulSoup(resp.text, "lxml")
table = convert_soup.find("table", {"class": "wikitable sortable"})
tickers = [rows.findAll("td")[0].text.strip() for rows in table.findAll("tr")[1:]]

all_data = pd.DataFrame()
test_data = pd.DataFrame()
no_data = []

# Extract data from Yahoo Finance
for i in tickers:
    try:
        test_data = pdr.get_data_yahoo(
            i, start=dt.datetime(1990, 1, 1), end=dt.date.today()
        )
        test_data["symbol"] = i
        all_data = all_data.append(test_data)
        clear_output(wait=True)
        print(i)
    except:
        no_data.append(i)
    clear_output(wait=True)

all_data["Close_Shifted"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.shift(-6)
)
all_data["Target"] = (
    (all_data["Close_Shifted"] - all_data["Open"]) / all_data["Open"] * 100
).shift(-1)
all_data["Target_Direction"] = np.where(all_data["Target"] > 0, 1, 0)
all_data = all_data.dropna().copy()

all_data["SMA_5"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.rolling(window=5).mean()
)
all_data["SMA_15"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.rolling(window=15).mean()
)
all_data["SMA_ratio"] = all_data["SMA_15"] / all_data["SMA_5"]
all_data["SMA5_Volume"] = all_data.groupby("symbol")["Volume"].transform(
    lambda x: x.rolling(window=5).mean()
)
all_data["SMA15_Volume"] = all_data.groupby("symbol")["Volume"].transform(
    lambda x: x.rolling(window=15).mean()
)
all_data["SMA_Volume_Ratio"] = all_data["SMA5_Volume"] / all_data["SMA15_Volume"]

# Wilder Smoothing
def Wilder(data, periods):
    start = np.where(~np.isnan(data))[0][0]  # Check if nans present in beginning
    Wilder = np.array([np.nan] * len(data))
    Wilder[start + periods - 1] = data[
        start : (start + periods)
    ].mean()  # Simple Moving Average
    for i in range(start + periods, len(data)):
        Wilder[i] = (
            Wilder[i - 1] * (periods - 1) + data[i]
        ) / periods  # Wilder Smoothing
    return Wilder


all_data["prev_close"] = all_data.groupby("symbol")["Close"].shift(1)
all_data["TR"] = np.maximum(
    (all_data["High"] - all_data["Low"]),
    np.maximum(
        abs(all_data["High"] - all_data["prev_close"]),
        abs(all_data["prev_close"] - all_data["Low"]),
    ),
)
for i in all_data["symbol"].unique():
    TR_data = all_data[all_data.symbol == i].copy()
    all_data.loc[all_data.symbol == i, "ATR_5"] = Wilder(TR_data["TR"], 5)
    all_data.loc[all_data.symbol == i, "ATR_15"] = Wilder(TR_data["TR"], 15)

all_data["ATR_Ratio"] = all_data["ATR_5"] / all_data["ATR_15"]
all_data["prev_high"] = all_data.groupby("symbol")["High"].shift(1)
all_data["prev_low"] = all_data.groupby("symbol")["Low"].shift(1)

all_data["+DM"] = np.where(
    ~np.isnan(all_data.prev_high),
    np.where(
        (all_data["High"] > all_data["prev_high"])
        & (
            (
                (all_data["High"] - all_data["prev_high"])
                > (all_data["prev_low"] - all_data["Low"])
            )
        ),
        all_data["High"] - all_data["prev_high"],
        0,
    ),
    np.nan,
)

all_data["-DM"] = np.where(
    ~np.isnan(all_data.prev_low),
    np.where(
        (all_data["prev_low"] > all_data["Low"])
        & (
            (
                (all_data["prev_low"] - all_data["Low"])
                > (all_data["High"] - all_data["prev_high"])
            )
        ),
        all_data["prev_low"] - all_data["Low"],
        0,
    ),
    np.nan,
)

for i in all_data["symbol"].unique():
    ADX_data = all_data[all_data.symbol == i].copy()
    all_data.loc[all_data.symbol == i, "+DM_5"] = Wilder(ADX_data["+DM"], 5)
    all_data.loc[all_data.symbol == i, "-DM_5"] = Wilder(ADX_data["-DM"], 5)
    all_data.loc[all_data.symbol == i, "+DM_15"] = Wilder(ADX_data["+DM"], 15)
    all_data.loc[all_data.symbol == i, "-DM_15"] = Wilder(ADX_data["-DM"], 15)

all_data["+DI_5"] = (all_data["+DM_5"] / all_data["ATR_5"]) * 100
all_data["-DI_5"] = (all_data["-DM_5"] / all_data["ATR_5"]) * 100
all_data["+DI_15"] = (all_data["+DM_15"] / all_data["ATR_15"]) * 100
all_data["-DI_15"] = (all_data["-DM_15"] / all_data["ATR_15"]) * 100

all_data["DX_5"] = np.round(
    abs(all_data["+DI_5"] - all_data["-DI_5"])
    / (all_data["+DI_5"] + all_data["-DI_5"])
    * 100
)

all_data["DX_15"] = np.round(
    abs(all_data["+DI_15"] - all_data["-DI_15"])
    / (all_data["+DI_15"] + all_data["-DI_15"])
    * 100
)

for i in all_data["symbol"].unique():
    ADX_data = all_data[all_data.symbol == i].copy()
    all_data.loc[all_data.symbol == i, "ADX_5"] = Wilder(ADX_data["DX_5"], 5)
    all_data.loc[all_data.symbol == i, "ADX_15"] = Wilder(ADX_data["DX_15"], 15)

all_data["Lowest_5D"] = all_data.groupby("symbol")["Low"].transform(
    lambda x: x.rolling(window=5).min()
)
all_data["High_5D"] = all_data.groupby("symbol")["High"].transform(
    lambda x: x.rolling(window=5).max()
)
all_data["Lowest_15D"] = all_data.groupby("symbol")["Low"].transform(
    lambda x: x.rolling(window=15).min()
)
all_data["High_15D"] = all_data.groupby("symbol")["High"].transform(
    lambda x: x.rolling(window=15).max()
)

all_data["Stochastic_5"] = (
    (all_data["Close"] - all_data["Lowest_5D"])
    / (all_data["High_5D"] - all_data["Lowest_5D"])
) * 100
all_data["Stochastic_15"] = (
    (all_data["Close"] - all_data["Lowest_15D"])
    / (all_data["High_15D"] - all_data["Lowest_15D"])
) * 100

all_data["Stochastic_%D_5"] = all_data["Stochastic_5"].rolling(window=5).mean()
all_data["Stochastic_%D_15"] = all_data["Stochastic_5"].rolling(window=15).mean()

all_data["Stochastic_Ratio"] = (
    all_data["Stochastic_%D_5"] / all_data["Stochastic_%D_15"]
)


all_data["Diff"] = all_data.groupby("symbol")["Close"].transform(lambda x: x.diff())
all_data["Up"] = all_data["Diff"]
all_data.loc[(all_data["Up"] < 0), "Up"] = 0

all_data["Down"] = all_data["Diff"]
all_data.loc[(all_data["Down"] > 0), "Down"] = 0
all_data["Down"] = abs(all_data["Down"])

all_data["avg_5up"] = all_data.groupby("symbol")["Up"].transform(
    lambda x: x.rolling(window=5).mean()
)
all_data["avg_5down"] = all_data.groupby("symbol")["Down"].transform(
    lambda x: x.rolling(window=5).mean()
)

all_data["avg_15up"] = all_data.groupby("symbol")["Up"].transform(
    lambda x: x.rolling(window=15).mean()
)
all_data["avg_15down"] = all_data.groupby("symbol")["Down"].transform(
    lambda x: x.rolling(window=15).mean()
)

all_data["RS_5"] = all_data["avg_5up"] / all_data["avg_5down"]
all_data["RS_15"] = all_data["avg_15up"] / all_data["avg_15down"]

all_data["RSI_5"] = 100 - (100 / (1 + all_data["RS_5"]))
all_data["RSI_15"] = 100 - (100 / (1 + all_data["RS_15"]))

all_data["RSI_ratio"] = all_data["RSI_5"] / all_data["RSI_15"]

all_data["5Ewm"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.ewm(span=5, adjust=False).mean()
)
all_data["15Ewm"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.ewm(span=15, adjust=False).mean()
)
all_data["MACD"] = all_data["15Ewm"] - all_data["5Ewm"]

all_data["15MA"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.rolling(window=15).mean()
)
all_data["SD"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.rolling(window=15).std()
)
all_data["upperband"] = all_data["15MA"] + 2 * all_data["SD"]
all_data["lowerband"] = all_data["15MA"] - 2 * all_data["SD"]

all_data["RC"] = all_data.groupby("symbol")["Close"].transform(
    lambda x: x.pct_change(periods=15)
)

Target_variables = [
    "SMA_ratio",
    "ATR_5",
    "ATR_15",
    "ATR_Ratio",
    "ADX_5",
    "ADX_15",
    "SMA_Volume_Ratio",
    "Stochastic_5",
    "Stochastic_15",
    "Stochastic_Ratio",
    "RSI_5",
    "RSI_15",
    "RSI_ratio",
    "MACD",
]
for variable in Target_variables:
    all_data.loc[:, variable] = mstats.winsorize(
        all_data.loc[:, variable], limits=[0.1, 0.1]
    )

# Extract the returns
all_data["return"] = all_data["Close"].pct_change()
returns = all_data[["symbol", "return"]].copy()
returns["Date"] = returns.index.copy()

# Pivot the returns to create series of returns for each stock
transposed = returns.pivot(index="Date", columns="symbol", values="return")

# Transpose the data to get companies on the index level and dates on the column level since clusters takes place on index level
X = transposed.dropna().transpose()

# Extract sum of squares for K-means clusters from 1 to 50 clusters
sum_of_sq = np.zeros([50, 1])
for k in range(1, 51):
    sum_of_sq[k - 1] = KMeans(n_clusters=k).fit(X).inertia_

plt.plot(range(1, 50), sum_of_sq[1:50])
plt.title("Elbow Method")
plt.xlabel("Number of Cluster")
plt.ylabel("Within-cluster Sum of Squares")
plt.show()

print(pd.DataFrame(sum_of_sq, columns=["Difference in SS"], index=range(1, 51)).diff())

# Get 17 clusters
gmm = GaussianMixture(n_components=17)
gmm.fit(transposed.dropna().transpose())

# Predict for each company
clusters = gmm.predict(transposed.dropna().transpose())
clusters_df = pd.DataFrame({"Cluster": clusters, "Companies": transposed.columns})

# Sort by Clusters
clusters_df = clusters_df.sort_values(["Cluster"]).reset_index(drop=True)
print(clusters_df)

# Save as csv
clusters_df.to_csv("clusters.csv")
clusters_df = pd.read_csv("clusters.csv", index_col=0)

all_data.index = pd.to_datetime(all_data.index)

train_data = all_data.loc[
    :"2018-12-31",
]
test_data = all_data.loc["2019-01-01":]

# Separate between X and Y
X_train = train_data.loc[:, Target_variables]

Y_train = train_data.loc[:, ["Target_Direction"]]

# Create validation curve for the Random Forest Classifier
rf = RandomForestClassifier()
train_scoreNum, test_scoreNum = validation_curve(
    rf,
    X=X_train["2010-01-01":],
    y=Y_train.loc["2010-01-01":, "Target_Direction"],
    param_name="n_estimators",
    param_range=[3, 4, 7, 10, 12, 15, 20, 25, 30],
    cv=TimeSeriesSplit(n_splits=3),
)

train_scores_mean = np.mean(train_scoreNum, axis=1)
train_scores_std = np.std(train_scoreNum, axis=1)
test_scores_mean = np.mean(test_scoreNum, axis=1)
test_scores_std = np.std(test_scoreNum, axis=1)

plt.figure(figsize=(15, 10))
plt.plot([3, 4, 7, 10, 12, 15, 20, 25, 30], train_scores_mean)
plt.plot([3, 4, 7, 10, 12, 15, 20, 25, 30], test_scores_mean)
plt.legend(["Train Score", "Test Score"], fontsize="large")
plt.title("Validation Curve Score for n_estimators", fontsize="large")
plt.show()

# Run the loop for every unique cluster - 17 loops
for cluster_selected in clusters_df.Cluster.unique():
    print(f"The current cluster running is : {cluster_selected}")

    # Get data for that cluster
    co_data = all_data[
        all_data.symbol.isin(
            clusters_df.loc[
                clusters_df.Cluster == cluster_selected, "Companies"
            ].tolist()
        )
    ].copy()
    co_train = co_data[:"2018-12-31"]
    co_train = co_train.dropna().copy()

    X_train = co_train.loc[:, Target_variables]

    Y_train = co_train.loc[:, ["Target_Direction"]]

    # Define paramters from Validation Curve
    params = {
        "max_depth": [5, 7],
        "max_features": ["sqrt"],
        "min_samples_leaf": [10, 15, 20],
        "n_estimators": [5, 7, 9],
        "min_samples_split": [20, 25, 30],
    }  # Using Validation Curves

    rf = RandomForestClassifier()

    # Perform a TimeSeriesSplit on the dataset
    time_series_split = TimeSeriesSplit(n_splits=3)

    rf_cv = GridSearchCV(rf, params, cv=time_series_split, n_jobs=-1, verbose=20)

    # Fit the random forest with our X_train and Y_train
    rf_cv.fit(X_train, Y_train)

    # Save the fited variable into a Pickle file
    file_loc = f"{os.getcwd()}\\Pickle_Files\\Cluster_{cluster_selected}"
    pickle.dump(rf_cv, open(file_loc, "wb"))

# Use Date
date = "2019-11-02"
day_data = test_data.loc[date]

pred_for_tomorrow = pd.DataFrame({"Date": [], "company": [], "prediction": []})

# Predict each stock using the 2nd January Data
for cluster_selected in clusters_df.Cluster.unique():
    rf_cv = pickle.load(
        open(os.getcwd() + f"\\Pickle_Files\\Cluster_{cluster_selected}", "rb")
    )
    best_rf = rf_cv.best_estimator_
    cluster_data = day_data.loc[
        day_data.symbol.isin(
            clusters_df.loc[
                clusters_df.Cluster == cluster_selected, "Companies"
            ].tolist()
        )
    ].copy()
    cluster_data = cluster_data.dropna()
    if cluster_data.shape[0] > 0:
        X_test = cluster_data.loc[:, Target_variables]

        pred_for_tomorrow = pred_for_tomorrow.append(
            pd.DataFrame(
                {
                    "Date": cluster_data.index,
                    "company": cluster_data["symbol"],
                    "prediction": best_rf.predict_proba(X_test)[:, 1],
                }
            ),
            ignore_index=True,
        )

top_10_pred = pred_for_tomorrow.sort_values(by=["prediction"], ascending=False).head(10)
print(top_10_pred)

for selected_company in top_10_pred["company"]:
    actual = all_data[all_data.symbol == selected_company].loc[date, "Target_Direction"]
    pct_change = all_data[all_data.symbol == selected_company].loc[date, "Target"]
    top_10_pred.loc[top_10_pred.company == selected_company, "actual"] = actual
    top_10_pred.loc[top_10_pred.company == selected_company, "pct_change"] = pct_change
print(top_10_pred)