<jupyter_start><jupyter_text>Predicting Interest Rates for Bank LoansIn this project, I use a dataset of loans from a bank to predict the interest rate an applicant can expect- first, I wrangle the data by cleansing it of missing values, useless variables, and leaky information- then, I feature engineer a few columns to hopefully improve my model's performance- finally, I iterate through various models, optimizing hyperparameters for each, to find the best model- I obtain a model which can account for 90% of the variance in the data and predict interest rate with a meanerror of 0.9%<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import shap
import seaborn as sns
import time

from pandas_profiling import ProfileReport

from sklearn.pipeline import make_pipeline
from category_encoders import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score, validation_curve # k-fold CV

from sklearn.linear_model import LinearRegression, Ridge
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor

from sklearn.metrics import mean_absolute_error
from sklearn.metrics import r2_score<jupyter_output><empty_output><jupyter_text>Data Wrangling and Feature EngineeringHere, I clean my data into a format that will be a suitable input for modeling. I get rid of columns with high cardinality or a cardinality of 1, excessive NaN values, as well as all the columns containing leaky information. The leaky variables were variables in the dataset available after the loan had been approved - had I not gotten rid of them, my model would have had much greater predictive power (I tested this for myself but running the same analyses without a fully wrangled dataframe - that analysis is not included here) but it would have been a *useless* model since the majority of the predictive power would have come from columns like the 'grade' and 'subgrade' columns describing how likely the bank thought the applicant was to pay back the loan, which is very tightly correlated with the interest rate. The more risky a loan is for the bank, the higher interest they will charge on it.Although it is not obvious from the code, I found the documentation to this dataset and meticulously examined each column to determine if it should be kept, removed, or changed to a more suitable format. I engineer a few features, several of which proved to have strong predictive power.<jupyter_code># my dataframe is quite large, starting at 150 columns and over 100,000 rows. 
# I will do some cleaning and feature engineering, but first I need to see all my columns
pd.set_option('display.max_columns', 150)
pd.set_option('display.max_rows', 50)
df = pd.read_csv('/Users/nicholashagemann/Lambda/Datasets/loan_data.csv')
df.head()<jupyter_output><empty_output><jupyter_text>Note: The following code block is the final result of much time and computing power spent analyzing the dataframe - I have removed the code it took to create the wrangle function as its length would only serve to clutter my notebook.<jupyter_code># each row in the data is a specific loan that was issued by a bank
# I choose the issue date of the loan to be the index
df = pd.read_csv('/Users/nicholashagemann/Lambda/Datasets/loan_data.csv', 
                 parse_dates = ['issue_d', 'earliest_cr_line'], 
                 index_col = 'issue_d').sort_values('issue_d')

# we will need to wrangle this data quite a bit before it's suitable for predicting interest rate
def wrangle(df):
    # Begin Data Cleaning
    
    # first I get rid of all columns with more than 20% NaN values
    # note that the necessary columns will be dropped at the same time at the end of wrangle()
    nan_cols = []
    for col in df.columns:
        if df[col].isna().sum()/len(df) > 0.2:
            nan_cols.append(col)

    # each of the following columns had data that was only available after the loan had been approved
    # we want to predict the interest rate for a loan before the bank decides to approve the loan
    # thus, these columns contain leakage and should be dropped
    leaky_cols = ['out_prncp', 'out_prncp_inv', 'total_pymnt', 'total_pymnt_inv', 'sub_grade',
                  'total_rec_int', 'total_rec_prncp', 'total_rec_late_fee', 'recoveries',
                  'collection_recovery_fee', 'last_pymnt_d', 'last_pymnt_amnt', 'initial_list_status',
                  'last_fico_range_high', 'last_fico_range_low', 'collections_12_mths_ex_med',
                  'funded_amnt', 'funded_amnt_inv', 'debt_settlement_flag', 'grade', 'loan_status']

    # the following columns did not contain any useful data or were duplicates of other columns
    cols_to_drop = ['title', 'zip_code', 'addr_state', 'fico_range_low',
                    'Unnamed: 0', 'emp_title', 'url', 'id', 'last_credit_pull_d',
                    'earliest_cr_line']
    
    # if a column only has one value it has no predictive power
    for col in df.columns:
        if df[col].nunique() == 1:
            cols_to_drop.append(col)

    # Begin Feature Engineering

    # first I create a feature which identifies how old the borrower's credit is in days
    df['credit_age'] = df.index - df['earliest_cr_line']
    df['credit_age'].astype('timedelta64[D]').astype(int)
    df['credit_age'] = df['credit_age'].dt.days.astype(int)

    # we want our interest rate to be a float, not a string
    # do the same for the amount of credit the borrower is using
    df['int_rate'] = df['int_rate'].str.replace('%', '').astype(float)
    df['revol_util'] = df['revol_util'].str.replace('%', '').astype(float)

    # we want to represent length of the loan as an integer
    df.loc[df['term'] == ' 36 months', 'term'] = 36
    df.loc[df['term'] == ' 60 months', 'term'] = 60
    df['term'] = df['term'].astype(int)

    # we want to represent employment length as an integer
    df.loc[df['emp_length'] == '< 1 year', 'emp_length'] = 0
    df['emp_length'] = df['emp_length'].str.replace(' ', '').str.replace('year', '').str.replace('s', '').str.replace('+', '')
    df['emp_length'] = pd.to_numeric(df['emp_length'], downcast = 'integer', errors = 'ignore')

    # we want to represent whether the application was individual or joint as binary
    df.loc[df['application_type'] == 'Individual', 'application_type'] = 1
    df.loc[df['application_type'] == 'Joint App', 'application_type'] = 0
    df['application_type'] = df['application_type'].astype(int)

    # drop columns from the data cleaning
    df.drop(columns = nan_cols + leaky_cols + cols_to_drop, inplace = True)

    return df

df = wrangle(df)
# check how many unique values we have for each of our object columns
# looks like it's small enough to use OneHotEncoder
df.select_dtypes('object').nunique()
# check shape and head of dataframe
print('Shape is: ', df.shape)
df.head()
df.describe()<jupyter_output><empty_output><jupyter_text>Creating Training and Validation DataHere, I create the X (predictor) and y (target) matrices which will then be split into training and testing sets for my models.<jupyter_code># we want to have ~80% of our data for training and ~20% for validation
# we see that the 80th percentile is December 2016
# we have a very large dataset, so we choose to start in January 2017
display(df.iloc[0:1], df.iloc[round(df.shape[0]*0.8):round(df.shape[0]*0.8)+1])
# our prediction target is the interest rate, our predictive matrix will be everything except the target
target = 'int_rate'
y = df[target]
X = df.drop(columns = target)
# separate our training and validation data by our cutoff of January 2017
cutoff = '2017-01-01'
mask = X.index < cutoff
X_train, y_train = X[mask], y[mask]
X_val, y_val = X[~mask], y[~mask]

# check shapes to make sure dimensions are correct
print(f'The shape of the training data is {X_train.shape} for X and {y_train.shape} for y.')
print(f'The shape of the validation data is {X_val.shape} for X and {y_val.shape} for y.')

# make sure no data was lost in the splice
assert len(X_train) + len(X_val) == len(X)
# our baseline mean absolute error will be the mean interest rate
# baseline R^2 score will of course be 0 or approximately 0
baseline_reg = [y_train.mean()]*len(y_val)
print('Baseline MAE:', mean_absolute_error(y_val, baseline_reg))
print('Baseline R2:', r2_score(y_val, baseline_reg))<jupyter_output>Baseline MAE: 3.360100268400727
Baseline R2: -0.009930713374478373<jupyter_text>Best Model - XGBoost For the sake of the reader's time, I include my final model here so that it is not necessary to scroll through the various models and tuning of hyperparameters.My final model predicts interest rate with a mean error below 0.9% and a cross-validated R^2 of 0.91.<jupyter_code>%%time
model_final = make_pipeline(OneHotEncoder(use_cat_names = True),
                           StandardScaler(),
                           SimpleImputer(strategy='mean'),
                           XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        reg_alpha = 0,
                                        reg_lambda = 1,
                                        gamma = 0))
model_final.fit(X_train, y_train)
print('Training MAE:', mean_absolute_error(y_train, model_final.predict(X_train)))
print('Validation MAE:', mean_absolute_error(y_val, model_final.predict(X_val)))
print('Validation R2:', model_boost.score(X_val, y_val))
cross_vals = cross_val_score(model_final, X, y)
print(f'Scores: {cross_vals}\nMean: {np.mean(cross_vals)}')<jupyter_output>Scores: [0.87971062 0.92950602 0.91982579 0.93255404 0.90019554]
Mean: 0.9123584023146762<jupyter_text>On Optimizing HyperparametersFor each model, I will optimize whatever hyperparameters are available for that model.In testing the hyperparameters, I will look at the Mean Absolute Error and R^2 (what percentage of the variance in the data my model accounts for) for the training and validation data, as well as the time it took - the time is necessary since when testing each hyperparameter I do not necessarily want the previously tested hyperparameters to be optimal since the model may take too much time.By having the MAE, R^2, and Time for each hyperparameter I can know which hyperparameter I will use in my final model while also learning which hyperparameter will be suitable for tuning other hyperparameters - i.e. for my final model I will use the hyperparameter with the lowest MAE and highest R^2, but when tuning the other hyperparameters I will use a value which has a MAE and R^2 *similar* to the best but not as time-expensive.I choose not to use a GridSearch for the hyperparameters since I am attempting to optimize every hyperparameter and the exponential increase in the time it would take to train my model is not worth the small amount of additional certainty I would gain from knowing I've found the absolute best combination - also, I am sure from past experience with exactly this problem that the increase in accuracy from a GridSearch would be very small compared to the overall increase in accuracy gained through my method (which itself is still small compared to the accuracy of the model itself). In going through each model, I also choose the order in which I optimize hyperparameters very deliberately. I start with hyperparameters that will have the greatest effect on the accuracy and then move to hyperparameters with smaller effects - in essence, ordering them so that the later hyperparameters do not affect the earlier hyperparameters.   Decision Tree RegressorHere, I will test the decision tree model and optimize hyperparameters. It is very unlikely that a decision tree will be the best model, as we naturally expect any random forest to outcompete a single tree.We find, after optimizing hyperparameters, the MAE is 2.28% and the R^2 is 0.4The only significant hyperparameter to optimize for a decision tree is the depth of the tree - since we have a large dataset, we can rely on the algorithm and not have to choose the minimum samples for a split or a leaf, the maximum features or leaves, or the minimum decrease in impurity. We see the depth of the tree significantly affects the accuracy and R^2 but the split criteria does not, and also that it is wise to stick with the default "best" splitter and not the "random" splitter.<jupyter_code># do a broad search first then a narrow search next to save time
depths = [1] + [i for i in range(4, 21, 4)]
train_acc_dt = []
val_acc_dt = []
train_r2_dt = []
val_r2_dt = []
times = []

# for each depth, train the model on the training data then add the accuracies, r^2, and time
for depth in depths:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         DecisionTreeRegressor(random_state = 42,
                                               max_depth = depth))
    model.fit(X_train, y_train)
    train_acc_dt.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_dt.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_dt.append(model.score(X_train, y_train))
    val_r2_dt.append(model.score(X_val, y_val))
    times.append(time.time() - start)

# display stored values to determine best parameter
pd.DataFrame(list(zip(train_acc_dt, val_acc_dt, train_r2_dt, val_r2_dt, times)), index = depths,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
# time for the narrow search
depths = range(8, 17)
train_acc_dt = []
val_acc_dt = []
train_r2_dt = []
val_r2_dt = []
times = []

# for each depth, train the model on the training data then add the accuracies, r^2, and time
for depth in depths:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         DecisionTreeRegressor(random_state = 42,
                                               max_depth = depth))

    model.fit(X_train, y_train)

    train_acc_dt.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_dt.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_dt.append(model.score(X_train, y_train))
    val_r2_dt.append(model.score(X_val, y_val))
    times.append(time.time() - start)

# display stored values to determine best parameter
pd.DataFrame(list(zip(train_acc_dt, val_acc_dt, train_r2_dt, val_r2_dt, times)), index = depths,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
# Here we decide if we should be using MSE or Friedman_MSE to split the cells in the tree - I ignore MAE
# and reducing poisson deviance because the models take significantly longer to train than the MSE variants
# and in attempting the model, my kernel crashes and must be restarted from the beginning :(
criteria = ["mse", "friedman_mse"]
train_acc_dt = []
val_acc_dt = []
train_r2_dt = []
val_r2_dt = []

for crit in criteria:
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         DecisionTreeRegressor(random_state = 42,
                                               criterion = crit,
                                               max_depth = 12))

    model.fit(X_train, y_train)

    train_acc_dt.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_dt.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_dt.append(model.score(X_train, y_train))
    val_r2_dt.append(model.score(X_val, y_val))

pd.DataFrame(list(zip(train_acc_dt, val_acc_dt, train_r2_dt, val_r2_dt)), index = criteria,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2'])
# was curious how the 'random' splitter would compare to the 'best' - as expected, 'best' is much better
splits = ['best', 'random']
train_acc_dt = []
val_acc_dt = []
train_r2_dt = []
val_r2_dt = []
times = []

for split in splits:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         DecisionTreeRegressor(random_state = 42,
                                               splitter = split,
                                               max_depth = 12))
    
    model.fit(X_train, y_train)

    train_acc_dt.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_dt.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_dt.append(model.score(X_train, y_train))
    val_r2_dt.append(model.score(X_val, y_val))
    times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_dt, val_acc_dt, train_r2_dt, val_r2_dt, times)), index = splits,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
# make our decision tree model here
model_dt = make_pipeline(OneHotEncoder(use_cat_names = True),
                        StandardScaler(),
                        SimpleImputer(strategy = 'mean'),
                        DecisionTreeRegressor(random_state = 42,
                                              max_depth = 12))
%%time
model_dt.fit(X_train, y_train)
print('Training MAE:', mean_absolute_error(y_train, model_dt.predict(X_train)))
print('Validation MAE:', mean_absolute_error(y_val, model_dt.predict(X_val)))
print('Validation R2:', model_dt.score(X_val,y_val))

# our MAE isn't bad for a weak model, but we can do better. R^2 is very low. 
# here I graph the most important features in our model

importances = model_dt.named_steps['decisiontreeregressor'].feature_importances_
features = model_dt.named_steps['onehotencoder'].get_feature_names()
feat_imp_dt = pd.Series(importances, index = features,
                       name = 'tree').abs().sort_values(ascending=False)
feat_imp_dt.head(10).plot(kind = 'barh', title = 'Most Important Features')<jupyter_output><empty_output><jupyter_text>Making Tuning DatasetFor the random forest model, my kernel kept crashing when optimizing hyperparameters by using the full dataset so I decide to tune the hyperparameters using a 'tuning' dataset which is 1/10 of the original dataset<jupyter_code># make tune data which is 1/10 of the train data
X_tune_train = X_train[:len(X_train)//10]
y_tune_train = y_train[:len(y_train)//10]
X_tune_val = X_val[:len(X_val)//10]
y_tune_val = y_val[:len(y_val)//10]

# check to make sure our tuning datasets are indeed 1/10 of our previous training sets
print(f'{X_tune_train.shape}, {X_train.shape}\n{y_tune_train.shape}, {y_train.shape}')
print(f'{X_tune_val.shape}, {X_val.shape}\n{y_tune_val.shape}, {y_val.shape}')<jupyter_output>(10364, 58), (103649, 58)
(10364,), (103649,)
(2468, 58), (24685, 58)
(2468,), (24685,)<jupyter_text>Random Forest Regressor<jupyter_code># first hyperparameter to optimize in a random forest is the number of trees
# obviously more trees will only increase the accuracy so we will instead examine when the accuracy approaches
# an asymptote - we will also learn how many estimators we should use for further tuning of hyperparameters, as
# we want enough so that the accuracy is similar to the best accuracy but also not too many for time constraints

estimators = range(10, 201, 10)
train_acc_rf = []
val_acc_rf = []
train_r2_rf = []
val_r2_rf = []
seconds = []

# we will use the previously determined max_depth hyperparameter
for estimator in estimators:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                        StandardScaler(),
                        SimpleImputer(strategy = 'mean'),
                        RandomForestRegressor(random_state = 42,
                                              n_jobs = -1,
                                              max_depth = 12,
                                              n_estimators = estimator))
    
    model.fit(X_tune_train, y_tune_train)
    
    train_acc_rf.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_rf.append(mean_absolute_error(y_train, model.predict(X_train)))
    train_r2_rf.append(model.score(X_train, y_train))
    val_r2_rf.append(model.score(X_val, y_val))
    seconds.append(time.time() - start)
    
pd.DataFrame(list(zip(train_acc_rf, val_acc_rf, train_r2_rf, val_r2_rf, seconds)), index = estimators,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
# determine what percentage of the data to sample for each decision tree
samples = np.arange(0.05, 1, 0.05)
train_acc_rf = []
val_acc_rf = []
train_r2_rf = []
val_r2_rf = []
seconds = []

for sample in samples:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                        StandardScaler(),
                        SimpleImputer(strategy = 'mean'),
                        RandomForestRegressor(random_state = 42,
                                              n_jobs = -1,
                                              max_depth = 12,
                                              n_estimators = 30,
                                              max_samples = sample))
    model.fit(X_tune_train, y_tune_train)
    
    train_acc_rf.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_rf.append(mean_absolute_error(y_train, model.predict(X_train)))
    train_r2_rf.append(model.score(X_train, y_train))
    val_r2_rf.append(model.score(X_val, y_val))
    seconds.append(time.time() - start)
    
pd.DataFrame(list(zip(train_acc_rf, val_acc_rf, train_r2_rf, val_r2_rf, seconds)), index = samples,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
model_rf = make_pipeline(OneHotEncoder(use_cat_names = True),
                        StandardScaler(),
                        SimpleImputer(strategy = 'mean'),
                        RandomForestRegressor(random_state = 42,
                                              n_jobs = -1,
                                              max_depth = 12,
                                              n_estimators = 180,
                                              max_samples = 0.4))
%%time
model_rf.fit(X_train, y_train)
print('Training MAE:', mean_absolute_error(y_train, model_rf.predict(X_train)))
print('Validation MAE:', mean_absolute_error(y_val, model_rf.predict(X_val)))
print('Validation R2:', model_rf.score(X_val,y_val))

# interestingly, our random forest model is not significantly better than our decision tree
# the top 3 importances are the same as the decision tree
importances = model_rf.named_steps['randomforestregressor'].feature_importances_
features = model_rf.named_steps['onehotencoder'].get_feature_names()
feat_imp_rf = pd.Series(importances, index = features,
                       name = 'forest').abs().sort_values(ascending=False)
feat_imp_rf.head(10).plot(kind = 'barh')<jupyter_output><empty_output><jupyter_text>Boost ModelAfter testing out decision trees, random forests, boosted random forests, linear regression, and ridge regression, I found that the boosted model gave the lowest MAE and highest R^2 by a significant amount. Given this information, I decided to retrain every hyperparameter instead of assuming that overlapping parameters with other models like max_depth would have the same optimized value.<jupyter_code># finding the optimal depth
depths = range(1, 20, 3)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
times = []

for depth in depths:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = depth))
    model.fit(X_train, y_train)

    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, times)), index = depths,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
# still finding the optimal depth
depths = range(7, 13)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
times = []

for depth in depths:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = depth))

    model.fit(X_train, y_train)

    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, times)), index = depths,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
# optimal depth is 9
# seeing how number of trees affects the model accuracy
estimators = [10] + [i for i in range(50, 300, 50)]
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
times = []

for estimator in estimators:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = estimator))

    model.fit(X_train, y_train)

    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, times)), index = estimators,
            columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 'Time'])
# 100 trees should be enough for optimizing hyperparameters, but use >200 for final model
# optimizing learning rate
learns = np.arange(0.01, 0.31, 0.05)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for learn in learns:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = learn))
    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = learns,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# need to expand learning rate parameter as accuracy kept increasing when learning rate hit 0.26
learns = np.arange(0.25, 0.51, 0.05)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for learn in learns:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = learn))
    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = learns,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# finding the true optimal learning rate
learns = np.arange(0.25, 0.35, 0.01)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for learn in learns:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = learn))
    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = learns,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# optimizing colsample_by tree hyperparameter, i.e. the subsample ratio for the columns when constructing
# the dataset for each tree
samples = np.arange(0.01, 1.00, 0.05)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for sample in samples:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = sample))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = samples,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# optimizing colsample_bylevel hyperparameter, i.e. the subsample ratio of the levels in each column
# when constructing each tree
samples = np.arange(0.01, 1.00, 0.05)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for sample in samples:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bylevel = sample))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = samples,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# optimizing colsample_bynode hyperparameter, i.e. the subsample ratio for each column when splitting
# a branch to make two new branches
samples = np.arange(0.01, 1.00, 0.05)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for sample in samples:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bynode = sample))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = samples,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# optimizing subsamples
subsamples = np.arange(0.01, 1.00, 0.05)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for sample in subsamples:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        subsample = sample))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = subsamples,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# we find that the best subsample value is the default value of 1
# learning if a different boosting method would provide greater accuracy
boosters = ['gbtree', 'gblinear', 'dart']
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for boost in boosters:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        booster = boost))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = boosters,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# optimizing alphas hyperparameter
alphas = range(0, 100, 5)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for alpha in alphas:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        reg_alpha = alpha))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = alphas,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# finish optimizing alphas hyperparameter
alphas = range(0, 5)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for alpha in alphas:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        reg_alpha = alpha))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_tune_train, model.predict(X_tune_train)))
    val_acc_boost.append(mean_absolute_error(y_tune_val, model.predict(X_tune_val)))
    train_r2_boost.append(model.score(X_tune_train, y_tune_train))
    val_r2_boost.append(model.score(X_tune_val, y_tune_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = alphas,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# optimizing lambda hyperparameter
lambdas = range(1, 6)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for l in lambdas:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        reg_lambda = l))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = lambdas,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])

# optimizing gamma hyperparameter
gammas = range(0, 6)
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for g in gammas:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = 'mean'),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        gamma = g))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = gammas,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])

# should we choose median instead of mean for SimpleImputer?
strats = ['mean', 'median']
train_acc_boost = []
val_acc_boost = []
train_r2_boost = []
val_r2_boost = []
train_times = []
score_times = []

for strat in strats:
    start = time.time()
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                         StandardScaler(),
                         SimpleImputer(strategy = strat),
                         XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        reg_alpha = 0,
                                        reg_lambda = 1,
                                        gamma = 0))

    model.fit(X_train, y_train)

    train_times.append(time.time() - start)
    train_acc_boost.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_boost.append(mean_absolute_error(y_val, model.predict(X_val)))
    train_r2_boost.append(model.score(X_train, y_train))
    val_r2_boost.append(model.score(X_val, y_val))
    score_times.append(time.time() - start)

pd.DataFrame(list(zip(train_acc_boost, val_acc_boost, train_r2_boost, val_r2_boost, train_times, score_times)), 
             index = strats,
             columns = ['Training Accuracy', 'Validation Accuracy', 'Training R^2', 'Validation R^2', 
                        'Train Time', 'Score Time'])
# mean works best, as expected - but median works surprisingly well
# time to build our final model with all our hyperparameters
model_boost = make_pipeline(OneHotEncoder(use_cat_names = True),
                           StandardScaler(),
                           SimpleImputer(strategy='mean'),
                           XGBRegressor(random_state = 42,
                                        n_jobs = -1,
                                        max_depth = 9,
                                        n_estimators = 100,
                                        learning_rate = 0.29,
                                        colsample_bytree = 1,
                                        reg_alpha = 0,
                                        reg_lambda = 1,
                                        gamma = 0))
%%time
model_boost.fit(X_train, y_train)
%%time
print('Training MAE:', mean_absolute_error(y_train, model_boost.predict(X_train)))
print('Validation MAE:', mean_absolute_error(y_val, model_boost.predict(X_val)))
print('Validation R2:', model_boost.score(X_val, y_val))
# check cross validation score
cross_vals = cross_val_score(model_boost, X, y)
print(np.mean(cross_vals))
# it's actually higher than our previous R^2
importances = model_boost.named_steps['xgbregressor'].feature_importances_
features = model_boost.named_steps['onehotencoder'].get_feature_names()
feat_imp_boost = pd.Series(importances, index = features,
                          name = 'boost').abs().sort_values(ascending=False)
feat_imp_boost.head(10).plot(kind = 'barh')<jupyter_output><empty_output><jupyter_text>Ridge Regression<jupyter_code>alphas = np.arange(0, 6)
train_acc_ridge = []
val_acc_ridge = []
for a in alphas:
    model = make_pipeline(OneHotEncoder(use_cat_names = True),
                            StandardScaler(),
                            SimpleImputer(strategy = 'mean'),
                            Ridge(alpha = a))
    model.fit(X_train, y_train)
    train_acc_ridge.append(mean_absolute_error(y_train, model.predict(X_train)))
    val_acc_ridge.append(mean_absolute_error(y_val, model.predict(X_val)))
    
pd.DataFrame(list(zip(train_acc_ridge, val_acc_ridge)), index = alphas, 
             columns = ['Training Accuracy', 'Validation Accuracy'])
model_ridge = make_pipeline(OneHotEncoder(use_cat_names = True),
                            StandardScaler(),
                            SimpleImputer(strategy = 'mean'),
                            Ridge(alpha = 1))
model_ridge.fit(X_train, y_train)
print('Training MAE:', mean_absolute_error(y_train, model_ridge.predict(X_train)))
print('Validation MAE:', mean_absolute_error(y_val, model_ridge.predict(X_val)))
print('Validation R2:', model_ridge.score(X_val,y_val))
importances = model_ridge.named_steps['ridge'].coef_
features = model_ridge.named_steps['onehotencoder'].get_feature_names()
feat_imp_ridge = pd.Series(importances, index = features,
                          name = 'ridge').abs().sort_values(ascending=False)
feat_imp_ridge.head(10).plot(kind = 'barh')<jupyter_output><empty_output><jupyter_text>Linear Regression<jupyter_code>model_linear = make_pipeline(OneHotEncoder(use_cat_names = True),
                            StandardScaler(),
                            SimpleImputer(strategy = 'mean'),
                            LinearRegression())
model_linear.fit(X_train, y_train)
print('Training MAE:', mean_absolute_error(y_train, model_linear.predict(X_train)))
print('Validation MAE:', mean_absolute_error(y_val, model_linear.predict(X_val)))
print('Validation R2:', model_linear.score(X_val,y_val))
importances = model_linear.named_steps['linearregression'].coef_
features = model_linear.named_steps['onehotencoder'].get_feature_names()
feat_imp_linear = pd.Series(importances, index = features,
                           name = 'linear').abs().sort_values(ascending=False)
feat_imp_linear.head(10).plot(kind = 'barh')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09548034805837638, 'fraction_numerical': 0.015576201606725463, 'mean_word_length': 2.568764895555867, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 130, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '12854', 'license': 'MIT', 'n_tokens_mistral': 15446, 'n_tokens_neox': 14570, 'n_words': 3850, 'path': 'InterestRatesProject-Final.ipynb', 'repo_name': 'nicholashagemann/InterestRatesProject', 'size': '50911'}