{ "cells": [ { "cell_type": "markdown", "id": "8ecbcb9f", "metadata": {}, "source": [ "# Chapter 3 \u2014 Strategy Diagnostics\n", "\n", "Chapter 2 built the research tables. Chapter 3 asks what those tables justify. A valid join, an available feature and a reconciled NAV are necessary foundations; the next task is to measure a signal, quantify its uncertainty and describe the losses behind its returns. The running case uses the actual public teaching outputs of Chapter 2, so students can trace a statistic to a key, a date and a denominator.\n", "\n", "There are **four classroom units**: distribution diagnosis, signal IC, backtest evidence and tail risk. The five existing section numbers and links are retained: take \u00a73.1 \u2192 \u00a73.4 \u2192 \u00a73.2\u20133.3 \u2192 \u00a73.5 in class. Section 3.3 supplies research discipline within the evidence unit. This is one three-hour session, not five lectures or a requirement to present every example.\n", "\n", "The classroom handoff contains 639 candidates, 630 valid labels, 315 keyed predictions and 34 realised monthly portfolio returns. Nine predictions have pending labels. The fixed Ridge baseline selects the three highest scores, targets 99% equal-weight long-only investment and charges 10 basis points per dollar traded. The small current-survivor sample uses assumed data clocks and static sectors. It illustrates diagnosis and reproducibility; it does not evaluate a full MyStrategy long-short implementation.\n", "\n", "Read the core subsections first. Collapsed advanced-reading panels preserve the original Apple, index, Dow, A/B and Bandit worked examples. They develop the theory and additional methods without expanding the classroom agenda. The full EVT fits, BCa details and experiment-design extensions are follow-up material, rather than prerequisites for completing the group report.\n", "\n", "### Class brief: a strategy diagnostic report\n", "\n", "**Decision question:** Does the supplied evidence justify further research into this fixed baseline, and what evidence is still missing?\n", "\n", "**Before class:** Review the core cells in \u00a7\u00a73.1, 3.4, 3.2\u20133.3 and 3.5. Download the [guided notebook](https://statpython.pages.dev/data/Chapter3_StrategyDiagnostics.ipynb). Predict one IC interpretation, one interval interpretation and one risk limitation; bring a compact note. Use the precomputed handoff rather than train or tune another model.\n", "\n", "**During class:** Work in groups of at most three; rotate analyst, validation reviewer and decision communicator. Predict, run supplied code, change one diagnostic setting, and challenge the interpretation. The times below are student activity blocks; teacher explanation and feedback are brief and targeted.\n", "\n", "| Minutes | Student activity | Evidence produced |\n", "|---|---|---|\n", "| 0\u201310 | Individual retrieval and initial judgement | One claim to test |\n", "| 10\u201340 | Distribution diagnosis (\u00a73.1) | Quantiles, ECDF, fixed training rule |\n", "| 40\u201380 | Signal diagnosis (\u00a73.4) | Date-level IC and coverage |\n", "| 80\u201390 | Break | \u2014 |\n", "| 90\u2013130 | Uncertainty and research discipline (\u00a73.2\u20133.3) | Block interval and search log |\n", "| 130\u2013155 | Risk diagnosis (\u00a73.5) | Tail count, ES, drawdown |\n", "| 155\u2013180 | Group report, peer challenge and revision | One runnable notebook and short memo |\n", "\n", "**Deliverable:** One runnable notebook with a short decision memo: question/data contract; IC and coverage; uncertainty and search history; risk horizon and tail count; recommendation, limitation and next test. Include each member's explanation of one code/output choice. Compare with another group and revise before submitting through the existing course submission route. This activity adds no new grading or attendance rules.\n", "\n", "**Check quality:** Keys and clocks are explicit; computations reproduce; units and denominators match the claim; limitations affect the recommendation. Follow-up work may explore alternative block lengths, nonlinear dependence or EVT, clearly labelled exploratory." ] }, { "cell_type": "markdown", "id": "d823002a", "metadata": {}, "source": [ "## Feature Distributions and Diagnostics" ] }, { "cell_type": "markdown", "id": "cc8d3e80", "metadata": {}, "source": [ "### Statistical foundations \u2014 explain before running\n", "\n", "#### Population, sample and the quantity we estimate\n", "\n", "A **population distribution** describes the random variable of interest; a **sample** supplies observations. Specify the unit: stock-month, portfolio-month or trading day.\n", "\n", "$$\\mu=E[X],\\qquad \\bar x=\\frac{1}{n}\\sum_i x_i,\\qquad s^2=\\frac{1}{n-1}\\sum_i(x_i-\\bar x)^2.$$\n", "\n", "The parameter $\\mu$ is unknown; $\\bar x$ is an estimator that changes across samples. Under i.i.d. sampling, $s^2$ is unbiased for variance; $s$ need not be unbiased for standard deviation.\n", "\n", "**Explain:** Why can 639 rows describe a distribution without providing 639 independent observations?\n", "\n", "#### CDF, ECDF and quantiles\n", "\n", "$$F(x)=P(X\\le x),\\qquad \\widehat F_n(x)=\\frac{1}{n}\\sum_i\\mathbf1\\{x_i\\le x\\}.$$\n", "\n", "The CDF maps a value to a cumulative probability. The ECDF replaces that probability with the observed fraction; its staircase jumps at sample values.\n", "\n", "$$q_p=\\inf\\{x:F(x)\\ge p\\}.$$\n", "\n", "The median is $q_{0.5}$. Software sample quantiles may interpolate between order statistics, so record the convention. A 95th percentile is a location in a distribution, not a confidence interval.\n", "\n", "#### Density estimation: histogram, KDE and bandwidth\n", "\n", "For a continuous variable, probability is **area**, not the height of its density:\n", "\n", "$$P(a2$.\n", "\n", "A **Q\u2013Q plot** compares sample and reference quantiles. Systematic tail departures question the reference model; a straight centre does not validate extreme tails.\n", "\n", "The KS distance is $D=\\sup_x|\\widehat F_n(x)-F_0(x)|$. Its standard calibration assumes independent observations and a specified continuous $F_0$. Fitting parameters on the same sample changes the null distribution; use a suitable calibrated test or simulation.\n", "\n", "**Check:** A descriptive distribution fit is not proof of independent or stationary returns." ] }, { "cell_type": "markdown", "id": "65e6e464", "metadata": {}, "source": [ "### The Chapter 2 handoff\n", "\n", "Chapter 2 supplied keyed features, mature labels, walk-forward predictions and a reconciled account. Read the exported tables rather than rebuild a model in this chapter. The nine-stock teaching sample contains 639 candidates, 630 valid labels and 315 predictions, including nine predictions whose labels are still pending. A missing future label must not remove a current candidate." ] }, { "cell_type": "code", "execution_count": null, "id": "aa1dac8f", "metadata": {}, "outputs": [], "source": [ "import numpy as np, pandas as pd\n", "import matplotlib.pyplot as plt\n", "data_url = \"https://statpython.pages.dev/data/\"\n", "panel = pd.read_csv(data_url + \"chapter2_scored.csv\",\n", " parse_dates=[\"date\", \"label_available_at\", \"trained_through\"])\n", "assert not panel.duplicated([\"date\", \"ticker\"]).any()\n", "assert panel.loc[panel.prediction.notna(), \"trained_through\"].le(\n", " panel.loc[panel.prediction.notna(), \"date\"]).all()\n", "print(len(panel), \"candidates;\", int(panel.label_valid.sum()), \"valid labels\")" ] }, { "cell_type": "markdown", "id": "eb7464c6", "metadata": {}, "source": [ "Keys and training cutoffs are checked before any statistic. The manifest records the source and file hashes. This public current-survivor sample uses assumed availability clocks and static sectors; it demonstrates methods, not historical investability." ] }, { "cell_type": "markdown", "id": "74a1ac0e", "metadata": {}, "source": [ "### Feature distributions across a fixed split\n", "\n", "Use 2019\u20132021 as the development period and 2022 onward as the held-out diagnostic period. Inspect the distribution before choosing a transform. The split is fixed for teaching; looking at the later period must not trigger repeated tuning presented as untouched OOS evidence." ] }, { "cell_type": "code", "execution_count": null, "id": "675afea9", "metadata": {}, "outputs": [], "source": [ "feature = \"ret_12_1\"\n", "panel[\"sample\"] = np.where(panel.date.lt(\"2022-01-01\"), \"development\", \"held-out\")\n", "distribution = panel.groupby(\"sample\")[feature].agg(n=\"count\", mean=\"mean\", sd=\"std\")\n", "quantiles = panel.groupby(\"sample\")[feature].quantile([0.05, 0.5, 0.95]).unstack()\n", "print(distribution.round(3))\n", "print(quantiles.round(3))" ] }, { "cell_type": "markdown", "id": "15275e84", "metadata": {}, "source": [ "The counts are stock-months, not independent experiments. A distribution difference is a prompt to investigate coverage, composition and market conditions; it does not itself identify a profitable signal or a causal mechanism." ] }, { "cell_type": "markdown", "id": "115fbf9c", "metadata": {}, "source": [ "### ECDF and a training-only clipping rule\n", "\n", "An ECDF answers what proportion of observations lie below a value without choosing histogram bins. Derive an illustrative clipping rule from development data only. Apply the same cutoffs to later observations; retain the raw column and report the fraction affected. This rule is a diagnostic experiment, not a recommendation to clip every feature." ] }, { "cell_type": "code", "execution_count": null, "id": "f22a197f", "metadata": {}, "outputs": [], "source": [ "development = panel.loc[panel[\"sample\"].eq(\"development\"), feature].dropna()\n", "held_out = panel.loc[panel[\"sample\"].eq(\"held-out\"), feature].dropna()\n", "lower, upper = development.quantile([0.01, 0.99])\n", "clipped = held_out.clip(lower=lower, upper=upper)\n", "print(\"training cutoffs:\", round(lower, 3), round(upper, 3))\n", "print(\"held-out fraction affected:\", round(float(clipped.ne(held_out).mean()), 3))\n", "fig, ax = plt.subplots(figsize=(9, 4))\n", "ax.step(np.sort(development), np.arange(1, len(development)+1)/len(development), label=\"development\")\n", "ax.step(np.sort(held_out), np.arange(1, len(held_out)+1)/len(held_out), label=\"held-out\")\n", "ax.set(xlabel=feature, ylabel=\"Empirical cumulative probability\")\n", "ax.legend(); plt.show()" ] }, { "cell_type": "markdown", "id": "2102a048", "metadata": {}, "source": [ "The staircase is descriptive even when observations are dependent. The classical i.i.d. DKW confidence bound in the advanced reading cannot automatically certify a panel of stock-months. KDE, fitted Student-t and Q-Q plots remain useful follow-up diagnostics." ] }, { "cell_type": "markdown", "id": "7b2b4a08", "metadata": {}, "source": [ "## Signal Association and Cross-sectional IC" ] }, { "cell_type": "markdown", "id": "a1b743bb", "metadata": {}, "source": [ "### Statistical foundations \u2014 explain before running\n", "\n", "#### Covariance, Pearson correlation and scale\n", "\n", "$$s_{xy}=\\frac{1}{n-1}\\sum_i(x_i-\\bar x)(y_i-\\bar y),\\qquad r_{xy}=\\frac{s_{xy}}{s_xs_y}.$$\n", "\n", "Covariance retains the variables\u2019 units. Pearson correlation removes their scale and summarises linear association between the same valid pairs; it lies between $-1$ and $1$.\n", "\n", "Constant inputs make correlation undefined. Outliers, missing-pair selection and pooling different periods can change the answer.\n", "\n", "**Explain:** A correlation of 0.6 does not mean a 60% probability of a correct trade or a 60% return.\n", "\n", "#### Rank association: Spearman and Kendall\n", "\n", "**Spearman** is Pearson correlation of both rank vectors:\n", "\n", "$$\\rho_s=\\operatorname{Corr}(\\operatorname{rank}X,\\operatorname{rank}Y).$$\n", "\n", "It measures monotonic association and is unchanged by strictly increasing transformations. Use average ranks for ties and a common valid-pair sample.\n", "\n", "**Kendall** compares concordant and discordant pairs; without ties, $\\tau=(C-D)/\\binom n2$. Tau-b adjusts for ties.\n", "\n", "An IC is an association statistic applied to predictions and realised targets. This chapter uses per-date Spearman IC; state that definition whenever reporting IC.\n", "\n", "#### Zero correlation, independence and an IC time series\n", "\n", "If $X$ is symmetric about zero and $Y=X^2$, their covariance can be zero even though $Y$ is determined by $X$. Zero Pearson or Spearman correlation does not establish independence.\n", "\n", "Distance correlation can detect broader dependence; it still does not establish causation or a tradable effect. Scatter plots and subject knowledge complement one-number summaries.\n", "\n", "For $T$ date-level IC values:\n", "\n", "$$\\overline{IC}=T^{-1}\\sum_t IC_t,\\qquad ICIR=\\overline{IC}/s_{IC}.$$\n", "\n", "This ICIR is not annualised. Dependence between dates affects uncertainty; inspect coverage, stability and a time-aware interval for mean IC." ] }, { "cell_type": "markdown", "id": "36d1bde2", "metadata": {}, "source": [ "### Define IC on one investable cross-section\n", "\n", "Information coefficient (IC) here means the association between a score known at the signal date and its later realised outcome, computed across stocks at that date. Pearson measures linear association; Spearman applies Pearson to ranks. The central MyStrategy task is to obtain one IC per signal date, then study the resulting time series." ] }, { "cell_type": "code", "execution_count": null, "id": "f67ec8e1", "metadata": {}, "outputs": [], "source": [ "import numpy as np, pandas as pd\n", "import matplotlib.pyplot as plt\n", "data_url = \"https://statpython.pages.dev/data/\"\n", "panel = pd.read_csv(data_url + \"chapter2_scored.csv\",\n", " parse_dates=[\"date\", \"label_available_at\", \"trained_through\"])\n", "assert not panel.duplicated([\"date\", \"ticker\"]).any()\n", "assert panel.loc[panel.prediction.notna(), \"trained_through\"].le(\n", " panel.loc[panel.prediction.notna(), \"date\"]).all()\n", "print(len(panel), \"candidates;\", int(panel.label_valid.sum()), \"valid labels\")" ] }, { "cell_type": "markdown", "id": "b20e7747", "metadata": {}, "source": [ "Raw correlation pooled across all dates mixes market periods and cross-sectional differences. It is not the same quantity as the mean of per-date IC. Neither a strong IC nor a rank transform establishes tradable profits after turnover and fees." ] }, { "cell_type": "markdown", "id": "7eb7f51b", "metadata": {}, "source": [ "### Keep prediction eligibility separate from evaluation\n", "\n", "Evaluate only prediction rows whose outcomes are mature and finite. Pending rows remain in the candidate/prediction table. For a sector-relative diagnostic, subtract each date-sector mean realised return after the outcome is known; that transformed outcome is evaluation data, never a feature available at the signal date." ] }, { "cell_type": "code", "execution_count": null, "id": "9f159fc5", "metadata": {}, "outputs": [], "source": [ "evaluation = panel.loc[panel.prediction.notna() & panel.label_valid].copy()\n", "assert np.isfinite(evaluation[[\"prediction\", \"ret_next\"]]).all().all()\n", "evaluation[\"sector_target\"] = evaluation.ret_next - evaluation.groupby(\n", " [\"date\", \"sector\"]).ret_next.transform(\"mean\")\n", "counts = evaluation.groupby(\"date\").size()\n", "print(len(evaluation), \"evaluated predictions;\", len(counts), \"signal dates\")\n", "print(\"stocks per date:\", int(counts.min()), \"to\", int(counts.max()))\n", "print(\"pending predictions retained:\", int((panel.prediction.notna() & ~panel.label_valid).sum()))" ] }, { "cell_type": "markdown", "id": "b28904bf", "metadata": {}, "source": [ "The handoff gives 306 evaluated predictions over 34 signal dates, nine stocks per date, and nine pending predictions retained. Sector-relative targets change the question; they do not prove neutral portfolio exposures. Three stocks per sector are a deliberately small classroom example." ] }, { "cell_type": "markdown", "id": "e028ee43", "metadata": {}, "source": [ "### Compute Spearman IC and ICIR by date\n", "\n", "Rank both the prediction and each target within the same valid date cross-section, using average ranks for ties. Correlate those ranks within the date. Use the dates as the observations when summarising stability; also retain coverage and investigate constant-score or constant-target dates rather than filling undefined correlations with zero." ] }, { "cell_type": "code", "execution_count": null, "id": "64b0d035", "metadata": {}, "outputs": [], "source": [ "ranked = evaluation.groupby(\"date\")[[\"prediction\", \"ret_next\", \"sector_target\"]].rank()\n", "date_correlations = ranked.groupby(evaluation.date).corr()\n", "ic = date_correlations.xs(\"prediction\", level=1)[[\"ret_next\", \"sector_target\"]]\n", "ic.columns = [\"raw_ic\", \"sector_ic\"]\n", "summary = ic.agg([\"count\", \"mean\", \"std\"]).T\n", "summary[\"icir\"] = summary[\"mean\"] / summary[\"std\"]\n", "print(summary.round(3))\n", "print(\"positive-date fraction:\", ic.gt(0).where(ic.notna()).mean().round(3).to_dict())" ] }, { "cell_type": "markdown", "id": "ab93faff", "metadata": {}, "source": [ "Mean raw IC is 0.044 and sector-relative IC 0.021 across 34 dates. ICIR is mean IC divided by its sample standard deviation, without annualisation. These are different targets; use consistent keys, valid pairs and tie handling when comparing versions." ] }, { "cell_type": "markdown", "id": "c462bed1", "metadata": {}, "source": [ "### Stability, redundancy and nonlinear dependence\n", "\n", "Inspect annual IC and the full dated series before reporting the average. A favourable year can conceal a later reversal. Feature redundancy should be investigated within development data, and distance correlation can screen non-monotonic dependence; it has no sign and does not by itself show a useful ranking direction." ] }, { "cell_type": "code", "execution_count": null, "id": "7c9ca636", "metadata": {}, "outputs": [], "source": [ "print(ic.groupby(ic.index.year).mean().round(3))\n", "fig, ax = plt.subplots(figsize=(9, 4))\n", "ic.plot(ax=ax, marker=\".\", linewidth=1)\n", "ax.axhline(0, color=\"black\", linewidth=1)\n", "ax.set(ylabel=\"Monthly cross-sectional Spearman IC\", xlabel=\"Signal date\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0ce2bdb3", "metadata": {}, "source": [ "An average is a summary of 34 noisy cross-sections. The uncertainty tools in \u00a73.2 apply to this IC time series too. The companion notebook estimates a block-bootstrap interval for mean raw IC and compares it with the net-return evidence." ] }, { "cell_type": "markdown", "id": "7538465c", "metadata": {}, "source": [ "## Uncertainty and Backtest Evidence" ] }, { "cell_type": "markdown", "id": "491c3137", "metadata": {}, "source": [ "### Statistical foundations \u2014 explain before running\n", "\n", "#### Sampling variation, bias and standard error\n", "\n", "The **sampling distribution** is the distribution of an estimator across repeated samples. It differs from the distribution of individual returns.\n", "\n", "$$\\operatorname{Bias}(\\widehat\\theta)=E[\\widehat\\theta]-\\theta,\\qquad \\operatorname{MSE}(\\widehat\\theta)=\\operatorname{Var}(\\widehat\\theta)+\\operatorname{Bias}(\\widehat\\theta)^2.$$\n", "\n", "Standard deviation describes observations; **standard error** describes uncertainty in an estimate. Under i.i.d. sampling, $\\widehat{SE}(\\bar X)=s/\\sqrt n$.\n", "\n", "**Explain:** More stock rows in one month need not reduce uncertainty as much as more independent months.\n", "\n", "#### Confidence intervals: the meaning of 95%\n", "\n", "For an i.i.d. Normal sample with unknown variance, a mean interval is\n", "\n", "$$\\bar x\\ \\pm\\ t_{n-1,\\,0.975}\\frac{s}{\\sqrt n}.$$\n", "\n", "In repeated sampling under the assumptions, this procedure covers the fixed population mean in 95% of samples. It does not assign a 95% probability to the mean being inside the realised interval.\n", "\n", "A confidence interval for the mean is not a prediction interval for a future return. Increasing sample size can reduce sampling error without removing selection bias or data leakage.\n", "\n", "#### Bootstrap: approximate an estimator\u2019s sampling distribution\n", "\n", "1. Choose a resampling unit consistent with the data-generating process.\n", "2. Draw $B$ resamples and recompute the statistic $\\widehat\\theta_b^*$.\n", "3. A percentile interval uses the 2.5th and 97.5th percentiles of these statistics.\n", "\n", "An i.i.d. bootstrap samples observations independently. A block bootstrap samples runs to preserve local time dependence; block length trades dependence preservation against the number of effectively distinct blocks.\n", "\n", "**Assumptions:** Enough data and a sufficiently stable process. More resamples reduce simulation noise; they do not create new historical information. BCa adds bias/acceleration adjustments and is developed in the Book.\n", "\n", "#### Why dependence changes a standard error\n", "\n", "For a stationary series with autocovariance $\\gamma_k$,\n", "\n", "$$\\operatorname{Var}(\\bar X)=\\frac{1}{n}\\left[\\gamma_0+2\\sum_{k=1}^{n-1}(1-k/n)\\gamma_k\\right].$$\n", "\n", "Positive and negative dependence can change uncertainty in different directions. The i.i.d. formula drops all non-zero-lag terms.\n", "\n", "**HAC** estimates uncertainty allowing heteroskedasticity and autocorrelation up to a chosen lag rule. In $r_t=\\mu+\\varepsilon_t$, the OLS intercept estimates mean return; HAC changes its standard error, not the fitted mean.\n", "\n", "With 34 months, compare lag/block sensitivity and avoid treating asymptotic inference as exact." ] }, { "cell_type": "markdown", "id": "ca6ea167", "metadata": {}, "source": [ "### The unit of evidence is a realised period\n", "\n", "A return estimate needs its frequency, cost basis and sampling unit. Load the 34 realised monthly net returns from Chapter 2. The period end is the outcome date; the separate date column is the signal/rebalance date. Never count a pending next-month outcome as zero." ] }, { "cell_type": "code", "execution_count": null, "id": "0efb4b35", "metadata": {}, "outputs": [], "source": [ "import numpy as np, pandas as pd\n", "import matplotlib.pyplot as plt\n", "data_url = \"https://statpython.pages.dev/data/\"\n", "ledger = pd.read_csv(data_url + \"chapter2_ledger.csv\",\n", " parse_dates=[\"date\", \"period_end\"])\n", "ledger = ledger.set_index(\"period_end\").sort_index()\n", "assert ledger.index.is_unique and ledger.index.is_monotonic_increasing\n", "assert np.allclose(ledger.nav_end - ledger.nav_start,\n", " ledger.gross_pnl - ledger.cost)\n", "assert np.allclose(ledger.net_return, ledger.nav_end / ledger.nav_start - 1)\n", "net = ledger.net_return\n", "print(len(net), \"realised monthly net returns\")" ] }, { "cell_type": "markdown", "id": "f60f8251", "metadata": {}, "source": [ "One portfolio has one realised return per month. Nine stocks observed in each month do not create nine independent portfolio histories. Reconciliation establishes accounting consistency; it does not establish alpha." ] }, { "cell_type": "markdown", "id": "e1250a9f", "metadata": {}, "source": [ "### Resample adjacent months\n", "\n", "An ordinary bootstrap samples individual observations independently. Trading returns and IC can have serial dependence, especially with overlapping horizons. A circular block bootstrap instead samples contiguous runs and wraps the last observation back to the first. The function makes one matrix of indices; no table-processing loop is hidden in it." ] }, { "cell_type": "code", "execution_count": null, "id": "0f672b34", "metadata": {}, "outputs": [], "source": [ "def circular_block_indices(n, length=3, draws=2000, seed=5650):\n", " \"\"\"Sample adjacent months, wrapping at the end; one row per replicate.\"\"\"\n", " assert 1 <= length <= n and draws >= 200\n", " rng = np.random.default_rng(seed)\n", " blocks = int(np.ceil(n / length))\n", " starts = rng.integers(0, n, size=(draws, blocks, 1))\n", " indices = (starts + np.arange(length)) % n\n", " return indices.reshape(draws, -1)[:, :n]" ] }, { "cell_type": "markdown", "id": "550bdf7a", "metadata": {}, "source": [ "Length 1 recovers independent resampling. Length 3 preserves local runs within each three-month block. Blocks are still resampled independently, and stationarity is an assumption: this does not recreate unseen regimes or correct strategy selection." ] }, { "cell_type": "markdown", "id": "6e0790d5", "metadata": {}, "source": [ "### A confidence interval and a sensitivity check\n", "\n", "Estimate uncertainty about the mean monthly net return with a fixed seed and 2,000 replicates. Compare independent resampling with three-month blocks, then change the block length to six. Neither method is required to give a wider interval on this particular sample." ] }, { "cell_type": "code", "execution_count": null, "id": "4c9b1611", "metadata": {}, "outputs": [], "source": [ "iid_means = net.to_numpy()[circular_block_indices(len(net), length=1)].mean(axis=1)\n", "block_means = net.to_numpy()[circular_block_indices(len(net), length=3)].mean(axis=1)\n", "print(\"mean monthly net return:\", round(float(net.mean()), 5))\n", "print(\"iid 95% interval:\", np.quantile(iid_means, [0.025, 0.975]).round(5))\n", "print(\"block-3 95% interval:\", np.quantile(block_means, [0.025, 0.975]).round(5))\n", "block6 = net.to_numpy()[circular_block_indices(len(net), length=6)].mean(axis=1)\n", "print(\"block-6 95% interval:\", np.quantile(block6, [0.025, 0.975]).round(5))" ] }, { "cell_type": "markdown", "id": "71c514d6", "metadata": {}, "source": [ "Mean monthly net return is 0.02494; the block-3 interval is [0.00792, 0.04238], narrower here than iid resampling. Report the 34-month denominator and block sensitivity. An interval excluding zero does not account for strategy search, sample selection or unverified execution assumptions." ] }, { "cell_type": "markdown", "id": "a09b80ee", "metadata": {}, "source": [ "### An intercept-only model with HAC uncertainty\n", "\n", "Reading a fitted constant is enough here: in a regression on a column of ones, the coefficient is the sample mean. HAC changes its estimated uncertainty to allow heteroskedasticity and autocorrelation. A three-lag choice is a teaching setting to challenge, not a validated optimum." ] }, { "cell_type": "code", "execution_count": null, "id": "f1029b19", "metadata": {}, "outputs": [], "source": [ "import statsmodels.api as sm\n", "mean_fit = sm.OLS(net.to_numpy(), np.ones((len(net), 1))).fit(\n", " cov_type=\"HAC\", cov_kwds={\"maxlags\": 3, \"use_correction\": True})\n", "naive_se = net.std(ddof=1) / np.sqrt(len(net))\n", "print(\"iid SE:\", round(float(naive_se), 5))\n", "print(\"HAC SE:\", round(float(mean_fit.bse[0]), 5))\n", "print(\"HAC mean / SE:\", round(float(mean_fit.tvalues[0]), 3))" ] }, { "cell_type": "markdown", "id": "07c93214", "metadata": {}, "source": [ "This is an approximate large-sample calculation with only 34 observations. Compare assumptions and sensitivity rather than treating a threshold such as |t| > 2 as an automatic acceptance rule. Regression modelling itself is developed in Chapter 4." ] }, { "cell_type": "markdown", "id": "a7d12fd6", "metadata": {}, "source": [ "### Apply the same resampling to mean IC\n", "\n", "The earlier IC table has one row per signal date. Resample those dates in blocks; do not treat each stock as an independent replication." ] }, { "cell_type": "code", "execution_count": null, "id": "4858609e", "metadata": {}, "outputs": [], "source": [ "raw_ic = ic.raw_ic.dropna().sort_index()\n", "ic_means = raw_ic.to_numpy()[circular_block_indices(len(raw_ic), length=3)].mean(axis=1)\n", "print(\"mean raw IC:\", round(float(raw_ic.mean()), 3))\n", "print(\"block-3 95% IC interval:\", np.quantile(ic_means, [0.025, 0.975]).round(3))" ] }, { "cell_type": "markdown", "id": "097bf1e2", "metadata": {}, "source": [ "## Research Design and Multiple Testing" ] }, { "cell_type": "markdown", "id": "2a9421ce", "metadata": {}, "source": [ "### Statistical foundations \u2014 explain before running\n", "\n", "#### Hypotheses, p-values and decision errors\n", "\n", "State $H_0$ and $H_1$ before testing. For example, $H_0:\\mu=0$ against $H_1:\\mu\\ne0$ for a fixed strategy and return definition.\n", "\n", "A **p-value** is the probability, under $H_0$ and the test assumptions, of a statistic at least as extreme as observed. It is not $P(H_0\\mid\\text{data})$ or an effect size.\n", "\n", "| Decision error | Meaning |\n", "|---|---|\n", "| Type I | Reject a true $H_0$; controlled by significance level $\\alpha$ for a valid test |\n", "| Type II | Fail to reject a specified false $H_0$; probability $\\beta$ |\n", "\n", "**Power** is $1-\\beta$ for a stated effect and sampling design. Failure to reject does not establish zero effect.\n", "\n", "#### Multiple comparisons: family-wise error and false discoveries\n", "\n", "For $m$ independent true-null tests at level $\\alpha$, the probability of at least one false rejection is $1-(1-\\alpha)^m$. At $m=20$, $\\alpha=0.05$, it is about 64%.\n", "\n", "**Bonferroni:** test each at $\\alpha/m$ to control family-wise error; independence is not required when individual tests are valid.\n", "\n", "**Benjamini\u2013Hochberg:** order the p-values; find the largest $k$ with $p_{(k)}\\le kq/m$, then reject through $k$. This controls expected false-discovery proportion under suitable independence/dependence conditions.\n", "\n", "Record every tried variant. Neither correction fixes invalid p-values, leakage or an undisclosed search.\n", "\n", "#### Research design: paired evidence and a fixed protocol\n", "\n", "Comparing strategies A and B on the same dates is a **paired** problem: analyse $d_t=r_{A,t}-r_{B,t}$, retaining time dependence in uncertainty estimation.\n", "\n", "Random assignment in an experiment helps identify effects; a historical market comparison does not inherit that property. Correlation, predictive improvement and a causal effect are distinct claims.\n", "\n", "Pre-specify the target, effect of interest, features, split, costs, metrics and stopping rule. Repeated peeking or tuning on a holdout changes the inference problem.\n", "\n", "Power planning, A/B testing, sequential testing and Bandits remain worked extensions; adaptive reward maximisation differs from evaluating a fixed rule." ] }, { "cell_type": "markdown", "id": "7f16062d", "metadata": {}, "source": [ "### Fix the research question before reading the answer\n", "\n", "The research protocol fixes a universe, timestamp contract, label horizon, primary metric, costs, model settings and evaluation dates. Record every candidate tried. Comparing strategies on the same dates is a paired comparison: analyse their net-return difference while preserving time dependence. Historical backtests are observational comparisons, not randomised A/B trials." ] }, { "cell_type": "markdown", "id": "b14e1e5b", "metadata": {}, "source": [ "An untouched holdout is used once for a stated decision. Repeatedly revising features after inspecting it turns that period into development data. Sample-size and sequential-testing ideas transfer from experiments, but randomisation assumptions do not automatically transfer to markets." ] }, { "cell_type": "markdown", "id": "5740d25a", "metadata": {}, "source": [ "### Two hundred searches, one attractive result\n", "\n", "This fixed-seed simulation isolates selection risk: 200 independent candidate signals, 150 observations each, ten with a non-zero mean. Run the same family once, then compare an unadjusted screen with Bonferroni and Benjamini\u2013Hochberg. It is a simulation of repeated research, not a measured result from MyStrategy." ] }, { "cell_type": "code", "execution_count": null, "id": "e10a6952", "metadata": {}, "outputs": [], "source": [ "import numpy as np, pandas as pd\n", "from scipy import stats\n", "from statsmodels.stats.multitest import multipletests\n", "rng = np.random.default_rng(9)\n", "signals = np.empty((200, 150))\n", "signals[:10] = rng.normal(0.4, 1.0, size=(10, 150))\n", "signals[10:] = rng.normal(0.0, 1.0, size=(190, 150))\n", "pvalues = stats.ttest_1samp(signals, 0, axis=1).pvalue\n", "naive = pvalues < 0.05\n", "bonferroni = multipletests(pvalues, method=\"bonferroni\")[0]\n", "bh = multipletests(pvalues, method=\"fdr_bh\")[0]\n", "print(\"naive / Bonferroni / BH discoveries:\", naive.sum(), bonferroni.sum(), bh.sum())\n", "print(\"false discoveries:\", naive[10:].sum(), bonferroni[10:].sum(), bh[10:].sum())" ] }, { "cell_type": "markdown", "id": "3595ee6e", "metadata": {}, "source": [ "The unadjusted screen finds 19 signals, nine false; Bonferroni finds nine, all real; BH finds ten, all real in this draw. Real features are dependent, and BH needs suitable dependence conditions. Adjustment requires valid base p-values and the full tested family; neither it nor a holdout repairs leakage." ] }, { "cell_type": "markdown", "id": "21630405", "metadata": {}, "source": [ "### Research log and stopping rule\n", "\n", "For each attempt, retain a dated hypothesis, configuration, source hashes, data cutoff and all evaluation metrics. State whether the outcome continues research, rejects the hypothesis or motivates a separately labelled exploratory change. Do not stop searching merely because the current winner crosses a significance threshold." ] }, { "cell_type": "markdown", "id": "96199a8f", "metadata": {}, "source": [ "The full A/B, power, peeking, alpha-spending and Bandit examples below are extension material. Bandit reward maximisation and unbiased evaluation of a fixed trading rule answer different questions; classroom familiarity is not an allocation policy." ] }, { "cell_type": "markdown", "id": "356d8e05", "metadata": {}, "source": [ "## Portfolio Losses and Tail Risk" ] }, { "cell_type": "markdown", "id": "b4e4f751", "metadata": {}, "source": [ "### Statistical foundations \u2014 explain before running\n", "\n", "#### Loss quantiles and Expected Shortfall\n", "\n", "Define loss as $L=-R$ for a stated horizon. A large positive value is a large loss.\n", "\n", "$$\\operatorname{VaR}_p=F_L^{-1}(p),\\qquad \\operatorname{ES}_p=\\frac{1}{1-p}\\int_p^1 F_L^{-1}(u)\\,du.$$\n", "\n", "For a continuous loss distribution, ES equals $E[L\\mid L\\ge\\operatorname{VaR}_p]$. VaR describes a cutoff; ES describes severity beyond it. Neither is a maximum possible loss.\n", "\n", "The classroom code uses an interpolated sample quantile and averages observations at/above it. With finite samples or ties this convention can differ from exact fractional-tail quantile integration; report the actual tail count.\n", "\n", "#### Drawdown is a path statistic\n", "\n", "For net portfolio value $V_t$ including the opening value,\n", "\n", "$$H_t=\\max_{s\\le t}V_s,\\qquad D_t=V_t/H_t-1,\\qquad D_{\\min}=\\min_t D_t.$$\n", "\n", "Drawdown measures loss from a running peak. We report its negative signed minimum; positive maximum-drawdown magnitude would be $-D_{\\min}$.\n", "\n", "Reordering identical period returns leaves their empirical VaR unchanged but can change the drawdown path. Risk comparisons must specify horizon, fees, opening NAV and sign convention.\n", "\n", "#### Extreme-value theory: model the tail, state assumptions\n", "\n", "**Block maxima:** under suitable conditions, maxima of large blocks are approximated by a generalised extreme-value (GEV) distribution. Apply to losses; worst single-day loss in a month is not the compounded monthly loss.\n", "\n", "**Peaks over threshold:** sufficiently high exceedances $L-u\\mid L>u$ can be approximated by a generalised Pareto distribution (GPD). Threshold choice trades approximation bias against a small tail sample.\n", "\n", "The tail-shape parameter $\\xi$ distinguishes heavy unbounded tails ($\\xi>0$), an exponential limit ($\\xi=0$) and a bounded upper tail ($\\xi<0$). Dependence may require declustering or adjusted inference.\n", "\n", "Tail fits need stability checks and uncertainty. Worked GEV/POT fitting follows; two tail months cannot validate an extreme-risk model." ] }, { "cell_type": "markdown", "id": "d956a0fa", "metadata": {}, "source": [ "### Measure risk on the reconciled net account\n", "\n", "State the return horizon before taking a quantile. Use the Chapter 2 ledger for monthly net portfolio losses. Daily stock returns, monthly portfolio returns and the worst single day inside a 20-day block are different random variables and require different labels." ] }, { "cell_type": "code", "execution_count": null, "id": "45f20a9d", "metadata": {}, "outputs": [], "source": [ "import numpy as np, pandas as pd\n", "import matplotlib.pyplot as plt\n", "data_url = \"https://statpython.pages.dev/data/\"\n", "ledger = pd.read_csv(data_url + \"chapter2_ledger.csv\",\n", " parse_dates=[\"date\", \"period_end\"])\n", "ledger = ledger.set_index(\"period_end\").sort_index()\n", "assert ledger.index.is_unique and ledger.index.is_monotonic_increasing\n", "assert np.allclose(ledger.nav_end - ledger.nav_start,\n", " ledger.gross_pnl - ledger.cost)\n", "assert np.allclose(ledger.net_return, ledger.nav_end / ledger.nav_start - 1)\n", "net = ledger.net_return\n", "print(len(net), \"realised monthly net returns\")" ] }, { "cell_type": "markdown", "id": "6609c06b", "metadata": {}, "source": [ "Fees are already in net_return; do not subtract them twice. Do not use gross returns for a net-risk report. A distributional risk estimate is evidence for a risk review, not a guarantee that losses stop at that level." ] }, { "cell_type": "markdown", "id": "e98a15c8", "metadata": {}, "source": [ "### Historical VaR and Expected Shortfall\n", "\n", "Define loss as minus the simple net return, so positive numbers mean losses. Historical 95% VaR is the empirical 95th loss percentile. The classroom tail-mean ES averages losses at or above that threshold; record the quantile convention, ties and number of tail observations." ] }, { "cell_type": "code", "execution_count": null, "id": "cc53d6fc", "metadata": {}, "outputs": [], "source": [ "loss = -net\n", "var95 = loss.quantile(0.95)\n", "tail = loss.loc[loss.ge(var95)]\n", "es95 = tail.mean()\n", "print(\"monthly loss VaR95:\", round(float(var95), 4))\n", "print(\"monthly loss ES95:\", round(float(es95), 4))\n", "print(\"tail observations:\", len(tail), \"of\", len(loss))" ] }, { "cell_type": "markdown", "id": "b028be16", "metadata": {}, "source": [ "Monthly loss VaR95 is 0.0823 and the tail-mean ES95 is 0.1345. Only two of 34 months enter that average. The sparse tail cannot support precise extreme-risk claims. An EVT fit adds a model and assumptions, not additional realised crises." ] }, { "cell_type": "markdown", "id": "0c284e9e", "metadata": {}, "source": [ "### Maximum drawdown depends on the path\n", "\n", "Quantiles ignore the order of returns; drawdown preserves it. Include the opening NAV before the first realised outcome so an initial loss is counted. A drawdown is a decline from the running peak, not a one-period loss quantile or a worst-day-in-block statistic." ] }, { "cell_type": "code", "execution_count": null, "id": "1da50037", "metadata": {}, "outputs": [], "source": [ "opening_date = ledger.index[0] - pd.offsets.MonthEnd(1)\n", "opening = pd.Series([ledger.nav_start.iloc[0]], index=[opening_date])\n", "wealth = pd.concat([opening, ledger.nav_end]).sort_index()\n", "drawdown = wealth / wealth.cummax() - 1\n", "print(\"maximum drawdown:\", round(float(drawdown.min()), 4))\n", "fig, ax = plt.subplots(figsize=(9, 4))\n", "drawdown.plot(ax=ax)\n", "ax.set(ylabel=\"Drawdown from running peak\", xlabel=\"Outcome date\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "bdf207a7", "metadata": {}, "source": [ "Maximum realised drawdown is \u22120.1833. This is the teaching sample path, not the worst possible future path. For long-short work, inspect each leg and their joint losses alongside portfolio NAV; gross exposure, concentration, borrow costs and execution also matter." ] }, { "cell_type": "markdown", "id": "c74abc33", "metadata": {}, "source": [ "### Tail co-exceedance and the limits of EVT\n", "\n", "Average correlation does not describe how often two assets fall together in their respective tails. Use the original S&P 500/Tesla daily-return example and keep it separate from the monthly ledger. Conditioning on a narrow market-return slice can reduce Pearson mechanically; report the tail event and its denominator." ] }, { "cell_type": "code", "execution_count": null, "id": "a4cd97b2", "metadata": {}, "outputs": [], "source": [ "daily = pd.read_csv(data_url + \"returns.csv\", parse_dates=[\"Date\"]).set_index(\"Date\").rename(columns={\"sp500\": \"tesla\", \"tesla\": \"sp500\"})\n", "daily = daily[[\"sp500\", \"tesla\"]].dropna()\n", "market_cut = daily.sp500.quantile(0.05)\n", "asset_cut = daily.tesla.quantile(0.05)\n", "market_tail = daily.loc[daily.sp500.le(market_cut)]\n", "co_exceedance = market_tail.tesla.le(asset_cut).mean()\n", "print(\"market-tail days:\", len(market_tail), \"of\", len(daily))\n", "print(\"P(Tesla tail | market tail):\", round(float(co_exceedance), 3))" ] }, { "cell_type": "markdown", "id": "52819b74", "metadata": {}, "source": [ "There are 126 market-tail days and conditional co-exceedance is about 0.325. These full-sample thresholds are descriptive; monitoring thresholds would be estimated using prior data. GEV block minima estimate the worst single daily return within a block, not the compounded return of that block. Tail fitting needs stability checks and uncertainty before any operational use." ] }, { "cell_type": "markdown", "id": "7542a939", "metadata": {}, "source": [ "### Why EVT belongs in a trading risk toolkit\n", "\n", "The centre of a return distribution describes ordinary days. A risk review also asks how severe the rare losses could be, how frequently a threshold might be exceeded and whether joint stress defeats diversification.\n", "\n", "**Extreme-value theory (EVT)** supplies limiting models for extremes and high-threshold excesses. It supports model-based extrapolation beyond observed quantiles, with assumptions and estimation uncertainty.\n", "\n", "The Central Limit Theorem concerns suitably normalised **sums or means**. EVT concerns suitably normalised **maxima or excesses**. Averaging and taking a maximum answer different questions.\n", "\n", "**Check:** A Normal-looking sample mean does not imply Normal one-day loss tails." ] }, { "cell_type": "markdown", "id": "d5d46fd6", "metadata": {}, "source": [ "### Two routes from the same daily losses\n", "\n", "**Two routes:** daily losses \u2192 block maxima \u2192 GEV; daily losses \u2192 high-threshold excesses \u2192 GPD \u2192 unconditional daily risk.\n", "**Read the diagram:** GEV keeps one extreme per block. POT retains all excesses over a threshold. The resulting probability statements have different denominators and horizons." ] }, { "cell_type": "markdown", "id": "d5f157df", "metadata": {}, "source": [ "### GEV: a distribution for block extremes\n", "\n", "For a block maximum loss $M$ and $z=(m-\\mu)/\\sigma$,\n", "\n", "$$G(m)=\\exp\\{-[1+\\xi z]^{-1/\\xi}\\},\\qquad 1+\\xi z>0,\\quad\\sigma>0.$$\n", "\n", "At $\\xi=0$, the limit is $G(m)=\\exp[-\\exp(-z)]$ (Gumbel).\n", "\n", "$\\mu$ is location, $\\sigma$ is scale, and $\\xi$ controls tail shape. Location is generally neither the mean nor the mode.\n", "\n", "**Python convention:** `scipy.stats.genextreme` returns `c = -xi`. A negative `c` indicates a positive heavy-tail index." ] }, { "cell_type": "markdown", "id": "92d9353c", "metadata": {}, "source": [ "### Tail shape controls support and moment existence\n", "\n", "| Shape | Tail family | Interpretation |\n", "|---|---|---|\n", "| $\\xi>0$ | Fr\u00e9chet | Heavy unbounded upper tail |\n", "| $\\xi=0$ | Gumbel | Exponential-type limit |\n", "| $\\xi<0$ | Weibull | Finite fitted upper endpoint |\n", "\n", "For a positive heavy-tail index, the fitted model has a finite moment of order $r$ only when $r\\xi<1$. Thus its variance requires $\\xi<1/2$.\n", "\n", "An estimated index does not prove the population has or lacks a moment. Support, confidence intervals, dependence and stability all matter." ] }, { "cell_type": "markdown", "id": "fba944b5", "metadata": {}, "source": [ "### Predict the tail family\n", "\n", "**Self-check:** `genextreme.fit` returns `c = -0.35`. Which interpretation follows from this fitted family?\n", "\n", "- A bounded upper loss tail because c is negative\n", "- A heavy unbounded upper loss tail because \u03be = 0.35\n", "- A Gaussian tail with standard deviation 0.35\n", "- A 35% probability of losing money\n", "\n", "**Feedback:** SciPy GEV uses c = \u2212\u03be. Therefore c = \u22120.35 means \u03be = +0.35, a heavy unbounded upper tail." ] }, { "cell_type": "markdown", "id": "3ea3b1d8", "metadata": {}, "source": [ "### Build daily losses and non-overlapping blocks\n", "\n", "**Predict:** How many complete 20-day blocks come from 1,496 daily returns?" ] }, { "cell_type": "code", "id": "d9ae3167", "metadata": {}, "source": [ "from scipy.stats import genextreme, genpareto\n", "dow = pd.read_csv(data_url + \"dji.csv\", parse_dates=[\"Date\"])\n", "dow = dow.set_index(\"Date\").sort_index()\n", "evt_loss = -np.log(dow.Close).diff().loc[\"1985\":\"1990\"].dropna()\n", "complete = len(evt_loss) // 20 * 20\n", "block_loss = evt_loss.iloc[:complete].to_numpy().reshape(-1, 20).max(axis=1)\n", "print(\"daily losses:\", len(evt_loss), \"complete blocks:\", len(block_loss))\n", "print(\"unused final days:\", len(evt_loss) - complete)\n", "print(\"largest daily log loss:\", round(float(evt_loss.max()), 4))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "ee38418c", "metadata": {}, "source": [ "There are **74 blocks**, with 16 final days excluded from block fitting. Positive loss means a fall. The largest daily log loss is **0.2563**, corresponding to a simple loss of about 22.6%. A block maximum remains a **single-day** loss." ] }, { "cell_type": "code", "id": "1d67372b", "metadata": {}, "source": [ "c_evt, mu_evt, sigma_evt = genextreme.fit(\n", " block_loss, loc=float(np.median(block_loss)), scale=float(block_loss.std()))\n", "xi_evt = -c_evt\n", "block_q95 = genextreme.ppf(0.95, c_evt, mu_evt, sigma_evt)\n", "assert np.isfinite(genextreme.logpdf(block_loss, c_evt, mu_evt, sigma_evt)).all()\n", "print(\"xi, location, scale:\", np.round([xi_evt, mu_evt, sigma_evt], 4))\n", "print(\"95th percentile of worst daily log loss in a block:\", round(block_q95, 4))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "1770c64a", "metadata": {}, "source": [ "The fit gives $\\xi\\approx0.354$, location 0.0136 and scale 0.0070. Its 95th block-extreme percentile is **0.0506**. This is the fitted worst-day threshold exceeded in 5% of 20-day blocks, rather than a compounded 20-day loss.\n", "\n", "Starting values seed numerical optimisation. Compare alternative starts and likelihoods before accepting a tail fit." ] }, { "cell_type": "markdown", "id": "de222d27", "metadata": {}, "source": [ "### POT: the generalised Pareto model\n", "\n", "Set a high loss threshold $u$. Conditional on $L>u$, define the positive excess $Y=L-u$.\n", "\n", "$$H(y)=1-(1+\\xi y/\\beta)^{-1/\\xi},\\qquad y\\ge0,\\quad1+\\xi y/\\beta>0.$$\n", "\n", "$\\beta>0$ is excess scale. At $\\xi=0$, $H(y)=1-e^{-y/\\beta}$. With $\\xi<0$, excesses end at $-\\beta/\\xi$.\n", "\n", "`genpareto` uses **`c = xi`**, unlike GEV. Fix excess location to zero with `floc=0`. The threshold $u$ belongs outside the excess distribution." ] }, { "cell_type": "markdown", "id": "07be9d80", "metadata": {}, "source": [ "### Explore the tail shape before fitting\n", "\n", "Change the shape in a GPD survival plot: compare \u03be = \u22120.25, 0 and +0.50 with \u03b2 = 0.01. The online slides include an interactive slider.\n", "**Try:** Set \u03be to \u22120.25, 0 and +0.50. Explain the endpoint and the change in the conditional 95th percentile while scale remains fixed." ] }, { "cell_type": "code", "id": "e51cc442", "metadata": {}, "source": [ "u_evt = 0.015\n", "excess = evt_loss.loc[evt_loss.gt(u_evt)] - u_evt\n", "xi_pot, _, beta_pot = genpareto.fit(excess, floc=0)\n", "exceedance_rate = len(excess) / len(evt_loss)\n", "print(\"exceedances:\", len(excess), \"of\", len(evt_loss))\n", "print(\"rate, xi, beta:\", np.round([exceedance_rate, xi_pot, beta_pot], 4))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "aa35d080", "metadata": {}, "source": [ "There are **91 exceedances out of 1,496 days**, a rate of **6.08%**. The fitted $\\xi$ is **0.4363**, with excess scale **0.00643**.\n", "\n", "This is an historical demonstration using the whole 1985\u20131990 window. A prospective risk forecast would fit only on observations already available before the forecast date." ] }, { "cell_type": "markdown", "id": "0ad0efb3", "metadata": {}, "source": [ "### From conditional excesses to unconditional daily risk\n", "\n", "Let $p_u=N_u/n$. For $q>1-p_u$ and $\\xi\\ne0$,\n", "\n", "$$\\operatorname{VaR}_q=u+\\frac{\\beta}{\\xi}\\left[\\left(\\frac{p_u}{1-q}\\right)^\\xi-1\\right].$$\n", "\n", "For the fitted continuous GPD tail and $\\xi<1$,\n", "\n", "$$\\operatorname{ES}_q=\\frac{\\operatorname{VaR}_q+\\beta-\\xi u}{1-\\xi}.$$\n", "\n", "At $\\xi=0$, VaR becomes $u+\\beta\\log[p_u/(1-q)]$, and ES becomes VaR $+\\beta$.\n", "\n", "The exceedance rate reconnects the conditional fit to **all trading days**. When $\\xi\\ge1$, this model has infinite ES, rather than a finite value obtained by ignoring the denominator." ] }, { "cell_type": "code", "id": "a4a617a3", "metadata": {}, "source": [ "q_evt = 0.99\n", "assert q_evt > 1 - exceedance_rate and 0 < xi_pot < 1\n", "var99_evt = u_evt + beta_pot / xi_pot * (\n", " (exceedance_rate / (1 - q_evt)) ** xi_pot - 1)\n", "es99_evt = (var99_evt + beta_pot - xi_pot * u_evt) / (1 - xi_pot)\n", "print(\"daily log-loss VaR99 / ES99:\", np.round([var99_evt, es99_evt], 4))\n", "print(\"VaR as a simple loss:\", round(float(1 - np.exp(-var99_evt)), 4))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "a36424c7", "metadata": {}, "source": [ "Daily **log-loss VaR99 = 0.0327** and **log-loss ES99 = 0.0577**. The VaR corresponds to a simple loss of about **3.21%**.\n", "\n", "Nonlinear transformation of an average differs from averaging transformed losses: `1 - exp(-ES_log)` does **not** equal simple-loss ES. Keep units explicit." ] }, { "cell_type": "markdown", "id": "5f23fcda", "metadata": {}, "source": [ "### Threshold selection trades bias against variance\n", "\n", "**Self-check:** Why should we compare several plausible POT thresholds?\n", "\n", "- The highest threshold always gives the most precise estimate\n", "- More excesses always mean a better tail approximation\n", "- Threshold choice changes approximation bias, sample size and uncertainty\n", "- We should select whichever threshold gives the smallest reported VaR\n", "\n", "**Feedback:** A lower threshold supplies more excesses but may include the body, where the tail approximation is poor. A higher threshold reduces approximation bias only at the cost of fewer observations and greater estimation noise." ] }, { "cell_type": "code", "id": "66807435", "metadata": {}, "source": [ "def tail_fit_at(threshold):\n", " e = evt_loss.loc[evt_loss.gt(threshold)] - threshold\n", " shape, _, scale = genpareto.fit(e, floc=0)\n", " rate = len(e) / len(evt_loss)\n", " quantile = threshold + genpareto.ppf(1 - 0.01 / rate, shape, scale=scale)\n", " return pd.Series({\"count\": len(e), \"xi\": shape, \"daily_var99\": quantile})\n", "thresholds = pd.Series([0.01, 0.015, 0.02], name=\"threshold\")\n", "tail_sensitivity = thresholds.apply(tail_fit_at)\n", "tail_sensitivity.index = thresholds\n", "print(tail_sensitivity.round(4))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "3f1c1810", "metadata": {}, "source": [ "Counts fall **159 \u2192 91 \u2192 45** as the threshold rises. Estimated \u03be changes **0.265 \u2192 0.436 \u2192 0.507**, while VaR99 is about **3.46%, 3.27%, 3.20%** in log-loss units.\n", "\n", "These differences reveal sensitivity. They do not constitute a calibrated confidence interval. Check mean-excess plots, parameter stability and dependence-aware uncertainty before extrapolating." ] }, { "cell_type": "code", "id": "dc82c78c", "metadata": {}, "source": [ "probabilities = (np.arange(1, len(excess) + 1) - 0.5) / len(excess)\n", "theoretical = genpareto.ppf(probabilities, xi_pot, scale=beta_pot)\n", "observed = np.sort(excess.to_numpy())\n", "limit = max(theoretical.max(), observed.max())\n", "fig, ax = plt.subplots(figsize=(8, 4))\n", "ax.scatter(theoretical, observed, s=18)\n", "ax.plot([0, limit], [0, limit], linestyle=\"--\", color=\"black\")\n", "ax.set(xlabel=\"Fitted GPD excess quantile\", ylabel=\"Observed log-loss excess\")\n", "plt.show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "9e686406", "metadata": {}, "source": [ "Compare the upper-right points with the diagonal. An apparently straight centre does not validate the far tail. The 1987 crash strongly influences extrapolation.\n", "\n", "Exceedances cluster during volatile periods. Consider declustering or conditional-volatility models, assess the changed target and denominator, and avoid treating all 91 excesses as independent evidence." ] }, { "cell_type": "markdown", "id": "095f0ed4", "metadata": {}, "source": [ "### Return levels extrapolate a probability, not a calendar promise\n", "\n", "A **1-in-T-day return level** $x_T$ satisfies $P(L>x_T)=1/T$ in the fitted stationary marginal model.\n", "\n", "$$x_T=u+\\frac{\\beta}{\\xi}\\left[(p_uT)^\\xi-1\\right].$$\n", "\n", "For this GPD fit, the 250-, 1,000- and 2,500-day levels are about **4.86%, 8.87% and 13.22% daily log loss**.\n", "\n", "An empirical interpolated quantile remains computable at those probabilities, but it cannot extend beyond the largest observation. The 2,500-day model estimate from 1,496 days is extrapolation.\n", "\n", "**Check:** \u201c1-in-1,000\u201d does not mean the next breach arrives exactly 1,000 days later. Clustering and changing regimes further weaken a simple waiting-time interpretation." ] }, { "cell_type": "markdown", "id": "f3ec6ada", "metadata": {}, "source": [ "### Hill estimation provides a second tail-shape diagnostic\n", "\n", "For descending positive losses $L_{(1)}\\ge\\cdots\\ge L_{(k+1)}>0$,\n", "\n", "$$\\widehat\\xi_k=\\frac1k\\sum_{i=1}^k\\log\\frac{L_{(i)}}{L_{(k+1)}}.$$\n", "\n", "Hill assumes a heavy, regularly varying tail. Plot $\\widehat\\xi_k$ against $k$: a stable region can support a threshold choice, but a plateau is not proof.\n", "\n", "Too few extremes make the estimate noisy. Too many admit the distribution\u2019s body. GEV, GPD and Hill need not agree exactly in a finite dependent sample.\n", "\n", "**Continue in the Book:** Run the Hill and return-level cells, then compare their assumptions and samples with the two fits above." ] }, { "cell_type": "markdown", "id": "91b73adb", "metadata": {}, "source": [ "### Self-check: write a defensible tail-risk statement\n", "\n", "**Self-check:** Which summary correctly reports the POT result?\n", "\n", "- The strategy cannot lose more than 3.27% in a month\n", "- A 1985\u20131990 Dow fit above 1.5% daily log loss uses 91 excesses and estimates daily log-loss VaR99 near 3.27%\n", "- Ninety-one exceedances guarantee a precise future tail probability\n", "- ES99 is a 99% confidence interval for mean return\n", "\n", "**Feedback:** State the daily log-loss target, threshold, number of exceedances and conditional model assumptions. A fitted VaR remains an estimated percentile with uncertainty.\n", "\n", "**Explain in your own words:** What changes if the threshold rises? Which units differ from the monthly portfolio report? What additional evidence would you need for a prospective risk decision?" ] }, { "cell_type": "markdown", "id": "746adaca", "metadata": {}, "source": [ "## Group investigation\n", "\n", "Choose one change: bootstrap block length 3 \u2192 6, feature ret_12_1 \u2192 rvol_21d, or tail percentile 95% \u2192 90%. Predict the result, rerun the affected cells and explain the changed denominator or assumption. Keep the model and evaluation dates fixed; label every additional comparison exploratory.\n", "\n", "Compare two diagnostic claims with another group before final revision." ] }, { "cell_type": "markdown", "id": "539a64ad", "metadata": {}, "source": [ "## Decision memo \u2014 complete during class\n", "\n", "**Question and data contract:** [your answer]\n", "\n", "**IC and coverage:** [your answer]\n", "\n", "**Uncertainty and search history:** [your answer]\n", "\n", "**Risk horizon and tail count:** [your answer]\n", "\n", "**Recommendation, strongest limitation and next test:** [your answer]\n", "\n", "**Peer challenge and revision:** [your answer]\n", "\n", "**Members and individual code/output explanations:** [your answer]" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.12" } }, "nbformat": 4, "nbformat_minor": 5 }