{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":84493,"databundleVersionId":9871156,"sourceType":"competition"}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport polars as pl\nfrom matplotlib import pyplot as plt\nimport seaborn as sns\nfrom pathlib import Path\nfrom tqdm import tqdm\nfrom functools import reduce\n\nfrom scipy.stats import pearsonr, spearmanr, kendalltau\n\nDATA_DIR = Path('/kaggle/input/jane-street-real-time-market-data-forecasting')\nN_PARTITION = 10\n\nfeature_cols = [f'feature_{x:02}' for x in range(79)]\nresponder_cols = [f'responder_{i}' for i in range(9)]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-11-16T23:52:43.115868Z","iopub.execute_input":"2024-11-16T23:52:43.116260Z","iopub.status.idle":"2024-11-16T23:52:46.261039Z","shell.execute_reply.started":"2024-11-16T23:52:43.116221Z","shell.execute_reply":"2024-11-16T23:52:46.259897Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Prepare data and have some basic EDA","metadata":{}},{"cell_type":"code","source":"train_parquets = [DATA_DIR / f\"train.parquet/partition_id={i}/part-0.parquet\" for i in range(N_PARTITION)]\n\npartition_dates = {}\nfor _f in train_parquets:\n    partition_id = int(_f.parents[0].stem.split('=')[-1])\n    _pl = pl.read_parquet(_f, columns=[\"date_id\", \"time_id\"])\n    partition_dates[partition_id] = (_pl['date_id'].min(), _pl['date_id'].max())\n\n# check date range in each partition\npartition_dates","metadata":{"execution":{"iopub.status.busy":"2024-11-16T23:52:46.262730Z","iopub.execute_input":"2024-11-16T23:52:46.263148Z","iopub.status.idle":"2024-11-16T23:52:46.839492Z","shell.execute_reply.started":"2024-11-16T23:52:46.263115Z","shell.execute_reply":"2024-11-16T23:52:46.838434Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pl_train = pl.read_parquet(train_parquets, columns=[\"date_id\", \"time_id\"])\n\ntime_id_per_day = pl_train.group_by(\"date_id\").agg([pl.col(\"time_id\").n_unique().alias(\"time_count\")])\ntime_id_per_day = time_id_per_day.sort('date_id')\n\nprint(\"Number of dates: \",  time_id_per_day.select('date_id').n_unique())\n\nprint(time_id_per_day.filter(pl.col('time_count')<900)['date_id'].max()) # before date_id 676, each day had 849 timestamps\nprint(time_id_per_day.filter(pl.col('time_count')>900)['date_id'].min()) # after date_id 677, each day had 968 timestamps","metadata":{"execution":{"iopub.status.busy":"2024-10-17T22:50:50.505062Z","iopub.execute_input":"2024-10-17T22:50:50.505907Z","iopub.status.idle":"2024-10-17T22:50:51.347505Z","shell.execute_reply.started":"2024-10-17T22:50:50.505864Z","shell.execute_reply":"2024-10-17T22:50:51.346446Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Check null values and trading symbols","metadata":{}},{"cell_type":"code","source":"# check distribution of Null values\nnull_counts_collect = []\n\nfor _col in tqdm(feature_cols):\n    pl_train = pl.read_parquet(train_parquets, columns=[\"date_id\", \"time_id\", \"symbol_id\", _col])\n    null_counts = pl_train.group_by([\"date_id\", \"symbol_id\"]).agg([pl.col(_col).is_null().sum().alias(f\"null_count_{_col}\")]).sort(['date_id', 'symbol_id'])\n    null_counts_collect.append(\n        null_counts.group_by('date_id').agg([\n            pl.col(f'null_count_{_col}').mean().alias(f\"mean_null_count_{_col}\"),\n            pl.col(f'null_count_{_col}').max().alias(f\"max_null_count_{_col}\"),\n            pl.col(f'null_count_{_col}').min().alias(f\"min_null_count_{_col}\"),\n            pl.col(f'null_count_{_col}').n_unique().alias(f\"nunique_null_count_{_col}\"),\n        ])\n    )\n    \npl_null_counts = reduce(lambda l, r: l.join(r, on=\"date_id\", how=\"inner\"), null_counts_collect)\n\npl_null_counts","metadata":{"execution":{"iopub.status.busy":"2024-10-17T22:52:12.159950Z","iopub.execute_input":"2024-10-17T22:52:12.160697Z","iopub.status.idle":"2024-10-17T22:52:25.504044Z","shell.execute_reply.started":"2024-10-17T22:52:12.160655Z","shell.execute_reply":"2024-10-17T22:52:25.502192Z"},"trusted":true,"collapsed":true,"jupyter":{"outputs_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot the null value distribution for features & dates \n# notice some features did not exist in the beginning dates \n\nmean_counts_col = [f'mean_null_count_{col}' for col in feature_cols]\nmean_counts_arr = pl_null_counts.select(mean_counts_col).to_numpy().T\n\nim = plt.pcolormesh(mean_counts_arr, cmap='Reds')\nplt.colorbar(im, label='Mean Null Counts')\nplt.xlabel('Date ID')\nplt.ylabel('Feature ID')\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot variation of traded symbols on each date. The early dates had less symbols than later. \n\npl_train = pl.concat([pl.read_parquet(_f, columns=[\"date_id\", \"symbol_id\"]) for _f in train_parquets])\nsymbol_count = pl_train.group_by(\"date_id\").agg([pl.col(\"symbol_id\").n_unique().alias('n_symbols')]).sort('date_id').to_pandas()\n\nfig, ax = plt.subplots(figsize=(10, 6))\nax.step(symbol_count['date_id'], symbol_count['n_symbols'])\nax.set_xlabel('Date ID')\nax.set_ylabel('Number of Symbols')\nax.grid(True, ls=\"--\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-17T22:52:30.313101Z","iopub.execute_input":"2024-10-17T22:52:30.313547Z","iopub.status.idle":"2024-10-17T22:52:31.490882Z","shell.execute_reply.started":"2024-10-17T22:52:30.313510Z","shell.execute_reply":"2024-10-17T22:52:31.489734Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot the start & end date of each symbol being traded\n\ndate_count = pl_train.group_by('symbol_id').agg([pl.col('date_id').min().alias('min_date'), pl.col('date_id').max().alias('max_date')]).sort('symbol_id')\n\nplt.figure(figsize=(10, 6))\nplt.barh(date_count['symbol_id'], date_count['max_date'], color='black', label='Traded')\nplt.barh(date_count['symbol_id'], date_count['min_date'], color='lightgray', label='Not Traded')\nplt.grid(axis='x')\nplt.legend(loc='center left', bbox_to_anchor=[1, 0.5])\nplt.xlabel('Date ID')\nplt.ylabel('Symbol ID')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-17T22:53:54.359596Z","iopub.execute_input":"2024-10-17T22:53:54.360024Z","iopub.status.idle":"2024-10-17T22:53:55.475601Z","shell.execute_reply.started":"2024-10-17T22:53:54.359984Z","shell.execute_reply":"2024-10-17T22:53:55.474499Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Weights ~ Responders\n\nIt seems that higher weights is related to lower dispersion of responders, especially for responder_6. Assume weight is the weight of instrument in a portfolio and responder_6 is some kind of return, this implies that a higher weight is applied if the expected return is more stable.","metadata":{}},{"cell_type":"code","source":"from scipy.stats import laplace\n\npl_train = pl.read_parquet(train_parquets, columns=[\"date_id\", \"time_id\", \"symbol_id\", 'weight']+responder_cols)\n\nhist, bin_edges = np.histogram(pl_train['weight'].to_numpy(), bins=100)\n\ncut_labels = (bin_edges[:-1] + bin_edges[1:]) / 2\ncut_labels = [\"0\"] + [str(x) for x in cut_labels] + [str(bin_edges[-1])]\n\npl_train = pl_train.with_columns(pl.col('weight').cut(bin_edges, labels=cut_labels).alias('weight_cut'))\n\nfor resp_col in responder_cols:\n\n    N, bins = np.histogram(pl_train[resp_col].to_numpy(), bins=100, density=True)\n    loc, scale = laplace.fit(pl_train[resp_col].to_numpy())\n    x_vals = (bins[:-1] + bins[1:]) / 2\n\n    resp_bin = pl_train.group_by('weight_cut').agg([\n        pl.col(resp_col).mean().alias(f'{resp_col}_mean'),\n        pl.col(resp_col).std().alias(f'{resp_col}_std'),\n        pl.col(resp_col).count().alias(f'{resp_col}_count'),\n        ((pl.col(resp_col) - loc)/scale).abs().mean().alias(f'{resp_col}_zscore')\n    ]).sort('weight_cut')\n    \n    x_ticks = [float(x) for x in resp_bin['weight_cut']]\n\n    fig, axes = plt.subplots(1, 4, figsize=(20, 4))\n    axes[0].scatter(x_ticks, resp_bin[f'{resp_col}_mean'] )\n    axes[0].set_ylim(-0.05, 0.05)\n    axes[0].set_title(f'{resp_col}_mean')\n    axes[0].set_ylabel(f'{resp_col}_mean')\n\n    axes[1].scatter(x_ticks, resp_bin[f'{resp_col}_std'] )\n    axes[1].set_title(f'{resp_col}_std')\n    axes[1].set_ylabel(f'{resp_col}_std')\n\n    axes[2].scatter(x_ticks, resp_bin[f'{resp_col}_zscore'] )\n    axes[2].set_title(f'{resp_col}_zscore')\n    axes[2].set_ylabel(f'{resp_col}_zscore')\n\n    axes[3].plot(x_vals, laplace.pdf(x_vals, loc=loc, scale=scale), 'r-', lw=2, alpha=0.6, label='laplace pdf')\n    axes[3].bar(x_vals, N, alpha=0.3, width=x_vals[2]-x_vals[1], color='b', label='Histogram')\n    axes[3].legend()\n    axes[3].set_title(f\"{resp_col} | loc={loc:.4f}; scale={scale:.4f}\")\n    axes[3].grid(True, ls=\"--\")\n    axes[3].set_xlabel(resp_col)\n    axes[3].set_ylabel('Density')\n\n    for ax in axes[:-1]:\n        ax.grid(True, ls=\"--\")\n        ax.set_xlabel(\"weights\")\n        ax_r = ax.twinx()\n        ax_r.bar(x_ticks, resp_bin[f'{resp_col}_count'], alpha=0.3, width=x_ticks[2]-x_ticks[1], color='k', label='Count')\n        ax_r.legend()\n        ax.set_xlim(0.1, 6)\n\n    fig.tight_layout()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-17T22:54:37.421679Z","iopub.execute_input":"2024-10-17T22:54:37.422165Z","iopub.status.idle":"2024-10-17T22:55:27.470277Z","shell.execute_reply.started":"2024-10-17T22:54:37.422123Z","shell.execute_reply":"2024-10-17T22:55:27.468080Z"},"trusted":true,"collapsed":true,"jupyter":{"outputs_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Responders day-to-day change","metadata":{}},{"cell_type":"code","source":"resp_col = 'responder_6'\n\ncheck_symbol = 38\n\ndaily_resp_group = pl_train.group_by(['date_id', 'symbol_id'])\n\ndaily_resp_on_symbol = [\n    (_date, batch.select('date_id', 'time_id', resp_col)) \n    for (_date, _symbol), batch in daily_resp_group if _symbol == check_symbol\n]\n\ndaily_resp_on_symbol = sorted(daily_resp_on_symbol, key=lambda x: x[0])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-17T23:01:43.037926Z","iopub.execute_input":"2024-10-17T23:01:43.038422Z","iopub.status.idle":"2024-10-17T23:01:55.219769Z","shell.execute_reply.started":"2024-10-17T23:01:43.038385Z","shell.execute_reply":"2024-10-17T23:01:55.218372Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, axes = plt.subplots(5, 1, figsize=(12, 15), sharex=True)\n\noffset = 1000\n\nfor i, ax in enumerate(axes.flatten()):\n    time_ids = daily_resp_on_symbol[i+offset][1]['time_id']\n    ts_resp = daily_resp_on_symbol[i+offset][1][resp_col]\n    ax.plot(time_ids[::9], ts_resp[::9])\n    ax.set_title(f\"Date ID = {daily_resp_on_symbol[i+offset][0]}\")\n    ax.grid(True, ls=\"--\")\n    \nfig.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-17T23:45:11.391193Z","iopub.execute_input":"2024-10-17T23:45:11.391632Z","iopub.status.idle":"2024-10-17T23:45:12.842184Z","shell.execute_reply.started":"2024-10-17T23:45:11.391596Z","shell.execute_reply":"2024-10-17T23:45:12.841019Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The plot below shows the daily change of the std of responder_6. Clearly, we can see the high volatility days and low volatility days. However, we can not expect the pattern of the time-series would look similar from day to day. ","metadata":{}},{"cell_type":"code","source":"resp_std = [x[1][resp_col].std() for x in daily_resp_on_symbol]\n\nresp_std_mean = pl.DataFrame({f'{resp_col}_std': resp_std}).with_columns(\n        pl.col(f'{resp_col}_std').rolling_mean(window_size=5).alias(\"r_mean\")\n    ).select('r_mean')\n\nfig = plt.figure(figsize=(15, 5))\nplt.plot(resp_std, label='std')\nplt.plot(resp_std_mean, label='rolling_mean')\nplt.xlabel('Date')\nplt.ylabel(f'{resp_col} std')\nplt.grid(ls=\"--\")\nplt.legend()\n# plt.xlim(1000, pl_train['date_id'].max())\nplt.xlim(1000, 1200)\nplt.show() ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-17T23:15:33.433470Z","iopub.execute_input":"2024-10-17T23:15:33.433894Z","iopub.status.idle":"2024-10-17T23:15:33.761089Z","shell.execute_reply.started":"2024-10-17T23:15:33.433858Z","shell.execute_reply":"2024-10-17T23:15:33.759940Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Let's use DTW to compare the difference between time-series of two days. Not really showing any pattern here.","metadata":{}},{"cell_type":"code","source":"# ! pip -q install dtw","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-17T23:52:02.374843Z","iopub.execute_input":"2024-10-17T23:52:02.375354Z","iopub.status.idle":"2024-10-17T23:52:14.442439Z","shell.execute_reply.started":"2024-10-17T23:52:02.375311Z","shell.execute_reply":"2024-10-17T23:52:14.441009Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# from dtw import *\n\n# offset = 1000\n# end = 1050\n# date_length = end - offset\n\n# dtw_dist_mat = np.zeros((date_length, date_length))\n\n# for i in range(date_length):\n#     for j in range(i+1, date_length):\n        \n#         query = daily_resp_on_symbol[j+offset][1][resp_col]\n#         template = daily_resp_on_symbol[i+offset][1][resp_col]\n        \n#         alignment = dtw(query[::9], template[::9], keep_internals=True)\n\n#         dtw_dist_mat[i, j] = alignment.distance\n#         dtw_dist_mat[j, i] = alignment.distance\n\n# sns.heatmap(dtw_dist_mat, cmap='viridis')\n# plt.title('DTW Distance Matrix')\n# plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-17T23:52:14.445546Z","iopub.execute_input":"2024-10-17T23:52:14.446089Z","iopub.status.idle":"2024-10-17T23:52:17.586471Z","shell.execute_reply.started":"2024-10-17T23:52:14.446036Z","shell.execute_reply":"2024-10-17T23:52:17.585061Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Maybe responders of the same day are correlated?","metadata":{}},{"cell_type":"code","source":"# plot responders of different symbols on the same day\npl_filter = pl_train.filter(pl.col('date_id')==500)\nunique_symbols = np.unique(np.vstack(pl_filter.group_by('time_id').agg(pl.col('symbol_id').unique())['symbol_id'].to_numpy()))\n\nfor col in responder_cols:\n    fig, ax = plt.subplots(1, 1, figsize=(14, 3))\n    for _id in unique_symbols:\n        ax.plot(pl_filter.filter(pl.col('symbol_id')==_id)[col], label=_id)\n    ax.set_title(col)\n#     ax.legend(loc='center left', bbox_to_anchor=[1, 0.5])\n    ax.grid(True, ls=\"--\")\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-17T23:52:39.226015Z","iopub.execute_input":"2024-10-17T23:52:39.226510Z","iopub.status.idle":"2024-10-17T23:52:42.417771Z","shell.execute_reply.started":"2024-10-17T23:52:39.226468Z","shell.execute_reply":"2024-10-17T23:52:42.416559Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot correlation of responders on the same day between symbols \nfor col in responder_cols:\n\n    corr_mat = np.ones((len(unique_symbols), len(unique_symbols)))\n    pval_mat = np.ones((len(unique_symbols), len(unique_symbols)))\n\n    for i in range(len(unique_symbols)):\n        for j in range(i+1, len(unique_symbols)):\n            ts_1 = pl_filter.filter(pl.col('symbol_id')==unique_symbols[i])[col]\n            ts_2 = pl_filter.filter(pl.col('symbol_id')==unique_symbols[j])[col]\n\n            corr, pval = spearmanr(ts_1, ts_2)\n\n            corr_mat[i, j] = corr; corr_mat[j, i] = corr\n            pval_mat[i, j] = pval; pval_mat[j, i] = pval\n\n\n    fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n\n    tick_kwargs = dict(xticklabels=unique_symbols, yticklabels=unique_symbols, annot=False, fmt=\".2f\", linewidth=.5,)\n\n    sns.heatmap(corr_mat, cmap='coolwarm', vmax=1, vmin=-1,  ax=axes[0], **tick_kwargs)\n    sns.heatmap(pval_mat, vmin=0, vmax=0.1, ax=axes[1], **tick_kwargs)\n\n    axes[0].set_title('Correlation')\n    axes[1].set_title('P-value')\n\n    fig.suptitle(col)\n    fig.tight_layout()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-17T23:54:16.465558Z","iopub.execute_input":"2024-10-17T23:54:16.466087Z","iopub.status.idle":"2024-10-17T23:54:43.020722Z","shell.execute_reply.started":"2024-10-17T23:54:16.466044Z","shell.execute_reply":"2024-10-17T23:54:43.019321Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Fix a specific date to plot intraday timeseries\n\ndate_choice = 1400\n\nfor k, (start, end) in partition_dates.items():\n    if start < date_choice < end:\n        load_partition = k\n        break\n\npl_train = pl.read_parquet(train_parquets[load_partition])\n\npl_train","metadata":{"execution":{"iopub.status.busy":"2024-11-16T23:53:06.928887Z","iopub.execute_input":"2024-11-16T23:53:06.929282Z","iopub.status.idle":"2024-11-16T23:53:12.682010Z","shell.execute_reply.started":"2024-11-16T23:53:06.929245Z","shell.execute_reply":"2024-11-16T23:53:12.680940Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot intraday ts of a symbol\n\nsymbol_choice = 37\n\npl_filter = pl_train.filter((pl.col('date_id')==date_choice) & (pl.col(\"symbol_id\")==symbol_choice))\ntime_series = pl_filter.select(\"time_id\").to_series()\nresponders = [pl_filter.select(r_col).to_series() for r_col in responder_cols]\n\nfig, axes = plt.subplots(3, 3, figsize=(14, 9), sharex=True)\nfor i, (ax, _responder) in enumerate(zip(axes.flatten(), responders)):\n    ax.plot(time_series, _responder)\n    ax.set_title(responder_cols[i])\n    ax.grid(True, linestyle=\"--\")\n    \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T23:53:19.155394Z","iopub.execute_input":"2024-11-16T23:53:19.155808Z","iopub.status.idle":"2024-11-16T23:53:21.009861Z","shell.execute_reply.started":"2024-11-16T23:53:19.155768Z","shell.execute_reply":"2024-11-16T23:53:21.008576Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, axes = plt.subplots(1, 5, figsize=(18, 3), sharey=True)\n\noffset = np.arange(-20, 30, 1)\n\nfor i, (ax, _off) in enumerate(zip(axes.flatten(), offset)):\n    _filter = pl_train.filter((pl.col('date_id')==date_choice+_off) & (pl.col(\"symbol_id\")==symbol_choice))\n    ax.plot(_filter.select(\"time_id\").to_series(), _filter.select(\"responder_6\").to_series() )\n    ax.set_title(f\"Date {date_choice+_off} | responder_6\")\n    ax.grid(True, linestyle=\"--\")\n\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T23:53:44.349890Z","iopub.execute_input":"2024-11-16T23:53:44.350298Z","iopub.status.idle":"2024-11-16T23:53:45.724776Z","shell.execute_reply.started":"2024-11-16T23:53:44.350263Z","shell.execute_reply":"2024-11-16T23:53:45.723769Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot intraday feature ts of the chosen symbol\nexclude_features = [f\"feature_{i:02d}\" for i in [9,10,11]+list(range(20, 32))]\nfeature_series = [pl_filter.select(col).to_series() for col in feature_cols if col not in exclude_features]\n\nfig, axes = plt.subplots(len(feature_series), 1, figsize=(15, 2*len(feature_series)), sharex=True)\nfor i, (ax, _series) in enumerate(zip(axes.flatten(), feature_series)):\n    ax.plot(time_series, _series)\n    ax.plot(time_series, _series.rolling_mean(5), color='orange')\n    ax.set_title(feature_cols[i])\n    ax.grid(True, linestyle=\"--\")\n\n    ax_r = ax.twinx()\n    ax_r.plot(responders[6], color='r', alpha=0.5)\n    ax_r.set_ylabel('Responder_6', color='r')\n    \nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-17T00:11:25.213009Z","iopub.execute_input":"2024-11-17T00:11:25.213402Z","iopub.status.idle":"2024-11-17T00:11:59.096229Z","shell.execute_reply.started":"2024-11-17T00:11:25.213368Z","shell.execute_reply":"2024-11-17T00:11:59.094451Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Check the predictive power of features using Information Coefficient ","metadata":{}},{"cell_type":"markdown","source":"**Information Coefficient (IC)**\n\nThe information coefficient measures the correlation between a stock's returns and the predicted returns from an alpha factor (i.e. a feature).","metadata":{}},{"cell_type":"code","source":"# calculate cross-section IC at each timestamp for all dates\n# very time-consuming...\n\ndates_arr = pl_train['date_id'].unique().sort().to_numpy()[-20:]\ntimestamps = pl_train['time_id'].unique().sort().to_numpy()\n\n# define the pair of feature-target for IC calculation\nfeat_cols = ['feature_78', 'feature_62', 'feature_53', 'feature_36']\nresp_col = 'responder_6'\n\ncorr_mat = np.zeros((len(feat_cols), len(dates_arr), len(timestamps)))\npval_mat = np.zeros_like(corr_mat)\n\nfor i, _feat in enumerate(feat_cols):\n    for j, _date in tqdm(enumerate(dates_arr), total=len(dates_arr)):\n        for k, _time in enumerate(timestamps):\n            _filter = pl_train.filter((pl.col('date_id')==_date) & (pl.col('time_id')==_time))\n            responder_arr = _filter[resp_col].to_numpy()\n            feature_arr = _filter[_feat].to_numpy()\n\n            # pearsonr for simple ic / spearmanr for ranked ic \n            corr, p_value = spearmanr(feature_arr, responder_arr)\n\n            corr_mat[i, j, k] = corr\n            pval_mat[i, j, k] = p_value","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import joblib\n\njoblib.dump(\n    {\n        \"corr_matrix\": corr_mat,\n        \"pval_matrix\": pval_mat\n    },\n    \"correlation_matrix.pkl\"\n)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from scipy.stats import chi2\n\nfor i, _feat in enumerate(feat_cols):\n\n    fig, axes = plt.subplots(2, 1, figsize=(12, 5), sharex=True)\n\n    corr_mean = np.nanmean(corr_mat[i], axis=0)\n    axes[0].plot(timestamps, corr_mean)\n    axes[0].axhline(np.nanmean(corr_mean))\n    axes[0].set_ylabel('Correlation')\n\n    # combine p-val using Fisher’s method\n    X = -2 * np.nansum(np.log(pval_mat[i]), axis=0)\n    p_combined = 1 - chi2.cdf(X, 2*pval_mat.shape[1])\n\n    axes[1].plot(timestamps, p_combined)\n    axes[1].set_xlabel('Time')\n    axes[1].set_ylabel('p-value')\n\n    corr_mean_pval = corr_mean.copy()\n    corr_mean_pval[p_combined >= 0.05] = np.nan \n    axes[0].scatter(timestamps, corr_mean_pval, s=10, c='r', marker='o', label='p<0.05')\n\n    for ax in axes:\n        ax.set_xticks(timestamps[::50])\n        ax.grid(True, ls=\"--\")\n        ax.legend()\n\n    axes[0].set_title(f\"{_feat} - {resp_col}\")\n\n    fig.tight_layout()\n    plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}