{"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":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Understanding Weighted $R^2$","metadata":{}},{"cell_type":"markdown","source":"![](https://img.memegenerator.net/instances/71735329.jpg)\n\n\nIn this kernel we take a deep dive into weighted $R^2$ for the [2024 Jane Street competition](https://www.kaggle.com/competitions/jane-street-real-time-market-data-forecasting).\n\nDue to the interesting structure of the competition weights in the training data we will also learn more about the [Gamma Distribution](https://en.wikipedia.org/wiki/Gamma_distribution).\n\nThis is the 4th notebook in a series on Kaggle competition metrics. \n\nPrevious Notebooks:\n\n[1. RMSLE](https://www.kaggle.com/code/carlolepelaars/understanding-the-metric-rmsle)\n\n[2. Quadratic Weighted Kappa](https://www.kaggle.com/code/carlolepelaars/understanding-the-metric-quadratic-weighted-kappa)\n\n[3. Spearman's Rho](https://www.kaggle.com/code/carlolepelaars/understanding-the-metric-spearman-s-rho)","metadata":{}},{"cell_type":"markdown","source":"# Preparation","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport polars as pl\nimport altair as alt\nfrom tqdm import tqdm\nfrom scipy.stats import gamma\nfrom collections import defaultdict\nfrom sklearn.metrics import r2_score\n\nimport kaggle_evaluation.jane_street_inference_server","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-10-23T11:07:49.300286Z","iopub.execute_input":"2024-10-23T11:07:49.301302Z","iopub.status.idle":"2024-10-23T11:07:51.813578Z","shell.execute_reply.started":"2024-10-23T11:07:49.301220Z","shell.execute_reply":"2024-10-23T11:07:51.812136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The full training data is too large to load into a Kaggle Notebook, so we load 1 partition. This provides us with sufficient data to explore weighted $R^2$.","metadata":{}},{"cell_type":"code","source":"BASE_PATH = \"/kaggle/input/jane-street-real-time-market-data-forecasting/\"\ntrain = pd.read_parquet(BASE_PATH + 'train.parquet/partition_id=0/part-0.parquet')","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:07:51.815924Z","iopub.execute_input":"2024-10-23T11:07:51.816579Z","iopub.status.idle":"2024-10-23T11:07:57.569719Z","shell.execute_reply.started":"2024-10-23T11:07:51.816530Z","shell.execute_reply":"2024-10-23T11:07:57.568267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def list_files_recursive(base_path):\n    print('\\n# Files and file sizes')\n    aggregated_sizes = defaultdict(float)\n    \n    for root, dirs, files in os.walk(base_path):\n        for file in files:\n            file_path = os.path.join(root, file)\n            size_mb = os.path.getsize(file_path) / 1000000\n            relative_path = os.path.relpath(file_path, base_path)\n            path_parts = relative_path.split(os.sep)\n            top_level = path_parts[0]\n            \n            if len(path_parts) > 1:\n                aggregated_sizes[top_level] += size_mb\n            else:\n                aggregated_sizes[relative_path] = size_mb\n    for path, size in sorted(aggregated_sizes.items()):\n        print(f'{path.ljust(50)}| {round(size, 2)} MB')\n\nlist_files_recursive(BASE_PATH)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-10-23T11:24:59.861200Z","iopub.execute_input":"2024-10-23T11:24:59.861828Z","iopub.status.idle":"2024-10-23T11:24:59.907340Z","shell.execute_reply.started":"2024-10-23T11:24:59.861771Z","shell.execute_reply":"2024-10-23T11:24:59.905927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The Metric","metadata":{}},{"cell_type":"markdown","source":"The weighted $R^2$ is [formally defined by Kaggle](https://www.kaggle.com/competitions/jane-street-real-time-market-data-forecasting) as:\n\n$$R^2 = 1 - \\frac{\\Sigma w_i (y_i - \\hat{y_i})^2}{\\Sigma w_i y^2_i}$$\n\nWe can see that in the nominator the squared error of each sample $(y_i - \\hat{y_i})^2$ is weighted by $w_i$. These weights are given with the training data and its distribution is explored below. In the denominator the error is normalized with the weight squared target value. \n\nIt is common to center around the mean in the denominator, but this is not done in the formula given by Kaggle. This is because the target distribution is already centered around $0$. \n\nIn scikit-learn we can calculate the uncentered $R^2$ using [explained_variance_score](https://scikit-learn.org/1.5/modules/generated/sklearn.metrics.explained_variance_score.html) and the centered version with [r2_score](https://scikit-learn.org/1.5/modules/generated/sklearn.metrics.r2_score.html). Due to the centering around $0$, for this competition the two metrics will yield very similar results. \n\nIf there is no error in our predictions then the score will be $1$. With random predictions the score is most likely negative. \n\nWhen we look at the formula closely we can see that the error can be potentially infinite so the minimum of $R^2$ is $-\\infty$. Thus the range of $R^2$ is $[-\\infty...1]$.","metadata":{}},{"cell_type":"markdown","source":"# Data Distributions","metadata":{}},{"cell_type":"markdown","source":"## Target","metadata":{}},{"cell_type":"markdown","source":"The main target `responder_6` is a continuous target in range $[-5...5]$. We can frame it as a regression problem and should clip predictions between $-5$ and $5$.","metadata":{}},{"cell_type":"code","source":"train['responder_6'].plot(kind=\"hist\", bins=100, title=\"Target distribution\", xlabel=\"responder_6\");","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:28:45.522260Z","iopub.execute_input":"2024-10-23T11:28:45.523994Z","iopub.status.idle":"2024-10-23T11:28:46.311431Z","shell.execute_reply.started":"2024-10-23T11:28:45.523932Z","shell.execute_reply":"2024-10-23T11:28:46.310021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Weight","metadata":{}},{"cell_type":"markdown","source":"Weights are given in the training data and are used for the evaluation.\n\nThe weight distribution seems to be in range $[0.4...6]$ and resembles a [Gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution). Let's fit this distribution to get good weight estimates for the full dataset.","metadata":{}},{"cell_type":"code","source":"train['weight'].plot(kind=\"hist\", bins=100, title=\"Training Weight Distribution\", \n                     xlabel=\"Weight\", ylabel=\"Probability Density\", density=True);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-10-23T11:07:58.457797Z","iopub.execute_input":"2024-10-23T11:07:58.458319Z","iopub.status.idle":"2024-10-23T11:07:59.308711Z","shell.execute_reply.started":"2024-10-23T11:07:58.458262Z","shell.execute_reply":"2024-10-23T11:07:59.307333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shape, loc, scale = gamma.fit(train['weight'])","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:07:59.310285Z","iopub.execute_input":"2024-10-23T11:07:59.310748Z","iopub.status.idle":"2024-10-23T11:08:07.979630Z","shell.execute_reply.started":"2024-10-23T11:07:59.310703Z","shell.execute_reply":"2024-10-23T11:08:07.978156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"shape (k) = {shape:.2f}\")\nprint(f\"loc (min) = {loc:.2f}\")\nprint(f\"scale (θ) = {scale:.2f}\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-10-23T11:08:07.981173Z","iopub.execute_input":"2024-10-23T11:08:07.981642Z","iopub.status.idle":"2024-10-23T11:08:07.989726Z","shell.execute_reply.started":"2024-10-23T11:08:07.981596Z","shell.execute_reply":"2024-10-23T11:08:07.988382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The `shape` ($k$) denotes the center of the distribution. It is similar to the $\\lambda$ parameter in a [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution). The `shape` in a Gamma distribution has a slightly different meaning from $\\lambda$, because there are also `scale` and `loc` parameters involved. \n\nThe `scale` ($\\theta$) determines how much of the probability density is centered around the mean. The larger the scale, the more the probability density will be spread out. The Jane Street weight distribution has a relatively small scale. This leads to most of weight being in range `[loc...3]`.\n\nThis brings us to the last parameter, `loc`. `loc` shifts the distribution and is most often $0$, which means the data minimum is $0$ and values are strictly positive. In our case its around $0.39$. Which is the minimum weight that the distribution will generate. The minimum weight in our training data is actually $0.44$. This shows that the Gamma distribution is not able to fit the weight distribution perfectly, but its close!","metadata":{}},{"cell_type":"code","source":"x = np.linspace(0, 6, 1000)\ny = gamma.pdf(x, a=shape, scale=scale, loc=loc)\ndf = pd.DataFrame({\"x\": x, \"probability\": y})\n\nalt.Chart(df).mark_area(\n    opacity=0.7,\n    color=\"#007bff\",\n    line={\"color\": \"#0056b3\"}\n).encode(\n    x=alt.X(\"x:Q\", title=\"Value\", axis=alt.Axis(labelFontSize=12, titleFontSize=14), scale=alt.Scale(domain=[0, 6])),\n    y=alt.Y(\"probability:Q\", title=\"Probability Density\", axis=alt.Axis(labelFontSize=12, titleFontSize=14)),\n    tooltip=[alt.Tooltip(\"x:Q\", title=\"Value\", format=\".2f\"), alt.Tooltip(\"probability:Q\", title=\"Probability Density\", format=\".4f\")]\n).properties(\n    width=700,\n    height=300,\n    title=alt.TitleParams(text=f\"Gamma Distribution (k = {shape:.2f}, θ = {scale:.2f})\", fontSize=20)\n).configure_view(\n    strokeWidth=0\n).configure_axis(\n    grid=False\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-10-23T11:08:07.991283Z","iopub.execute_input":"2024-10-23T11:08:07.991984Z","iopub.status.idle":"2024-10-23T11:08:08.112556Z","shell.execute_reply.started":"2024-10-23T11:08:07.991936Z","shell.execute_reply":"2024-10-23T11:08:08.111133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Baselines","metadata":{}},{"cell_type":"markdown","source":"## Random Scores","metadata":{}},{"cell_type":"markdown","source":"Most of the time $R^2$ scores will be in range $[0...1]$, but can be negative. Random scores will generally score below $0$, which is worse than predicting a constant. The score for random predictions in range $[-5...5]$ using real targets and weights is around $-11.5$. The score can potentially be $-\\infty$ if we don't clip our predictions.","metadata":{}},{"cell_type":"code","source":"size = 1000\nn = 10000\nsample = train.sample(size)\ny_true = sample[\"responder_6\"]\nweights = sample[\"weight\"]\n\nr2s = []\nfor i in tqdm(range(n)):\n    y_pred = np.random.uniform(-5, 5, size=size)\n    score = r2_score(y_true, y_pred, sample_weight=weights)\n    r2s.append(score)","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:08:08.114369Z","iopub.execute_input":"2024-10-23T11:08:08.114854Z","iopub.status.idle":"2024-10-23T11:08:15.235677Z","shell.execute_reply.started":"2024-10-23T11:08:08.114810Z","shell.execute_reply":"2024-10-23T11:08:15.234504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.Series(r2s).plot(kind=\"hist\", bins=100, title=\"$R^2$ For Random Predictions\", xlabel=\"Score\");","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:08:15.240037Z","iopub.execute_input":"2024-10-23T11:08:15.240621Z","iopub.status.idle":"2024-10-23T11:08:15.862227Z","shell.execute_reply.started":"2024-10-23T11:08:15.240563Z","shell.execute_reply":"2024-10-23T11:08:15.860677Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"R^2 scores from random targets and predictions in range [-5...5]:\")\nfor _ in range(5):\n    print(r2_score(np.random.uniform(-5, 5, size=5), np.random.uniform(-5, 5, size=5), sample_weight=[0.1, 0.2, 0.3, 0.4, 0.5]))","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:08:15.864166Z","iopub.execute_input":"2024-10-23T11:08:15.864752Z","iopub.status.idle":"2024-10-23T11:08:15.878107Z","shell.execute_reply.started":"2024-10-23T11:08:15.864682Z","shell.execute_reply":"2024-10-23T11:08:15.876576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Constant Predictions","metadata":{}},{"cell_type":"markdown","source":"Because the Jane Street data is centered around $0$ the best constant prediction we can make is $0$. This will result in an $R^2$ score of $0$. Any constant prediction that deviates from this will result in a negative score.","metadata":{"_kg_hide-input":true}},{"cell_type":"code","source":"size = 1000\nn = 1000\nsample = train.sample(size)\ny_true = sample[\"responder_6\"]\nweights = sample[\"weight\"]\n\nr2s = []\nfor i in tqdm(np.linspace(-5, 5, n)):\n    y_pred = np.array([i] * size)\n    score = r2_score(y_true, y_pred, sample_weight=weights)\n    r2s.append(score)","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:14:14.455360Z","iopub.execute_input":"2024-10-23T11:14:14.456569Z","iopub.status.idle":"2024-10-23T11:14:15.311601Z","shell.execute_reply.started":"2024-10-23T11:14:14.456509Z","shell.execute_reply":"2024-10-23T11:14:15.310218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = pd.DataFrame({\n    'x': np.linspace(-5, 5, len(r2s)),\n    'R-squared': r2s\n})\n\n# Create the Altair chart\nalt.Chart(data).mark_line().encode(\n    x=alt.X('x', title='Constant', scale=alt.Scale(domain=[-5, 5])),\n    y=alt.Y('R-squared', title='R-squared')\n).properties(\n    width=600,\n    height=400,\n    title='R-squared scores for constants'\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:21:24.065174Z","iopub.execute_input":"2024-10-23T11:21:24.066601Z","iopub.status.idle":"2024-10-23T11:21:24.125524Z","shell.execute_reply.started":"2024-10-23T11:21:24.066542Z","shell.execute_reply":"2024-10-23T11:21:24.123648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Perfect Score","metadata":{}},{"cell_type":"markdown","source":"![](https://lh3.googleusercontent.com/proxy/7FRwnnqtUtiB9jIbkw5147Cenukdi1HtYCNCc3h_fjKe_18NB2kyb1CdaPkxntBsokhphk9b6XE-XgzTdXsi968y56h5sOX39w8TNnUnhA4MtIVMWsEFZFTXg-tXf0PNAXhWo3_6j3E0pC_9mrS5gij8lADH)\n\n\n\nA perfect $R^2$ score is $1$. In this case the model gives perfect predicting and explains `100%` of the variance in the data.","metadata":{}},{"cell_type":"code","source":"r2_score([1,0,1,0,1], [1,0,1,0,1], sample_weight=[0.1, 0.2, 0.3, 0.4, 0.5])","metadata":{"execution":{"iopub.status.busy":"2024-10-23T11:15:11.013762Z","iopub.execute_input":"2024-10-23T11:15:11.014469Z","iopub.status.idle":"2024-10-23T11:15:11.027467Z","shell.execute_reply.started":"2024-10-23T11:15:11.014382Z","shell.execute_reply":"2024-10-23T11:15:11.026005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I expect the top $R^2$ scores for the Jane Street competition to be somewhere between $0.01$ and $0.015$. This is low for ML models in general, but really high for financial models due to low [signal-to-noise](https://en.wikipedia.org/wiki/Signal-to-noise_ratio) in the data.","metadata":{}},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"markdown","source":"For a test submission we use sample code given in the [Demo Submission Notebook](https://www.kaggle.com/code/ryanholbrook/jane-street-rmf-demo-submission). It creates a `predict` function that is submitted to the Kaggle evaluation server. \n\nAs we have seen in this notebook, random predictions perform worse compared to picking the best constant. Our baseline will therefore be a prediction of all $0$. This will result in an $R^2$ score of $0$.","metadata":{}},{"cell_type":"code","source":"lags_ : pl.DataFrame | None = None\n\n# Replace this function with your inference code.\n# You can return either a Pandas or Polars dataframe, though Polars is recommended.\n# Each batch of predictions (except the very first) must be returned within 10 minutes of the batch features being provided.\ndef predict(test: pl.DataFrame, lags: pl.DataFrame | None) -> pl.DataFrame | pd.DataFrame:\n    \"\"\"Make a prediction.\"\"\"\n    # All the responders from the previous day are passed in at time_id == 0. We save them in a global variable for access at every time_id.\n    # Use them as extra features, if you like.\n    global lags_\n    if lags is not None:\n        lags_ = lags\n\n    # Predictions are clipped between -5 and 5.\n    predictions = test.select(\n        'row_id',\n        pl.lit(0.0).clip(-5, 5).alias('responder_6'),\n    )\n\n    # The predict function must return a DataFrame\n    assert isinstance(predictions, pl.DataFrame | pd.DataFrame)\n    # with columns 'row_id', 'responder_6'\n    assert predictions.columns == ['row_id', 'responder_6']\n    # and as many rows as the test data.\n    assert len(predictions) == len(test)\n\n    return predictions","metadata":{"execution":{"iopub.status.busy":"2024-10-16T15:23:34.355150Z","iopub.execute_input":"2024-10-16T15:23:34.355624Z","iopub.status.idle":"2024-10-16T15:23:34.363282Z","shell.execute_reply.started":"2024-10-16T15:23:34.355578Z","shell.execute_reply":"2024-10-16T15:23:34.362035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inference_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":{"execution":{"iopub.status.busy":"2024-10-16T15:23:44.515979Z","iopub.execute_input":"2024-10-16T15:23:44.516444Z","iopub.status.idle":"2024-10-16T15:23:44.826826Z","shell.execute_reply.started":"2024-10-16T15:23:44.516399Z","shell.execute_reply":"2024-10-16T15:23:44.825664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"![](https://media.makeameme.org/created/when-you-finish-b158318ec5.jpg)","metadata":{}},{"cell_type":"markdown","source":"**That's it! If you like this Kaggle kernel, feel free to give an upvote and leave a comment!**","metadata":{}}]}