{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-08-19T19:17:47.23368Z","iopub.execute_input":"2021-08-19T19:17:47.234007Z","iopub.status.idle":"2021-08-19T19:17:47.384747Z","shell.execute_reply.started":"2021-08-19T19:17:47.233979Z","shell.execute_reply":"2021-08-19T19:17:47.383818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Import relevant libraries\nimport pandas as pd\nimport numpy as np\nimport os\nimport glob\nfrom functools import reduce\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport xgboost as xgb\nfrom sklearn.datasets import load_boston\nfrom sklearn.model_selection import train_test_split\nimport matplotlib.pyplot as plt\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import GridSearchCV\nimport keras\nfrom keras.callbacks import ModelCheckpoint\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Activation, Flatten\nfrom keras import backend as K\nimport tensorflow as tf\n\n# Import training data\ntrain = pd.read_csv('../input/optiver-realized-volatility-prediction/train.csv')\n\n# Grab Unique Array IDs\nstock_id_array = train.stock_id.unique()\n\n","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:52.269201Z","iopub.execute_input":"2021-08-19T19:17:52.269546Z","iopub.status.idle":"2021-08-19T19:17:57.017104Z","shell.execute_reply.started":"2021-08-19T19:17:52.269517Z","shell.execute_reply":"2021-08-19T19:17:57.015888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will now examine the structure of the data that we were provided.","metadata":{}},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:18:47.266825Z","iopub.execute_input":"2021-08-19T19:18:47.267184Z","iopub.status.idle":"2021-08-19T19:18:47.277227Z","shell.execute_reply.started":"2021-08-19T19:18:47.267148Z","shell.execute_reply":"2021-08-19T19:18:47.276575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.shape","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:18:49.353764Z","iopub.execute_input":"2021-08-19T19:18:49.35435Z","iopub.status.idle":"2021-08-19T19:18:49.35942Z","shell.execute_reply.started":"2021-08-19T19:18:49.354301Z","shell.execute_reply":"2021-08-19T19:18:49.358269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Take a sample from the book data\n\nsample = pd.read_parquet(f'../input/optiver-realized-volatility-prediction/book_train.parquet/stock_id=1')","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:19:23.221841Z","iopub.execute_input":"2021-08-19T19:19:23.222548Z","iopub.status.idle":"2021-08-19T19:19:24.129318Z","shell.execute_reply.started":"2021-08-19T19:19:23.222508Z","shell.execute_reply":"2021-08-19T19:19:24.128356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample.head()","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:19:34.866405Z","iopub.execute_input":"2021-08-19T19:19:34.866794Z","iopub.status.idle":"2021-08-19T19:19:34.881493Z","shell.execute_reply.started":"2021-08-19T19:19:34.866761Z","shell.execute_reply":"2021-08-19T19:19:34.880448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample.shape","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:19:49.216645Z","iopub.execute_input":"2021-08-19T19:19:49.216996Z","iopub.status.idle":"2021-08-19T19:19:49.223675Z","shell.execute_reply.started":"2021-08-19T19:19:49.216967Z","shell.execute_reply":"2021-08-19T19:19:49.222478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Feature Gathering**\nAside from the other features already given to us, we want to gather the 4 statistics that we can develop with our given data and that are most likely related to future realized volatility: Bid/Ask Spread (BAS), Weighted Average Price (WAP), Log Return, and Past Realized Volatility. Since BAS, WAP, and Log Returns are calculated upon each second in bucket, we will have to aggregate these so that there is only one BAS and WAP per time ID and stock ID pair. We will also need to aggregate the other features we are already given from the competition repository.","metadata":{}},{"cell_type":"code","source":"# Calculate Bid Ask Spread\n\n# BAS calculation for a stock\ndef bas_calculation_per_id(stock_id, data_type):\n    df_book_data = pd.read_parquet(f'../input/optiver-realized-volatility-prediction/book_{data_type}.parquet/stock_id={stock_id}')\n    df_book_data['bas'] = df_book_data[['ask_price1', 'ask_price2']].min(axis=1)/df_book_data[['bid_price1', 'bid_price2']].max(axis=1) - 1\n    df_book_data['stock_id'] = stock_id\n    return df_book_data\n\n# Loop through each stock\ndef bas_calculation(stock_id_array, data_type):\n    df_bas = pd.DataFrame()\n    for stock_id in stock_id_array:\n        df_bas_id = bas_calculation_per_id(stock_id, data_type).groupby(by = ['stock_id', 'time_id'], as_index = False)['bas'].mean()\n        df_bas = pd.concat([df_bas,df_bas_id])\n    return df_bas\n\ndf_bas = bas_calculation(stock_id_array, 'train')\nprint(df_bas)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:20:29.5752Z","iopub.execute_input":"2021-08-19T19:20:29.575535Z","iopub.status.idle":"2021-08-19T19:21:36.352095Z","shell.execute_reply.started":"2021-08-19T19:20:29.575507Z","shell.execute_reply":"2021-08-19T19:21:36.350562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate Weighted Average Price\n\n# WAP calculation for a stock\ndef wap_calculation_per_id(stock_id, data_type):\n    df_book_data = pd.read_parquet(f'../input/optiver-realized-volatility-prediction/book_{data_type}.parquet/stock_id={stock_id}')\n    df_book_data['stock_id'] = stock_id\n    df_book_data['wap'] = (df_book_data['bid_price1'] * df_book_data['ask_size1'] + df_book_data['ask_price1'] \n                           * df_book_data['bid_size1']) / (df_book_data['bid_size1'] + df_book_data['ask_size1'])           \n    return df_book_data[['stock_id', 'time_id', 'seconds_in_bucket', 'wap']]\n\n# Loop through each stock\ndef wap_calculation(stock_id_array, data_type):\n    df_wap = pd.DataFrame()\n    for stock_id in stock_id_array:\n        df_wap_id = wap_calculation_per_id(stock_id, data_type).groupby(by = ['stock_id', 'time_id'], as_index = False)['wap'].mean()\n        df_wap = pd.concat([df_wap,df_wap_id])\n    return df_wap\n\ndf_wap = wap_calculation(stock_id_array, 'train')\nprint(df_wap)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:21:36.354559Z","iopub.execute_input":"2021-08-19T19:21:36.354908Z","iopub.status.idle":"2021-08-19T19:22:18.03377Z","shell.execute_reply.started":"2021-08-19T19:21:36.354877Z","shell.execute_reply":"2021-08-19T19:22:18.032878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate Log Returns (LR)\n\n# LR calculation for a list of stock prices\ndef lr(list_stock_prices):\n    list_stock_prices['lr'] = np.log(list_stock_prices['wap']).diff()\n    return list_stock_prices\n\n# LR calculation for a stock\ndef lr_calculation_per_id(stock_id, data_type):\n    df_book_data = wap_calculation_per_id(stock_id, data_type)\n    df_lr_per_id = lr(df_book_data)\n    return df_lr_per_id[['stock_id', 'time_id', 'seconds_in_bucket', 'lr']]\n\n# Loop through each stock    \ndef lr_calculation(stock_id_array, data_type):\n    df_lr = pd.DataFrame()\n    for stock_id in stock_id_array:\n        df_lr_id = lr_calculation_per_id(stock_id, data_type).groupby(by = ['stock_id', 'time_id'], as_index = False)['lr'].mean()\n        df_lr = pd.concat([df_lr,df_lr_id])\n    return df_lr\n\ndf_lr = lr_calculation(stock_id_array, 'train')\nprint(df_lr)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:22:18.035872Z","iopub.execute_input":"2021-08-19T19:22:18.036265Z","iopub.status.idle":"2021-08-19T19:22:55.708592Z","shell.execute_reply.started":"2021-08-19T19:22:18.036224Z","shell.execute_reply":"2021-08-19T19:22:55.707499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate Past Realized Volatility (RV)\n\n# RV calculation for a series of LRs\ndef realized_volatility(series_log_return):\n    return np.sqrt(np.sum(series_log_return['lr']**2))\n    \n# Loop through each stock\ndef rv_calculation(stock_id_array, data_type):\n    df_rv = pd.DataFrame()\n    for stock_id in stock_id_array:\n        df_book_data = lr_calculation_per_id(stock_id, data_type)\n        df_rv_id = df_book_data.groupby(by = ['stock_id', 'time_id'], as_index = False).apply(realized_volatility)\n        df_rv_id.columns = ['stock_id', 'time_id', 'rv']\n        df_rv = pd.concat([df_rv,df_rv_id])\n    return df_rv\n\ndf_rv = rv_calculation(stock_id_array, 'train')\nprint(df_rv)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:22:55.710238Z","iopub.execute_input":"2021-08-19T19:22:55.710572Z","iopub.status.idle":"2021-08-19T19:25:52.825517Z","shell.execute_reply.started":"2021-08-19T19:22:55.710541Z","shell.execute_reply":"2021-08-19T19:25:52.824491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Grab and Aggregate other given Relevant Features\n\n# WAP calculation for a stock\ndef wap_calculation_per_id(stock_id, data_type):\n    df_book_data = pd.read_parquet(f'../input/optiver-realized-volatility-prediction/book_{data_type}.parquet/stock_id={stock_id}')\n    df_book_data['stock_id'] = stock_id\n    df_book_data['wap'] = (df_book_data['bid_price1'] * df_book_data['ask_size1'] + df_book_data['ask_price1'] \n                           * df_book_data['bid_size1']) / (df_book_data['bid_size1'] + df_book_data['ask_size1'])           \n    return df_book_data[['stock_id', 'time_id', 'seconds_in_bucket', 'wap']]\n\n# Loop through each stock\ndef wap_calculation(stock_id_array, data_type):\n    df_wap = pd.DataFrame()\n    for stock_id in stock_id_array:\n        df_wap_id = wap_calculation_per_id(stock_id, data_type).groupby(by = ['stock_id', 'time_id'], as_index = False)['wap'].mean()\n        df_wap = pd.concat([df_wap,df_wap_id])\n    return df_wap\n\ndf_wap = wap_calculation(stock_id_array, 'train')\nprint(df_wap)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:25:52.826933Z","iopub.execute_input":"2021-08-19T19:25:52.827218Z","iopub.status.idle":"2021-08-19T19:26:27.104532Z","shell.execute_reply.started":"2021-08-19T19:25:52.82719Z","shell.execute_reply":"2021-08-19T19:26:27.103569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Gather and Aggregate other given features\n\n# Grab features for each stock\ndef other_calculation_per_id(stock_id, data_type):\n    df_book_data = pd.read_parquet(f'../input/optiver-realized-volatility-prediction/book_{data_type}.parquet/stock_id={stock_id}')\n    df_book_data['stock_id'] = stock_id\n    df = df_book_data.groupby(by = ['stock_id', 'time_id'], \n                      as_index = False)[['bid_price1', 'ask_price1', 'bid_price2', \n                      'ask_price2', 'bid_size1', 'ask_size1', 'bid_size2', 'ask_size2']].mean()\n    df_two =  df_book_data.groupby(by = ['stock_id', 'time_id'],as_index = False)[['seconds_in_bucket']].max()\n    return pd.merge(df, df_two,  how='left', left_on=['stock_id','time_id'], right_on = ['stock_id','time_id'])\n\n# Loop through each stock\ndef other_calculation(stock_id_array, data_type):\n    df_other = pd.DataFrame()\n    for stock_id in stock_id_array:\n        df_other_id = other_calculation_per_id(stock_id, data_type)\n        df_other = pd.concat([df_other,df_other_id])\n    return df_other\n\ndf_other = other_calculation(stock_id_array, 'train')\nprint(df_other)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:26:27.106008Z","iopub.execute_input":"2021-08-19T19:26:27.106419Z","iopub.status.idle":"2021-08-19T19:27:19.1137Z","shell.execute_reply.started":"2021-08-19T19:26:27.10638Z","shell.execute_reply":"2021-08-19T19:27:19.112658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Concatonate each statistic into one dataframe\n\ndf_features = df_bas.merge(df_wap,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(df_lr,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(df_rv,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(df_other,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(train,on=['stock_id', 'time_id'])\n\n# Delete rows with NaN or missing values\ndf_features = df_features.dropna()\n\n# Make feature dataframe copy\ndf_features_copy = df_features\nprint(df_features)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:27:19.1151Z","iopub.execute_input":"2021-08-19T19:27:19.115448Z","iopub.status.idle":"2021-08-19T19:27:19.731513Z","shell.execute_reply.started":"2021-08-19T19:27:19.115416Z","shell.execute_reply":"2021-08-19T19:27:19.7305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let us now explore stock ID 0 at all time IDs and try to spot the patterns within our features, if there exists any in the first place.","metadata":{}},{"cell_type":"code","source":"df_example = df_features.loc[(df_features['stock_id'] == 0)]\n\ndf_example.plot.scatter(x = \"time_id\", y = 'target')","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:27:19.733621Z","iopub.execute_input":"2021-08-19T19:27:19.733911Z","iopub.status.idle":"2021-08-19T19:27:20.005956Z","shell.execute_reply.started":"2021-08-19T19:27:19.733882Z","shell.execute_reply":"2021-08-19T19:27:20.004793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Appears as though there is no correlation between time ID and realized volatility, so we will go on looking at the more relevant features we have gathered.","metadata":{}},{"cell_type":"code","source":"feature_array = ['bas', 'wap', 'lr', 'rv', 'seconds_in_bucket', \n                 'bid_price1', 'ask_price1', 'bid_price2', 'ask_price2', \n                 'bid_size1', 'ask_size1', 'bid_size2', 'ask_size2']\n\nfor feature in feature_array:\n    df_example.plot.scatter(x = feature, y = 'target')\n\ncorr = df_example[feature_array].corr()\nfig, ax = plt.subplots(figsize=(15, 15))\nax.matshow(corr)\nplt.xticks(range(len(corr.columns)), corr.columns)\nplt.yticks(range(len(corr.columns)), corr.columns)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:27:20.007494Z","iopub.execute_input":"2021-08-19T19:27:20.007785Z","iopub.status.idle":"2021-08-19T19:27:22.657837Z","shell.execute_reply.started":"2021-08-19T19:27:20.007756Z","shell.execute_reply":"2021-08-19T19:27:22.657126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There appears to be a relationship that we can capture with regression models between realized volatility, bid ask spread (and associated features), and past realized volatility. We will add the other features, which appear to have little correlation with realized volatility, to see if our later models can catch a pattern that we cannot see. The data also appears to need standardization before we train our models on them because of how small and spread out the feature data is. We will now apply standardization to all relevant features.","metadata":{}},{"cell_type":"markdown","source":"**Apply Standard Scaling**","metadata":{}},{"cell_type":"code","source":"scalers = []\nfor i in range(len(feature_array)):\n    scaler = StandardScaler()\n    feature = feature_array[i]\n    df_features[feature] = scaler.fit_transform(np.asarray(df_features[feature]).reshape(-1,1))\n    scalers.append(scaler)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:27:22.658794Z","iopub.execute_input":"2021-08-19T19:27:22.659151Z","iopub.status.idle":"2021-08-19T19:27:22.759027Z","shell.execute_reply.started":"2021-08-19T19:27:22.659124Z","shell.execute_reply":"2021-08-19T19:27:22.758197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Examine Data spread for stock ID 0 features again\ndf_example = df_features.loc[(df_features['stock_id'] == 0)]\n\nfor feature in feature_array:\n    df_example.plot.scatter(x = feature, y = 'target')\n\ncorr = df_example[feature_array].corr()\nfig, ax = plt.subplots(figsize=(15, 15))\nax.matshow(corr)\nplt.xticks(range(len(corr.columns)), corr.columns)\nplt.yticks(range(len(corr.columns)), corr.columns)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:27:22.759976Z","iopub.execute_input":"2021-08-19T19:27:22.760364Z","iopub.status.idle":"2021-08-19T19:27:25.897912Z","shell.execute_reply.started":"2021-08-19T19:27:22.760323Z","shell.execute_reply":"2021-08-19T19:27:25.896501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Even with standardization, our features appear to follow the same pattern in relation to the target as they did pre-standardization. We will try out both standardized and pre-standardized training data in our following models.","metadata":{}},{"cell_type":"markdown","source":"**Linear Regression and XGBoost Regression**","metadata":{}},{"cell_type":"code","source":"# Train Test Splitting\n\nX = df_features[feature_array]\nY = df_features['target']\n\nX_train, X_test, Y_train, Y_test = train_test_split(X, Y, train_size = .9)\n\n# Train Test Splitting for pre-standardized data\n\nX_copy = df_features_copy[feature_array]\nY_copy = df_features_copy['target']\n\nX_train_copy, X_test_copy, Y_train_copy, Y_test_copy = train_test_split(X_copy, Y_copy, train_size = .9)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.319087Z","iopub.status.idle":"2021-08-19T19:17:36.319694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Linear Regression Model Fitting \n\nlinreg_copy = LinearRegression().fit(X_train_copy, Y_train_copy)\nprint(\"Pre Standardization\")\nprint(\"R2 Score: \")\nprint(linreg_copy.score(X_train_copy, Y_train_copy))\nprint(\"Coefficients: \")\nprint(linreg_copy.coef_)\nprint(\"Intercepts: \")\nprint(linreg_copy.intercept_)\n\nlinreg = LinearRegression().fit(X_train, Y_train)\nprint(\"Post Standardization\")\nprint(\"R2 Score: \")\nprint(linreg.score(X_train, Y_train))\nprint(\"Coefficients: \")\nprint(linreg.coef_)\nprint(\"Intercepts: \")\nprint(linreg.intercept_)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.320819Z","iopub.status.idle":"2021-08-19T19:17:36.321187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# XGBoost Model Fitting \n\nregressor_copy = xgb.XGBRegressor(\n    n_estimators=50,\n    reg_lambda=1,\n    gamma=0,\n    max_depth=20\n)\n\nregressor_copy.fit(X_train_copy, Y_train_copy)\n\n\nregressor = xgb.XGBRegressor(\n    n_estimators=50,\n    reg_lambda=1,\n    gamma=0,\n    max_depth=20\n)\n\nregressor.fit(X_train, Y_train)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.322182Z","iopub.status.idle":"2021-08-19T19:17:36.32263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Examine Feature Importance\nprint(\"Pre Standardization Feature Importance\")\nprint(pd.DataFrame(regressor_copy.feature_importances_.reshape(1, -1), columns=feature_array))\n\nprint(\"Post Standardization Feature Importance\")\nprint(pd.DataFrame(regressor.feature_importances_.reshape(1, -1), columns=feature_array))","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.323991Z","iopub.status.idle":"2021-08-19T19:17:36.324531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As expected, past bid ask spread and realized volatility had the most importance in determining future, followed by Log Return, and our other features with much less importance.","metadata":{}},{"cell_type":"code","source":"# Calculate Root Mean Squared Prediction Error (RMSPE)\n\ndef RMSPE(actual, predict):\n    return (np.sqrt(np.mean(np.square((actual - predict) / actual))))\n\nY_predict_copy = linreg_copy.predict(X_test_copy)\nerror = RMSPE(Y_test_copy, Y_predict_copy)\nprint(\"Pre Standardization\")\nprint(\"Linear Regression Error:\")\nprint(error)\nY_predict_copy = regressor_copy.predict(X_test_copy)\nerror = RMSPE(Y_test_copy, Y_predict_copy)\nprint(\"XGBoost Error:\")\nprint(error)\n\n\nY_predict = linreg.predict(X_test)\nerror = RMSPE(Y_test, Y_predict)\nprint(\"Post Standardization\")\nprint(\"Linear Regression Error:\")\nprint(error)\nY_predict = regressor.predict(X_test)\nerror = RMSPE(Y_test, Y_predict)\nprint(\"XGBoost Error:\")\nprint(error)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.325813Z","iopub.status.idle":"2021-08-19T19:17:36.326173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Higher than the error from the other naive approach given to us in Optiver's tutorial, but not terrible for a first approach. There are still quite a few things we could obviously optimize in terms of our XGBoost Regression model, which has almost default settings right now. Additionally, standardization appeared to help in our models (the results above may not show it, but the majority of times that this has run, standardization has helped), so we will use that for our predictions.","metadata":{}},{"cell_type":"markdown","source":"Since our XGBRegressor model with vanilla parameters did fairly well (below .4 RMSPE for both standardized and non-standardized data), we will be using the same model with tuned hyperparameters to make our final test predictions.","metadata":{}},{"cell_type":"code","source":"# XGBRegressor Hyperparameter Tuning\n\n\"\"\"params = {\n    'reg_lambda': [1],\n    'gamma': [0],\n    'learning_rate': [0.01, 0.1],\n    'min_child_weight': [1, 3, 5],\n    'n_estimators' : [50 ,100, 200],\n    'max_depth': [20, 40, 60],\n    'objective': ['reg:squarederror']\n}\n\nxgb_regressor = xgb.XGBRegressor()\n\ngsearch = GridSearchCV(estimator = xgb_regressor,\n                       param_grid = params,\n                       # 2 Stratified Kfold\n                       cv = 2,\n                       # Use All Processors\n                       n_jobs = -1,\n                       # Display time taken per fold and paramater candidate is displayed\n                       # Display score\n                       verbose = 4)\n\ngsearch.fit(X_train,Y_train)\n\nprint(gsearch.best_params_)\"\"\"\n\n# May take 5+ hours to run. Run this only once if needed.","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.327356Z","iopub.status.idle":"2021-08-19T19:17:36.327707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Test Best Hyperparameter Configuration on Validation Data\n\nregressor = xgb.XGBRegressor(\n    reg_lambda=1,\n    gamma=0,\n    max_depth=20,\n    learning_rate=0.1,  \n    min_child_weight=3, \n    n_estimators=200,\n    objective='reg:squarederror'\n)\n\nregressor.fit(X_train, Y_train)\n\nY_predict = regressor.predict(X_test)\nerror = RMSPE(Y_test, Y_predict)\nprint(\"Hyperparameter Tuned XGBoost Error:\")\nprint(error)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.328899Z","iopub.status.idle":"2021-08-19T19:17:36.329255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will now attempt to create a regression model from a Deep Neural Network (DNN) that aims to lower the best root mean squared error from our previous models. Due to computing constraints, we will only use standardized data on a single DNN.\n\nThe model will be as follows:\n\n1 Input Layer with Relu Activation, 128 Units, and a input dimension that matches our training data shape\n \n1 Hidden Layer with Relu Activation and 256 Units\n\n1 Hidden Layer with Relu Activation and 128 Units\n\n1 Hidden Layer with Relu Activation and 64 Units\n\n1 Hidden Layer with Relu Activation and 32 Units\n\n1 Hidden Layer with Relu Activation and 16 Units\n\n1 Output Layer with Linear Activation and 1 Unit to match our output data shape\n\nAll layers will be kernel initalized to normal to match our standardized data and we will use the Adam optimizer as it appears to have give us the smoothest downward trend in our loss and validation loss history plot.","metadata":{}},{"cell_type":"code","source":"# Create RMSPE Function for the keras model\ndef rmspe(y_true, y_pred):\n    return K.sqrt(K.mean(K.square((y_true - y_pred) / y_true)))\n\n# Load Sequential to start forming our DNN\nDNN = Sequential()\n\n# The 128 Unit Input Layer:\nDNN.add(Dense(128, kernel_initializer='normal',input_dim = X_train.shape[1], activation='relu'))\n\n# The 256 Unit Hidden Layer:\nDNN.add(Dense(256, kernel_initializer='normal',activation='relu'))\n\n# The 128 Unit Hidden Layer:\nDNN.add(Dense(128, kernel_initializer='normal',activation='relu'))\n\n# The 64 Unit Hidden Layer:\nDNN.add(Dense(64, kernel_initializer='normal',activation='relu'))\n\n# The 32 Unit Hidden Layer\nDNN.add(Dense(32, kernel_initializer='normal',activation='relu'))\n\n# The 16 Unit Hidden Layer\nDNN.add(Dense(16, kernel_initializer='normal',activation='relu'))\n\n# The 1 Unit Outer Layer:\nDNN.add(Dense(1, kernel_initializer='normal',activation='linear'))\n\n# Compile the network:\nDNN.compile(loss=rmspe, optimizer='adam')\nDNN.summary()","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.330394Z","iopub.status.idle":"2021-08-19T19:17:36.330848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Initialize Callbacks\n\"\"\"my_callbacks = [\n    tf.keras.callbacks.EarlyStopping(patience=2),\n    tf.keras.callbacks.ModelCheckpoint(filepath='model.{epoch:02d}-{val_loss:.2f}.h5'),\n    tf.keras.callbacks.TensorBoard(log_dir='./logs'),\n]\"\"\"\n\n# Fit the DNN to our training data\nhistory = DNN.fit(X, Y, epochs=50, batch_size=128, validation_split = 0.2, verbose=1)","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.331873Z","iopub.status.idle":"2021-08-19T19:17:36.332221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot loss and validation loss per epoch\npd.DataFrame(history.history).plot(figsize=(8,5))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.33496Z","iopub.status.idle":"2021-08-19T19:17:36.335328Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our DNN seems to have beaten our previous regression models in RMSPE by quite a bit, so we will be using it in our final predictions.","metadata":{}},{"cell_type":"markdown","source":" **Predictions for the true test set**","metadata":{}},{"cell_type":"code","source":"# Import test Data\n\ntest = pd.read_csv('../input/optiver-realized-volatility-prediction/test.csv')\n\n# Grab Unique Array IDs\nstock_id_array = test.stock_id.unique()\n\ndf_bas = bas_calculation(stock_id_array, 'test')\ndf_wap = wap_calculation(stock_id_array, 'test')\ndf_lr = lr_calculation(stock_id_array, 'test')\ndf_rv = rv_calculation(stock_id_array, 'test')\ndf_other = other_calculation(stock_id_array, 'test')\n\n# Concatonate each statistic into one dataframe\n\ndf_features = df_bas.merge(df_wap,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(df_lr,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(df_rv,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(df_other,on=['stock_id', 'time_id'])\ndf_features = df_features.merge(test,on=['stock_id', 'time_id'])\n\n# Standardize the feature data\n\nfor i in range(len(feature_array)):\n    scaler = scalers[i]\n    feature = feature_array[i]\n    df_features[feature] = scaler.transform(np.asarray(df_features[feature]).reshape(-1,1))\n\n# Make Predictions\nX_predict = df_features[feature_array]\nY_predict = DNN.predict(X_predict)\n\ndata = []\nfor i in range(len(Y_predict)):\n    data.append([str(df_features['stock_id'].iloc[i]) + \"-\" + str(df_features['time_id'].iloc[i]), Y_predict[i][0]])\ndf = pd.DataFrame(data, columns = ['row_id', 'target'])\ndf.to_csv('submission.csv', index=False)\npd.read_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2021-08-19T19:17:36.336178Z","iopub.status.idle":"2021-08-19T19:17:36.336543Z"},"trusted":true},"execution_count":null,"outputs":[]}]}