{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":84493,"databundleVersionId":9871156,"sourceType":"competition"},{"sourceId":216354576,"sourceType":"kernelVersion"}],"dockerImageVersionId":30823,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What are the significant relationships between features and responder_6, and how can statistical tests validate or invalidate these relationships? #","metadata":{}},{"cell_type":"markdown","source":"## Dataset Description\n\n\nThanks to Jane Street 2024 competition, the dataset consists of timeseries data with 79 anonymized features and 9 responders, derived from real market data. The task is to forecast responder_6 for up to six months. \n\nTrain Data: Historical data partitioned into ten parts, with columns for date_id, time_id, symbol_id, 79 feature_{00...78} columns, and 9 responder_{0...8} columns.\nTest Data: A mock test set structured similarly to the training data, provided in batches via an evaluation API.\nLags: Previous day’s responder values, available at the first time_id of each new day.\n\n\nFor this EDA and Statistical analysis i will only be considering the last partition, which might contain more recent data that may align better with the patterns in the unseen test set. ","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"This notebook forked from: \nhttps://www.kaggle.com/code/yanisbelami/jane-street-real-time-market-data-forecasting-eda\nplease Upvote me and it.","metadata":{}},{"cell_type":"code","source":"%%time\nimport os\nimport gc\nimport random\nimport pandas as pd\nimport polars as pl\nimport numpy as np\nfrom matplotlib import pyplot as plt\nimport seaborn as sns\nfrom math import log, sqrt\nfrom scipy.stats import (ttest_1samp, pearsonr, spearmanr, f_oneway, shapiro, kstest, norm, chi2)\nimport lightgbm as lgb\nfrom lightgbm import LGBMRegressor\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.inspection import permutation_importance\nfrom sklearn.model_selection import train_test_split\nimport statsmodels.api as sm\nimport kaggle_evaluation.jane_street_inference_server\nimport warnings\n\nwarnings.filterwarnings('ignore')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:07:36.099247Z","iopub.execute_input":"2025-01-06T23:07:36.099570Z","iopub.status.idle":"2025-01-06T23:07:41.500623Z","shell.execute_reply.started":"2025-01-06T23:07:36.099545Z","shell.execute_reply":"2025-01-06T23:07:41.499860Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"feature_cols = [f\"feature_{idx:02d}\" for idx in range(79)]+ [f\"responder_{idx}_lag_1\" for idx in range(9)]\ntarget_col = \"responder_6\"\nselected_features = [\"symbol_id\", \"time_id\"] + feature_cols+[target_col]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:13:42.580356Z","iopub.execute_input":"2025-01-06T23:13:42.580687Z","iopub.status.idle":"2025-01-06T23:13:42.584759Z","shell.execute_reply.started":"2025-01-06T23:13:42.580660Z","shell.execute_reply":"2025-01-06T23:13:42.583904Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"responders = pd.read_csv('/kaggle/input/jane-street-real-time-market-data-forecasting/responders.csv')\nfeatures = pd.read_csv('/kaggle/input/jane-street-real-time-market-data-forecasting/features.csv')\nsample_submission = pd.read_csv('/kaggle/input/jane-street-real-time-market-data-forecasting/sample_submission.csv')\n\ntrain = (\n    pl.read_parquet('/kaggle/input/jane-street-real-time-market-data-forecasting/train.parquet/partition_id=9/part-0.parquet')\n)\n\ntrain = train.to_pandas()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:05.020781Z","iopub.execute_input":"2025-01-06T16:23:05.021089Z","iopub.status.idle":"2025-01-06T16:23:20.466162Z","shell.execute_reply.started":"2025-01-06T16:23:05.021059Z","shell.execute_reply":"2025-01-06T16:23:20.465405Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Basic introduction ###\nPartition 9 consists of 6274576 rows and 92 columns","metadata":{}},{"cell_type":"code","source":"print(\"Shape of the dataset:\", train.shape)\nprint(train.info())\n#print(train.describe())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:20.466898Z","iopub.execute_input":"2025-01-06T16:23:20.467158Z","iopub.status.idle":"2025-01-06T16:23:20.496128Z","shell.execute_reply.started":"2025-01-06T16:23:20.467136Z","shell.execute_reply":"2025-01-06T16:23:20.495111Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ntrain_data = pl.scan_parquet(f\"/kaggle/input/js24-dataset-stats-with-lags/training.parquet/\").collect().to_pandas()\nval_data = pl.scan_parquet(f\"/kaggle/input/js24-dataset-stats-with-lags/validation.parquet/\").collect().to_pandas()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:07:57.836077Z","iopub.execute_input":"2025-01-06T23:07:57.836710Z","iopub.status.idle":"2025-01-06T23:08:37.526250Z","shell.execute_reply.started":"2025-01-06T23:07:57.836683Z","shell.execute_reply":"2025-01-06T23:08:37.525437Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def reduce_memory_usage(df):\n    \"\"\" \n    iterate through all the columns of a dataframe and \n    modify the data type to reduce memory usage.        \n    \"\"\"\n    start_mem = df.memory_usage().sum() / 1024**3\n    print(('Memory usage of dataframe is {:.2f}' \n                     'GB').format(start_mem))\n    \n    for col in df.columns:\n        col_type = df[col].dtype\n        \n        if col_type != object:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max <\\\n                  np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max <\\\n                   np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max <\\\n                   np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max <\\\n                   np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)  \n            else:\n                if c_min > np.finfo(np.float16).min and c_max <\\\n                   np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max <\\\n                   np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)\n        else:\n            df[col] = df[col].astype('category')\n    end_mem = df.memory_usage().sum() / 1024**3\n    print(('Memory usage after optimization is: {:.2f}' \n                              'GB').format(end_mem))\n    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) \n                                             / start_mem))\n    \n    return df\n\ndef percentage_missing_values(df):\n    missing_values_count = df.isnull().sum()\n    total_cells = np.product(df.shape)\n    total_missing = missing_values_count.sum()\n    print (\"Percentage of Missing Data = \",(total_missing/total_cells) * 100,\"%\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:09:04.117701Z","iopub.execute_input":"2025-01-06T23:09:04.117997Z","iopub.status.idle":"2025-01-06T23:09:04.127610Z","shell.execute_reply.started":"2025-01-06T23:09:04.117975Z","shell.execute_reply":"2025-01-06T23:09:04.126679Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ntrain_data = reduce_memory_usage(train_data)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:09:10.073948Z","iopub.execute_input":"2025-01-06T23:09:10.074286Z","iopub.status.idle":"2025-01-06T23:09:31.834949Z","shell.execute_reply.started":"2025-01-06T23:09:10.074259Z","shell.execute_reply":"2025-01-06T23:09:31.834199Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"features_has_nan=[]\nfeatures_not_nan=[]\nfor col in train_data.columns:\n    if train_data[col].isna().sum()>0:\n        features_has_nan+=[col]\n        print(f'{col}: {train_data[col].isna().sum()}')\n    else:\n        features_not_nan+=[col]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:10:50.488279Z","iopub.execute_input":"2025-01-06T23:10:50.488622Z","iopub.status.idle":"2025-01-06T23:11:01.046820Z","shell.execute_reply.started":"2025-01-06T23:10:50.488599Z","shell.execute_reply":"2025-01-06T23:11:01.045976Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(15, 5))\nbalance= pd.Series(train_data[target_col]).cumsum()\nax.set_xlabel(\"Trade\", fontsize=18)\nax.set_ylabel(\"Cumulative resp\", fontsize=18);\nbalance.plot(lw=3);\ndel balance\ngc.collect();","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:14:06.979320Z","iopub.execute_input":"2025-01-06T23:14:06.979668Z","iopub.status.idle":"2025-01-06T23:14:10.715969Z","shell.execute_reply.started":"2025-01-06T23:14:06.979641Z","shell.execute_reply":"2025-01-06T23:14:10.715128Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(15, 5))\nbalance= pd.Series(val_data[target_col]).cumsum()\nax.set_xlabel(\"Trade\", fontsize=18)\nax.set_ylabel(\"Cumulative resp\", fontsize=18);\nbalance.plot(lw=3);\ndel balance\ngc.collect();","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:14:20.788285Z","iopub.execute_input":"2025-01-06T23:14:20.788593Z","iopub.status.idle":"2025-01-06T23:14:21.360555Z","shell.execute_reply.started":"2025-01-06T23:14:20.788563Z","shell.execute_reply":"2025-01-06T23:14:21.359681Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport pandas as pd\nimport seaborn as sns\nimport gc\n\ndef plot_weighted_resps(df, cols=['responder_6'], marker_day=85, start_date=None, use_weight=True):\n    \"\"\"\n    Plots the cumulative daily return for weighted or unweighted responses over specified time horizons.\n\n    Parameters:\n        df (pd.DataFrame): DataFrame containing the data with 'weight' and specified horizon columns.\n        cols (list): List of horizon column names (e.g., 'responder_6', etc.).\n        marker_day (int): Day to mark with a vertical line and shaded region.\n        start_date (int or None): Start date for the x-axis. If None, defaults to the minimum date in the data.\n        use_weight (bool): If True, uses weight in the calculation; otherwise, plots unweighted cumulative return.\n\n    Returns:\n        None: Displays the plot.\n    \"\"\"\n    if start_date is None:\n        start_date = df['date_id'].min()  # Default to the minimum date in the dataset\n        marker_day+=start_date\n\n    fig, ax = plt.subplots(figsize=(15, 5))\n    \n    # Compute and plot cumulative daily returns for each column in `cols`\n    for horizon in cols:\n        if use_weight:\n            weighted_resp_col = f'weight_{horizon}'\n            if weighted_resp_col not in df.columns:\n                # Compute weighted response if not already in DataFrame\n                df[weighted_resp_col] = df['weight'] * df[horizon]\n            col_to_plot = weighted_resp_col\n        else:\n            col_to_plot = horizon\n        \n        # Group by date, compute mean, and calculate cumulative product\n        cumulative_return = pd.Series(1 + df.groupby('date_id')[col_to_plot].mean()).cumprod()\n        cumulative_return.index -= start_date  # Adjust the x-axis to start from `start_date`\n        label = f\"{horizon} {'x weight' if use_weight else ''}\"\n        cumulative_return.plot(lw=3, label=label, ax=ax)\n        \n        # Clean up memory\n        del cumulative_return\n        gc.collect()\n    \n    # Add labels and title\n    ax.set_xlabel(\"Day\", fontsize=18)\n    ax.set_title(f\"Cumulative daily return for {'weighted' if use_weight else 'unweighted'} responses over time horizons\", fontsize=18)\n    \n    # Adjust marker day to align with the new x-axis starting point\n    adjusted_marker_day = marker_day - start_date\n    ax.axvline(x=adjusted_marker_day, linestyle='--', alpha=0.3, c='red', lw=1, label=f'Marker Day {marker_day}')\n    ax.axvspan(0, adjusted_marker_day, color=sns.xkcd_rgb['grey'], alpha=0.1)\n    \n    # Add legend\n    plt.legend(loc=\"lower left\")\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:16:01.258762Z","iopub.execute_input":"2025-01-06T23:16:01.259129Z","iopub.status.idle":"2025-01-06T23:16:01.266657Z","shell.execute_reply.started":"2025-01-06T23:16:01.259093Z","shell.execute_reply":"2025-01-06T23:16:01.265793Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_weighted_resps(train_data,cols=[target_col],use_weight=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:16:52.524091Z","iopub.execute_input":"2025-01-06T23:16:52.524451Z","iopub.status.idle":"2025-01-06T23:16:53.713324Z","shell.execute_reply.started":"2025-01-06T23:16:52.524424Z","shell.execute_reply":"2025-01-06T23:16:53.712472Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_weighted_resps(train_data,cols=[target_col],use_weight=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:17:10.094151Z","iopub.execute_input":"2025-01-06T23:17:10.094496Z","iopub.status.idle":"2025-01-06T23:17:11.035018Z","shell.execute_reply.started":"2025-01-06T23:17:10.094471Z","shell.execute_reply":"2025-01-06T23:17:11.034008Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_weighted_resps(val_data,marker_day=10,use_weight=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:19:02.622103Z","iopub.execute_input":"2025-01-06T23:19:02.622443Z","iopub.status.idle":"2025-01-06T23:19:02.960589Z","shell.execute_reply.started":"2025-01-06T23:19:02.622418Z","shell.execute_reply":"2025-01-06T23:19:02.959780Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_weighted_resps(val_data,marker_day=10,use_weight=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T23:19:41.298145Z","iopub.execute_input":"2025-01-06T23:19:41.298465Z","iopub.status.idle":"2025-01-06T23:19:41.773418Z","shell.execute_reply.started":"2025-01-06T23:19:41.298438Z","shell.execute_reply":"2025-01-06T23:19:41.772520Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Missing values\n\n33 features have missing values, with 4 features having more than 7% missing values","metadata":{}},{"cell_type":"code","source":"missing_values = train.isnull().sum()\nmissing_percentage = (missing_values / train.shape[0]) * 100\n\nmissing_data = pd.DataFrame({\n    \"Feature\": train.columns,\n    \"MissingCount\": missing_values,\n    \"MissingPercentage\": missing_percentage\n}).sort_values(by=\"MissingPercentage\", ascending=False)\n\n# threshold > 0% 33 features\n# threshold > 1% 13 features\n# threshold > 5% 4 features\nmissing_data_filtered = missing_data[missing_data[\"MissingPercentage\"] > 0]\n\nplt.figure(figsize=(14, 8))\nsns.barplot(\n    x=\"MissingPercentage\",\n    y=\"Feature\",\n    data=missing_data_filtered,\n    palette=\"vlag\"\n)\nplt.title(\"Missing Values Analysis\", fontsize=16)\nplt.xlabel(\"Percentage of Missing Values\", fontsize=14)\nplt.ylabel(\"Features\", fontsize=14)\nplt.show()\n\nprint(\"Number of features with missing values : \", len(missing_data_filtered))\nmissing_data_filtered","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:20.496842Z","iopub.execute_input":"2025-01-06T16:23:20.497147Z","iopub.status.idle":"2025-01-06T16:23:21.624259Z","shell.execute_reply.started":"2025-01-06T16:23:20.497126Z","shell.execute_reply":"2025-01-06T16:23:21.623563Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Responder_6 distribution","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(8, 5))\nsns.histplot(train['responder_6'], bins=100, kde=True,fill=False)\nplt.title(\"Distribution of Responder_6\")\nplt.xlabel(\"Responder_6\")\nplt.ylabel(\"Frequency\")\nplt.show()\n\nall_response = train['responder_6'].mean() * 100\nprint(f\"Proportion of all Responder_6: {all_response:.2f}%\")\n\nzero_response = (train['responder_6'] == 0).mean() * 100\nprint(f\"Proportion of exact zeros in Responder_6: {zero_response:.2f}%\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:21.625065Z","iopub.execute_input":"2025-01-06T16:23:21.625334Z","iopub.status.idle":"2025-01-06T16:23:45.374549Z","shell.execute_reply.started":"2025-01-06T16:23:21.625312Z","shell.execute_reply":"2025-01-06T16:23:45.373668Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The mean of responder_6 is -0.38%, indicating a slight negative bias, and also no exact zeros are present in the data.\nThe distribution of responder_6 looks roughly symmetric and bell-shaped, kind of like a normal distribution. But to be sure, it could be interesting to run some normality test, to check later if it really follows a normal law, even if it seems useless in our real world, https://stats.stackexchange.com/questions/2492/is-normality-testing-essentially-useless.","metadata":{}},{"cell_type":"markdown","source":"## Date_id and Time_id","metadata":{}},{"cell_type":"code","source":"daily_mean = train.groupby('date_id')['responder_6'].mean()\nplt.figure(figsize=(12, 6))\nplt.plot(daily_mean)\nplt.title(\"Date_id Mean of Responder_6 Over Time\")\nplt.xlabel(\"Date ID\")\nplt.ylabel(\"Mean Responder_6\")\nplt.grid()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:45.375227Z","iopub.execute_input":"2025-01-06T16:23:45.375445Z","iopub.status.idle":"2025-01-06T16:23:45.74648Z","shell.execute_reply.started":"2025-01-06T16:23:45.375426Z","shell.execute_reply":"2025-01-06T16:23:45.745563Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"daily_std = train.groupby('date_id')['responder_6'].std()\nplt.figure(figsize=(12, 6))\nplt.plot(daily_std)\nplt.title(\"Date_id Std of Responder_6 Over Time\")\nplt.xlabel(\"Date ID\")\nplt.ylabel(\"Mean Responder_6\")\nplt.grid()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:45.748447Z","iopub.execute_input":"2025-01-06T16:23:45.748668Z","iopub.status.idle":"2025-01-06T16:23:46.124917Z","shell.execute_reply.started":"2025-01-06T16:23:45.748649Z","shell.execute_reply":"2025-01-06T16:23:46.124118Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"intraday_mean = train.groupby('time_id')['responder_6'].mean()\nplt.figure(figsize=(12, 6))\nplt.plot(intraday_mean)\nplt.title(\"Time_id Mean of Responder_6\")\nplt.xlabel(\"Time ID\")\nplt.ylabel(\"Mean Responder_6\")\nplt.grid()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:46.126214Z","iopub.execute_input":"2025-01-06T16:23:46.126472Z","iopub.status.idle":"2025-01-06T16:23:46.469435Z","shell.execute_reply.started":"2025-01-06T16:23:46.12645Z","shell.execute_reply":"2025-01-06T16:23:46.468605Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The dataset has in total for partition_6 169 days and 968 time steps per day every day.","metadata":{}},{"cell_type":"code","source":"print(\"Unique date-id:\", train['date_id'].nunique())\nprint(\"Total Unique time_id steps:\", train['time_id'].nunique())\n\n\ntime_per_day = train.groupby('date_id')['time_id'].nunique()\nprint(\"\\nTime Steps per Day:\")\nprint(time_per_day.describe()) \n\n\n### time_id per date_id\nrows_per_day = train['date_id'].value_counts()\nprint(\"\\nRows per Day:\")\nprint(rows_per_day.describe())\n\nplt.figure(figsize=(10, 6))\nsns.histplot(rows_per_day, bins=100, kde=True)\nplt.title(\"Number of Tick per Date_id\")\nplt.xlabel(\"Number of Rows\")\nplt.ylabel(\"Frequency\")\nplt.grid()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:46.470232Z","iopub.execute_input":"2025-01-06T16:23:46.470441Z","iopub.status.idle":"2025-01-06T16:23:47.508963Z","shell.execute_reply.started":"2025-01-06T16:23:46.470422Z","shell.execute_reply":"2025-01-06T16:23:47.508107Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Checking for a specific date_id 1530","metadata":{}},{"cell_type":"code","source":"specific_date = 1530\ntrain_one_day = train[train['date_id'] == specific_date]\nprint(f\"Number of time_id for date_id 1530: {train_one_day.shape[0]}\")\n\nprint(\"from \", train_one_day['time_id'].min(), \"to\", train_one_day['time_id'].max())\n\nplt.figure(figsize=(12, 6))\nplt.plot(train_one_day['time_id'], train_one_day['responder_6'])\nplt.title(f\"Responder_6 Over Time for date_id {specific_date}\")\nplt.xlabel(\"Time ID\")\nplt.ylabel(\"Responder_6\")\nplt.grid()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:47.510068Z","iopub.execute_input":"2025-01-06T16:23:47.510406Z","iopub.status.idle":"2025-01-06T16:23:47.733063Z","shell.execute_reply.started":"2025-01-06T16:23:47.51037Z","shell.execute_reply":"2025-01-06T16:23:47.732201Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The weight distribution is right-skewed, with most values between 1 and 3, which should be analyzed further.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(8, 5))\nsns.histplot(train['weight'], bins=100, kde=True,fill=False)\nplt.title(\"Distribution of Weights\")\nplt.xlabel(\"Weight\")\nplt.ylabel(\"Frequency\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:23:47.733861Z","iopub.execute_input":"2025-01-06T16:23:47.734196Z","iopub.status.idle":"2025-01-06T16:24:11.360508Z","shell.execute_reply.started":"2025-01-06T16:23:47.734164Z","shell.execute_reply":"2025-01-06T16:24:11.358192Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Symbol_id \n\nThe cumulative returns of responder_6 show significant divergence across all symbol_id values. Symbols exhibit varying levels of volatility and form clusters of similar behaviors. This suggests that symbol_id influences the performance of responder_6 over time and warrants further analysis. And the big peak at circa 1.2 seems to be correlated with the same high std seen in date_id plot. Finding this exact value and delete it may help.","metadata":{}},{"cell_type":"code","source":"train['id'] = train.index.values\n\nplt.figure(figsize=(16, 8))\nfor symbol_id in train['symbol_id'].unique():\n    xx = train[train['symbol_id'] == symbol_id]['id']\n    yy = train[train['symbol_id'] == symbol_id]['responder_6']\n    plt.plot(xx, yy.cumsum(), label=f'Symbol ID {symbol_id}', linewidth=0.5)\n\nplt.title('Cumulative responder_6 for All Symbol IDs', fontsize=16)\nplt.xlabel(\"Time\", fontsize=12)\nplt.ylabel(\"Cumulative Returns\", fontsize=12)\nplt.grid(color='lightgray', linewidth=0.5)\nplt.axhline(0, color='red', linestyle='-', linewidth=0.7)\nplt.legend(loc='upper left', fontsize=8)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.36117Z","iopub.status.idle":"2025-01-06T16:24:11.361457Z","shell.execute_reply":"2025-01-06T16:24:11.361345Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"All responders cumulative sum comparison, Responder_6 has the same pattern as Responder_7 and Responder_8.","metadata":{}},{"cell_type":"code","source":"responders = ['responder_0', 'responder_1', 'responder_2', 'responder_3', 'responder_4', 'responder_5', 'responder_6', 'responder_7', 'responder_8']\nplt.figure(figsize=(16, 8))\n\nfor responder in responders:\n    cumulative_sum = train.groupby('id')[responder].sum().cumsum()\n    plt.plot(cumulative_sum, label=f'{responder}', linewidth=0.8)\n\nplt.title('Cumulative Sum of All Responders', fontsize=16)\nplt.xlabel(\"Time\", fontsize=12)\nplt.ylabel(\"Cumulative Returns\", fontsize=12)\nplt.grid(color='lightgray', linewidth=0.5)\nplt.axhline(0, color='red', linestyle='-', linewidth=0.7)\nplt.legend(loc='upper left', fontsize=8)\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.362214Z","iopub.status.idle":"2025-01-06T16:24:11.36264Z","shell.execute_reply":"2025-01-06T16:24:11.362456Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Check whether Responder_6 follows a normal distribution using Kolmogorov-Smirnov normality test:","metadata":{}},{"cell_type":"markdown","source":"The **Kolmogorov-Smirnov** test compares the Empirical Cumulative Distribution Function of the sample data to the Cumulative Distribution Function of a reference distribution. The test statistic $ D $ is defined as:\n$$D = \\sup_x \\left| F_n(x) - F(x) \\right|$$\n\nWhere:\n- $F_n(x) $ is the Empirical CDF:\n$\nF_n(x) = \\frac{1}{n} \\sum_{i=1}^n \\mathbb{I}(x_i \\leq x)\n$\n- $ F(x) $ is the theoretical CDF  \n- $ \\sup_x $ is the supremum over all values of $ x $.\n\n**Hypotheses**:\n- $ H_0: F_n(x) = F(x) $ (Normal distribution)\n- $ H_1: F_n(x) \\neq F(x) $ (Not Normal distribution)\n\n**Interpretation**:\n- If $ p \\leq \\alpha $, reject $ H_0 $: The data does not follow a Normal distribution.\n- If $ p > \\alpha $, fail to reject $ H_0 $: The data may follow a Normal distribution.\n","metadata":{}},{"cell_type":"code","source":"responder_data = train['responder_6'].dropna()\nmean_responder = responder_data.mean()\nstd_responder = responder_data.std()\n\n\nks_statistic, p_value = kstest(responder_data, 'norm', args=(mean_responder, std_responder))\nprint(f\"KS Test Statistic: {ks_statistic:.9f}\")\nprint(f\"P-value: {p_value:.9e}\")\n\n\nif p_value > 0.05:\n    print(\"Fail to reject H0: The data appears to follow a normal distribution.\")\nelse:\n    print(\"Reject H0: The data does not follow a normal distribution.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.364202Z","iopub.status.idle":"2025-01-06T16:24:11.364595Z","shell.execute_reply":"2025-01-06T16:24:11.364421Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":" The test statistic was about 0.1, and the p-value was approximately 0. This result strongly rejects the null hypothesis, indicating that the responder_6 data does **not follow a normal distribution**.","metadata":{}},{"cell_type":"markdown","source":"## Features EDA and Correlogram","metadata":{}},{"cell_type":"code","source":"feature_columns = [col for col in train.columns if 'feature' in col]\n\nprint(\"Total Features:\", len(feature_columns))\nprint(\"Feature Columns:\", feature_columns)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.365186Z","iopub.status.idle":"2025-01-06T16:24:11.365564Z","shell.execute_reply":"2025-01-06T16:24:11.365388Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Bravais-Pearson correlation:**\nThe Pearson correlation measures the linear relationship between two variables $X$ and $Y$.\n\nIn the context of your dataset:\n$X$ represents one feature column\n$Y$ represents the target variable responder_6\n\nThe formula for the Pearson correlation coefficient is:\n$$\\rho_{xy} = \\frac{\\text{Cov}(X,Y)}{\\sigma_x \\sigma_y}  = \\frac{E[(X-\\mu_x)(Y-\\mu_y)]}{\\sigma_x \\sigma_y} = \\frac{\\sum_{i=1}^{n} (x_i - \\bar{x})(y_i - \\bar{y})}{\\sqrt{\\sum_{i=1}^{n} (x_i - \\bar{x})^2} \\sqrt{\\sum_{i=1}^{n} (y_i - \\bar{y})^2}} $$\n\nWhere:\n- $\\text{Cov}(X,Y)$: Covariance between $X$ and $Y$\n- $\\sigma_x, \\sigma_y$: Standard deviations of $X$ and $Y$\n- $\\mu_x, \\mu_y$: Means of $X$ and $Y$\n- $E$: mathematical expectation\n- $x_i$: Values of feature $X$ \n- $y_i$: Values of responder_6\n- $\\bar{x}, \\bar{y}$: Mean values of $X$ and $Y$\n- $n$: Total number of data points\n\n\nIn this analysis, I use the Pearson correlation by default. But for rank-based relationships, **Spearman correlation** would be considered, which measures monotonic relationships based on **rank**.\nThe formula is based on the ranks of $X$ and $Y$ instead of their raw values:\n$$\\rho = 1 - \\frac{6\\sum d_i^2}{n(n^2-1)} = 1 - \\frac{6\\sum_{i=1}^{n} (R(x_i) - R(y_i))^2}{n(n^2-1)}$$\n\nWhere:\n- $d_i = R(x_i) - R(y_i)$: Difference between the ranks of $x_i$ and $y_i$\n- $R(x_i)$: Rank of $x_i$\n- $R(y_i)$: Rank of $y_i$\n- $n$: Total number of observations","metadata":{}},{"cell_type":"code","source":"%%time\nfeature_corr = train[feature_columns].corr(method ='pearson')\n\nplt.figure(figsize=(24, 20))\nsns.heatmap(feature_corr, cmap=\"coolwarm\", annot=False, center=0)\nplt.title(\"Feature-to-Feature Correlation Heatmap\", fontsize=16)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.366098Z","iopub.status.idle":"2025-01-06T16:24:11.366487Z","shell.execute_reply":"2025-01-06T16:24:11.366305Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nfull_corr_s = train.corr(method='spearman')\n\nplt.figure(figsize=(24, 20))\nsns.heatmap(full_corr_s, cmap=\"coolwarm\", annot=False, center=0)\nplt.title(\"Full Correlation Heatmap\", fontsize=16)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.367071Z","iopub.status.idle":"2025-01-06T16:24:11.367421Z","shell.execute_reply":"2025-01-06T16:24:11.367265Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Their is no correlation between features and Responder_6, 0.03 being pretty low ","metadata":{}},{"cell_type":"code","source":"correlation_with_responder = train[feature_columns + ['responder_6']].corr()\n\nresponder_corr = correlation_with_responder['responder_6'].drop('responder_6').sort_values(ascending=True)\n\nplt.figure(figsize=(12, 6))\nsns.barplot(x=responder_corr.index, y=responder_corr.values)\nplt.xticks(rotation=90)\nplt.title(\"Correlation of Features with Responder_6\", fontsize=16)\nplt.xlabel(\"Features\")\nplt.ylabel(\"Correlation Coefficient\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.367934Z","iopub.status.idle":"2025-01-06T16:24:11.368347Z","shell.execute_reply":"2025-01-06T16:24:11.368174Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Full Heatmap just to see interactions between date_id, time_id, symbol_id and other responders, using pearson method ","metadata":{}},{"cell_type":"code","source":"%%time\nfull_corr_p = train.corr(method ='pearson')\n\nplt.figure(figsize=(24, 20))\nsns.heatmap(full_corr_p, cmap=\"coolwarm\", annot=False, center=0)\nplt.title(\"Full Correlation Heatmap\", fontsize=16)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.369657Z","iopub.status.idle":"2025-01-06T16:24:11.370062Z","shell.execute_reply":"2025-01-06T16:24:11.369882Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Responder Correlations:** \nresponder_6 is strongly correlated with responder_3, responder_7, and responder_8.\n\n**Feature Clusters:**-\nfeature_73 to feature_78 are highly correlated.\nfeature_50-60 each correlate with feature_39-49 respectively.\n\n**Moderate Correlations:**\nfeature_72 to feature_67 are all moderately correlated and also with feature_12, feature_13, and feature_14.\nfeature_65 with feature_18, feature_66 with feature_19.\n\n**Other Pairs:**\nfeature_02 and feature_03 with feature_00.\nfeature_22 with weight.","metadata":{}},{"cell_type":"markdown","source":"As found previously, responders 3, 8 and 7 may share a similar pattern with responder_6.\nAlso the negative correlations are too weak to indicate strong inverse relationships with responder_6. These columns are unlikely to play a significant role so they may have limited utility for the final model.","metadata":{}},{"cell_type":"code","source":"combined_corr =  full_corr_s\n#combined_corr = (full_corr_p + full_corr_p) / 2\n\ncombined_corr_responder = combined_corr['responder_6'].drop('responder_6').abs().sort_values(ascending=False)\ntop_10_features = combined_corr_responder.head(15)\n\nprint(\"Top 15 Features with highest mean (spearman) correlation:\")\nprint(top_10_features)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.371013Z","iopub.status.idle":"2025-01-06T16:24:11.371388Z","shell.execute_reply":"2025-01-06T16:24:11.371231Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"combined_corr = full_corr_p \n#combined_corr = (full_corr_p + full_corr_p) / 2\n\ncombined_corr_responder = combined_corr['responder_6'].drop('responder_6').abs().sort_values(ascending=False)\ntop_10_features = combined_corr_responder.head(15)\n\nprint(\"Top 15 Features with highest mean (pearson) correlation:\")\nprint(top_10_features)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.373438Z","iopub.status.idle":"2025-01-06T16:24:11.373814Z","shell.execute_reply":"2025-01-06T16:24:11.373655Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Could be a good idea to perform feature selection, and delete the ones that have very low correlations, to avoid adding unnecessary noise to the model.","metadata":{}},{"cell_type":"code","source":"combined_corr = (full_corr_p + full_corr_s) / 2\n#combined_corr = (full_corr_p + full_corr_p) / 2\n\ncombined_corr_responder = combined_corr['responder_6'].drop('responder_6').abs().sort_values(ascending=False)\ntop_10_features = combined_corr_responder.head(10)\nbottom_10_features = combined_corr_responder.tail(10)\n\nprint(\"Top 10 Features with highest mean (spearman + pearson) correlation:\")\nprint(top_10_features)\nprint(\"Bottom 10 Features with lowest mean (spearman + pearson) correlation:\")\nprint(bottom_10_features)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.374446Z","iopub.status.idle":"2025-01-06T16:24:11.375256Z","shell.execute_reply":"2025-01-06T16:24:11.375071Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"These highly correlated features may contain redundant information, which can impact modeling performance. Dimensionality reduction techniques like PCA or even feature deletion/selection may help improving the final model.","metadata":{}},{"cell_type":"code","source":"threshold = 0.9\nhigh_corr_pairs = []\n\nfor i in range(len(combined_corr.columns)):\n    for j in range(i):\n        if combined_corr.iloc[i, j] > threshold:\n            high_corr_pairs.append((combined_corr.columns[i], combined_corr.columns[j], combined_corr.iloc[i, j]))\n\n\nprint(f\"Highly Correlated Column Pairs (|corr| > {threshold}:\")\nfor pair in high_corr_pairs:\n    print(f\"{pair[0]} - {pair[1]}: {pair[2]:.2f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.376285Z","iopub.status.idle":"2025-01-06T16:24:11.376725Z","shell.execute_reply":"2025-01-06T16:24:11.376476Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We can talk about denser central areas where values cluster tightly near the mean, or even the present vertical symmetry, but as observed earlier, there is no strong correlation between the features and responder_6. Even for the features with relatively higher (but still weak) correlations, plotting them does not reveal any clear patterns or significant relationships.\n","metadata":{}},{"cell_type":"code","source":"combined_corr_responder = combined_corr['responder_6'].drop('responder_6').abs().sort_values(ascending=False)\ncombined_corr_features = combined_corr_responder[combined_corr_responder.index.str.contains('feature')].abs().sort_values(ascending=False)\n\ntop_9_features = combined_corr_features.head(9).index.tolist()\nprint(\"Top 9 Features Most Correlated with Responder_6 (features only):\", top_9_features)\n\nplt.figure(figsize=(16, 16))\nfor i, feature in enumerate(top_9_features, 1):\n    plt.subplot(3, 3, i)\n    plt.hexbin(train['responder_6'], train[feature],gridsize=1000, bins='log', cmap='inferno')\n    plt.xlabel(feature, fontsize=10)\n    plt.ylabel('Responder 6', fontsize=10)\n    plt.tick_params(axis='x', labelsize=8)\n    plt.tick_params(axis='y', labelsize=8)\n    plt.title(f\"Responder_6 vs {feature}\", fontsize=12)\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.377766Z","iopub.status.idle":"2025-01-06T16:24:11.378346Z","shell.execute_reply":"2025-01-06T16:24:11.378184Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Statistical Tests","metadata":{}},{"cell_type":"markdown","source":"**Neyman-Pearson Test**\n\nThis example illustrates a Neyman-Pearson test for correlation in the classical theoretical sense. This involves several strong/false assumptions and limitations that differ from the practical implementation.\n\n- Simple Hypotheses\n\nThe Neyman-Pearson lemma applies only to simple vs. simple hypotheses. We must specify:\n$$ H_0: \\rho = \\rho_1 \\text{ and } H_1: \\rho = \\rho_1 $$\n\nwhere $\\rho_0$ and $\\rho_1$ are both known, fixed constants, and $\\rho_0 \\neq \\rho_1$. For example, we might choose $\\rho_0 = 0$ and $\\rho_1 = 0.1$.\n\n- Model Assumptions ( totally biased and in opposition to the previous analysis )\n\nAssume $(X_i, Y_i)$ are i.i.d. from a bivariate normal distribution.\nUnder $H_0$: $\\rho = 0$, so $(X_i, Y_i)$ are independent standard normals. Under $H_1$: $\\rho = \\rho_1$, known and fixed.\n\n- Neyman-Pearson Lemma and Likelihood Ratio Test\n\nThe Neyman-Pearson lemma states that the most powerful test for these two simple hypotheses is a likelihood ratio test:\n\n$$ \\Lambda = \\frac{L(X_1,\\ldots,X_n,Y_1,\\ldots,Y_n|\\rho=\\rho_1)}{L(X_1,\\ldots,X_n,Y_1,\\ldots,Y_n|\\rho=0)}$$\n\nReject $H_0$ if $\\Lambda > k$ for some critical value $k$ determined by the size $\\alpha$.\n\nFor a bivariate normal model:\nUnder $H_0$, the joint density of each pair is:\n       $\n            f_0(x,y) = \\frac{1}{2\\pi}e^{-(x^2+y^2)/2}\n        $\n    Under $H_1$, the joint density is:\n    $$ f_1(x,y) = \\frac{1}{2\\pi\\sqrt{1-\\rho_1^2}}\\exp\\left(-\\frac{x^2-2\\rho_1xy+y^2}{2(1-\\rho_1^2)}\\right)$$\n\n\nThe likelihood ratio for the entire sample is:\n$$\\Lambda = \\left(\\frac{1}{1-\\rho_1^2}\\right)^{n/2}\\exp\\left(-\\frac{\\sum(x_i^2-2\\rho_1x_iy_i+y_i^2)}{2(1-\\rho_1^2)}+\\frac{\\sum(x_i^2+y_i^2)}{2}\\right)$$\n\nThis simplifies and can be expressed in terms of the sample correlation:\n$r = \\frac{\\sum x_iy_i}{\\sqrt{\\sum x_i^2\\sum y_i^2}}$\n\n- Determining the Critical Value\n\nTo have a size $\\alpha$ test, we solve:\n$ P_{H_0}(\\Lambda > k) = \\alpha $\n\nThis can be done using Fisher's z-transform:\n$$Z = \\frac{1}{2}\\ln\\left(\\frac{1+r}{1-r}\\right) \\sim N\\left(0,\\frac{1}{n-3}\\right)$$\n\nunder $H_0$. For a chosen $\\alpha$, find the critical z-value $z_\\alpha$ and translate that back to a threshold on $r$. Under $H_1$, the distribution of $Z$ shifts, ensuring this test is the most powerful for detecting $\\rho_1$ specifically.","metadata":{}},{"cell_type":"markdown","source":"Based on previous correlation analysis, i will only experiment on 10 best ","metadata":{}},{"cell_type":"code","source":"top_10_features.index.tolist()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.378957Z","iopub.status.idle":"2025-01-06T16:24:11.379626Z","shell.execute_reply":"2025-01-06T16:24:11.37945Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"With so much data, the test finds even tiny positive correlations significant, while negative ones do not pass the threshold. Since the sample is huge, the test is incredibly sensitive and has essentially perfect power. However, this also means that even very small effects appear significant, so the practical importance of these correlations might be questionable to not say useless.","metadata":{}},{"cell_type":"code","source":"alpha = 0.05\nrho_1 = 0.1\nresponder_name = \"responder_6\"\n\nresponders_to_test = top_10_features.index.tolist()\n\nresults_list = []\n\nfor feature_name in responders_to_test:\n    df_pair = train[[responder_name, feature_name]].dropna()\n    X = df_pair[feature_name].values\n    Y = df_pair[responder_name].values\n    n = len(X)\n\n    r = np.corrcoef(X, Y)[0,1]\n    Z_obs = 0.5 * log((1+r)/(1-r))\n    z_alpha_standard = norm.ppf(1 - alpha)\n    z_alpha = z_alpha_standard * sqrt(1/(n-3))\n    reject_H0 = (Z_obs > z_alpha)\n    mean_Z_h1 = 0.5 * log((1+rho_1)/(1-rho_1))\n    power = 1 - norm.cdf(z_alpha, loc=mean_Z_h1, scale=sqrt(1/(n-3)))\n    results_list.append({\n        'Feature': feature_name,\n        'N': n,\n        'Correlation': r,\n        'Z_Obs': Z_obs,\n        'Z_Critical': z_alpha,\n        'Reject_H0': reject_H0,\n        'Power_at_rho1': power\n    })\n\nnp_results = pd.DataFrame(results_list)\nnp_results","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.380176Z","iopub.status.idle":"2025-01-06T16:24:11.381075Z","shell.execute_reply":"2025-01-06T16:24:11.380892Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Wald Test**\n\nThe Wald test assesses whether an estimated parameter $\\hat{\\theta}$ differs significantly from a hypothesized value $\\theta_0$.\n\nSuppose we have:\n\\begin{align*}\n   H_0&: \\theta = \\theta_0 \\text{ vs. } H_1: \\theta \\neq \\theta_0\n\\end{align*}\n\nWe estimate $\\theta$ from the data to get $\\hat{\\theta}$ and its standard error $SE(\\hat{\\theta})$.\n\nThe Wald statistic is:\n\\begin{equation}\n   W_n = \\frac{(\\hat{\\theta}_n - \\theta_0)^2}{V(\\hat{\\theta}_n)}\n\\end{equation}\n\nUnder $H_0$, $W_n$ converges in distribution to a $\\chi^2(1)$.\nThe Wald test thus follows an asymptotic $\\chi^2(1)$ distribution under $H_0$.","metadata":{}},{"cell_type":"markdown","source":"The Wald test results show that all features have statistically significant relationships with responder_6, as all p-values are 0, and Hypothesis0 is rejected for every feature. Features like responder_3, responder_8, and responder_7 exhibit particularly strong effects as before, and as reflected by their large slopes and high Wald statistics. As of features_16 and features_17, who seems to be most significant features.","metadata":{}},{"cell_type":"code","source":"responder_name = \"responder_6\"\nresponders_to_test = top_10_features.index.tolist()\n\nresults_list = []\n\nfor feature_name in responders_to_test:\n    df_pair = train[[responder_name, feature_name]].dropna()\n    X = df_pair[feature_name].values\n    Y = df_pair[responder_name].values\n    n = len(X)\n    \n    X_with_const = sm.add_constant(X)\n    model = sm.OLS(Y, X_with_const).fit()\n    \n    slope = model.params[1]\n    se = model.bse[1]\n    wald_stat = (slope ** 2) / (se ** 2)  \n    p_val = 1 - chi2.cdf(wald_stat, df=1)  \n    \n    alpha = 0.05\n    reject_H0 = p_val < alpha\n\n    results_list.append({\n        'Feature': feature_name,\n        'N': n,\n        'Slope': slope,\n        'SE': se,\n        'Wald_Stat': wald_stat,\n        'P_Value': p_val,\n        'Reject_H0': reject_H0\n    })\n\nwald_results = pd.DataFrame(results_list)\nwald_results","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.38186Z","iopub.status.idle":"2025-01-06T16:24:11.382447Z","shell.execute_reply":"2025-01-06T16:24:11.382282Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Permutation Importance**\n\nPermutation importance is a model-agnostic technique to assess the significance of features in a predictive model. It measures how much the model's performance deteriorates when the values of a particular feature are randomly shuffled, breaking its relationship with the target variable.\n\n- Compute Baseline Error\nEvaluate the model's performance using the unaltered dataset to get the baseline error:\n\\begin{equation}\n   \\text{Baseline Error} = \\text{Metric}(y, \\hat{y})\n\\end{equation}\nwhere $y$ are the true values and $\\hat{y}$ are the model's predictions.\n\n- Permute Feature Values\nFor each feature $X_j$, shuffle its values randomly to create a modified dataset $\\tilde{X}_j$. This destroys any relationship between $X_j$ and the target $y$.\n\n- Compute Permutation Error\nEvaluate the model's performance using the dataset with permuted $X_j$:\n\\begin{equation}\n   \\text{Permuted Error}_j = \\text{Metric}(y, \\hat{y}^{\\text{perm}})\n\\end{equation}\nwhere $\\hat{y}^{\\text{perm}}$ are the model's predictions on the permuted dataset.\n\n- Calculate Importance\nThe importance of feature $X_j$ is quantified by the increase in error:\n\\begin{equation}\n   \\Delta_j = \\text{Permuted Error}_j - \\text{Baseline Error}\n\\end{equation}\n\n- Repeat for Stability\nRepeat the permutation process multiple times to estimate the mean and standard deviation of $\\Delta_j$.","metadata":{}},{"cell_type":"markdown","source":"The permutation importance analysis shows that responder_3 is the most critical feature, with a large impact on model performance $(ΔMSE=1.45)$. responder_5, responder_8, and responder_7 also significantly contribute, while features like feature_16 and feature_17 have minimal impact $(ΔMSE≈0)$. The consistency between both results confirms the robustness of these findings.","metadata":{}},{"cell_type":"code","source":"X = train[top_10_features.index.tolist()]\ny = train['responder_6']\n\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n\nlgb_model = lgb.LGBMRegressor(random_state=42)\nlgb_model.fit(X_train, y_train)\n\nbaseline_mse = mean_squared_error(y_test, lgb_model.predict(X_test))\nprint(f\"Baseline MSE: {baseline_mse}\")\n\nperm_importance = permutation_importance(lgb_model, X_test, y_test, n_repeats=10, random_state=42)\nperm_importance_df = pd.DataFrame({\n    'Feature': X_test.columns,\n    'Importance': perm_importance.importances_mean,\n    'Std_Dev': perm_importance.importances_std\n}).sort_values(by='Importance', ascending=False)\nprint(\"Permutation Importance sklearn:\")\nprint(perm_importance_df)\n\nresults_list = []\nfor feature_name in top_10_features.index.tolist():\n\n    X_test_permuted = X_test.copy()\n    X_test_permuted[feature_name] = np.random.permutation(X_test_permuted[feature_name])\n    \n    y_pred_permuted = lgb_model.predict(X_test_permuted)\n    permuted_mse = mean_squared_error(y_test, y_pred_permuted)\n    delta_mse = permuted_mse - baseline_mse\n    \n    results_list.append({\n        'Feature': feature_name,\n        'Baseline_MSE': baseline_mse,\n        'Permuted_MSE': permuted_mse,\n        'Delta_MSE': delta_mse\n    })\n\ncustom_perm_df = pd.DataFrame(results_list).sort_values(by='Delta_MSE', ascending=False)\nprint(\"\\nCustom Permutation Importance:\")\nprint(custom_perm_df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:24:11.383262Z","iopub.status.idle":"2025-01-06T16:24:11.383638Z","shell.execute_reply":"2025-01-06T16:24:11.383469Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Prediction/ Conclusion\n\nIn this notebook, we explored the relationships between features and responder_6 using different methods. Features like responder_3, responder_8, and feature_36, feature_16 and feature_17 were identified as highly correlated with responder_6, confirmed through both Pearson and Spearman correlations. Statistical tests (which have limitations due to incorrect assumptions being made), including Neyman-Pearson and permutation importance, validated these findings and demonstrated the significance of the selected features. Feature engineering, such as interaction terms and time-based sinusoidal features, further enhanced the model's performance. A LightGBM model successfully predicted responder_6, confirming the relevance of the chosen features. Overall, we addressed the initial hypothesis by identifying significant relationships between certain features and responder_6 and validating them with statistical methods.\n\n\nLet's perform predictions on the responder_6 target variable using a LGBM model. It incorporates feature engineering techniques and important features identified earlier.","metadata":{}},{"cell_type":"code","source":"%%time\n\nresponders = pd.read_csv('/kaggle/input/jane-street-real-time-market-data-forecasting/responders.csv')\nfeatures = pd.read_csv('/kaggle/input/jane-street-real-time-market-data-forecasting/features.csv')\nsample_submission = pd.read_csv('/kaggle/input/jane-street-real-time-market-data-forecasting/sample_submission.csv')\n\ntrain = (\n    pl.read_parquet('/kaggle/input/jane-street-real-time-market-data-forecasting/train.parquet/partition_id=9/part-0.parquet')\n).to_pandas()\n\nlgb_params = {\n    \"boosting_type\": \"gbdt\",\n    \"metric\": \"rmse\",\n    \"random_state\": 9,\n    \"learning_rate\": 0.05,\n    \"n_estimators\": 1000,\n    \"num_leaves\": 64,\n    \"colsample_bytree\": 0.8,\n    \"subsample\": 0.8,\n    \"reg_alpha\": 0.1,\n    \"reg_lambda\": 10,\n    \"min_child_weight\": 10,\n    \"device\": \"gpu\", \n}\n\nimportant_features = [\n    'symbol_id', 'feature_16', 'feature_17', 'feature_36', \n    'responder_3', 'responder_7', 'responder_8'\n]\n\ntrain['feature_16_17_interaction'] = train['feature_16'] * train['feature_17']\ntrain['responder_avg'] = (train['responder_3'] + train['responder_7'] + train['responder_8']) / 3\ntrain['sin_time_id'] = np.sin(2 * np.pi * train['time_id'] / 967)\ntrain['cos_time_id'] = np.cos(2 * np.pi * train['time_id'] / 967)\ntrain['feature_36_squared'] = train['feature_36'] ** 2\ntrain['feature_16_36_interaction'] = train['feature_16'] * train['feature_36']\ntrain['feature_ratio'] = train['feature_16'] / (train['feature_17'] + 1e-9)\ntrain['responder_sum'] = train[['responder_3', 'responder_7', 'responder_8']].sum(axis=1)\ntrain['feature_16_rolling_mean'] = train['feature_16'].rolling(window=5, min_periods=1).mean()\ntrain['feature_16_rolling_std'] = train['feature_16'].rolling(window=5, min_periods=1).std()\n\nfinal_feature = important_features + [\n    'feature_16_17_interaction', 'responder_avg', 'sin_time_id', 'cos_time_id',\n    'feature_36_squared', 'feature_16_36_interaction', 'feature_ratio', 'responder_sum',\n    'feature_16_rolling_mean', 'feature_16_rolling_std'\n]\ntrain = train[['responder_6'] + final_feature]\n\nlgb_model = LGBMRegressor(**lgb_params)\nlgb_model.fit(train[final_feature], train['responder_6'])\n\ndef predict(test: pl.DataFrame, lags: pl.DataFrame | None) -> pl.DataFrame:\n    global lags_\n    if lags is not None:\n        lags_ = lags\n\n    test_df = test.to_pandas()\n\n    for col in important_features + ['time_id']:\n        if col not in test_df.columns:\n            test_df[col] = 0\n\n    test_df['feature_16_17_interaction'] = test_df['feature_16'] * test_df['feature_17']\n    test_df['responder_avg'] = (test_df['responder_3'] + test_df['responder_7'] + test_df['responder_8']) / 3\n    test_df['sin_time_id'] = np.sin(2 * np.pi * test_df['time_id'] / 967)\n    test_df['cos_time_id'] = np.cos(2 * np.pi * test_df['time_id'] / 967)\n    test_df['feature_36_squared'] = test_df['feature_36'] ** 2\n    test_df['feature_16_36_interaction'] = test_df['feature_16'] * test_df['feature_36']\n    test_df['feature_ratio'] = test_df['feature_16'] / (test_df['feature_17'] + 1e-9)\n    test_df['responder_sum'] = test_df[['responder_3', 'responder_7', 'responder_8']].sum(axis=1)\n    test_df['feature_16_rolling_mean'] = test_df['feature_16'].rolling(window=5, min_periods=1).mean()\n    test_df['feature_16_rolling_std'] = test_df['feature_16'].rolling(window=5, min_periods=1).std()\n\n    test_df = test_df[final_feature].fillna(-1)\n    predictions = lgb_model.predict(test_df)\n    eps = 1e-10\n    predictions = np.clip(predictions, -5+eps, 5-eps)\n\n    predictions_df = test.select(\n        'row_id',\n        pl.Series('responder_6', predictions),\n    )\n\n    assert isinstance(predictions_df, (pl.DataFrame, pd.DataFrame)), \"Predictions must be a Polars or Pandas DataFrame.\"\n    print(\"Predictions are a Polars or Pandas DataFrame.\")\n\n    assert list(predictions_df.columns) == ['row_id', 'responder_6'], \"Predictions must have columns ['row_id', 'responder_6'].\"\n    print(\"Predictions have the correct columns ['row_id', 'responder_6'].\")\n\n    assert len(predictions_df) == len(test), \"Predictions must have the same number of rows as the test data.\"\n    print(\"Predictions have the same number of rows as the test data.\")\n\n    return predictions_df\n\ninference_server = kaggle_evaluation.jane_street_inference_server.JSInferenceServer(predict)\n\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    inference_server.serve()\nelse:\n    inference_server.run_local_gateway(\n        (\n            '/kaggle/input/jane-street-real-time-market-data-forecasting/test.parquet',\n            '/kaggle/input/jane-street-real-time-market-data-forecasting/lags.parquet',\n        )\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:30:05.796254Z","iopub.execute_input":"2025-01-06T16:30:05.796599Z","iopub.status.idle":"2025-01-06T16:31:31.069117Z","shell.execute_reply.started":"2025-01-06T16:30:05.796572Z","shell.execute_reply":"2025-01-06T16:31:31.068236Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('EOF')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-06T16:31:35.990388Z","iopub.execute_input":"2025-01-06T16:31:35.990693Z","iopub.status.idle":"2025-01-06T16:31:35.99509Z","shell.execute_reply.started":"2025-01-06T16:31:35.99067Z","shell.execute_reply":"2025-01-06T16:31:35.994228Z"}},"outputs":[],"execution_count":null}]}