# linear algebra
import numpy as np     

# data processing
import pandas as pd    
#display max columns and rows, to be able to inspect big data tables
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)

#data visualization
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()

# conntects to google drive if ran in colab
# ignored if ran outside of colab

if 'google.colab' in str(get_ipython()):
  from google.colab import drive
  drive.mount('/content/drive')

# load file from google drive path, if ran in colab 
# load file from notebook folder, if ran outside colab

if 'google.colab' in str(get_ipython()):
    path = '/content/drive/My Drive/Colab Notebooks/stock_pricing/'
else:
    path = str()
file = path + 'large_Financials_unlinked.xlsm'

company_data = pd.read_excel(file)

#original data structure
company_data.head(5)

#original data shape
original_data_shape = company_data.shape
print('The original data shape is', original_data_shape)

#drop duplicates
company_data = company_data.drop_duplicates('Name' , keep='first')
print(original_data_shape[0] - company_data.shape[0], 'companies were removed due to duplication.')
print('The new data shape is', company_data.shape)

#remove not meaningful samples
current_num_samples = company_data.shape[0]
company_data = company_data[(company_data.EV > 0)]
company_data = company_data[(company_data.Market_Cap > 0)]
company_data = company_data[(company_data.EBITDA > 0)]
company_data = company_data[(company_data.Sector != 'Financials')]
new_num_samples = company_data.shape[0]
print(current_num_samples - new_num_samples, 'companies were removed due to not meaningful data.')
print('The new data shape is', company_data.shape)

#filter data for required columns
data = company_data[['Sector', 
                        'Name', 
                        'EV_TTM_EBITDA', 
                        'EBITDA_Forecast_CAGR_10y',
                        'Avg_EBITDA_Margin_Forecast_10y', 
                        'Beta_5y', 
                        'Net_Debt_perc_Market_Cap'
                       ]].set_index('Name')
data.head(5)

#check column data types
data.info()

#ensure columns have the correct data type
data['EV_TTM_EBITDA'] = data['EV_TTM_EBITDA'].apply(pd.to_numeric, errors='coerce')
data['EBITDA_Forecast_CAGR_10y'] = data['EBITDA_Forecast_CAGR_10y'].apply(pd.to_numeric, errors='coerce')
data['Avg_EBITDA_Margin_Forecast_10y'] = data['Avg_EBITDA_Margin_Forecast_10y'].apply(pd.to_numeric, errors='coerce')
data['Beta_5y'] = data['Beta_5y'].apply(pd.to_numeric, errors='coerce')
data['Net_Debt_perc_Market_Cap'] = data['Net_Debt_perc_Market_Cap'].apply(pd.to_numeric, errors='coerce')

data.info()

#check for missing values
data.isna().sum()

#remove missing values
current_num_samples = data.shape[0]
data = data.dropna()
data.isna().sum()
new_num_samples = data.shape[0]
print(current_num_samples - new_num_samples, 'companies were removed due to missing data.')
print('The new data shape is', data.shape)

#summary statistics before outlier removal
data.describe()

sns.histplot(data[['EV_TTM_EBITDA']])

sns.histplot(np.log(data[['EV_TTM_EBITDA']]))

sns.histplot(data[['EBITDA_Forecast_CAGR_10y']])

sns.histplot(data[['Avg_EBITDA_Margin_Forecast_10y']])

sns.histplot(data[['Beta_5y']])

sns.histplot(data[['Net_Debt_perc_Market_Cap']])

#check for outliers with multivariate Guassian
from sklearn.covariance import EllipticEnvelope

X = data[['EV_TTM_EBITDA', 
             'EBITDA_Forecast_CAGR_10y',
             'Avg_EBITDA_Margin_Forecast_10y',
             'Beta_5y',
             'Net_Debt_perc_Market_Cap']]
cov = EllipticEnvelope(random_state=0, contamination=0.1).fit(X)
# predict returns 1 for an inlier and -1 for an outlier
data[cov.predict(X)== -1]

#outlier removal
current_num_samples = data.shape[0]
data_new = data[cov.predict(X)== 1]
new_num_samples = data_new.shape[0]
print(current_num_samples - new_num_samples, 'companies were removed due to outliers.')
print('The new data shape is', data_new.shape)

#summary statistics after outlier removal
data_new.describe()

sns.histplot(data_new[['EV_TTM_EBITDA']])

sns.histplot(np.log(data_new[['EV_TTM_EBITDA']]))

sns.histplot(data_new[['EBITDA_Forecast_CAGR_10y']])

sns.histplot(data[['Avg_EBITDA_Margin_Forecast_10y']])

sns.histplot(data_new[['Beta_5y']])

sns.histplot(data_new[['Net_Debt_perc_Market_Cap']])

# split data into train and test splits
from sklearn.model_selection import train_test_split

X = data_new[['Sector', 'EBITDA_Forecast_CAGR_10y','Avg_EBITDA_Margin_Forecast_10y', 'Beta_5y', 'Net_Debt_perc_Market_Cap']];
Y = data_new[['EV_TTM_EBITDA']];
Y['Log_EV_TTM_EBITDA'] = np.log(Y['EV_TTM_EBITDA']);

X_train, X_test, Y_train, Y_test = train_test_split(X, Y, shuffle=True, test_size = 0.2, random_state=1);

# check training data shape
print('X_train has shape', X_train.shape)
print('y_train has shape', Y_train.shape)

#check X_train
X_train.head(5)

#check y_train
Y_train.head(5)

train_data = Y_train.merge(X_train, how='outer', left_index=True, right_index=True)
test_data = Y_test.merge(X_test, how='outer', left_index=True, right_index=True)

train_data.describe().round(2)

# correlation matrix for outcome variable

from string import ascii_letters
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

sns.set_theme(style="white")

# Compute the correlation matrix
corr = train_data.drop(columns="Log_EV_TTM_EBITDA").corr()

# Generate a mask for the upper triangle
mask = np.triu(np.ones_like(corr, dtype=bool))

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(11, 9))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(20, 230, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr, mask=mask, annot= True, cmap=cmap, vmax=.3, center=0,
            square=True, linewidths=.5, cbar_kws={"shrink": .5})

# correlation matrix for log of outcome variable

# Compute the correlation matrix
corr = train_data.drop(columns="EV_TTM_EBITDA").corr()

# Generate a mask for the upper triangle
mask = np.triu(np.ones_like(corr, dtype=bool))

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(11, 9))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(20, 230, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr, mask=mask, annot= True, cmap=cmap, vmax=.3, center=0,
            square=True, linewidths=.5, cbar_kws={"shrink": .5})

sns.set_theme()
sns.jointplot(data=train_data, x="EBITDA_Forecast_CAGR_10y", y="EV_TTM_EBITDA")
sns.jointplot(data=train_data, x="EBITDA_Forecast_CAGR_10y", y="Log_EV_TTM_EBITDA")

sns.jointplot(data=train_data, x="Avg_EBITDA_Margin_Forecast_10y", y="EV_TTM_EBITDA")
sns.jointplot(data=train_data, x="Avg_EBITDA_Margin_Forecast_10y", y="Log_EV_TTM_EBITDA")

sns.jointplot(data=train_data, x="Beta_5y", y="EV_TTM_EBITDA")
sns.jointplot(data=train_data, x="Beta_5y", y="Log_EV_TTM_EBITDA")

sns.jointplot(data=train_data, x="Net_Debt_perc_Market_Cap", y="EV_TTM_EBITDA")
sns.jointplot(data=train_data, x="Net_Debt_perc_Market_Cap", y="Log_EV_TTM_EBITDA")

#ML dependancies

#data preprocessing pipeline to preprocess data table
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import StandardScaler

#Grid search cross validation to find best hyper parameters
from sklearn.model_selection import GridSearchCV

#models considered
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.neighbors import KNeighborsRegressor
from sklearn.svm import SVR
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import VotingRegressor

#evaluation
from sklearn.model_selection import cross_val_score

#split data back into X and y
X_train = train_data[['EBITDA_Forecast_CAGR_10y','Avg_EBITDA_Margin_Forecast_10y', 'Beta_5y', 'Net_Debt_perc_Market_Cap', 'Sector']]
y_train = train_data['Log_EV_TTM_EBITDA']

X_test = test_data[['EBITDA_Forecast_CAGR_10y','Avg_EBITDA_Margin_Forecast_10y', 'Beta_5y', 'Net_Debt_perc_Market_Cap', 'Sector']]
y_test = test_data['Log_EV_TTM_EBITDA']

# check training and test data shape
print('X_train has shape', X_train.shape)
print('y_train has shape', y_train.shape)
print('X_test has shape', X_test.shape)
print('y_test has shape', y_test.shape)

# define categorical columns for transformation
categorical_variables = ['Sector']

# define numerical columns for transformation
numerical_variables = list(X_train.drop(columns=categorical_variables).columns)

# define column transformation for data preprocessing pipeline
column_transformer = ColumnTransformer(transformers=[
                                                 ('categorical_columns', OneHotEncoder(handle_unknown='ignore'), categorical_variables),
                                                 ('numeric_columns', StandardScaler(), numerical_variables)])

#ridge regression

grid_values = {'model__alpha': [0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 200, 300, 400, 500]}

model_1 = Pipeline([('column_transformer', column_transformer),
                    ('model', Ridge(max_iter = 100000))])

# default metric to optimize over grid parameters: accuracy
grid_model_1 = GridSearchCV(model_1, param_grid = grid_values, cv=5, scoring = 'r2')
grid_model_1 = grid_model_1.fit(X_train, y_train)

print('Grid best parameter: ', grid_model_1.best_params_)
print('Grid best score: ', grid_model_1.best_score_.round(3))

#lasso regression

grid_values = {'model__alpha': [0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 200, 300, 400, 500]}

model_2 = Pipeline([('column_transformer', column_transformer),
                    ('model', Lasso(max_iter = 100000))])

# default metric to optimize over grid parameters: accuracy
grid_model_2 = GridSearchCV(model_2, param_grid = grid_values, cv=5, scoring = 'r2')
grid_model_2 = grid_model_2.fit(X_train, y_train)

print('Grid best parameter: ', grid_model_2.best_params_)
print('Grid best score: ', grid_model_2.best_score_.round(3))

#KNeighborsRegressor

grid_values = {'model__n_neighbors': [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25]}

model_3 = Pipeline([('column_transformer', column_transformer),
                    ('model', KNeighborsRegressor())])

# default metric to optimize over grid parameters: accuracy
grid_model_3 = GridSearchCV(model_3, param_grid = grid_values, cv=5, scoring = 'r2')
grid_model_3 = grid_model_3.fit(X_train, y_train)

print('Grid best parameter: ', grid_model_3.best_params_)
print('Grid best score: ', grid_model_3.best_score_.round(3))

#SVR linear

grid_values = {'model__C': [0.01, 0.1, 1, 10, 100]}

model_4 = Pipeline([('column_transformer', column_transformer),
                    ('model', SVR(kernel='linear'))])

# default metric to optimize over grid parameters: accuracy
grid_model_4 = GridSearchCV(model_4, param_grid = grid_values, cv=5, scoring = 'r2')
grid_model_4 = grid_model_4.fit(X_train, y_train)

print('Grid best parameter: ', grid_model_4.best_params_)
print('Grid best score: ', grid_model_4.best_score_.round(3))

#SVR rbf

grid_values = {'model__gamma': [0.01, 0.1, 1, 10, 100], 
               'model__C': [0.01, 0.1, 1, 10, 100]}

model_5 = Pipeline([('column_transformer', column_transformer),
                    ('model', SVR(kernel='rbf'))])

# default metric to optimize over grid parameters: accuracy
grid_model_5 = GridSearchCV(model_5, param_grid = grid_values, cv=5, scoring = 'r2')
grid_model_5 = grid_model_5.fit(X_train, y_train)

print('Grid best parameter: ', grid_model_5.best_params_)
print('Grid best score: ', grid_model_5.best_score_.round(3))

#DecisionTreeRegressor

grid_values = {'model__max_depth': [2, 3, 4, 5, 6, 7, 8]}

model_6 = Pipeline([('column_transformer', column_transformer),
                    ('model', DecisionTreeRegressor())])

# default metric to optimize over grid parameters: accuracy
grid_model_6 = GridSearchCV(model_6, param_grid = grid_values, cv=5, scoring = 'r2')
grid_model_6 = grid_model_6.fit(X_train, y_train)

print('Grid best parameter: ', grid_model_6.best_params_)
print('Grid best score: ', grid_model_6.best_score_.round(3))

#RandomForestRegressor

grid_values = {'model__max_depth': [2, 3, 4, 5, 6, 7, 8, 9, 10], 
               'model__n_estimators': [50, 75, 100, 125, 150, 200, 250]}

model_7 = Pipeline([('column_transformer', column_transformer),
                    ('model', RandomForestRegressor())])

# default metric to optimize over grid parameters: accuracy
grid_model_7 = GridSearchCV(model_7, param_grid = grid_values, cv=5, scoring = 'r2')
grid_model_7 = grid_model_7.fit(X_train, y_train)

print('Grid best parameter: ', grid_model_7.best_params_)
print('Grid best score: ', grid_model_7.best_score_.round(3))

#VotingRegressor

model_8 = VotingRegressor([('ridge', grid_model_1), ('lasso', grid_model_2), ('knn', grid_model_3), ('svr_linear', grid_model_4), 
                                     ('svr_rbf', grid_model_5), ('tree', grid_model_6), ('rf', grid_model_7)])
model_8 = model_8.fit(X_train, y_train)
cv_score = cross_val_score(model_8, X_train, y_train, cv=5, scoring = 'r2').mean().round(3)
print('CV score', cv_score)

#model performance summary on cross validation folds

models = ['Ridge', 
          'Lasso',
          'KNeighborsRegressor',
          'SVR linear',
          'SVR rbf',
          'Decision Tree',
          'RandomForestRegressor', 
          'VotingRegressor']

scores = [grid_model_1.best_score_.round(3),
          grid_model_2.best_score_.round(3),
          grid_model_3.best_score_.round(3),
          grid_model_4.best_score_.round(3),
          grid_model_5.best_score_.round(3),
          grid_model_6.best_score_.round(3),
          grid_model_7.best_score_.round(3), 
          cv_score.round(3)]

model_scores = pd.DataFrame(data= scores, columns = ['CV_r2'], index = models)
print('The average cross validation performance of all trained models was as below')
model_scores

#get the best model as per cross validation performance
print('The model selected with the best Cross Validation performance was:', model_scores.index[np.argmax(model_scores['CV_r2'])])
model_scores.iloc[[np.argmax(model_scores['CV_r2'])]]

#model dictionary
model_dict = {'Ridge': grid_model_1,
              'Lasso': grid_model_2,
              'KNeighborsRegressor': grid_model_3,
              'SVR linear': grid_model_4,
              'SVR rbf': grid_model_5,
              'Decision Tree': grid_model_6,
              'RandomForestRegressor': grid_model_7, 
              'VotingRegressor': model_8
              }

#test set performance best models
print('The best model reaches below out of sample test set performance:')
model_dict[model_scores.index[np.argmax(model_scores['CV_r2'])]].score(X_test, y_test).round(3)