{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ada17624",
   "metadata": {},
   "source": [
    "# Chapter 5 — Bayesian Investment Research\n",
    "\n",
    "This tutorial mirrors Book and Slides §§5.5–5.6. Run cells in order. Predict an output, run it, change one declared assumption and explain the difference. The two repair exercises contain labelled worked solutions; try the online starters first.\n",
    "\n",
    "**Prerequisites:** Chapter 2 keyed features, Chapter 3 uncertainty, Chapter 4 chronological validation, and Chapter 5 priors, likelihoods and conjugate regression. Packages: numpy, pandas, scipy and matplotlib. The first data cell downloads the public Chapter 2 CSV; no private MyStrategy data or credentials are used.\n",
    "\n",
    "**Scope:** a nine-stock survivor sample with static sectors and assumed availability, already examined in earlier chapters. Fixed residual and hierarchy scales are instructional assumptions. This historical teaching experiment is not a fresh validation of a trading strategy.\n",
    "\n",
    "**Learning goals:** compute partial pooling, fit a hierarchical normal signal, distinguish mean and predictive uncertainty, compare prior models, generate three annual forward blocks, retain pending forecasts and reconcile cash, drift and fees against simple portfolios.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e28bcca4",
   "metadata": {},
   "source": [
    "## Bayesian Signals and Partial Pooling Across Stocks\n",
    "\n",
    "Use the keyed features built in Chapter 2 to predict next-month returns. Learn what a hierarchical prior shares across stocks, how it changes noisy effects, and which probability is relevant to a decision.\n",
    "\n",
    "### An investment question needs a target and a clock\n",
    "\n",
    "The question is: **which candidates have favourable expected next-month returns given their available characteristics?**\n",
    "\n",
    "Use three features: momentum, realised volatility and sector-relative momentum. Preserve stock/date keys. Keep future returns outside the feature matrix and retain pending candidates.\n",
    "\n",
    "The nine-stock classroom panel has survivor membership, static sectors and assumed monthly availability. Its previous predictions are ignored. Its realised results have already been explored in earlier chapters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2735fe24",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, pandas as pd, matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from scipy.special import logsumexp\n",
    "investment = pd.read_csv(\"https://statpython-book.pages.dev/data/chapter2_scored.csv\",\n",
    "                         parse_dates=[\"date\", \"label_available_at\"])\n",
    "investment = investment.sort_values([\"date\", \"ticker\"]).reset_index(drop=True)\n",
    "investment_features = [\"ret_12_1\", \"rvol_21d\", \"momentum_sector_z\"]\n",
    "investment_tickers = [\"AAPL\", \"BAC\", \"COP\", \"CVX\", \"GS\", \"JPM\", \"MSFT\", \"NVDA\", \"XOM\"]\n",
    "investment_sectors = [\"Financial Services\", \"Technology\"]\n",
    "assert not investment.duplicated([\"date\", \"ticker\"]).any()\n",
    "assert np.isfinite(investment[investment_features].to_numpy()).all()\n",
    "print(\"candidates / mature / pending:\", len(investment), int(investment.label_valid.sum()),\n",
    "      int((~investment.label_valid).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51e38703",
   "metadata": {},
   "source": [
    "A target is an outcome over a declared horizon. Here it is next-month simple return, rather than price level, a cumulative multi-month return or a causal treatment effect. A stock can be eligible for prediction while its outcome remains unavailable. That distinction carries through fitting, scoring and accounting.\n",
    "\n",
    "The public panel has 639 candidate rows, 630 valid mature labels and nine pending labels. Some raw pending return values happen to be present in the file; the declared validity flags still exclude them from fitting and evaluation. Treating those raw numbers as realised outcomes would violate the exercise contract. The earlier stored prediction column is also excluded: this chapter constructs new predictions from its own training samples.\n",
    "\n",
    "### Partial pooling sits between two strong assumptions\n",
    "\n",
    "No pooling estimates each stock separately. Complete pooling assumes a common effect. Partial pooling lets stocks differ while sharing a distribution:\n",
    "\n",
    "$$a_i\\sim N(0,\\tau^2),\\qquad \\bar e_i\\mid a_i\\sim N(a_i,\\sigma^2/n_i).$$\n",
    "\n",
    "$$E[a_i\\mid\\bar e_i]=w_i\\bar e_i,\\qquad w_i=\\frac{\\tau^2}{\\tau^2+\\sigma^2/n_i}.$$\n",
    "\n",
    "Here $\\bar e_i$ is a residual mean after a declared baseline. A large noisy residual mean need not imply a large persistent stock effect.\n",
    "\n",
    "A stock-specific effect represents the part of the conditional mean not described by the supplied characteristics and sector effects. It is not automatically a risk-adjusted alpha: that interpretation requires an appropriate factor benchmark and return definition. The exercise calls it a stock deviation to avoid borrowing an economic interpretation that the model has not earned.\n",
    "\n",
    "The prior standard deviation $\\tau$ describes heterogeneity across stocks. A small value expresses a belief that persistent differences are small; a large value permits more dispersion. The sampling variance $\\sigma^2/n_i$ describes how noisy each individual estimate is. Their ratio determines the weight on that stock's own evidence. This is the same precision-weighting used in the conjugate normal update, now applied across assets.\n",
    "\n",
    "### Predict which apparent winner shrinks most\n",
    "\n",
    "\n",
    "**Self-check:** Stock A has only three monthly residuals averaging +8%. Stock C has sixty averaging −1%. Under the same prior and residual variance, which receives less weight on its own sample mean?\n",
    "\n",
    "- A, because three observations carry less precision\n",
    "- C, because negative estimates must be removed\n",
    "- Both receive the same weight\n",
    "- Neither: Bayesian estimates must equal sample means\n",
    "\n",
    "**Feedback:** The weight depends on sampling precision and prior variance. Fewer observations imply greater pooling; the sign alone does not determine the weight.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Compute the shrinkage instead of declaring a winner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37b787a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "toy = pd.DataFrame({\"stock\": [\"A\", \"B\", \"C\"], \"n\": [3, 12, 60],\n",
    "                    \"observed_residual_mean\": [0.08, 0.035, -0.01]})\n",
    "tau, residual_sd = 0.02, 0.10\n",
    "toy[\"sampling_var\"] = residual_sd**2 / toy.n\n",
    "toy[\"own_weight\"] = tau**2 / (tau**2 + toy.sampling_var)\n",
    "toy[\"posterior_mean\"] = toy.own_weight * toy.observed_residual_mean\n",
    "toy[\"posterior_sd\"] = np.sqrt(1 / (1/tau**2 + 1/toy.sampling_var))\n",
    "print(toy[[\"stock\", \"own_weight\", \"posterior_mean\", \"posterior_sd\"]].round(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "270593c2",
   "metadata": {},
   "source": [
    "Stock A shrinks from **8% to 0.86%**, with own-data weight 0.1071. Stock B becomes 1.14%; C becomes −0.71%. The largest raw sample mean is no longer the largest posterior mean.\n",
    "\n",
    "The fixture is simulated to vary sample sizes transparently. It demonstrates a mechanism, rather than proving that the pooled ranking will earn money.\n",
    "\n",
    "\n",
    "\n",
    "### Change the prior and explain what moves\n",
    "\n",
    "*Use the online Book’s prior slider to explore the same shrinkage calculation interactively.*\n",
    "**Try:** compare τ = 0.5%, 2% and 10%. Explain both the posterior mean and interval. Prior sensitivity is part of the result.\n",
    "\n",
    "\n",
    "\n",
    "### A feature model and stock effects can be estimated together\n",
    "\n",
    "$$r_{i,t+1}=b_0+z_{it}^{\\top}\\beta+s_i^{\\top}\\gamma+a_i+\\varepsilon_{it},\\qquad \\varepsilon_{it}\\sim N(0,\\sigma^2).$$\n",
    "\n",
    "$z$ contains three train-standardised features; $s$ contains two sector indicators with Energy as reference. Nine stock effects share $a_i\\sim N(0,\\tau^2)$.\n",
    "\n",
    "Fixed exercise scales: intercept/features 5%, sector effects 3%, stock effects τ = 2%, residual sd σ = 10%. These are assumptions for sensitivity analysis, rather than universally appropriate market values.\n",
    "\n",
    "This is a hierarchical normal model conditional on fixed hyperparameters. There are two levels: the observed returns depend on stock effects, and those effects share a prior distribution. A fully Bayesian extension would assign hyperpriors to $\\tau$ and $\\sigma$, estimate them and integrate over their uncertainty. The browser example fixes them so that the posterior is available in closed form and every matrix operation can be inspected.\n",
    "\n",
    "The design contains an intercept, sector indicators and all nine stock indicators. The corresponding unpenalised design is redundant, but the proper Gaussian priors give a positive-definite posterior precision matrix. Consequently a conventional OLS decomposition cannot be copied over mechanically. Individual stock deviations are interpreted conditional on this prior and baseline. Predictions combine all terms. New stocks would need an explicit unseen-identity policy, such as integrating a new effect drawn from the shared prior.\n",
    "\n",
    "### Standardisation belongs inside the fitting boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94652626",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_design(frame, centre, spread):\n",
    "    scaled = (frame[investment_features] - centre) / spread\n",
    "    sector = pd.get_dummies(frame.sector).reindex(columns=investment_sectors, fill_value=0)\n",
    "    stock = pd.get_dummies(frame.ticker).reindex(columns=investment_tickers, fill_value=0)\n",
    "    return np.column_stack([np.ones(len(frame)), scaled, sector, stock]).astype(float)\n",
    "\n",
    "investment_names = [\"intercept\"] + investment_features + investment_sectors + investment_tickers\n",
    "print(\"design columns:\", len(investment_names), \"stock deviations:\", len(investment_tickers))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af51ee1f",
   "metadata": {},
   "source": [
    "The design has **15 columns**. The function receives the training mean and standard deviation; it never estimates a new scale from prediction rows.\n",
    "\n",
    "Ticker indicators represent this declared teaching universe. They do not establish historical membership. A production universe and sector history need separate point-in-time evidence.\n",
    "\n",
    "\n",
    "\n",
    "### The normal posterior makes the calculation transparent\n",
    "\n",
    "For $\\theta\\sim N(0,V_0)$ and $y\\mid\\theta\\sim N(X\\theta,\\sigma^2I)$:\n",
    "\n",
    "$$V_n=(V_0^{-1}+X^{\\top}X/\\sigma^2)^{-1},\\qquad m_n=V_nX^{\\top}y/\\sigma^2.$$\n",
    "\n",
    "The same posterior mean minimises squared error plus a quadratic penalty. **A prior changes parameter estimates and predictions**, while the posterior covariance adds uncertainty around them.\n",
    "\n",
    "The next function also calculates the normalised marginal likelihood needed in §5.6.\n",
    "\n",
    "Prior precision and data precision add. The stock-effect block of $V_0$ determines pooling strength; the feature and sector blocks express separate prior scales. Standardising features on training rows makes a one-unit coefficient correspond to one training standard deviation of that characteristic. It does not make economically different prior choices equivalent.\n",
    "\n",
    "The likelihood assumes constant residual variance and independent rows. Stocks observed in the same month often share shocks; monthly returns can also exhibit temporal dependence. The posterior intervals in this exercise therefore describe a simplified model, not dependence-robust guarantees. A richer model can include common factors or structured residual covariance. Checking this assumption is as necessary as choosing a prior.\n",
    "\n",
    "### Fit parameters and preserve their covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1431b9dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_posterior(training, stock_tau=0.02, residual_sd=0.10):\n",
    "    centre = training[investment_features].mean()\n",
    "    spread = training[investment_features].std(ddof=0).replace(0, 1)\n",
    "    matrix = investment_design(training, centre, spread)\n",
    "    outcome = training.ret_next.to_numpy()\n",
    "    prior_var = np.r_[0.05**2, np.repeat(0.05**2, 3), np.repeat(0.03**2, 2),\n",
    "                      np.repeat(stock_tau**2, 9)]\n",
    "    precision = np.diag(1/prior_var) + matrix.T @ matrix / residual_sd**2\n",
    "    covariance = np.linalg.solve(precision, np.eye(len(prior_var)))\n",
    "    mean = np.linalg.solve(precision, matrix.T @ outcome / residual_sd**2)\n",
    "    log_evidence = (-len(outcome)/2 * np.log(2*np.pi*residual_sd**2)\n",
    "        + 0.5 * (np.linalg.slogdet(covariance)[1] - np.log(prior_var).sum())\n",
    "        - 0.5 * (outcome @ outcome / residual_sd**2 - mean @ precision @ mean))\n",
    "    return dict(mean=mean, covariance=covariance, centre=centre, spread=spread,\n",
    "                residual_sd=residual_sd, log_evidence=float(log_evidence))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32fd1fc7",
   "metadata": {},
   "source": [
    "The function uses linear solves rather than explicitly inverting the precision matrix. All prior variances are strictly positive. Its evidence calculation includes normalising constants; an unnormalised posterior density cannot replace it.\n",
    "\n",
    "\n",
    "\n",
    "### Read coefficients in their transformed units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c8a1d44",
   "metadata": {},
   "outputs": [],
   "source": [
    "investment_freeze = pd.Timestamp(\"2021-01-01\")\n",
    "investment_train = investment.loc[investment.date.lt(investment_freeze)\n",
    "    & investment.label_available_at.lt(investment_freeze) & investment.label_valid].copy()\n",
    "investment_model = investment_posterior(investment_train)\n",
    "coefficient_table = pd.DataFrame({\"mean\": investment_model[\"mean\"],\n",
    "    \"sd\": np.sqrt(np.diag(investment_model[\"covariance\"]))}, index=investment_names)\n",
    "print(\"training rows / outcome dates:\", len(investment_train), investment_train.date.nunique())\n",
    "print(coefficient_table.loc[investment_features].round(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9db01f1b",
   "metadata": {},
   "source": [
    "There are **207 training rows over 23 outcome periods** before the 1 January 2021 freeze. Posterior feature means are −0.0129, +0.0331 and +0.0050 per training standard deviation.\n",
    "\n",
    "A coefficient's sign is conditional on the included predictors, sector effects, stock effects and prior. It does not establish a causal trading mechanism.\n",
    "\n",
    "The difference between 207 rows and 23 monthly snapshots matters. Nine stocks in one month do not necessarily supply nine independent market experiments. The code's numerical solution is exact for its Gaussian assumptions; the economic evidence still depends on those assumptions and the source data.\n",
    "\n",
    "An estimate of −0.0129 on standardised momentum means that, within this model, increasing momentum by one training standard deviation changes expected next-month simple return by about −1.29 percentage points, holding the other columns fixed. It is neither a raw-unit slope nor a statement that increasing historical momentum causes returns to fall.\n",
    "\n",
    "### Parameter uncertainty and future-return uncertainty differ\n",
    "\n",
    "For a candidate row $x_*$:\n",
    "\n",
    "$$E[r_*\\mid D]=x_*^\\top m_n,\\quad \\operatorname{Var}(E[r_*\\mid\\theta]\\mid D)=x_*^\\top V_nx_*.$$\n",
    "$$\\operatorname{Var}(r_*\\mid D)=\\sigma^2+x_*^\\top V_nx_*.$$\n",
    "\n",
    "A credible interval for the conditional mean leaves out the new observation's residual shock. A predictive interval includes it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2e5c4f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_prediction(model, candidates):\n",
    "    matrix = investment_design(candidates, model[\"centre\"], model[\"spread\"])\n",
    "    mean = matrix @ model[\"mean\"]\n",
    "    mean_var = np.einsum(\"ij,jk,ik->i\", matrix, model[\"covariance\"], matrix)\n",
    "    mean_var = np.maximum(mean_var, 1e-15)\n",
    "    predictive_var = mean_var + model[\"residual_sd\"]**2\n",
    "    return mean, mean_var, predictive_var\n",
    "\n",
    "investment_snapshot = investment.loc[investment.date.eq(\"2021-01-31\")].copy()\n",
    "conditional_mean, parameter_var, predictive_var = investment_prediction(investment_model, investment_snapshot)\n",
    "print(\"median sd of conditional mean / future return:\",\n",
    "      round(float(np.median(np.sqrt(parameter_var))), 4),\n",
    "      round(float(np.median(np.sqrt(predictive_var))), 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ca3f036",
   "metadata": {},
   "source": [
    "For the January 2021 candidates, the median posterior standard deviation of the conditional mean is 1.83%, while the predictive standard deviation is 10.17%. Most of the next month's uncertainty comes from the residual return, not from uncertainty about the fitted mean. A narrow interval for expected return is compatible with a wide range of realised returns.\n",
    "\n",
    "This is why a probability that the conditional mean is positive cannot be described as the probability that next month's trade will win. Both can be computed, but they answer different questions and use different variances. A portfolio-level prediction also needs cross-stock covariance; adding individual predictive variances while assuming zero covariance can severely understate aggregate risk.\n",
    "\n",
    "### A probability of positive mean is not a win probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "640b9d23",
   "metadata": {},
   "outputs": [],
   "source": [
    "investment_snapshot[\"forecast\"] = conditional_mean\n",
    "investment_snapshot[\"p_mean_above_cost\"] = norm.sf((0.002 - conditional_mean) / np.sqrt(parameter_var))\n",
    "investment_snapshot[\"p_return_above_cost\"] = norm.sf((0.002 - conditional_mean) / np.sqrt(predictive_var))\n",
    "print(investment_snapshot[[\"ticker\", \"forecast\", \"p_mean_above_cost\",\n",
    "                           \"p_return_above_cost\"]].round(3).head(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "197a630f",
   "metadata": {},
   "source": [
    "For AAPL, `P(conditional mean > 0.2%)` is **0.974**, whereas `P(next return > 0.2%)` is **0.638**.\n",
    "\n",
    "The 0.2% screen is an assumed round-trip hurdle for this experiment, rather than a prediction of actual turnover costs. Cash earns zero. Both probabilities are conditional on the declared likelihood and prior.\n",
    "\n",
    "\n",
    "\n",
    "### Self-check: distinguish the two probabilities\n",
    "\n",
    "\n",
    "**Self-check:** A posterior gives P(conditional mean > cost) = 0.97. What follows?\n",
    "\n",
    "- The next trade has a 97% probability of making money\n",
    "- The conditional mean exceeds the hurdle with 97% posterior probability under this model\n",
    "- The stock has a causal alpha of 97%\n",
    "- The strategy will be profitable after every rebalance\n",
    "\n",
    "**Feedback:** Realised return includes a new residual shock. Its posterior predictive probability generally differs from the probability about the conditional mean.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Check whether a stronger prior changes the decision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33f6e672",
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_sensitivity = []\n",
    "for stock_tau in [0.005, 0.02, 0.08]:\n",
    "    model = investment_posterior(investment_train, stock_tau)\n",
    "    means, variances, _ = investment_prediction(model, investment_snapshot)\n",
    "    prior_sensitivity.append({\"stock_tau\": stock_tau,\n",
    "        \"max_stock_deviation\": np.max(np.abs(model[\"mean\"][-9:])),\n",
    "        \"median_mean_sd\": np.median(np.sqrt(variances)),\n",
    "        \"top_stock\": investment_snapshot.iloc[int(np.argmax(means))].ticker})\n",
    "print(pd.DataFrame(prior_sensitivity).round(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe2248b5",
   "metadata": {},
   "source": [
    "Across τ = 0.5%, 2% and 8%, the maximum absolute stock deviation changes from **0.11% to 1.26% to 4.35%**. AAPL remains the highest forecast in this snapshot.\n",
    "\n",
    "An unchanged top name does not imply that predictions, confidence or portfolio weights are insensitive. Report what changed and what stayed stable.\n",
    "\n",
    "\n",
    "\n",
    "### Repair exercise: uncertainty needs the observation noise\n",
    "\n",
    "Predictive variance must include both parameter uncertainty and residual variation. Fix the assignment; keep the conditional-mean variance separately.\n",
    "**Repair task:** Repair the omitted uncertainty term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e7a0074",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Worked repair: observation noise belongs in future-return variance.\n",
    "reported_predictive_var = parameter_var + investment_model[\"residual_sd\"]**2\n",
    "assert np.allclose(reported_predictive_var, predictive_var)\n",
    "print(\"reported median predictive sd:\", round(float(np.median(np.sqrt(reported_predictive_var))), 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69db3fa7",
   "metadata": {},
   "source": [
    "**Worked repair:** `reported_predictive_var = parameter_var + investment_model[\"residual_sd\"]**2`. The online starter reports the wrong information set; the repair preserves the separate quantities.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### What you discovered\n",
    "\n",
    "- Partial pooling changes a noisy signal rather than merely attaching an error bar.\n",
    "- A hierarchical prior shares information across stocks; fixed hyperparameters keep the exercise transparent.\n",
    "- Posterior covariance describes parameter uncertainty. Predictive variance also includes future residual noise.\n",
    "- Prior sensitivity, dependence and timing remain part of the research result.\n",
    "\n",
    "**Next:** §5.6 combines competing models and takes keyed forecasts through a declared portfolio policy.\n",
    "\n",
    "\n",
    "\n",
    "## Model Uncertainty, Portfolio Decisions and Walk-Forward Evidence\n",
    "\n",
    "Bayesian modelling can produce trading signals and portfolio decisions. This section compares prior models, builds predictive mixtures and evaluates complete net accounts. Walk-forward supplies the time boundary; the Bayesian model supplies the forecast distribution.\n",
    "\n",
    "### One research chain connects the chapters\n",
    "\n",
    "**Research flow:** keyed features → eligible historical labels → posterior models → predictive mixture → fixed portfolio rule → realised net accounts.\n",
    "\n",
    "The process must be auditable at each fitting cutoff. Posterior calculations do not repair future information in an input, an overlapping label or an incorrectly reconstructed trading account.\n",
    "\n",
    "\n",
    "\n",
    "### Model probabilities use integrated evidence\n",
    "\n",
    "For a declared candidate set $M_1,\\ldots,M_K$:\n",
    "\n",
    "$$p(M_k\\mid D)=\\frac{p(D\\mid M_k)p(M_k)}{\\sum_jp(D\\mid M_j)p(M_j)},\\qquad p(D\\mid M_k)=\\int p(D\\mid\\theta_k,M_k)p(\\theta_k\\mid M_k)d\\theta_k.$$\n",
    "\n",
    "The evidence averages fit over the prior; it is not the likelihood at the best-fitting parameter. Proper priors and their normalisation matter.\n",
    "\n",
    "Our three models differ only in stock-effect prior sd τ = 0.5%, 2%, 8%, with equal prior model probabilities.\n",
    "\n",
    "A maximised likelihood rewards fit without integrating parameter uncertainty. Marginal likelihood instead accounts for the probability of the observations across the model's prior parameter distribution. Priors that spread mass over many implausible parameter values can have low evidence even if a small region fits well. Consequently an arbitrary change in prior width can change model weights. Bayes factors are ratios of evidences; they inherit this sensitivity.\n",
    "\n",
    "All three models here condition on the same training observations, feature construction and fixed residual standard deviation. Their evidence weights are probabilities within this finite declared model set. They are not probabilities that any model is universally true, or probabilities of a profitable backtest. The candidate set omits heavy tails and residual correlation, so a dominant weight does not establish that those omissions are harmless.\n",
    "\n",
    "### Compute an evidence-weighted predictive mixture\n",
    "\n",
    "This section is standalone: run the setup functions in order. Models are recomputed from the current training sample; future portfolio performance never determines their weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64c4dc9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, pandas as pd, matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from scipy.special import logsumexp\n",
    "investment = pd.read_csv(\"https://statpython-book.pages.dev/data/chapter2_scored.csv\",\n",
    "                         parse_dates=[\"date\", \"label_available_at\"])\n",
    "investment = investment.sort_values([\"date\", \"ticker\"]).reset_index(drop=True)\n",
    "investment_features = [\"ret_12_1\", \"rvol_21d\", \"momentum_sector_z\"]\n",
    "investment_tickers = [\"AAPL\", \"BAC\", \"COP\", \"CVX\", \"GS\", \"JPM\", \"MSFT\", \"NVDA\", \"XOM\"]\n",
    "investment_sectors = [\"Financial Services\", \"Technology\"]\n",
    "assert not investment.duplicated([\"date\", \"ticker\"]).any()\n",
    "assert np.isfinite(investment[investment_features].to_numpy()).all()\n",
    "print(\"candidates / mature / pending:\", len(investment), int(investment.label_valid.sum()),\n",
    "      int((~investment.label_valid).sum()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b0fa98",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_design(frame, centre, spread):\n",
    "    scaled = (frame[investment_features] - centre) / spread\n",
    "    sector = pd.get_dummies(frame.sector).reindex(columns=investment_sectors, fill_value=0)\n",
    "    stock = pd.get_dummies(frame.ticker).reindex(columns=investment_tickers, fill_value=0)\n",
    "    return np.column_stack([np.ones(len(frame)), scaled, sector, stock]).astype(float)\n",
    "\n",
    "investment_names = [\"intercept\"] + investment_features + investment_sectors + investment_tickers\n",
    "print(\"design columns:\", len(investment_names), \"stock deviations:\", len(investment_tickers))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "828ca9be",
   "metadata": {},
   "source": [
    "### Refit the same normal posterior inside each training window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2e57a1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_posterior(training, stock_tau=0.02, residual_sd=0.10):\n",
    "    centre = training[investment_features].mean()\n",
    "    spread = training[investment_features].std(ddof=0).replace(0, 1)\n",
    "    matrix = investment_design(training, centre, spread)\n",
    "    outcome = training.ret_next.to_numpy()\n",
    "    prior_var = np.r_[0.05**2, np.repeat(0.05**2, 3), np.repeat(0.03**2, 2),\n",
    "                      np.repeat(stock_tau**2, 9)]\n",
    "    precision = np.diag(1/prior_var) + matrix.T @ matrix / residual_sd**2\n",
    "    covariance = np.linalg.solve(precision, np.eye(len(prior_var)))\n",
    "    mean = np.linalg.solve(precision, matrix.T @ outcome / residual_sd**2)\n",
    "    log_evidence = (-len(outcome)/2 * np.log(2*np.pi*residual_sd**2)\n",
    "        + 0.5 * (np.linalg.slogdet(covariance)[1] - np.log(prior_var).sum())\n",
    "        - 0.5 * (outcome @ outcome / residual_sd**2 - mean @ precision @ mean))\n",
    "    return dict(mean=mean, covariance=covariance, centre=centre, spread=spread,\n",
    "                residual_sd=residual_sd, log_evidence=float(log_evidence))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6694e0cd",
   "metadata": {},
   "source": [
    "The next cell defines prediction without relying on a previous section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d504daf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_prediction(model, candidates):\n",
    "    matrix = investment_design(candidates, model[\"centre\"], model[\"spread\"])\n",
    "    mean = matrix @ model[\"mean\"]\n",
    "    mean_var = np.einsum(\"ij,jk,ik->i\", matrix, model[\"covariance\"], matrix)\n",
    "    mean_var = np.maximum(mean_var, 1e-15)\n",
    "    predictive_var = mean_var + model[\"residual_sd\"]**2\n",
    "    return mean, mean_var, predictive_var"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cce7e40d",
   "metadata": {},
   "source": [
    "### Mixture uncertainty includes model disagreement\n",
    "\n",
    "$$p(r_*\\mid D)=\\sum_kw_k\\,p(r_*\\mid D,M_k),\\qquad \\bar m=\\sum_kw_km_k.$$\n",
    "$$\\operatorname{Var}(r_*\\mid D)=\\sum_kw_k(v_k+m_k^2)-\\bar m^2.$$\n",
    "\n",
    "A weighted average of standard deviations is not this variance. Averaging only $v_k$ also misses disagreement between the component means.\n",
    "\n",
    "Mixture tail probabilities average component tail probabilities; the mixture is generally not Normal.\n",
    "\n",
    "There are now three uncertainty sources: uncertain fitted parameters within each model, new return noise and disagreement across models. Each component predictive variance already includes the first two; the law of total variance adds the third through the squared component means. This decomposition prevents the false precision that arises when the best-looking model is treated as certain.\n",
    "\n",
    "The mean and variance summarise a mixture without determining its complete distribution. In particular, a Normal density with these matching moments is an approximation. For log predictive scoring we evaluate the actual weighted component densities. Predictive stacking, which estimates combination weights for predictive performance, is a different procedure from the evidence weights used here; random-row leave-one-out is also not automatically appropriate for a chronological investment question.\n",
    "\n",
    "### Use logs to normalise model weights safely"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "419f7db7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_mixture(training, candidates):\n",
    "    models = [investment_posterior(training, tau) for tau in [0.005, 0.02, 0.08]]\n",
    "    log_weight = np.array([m[\"log_evidence\"] for m in models]) + np.log(1/3)\n",
    "    weights = np.exp(log_weight - logsumexp(log_weight))\n",
    "    components = [investment_prediction(m, candidates) for m in models]\n",
    "    means = np.stack([c[0] for c in components])\n",
    "    mean_vars = np.stack([c[1] for c in components])\n",
    "    return_vars = np.stack([c[2] for c in components])\n",
    "    mean = weights @ means\n",
    "    predictive_var = weights @ (return_vars + means**2) - mean**2\n",
    "    p_mean = weights @ norm.sf((0.002 - means) / np.sqrt(mean_vars))\n",
    "    p_return = weights @ norm.sf((0.002 - means) / np.sqrt(return_vars))\n",
    "    return mean, predictive_var, p_mean, p_return, weights"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "771808cb",
   "metadata": {},
   "source": [
    "`logsumexp` prevents underflow when multiplying many small densities. The returned probabilities distinguish a conditional-mean hurdle from an actual return hurdle.\n",
    "\n",
    "\n",
    "\n",
    "### Read model weights within their candidate set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbc81298",
   "metadata": {},
   "outputs": [],
   "source": [
    "investment_train = investment.loc[investment.date.lt(\"2021-01-01\")\n",
    "    & investment.label_available_at.lt(\"2021-01-01\") & investment.label_valid].copy()\n",
    "investment_snapshot = investment.loc[investment.date.eq(\"2021-01-31\")].copy()\n",
    "mix_mean, mix_var, mix_p_mean, mix_p_return, mix_weights = investment_mixture(\n",
    "    investment_train, investment_snapshot)\n",
    "print(\"model weights for stock tau 0.005 / 0.02 / 0.08:\", np.round(mix_weights, 4))\n",
    "print(\"sum of model probabilities:\", round(float(mix_weights.sum()), 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f62d168",
   "metadata": {},
   "source": [
    "The model weights are **0.7728, 0.2264 and 0.0008**. The training sample favours stronger pooling within these assumptions.\n",
    "\n",
    "This does not mean a 77.28% probability that the tightest-prior strategy will make money.\n",
    "\n",
    "\n",
    "\n",
    "### Self-check: the highest evidence is not a profit guarantee\n",
    "\n",
    "\n",
    "**Self-check:** A model has posterior model probability 0.90 within a three-model comparison. What does 0.90 describe?\n",
    "\n",
    "- Its probability of positive net profit next year\n",
    "- Its relative posterior probability within the declared model set and priors\n",
    "- Its probability of being the true causal mechanism\n",
    "- Its expected annual return\n",
    "\n",
    "**Feedback:** Model probability depends on the candidate set, prior probabilities and integrated likelihood. Future trading profit is a different random quantity.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### A posterior supports a decision only through an objective\n",
    "\n",
    "A general Bayesian decision chooses an action $w$ to maximise:\n",
    "\n",
    "$$E[U(w^{\\top}r_{t+1}-C(w,w_{t^-}))\\mid D_t].$$\n",
    "\n",
    "Utility $U$, predictive joint returns, transaction costs $C$, constraints and drifted holdings all matter.\n",
    "\n",
    "This exercise uses a simpler declared rule: rank posterior mean returns; reserve three 0.33 slots; retain a slot only if **P(conditional mean > 0.2%) ≥ 0.95**. Otherwise keep cash. This is an illustrative policy, not a solved optimal portfolio.\n",
    "\n",
    "The probability threshold is an assumption to investigate, not a result delivered by Bayes' theorem. A mean-return screen ignores risk aversion unless risk enters elsewhere in the policy. A portfolio optimiser would require the joint predictive covariance and explicit exposure, turnover and concentration constraints. The small normal model here is designed to teach the chain from inference to action without claiming to solve that larger allocation problem.\n",
    "\n",
    "The 95% threshold and 0.2% hurdle are fixed for the reported classroom comparison; neither is presented as a tuned optimum. The supplied data have already been explored in the course, so this exercise is an historical illustration rather than a newly untouched performance test. Changing the threshold after viewing its account is a sensitivity experiment and must be labelled as such. Repeated strategy searches need an additional evaluation design and a record of the attempted rules.\n",
    "\n",
    "### Each annual refit sees only labels already available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38190d32",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_walk_forward(frame):\n",
    "    forecasts, audit = [], []\n",
    "    for year in [2022, 2023, 2024]:\n",
    "        freeze = pd.Timestamp(year=year, month=1, day=1)\n",
    "        train = frame.loc[frame.date.lt(freeze) & frame.label_available_at.lt(freeze) & frame.label_valid]\n",
    "        candidates = frame.loc[frame.date.dt.year.eq(year)].copy()\n",
    "        mean, variance, p_mean, p_return, model_weights = investment_mixture(train, candidates)\n",
    "        candidates[\"forecast\"], candidates[\"predictive_sd\"] = mean, np.sqrt(variance)\n",
    "        candidates[\"p_mean_above_cost\"], candidates[\"p_return_above_cost\"] = p_mean, p_return\n",
    "        candidates[\"baseline_forecast\"] = train.ret_next.mean()\n",
    "        candidates[\"trained_through\"] = train.label_available_at.max()\n",
    "        assert candidates.trained_through.lt(freeze).all()\n",
    "        forecasts.append(candidates)\n",
    "        audit.append({\"year\": year, \"train_rows\": len(train), \"predictions\": len(candidates),\n",
    "                      \"small_tau_weight\": model_weights[0], \"large_tau_weight\": model_weights[-1]})\n",
    "    return pd.concat(forecasts).sort_values([\"date\", \"ticker\"]), pd.DataFrame(audit)\n",
    "investment_walk, investment_audit = investment_walk_forward(investment)\n",
    "print(investment_audit.round(4))\n",
    "print(\"predicted / evaluated / pending:\", len(investment_walk),\n",
    "      int(investment_walk.label_valid.sum()), int((~investment_walk.label_valid).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "723f80ee",
   "metadata": {},
   "source": [
    "Training row counts are **315, 423 and 531**. The three annual blocks produce **315 predictions: 306 evaluable, nine pending**.\n",
    "\n",
    "Features and priors are transformed inside the training window. January's freeze uses labels strictly before 1 January. Refit from the original prior on the complete expanding history, avoiding double-counting.\n",
    "\n",
    "\n",
    "\n",
    "### Refitting and sequential updating count evidence differently\n",
    "\n",
    "\n",
    "**Self-check:** You retain last year’s posterior and also refit on every historical observation, including last year’s data. What is the problem?\n",
    "\n",
    "- No problem: more repetitions always improve precision\n",
    "- Old observations are counted again through both the prior and likelihood\n",
    "- The posterior must become uniform\n",
    "- Walk-forward prohibits any refitting\n",
    "\n",
    "**Feedback:** Use the original prior with all eligible history, or the previous posterior with genuinely new likelihood contributions. Combining the previous posterior with the same old data counts evidence twice.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Score predictions before reading the portfolio curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d32caba",
   "metadata": {},
   "outputs": [],
   "source": [
    "investment_realised = investment_walk.loc[investment_walk.label_valid].copy()\n",
    "log_score = norm.logpdf(investment_realised.ret_next, investment_realised.forecast,\n",
    "                       investment_realised.predictive_sd)\n",
    "# Score the actual mixture, not a Normal approximation with matching moments.\n",
    "mixture_scores = []\n",
    "for year in [2022, 2023, 2024]:\n",
    "    freeze = pd.Timestamp(year=year, month=1, day=1)\n",
    "    training = investment.loc[investment.date.lt(freeze) & investment.label_available_at.lt(freeze)\n",
    "                              & investment.label_valid]\n",
    "    observed = investment_realised.loc[investment_realised.date.dt.year.eq(year)]\n",
    "    models = [investment_posterior(training, tau) for tau in [0.005, 0.02, 0.08]]\n",
    "    lw = np.array([m[\"log_evidence\"] for m in models]); lw -= logsumexp(lw)\n",
    "    component_logpdf = []\n",
    "    for model in models:\n",
    "        mean, _, variance = investment_prediction(model, observed)\n",
    "        component_logpdf.append(norm.logpdf(observed.ret_next, mean, np.sqrt(variance)))\n",
    "    mixture_scores.extend(logsumexp(np.stack(component_logpdf) + lw[:, None], axis=0))\n",
    "rmse = np.sqrt(np.mean((investment_realised.ret_next - investment_realised.forecast)**2))\n",
    "print(\"RMSE / exact mixture mean log score:\", round(float(rmse), 4),\n",
    "      round(float(np.mean(mixture_scores)), 4))\n",
    "print(\"training-mean baseline RMSE:\", round(float(np.sqrt(np.mean(\n",
    "    (investment_realised.ret_next - investment_realised.baseline_forecast)**2))), 4))\n",
    "print(\"pending values excluded:\", int((~investment_walk.label_valid).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f11a1e2",
   "metadata": {},
   "source": [
    "The mixture has **RMSE 0.0952**, versus **0.0930** for a training-mean baseline. Its exact-mixture mean log predictive density is **0.9315** in decimal-return units; higher is better for comparable models on the same outcomes.\n",
    "\n",
    "This forecast does not beat the simple baseline on RMSE. A favourable portfolio ranking would not reverse that finding.\n",
    "\n",
    "RMSE evaluates the posterior predictive mean under squared-error loss. Log predictive density evaluates the complete distribution, including dispersion and tails. Density scores depend on measurement units, so the value 0.9315 is not an accuracy percentage and cannot be compared with a score computed after rescaling returns into percentages. A baseline log score under the same units and observations is needed to interpret predictive-density improvements.\n",
    "\n",
    "An individual stock prediction is not a portfolio allocation. A model can lose on average squared error while a ranking-based rule happens to select a profitable subset, and the converse is also possible. Report the forecast metrics and the economic policy metrics separately. Excluding pending labels is necessary for either evaluation; excluding candidates before prediction because their future outcome is missing would be a different error.\n",
    "\n",
    "### Build a holdings rule without changing the candidate spine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93670428",
   "metadata": {},
   "outputs": [],
   "source": [
    "investment_walk[\"rank\"] = investment_walk.groupby(\"date\").forecast.rank(ascending=False, method=\"first\")\n",
    "investment_walk[\"rank_weight\"] = np.where(investment_walk[\"rank\"].le(3), 0.33, 0.0)\n",
    "investment_walk[\"bayes_weight\"] = investment_walk.rank_weight.where(\n",
    "    investment_walk.p_mean_above_cost.ge(0.95), 0.0)\n",
    "investment_walk[\"equal_weight\"] = 0.99 / 9\n",
    "investment_exposure = investment_walk.groupby(\"date\")[[\"rank_weight\", \"bayes_weight\", \"equal_weight\"]].sum()\n",
    "assert investment_exposure.le(0.99 + 1e-12).all().all()\n",
    "print(\"mean exposure:\\n\", investment_exposure.mean().round(4))\n",
    "print(\"blocked selected positions:\", int((investment_walk.rank_weight.gt(0)\n",
    "    & investment_walk.bayes_weight.eq(0)).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dcf37ae",
   "metadata": {},
   "source": [
    "A permitted name retains its original 0.33 allocation. Blocked slots stay in cash. The equal-weight comparator invests 0.99 across all nine names.\n",
    "\n",
    "Exposure, costs and drawdown must accompany returns. A lower-exposure portfolio is not directly evidence of better stock selection.\n",
    "\n",
    "\n",
    "\n",
    "### Carry drifted holdings and cash through the ledger\n",
    "\n",
    "The loop below represents changing portfolio state. Cash earns zero, trades cost 10 bps per dollar bought or sold, opening NAV is 100, and there is no terminal liquidation. The 1% target cash reserve funds the specified fees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00052b6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def investment_account(weights, returns, cost_bps=10):\n",
    "    holdings = pd.Series(0.0, index=weights.columns)\n",
    "    cash, rows = 100.0, []\n",
    "    for date in weights.index:\n",
    "        opening_nav = cash + holdings.sum()\n",
    "        target = weights.loc[date] * opening_nav\n",
    "        fee = (target - holdings).abs().sum() * cost_bps / 10000\n",
    "        cash = opening_nav - target.sum() - fee\n",
    "        assert cash >= -1e-10\n",
    "        pnl = target * returns.loc[date]\n",
    "        holdings = target + pnl\n",
    "        closing_nav = cash + holdings.sum()\n",
    "        rows.append({\"date\": date, \"opening_nav\": opening_nav, \"nav\": closing_nav,\n",
    "                     \"pnl\": pnl.sum(), \"fee\": fee})\n",
    "    ledger = pd.DataFrame(rows).set_index(\"date\")\n",
    "    assert np.allclose(ledger.nav - ledger.opening_nav, ledger.pnl - ledger.fee)\n",
    "    ledger[\"net_return\"] = ledger.nav / ledger.opening_nav - 1\n",
    "    return ledger"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "384796be",
   "metadata": {},
   "source": [
    "At every month: **closing NAV − opening NAV = security P&L − fees**. A positive prediction never substitutes for a realised return.\n",
    "\n",
    "Turnover is calculated from current holdings after the previous period's gains and losses. Subtracting the previous target weights would omit that drift. Targets in this exercise are fractions of opening NAV before the current trade fees; cash then absorbs fees and the unused budget. The account assumes the supplied simple returns can be earned on those target positions. Execution timing, spread, capacity and source accuracy need evidence beyond the accounting identity.\n",
    "\n",
    "Only mature, valid monthly outcomes enter the realised ledger. Nine pending predictions remain in the forecast output, while their unknown holding period stays outside realised accounting. The final NAV is marked at the end of the last observed holding period without an assumed sale of all holdings. That endpoint convention matters if a different study includes terminal transaction costs.\n",
    "\n",
    "### Compare all three complete accounts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11f64f38",
   "metadata": {},
   "outputs": [],
   "source": [
    "investment_realised = investment_walk.loc[investment_walk.label_valid].copy()\n",
    "realised_returns = investment_realised.pivot(index=\"date\", columns=\"ticker\", values=\"ret_next\")\n",
    "assert realised_returns.notna().all().all()\n",
    "investment_accounts, account_rows = {}, []\n",
    "for name in [\"equal_weight\", \"rank_weight\", \"bayes_weight\"]:\n",
    "    weights = investment_realised.pivot(index=\"date\", columns=\"ticker\", values=name)\n",
    "    ledger = investment_account(weights, realised_returns)\n",
    "    investment_accounts[name] = ledger\n",
    "    wealth = np.r_[100.0, ledger.nav.to_numpy()]\n",
    "    account_rows.append({\"policy\": name, \"ending_nav\": ledger.nav.iloc[-1],\n",
    "        \"max_drawdown\": (wealth / np.maximum.accumulate(wealth) - 1).min(),\n",
    "        \"fees\": ledger.fee.sum(), \"mean_exposure\": weights.sum(axis=1).mean()})\n",
    "print(pd.DataFrame(account_rows).round(4))\n",
    "print(\"realised monthly accounts:\", len(ledger))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7e635a9",
   "metadata": {},
   "source": [
    "With opening NAV 100, the net ending values are **183.90 equal-weight, 204.77 ranking alone and 171.27 Bayesian screening**. Maximum drawdowns are **18.15%, 34.26% and 34.26%**, respectively. Exposure in this table averages the same 34 realised months as the ledger; the earlier signal summary covers all 35 prediction months.\n",
    "\n",
    "The Bayesian threshold is evaluated as part of the whole policy. It does not inherit a performance claim merely from using posterior probabilities.\n",
    "\n",
    "\n",
    "\n",
    "### Read the incremental result, including a null result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dc90612",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "for name, ledger in investment_accounts.items():\n",
    "    ax.plot(ledger.index + pd.offsets.MonthEnd(1), ledger.nav, label=name)\n",
    "ax.set(xlabel=\"Outcome month-end\", ylabel=\"Net NAV; opening 100\")\n",
    "ax.legend(); plt.show()\n",
    "paired_difference = (investment_accounts[\"bayes_weight\"].net_return\n",
    "                     - investment_accounts[\"rank_weight\"].net_return)\n",
    "print(\"mean paired monthly return difference:\", round(float(paired_difference.mean()), 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "466b2203",
   "metadata": {},
   "source": [
    "Compare the threshold policy with ranking alone using **paired monthly net returns**, then apply Chapter 3's dependence-aware uncertainty methods. Thirty-four periods and a survivor sample provide limited evidence.\n",
    "\n",
    "Here the screen subtracts about **0.55 percentage points per month** from ranking alone on paired net returns. It also changes exposure. Forecast RMSE 0.0952 is worse than the training-mean baseline 0.0930. These results support further investigation or rejection of this rule, rather than a claim that the Bayesian screen adds value. The examined historical sample is not untouched evidence.\n",
    "\n",
    "\n",
    "\n",
    "### Changing future outcomes must not change the frozen forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27a76189",
   "metadata": {},
   "outputs": [],
   "source": [
    "altered_future = investment.copy()\n",
    "altered_future.loc[altered_future.date.ge(\"2022-01-01\"), \"ret_next\"] = 0.90\n",
    "original_first = investment_walk.loc[investment_walk.date.dt.year.eq(2022), \"forecast\"]\n",
    "altered_walk, _ = investment_walk_forward(altered_future)\n",
    "np.testing.assert_allclose(original_first,\n",
    "    altered_walk.loc[altered_walk.date.dt.year.eq(2022), \"forecast\"])\n",
    "print(\"PASS: changing 2022+ outcomes cannot change the 2022 frozen forecasts\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c017f0a3",
   "metadata": {},
   "source": [
    "The altered future does change later training samples when those outcomes become eligible. It must not change predictions from the earlier frozen fit.\n",
    "\n",
    "This implementation check does not verify historical source availability or validate the chosen likelihood.\n",
    "\n",
    "\n",
    "\n",
    "### Repair exercise: keep pending labels outside the score\n",
    "\n",
    "The online starter includes every prediction in an RMSE denominator. The worked repair below fixes evaluation rows without removing pending prediction candidates.\n",
    "**Repair task:** Keep all forecasts; restrict evaluation to valid outcomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54c0dfa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Worked repair: retain pending forecasts, but score only valid labels.\n",
    "evaluation_rows = investment_walk.loc[investment_walk.label_valid].copy()\n",
    "assert len(evaluation_rows) == 306 and evaluation_rows.label_valid.all()\n",
    "print(\"evaluation rows:\", len(evaluation_rows), \"pending included:\", int((~evaluation_rows.label_valid).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ad0d8c7",
   "metadata": {},
   "source": [
    "**Worked repair:** `evaluation_rows = investment_walk.loc[investment_walk.label_valid].copy()`. The online starter reports the wrong information set; the repair preserves the separate quantities.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### A research memo must separate signal, model and policy\n",
    "\n",
    "Report the **target and clock; likelihood and priors; pooling assumptions; model set and evidence weights; parameter versus predictive uncertainty; forecast baseline; frozen holding rule; exposure, costs and net accounts**.\n",
    "\n",
    "State which result would change your decision. Include prior sensitivity and unsuccessful comparisons. A posterior probability describes the specified model, rather than certifying economic profitability.\n",
    "\n",
    "\n",
    "\n",
    "### What you discovered\n",
    "\n",
    "- Bayesian inference changes signals through pooling and combines uncertainty through predictive mixtures.\n",
    "- Evidence weights are conditional on the declared model set and proper priors.\n",
    "- Walk-forward trains and evaluates that model process at explicit time boundaries.\n",
    "- Prediction quality, portfolio profit and risk are separate results to report.\n",
    "\n",
    "**Next:** Chapter 6 develops time-series dependence and dynamic risk models under the same availability discipline."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfb0f95d",
   "metadata": {},
   "source": [
    "## Explain your decision\n",
    "\n",
    "The screening rule ends below both comparator portfolios in this example, and its RMSE is worse than the training-mean baseline. Explain why posterior model probabilities do not certify profitability. Report dates, sample, exposure, costs, drawdown and the paired return difference before proposing an extension.\n",
    "\n",
    "Try a declared prior sensitivity or dependence-aware likelihood, preserve the original result and evaluate any change on a genuinely fresh period or universe. Do not rename this explored sample as untouched validation.\n",
    "\n",
    "See [the companion Book](https://statpython-book.pages.dev/topic5.html) for Bayesian foundations, dynamic risk and the interactive prior explorer, and [Feng and He](https://arxiv.org/abs/1902.01015) for a fuller investment hierarchy.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
