{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "55201f8b",
   "metadata": {},
   "source": [
    "### Chapter 4 \u2014 Trading and Risk Models\n",
    "\n",
    "Run from a fresh session, top to bottom. This notebook mirrors the public self-study slides. Repair exercises contain worked solutions; predict and attempt them before reading the solution. All data are public classroom artifacts, with survivor membership and assumed clocks. No broker or private research system is accessed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21ce1e18",
   "metadata": {},
   "source": [
    "### What This Chapter Builds\n",
    "\n",
    "- A continuous next-month return target and a fixed adverse-event label from the Chapter 2 handoff.\n",
    "- Multiple regression: coefficients, residuals, assumptions, inference and forecast errors.\n",
    "- Defensible model selection: interactions, best subset, VIF and regularisation.\n",
    "- Trees, forests, boosting and logistic probabilities with separate ranking, calibration and decision checks.\n",
    "- Causal designs and attention features, with explicit assumptions and limitations.\n",
    "- A walk-forward return/risk pipeline and a cash-and-fee account."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ded95c44",
   "metadata": {},
   "source": [
    "### How to study this chapter independently\n",
    "\n",
    "Read each section\u2019s concept and formula before running Python. Predict the result, run cells from top to bottom, modify one input, then explain what changed and why. A new section deck starts a fresh Python session.\n",
    "\n",
    "The slides are the complete self-study route. The Book expands intuition, derivations and failure modes. Classroom time selects discussions and experiments from this material, rather than limiting its coverage.\n",
    "\n",
    "Bring one coefficient interpretation, one model-comparison question and one data/risk limitation. The case demonstrates statistical methods and can transfer to forecasting and adverse-event decisions outside finance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a7dc1c9",
   "metadata": {},
   "source": [
    "### Roadmap\n",
    "\n",
    "| Section | Statistical question | Evidence |\n",
    "|---|---|---|\n",
    "| \u00a74.1 | What does a multiple regression estimate? | Coefficients, residuals and forecast errors |\n",
    "| \u00a74.2 | Which specification and uncertainty claims are defensible? | Inference, validation, regularisation and VIF |\n",
    "| \u00a74.3 | How do nonlinear models and event probabilities work? | Trees, ensembles, ROC, Brier and thresholds |\n",
    "| \u00a74.4 | What would change under an intervention? | Confounding, DiD and IV assumptions |\n",
    "| \u00a74.5 | How can cross-sectional context form features? | Fixed attention and walk-forward comparison |\n",
    "| \u00a74.6 | How do forecasts enter a risk decision? | Frozen rules, exposure, cash, costs and NAV |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3c04340",
   "metadata": {},
   "source": [
    "### Three-hour classroom application\n",
    "\n",
    "| Minutes | Activity |\n",
    "|---|---|\n",
    "| 0\u201315 | Readiness discussion: target, coefficient and timing |\n",
    "| 15\u201335 | Guided regression / classification interpretation |\n",
    "| 35\u201365 | Group experiment on a fixed validation sample |\n",
    "| 65\u201375 | Break |\n",
    "| 75\u2013135 | Walk-forward risk-rule and accounting investigation |\n",
    "| 135\u2013165 | Peer challenge: statistics, causality, costs and limitations |\n",
    "| 165\u2013180 | Revise the decision memo and individual explanation |\n",
    "\n",
    "Prepare all assigned self-study sections. Each group submits one runnable notebook and a short memo through the existing course route. This schedule adds no new grading rule."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12108593",
   "metadata": {},
   "source": [
    "### Return and Loss Targets: Data Contracts, Multiple Regression and OLS\n",
    "\n",
    "Study the definitions, then predict, run, modify and explain the worked results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08644fa5",
   "metadata": {},
   "source": [
    "### One table, two statistical questions\n",
    "\n",
    "**Regression:** What is the next-month simple return $Y_{i,t+1}$ conditional on today\u2019s features?\n",
    "\n",
    "**Classification:** What is the probability of $Y_{i,t+1}\\le-0.05$? This event means a next-month closing-to-closing loss of at least 5%. It does not mean an intramonth 5% drawdown.\n",
    "\n",
    "The public Chapter 2 handoff supplies nine stocks, 71 signal dates and 639 candidates. This small current-survivor sample teaches statistical reasoning, rather than validating a profitable strategy.\n",
    "\n",
    "The same questions arise in demand forecasting, equipment failures and insurance: predict a continuous outcome or the probability of a defined adverse event."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1a1a1d0",
   "metadata": {},
   "source": [
    "### The modelling path preserves the data contract\n",
    "\n",
    "**Diagram concepts:** Chapter 2: keyed available features; Chronological train / validation / test; OLS: next-month return; Logit: probability of loss \u2265 5%; Evaluate prediction and uncertainty; Freeze a decision rule; Walk forward and reconcile net accounts.\n",
    "\n",
    "Keep keys, feature availability and label maturity visible through every stage. A missing future label prevents fitting or evaluation, but does not remove the current prediction candidate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ea957cc",
   "metadata": {},
   "source": [
    "### Load the case and freeze the three clocks\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "6a925616",
   "metadata": {},
   "source": [
    "import numpy as np, pandas as pd, matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import Ridge, LogisticRegression\n",
    "from sklearn.metrics import mean_squared_error, roc_auc_score, brier_score_loss\n",
    "case = pd.read_csv(\"https://statpython.pages.dev/data/chapter2_scored.csv\",\n",
    "                   parse_dates=[\"date\", \"label_available_at\"])\n",
    "case = case.sort_values([\"date\", \"ticker\"]).reset_index(drop=True)\n",
    "model_keys = [\"date\", \"ticker\"]\n",
    "assert not case.duplicated(model_keys).any()\n",
    "numeric_features = [\"ret_12_1\", \"rvol_21d\", \"momentum_sector_z\"]\n",
    "sector_dummies = pd.get_dummies(case.sector, prefix=\"sector\", drop_first=True, dtype=float)\n",
    "X_case = pd.concat([case[numeric_features], sector_dummies], axis=1)\n",
    "assert np.isfinite(X_case.to_numpy()).all()\n",
    "y_case = case.ret_next.where(case.label_valid)\n",
    "loss_event = y_case.le(-0.05).where(case.label_valid)\n",
    "train_case = case.date.lt(\"2021-01-01\") & case.label_available_at.lt(\"2021-01-01\") & case.label_valid\n",
    "valid_case = (case.date.ge(\"2021-01-01\") & case.date.lt(\"2022-01-01\")\n",
    "              & case.label_available_at.lt(\"2022-01-01\") & case.label_valid)\n",
    "test_case = case.date.ge(\"2022-01-01\") & case.label_valid\n",
    "print(\"train / validation / test / pending:\",\n",
    "      int(train_case.sum()), int(valid_case.sum()), int(test_case.sum()), int((~case.label_valid).sum()))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "9893dc65",
   "metadata": {},
   "source": [
    "The split prints **207 training, 99 validation, 306 test and 9 pending rows**. Freeze fitting at 1 January 2021 and rule selection at 1 January 2022. December 2020 and December 2021 each contribute nine boundary rows whose labels arrive after the respective cutoff.\n",
    "\n",
    "The table\u2019s stored predictions are ignored. We fit new models using only the three named features and two sector dummies. Future returns, label status and prior predictions never enter X."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af185953",
   "metadata": {},
   "source": [
    "### Read the feature matrix before fitting\n",
    "\n",
    "| Input | Meaning | Unit |\n",
    "|---|---|---|\n",
    "| `ret_12_1` | Twelve-to-one momentum | Simple cumulative return |\n",
    "| `rvol_21d` | Trailing realised volatility | Source volatility units; verify annualisation metadata |\n",
    "| `momentum_sector_z` | Momentum relative to same-date sector peers | Standard deviations |\n",
    "| Sector dummies | Financial Services and Technology versus Energy | 0/1 indicators |\n",
    "\n",
    "An intercept plus every sector dummy is collinear. `drop_first=True` supplies a reference group. Static sectors and assumed month-end availability remain classroom limitations.\n",
    "\n",
    "**Predict:** Why would fitting on every numeric column accidentally admit future information?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df941cef",
   "metadata": {},
   "source": [
    "### Multiple regression separates conditional associations\n",
    "\n",
    "$$Y_i=\\beta_0+\\beta_1 x_{i1}+\\cdots+\\beta_p x_{ip}+\\varepsilon_i,\\qquad E[\\varepsilon_i\\mid X_i]=0.$$\n",
    "\n",
    "$\\beta_j$ is the conditional change in expected Y for a one-unit increase in $x_j$, holding the other regressors fixed within the model. Correlated features can have very different marginal and conditional associations.\n",
    "\n",
    "The residual $e_i=y_i-\\widehat y_i$ is observed after fitting. The unobserved error $\\varepsilon_i$ belongs to the population model.\n",
    "\n",
    "A predictive coefficient requires stability to generalise. A causal interpretation requires additional identification assumptions (\u00a74.4)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c9b40ef",
   "metadata": {},
   "source": [
    "### OLS chooses coefficients by squared residuals\n",
    "\n",
    "$$\\widehat\\beta=\\arg\\min_b\\sum_i(y_i-x_i^\\top b)^2.$$\n",
    "\n",
    "With full column rank, $\\widehat\\beta=(X^\\top X)^{-1}X^\\top y$. Numerical software uses more stable linear algebra rather than requiring you to form the inverse.\n",
    "\n",
    "Squared loss heavily penalises large errors. An extreme return can materially influence coefficients. Check leverage, residuals and sensitivity before deleting an observation.\n",
    "\n",
    "**Check:** A low training sum of squares establishes fit to observed outcomes; it does not establish future prediction quality."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dceee2f0",
   "metadata": {},
   "source": [
    "### Fit OLS and compare training with validation\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "502c0e5f",
   "metadata": {},
   "source": [
    "design = sm.add_constant(X_case, has_constant=\"add\")\n",
    "ols_case = sm.OLS(y_case.loc[train_case], design.loc[train_case]).fit()\n",
    "robust_case = sm.OLS(y_case.loc[train_case], design.loc[train_case]).fit(\n",
    "    cov_type=\"cluster\", cov_kwds={\"groups\": case.loc[train_case, \"date\"]})\n",
    "validation_prediction = ols_case.predict(design.loc[valid_case])\n",
    "print(\"training R2 / adjusted R2:\", round(ols_case.rsquared, 4), round(ols_case.rsquared_adj, 4))\n",
    "print(\"validation RMSE:\", round(np.sqrt(mean_squared_error(y_case.loc[valid_case], validation_prediction)), 4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "49083976",
   "metadata": {},
   "source": [
    "Training R\u00b2 is **0.1766**, adjusted R\u00b2 **0.1561**, while validation RMSE is **0.1014** in return units (10.14 percentage points).\n",
    "\n",
    "The clustered fit has identical point estimates but different standard errors. It allows within-date correlation across stocks. Only 23 training dates supply clusters, so large-sample inference remains approximate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7deb2bde",
   "metadata": {},
   "source": [
    "### Interpret a coefficient in the right units\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "e543bcfa",
   "metadata": {},
   "source": [
    "print(ols_case.params.round(4))\n",
    "print(\"predicted return change for +10 percentage points momentum:\",\n",
    "      round(float(0.10 * ols_case.params[\"ret_12_1\"]), 4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "43346a45",
   "metadata": {},
   "source": [
    "The momentum coefficient is about **\u22120.0484**. A 0.10 increase in momentum corresponds to about **\u22120.00484** predicted return, conditional on the other inputs in this fitted sample.\n",
    "\n",
    "The Technology coefficient is about **+0.0880** relative to Energy, conditional on the numeric features. This is a sample association, rather than a prescription to change a stock\u2019s sector."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b343a7e9",
   "metadata": {},
   "source": [
    "### Assumptions have different jobs\n",
    "\n",
    "| Assumption | What it supports |\n",
    "|---|---|\n",
    "| Full rank | Identifiable coefficients |\n",
    "| Correct conditional mean / exogeneity | Unbiased or consistent OLS under suitable sampling conditions |\n",
    "| Constant error variance and uncorrelated errors | Classical OLS variance formula |\n",
    "| Normal errors with fixed-X classical assumptions | Exact finite-sample t/F inference |\n",
    "| Suitable temporal stability | Generalisation to future periods |\n",
    "\n",
    "OLS can produce predictions without Normal errors. Robust covariance changes inference, rather than repairing a wrong target, leakage or an unstable conditional mean."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5db83ea9",
   "metadata": {},
   "source": [
    "### R\u00b2, adjusted R\u00b2 and forecast errors\n",
    "\n",
    "$$R^2=1-\\frac{\\sum_i e_i^2}{\\sum_i(y_i-\\bar y)^2},\\qquad\n",
    "\\bar R^2=1-(1-R^2)\\frac{n-1}{n-p-1}.$$\n",
    "\n",
    "These formulas assume a model with an intercept. Adding regressors cannot reduce training R\u00b2. Adjusted R\u00b2 penalises their number, but remains a training criterion.\n",
    "\n",
    "$$\\operatorname{RMSE}=\\sqrt{n^{-1}\\sum_i(y_i-\\widehat y_i)^2},\\qquad\n",
    "\\operatorname{MAE}=n^{-1}\\sum_i|y_i-\\widehat y_i|.$$\n",
    "\n",
    "Evaluate future predictions against a simple training-fitted baseline. Out-of-sample R\u00b2 can be negative."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4967a7e4",
   "metadata": {},
   "source": [
    "### A baseline makes model error interpretable\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "f1fccb36",
   "metadata": {},
   "source": [
    "baseline_return = pd.Series(y_case.loc[train_case].mean(), index=case.index[valid_case])\n",
    "comparison = pd.Series({\n",
    "    \"training-mean baseline RMSE\": np.sqrt(mean_squared_error(y_case.loc[valid_case], baseline_return)),\n",
    "    \"OLS validation RMSE\": np.sqrt(mean_squared_error(y_case.loc[valid_case], validation_prediction))})\n",
    "print(comparison.round(4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "aaa53ce2",
   "metadata": {},
   "source": [
    "Compare both errors on exactly the same validation rows. A flexible model is useful only if the improvement serves the intended decision and survives future evaluation.\n",
    "\n",
    "**Modify:** Restrict the numeric feature list, rerun fitting, and explain which quantities may improve mechanically and which need validation evidence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b8ce542",
   "metadata": {},
   "source": [
    "### Self-check: inference, prediction and timing\n",
    "\n",
    "**Predict:** Which statement is justified by the OLS example?\n",
    "\n",
    "- Training R\u00b2 proves the signal will make money\n",
    "- A coefficient measures a conditional fitted association, and future forecast quality needs chronological evaluation\n",
    "- Normal-looking features guarantee independent errors\n",
    "- Pending outcomes should be filled with zero before training\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "A coefficient measures a conditional fitted association, and future forecast quality needs chronological evaluation. OLS coefficients describe conditional fitted associations. Prediction, causal interpretation and portfolio profitability require separate evidence.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92dce557",
   "metadata": {},
   "source": [
    "### Repair the regression sample\n",
    "\n",
    "Select only the declared eligible training rows. The starter currently includes future validation/test dates and pending labels."
   ]
  },
  {
   "cell_type": "code",
   "id": "cba08750",
   "metadata": {},
   "source": [
    "regression_rows = case.loc[train_case].copy()  # replace with case.loc[...]\n",
    "print(len(regression_rows))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "88ac9a1c",
   "metadata": {},
   "source": [
    "### OLS does not require Normal predictors\n",
    "\n",
    "**Predict:** Which condition is needed to define an identifiable OLS coefficient vector?\n",
    "\n",
    "- Every predictor is normally distributed\n",
    "- The design matrix has full column rank\n",
    "- The training R\u00b2 exceeds 0.5\n",
    "- Every outcome is positive\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "The design matrix has full column rank. Full column rank prevents redundant coefficient combinations. Normal predictors are not an OLS requirement.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16bb68fc",
   "metadata": {},
   "source": [
    "### Training and validation measure different things\n",
    "\n",
    "**Predict:** Adding an irrelevant regressor to OLS with an intercept can do what?\n",
    "\n",
    "- Reduce training R\u00b2 mechanically\n",
    "- Leave or increase training R\u00b2 while worsening validation error\n",
    "- Guarantee a lower validation RMSE\n",
    "- Eliminate residual variance\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Leave or increase training R\u00b2 while worsening validation error. The larger model includes the smaller model as a special case, so its training SSE cannot rise. Future error can worsen.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a268d089",
   "metadata": {},
   "source": [
    "### What you discovered\n",
    "\n",
    "- Define the return and loss-event targets with a horizon.\n",
    "- Preserve candidate keys and available labels through chronological splits.\n",
    "- Interpret OLS coefficients, residuals and assumptions.\n",
    "- Compare forecast errors against a training-fitted baseline."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6af3e2fc",
   "metadata": {},
   "source": [
    "### Model Selection and Inference: Interactions, Subsets, Regularisation and VIF\n",
    "\n",
    "Study the definitions, then predict, run, modify and explain the worked results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddb24701",
   "metadata": {},
   "source": [
    "### Rebuild the case before model selection\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "c5cda2e1",
   "metadata": {},
   "source": [
    "import numpy as np, pandas as pd, matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import Ridge, LogisticRegression\n",
    "from sklearn.metrics import mean_squared_error, roc_auc_score, brier_score_loss\n",
    "case = pd.read_csv(\"https://statpython.pages.dev/data/chapter2_scored.csv\",\n",
    "                   parse_dates=[\"date\", \"label_available_at\"])\n",
    "case = case.sort_values([\"date\", \"ticker\"]).reset_index(drop=True)\n",
    "model_keys = [\"date\", \"ticker\"]\n",
    "assert not case.duplicated(model_keys).any()\n",
    "numeric_features = [\"ret_12_1\", \"rvol_21d\", \"momentum_sector_z\"]\n",
    "sector_dummies = pd.get_dummies(case.sector, prefix=\"sector\", drop_first=True, dtype=float)\n",
    "X_case = pd.concat([case[numeric_features], sector_dummies], axis=1)\n",
    "assert np.isfinite(X_case.to_numpy()).all()\n",
    "y_case = case.ret_next.where(case.label_valid)\n",
    "loss_event = y_case.le(-0.05).where(case.label_valid)\n",
    "train_case = case.date.lt(\"2021-01-01\") & case.label_available_at.lt(\"2021-01-01\") & case.label_valid\n",
    "valid_case = (case.date.ge(\"2021-01-01\") & case.date.lt(\"2022-01-01\")\n",
    "              & case.label_available_at.lt(\"2022-01-01\") & case.label_valid)\n",
    "test_case = case.date.ge(\"2022-01-01\") & case.label_valid\n",
    "print(\"train / validation / test / pending:\",\n",
    "      int(train_case.sum()), int(valid_case.sum()), int(test_case.sum()), int((~case.label_valid).sum()))\n",
    "\n",
    "\n",
    "design = sm.add_constant(X_case, has_constant=\"add\")\n",
    "ols_case = sm.OLS(y_case.loc[train_case], design.loc[train_case]).fit()\n",
    "robust_case = sm.OLS(y_case.loc[train_case], design.loc[train_case]).fit(\n",
    "    cov_type=\"cluster\", cov_kwds={\"groups\": case.loc[train_case, \"date\"]})\n",
    "validation_prediction = ols_case.predict(design.loc[valid_case])\n",
    "print(\"training R2 / adjusted R2:\", round(ols_case.rsquared, 4), round(ols_case.rsquared_adj, 4))\n",
    "print(\"validation RMSE:\", round(np.sqrt(mean_squared_error(y_case.loc[valid_case], validation_prediction)), 4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "6c45e238",
   "metadata": {},
   "source": [
    "This section runs independently. Preserve the same training and validation clocks. Model selection uses validation data; the 2022\u20132024 test remains unopened until the rule is fixed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dfb96f9",
   "metadata": {},
   "source": [
    "### Confidence intervals describe coefficient uncertainty\n",
    "\n",
    "For a classical homoskedastic linear model,\n",
    "\n",
    "$$\\widehat{\\operatorname{Var}}(\\widehat\\beta)=s^2(X^\\top X)^{-1},\\qquad s^2=\\frac{\\sum_i e_i^2}{n-p-1}.$$\n",
    "\n",
    "A coefficient interval uses estimate \u00b1 critical value \u00d7 standard error. A t statistic tests a specified coefficient value, often $H_0:\\beta_j=0$.\n",
    "\n",
    "A confidence interval for the conditional mean differs from a prediction interval for one future return, which also includes outcome noise. Clustered standard errors address a different error structure and typically use approximate inference."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a337e304",
   "metadata": {},
   "source": [
    "### Read robust coefficient output\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "ac4152ca",
   "metadata": {},
   "source": [
    "interval = robust_case.conf_int()\n",
    "coefficient_report = pd.DataFrame({\"estimate\": robust_case.params,\n",
    "    \"cluster_se\": robust_case.bse, \"statistic\": robust_case.tvalues,\n",
    "    \"p_value\": robust_case.pvalues, \"lower\": interval[0], \"upper\": interval[1]})\n",
    "print(coefficient_report.round(4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "9661c179",
   "metadata": {},
   "source": [
    "State the covariance estimator, clustering unit and number of dates. Robust standard errors do not make the coefficient causal.\n",
    "\n",
    "Repeatedly testing coefficients after searching many subsets changes the inference problem. A nominal p-value from the selected model ignores selection unless the procedure accounts for it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27720d40",
   "metadata": {},
   "source": [
    "### An F test asks about several coefficients together\n",
    "\n",
    "For q restrictions in nested classical OLS models,\n",
    "\n",
    "$$F=\\frac{(SSE_R-SSE_U)/q}{SSE_U/(n-p_U-1)}.$$\n",
    "\n",
    "The null says all specified restrictions hold. A significant joint test need not imply every individual coefficient is significant, particularly with correlated regressors."
   ]
  },
  {
   "cell_type": "code",
   "id": "5e720473",
   "metadata": {},
   "source": [
    "joint_test = robust_case.wald_test(\"ret_12_1 = 0, momentum_sector_z = 0\", scalar=True)\n",
    "print(\"cluster-robust Wald statistic / p-value:\", round(float(joint_test.statistic), 3), round(float(joint_test.pvalue), 4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "3d436127",
   "metadata": {},
   "source": [
    "The code uses a cluster-robust Wald test with asymptotic calibration, rather than claiming the exact classical F distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "155b890b",
   "metadata": {},
   "source": [
    "### Residual plots test the story behind the fit\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "aff119b2",
   "metadata": {},
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.scatter(ols_case.fittedvalues, ols_case.resid, s=15, alpha=0.6)\n",
    "ax.axhline(0, color=\"black\", linestyle=\"--\")\n",
    "ax.set(xlabel=\"Fitted next-month return\", ylabel=\"Observed minus fitted return\")\n",
    "plt.show()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "17650626",
   "metadata": {},
   "source": [
    "A curve suggests a missing nonlinear mean relationship. A widening fan suggests changing error variance. Clusters by date may reflect common shocks. None can be ruled out by a single global residual mean near zero, which OLS with an intercept enforces mechanically.\n",
    "\n",
    "**Modify:** Colour the plot by year or sector. Explain which assumption the visible pattern questions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9002f39a",
   "metadata": {},
   "source": [
    "### Interactions and nonlinear terms change the question\n",
    "\n",
    "$$E[Y\\mid X]=\\beta_0+\\beta_1m+\\beta_2v+\\beta_3mv+\\beta_4m^2.$$\n",
    "\n",
    "The momentum slope becomes $\\beta_1+\\beta_3v+2\\beta_4m$. With an interaction, its main coefficient is the slope at the other interacting variable\u2019s zero value."
   ]
  },
  {
   "cell_type": "code",
   "id": "f9e5c8ae",
   "metadata": {},
   "source": [
    "X_expanded = X_case.assign(momentum_x_vol=X_case.ret_12_1 * X_case.rvol_21d,\n",
    "                           momentum_squared=X_case.ret_12_1 ** 2)\n",
    "expanded_design = sm.add_constant(X_expanded, has_constant=\"add\")\n",
    "expanded_fit = sm.OLS(y_case.loc[train_case], expanded_design.loc[train_case]).fit()\n",
    "expanded_valid = expanded_fit.predict(expanded_design.loc[valid_case])\n",
    "print(\"expanded train adjusted R2:\", round(expanded_fit.rsquared_adj, 4))\n",
    "print(\"expanded validation RMSE:\", round(np.sqrt(mean_squared_error(y_case.loc[valid_case], expanded_valid)), 4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "c1625a45",
   "metadata": {},
   "source": [
    "Feature construction uses information available at each row. Select transformations using training/validation, then freeze them before test evaluation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "459dd6d7",
   "metadata": {},
   "source": [
    "### Best subset is a search with a statistical cost\n",
    "\n",
    "Five candidate regressors produce $2^5-1=31$ non-empty subsets. Searching them can overfit the validation set even though each OLS fit is simple.\n",
    "\n",
    "The helper below iterates over **model specifications**, rather than processing data rows. It records every attempted subset. Adjusted R\u00b2 is a useful diagnostic, but we select by the declared validation RMSE."
   ]
  },
  {
   "cell_type": "code",
   "id": "a7cb40de",
   "metadata": {},
   "source": [
    "from itertools import combinations\n",
    "subset_records = []\n",
    "for size in range(1, X_case.shape[1] + 1):\n",
    "    for columns in combinations(X_case.columns, size):\n",
    "        matrix = sm.add_constant(X_case[list(columns)], has_constant=\"add\")\n",
    "        fitted = sm.OLS(y_case.loc[train_case], matrix.loc[train_case]).fit()\n",
    "        error = np.sqrt(mean_squared_error(y_case.loc[valid_case], fitted.predict(matrix.loc[valid_case])))\n",
    "        subset_records.append({\"columns\": columns, \"size\": size,\n",
    "                               \"train_adj_r2\": fitted.rsquared_adj, \"validation_rmse\": error})\n",
    "subset_report = pd.DataFrame(subset_records).sort_values(\"validation_rmse\")\n",
    "print(\"models tried:\", len(subset_report))\n",
    "print(subset_report.head(3).to_string(index=False))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "6203a719",
   "metadata": {},
   "source": [
    "### Multicollinearity inflates coefficient uncertainty\n",
    "\n",
    "Regress feature j on the remaining regressors to obtain $R_j^2$:\n",
    "\n",
    "$$VIF_j=\\frac{1}{1-R_j^2}.$$\n",
    "\n",
    "Near-collinearity makes separate coefficients unstable. Exact collinearity prevents unique identification. Prediction may remain stable within the observed feature region even while coefficients change sharply."
   ]
  },
  {
   "cell_type": "code",
   "id": "fd46c5b4",
   "metadata": {},
   "source": [
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "vif_matrix = design.loc[train_case]\n",
    "vif_values = pd.Series({name: variance_inflation_factor(vif_matrix.to_numpy(), j)\n",
    "                       for j, name in enumerate(vif_matrix.columns) if name != \"const\"})\n",
    "print(vif_values.round(2))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "2f04f31f",
   "metadata": {},
   "source": [
    "A VIF cutoff is a diagnostic convention, rather than a theorem that a feature should be deleted. Consider interpretation, redundancy, sample size and validation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d760ecda",
   "metadata": {},
   "source": [
    "### Ridge and lasso trade variance for bias\n",
    "\n",
    "Ridge minimises $SSE+\\lambda\\sum_j\\beta_j^2$. Lasso uses $SSE+\\lambda\\sum_j|\\beta_j|$ and can set some coefficients to zero. The intercept is normally unpenalised.\n",
    "\n",
    "Scaling matters because penalties act on coefficient magnitudes. Fit the scaler on training rows only. Choose penalty strength on validation data or time-aware folds."
   ]
  },
  {
   "cell_type": "code",
   "id": "8ad180e5",
   "metadata": {},
   "source": [
    "ridge_errors = {}\n",
    "for penalty in [0.1, 1.0, 10.0, 100.0]:\n",
    "    candidate = make_pipeline(StandardScaler(), Ridge(alpha=penalty))\n",
    "    candidate.fit(X_case.loc[train_case], y_case.loc[train_case])\n",
    "    estimate = candidate.predict(X_case.loc[valid_case])\n",
    "    ridge_errors[penalty] = np.sqrt(mean_squared_error(y_case.loc[valid_case], estimate))\n",
    "ridge_alpha = min(ridge_errors, key=ridge_errors.get)\n",
    "print(\"validation RMSE by alpha:\", pd.Series(ridge_errors).round(4).to_dict())\n",
    "print(\"frozen alpha:\", ridge_alpha)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "7b408442",
   "metadata": {},
   "source": [
    "### Bias, variance and the role of the untouched test\n",
    "\n",
    "For squared prediction loss, expected error can be decomposed into irreducible noise, squared bias and variance under a specified data-generating model.\n",
    "\n",
    "More flexible fits often reduce training bias while increasing estimation variance. Regularisation and averaging can reduce variance. The trade-off depends on the sample and model, rather than a universal algorithm ranking.\n",
    "\n",
    "**Diagram concepts:** Training: estimate parameters and preprocessing; Validation: choose a declared finite set of rules; Freeze features, model settings and thresholds; Test: evaluate once and report limitations; New ideas need a new evaluation protocol."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fc2b0ab",
   "metadata": {},
   "source": [
    "### Self-check: selection does not establish causation\n",
    "\n",
    "**Predict:** Which is a defensible model-selection claim?\n",
    "\n",
    "- The subset with the smallest p-value must have the best future return\n",
    "- A higher training R\u00b2 guarantees lower future RMSE\n",
    "- Validation can guide a fixed search, but repeated tuning consumes its independence and coefficients still need identification for causal claims\n",
    "- A large VIF proves every prediction is wrong\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Validation can guide a fixed search, but repeated tuning consumes its independence and coefficients still need identification for causal claims. Validation guides a declared search. Selection uncertainty, generalisation and causal identification remain distinct issues.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91d62b6d",
   "metadata": {},
   "source": [
    "### A coefficient interval is not a future-return interval\n",
    "\n",
    "**Predict:** What does a coefficient confidence interval describe?\n",
    "\n",
    "- A band containing 95% of future stock returns\n",
    "- Uncertainty in a specified coefficient under the estimation and sampling assumptions\n",
    "- The probability that every strategy trade succeeds\n",
    "- The maximum possible drawdown\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Uncertainty in a specified coefficient under the estimation and sampling assumptions. Coefficient uncertainty concerns the estimated model parameter. A future-outcome prediction interval also includes outcome noise.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "addf1e21",
   "metadata": {},
   "source": [
    "### VIF diagnoses coefficient redundancy\n",
    "\n",
    "**Predict:** A large VIF most directly indicates what?\n",
    "\n",
    "- A proven causal relationship\n",
    "- Strong linear predictability of one regressor from the others\n",
    "- Guaranteed poor future portfolio performance\n",
    "- A calibrated classifier\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Strong linear predictability of one regressor from the others. VIF reflects collinearity in X and potential inflation of coefficient variance. It does not independently determine prediction quality.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "039bc727",
   "metadata": {},
   "source": [
    "### What you discovered\n",
    "\n",
    "- Match uncertainty estimates and joint tests to the sampling structure.\n",
    "- Diagnose residuals, interactions and collinearity.\n",
    "- Record subset and penalty searches using validation.\n",
    "- Freeze the selected rule before examining test outcomes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1ae47c0",
   "metadata": {},
   "source": [
    "### Nonlinear Regression and Loss Classification: Trees, Ensembles and Logistic Models\n",
    "\n",
    "Study the definitions, then predict, run, modify and explain the worked results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cec777ec",
   "metadata": {},
   "source": [
    "### Rebuild the two targets for nonlinear learning\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "a406dbf0",
   "metadata": {},
   "source": [
    "import numpy as np, pandas as pd, matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import Ridge, LogisticRegression\n",
    "from sklearn.metrics import mean_squared_error, roc_auc_score, brier_score_loss\n",
    "case = pd.read_csv(\"https://statpython.pages.dev/data/chapter2_scored.csv\",\n",
    "                   parse_dates=[\"date\", \"label_available_at\"])\n",
    "case = case.sort_values([\"date\", \"ticker\"]).reset_index(drop=True)\n",
    "model_keys = [\"date\", \"ticker\"]\n",
    "assert not case.duplicated(model_keys).any()\n",
    "numeric_features = [\"ret_12_1\", \"rvol_21d\", \"momentum_sector_z\"]\n",
    "sector_dummies = pd.get_dummies(case.sector, prefix=\"sector\", drop_first=True, dtype=float)\n",
    "X_case = pd.concat([case[numeric_features], sector_dummies], axis=1)\n",
    "assert np.isfinite(X_case.to_numpy()).all()\n",
    "y_case = case.ret_next.where(case.label_valid)\n",
    "loss_event = y_case.le(-0.05).where(case.label_valid)\n",
    "train_case = case.date.lt(\"2021-01-01\") & case.label_available_at.lt(\"2021-01-01\") & case.label_valid\n",
    "valid_case = (case.date.ge(\"2021-01-01\") & case.date.lt(\"2022-01-01\")\n",
    "              & case.label_available_at.lt(\"2022-01-01\") & case.label_valid)\n",
    "test_case = case.date.ge(\"2022-01-01\") & case.label_valid\n",
    "print(\"train / validation / test / pending:\",\n",
    "      int(train_case.sum()), int(valid_case.sum()), int(test_case.sum()), int((~case.label_valid).sum()))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "51324f85",
   "metadata": {},
   "source": [
    "Regression targets a continuous return. Classification targets the fixed 5% loss event. Keep the split and candidate spine unchanged when comparing models."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "044fbebc",
   "metadata": {},
   "source": [
    "### A regression tree chooses a split by loss reduction\n",
    "\n",
    "For a candidate split s, minimise the two child sums of squared errors:\n",
    "\n",
    "$$\\sum_{i\\in L(s)}(y_i-\\bar y_L)^2+\\sum_{i\\in R(s)}(y_i-\\bar y_R)^2.$$\n",
    "\n",
    "Each leaf predicts its training mean. Recursive splits capture thresholds and interactions automatically, but make piecewise-constant predictions.\n",
    "\n",
    "**Diagram concepts:** Is volatility below the fitted threshold?; Is momentum below a second threshold?; Leaf: mean training return; Leaf: mean training return; Leaf: mean training return.\n",
    "\n",
    "This diagram illustrates the mechanism. The fitted thresholds below come from training data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70c15f14",
   "metadata": {},
   "source": [
    "### Fit a shallow tree and inspect its rule\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "f85e106b",
   "metadata": {},
   "source": [
    "from sklearn.tree import DecisionTreeRegressor, export_text\n",
    "tree_case = DecisionTreeRegressor(max_depth=2, min_samples_leaf=15, random_state=5650)\n",
    "tree_case.fit(X_case.loc[train_case], y_case.loc[train_case])\n",
    "print(export_text(tree_case, feature_names=list(X_case.columns)))\n",
    "print(\"validation RMSE:\", round(np.sqrt(mean_squared_error(\n",
    "    y_case.loc[valid_case], tree_case.predict(X_case.loc[valid_case]))), 4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "2887a38e",
   "metadata": {},
   "source": [
    "Read each path as a sequence of conditional comparisons. The leaf\u2019s predicted return comes from outcomes of training rows reaching that leaf.\n",
    "\n",
    "**Modify:** Increase depth or reduce minimum leaf size. Compare training and validation error; do not select settings using test returns."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca478890",
   "metadata": {},
   "source": [
    "### Classification trees use impurity or log loss\n",
    "\n",
    "For class proportions $p_k$ in a node, Gini impurity is $1-\\sum_k p_k^2$ and entropy is $-\\sum_k p_k\\log p_k$.\n",
    "\n",
    "A split maximises the weighted reduction in impurity. A leaf estimates class probability from its training observations. A pure leaf can be a tiny, unreliable sample.\n",
    "\n",
    "Classification impurity and regression SSE answer different loss functions. Neither automatically reflects the economic cost of a missed loss event.\n",
    "\n",
    "**Check:** Why can a deep classifier attain perfect training accuracy yet give poor probability forecasts?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74c0d5c4",
   "metadata": {},
   "source": [
    "### Forests average; boosting adds corrective learners\n",
    "\n",
    "**Bagging** fits models on bootstrap samples and averages their predictions. A random forest also samples features at candidate splits, helping decorrelate trees. Averaging reduces variance most when errors are not perfectly correlated.\n",
    "\n",
    "**Gradient boosting** adds a small learner in the direction that reduces the chosen loss. Under squared error, the negative gradient is the current residual.\n",
    "\n",
    "$$F_m(x)=F_{m-1}(x)+\\eta h_m(x).$$\n",
    "\n",
    "Learning rate \u03b7, tree depth, leaf size and number of stages jointly control flexibility. Tune them on time-aware validation, rather than assuming one ensemble dominates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f2ff173",
   "metadata": {},
   "source": [
    "### Compare regression learners on the same validation rows\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "6e8e83b1",
   "metadata": {},
   "source": [
    "from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor\n",
    "regression_models = {\n",
    "    \"tree\": DecisionTreeRegressor(max_depth=2, min_samples_leaf=15, random_state=5650),\n",
    "    \"forest\": RandomForestRegressor(n_estimators=60, max_depth=3, min_samples_leaf=10, random_state=5650, n_jobs=1),\n",
    "    \"boosting\": GradientBoostingRegressor(n_estimators=60, max_depth=1, learning_rate=0.05, random_state=5650)}\n",
    "regression_scores = []\n",
    "for name, learner in regression_models.items():\n",
    "    learner.fit(X_case.loc[train_case], y_case.loc[train_case])\n",
    "    prediction = learner.predict(X_case.loc[valid_case])\n",
    "    regression_scores.append({\"model\": name, \"validation_rmse\": np.sqrt(mean_squared_error(y_case.loc[valid_case], prediction))})\n",
    "print(pd.DataFrame(regression_scores).round(4))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "cafd5143",
   "metadata": {},
   "source": [
    "These settings are a small fixed demonstration, not an exhaustive optimisation. Compare with the OLS and constant baselines from \u00a74.1. More models tried means more selection uncertainty."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b5d5f35",
   "metadata": {},
   "source": [
    "### Importance describes a model\u2019s use of a feature\n",
    "\n",
    "Impurity importance totals the training loss reductions credited to each feature. It can favour features with many possible split points.\n",
    "\n",
    "Permutation importance measures how a declared evaluation score changes when a feature is shuffled. Correlated substitutes can obscure the contribution of either variable. Arbitrary shuffling can create implausible market observations or break time structure.\n",
    "\n",
    "Neither measure is a causal effect. Report the model, sample, score and perturbation scheme. A stable predictive feature still needs an availability and implementation check."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23ca512a",
   "metadata": {},
   "source": [
    "### Logistic regression models an event probability\n",
    "\n",
    "$$p_i=P(Y_i\\le-0.05\\mid X_i),\\qquad\n",
    "\\log\\frac{p_i}{1-p_i}=\\beta_0+X_i^\\top\\beta,\\qquad\n",
    "p_i=\\frac1{1+e^{-(\\beta_0+X_i^\\top\\beta)}}.$$\n",
    "\n",
    "A coefficient changes log odds. $e^{\\beta_j}$ is the conditional odds multiplier for a one-unit feature increase. Odds and probabilities differ.\n",
    "\n",
    "Logit estimates coefficients by maximising Bernoulli likelihood, equivalently minimising log loss. Regularisation adds a penalty. A probability model requires calibration checks as well as ranking scores."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b0e4e35",
   "metadata": {},
   "source": [
    "### Fit the loss classifier with training-only scaling\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "b3d2b7d4",
   "metadata": {},
   "source": [
    "risk_logit = make_pipeline(StandardScaler(), LogisticRegression(C=1, max_iter=1000))\n",
    "risk_logit.fit(X_case.loc[train_case], loss_event.loc[train_case].astype(int))\n",
    "validation_probability = risk_logit.predict_proba(X_case.loc[valid_case])[:, 1]\n",
    "print(\"event rate train / validation:\",\n",
    "      round(loss_event.loc[train_case].mean(), 3), round(loss_event.loc[valid_case].mean(), 3))\n",
    "print(\"validation AUC:\", round(roc_auc_score(loss_event.loc[valid_case].astype(int), validation_probability), 3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "7c8d8f6f",
   "metadata": {},
   "source": [
    "The event rates are about **21.7% in training** and **16.2% in validation**. Validation AUC is about **0.544**.\n",
    "\n",
    "`StandardScaler` inside the pipeline learns its means and scales only during training. A numeric coefficient here measures a one-training-standard-deviation change; scaling dummy variables also changes their coefficient units.\n",
    "\n",
    "**Check:** AUC above 0.5 in one validation sample does not establish a calibrated probability or a net economic benefit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65c60b37",
   "metadata": {},
   "source": [
    "### Confusion counts depend on a threshold\n",
    "\n",
    "Predict an adverse event when $p\\ge c$.\n",
    "\n",
    "| Actual / prediction | Predict event | Predict no event |\n",
    "|---|---|---|\n",
    "| Event | True positive (TP) | False negative (FN) |\n",
    "| No event | False positive (FP) | True negative (TN) |\n",
    "\n",
    "Precision = TP/(TP+FP), recall = TP/(TP+FN), specificity = TN/(TN+FP). Accuracy can look high when events are rare.\n",
    "\n",
    "Lowering c tends to increase both true positives and false positives. State how undefined ratios are handled when a denominator is zero."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9bd0d6b",
   "metadata": {},
   "source": [
    "### Ranking and probability quality need separate scores\n",
    "\n",
    "**ROC AUC** measures ranking across event/non-event pairs. ROC traces recall against false-positive rate as the threshold changes. Precision\u2013recall curves are useful when positives are rare, and depend strongly on prevalence.\n",
    "\n",
    "**Brier score** is $n^{-1}\\sum_i(p_i-y_i)^2$. **Log loss** is $-n^{-1}\\sum_i[y_i\\log p_i+(1-y_i)\\log(1-p_i)]$.\n",
    "\n",
    "Calibration compares predicted probabilities with observed event frequencies. Small bins and dependent observations produce uncertain calibration estimates.\n",
    "\n",
    "Use a training-prevalence baseline. AUC ignores threshold costs and probability calibration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1243553",
   "metadata": {},
   "source": [
    "### Interactive threshold and error costs\n",
    "\n",
    "Interactive diagram: open the corresponding [online slide deck](https://statpython.pages.dev/). Change the control, predict the outcome and explain it.\n",
    "\n",
    "**Predict, then move:** Which cases change at thresholds 0.20 and 0.80? Explain why the minimum toy cost does not identify the correct threshold for the real validation sample."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "154b8e8b",
   "metadata": {},
   "source": [
    "### Choose the decision threshold on validation\n",
    "\n",
    "Suppose a false alarm costs **1** and a missed adverse event costs **4**, in declared abstract decision-cost units. With calibrated probabilities and constant costs, the expected-cost threshold is $1/(1+4)=0.20$.\n",
    "\n",
    "The model may be miscalibrated. This exercise also evaluates a fixed threshold grid on validation only."
   ]
  },
  {
   "cell_type": "code",
   "id": "1d7be12b",
   "metadata": {},
   "source": [
    "threshold_records = []\n",
    "actual_validation = loss_event.loc[valid_case].astype(bool).to_numpy()\n",
    "for cutoff in [0.20, 0.35, 0.50, 0.65]:\n",
    "    alarm = validation_probability >= cutoff\n",
    "    false_positive = int((alarm & ~actual_validation).sum())\n",
    "    false_negative = int((~alarm & actual_validation).sum())\n",
    "    threshold_records.append({\"cutoff\": cutoff, \"FP\": false_positive, \"FN\": false_negative,\n",
    "                              \"cost\": false_positive + 4 * false_negative})\n",
    "threshold_report = pd.DataFrame(threshold_records)\n",
    "risk_cutoff = float(threshold_report.loc[threshold_report.cost.idxmin(), \"cutoff\"])\n",
    "print(threshold_report)\n",
    "print(\"frozen cutoff:\", risk_cutoff)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "0c5c41df",
   "metadata": {},
   "source": [
    "The grid selects **0.65**, costing **56** units versus 64 for never alarming and 83 for always alarming. This is a validation result, with tuning uncertainty. These abstract costs are not realised trading P&L."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22593014",
   "metadata": {},
   "source": [
    "### Open the held-out test after freezing the rule\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "9e36cbf9",
   "metadata": {},
   "source": [
    "from sklearn.metrics import confusion_matrix, log_loss\n",
    "test_probability = risk_logit.predict_proba(X_case.loc[test_case])[:, 1]\n",
    "test_actual = loss_event.loc[test_case].astype(int)\n",
    "constant_probability = np.repeat(loss_event.loc[train_case].mean(), len(test_actual))\n",
    "print(\"test AUC:\", round(roc_auc_score(test_actual, test_probability), 3))\n",
    "print(\"test Brier / baseline Brier:\", round(brier_score_loss(test_actual, test_probability), 4),\n",
    "      round(brier_score_loss(test_actual, constant_probability), 4))\n",
    "print(\"test log loss:\", round(log_loss(test_actual, test_probability), 4))\n",
    "print(\"confusion matrix, rows actual 0/1, columns predicted 0/1:\\n\",\n",
    "      confusion_matrix(test_actual, test_probability >= risk_cutoff, labels=[0, 1]))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "e3049c30",
   "metadata": {},
   "source": [
    "Test AUC is about **0.431**, and Brier score about **0.2176**. This weak result is useful evidence against assuming the fitted risk model generalises.\n",
    "\n",
    "Inspecting this test to invent another feature or threshold would turn it into development data. A revised model needs a new evaluation protocol."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98ba5d82",
   "metadata": {},
   "source": [
    "### Self-check: probabilities and decisions\n",
    "\n",
    "**Predict:** Which claim is supported by a validation AUC of 0.544?\n",
    "\n",
    "- The model predicts a 54.4% return\n",
    "- Its default cutoff must be 0.5\n",
    "- The score showed modest event ranking in that sample, while calibration, costs and future performance need separate checks\n",
    "- The model should automatically block every trade\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "The score showed modest event ranking in that sample, while calibration, costs and future performance need separate checks. AUC measures event ranking in a stated sample. It does not supply return, calibration or a decision threshold.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f719f38",
   "metadata": {},
   "source": [
    "### Repair the adverse-event definition\n",
    "\n",
    "The event is a next-month simple return of \u22125% or worse. Fix the comparison without inventing labels for pending rows."
   ]
  },
  {
   "cell_type": "code",
   "id": "05bd020b",
   "metadata": {},
   "source": [
    "event_check = y_case.le(-0.05).where(case.label_valid)  # fix ge / le\n",
    "print(event_check.value_counts(dropna=False))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "a33837ce",
   "metadata": {},
   "source": [
    "### Log odds differ from probability\n",
    "\n",
    "**Predict:** A logit coefficient is log(2) per unit. What doubles for a one-unit increase, holding other inputs fixed?\n",
    "\n",
    "- The event probability at every starting value\n",
    "- The conditional odds\n",
    "- The expected next-month return\n",
    "- The classifier AUC\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "The conditional odds. Exponentiating a logit coefficient gives an odds multiplier. A probability change depends on the starting probability.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "974972db",
   "metadata": {},
   "source": [
    "### Ensemble methods change the estimation trade-off\n",
    "\n",
    "**Predict:** Under squared-error gradient boosting, the next learner fits what?\n",
    "\n",
    "- The current residual direction\n",
    "- The future test returns before they arrive\n",
    "- Only the largest feature values\n",
    "- A causal treatment effect\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "The current residual direction. The negative gradient of squared loss is proportional to observed minus current fitted value. General losses have different pseudo-residuals.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79e715c4",
   "metadata": {},
   "source": [
    "### What you discovered\n",
    "\n",
    "- Explain tree splitting, forest averaging and boosting updates.\n",
    "- Interpret logistic probabilities and odds.\n",
    "- Separate ranking, calibration and decision costs.\n",
    "- Report the weak fixed-classifier test result without further test tuning."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe5c6a11",
   "metadata": {},
   "source": [
    "### Causal Analysis: Prediction vs Causation, Confounders, Difference-in-Differences, Instrumental Variables\n",
    "\n",
    "Everything so far answered \"given $x$, what is $y$?\". A manager asks \"if I *change* $x$, what happens to $y$?\". You will build a case where the best predictor is a useless lever, fix it with a control, and then meet two designs for when the control is not observed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c221b00",
   "metadata": {},
   "source": [
    "### Two questions that look alike\n",
    "\n",
    "**Predict:** \"Traders who use a risk dashboard have higher risk-control scores.\" Which action needs a <em>causal</em> answer?\n",
    "\n",
    "- Use adoption as a predictor of a later risk-control score\n",
    "- Report the adoption\u2013score association\n",
    "- Require every trader to adopt the dashboard, expecting the observed improvement\n",
    "- Report the correlation in the annual review\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Require every trader to adopt the dashboard, expecting the observed improvement. A prediction only needs association: any variable correlated with y helps. A causal claim needs that nothing else drives both x and y. The same regression coefficient serves the first question and, in general, answers the second wrongly.\n",
    "\n",
    "</details>\n",
    "\n",
    "Potential outcomes: each trader has $Y(1)$ with the tool and $Y(0)$ without; we see one. The observed gap is $\\underbrace{E[Y(1)-Y(0)\\mid T=1]}_{\\text{causal}} + \\underbrace{E[Y(0)\\mid T=1] - E[Y(0)\\mid T=0]}_{\\text{selection bias}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40c1c84b",
   "metadata": {},
   "source": [
    "### A great predictor that is a useless lever\n",
    "\n",
    "Simulate a hidden trait $U$ (risk discipline) that drives *both* tool adoption $T$ and risk-control score $Y$. By construction the tool has **zero** effect."
   ]
  },
  {
   "cell_type": "code",
   "id": "c3ad7cf6",
   "metadata": {},
   "source": [
    "import numpy as np, pandas as pd, statsmodels.api as sm, matplotlib.pyplot as plt\n",
    "rng = np.random.default_rng(44); n = 2000\n",
    "U = rng.normal(size=n)                                # discipline: not in the naive analyst's data\n",
    "T = (U + rng.normal(size=n) > 0.5).astype(int)        # disciplined traders adopt the tool\n",
    "Y = 2.0 * U + 0.0 * T + rng.normal(size=n)            # risk-control score: the tool does NOTHING\n",
    "naive = sm.OLS(Y, sm.add_constant(T)).fit()\n",
    "print(\"share using the tool:\", round(T.mean(), 3), \"  corr(T, Y) =\", round(np.corrcoef(T, Y)[0, 1], 3))\n",
    "print(\"naive coefficient on T =\", round(naive.params[1], 3), \"  t =\", round(naive.tvalues[1], 1), \"  R2 =\", round(naive.rsquared, 3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "e82a2c30",
   "metadata": {},
   "source": [
    "Coefficient **2.17** with t = 24 and R\u00b2 = 0.23. As a *predictor* of risk-control score, tool usage is excellent. As a *policy*, pushing the tool would move risk-control score by exactly **0**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da8b9994",
   "metadata": {},
   "source": [
    "### Regression adjustment: control for the confounder\n",
    "\n",
    "**Predict:** Add $U$ as a second regressor. The coefficient on $T$ will:\n",
    "\n",
    "- Stay near 2.17 \u2014 $T$ is still strongly correlated with $Y$\n",
    "- Collapse toward 0 \u2014 the true effect\n",
    "- Flip to about \u22122\n",
    "- Become undefined because $T$ and $U$ are correlated\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Collapse toward 0 \u2014 the true effect. Once U is in the regression, the comparison is between tool users and non-users with the same discipline. The backdoor path T <- U -> Y is closed and the coefficient on T collapses to the true effect, zero, up to noise.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "id": "33b83813",
   "metadata": {},
   "source": [
    "adj = sm.OLS(Y, sm.add_constant(np.column_stack([T, U]))).fit()\n",
    "print(\"adjusted coefficient on T =\", round(adj.params[1], 3), \"  on U =\", round(adj.params[2], 3))\n",
    "Y2 = 2.0 * U + 0.5 * T + rng.normal(size=n)           # same world, but now the tool adds 0.5\n",
    "print(\"true effect 0.5:  naive\", round(sm.OLS(Y2, sm.add_constant(T)).fit().params[1], 3),\n",
    "      \"  adjusted\", round(sm.OLS(Y2, sm.add_constant(np.column_stack([T, U]))).fit().params[1], 3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "64d33922",
   "metadata": {},
   "source": [
    "Adjusted: **0.066** (truth 0). With a real effect of 0.5, naive says 2.67, adjusted says **0.45**. This works only because $U$ was *measured*. Draw the graph $T \\leftarrow U \\rightarrow Y$ first; choose controls from the causal graph. Conditioning on mediators or colliders can change or bias the target effect."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65214696",
   "metadata": {},
   "source": [
    "### When the confounder is not in your data\n",
    "\n",
    "**Diagram concepts:** Unobserved risk discipline; Dashboard adoption; Risk-control score.\n",
    "\n",
    "- Regression adjustment cannot close a backdoor through a variable you do not have. No amount of data on $T$ and $Y$ alone recovers the effect.\n",
    "- **Randomise** if you can: an A/B test makes $T$ independent of $U$ and the difference of means is causal.\n",
    "- If you cannot, two observational designs replace the missing control with a *structure*:\n",
    "  - **Difference-in-differences** \u2014 a control group that shares the same time trend.\n",
    "  - **Instrumental variables** \u2014 a source of variation in $T$ that has nothing to do with $U$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96ec9ec6",
   "metadata": {},
   "source": [
    "### Difference-in-differences: the gap between gaps\n",
    "\n",
    "Treated group gets a policy between period 0 and 1; the control group does not.\n",
    "$$\\hat\\tau_{\\text{DiD}} = (\\bar Y^{\\text{tr}}_{1} - \\bar Y^{\\text{tr}}_{0}) - (\\bar Y^{\\text{ctl}}_{1} - \\bar Y^{\\text{ctl}}_{0})$$\n",
    "\n",
    "**Predict:** The identifying assumption of DiD is:\n",
    "\n",
    "- The two groups have the same level of $Y$ before treatment\n",
    "- Treatment was randomly assigned\n",
    "- Without treatment, both groups would have followed parallel trends\n",
    "- The control group is larger than the treated group\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Without treatment, both groups would have followed parallel trends. The control group's change stands in for what the treated group would have done without the policy. That only works if, absent treatment, both groups would have moved in parallel. Level differences are fine; different trends are fatal.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2638e60",
   "metadata": {},
   "source": [
    "### Simulate a two-group, two-period panel\n",
    "\n",
    "Treated group starts 1.0 higher, everyone drifts up 0.5, and the policy adds **0.8** to the treated group in period 1."
   ]
  },
  {
   "cell_type": "code",
   "id": "9da25f25",
   "metadata": {},
   "source": [
    "rng = np.random.default_rng(104); rows = []\n",
    "for grp in [0, 1]:\n",
    "    for per in [0, 1]:\n",
    "        y = 5 + 1.0 * grp + 0.5 * per + 0.8 * grp * per + rng.normal(scale=0.5, size=400)\n",
    "        rows.append(pd.DataFrame({\"group\": grp, \"period\": per, \"y\": y}))\n",
    "panel = pd.concat(rows, ignore_index=True)\n",
    "means = panel.groupby([\"group\", \"period\"])[\"y\"].mean().unstack()\n",
    "print(means.round(3))\n",
    "print(\"naive post-period gap      :\", round(means.loc[1, 1] - means.loc[0, 1], 3))\n",
    "print(\"naive before-after, treated:\", round(means.loc[1, 1] - means.loc[1, 0], 3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "8b17b9d6",
   "metadata": {},
   "source": [
    "Both naive numbers are wrong: the post gap (1.85) includes the pre-existing level difference; the before-after (1.31) includes the common drift."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80218337",
   "metadata": {},
   "source": [
    "### Your turn: compute the DiD by hand\n",
    "\n",
    "`did` currently holds the treated group's before-after change. Subtract the control group's change so that `did` is the difference-in-differences."
   ]
  },
  {
   "cell_type": "code",
   "id": "d513d3b3",
   "metadata": {},
   "source": [
    "did = (means.loc[1, 1] - means.loc[1, 0]) - (means.loc[0, 1] - means.loc[0, 0])  # solution: subtract (means.loc[0, 1] - means.loc[0, 0])\n",
    "print(\"DiD estimate:\", round(did, 3), \"  (true effect 0.8)\")"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "id": "bbd02043",
   "metadata": {},
   "source": [
    "X = sm.add_constant(pd.DataFrame({\"group\": panel[\"group\"], \"period\": panel[\"period\"],\n",
    "                                  \"group_x_period\": panel[\"group\"] * panel[\"period\"]}))\n",
    "r = sm.OLS(panel[\"y\"], X).fit()\n",
    "print(pd.DataFrame({\"coef\": r.params, \"se\": r.bse}).round(3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "81b0973a",
   "metadata": {},
   "source": [
    "The regression $y = \\alpha + \\beta\\,\\text{group} + \\gamma\\,\\text{period} + \\delta\\,(\\text{group} \\times \\text{period})$ gives $\\hat\\delta$ = **0.861** (se 0.051) \u2014 identical to the hand calculation, now with a standard error. Counterfactual parallel trends cannot be verified directly. Pre-treatment trends can challenge the assumption but cannot prove it. This independently simulated cross-section uses an ordinary OLS standard error; real repeated-unit designs need design-appropriate inference."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95438663",
   "metadata": {},
   "source": [
    "### Instrumental variables: borrow a lottery\n",
    "\n",
    "An instrument $Z$ moves $T$ but touches $Y$ only *through* $T$. Three conditions: **relevance** ($Z \\to T$), **exclusion** (no $Z \\to Y$ path except via $T$), **independence** ($Z \\perp U$).\n",
    "\n",
    "**Predict:** Which condition can <em>never</em> be verified from the data alone?\n",
    "\n",
    "- Relevance \u2014 $Z$ predicts $T$\n",
    "- Exclusion \u2014 $Z$ affects $Y$ only through $T$\n",
    "- That $n$ is large enough\n",
    "- That $T$ is binary\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Exclusion \u2014 $Z$ affects $Y$ only through $T$. Relevance is testable (first-stage F). Independence is often argued from the design (a lottery). Exclusion, that the instrument has no other channel to Y, cannot be tested from data; it has to be defended on substantive grounds and is where most IV papers are attacked.\n",
    "\n",
    "</details>\n",
    "\n",
    "Classic instruments: distance to college for schooling (Card 1995), draft-lottery number for military service (Angrist 1990), rainfall for agricultural income. Two-stage least squares: regress $T$ on $Z$, then $Y$ on $\\hat T$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b81d71a",
   "metadata": {},
   "source": [
    "### 2SLS by hand with numpy\n",
    "\n",
    "Hidden $U$ drives $T$ and $Y$; true effect of $T$ on $Y$ is **1.0**; $Z$ is a clean shock to $T$."
   ]
  },
  {
   "cell_type": "code",
   "id": "8e0b23a3",
   "metadata": {},
   "source": [
    "rng = np.random.default_rng(102); n = 2000\n",
    "U, Z = rng.normal(size=n), rng.normal(size=n)\n",
    "T = 0.5 * Z + 0.8 * U + rng.normal(scale=0.3, size=n)\n",
    "Y = 1.0 * T + 1.2 * U + rng.normal(scale=0.5, size=n)\n",
    "one = np.ones(n)\n",
    "b_ols = np.linalg.lstsq(np.column_stack([one, T]), Y, rcond=None)[0]\n",
    "g1 = np.linalg.lstsq(np.column_stack([one, Z]), T, rcond=None)[0]         # stage 1: T on Z\n",
    "T_hat = g1[0] + g1[1] * Z\n",
    "b_iv = np.linalg.lstsq(np.column_stack([one, T_hat]), Y, rcond=None)[0]   # stage 2: Y on T_hat\n",
    "print(\"OLS slope  :\", round(b_ols[1], 3), \"  (true 1.0)\")\n",
    "print(\"2SLS slope :\", round(b_iv[1], 3))\n",
    "print(\"cov(Z,Y)/cov(Z,T):\", round(np.cov(Z, Y)[0, 1] / np.cov(Z, T)[0, 1], 3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "358d82b4",
   "metadata": {},
   "source": [
    "OLS says **2.013** \u2014 twice the truth, because $U$ pushes $T$ and $Y$ the same way. 2SLS uses only the part of $T$ that $Z$ explains and lands on **1.027**. The one-instrument case reduces to the ratio $\\text{cov}(Z,Y)/\\text{cov}(Z,T)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b378b44",
   "metadata": {},
   "source": [
    "### Weak instruments: when the lottery barely moves T\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "2e8f6027",
   "metadata": {},
   "source": [
    "print(\"first-stage F, strong Z:\", round(sm.OLS(T, sm.add_constant(Z)).fit().fvalue, 1))\n",
    "T_w = 0.05 * Z + 0.8 * U + rng.normal(scale=0.3, size=n)                 # Z now barely moves T\n",
    "Y_w = 1.0 * T_w + 1.2 * U + rng.normal(scale=0.5, size=n)\n",
    "g1 = np.linalg.lstsq(np.column_stack([one, Z]), T_w, rcond=None)[0]\n",
    "b_w = np.linalg.lstsq(np.column_stack([one, g1[0] + g1[1] * Z]), Y_w, rcond=None)[0]\n",
    "print(\"first-stage F, weak Z  :\", round(sm.OLS(T_w, sm.add_constant(Z)).fit().fvalue, 1), \"  2SLS slope:\", round(b_w[1], 3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "8fe3ee71",
   "metadata": {},
   "source": [
    "Strong instrument: F = **637**. Weak one: F = **7.0** and the 2SLS estimate drifts to 0.72 and weak identification. The first-stage F > 10 rule is a context-dependent heuristic, not proof of valid identification. The hand calculation reports point estimates only; second-stage OLS standard errors are not valid IV standard errors.\n",
    "\n",
    "| Design | Replaces the missing control with | Fails when |\n",
    "|---|---|---|\n",
    "| Regression adjustment | the measured confounder | a confounder is unmeasured |\n",
    "| Difference-in-differences | a control group's time trend | trends are not parallel |\n",
    "| Instrumental variables | an exogenous shock to $T$ | $Z$ is weak or has its own path to $Y$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21528bae",
   "metadata": {},
   "source": [
    "### What you discovered\n",
    "\n",
    "- A predictor with t = 24 and R\u00b2 = 0.23 had a causal effect of **zero**; the correlation came entirely from the hidden trait behind both $T$ and $Y$.\n",
    "\n",
    "- Regression adjustment recovers the effect (0.066 \u2248 0; 0.45 \u2248 0.5) \u2014 but only for confounders you measured.\n",
    "\n",
    "- DiD nets out level differences and common trends: 0.861 for a true 0.8. Its price is the parallel-trends assumption.\n",
    "\n",
    "- 2SLS turned an OLS estimate of 2.01 into 1.03 using a clean instrument; with F = 7 the instrument was too weak to trust.\n",
    "\n",
    "- Never report a causal number without naming the assumption it rests on."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a96c14",
   "metadata": {},
   "source": [
    "### Cross-Sectional Attention Features: Transformers for the Stock Cross-Section\n",
    "\n",
    "Every month a fund sees the same object: a matrix of stocks by characteristics, and next month's returns as the target. You will standardise it without peeking at the future, replace the sector as \"peer group\" with one attention layer written in numpy, and put both through the walk-forward protocol of \u00a74.3 \u2014 and discover that the protocol, not the architecture, is the lesson."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9972eca",
   "metadata": {},
   "source": [
    "### One month, one matrix\n",
    "\n",
    "At month-end $t$, $X_t \\in \\mathbb{R}^{N \\times K}$ holds $N$ stocks \u00d7 $K$ characteristics and $y_t \\in \\mathbb{R}^N$ holds next month's returns. The prediction problem is *cross-sectional*: rank the $N$ stocks, month by month."
   ]
  },
  {
   "cell_type": "code",
   "id": "564de20c",
   "metadata": {},
   "source": [
    "import numpy as np, pandas as pd\n",
    "xs = pd.read_csv(\"https://statpython.pages.dev/data/xs_panel_monthly.csv\", parse_dates=[\"date\"])\n",
    "chars = [\"ret_1_0\", \"ret_12_1\", \"rvol_21d\", \"log_dvol_126d\", \"prc_highprc_252d\"]\n",
    "print(xs.shape, \" months:\", xs[\"date\"].nunique(), \" stocks:\", xs[\"ticker\"].nunique(), \" sectors:\", xs[\"sector\"].nunique())\n",
    "print(xs[\"date\"].min().date(), \"\u2026\", xs[\"date\"].max().date())\n",
    "a = xs[xs[\"ticker\"] == \"AAPL\"].set_index(\"date\")\n",
    "print(\"AAPL ret_next at 2015-01:\", a.loc[\"2015-01-31\", \"ret_next\"], \" ret_1_0 at 2015-02:\", a.loc[\"2015-02-28\", \"ret_1_0\"])\n",
    "print(xs[chars + [\"ret_next\"]].describe().loc[[\"mean\", \"std\", \"min\", \"max\"]].round(3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "634b4a07",
   "metadata": {},
   "source": [
    "25 335 rows = 119 months \u00d7 213 stocks in 11 sectors, 2015-01 \u2026 2024-11. Five price-based characteristics: last-month return, 12-1 momentum, 21-day realised volatility, log dollar volume, price / 252-day high. `ret_next` at month $t$ is exactly `ret_1_0` at month $t+1$ (0.100776 = 0.100776) \u2014 the target is the *next* row's return."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbb76a97",
   "metadata": {},
   "source": [
    "### Whose 213 stocks are these?\n",
    "\n",
    "The universe is **today's** large caps, tracked back to 2015. Nothing that was large in 2015 and has since failed, shrunk or been delisted is in the file.\n",
    "\n",
    "**Predict:** What does this <em>survivorship bias</em> do to a backtest run on the file?\n",
    "\n",
    "- Nothing \u2014 every price is the real price on that day\n",
    "- Flatters it: firms that failed are absent, so returns are biased upward and losers are under-represented\n",
    "- Hurts it: today's winners were expensive in 2015\n",
    "- Only changes the sector mix\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Flatters it: firms that failed are absent, so returns are biased upward and losers are under-represented. The universe was chosen in 2024 on the condition of having survived and grown. Stocks that collapsed are missing, so the average stock in the file did better than the average stock a 2015 investor could have bought. Any long-only or long-tilted signal inherits that upward bias; a long-short spread partly cancels it, but the losers that would have populated the short leg are also gone.\n",
    "\n",
    "</details>\n",
    "\n",
    "Production rule: build the universe *as of each month* from a point-in-time index membership file. This file is fine for learning the mechanics \u2014 never quote its Sharpe ratio as evidence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d4221e1",
   "metadata": {},
   "source": [
    "### Standardise per snapshot, never across time\n",
    "\n",
    "**Predict:** Why is <code>(rvol_21d \u2212 rvol_21d.mean()) / rvol_21d.std()</code> over the whole 2015\u20132024 file a look-ahead?\n",
    "\n",
    "- Volatility is not normally distributed\n",
    "- The mean is dominated by 2020\n",
    "- The 2015 z-score depends on 2024's mean and std, which a live model in 2015 could not have seen\n",
    "- It is not: the transform is monotone\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "The 2015 z-score depends on 2024's mean and std, which a live model in 2015 could not have seen. A full-sample z-score of the 2015 row uses the mean and standard deviation of 2015\u20132024, including nine years that had not happened. The live model in 2015 would have seen a different number \u2014 the backtest is not reproducing what it could have known. Per-month standardisation uses only the cross-section observed at t.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "id": "56c62301",
   "metadata": {},
   "source": [
    "def zscore(df, by):\n",
    "    g = df.groupby(by)[chars]\n",
    "    return ((df[chars] - g.transform(\"median\")) / g.transform(\"std\")).clip(-3, 3).fillna(0.0)\n",
    "zc = [c + \"_z\" for c in chars]; sc = [c + \"_s\" for c in chars]\n",
    "xs[zc] = zscore(xs, \"date\").values              # market-relative: median / std across all stocks in month t\n",
    "xs[sc] = zscore(xs, [\"date\", \"sector\"]).values  # sector-relative: median / std within sector in month t\n",
    "chk = xs.groupby([\"date\", \"sector\"])[\"ret_12_1_s\"].median().abs().max()\n",
    "print(\"largest |sector median| after standardising:\", round(chk, 6))\n",
    "print(xs[[\"ret_12_1_z\", \"ret_12_1_s\"]].describe().loc[[\"mean\", \"std\", \"min\", \"max\"]].round(2).T)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "8367b4fe",
   "metadata": {},
   "source": [
    "Median and std are computed **within each month** (and, for the `_s` columns, within each month \u00d7 sector \u2014 the within-date transformation rule: a stock is judged against the peers it is compared with on the day). Sector medians are exactly 0 afterwards; values are clipped at \u00b13 to reduce the numerical influence of extremes; clipping does not guarantee robust estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4752b5b1",
   "metadata": {},
   "source": [
    "### Why is a sector median a crude peer group?\n",
    "\n",
    "\"Sector\" is a label assigned by an index provider. In October 2022 NVDA's *behavioural* peers \u2014 high volume, high volatility, far below the 52-week high \u2014 included AMD, but also Netflix, Alibaba and Tesla, none of them \"Technology\".\n",
    "\n",
    "Attention builds the peer group from the data. With queries $Q = XW_q$, keys $K = XW_k$, values $V = XW_v$:\n",
    "\n",
    "$$A = \\text{softmax}_{\\text{rows}}\\!\\left(\\frac{QK^\\top}{\\sqrt{d}\\,\\tau}\\right), \\qquad C = AV$$\n",
    "\n",
    "**Predict:** With $W_q = W_k = I$, the weight stock $i$ gives stock $j$ is proportional to\n",
    "\n",
    "- $1/N$ \u2014 every stock equally\n",
    "- $\\exp(x_i \\cdot x_j / \\sqrt{K})$ \u2014 a similarity kernel in characteristic space\n",
    "- the sector indicator $\\mathbb{1}[s_i = s_j]$\n",
    "- the Euclidean distance $\\lVert x_i - x_j \\rVert$\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "$\\exp(x_i \\cdot x_j / \\sqrt{K})$ \u2014 a similarity kernel in characteristic space. With identity weights the score is the inner product x_i\u00b7x_j/\u221aK and the softmax turns it into exp(x_i\u00b7x_j/\u221aK) normalised over j: a similarity kernel in characteristic space. Stocks that look like i get high weight, whatever sector label they carry.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8111e001",
   "metadata": {},
   "source": [
    "### One attention layer in twelve lines of numpy\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "51fb0c8b",
   "metadata": {},
   "source": [
    "def attention(X, Wq, Wk, Wv, tau=1.0):\n",
    "    Q, K, V = X @ Wq, X @ Wk, X @ Wv                 # queries, keys, values\n",
    "    S = Q @ K.T / np.sqrt(K.shape[1]) / tau          # N \u00d7 N scores\n",
    "    A = np.exp(S - S.max(axis=1, keepdims=True))     # softmax across stocks, row by row\n",
    "    A /= A.sum(axis=1, keepdims=True)\n",
    "    return A, A @ V                                  # weights, context\n",
    "I = np.eye(5)\n",
    "m = xs[xs[\"date\"] == \"2022-10-31\"].reset_index(drop=True)\n",
    "X = m[zc].to_numpy(); i = int(np.where(m[\"ticker\"] == \"NVDA\")[0][0])\n",
    "A, C = attention(X, I, I, I)\n",
    "row = pd.Series(A[i], index=m[\"ticker\"]).sort_values(ascending=False)\n",
    "print(row.head(6).round(3).to_dict(), \" stocks with weight > 1 %:\", int((A[i] > 0.01).sum()))\n",
    "print(m.set_index(\"ticker\").loc[row.index[:6], \"sector\"].to_dict())\n",
    "print(\"NVDA's row of X:\", X[i].round(2), \"  context C:\", C[i].round(2))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "8af0f00b",
   "metadata": {},
   "source": [
    "NVDA's attention row in October 2022: **AMD 0.353**, itself 0.263, **NFLX 0.154**, **BABA 0.077**, **TSLA 0.055**, AMZN 0.028 \u2014 six stocks carry more than 1 %, from three sectors. The context $C$ is the peer-weighted characteristic vector: 12-1 momentum \u22121.49 for NVDA against \u22121.38 for its peers, dollar volume 3.0 against 2.52."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dac51139",
   "metadata": {},
   "source": [
    "### Predict: what happens as the temperature grows?\n",
    "\n",
    "With identity weights and $\\tau \\to \\infty$ every score goes to 0. What does NVDA's context vector collapse to?\n",
    "\n",
    "```{.p-code}\n",
    "A, C = attention(X, I, I, I, tau=1e6)\n",
    "print(C[i].round(2))\n",
    "```\n",
    "[`the cross-sectional mean of X \u2014 every weight is 1/N`]{.p-answer}"
   ]
  },
  {
   "cell_type": "code",
   "id": "1f2979ff",
   "metadata": {},
   "source": [
    "for tau in [0.05, 1, 5, 1e6]:\n",
    "    A, C = attention(X, I, I, I, tau)\n",
    "    print(f\"tau={tau:<8g} self-weight {A[i, i]:.3f}  largest other {np.delete(A[i], i).max():.3f}  context {C[i].round(2)}\")\n",
    "print(\"cross-sectional mean of X:      \", X.mean(axis=0).round(2))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "45df9406",
   "metadata": {},
   "source": [
    "Temperature is the dial between two crude peer groups: $\\tau \\to 0$ attends only to the most similar stock (at $\\tau = 0.05$ AMD carries 0.997 \u2014 NVDA is not even its own nearest point), $\\tau \\to \\infty$ gives the cross-sectional mean (0.02, 0.10, 0.19, 0.00, \u22120.14), while median-centred per-month features need not have exactly zero means. $\\tau = 1$ sits in between: a handful of behavioural peers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61a48b1f",
   "metadata": {},
   "source": [
    "### What Kelly, Kuznetsov, Malamud and Xu actually train\n",
    "\n",
    "**\"Artificial Intelligence Asset Pricing Models\"** (KKMX, 2025) puts a transformer over the stock cross-section, not over words:\n",
    "\n",
    "- Each layer: **linear attention** $A = XW_qW_k^\\top X^\\top$ \u2014 no softmax \u2014 then a **feed-forward block** and a **residual** connection, $X^{(\\ell+1)} = X^{(\\ell)} + \\text{FFN}(A X^{(\\ell)} W_v)$, stacked $L$ times.\n",
    "- The final layer maps each stock's row to a **portfolio weight**; all $W$'s are trained end-to-end to **maximise the Sharpe ratio** of that portfolio, not to minimise a squared error.\n",
    "- Attention lets a stock's weight depend on the whole cross-section \u2014 the \"peer group\" is learned.\n",
    "\n",
    "Weights stay **fixed** ($W = I$; a seeded random $W$ works the same way and is the random-features view). The attention output is used only as **feature engineering**: for every stock, the context $C$ (peer-weighted characteristics) and the deviation $X - C$ (\"how I differ from my attention-weighted peers\"). The learner on top is the boosting machine from \u00a74.3."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90523509",
   "metadata": {},
   "source": [
    "### From attention to features: $C$ and $X - C$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "942dc3d2",
   "metadata": {},
   "source": [
    "cc = [c + \"_c\" for c in chars]; dc = [c + \"_d\" for c in chars]\n",
    "def attention_features(tau=1.0, keys=zc):\n",
    "    parts = []\n",
    "    for d, g in xs.groupby(\"date\"):                  # one attention pass per month\n",
    "        Xk, Xv = g[keys].to_numpy(), g[zc].to_numpy()\n",
    "        S = Xk @ Xk.T / np.sqrt(Xk.shape[1]) / tau\n",
    "        A = np.exp(S - S.max(axis=1, keepdims=True)); A /= A.sum(axis=1, keepdims=True)\n",
    "        C = A @ Xv\n",
    "        parts.append(pd.DataFrame(np.hstack([C, Xv - C]), index=g.index, columns=cc + dc))\n",
    "    return pd.concat(parts)\n",
    "xs[cc + dc] = attention_features()            # aligned by index, not by position\n",
    "nv = xs[(xs[\"ticker\"] == \"NVDA\") & (xs[\"date\"] == \"2022-10-31\")]\n",
    "print(nv[[\"ret_12_1_z\", \"ret_12_1_c\", \"ret_12_1_d\", \"prc_highprc_252d_z\", \"prc_highprc_252d_c\", \"prc_highprc_252d_d\"]].round(2).T)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "a5e85246",
   "metadata": {},
   "source": [
    "Fifteen columns per stock-month: 5 market-relative characteristics, 5 peer contexts, 5 deviations. NVDA in October 2022: 12-1 momentum \u22121.49 against \u22121.38 for its peers, price / 52-week high \u22122.67 against \u22122.59 \u2014 deviations of \u22120.10 and \u22120.08. With a self-weight of 0.26 and AMD at 0.35, NVDA largely *is* its peer group that month. Only month-$t$ rows enter the transformation. This prevents future cross-sections entering this step, but does not establish historical source availability."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73bc3f3a",
   "metadata": {},
   "source": [
    "### The walk-forward protocol\n",
    "\n",
    "**Predict:** The row dated 2018-12-31 has <code>ret_next</code> = January 2019's return. May it sit in the training set when we predict the 2019-01-31 row?\n",
    "\n",
    "- No \u2014 its target overlaps the first test month\n",
    "- Yes \u2014 its target is realised on 31 Jan 2019, exactly when the 2019-01-31 features are observed\n",
    "- Only if we drop one month between train and test\n",
    "- Yes \u2014 targets can never leak\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Yes \u2014 its target is realised on 31 Jan 2019, exactly when the 2019-01-31 features are observed. The 2018-12-31 row's target is the return over January 2019, realised on 31 January 2019 \u2014 the same instant at which the 2019-01-31 features are observed. So it is known when we predict the 2019-01-31 row, and training on all rows dated strictly before t is legal. A two-month-ahead target would need a one-month gap.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "id": "c4fb060a",
   "metadata": {},
   "source": [
    "from scipy.stats import spearmanr\n",
    "from sklearn.ensemble import HistGradientBoostingRegressor\n",
    "months = np.sort(xs[\"date\"].unique()); test_months = months[months >= np.datetime64(\"2019-01-31\")]\n",
    "def walk_forward(cols):\n",
    "    preds = []\n",
    "    for j in range(0, len(test_months), 12):             # refit every 12 months, predict the next 12\n",
    "        block = test_months[j:j + 12]\n",
    "        trn, tst = xs[xs[\"date\"] < block[0]], xs[xs[\"date\"].isin(block)]\n",
    "        hgb = HistGradientBoostingRegressor(max_depth=2, learning_rate=0.05, max_iter=200, min_samples_leaf=50, random_state=0)\n",
    "        preds.append(pd.Series(hgb.fit(trn[cols], trn[\"ret_next\"]).predict(tst[cols]), index=tst.index))\n",
    "    return pd.concat(preds)\n",
    "def evaluate(pred):\n",
    "    df = xs.loc[pred.index, [\"date\", \"ret_next\"]].assign(pred=pred)\n",
    "    ic = df.groupby(\"date\")[[\"pred\", \"ret_next\"]].apply(lambda g: spearmanr(g[\"pred\"], g[\"ret_next\"])[0])\n",
    "    sp = df.groupby(\"date\")[[\"pred\", \"ret_next\"]].apply(lambda g: g.nlargest(20, \"pred\")[\"ret_next\"].mean() - g.nsmallest(20, \"pred\")[\"ret_next\"].mean())\n",
    "    return ic, sp\n",
    "print(len(test_months), \"test months;\", int(np.ceil(len(test_months) / 12)), \"refits; first training set:\", int((xs[\"date\"] < test_months[0]).sum()), \"rows\")"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "82c153fd",
   "metadata": {},
   "source": [
    "Expanding window from 2019-01: 71 test months, 6 refits (every 12 months, predicting the next 12 \u2014 a monthly refit would be 71 fits, same protocol, ten times the runtime). Two scores per month: **rank IC** (Spearman between prediction and realised return across the 213 stocks) and the **top-20 minus bottom-20** equal-weight spread."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72f60b5a",
   "metadata": {},
   "source": [
    "### Three feature sets, one protocol\n",
    "\n",
    "**Predict:** Which feature set will have the highest mean rank IC \u2014 and will the difference be significant?\n",
    "\n",
    "- Raw standardised $X$: fewer features, less over-fitting\n",
    "- $X$ + sector-relative: sector-relative controls win\n",
    "- $X$ + attention context and deviations: data-driven peers win\n",
    "- No approximate single-set test rejects zero here; pairwise superiority requires a separate paired comparison\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "No approximate single-set test rejects zero here; pairwise superiority requires a separate paired comparison. Five price-based characteristics on 213 large caps over 71 months is a small, noisy problem. Mean IC is within one to two standard errors of zero for all three sets, and the ranking between them flips with the temperature. The honest reading is 'no detectable difference'.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "id": "06cad829",
   "metadata": {},
   "source": [
    "sets = {\"raw z\": zc, \"z + sector\": zc + sc, \"z + attention\": zc + cc + dc}\n",
    "ics, spreads = {}, {}\n",
    "for name, cols in sets.items():\n",
    "    ics[name], spreads[name] = evaluate(walk_forward(cols))\n",
    "    ic, sp = ics[name], spreads[name]\n",
    "    print(f\"{name:14s} {len(cols):2d} features  mean IC {ic.mean():+.4f}  t-stat {ic.mean() / ic.std() * np.sqrt(len(ic)):+.2f}\"\n",
    "          f\"   top20 \u2212 bottom20 {sp.mean() * 12 * 100:+.1f} %/yr  Sharpe {sp.mean() / sp.std() * np.sqrt(12):+.2f}\")"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "fe79986f",
   "metadata": {},
   "source": [
    "Mean rank IC: **\u22120.015** (t = \u22121.04) raw, **\u22120.014** (t = \u22120.99) with sector features, **\u22120.010** (t = \u22120.67) with attention. None of these approximate unadjusted tests rejects zero at the usual 5% level. Failure to reject does not establish zero effect or equivalence. The long-short spread is positive (+7.7 %, +9.3 %, +2.8 % a year; Sharpe 0.47, 0.62, 0.17, before costs) \u2014 the extremes rank slightly better than the middle, but 71 months cannot separate a Sharpe of 0.6 from 0. Attention features did **not** help here; the protocol that tells you so is the deliverable."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26b6d7a7",
   "metadata": {},
   "source": [
    "### Is a Sharpe of 0.5 evidence?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "bcd9820d",
   "metadata": {},
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(9, 4))\n",
    "for name, sp in spreads.items():\n",
    "    plt.plot(sp.index, sp.cumsum() * 100, label=f\"{name}  (t = {sp.mean() / sp.std() * np.sqrt(len(sp)):+.2f})\")\n",
    "plt.axhline(0, color=\"k\", lw=0.5); plt.ylabel(\"sum of monthly spreads (%)\"); plt.legend(); plt.show()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "aeb97719",
   "metadata": {},
   "source": [
    "Under an i.i.d. standard-error approximation, the t-statistic of a spread with Sharpe $S$ over $T$ months is $S\\sqrt{T/12}$: with $T = 71$, the best set's Sharpe of 0.62 gives t = **1.51**; the raw set's 0.47 gives 1.15; attention's 0.17 gives 0.40. To reach t = 2 the sector set would need 124 months of the same performance (10 years), the raw set 18 years. Every gain in the plot comes with an equally plausible zero \u2014 and the file is survivorship-biased in the spread's favour."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ac334c2",
   "metadata": {},
   "source": [
    "### Your turn: change the peer group\n",
    "\n",
    "Change the temperature (`tau_new = 3.0`) **or** build the keys from a subset of characteristics (`keys_new = [\"ret_12_1_z\", \"rvol_21d_z\"]`) so that `ic_new` differs from the baseline attention IC. Does the peer-group definition move the IC by more than its standard error (about 0.015)?"
   ]
  },
  {
   "cell_type": "code",
   "id": "8c5d9ccc",
   "metadata": {},
   "source": [
    "tau_new, keys_new = 1.0, [\"ret_12_1_z\", \"rvol_21d_z\"]  # worked solution\n",
    "xs[cc + dc] = attention_features(tau_new, keys_new)\n",
    "ic_new, sp_new = evaluate(walk_forward(zc + cc + dc))\n",
    "print(\"mean IC:\", round(ic_new.mean(), 4), \" baseline:\", round(ics[\"z + attention\"].mean(), 4),\n",
    "      \" Sharpe:\", round(sp_new.mean() / sp_new.std() * np.sqrt(12), 2))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "d158d7a1",
   "metadata": {},
   "source": [
    "$\\tau = 3$: IC \u22120.010 (unchanged); keys from momentum and volatility only: \u22120.025 (t = \u22121.64); $\\tau = 0.3$: \u22120.041 (t = \u22122.60); $\\tau = 10$: \u22120.008. The sign of the *change* flips with the setting, and the one \"significant\" number appears after four tries \u2014 that is multiple testing, not a signal. In KKMX the peer group is *learned* on a Sharpe objective with thousands of stocks and six decades; here it is a fixed kernel on 213 survivors."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "658e54ea",
   "metadata": {},
   "source": [
    "### What you discovered\n",
    "\n",
    "- The cross-sectional frame: one matrix $X_t$ per month, standardised **per month** (median / std, within sector in production) so that no row uses a number from its own future.\n",
    "\n",
    "- One attention layer is twelve lines of numpy; with $W = I$ its weights are a similarity kernel, and the temperature interpolates between \"nearest stock\" ($\\tau \\to 0$) and \"cross-sectional mean\" ($\\tau \\to \\infty$).\n",
    "\n",
    "- NVDA's data-driven peers in October 2022 were AMD, Netflix, Alibaba and Tesla \u2014 three sectors; a sector median would have seen none of that.\n",
    "\n",
    "- KKMX learn $W$ end-to-end on a Sharpe objective; we froze it and used $C$ and $X - C$ as features for the \u00a74.3 boosting machine.\n",
    "\n",
    "- Walk-forward, 71 months, 6 refits: mean rank IC \u22120.015 / \u22120.014 / \u22120.010 for raw, sector and attention features \u2014 none distinguishable from zero. Spread Sharpe 0.17\u20130.62 needs a decade or more to become evidence.\n",
    "\n",
    "- Next: the Mistakes Library \u2014 what happens when a good predictor becomes an action."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49fa563b",
   "metadata": {},
   "source": [
    "### From Predictions to Risk Decisions: Walk-Forward Models and Portfolio Accounts\n",
    "\n",
    "Bring return forecasts and event probabilities into a declared, testable risk policy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cea0a370",
   "metadata": {},
   "source": [
    "### A forecast becomes a portfolio only through an explicit rule\n",
    "\n",
    "**Diagram concepts:** Return forecast; Rank stocks and select top three; Loss probability; Apply frozen risk threshold; Target weights, retained cash; Trades from drifted holdings; Costs, realised security P&L and NAV.\n",
    "\n",
    "The return model ranks candidates. The classifier can reduce exposure under a declared policy. Evaluate each statistical model separately, then evaluate the complete net account. Weak classification evidence remains a limitation even if one portfolio path looks good."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bd67029",
   "metadata": {},
   "source": [
    "### Rebuild an availability-aware walk-forward sample\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "c26093bf",
   "metadata": {},
   "source": [
    "import numpy as np, pandas as pd, matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import Ridge, LogisticRegression\n",
    "from sklearn.metrics import mean_squared_error, roc_auc_score, brier_score_loss\n",
    "case = pd.read_csv(\"https://statpython.pages.dev/data/chapter2_scored.csv\",\n",
    "                   parse_dates=[\"date\", \"label_available_at\"])\n",
    "case = case.sort_values([\"date\", \"ticker\"]).reset_index(drop=True)\n",
    "model_keys = [\"date\", \"ticker\"]\n",
    "assert not case.duplicated(model_keys).any()\n",
    "numeric_features = [\"ret_12_1\", \"rvol_21d\", \"momentum_sector_z\"]\n",
    "sector_dummies = pd.get_dummies(case.sector, prefix=\"sector\", drop_first=True, dtype=float)\n",
    "X_case = pd.concat([case[numeric_features], sector_dummies], axis=1)\n",
    "assert np.isfinite(X_case.to_numpy()).all()\n",
    "y_case = case.ret_next.where(case.label_valid)\n",
    "loss_event = y_case.le(-0.05).where(case.label_valid)\n",
    "train_case = case.date.lt(\"2021-01-01\") & case.label_available_at.lt(\"2021-01-01\") & case.label_valid\n",
    "valid_case = (case.date.ge(\"2021-01-01\") & case.date.lt(\"2022-01-01\")\n",
    "              & case.label_available_at.lt(\"2022-01-01\") & case.label_valid)\n",
    "test_case = case.date.ge(\"2022-01-01\") & case.label_valid\n",
    "print(\"train / validation / test / pending:\",\n",
    "      int(train_case.sum()), int(valid_case.sum()), int(test_case.sum()), int((~case.label_valid).sum()))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "ee25c264",
   "metadata": {},
   "source": [
    "For each calendar year, refit using labels available **strictly before 1 January**. Predict every month in that year using its available features. This conservative annual protocol differs from refitting at each month-end.\n",
    "\n",
    "Ridge alpha 100 (selected by the fixed validation grid), logit C 1, the 5% event definition and risk threshold 0.65 are frozen settings. The threshold grid was studied using 2021 validation. Do not retune them using the subsequent portfolio outcomes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3b8a20a",
   "metadata": {},
   "source": [
    "### A function isolates the stateful refitting logic\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "815abfba",
   "metadata": {},
   "source": [
    "def annual_predictions(frame, matrix, alpha=100.0):\n",
    "    records = []\n",
    "    for year in [2022, 2023, 2024]:\n",
    "        freeze = pd.Timestamp(year=year, month=1, day=1)\n",
    "        train = frame.date.lt(freeze) & frame.label_available_at.lt(freeze) & frame.label_valid\n",
    "        predict = frame.date.dt.year.eq(year)\n",
    "        return_model = make_pipeline(StandardScaler(), Ridge(alpha=alpha))\n",
    "        event_model = make_pipeline(StandardScaler(), LogisticRegression(C=1, max_iter=1000))\n",
    "        return_model.fit(matrix.loc[train], frame.loc[train, \"ret_next\"])\n",
    "        event_model.fit(matrix.loc[train], frame.loc[train, \"ret_next\"].le(-0.05).astype(int))\n",
    "        part = frame.loc[predict, model_keys + [\"ret_next\", \"label_valid\"]].copy()\n",
    "        part[\"return_forecast\"] = return_model.predict(matrix.loc[predict])\n",
    "        part[\"loss_probability\"] = event_model.predict_proba(matrix.loc[predict])[:, 1]\n",
    "        part[\"trained_through\"] = frame.loc[train, \"label_available_at\"].max()\n",
    "        assert part.trained_through.lt(freeze).all()\n",
    "        records.append(part)\n",
    "    return pd.concat(records).sort_values(model_keys)\n",
    "walk = annual_predictions(case, X_case, alpha=100.0)\n",
    "assert not walk.duplicated(model_keys).any()\n",
    "print(\"predictions / evaluated / pending:\", len(walk), int(walk.label_valid.sum()), int((~walk.label_valid).sum()))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "5587fbdd",
   "metadata": {},
   "source": [
    "There are **315 predictions, 306 evaluable outcomes and 9 pending predictions**. Preprocessing is refitted within each training sample. Prediction keys survive concatenation.\n",
    "\n",
    "The loop represents three successive fitting clocks. It does not mix future outcomes into an earlier fit. Exact execution, historical membership and vendor timestamps remain unverified classroom assumptions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cced1cb6",
   "metadata": {},
   "source": [
    "### Risk filtering retains cash instead of silently releveraging\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "bafc7e5e",
   "metadata": {},
   "source": [
    "walk[\"rank\"] = walk.groupby(\"date\").return_forecast.rank(ascending=False, method=\"first\")\n",
    "walk[\"base_weight\"] = np.where(walk[\"rank\"].le(3), 0.99 / 3, 0.0)\n",
    "walk[\"risk_weight\"] = walk.base_weight.where(walk.loss_probability.lt(0.65), 0.0)\n",
    "exposure = walk.groupby(\"date\")[[\"base_weight\", \"risk_weight\"]].sum()\n",
    "assert exposure.le(0.99 + 1e-12).all().all()\n",
    "print(\"mean invested share:\\n\", exposure.mean().round(3))\n",
    "print(\"mean retained cash share:\\n\", (1 - exposure).mean().round(3))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "05b69053",
   "metadata": {},
   "source": [
    "At the frozen 0.65 threshold, this sample blocks none of the selected positions, so the baseline and filtered accounts are identical. This rule supplies no realised risk reduction here. An allowed selected stock keeps weight 0.33. A blocked stock receives zero, and its allocation stays in cash. Rescaling the survivors to 99% would define a different risk policy.\n",
    "\n",
    "The row order provides a deterministic tie rule in this exercise. A production specification must state tie-breaking, capacity, borrow and concentration constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe807333",
   "metadata": {},
   "source": [
    "### Accounting follows holdings, cash and fees through time\n",
    "\n",
    "This function carries holdings and cash across periods, so a short stateful loop is appropriate. Trades start from **drifted holdings**, rather than last month\u2019s target weights. Cash earns zero in this exercise."
   ]
  },
  {
   "cell_type": "code",
   "id": "13393f74",
   "metadata": {},
   "source": [
    "def reconcile_account(weights, returns, cost_bps=10):\n",
    "    holdings = pd.Series(0.0, index=weights.columns)\n",
    "    cash = 100.0\n",
    "    rows = []\n",
    "    for date in weights.index:\n",
    "        opening_nav = float(cash + holdings.sum())\n",
    "        target = weights.loc[date] * opening_nav\n",
    "        traded = float((target - holdings).abs().sum())\n",
    "        fee = traded * cost_bps / 10000\n",
    "        cash = opening_nav - float(target.sum()) - fee\n",
    "        assert cash >= -1e-10\n",
    "        security_pnl = target * returns.loc[date]\n",
    "        holdings = target + security_pnl\n",
    "        closing_nav = float(cash + holdings.sum())\n",
    "        rows.append({\"date\": date, \"nav_start\": opening_nav,\n",
    "                     \"gross_pnl\": float(security_pnl.sum()), \"cost\": fee,\n",
    "                     \"cash\": cash, \"nav_end\": closing_nav, \"traded_notional\": traded})\n",
    "    ledger = pd.DataFrame(rows).set_index(\"date\")\n",
    "    assert np.allclose(ledger.nav_end - ledger.nav_start, ledger.gross_pnl - ledger.cost)\n",
    "    ledger[\"net_return\"] = ledger.nav_end / ledger.nav_start - 1\n",
    "    return ledger\n",
    "realised_walk = walk.loc[walk.label_valid].copy()\n",
    "return_matrix = realised_walk.pivot(index=\"date\", columns=\"ticker\", values=\"ret_next\")\n",
    "base_weights = realised_walk.pivot(index=\"date\", columns=\"ticker\", values=\"base_weight\")\n",
    "risk_weights = realised_walk.pivot(index=\"date\", columns=\"ticker\", values=\"risk_weight\")\n",
    "assert return_matrix.notna().all().all()\n",
    "base_account = reconcile_account(base_weights, return_matrix)\n",
    "risk_account = reconcile_account(risk_weights, return_matrix)\n",
    "print(\"realised monthly accounts:\", len(base_account), len(risk_account))"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "23111fd2",
   "metadata": {},
   "source": [
    "Both accounts contain **34 realised monthly holding periods**. Initial purchases incur fees. NAV is measured after fees and realised simple-return P&L. The identity `NAV change = gross P&L \u2212 fees` must hold at every row."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14aacf2a",
   "metadata": {},
   "source": [
    "### Compare net outcomes with their exposure and uncertainty\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "b8e0add2",
   "metadata": {},
   "source": [
    "def account_summary(account):\n",
    "    wealth = np.r_[100.0, account.nav_end.to_numpy()]\n",
    "    drawdown = wealth / np.maximum.accumulate(wealth) - 1\n",
    "    return pd.Series({\"ending_nav\": account.nav_end.iloc[-1],\n",
    "        \"mean_monthly_net_return\": account.net_return.mean(),\n",
    "        \"max_drawdown\": drawdown.min(), \"total_fees\": account.cost.sum()})\n",
    "account_comparison = pd.DataFrame({\"baseline\": account_summary(base_account),\n",
    "                                   \"risk_filter\": account_summary(risk_account)})\n",
    "print(account_comparison.round(4))\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.plot(base_account.index + pd.offsets.MonthEnd(1), base_account.nav_end, label=\"baseline\")\n",
    "ax.plot(risk_account.index + pd.offsets.MonthEnd(1), risk_account.nav_end, label=\"risk filter\")\n",
    "ax.set(xlabel=\"Holding-period outcome date\", ylabel=\"Net NAV from opening 100\")\n",
    "ax.legend(); plt.show()"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "2132bd94",
   "metadata": {},
   "source": [
    "Report the actual values produced, exposure differences, turnover and fee assumptions. Apply Chapter 3\u2019s paired, dependence-aware uncertainty to the monthly return difference.\n",
    "\n",
    "Thirty-four months and nine surviving stocks cannot establish a reliable risk reduction. A net account can reconcile perfectly while its model generalises poorly."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b425ae0",
   "metadata": {},
   "source": [
    "### Self-check: cash is part of the policy\n",
    "\n",
    "**Predict:** One of three 0.33-weight stocks is blocked by the risk rule. What follows from this policy?\n",
    "\n",
    "- The other stocks automatically rise to 0.495 each\n",
    "- The missing allocation stays in cash, and exposure falls from 0.99 to 0.66 before fees\n",
    "- The blocked stock\u2019s return becomes zero in the source file\n",
    "- No transaction can occur\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "The missing allocation stays in cash, and exposure falls from 0.99 to 0.66 before fees. The declared filter keeps the original weights on allowed names and sends blocked allocations to cash. Renormalising creates a different policy.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e47bbe0",
   "metadata": {},
   "source": [
    "### The decision memo links statistical evidence to action\n",
    "\n",
    "Report **target and horizon; candidate/label counts; feature availability; training/validation/test clocks; fitted models; baseline metrics; frozen decision costs and threshold; exposure, net fees and drawdown; uncertainty and limitations**.\n",
    "\n",
    "The fixed classifier\u2019s test AUC of 0.431 does not support claiming a dependable loss detector. The attention extension does not establish that a transformer learns a profitable strategy. Keep negative findings visible.\n",
    "\n",
    "**Your investigation:** Change one declared assumption on a development sample, record the attempt, explain the changed predictions and account, and propose fresh evidence for evaluation. Do not rename an examined test \u201cuntouched\u201d."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b49cc4a2",
   "metadata": {},
   "source": [
    "### Pending predictions stay outside realised accounting\n",
    "\n",
    "**Predict:** The last nine candidates have pending outcomes. What is appropriate?\n",
    "\n",
    "- Delete them before predicting\n",
    "- Predict them if features are ready, but omit their unobserved holding periods from realised evaluation\n",
    "- Fill future returns with the training mean\n",
    "- Use their label-status flags as return predictors\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "Predict them if features are ready, but omit their unobserved holding periods from realised evaluation. Prediction eligibility and outcome availability are separate. Realised accounting needs actual mature outcomes.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec90f2d",
   "metadata": {},
   "source": [
    "### Reconciliation is necessary evidence\n",
    "\n",
    "**Predict:** If NAV change equals P&L minus fees at every date, what has been established?\n",
    "\n",
    "- The strategy is profitable out of sample\n",
    "- The specified accounting identity holds for the supplied inputs\n",
    "- All source records were known historically\n",
    "- The model has a causal effect\n",
    "\n",
    "<details><summary>Feedback</summary>\n",
    "\n",
    "The specified accounting identity holds for the supplied inputs. Reconciliation checks accounting consistency. Source availability, generalisation, causality and execution require separate evidence.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5581d06d",
   "metadata": {},
   "source": [
    "### Chapter Summary\n",
    "\n",
    "| Concept | What to explain |\n",
    "|---|---|\n",
    "| Regression | Conditional mean, residuals, assumptions and forecast errors |\n",
    "| Inference | Covariance estimator, confidence interval, tests and selection uncertainty |\n",
    "| Classification | Event horizon, odds, calibration, ranking and threshold costs |\n",
    "| Causal reasoning | Identification assumptions rather than predictive association |\n",
    "| Feature context | Fixed attention as engineering, with availability and survival limits |\n",
    "| Risk decisions | Exposure, pending labels, drifted holdings, fees and NAV reconciliation |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ea847c2",
   "metadata": {},
   "source": [
    "### Discussion Questions\n",
    "\n",
    "1. Why does a low training RMSE not establish a useful future model?\n",
    "2. What does the fixed classifier\u2019s weak test result imply for a risk decision?\n",
    "3. Which assumptions support a causal effect rather than a prediction?\n",
    "4. Why can a correctly reconciled account still provide weak evidence for a strategy?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
