MATH307Su26 - Fat Tails and Black Swans

# Fat Tails and Black Swans > [!abstract] The project in one line > Work the normal and Cauchy probability integrals by hand (one needs a table, the other collapses to an arctangent), verify both with your own trapezoid and Simpson rules, watch every moment of the Cauchy distribution fail to exist while the software hands you confident nonsense, then hold both models against eight years of real stock returns and finish by manufacturing fat tails yourself with a two-dial contamination model whose every target you can compute by hand. ## Introduction Consider [André René Roussimoff](https://en.wikipedia.org/wiki/Andr%C3%A9_the_Giant), better known as André the Giant, who stood 7'4" (88 inches) and weighed 520 pounds. Under a textbook [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution) model of adult male height (mean $69.3$ in, standard deviation $3.0$ in), André sits $6.23$ standard deviations above the mean, an event with probability $2.3\times10^{-10}$ (about 1 in 4 billion). For weight (a normal model with mean $198$ lb, standard deviation $40$ lb), he sits $8.05$ standard deviations out: probability $4.1\times10^{-16}$, about 1 in 2.4 quadrillion. Yet André existed. He was not a statistical impossibility; he had [gigantism](https://en.wikipedia.org/wiki/Gigantism), a rare but real condition (on the order of a few cases per million people). The normal model was catastrophically wrong because human variation has [fat tails](https://en.wikipedia.org/wiki/Fat-tailed_distribution): rare mechanisms create outliers that "should not exist" but regularly do. The same mathematical blindness devastated the financial world on October 19, 1987. Professor [Robert Shiller](https://en.wikipedia.org/wiki/Robert_J._Shiller) was teaching his finance class at Yale when a student with a transistor radio announced that the stock market was crashing during the lecture. By the close, the S&P 500 had fallen $20.47\%$ (the Dow fell $22.6\%$) in a single day, the event now called [Black Monday](https://en.wikipedia.org/wiki/Black_Monday_(1987)). Under the normal-distribution models used across the industry, a one-day drop that size was, for every practical purpose, impossible (Shiller quotes a probability of $3\times10^{-71}$; you will compute the number yourself in Step 5 and see that "impossible" is the right summary however you slice the assumptions). People lost fortunes because they trusted models that systematically ignore the tails where extreme events actually live. If you want the story firsthand, Shiller tells it in his open Yale lecture: [The Problem of Outliers, Crash of 1987 (11:30 to 22:00)](https://youtu.be/XmyZdcvpFgs?si=eRtxm3QquW4b8oIb&t=690). This project explores the mathematics of that model failure with the numerical integration machinery from class. The tools that verify the normal distribution's famous properties to eight digits are the same tools that expose, concretely and reproducibly, why a fat-tailed distribution has no mean for a portfolio optimizer to optimize. The failure is not a software bug; it is an [improper integral](https://en.wikipedia.org/wiki/Improper_integral) that does not converge, and you will hold the divergence in your hands. ![Standard normal and standard Cauchy densities on linear and log scales; near twins in the center, entirely different animals in the tails](Media/fattails_densities.png) ## Project Description The project moves in four passes. First, work both distributions by hand: the Cauchy probability $P(|X|\le1)$ collapses to an exact $\tfrac12$ by antiderivative, the normal counterpart needs the [error function](https://en.wikipedia.org/wiki/Error_function), and the Cauchy mean integral diverges in a way you can compute exactly (moving windows give any answer you like). Second, mirror the hand work in code: build the trapezoid and Simpson rules from scratch, confirm the [68-95-99.7 rule](https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule) and both convergence orders, and reproduce every hand number through the verification helper. Third, run the instrument on reality: the 1987 crash probabilities under both models, then eight years of real DAX daily returns (the shipped dataset) held against both models, plus a generalization that tracks exactly how fast the two tails part ways. Finally, manufacture fat tails yourself: a two-dial contamination model (how often a crisis day strikes, and how hard) whose every diagnostic has an exact hand-computable target, run at three dial settings and then tuned to match the real data. **Foundation task:** Perform the hand calculations in the Mathematical Background (the arctangent probability, the erf probability, the divergent mean with its window dependence, the truncated second moments), then confirm every number computationally: quadrature rules built from scratch, the 68-95-99.7 rule with measured convergence orders, and the Cauchy checks including two machine-precision identities. (This is Steps 1 to 4 below.) **Application task:** Compute the Black Monday probability under the normal and Cauchy models from the classic market statistics and quantify the $10^{60}$-fold disagreement; then load `eustocks_returns.csv` and count the DAX's extreme days against both models' predictions. Generalize by computing the tail probability as a function of threshold and watching the two models' ratio explode. (This is Steps 5 and 6 plus the Generalization section below.) **Key deliverable:** The contamination study at three dial settings (pure, mild, wild): the excess kurtosis and tail table computed BY HAND from the mixture formulas before any code, the formulas verified mechanically by your own Simpson rule at machine precision, simulation held against the exact targets with sampling scatter understood, and a closing dial-tuning exercise that matches the model to the DAX and says what a two-dial world captures and misses. Two companions support the milestone forms: the [[MATH307Su26 - Fat Tails and Black Swans (Milestone Map)|Milestone Map]] says which artifact from this handout answers each form field, and the [[MATH307Su26 - Fat Tails and Black Swans (Verification Table)|Verification Table]] is filled in as you verify and submitted with Milestone 2. ## Mathematical Background: two bells, only one behaves **Idea.** Both densities below are symmetric bells, and near the center they look interchangeable. Everything that matters happens in the tails: the normal density dies *exponentially*, so every probability and every [moment](https://en.wikipedia.org/wiki/Moment_(mathematics)) is a convergent integral; the Cauchy density dies only *polynomially* (like $1/x^2$), so the very first moment already diverges. Numerical integration is our microscope: the same quadrature rules that confirm the normal facts to eight digits will show you, window by window, exactly how the Cauchy integrals fail. ### The normal distribution's mathematical paradise **Recall.** The normal distribution $\mathcal{N}(\mu,\sigma^{2})$ has the [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) $f(x) = \frac{1}{\sigma\sqrt{2\pi}}\, e^{-\frac{(x-\mu)^2}{2\sigma^2}}.$ The exponential factor $e^{-x^{2}/2}$ crushes any polynomial: for every $k$, $x^{k}e^{-x^{2}/2}\to0$ as $|x|\to\infty$, fast enough that all moments exist and are finite, $E[X^k] = \int_{-\infty}^{\infty} x^k\, \frac{1}{\sqrt{2\pi}}\, e^{-x^2/2}\, dx < \infty \quad \text{for all } k.$ This mathematical ideal is what the [68-95-99.7 rule](https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule) summarizes: - 68% of values within 1 standard deviation: $P(-\sigma \leq X \leq \sigma) \approx 0.6827$ - 95% within 2 standard deviations: $P(-2\sigma \leq X \leq 2\sigma) \approx 0.9545$ - 99.7% within 3 standard deviations: $P(-3\sigma \leq X \leq 3\sigma) \approx 0.9973$ There is a catch worth naming: $e^{-x^{2}/2}$ has no elementary antiderivative, so these probabilities are *defined* through a special function, the [error function](https://en.wikipedia.org/wiki/Error_function) $P(|X|\le t)=\int_{-t}^{t}\frac{e^{-x^{2}/2}}{\sqrt{2\pi}}\,dx=\operatorname{erf}\!\Bigl(\frac{t}{\sqrt2}\Bigr),$ which is itself computed numerically (that is what the table in a statistics book is). Quadrature is not a toy here; it is how these numbers exist at all. For the warm-up, $\operatorname{erf}(1/\sqrt2)=0.682689$ and $\operatorname{erf}(3/\sqrt2)=0.997300$ to six digits. ### The Cauchy distribution's mathematical nightmare The [Cauchy distribution](https://en.wikipedia.org/wiki/Cauchy_distribution) appears deceptively similar, $f(x; x_0, \gamma) = \frac{1}{\pi}\, \frac{\gamma}{(x-x_0)^2 + \gamma^2},$ with center $x_0$ and scale $\gamma$ (the standard case is $x_0=0$, $\gamma=1$). Unlike the normal density it *does* have an elementary antiderivative, which is why we can work it exactly by hand: $\int \frac{dx}{\pi(1+x^{2})} = \frac{\arctan x}{\pi} + C \qquad\Longrightarrow\qquad \boxed{\,P(|X|\le t) = \frac{2}{\pi}\arctan t.\,}$ **Worked probability (run to the checkable number).** At $t=1$: $P(|X|\le1)=\tfrac{2}{\pi}\arctan 1=\tfrac{2}{\pi}\cdot\tfrac{\pi}{4}=\tfrac12$, *exactly*. Compare the two bells at matched scale points: the normal packs $68.27\%$ of its probability within one scale unit; the Cauchy packs only $50\%$. At $t=3$ the gap is starker: $\tfrac{2}{\pi}\arctan 3 = 0.795167$ versus $0.997300$. Read it from the tails: beyond three scale units the Cauchy still holds $20.48\%$ of its probability where the normal holds $0.27\%$, a factor of $75.9$. These four numbers ($\tfrac12$ exactly, $0.682689$, $0.795167$ vs $0.997300$, ratio $75.9$) are what your warm-up code must reproduce. **The mean does not exist (and you can compute exactly how).** The polynomial decay $f(x)\sim \tfrac{1}{\pi x^{2}}$ is fatal to the moments. Formally the mean would be $E[X] = \int_{-\infty}^{\infty} x \cdot \frac{1}{\pi(1+x^2)}\, dx,$ and this integrand also has an elementary antiderivative, $\int \frac{x\,dx}{\pi(1+x^{2})}=\frac{\ln(1+x^{2})}{2\pi}+C,$ which grows like $\tfrac{\ln|x|}{\pi}$ without bound. So the integral is not [absolutely convergent](https://en.wikipedia.org/wiki/Absolute_convergence), and the value of a truncated version depends on *how* you open the window: $\int_{-R}^{R} \frac{x\,dx}{\pi(1+x^{2})} = 0 \quad\text{(odd integrand, symmetric window)},$ $\boxed{\,\int_{-R}^{2R} \frac{x\,dx}{\pi(1+x^{2})} = \frac{1}{2\pi}\ln\!\frac{1+4R^{2}}{1+R^{2}} \;\xrightarrow[R\to\infty]{}\; \frac{\ln 2}{\pi} = 0.220636.\,}$ Same distribution, same "$\to\inftyquot;, two different limits (and a window $[-R, cR]$ can be tuned to hit *any* value $\tfrac{\ln c}{\pi}$). A quantity whose value depends on the bookkeeping of the limit is not a number; the mean is undefined. The symmetric-window $0$ is called the [Cauchy principal value](https://en.wikipedia.org/wiki/Cauchy_principal_value), and mistaking it for "the mean" is exactly the trap your Step 3 code will spring on purpose. **The second moment diverges in slow motion.** For the variance we would need $E[X^{2}]$. Truncate and compute exactly (divide $\tfrac{x^2}{1+x^2}=1-\tfrac{1}{1+x^2}$): $\int_{-R}^{R} \frac{x^{2}\,dx}{\pi(1+x^{2})} = \frac{2}{\pi}\bigl(R-\arctan R\bigr) \;\approx\; \frac{2R}{\pi} \quad\text{for large } R,$ growing *linearly* with the window, forever. The normal counterpart, by an integration by parts (use $\varphi'(x)=-x\varphi(x)$ for the standard normal density $\varphi$), is $\int_{-R}^{R} x^{2}\varphi(x)\,dx = \operatorname{erf}\!\Bigl(\frac{R}{\sqrt2}\Bigr) - 2R\,\varphi(R) \;\xrightarrow[R\to\infty]{}\; 1,$ already $0.999985$ at $R=5$. Both closed forms are hand results your Step 4 code checks by quadrature, window by window. > [!note] Where this comes from > Nothing here is new machinery: the antiderivatives ($\arctan$, $\ln$, parts) are Calc 2, the trapezoid and Simpson rules with their $O(h^{2})$/$O(h^{4})$ error orders are the Taylor-series analysis from class, and the probability content is supplied. The project is the assembly, plus one hard lesson: a numerical answer is only as meaningful as the convergence of the integral behind it. ### Worked example: the 1987 crash prediction Consider daily stock returns modeled with the classic historical statistics $\mu = 0.05\%$ and $\sigma = 1.2\%$ (these are the Su25 anchor numbers; the real dataset in Step 6 will land remarkably close to them). A $20\%$ single-day drop sits at $z = \frac{-0.20 - 0.0005}{0.012} = -16.7083$ standard deviations. The normal model's probability of a day at least that bad is $P(X \leq -20\%) = \int_{-\infty}^{-0.20} \frac{1}{0.012\sqrt{2\pi}}\, e^{-\frac{(x-0.0005)^2}{2(0.012)^2}}\, dx = \tfrac12\operatorname{erfc}\!\Bigl(\tfrac{16.7083}{\sqrt{2}}\Bigr) \approx 5.7\times10^{-63},$ (here $\operatorname{erfc}(x)=1-\operatorname{erf}(x)$ is the [complementary error function](https://en.wikipedia.org/wiki/Error_function#Complementary_error_function), the standard way to compute deep normal tails without cancellation), i.e., one such day per $7.0\times10^{59}$ years of trading. The age of the universe is about $1.4\times10^{10}$ years. Shiller quotes $3\times10^{-71}$ for this event (slightly different assumed statistics); at this depth in the tail the eighth decimal of $\sigma$ moves the answer by decades of orders of magnitude, and every version of the number means the same thing: *the model calls the event impossible*. Reality check: it happened, on October 19, 1987, while the models were running. The same $z$-score under a Cauchy model with the same center and scale ($x_0=0.0005$, $\gamma=0.012$) is an arctangent away: $P(X\le-20\%) = \frac{1}{\pi}\arctan\!\frac{\gamma}{0.20+x_0} = 0.019028,$ about once per $53$ trading days. That is absurd in the *other* direction (roughly five crashes a year), which is its own lesson: $\gamma$ is not a standard deviation, and matching "scale" across families is a modeling choice with enormous consequences. The two models, fed the same two parameters, disagree about Black Monday by a factor of $3.3\times10^{60}$. Model choice is not a detail; in the tails, model choice is everything. ## Implementation Guidelines The workflow is the same in every language, so we state it once and then implement it step by step. > [!abstract] The workflow (state it once, run it in any language) > Given the two densities, the market statistics $(\mu,\sigma)=(0.0005, 0.012)$, and the returns file: > 1. **Warm up**: define the verification helper, then reproduce the hand numbers ($\tfrac{2}{\pi}\arctan 1=\tfrac12$ exactly; $\operatorname{erf}(1/\sqrt2)=0.682689$; within-3 values $0.795167$ vs $0.997300$, tail ratio $75.9$; crash $z=-16.7083$). > 2. **Build the rules**: trapezoid $\;h\bigl(\tfrac{f_0}{2}+f_1+\cdots+f_{n-1}+\tfrac{f_n}{2}\bigr)$ and Simpson $\;\tfrac{h}{3}\bigl(f_0+4f_{\text{odd}}+2f_{\text{even}}+f_n\bigr)$; verify the 68-95-99.7 rule against the erf values and measure the orders (error ratios $4$ and $16$ per doubling of $n$). > 3. **Spring the Cauchy trap**: quadrature reproduces $P(|X|\le1)=\tfrac12$ to machine precision; a symmetric window makes $\int x f\,dx$ cancel to machine zero; windows $[-R,2R]$ converge to $\tfrac{\ln 2}{\pi}=0.220636$ instead. The "mean" depends on the window, so it does not exist. > 4. **Watch the moments**: truncated $E[X^{2}]$ on $[-R,R]$ matches $\operatorname{erf}(R/\sqrt2)-2R\varphi(R)\to1$ (normal) and $\tfrac{2}{\pi}(R-\arctan R)\approx\tfrac{2R}{\pi}$ (Cauchy, diverging); then ask your library integrator for the Cauchy moments and record what it does. > 5. **Price Black Monday**: $P(X\le-20\%)$ under both models; normal $5.7\times10^{-63}$ versus Cauchy $0.019028$; quantify the $3.3\times10^{60}$ disagreement. > 6. **Meet reality**: load `eustocks_returns.csv`, summarize the DAX column, and count days beyond $t\sigma$ for $t=2,3,4,5$ against both models' predicted counts. > 7. **Generalize**: tail probability versus threshold, $\operatorname{erfc}(t/\sqrt2)$ versus $\tfrac{2}{\pi}\arctan\tfrac1t$, and their exploding ratio (its own section below). > 8. **Your own tails**: three contamination-dial settings, hand targets first (the finale). Pick your language below; each callout contains Steps 1 to 6 in order, run top to bottom. Code is four-language and portable (Python is the run-verified reference; Mathematica is numerical only, i.e., no symbolic `Integrate`). Download the shipped dataset from the course website (link on the Assignments page) and put `eustocks_returns.csv` next to your script; it is a plain CSV with a header row (`day,DAX,SMI,CAC,FTSE`) and 1859 rows of daily [log returns](https://en.wikipedia.org/wiki/Rate_of_return#Logarithmic_or_continuously_compounded_return) derived from the [EuStockMarkets](https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/EuStockMarkets.html) prices (four major European indices, 1991 to 1998). > [!example]- Python (reference) > > **Step 1: define the helper, then warm up on the hand numbers.** The helper compares two values (or arrays) and reports the largest difference; every check in this project runs through it. > ```python > import numpy as np > import matplotlib.pyplot as plt > from math import erf, erfc, atan, sqrt, pi, log > def same(name, X, Y, tol=1e-9): # verification helper > d = np.max(np.abs(np.asarray(X) - np.asarray(Y))) > print(f"{name}: max difference {d:.1e}", "OK" if d < tol else "CHECK THIS") > normal_pdf = lambda x: np.exp(-x**2/2)/np.sqrt(2*np.pi) > cauchy_pdf = lambda x: 1/(np.pi*(1 + x**2)) > c_within = lambda t: 2*np.arctan(t)/np.pi # Cauchy P(|X| <= t), exact > n_within = lambda t: erf(t/sqrt(2)) # normal P(|X| <= t), exact > same("Cauchy P(|X| <= 1) equals 1/2 exactly", c_within(1.0), 0.5) > print("normal P(|X| <= 1) =", round(n_within(1), 6)) # 0.682689 > print("within 3:", round(c_within(3), 6), round(n_within(3), 6)) # 0.795167 0.997300 > print("beyond-3 ratio:", round((1 - c_within(3))/(1 - n_within(3)), 1)) # 75.9 > mu, sd = 0.0005, 0.012 # the 1987 market statistics > print("crash z-score:", round((-0.20 - mu)/sd, 4)) # -16.7083 > ``` > *Look for:* the exact-$\tfrac12$ helper check at $0.0$ (this identity is exact, not approximate), $0.682689$, $0.795167$ vs $0.997300$, ratio $75.9$, and $z=-16.7083$. > > **Step 2: build the quadrature rules and verify the 68-95-99.7 rule.** Both rules from scratch; the erf values are the ground truth, and doubling $n$ must cut the errors by $4$ (trapezoid) and $16$ (Simpson). > ```python > def my_trap(f, a, b, n): > x = np.linspace(a, b, n+1); y = f(x); h = (b - a)/n > return h*(np.sum(y) - 0.5*(y[0] + y[-1])) > def my_simp(f, a, b, n): # n must be even > if n % 2: raise ValueError("n must be even for Simpson's rule") > x = np.linspace(a, b, n+1); y = f(x); h = (b - a)/n > return h/3*(y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]) + y[-1]) > for t in (1, 2, 3): > exact = n_within(t) > pt, ps = my_trap(normal_pdf, -t, t, 100), my_simp(normal_pdf, -t, t, 100) > print(f"[-{t},{t}] exact {exact:.6f} trap err {abs(pt-exact):.1e} simp err {abs(ps-exact):.1e}") > errs = [(abs(my_trap(normal_pdf, -1, 1, n) - n_within(1)), > abs(my_simp(normal_pdf, -1, 1, n) - n_within(1))) for n in (4, 8, 16, 32, 64)] > for i in range(1, 5): > print(f"doubling ratio: trap {errs[i-1][0]/errs[i][0]:.2f} simp {errs[i-1][1]/errs[i][1]:.2f}") > ``` > *Look for:* at $n=100$ on $[-1,1]$, trapezoid error $1.6\times10^{-5}$ and Simpson error $8.6\times10^{-10}$; doubling ratios settling to $4.00$ and $16.02$ (orders $2$ and $4$, the Taylor analysis made visible). > > **Step 3: spring the Cauchy trap.** The rules still work on any finite window (the integrand is smooth); what fails is the *concept* of the infinite-window value. The symmetric window is the machine-precision exactness check: an odd integrand on a symmetric grid cancels term by term. > ```python > print("Simpson Cauchy [-1,1] err:", abs(my_simp(cauchy_pdf, -1, 1, 2000) - 0.5)) > mean_integrand = lambda x: x*cauchy_pdf(x) > sym = my_trap(mean_integrand, -100, 100, 10000) > same("odd integrand on a symmetric window cancels", sym, 0.0, tol=1e-14) > for R in (10.0, 100.0, 1000.0, 10000.0): # now open the window asymmetrically > q = my_simp(mean_integrand, -R, 2*R, 400000) > exact = log((1 + 4*R**2)/(1 + R**2))/(2*pi) > print(f"R={R:8.0f} window integral {q:.6f} hand value {exact:.6f}") > print("limit ln(2)/pi =", round(log(2)/pi, 6)) # 0.220636 > ``` > *Look for:* Simpson error $0.0$ at $n=2000$ (machine precision) and the symmetric cancellation at $\sim10^{-16}$; then the asymmetric windows marching to $0.220636$, matching the hand formula every step. Two different "answers" for the same integral is the whole point: the mean does not exist. > > **Step 4: watch the moments on growing windows, then ask the library.** The truncated second moments have exact hand values; check quadrature against them, then hand the divergent integrals to `scipy.integrate.quad` and record what comes back. > ```python > from scipy import integrate > for R in (1.0, 2.0, 3.0, 5.0, 10.0): # normal: converges to 1 > q = my_simp(lambda x: x**2*normal_pdf(x), -R, R, 20000) > print(f"normal R={R:4.1f} {q:.6f} exact {erf(R/sqrt(2)) - 2*R*float(normal_pdf(R)):.6f}") > for R in (5.0, 50.0, 500.0, 5000.0): # Cauchy: grows like 2R/pi forever > q = my_simp(lambda x: x**2*cauchy_pdf(x), -R, R, 400000) > print(f"Cauchy R={R:6.0f} {q:10.4f} exact {(2/pi)*(R - atan(R)):10.4f} 2R/pi {2*R/pi:10.4f}") > print(integrate.quad(cauchy_pdf, -np.inf, np.inf)) # total probability: fine > print(integrate.quad(mean_integrand, -np.inf, np.inf)) # the "mean": silent nonsense > print(integrate.quad(mean_integrand, 0, np.inf)) # half-line: warning + garbage > print(integrate.quad(lambda x: x**2*cauchy_pdf(x), -np.inf, np.inf)) # negative "moment"! > ``` > *Look for:* normal $0.198748, 0.738536, 0.970709, 0.999985, 1.000000$; Cauchy $2.3088, 30.8437, 317.3112, 3182.0990$ tracking $2R/\pi$. Then the library: total probability $1.00000000$ (fine); the mean comes back as $0.0$ with error estimate $0.0$ and *no warning at all* (it silently computed the principal value); the half-line mean returns $13.05$ with error estimate $2.6$ and a max-subdivisions warning; the second moment returns $-1.63662$, a *negative* value for a positive integrand, with a "probably divergent" warning. Confident, wrong, and only sometimes flagged: record all four printouts verbatim. > > **Step 5: price Black Monday under both models.** The normal tail through `erfc` (this is how the "impossible" number is actually computed), a Simpson check on the $[-100\%, -20\%]$ band, and the Cauchy tail through the arctangent. > ```python > p_norm = 0.5*erfc(((0.20 + mu)/sd)/sqrt(2)) # P(X <= -20%), exact > print(f"normal P = {p_norm:.3e} -> once per {1/p_norm/252:.2e} years") > norm_market = lambda x: np.exp(-((x - mu)/sd)**2/2)/(sd*np.sqrt(2*np.pi)) > print(f"Simpson on [-1, -0.2]: {my_simp(norm_market, -1.0, -0.20, 20000):.3e}") > gam = 0.012 # Cauchy, same center and scale > p_cau = atan(gam/(0.20 + mu))/pi > print(f"Cauchy P = {p_cau:.6f} -> once per {1/p_cau:.1f} trading days") > cauchy_market = lambda x: (gam/np.pi)/((x - mu)**2 + gam**2) > band = (atan((-0.20 - mu)/gam) - atan((-1.0 - mu)/gam))/pi > same("Cauchy band [-100%, -20%] quadrature vs arctan", > my_simp(cauchy_market, -1.0, -0.20, 20000), band) > print(f"model disagreement: {p_cau/p_norm:.1e}") > ``` > *Look for:* normal $5.699\times10^{-63}$ (once per $7.0\times10^{59}$ years) with Simpson agreeing on the band; Cauchy $0.019028$ (once per $52.6$ trading days) with the band check at $\sim10^{-17}$; disagreement $3.3\times10^{60}$. > > **Step 6: meet reality (the shipped returns).** Load the CSV, summarize the DAX column, fit the normal by moments and the Cauchy by the middle half (median and half the [interquartile range](https://en.wikipedia.org/wiki/Interquartile_range); you *cannot* fit a Cauchy by moments, it has none), then count extreme days against both models. > ```python > data = np.genfromtxt("eustocks_returns.csv", delimiter=",", names=True) > r = data["DAX"]; n = r.size > m, s = r.mean(), r.std(ddof=1) > med = np.median(r); q1, q3 = np.percentile(r, [25, 75]); gam_hat = (q3 - q1)/2 > print(f"n={n} mean {100*m:.4f}% sd {100*s:.4f}% min {100*r.min():.3f}% max {100*r.max():.3f}%") > m2, m4 = np.mean((r - m)**2), np.mean((r - m)**4) > print(f"median {100*med:.4f}% IQR/2 {100*gam_hat:.4f}% excess kurtosis {m4/m2**2 - 3:.2f}") > z_worst = (r.min() - m)/s > print(f"worst day z = {z_worst:.2f}, normal P = {0.5*erfc(-z_worst/sqrt(2)):.2e}") > for t in (2, 3, 4, 5): > obs = int(np.sum(np.abs(r - m) > t*s)) > e_n = n*erfc(t/sqrt(2)); e_c = n*(2/pi)*atan(gam_hat/(t*s)) > print(f"t={t}: observed {obs:3d} normal {e_n:7.2f} Cauchy {e_c:7.2f}") > fw = np.mean(np.abs(r - m) <= 3*s); fb = np.mean(np.abs(r - m) > 3*s) > same("within + beyond fractions sum to 1", fw + fb, 1.0, tol=1e-15) > bins = np.linspace(r.min(), r.max(), 80); grid = np.linspace(r.min(), r.max(), 800) > plt.hist(r, bins=bins, density=True, alpha=0.45, label="DAX daily log returns") > plt.plot(grid, np.exp(-((grid-m)/s)**2/2)/(s*np.sqrt(2*np.pi)), label="normal (mean, sd)") > plt.plot(grid, (gam_hat/np.pi)/((grid-med)**2 + gam_hat**2), label="Cauchy (median, IQR/2)") > plt.yscale("log"); plt.ylim(1e-3, 200); plt.xlabel("daily log return"); plt.legend(); plt.show() > ``` > *Look for:* $n=1859$, mean $0.0652\%$/day, sd $1.0301\%$/day (eerily close to the Su25 anchor $0.05\%$/$1.2\%$), min $-9.628\%$, excess [kurtosis](https://en.wikipedia.org/wiki/Kurtosis) $6.28$ (normal would be $0$); worst day $z=-9.41$ with normal $P=2.5\times10^{-21}$ (the model expected $5\times10^{-18}$ such days in the sample; the sample contains one); the count table $90/24/6/2$ observed versus $84.6/5.0/0.12/0.001$ (normal) and $309.8/209.2/157.6/126.4$ (Cauchy); the fraction check exact at $0$. The normal nails $t=2$ and dies in the tails; the Cauchy overshoots everywhere. Reality lives between the two models. > [!example]- MATLAB > > **Step 1: define the helper, then warm up on the hand numbers.** The helper reports the largest difference; every check runs through it. (Careful with `sprintf`: avoid stray backslashes; `\n` is the only escape used here.) > ```matlab > same = @(name, X, Y) fprintf('%s: max difference %.1e\n', name, max(abs(X - Y), [], 'all')); > normal_pdf = @(x) exp(-x.^2/2)/sqrt(2*pi); > cauchy_pdf = @(x) 1./(pi*(1 + x.^2)); > c_within = @(t) 2*atan(t)/pi; % Cauchy P(|X| <= t), exact > n_within = @(t) erf(t/sqrt(2)); % normal P(|X| <= t), exact > same('Cauchy P(|X| <= 1) equals 1/2 exactly', c_within(1), 0.5) > fprintf('normal P(|X| <= 1) = %.6f\n', n_within(1)) % 0.682689 > fprintf('within 3: %.6f %.6f\n', c_within(3), n_within(3)) % 0.795167 0.997300 > fprintf('beyond-3 ratio: %.1f\n', (1 - c_within(3))/(1 - n_within(3))) % 75.9 > mu = 0.0005; sd = 0.012; % the 1987 market statistics > fprintf('crash z-score: %.4f\n', (-0.20 - mu)/sd) % -16.7083 > ``` > *Look for:* the exact-$\tfrac12$ check at $0.0$, then $0.682689$, $0.795167$ vs $0.997300$, ratio $75.9$, $z=-16.7083$. > > **Step 2: build the quadrature rules and verify the 68-95-99.7 rule.** Written as weight vectors dotted with the samples ($[\tfrac12,1,\dots,1,\tfrac12]$ and $\tfrac13[1,4,2,\dots,2,4,1]$), which keeps both rules anonymous and script-safe (MATLAB scripts only allow `function` definitions at the very end of the file). > ```matlab > my_trap = @(f, a, b, n) ((b-a)/n) * (f(linspace(a, b, n+1)) * [0.5, ones(1, n-1), 0.5].'); > my_simp = @(f, a, b, n) ((b-a)/n)/3 * ... % n must be even > (f(linspace(a, b, n+1)) * [1, 3 - (-1).^(1:n-1), 1].'); > for t = 1:3 > exact = n_within(t); > fprintf('[-%d,%d] exact %.6f trap err %.1e simp err %.1e\n', t, t, exact, ... > abs(my_trap(normal_pdf, -t, t, 100) - exact), abs(my_simp(normal_pdf, -t, t, 100) - exact)); > end > ns = [4 8 16 32 64]; et = zeros(1,5); es = et; > for i = 1:5 > et(i) = abs(my_trap(normal_pdf, -1, 1, ns(i)) - n_within(1)); > es(i) = abs(my_simp(normal_pdf, -1, 1, ns(i)) - n_within(1)); > end > fprintf('doubling ratios: trap %.2f simp %.2f\n', et(4)/et(5), es(4)/es(5)) > ``` > *Look for:* trapezoid error $1.6\times10^{-5}$ and Simpson error $8.6\times10^{-10}$ at $n=100$ on $[-1,1]$; final doubling ratios $4.00$ and $16.02$. > > **Step 3: spring the Cauchy trap.** > ```matlab > fprintf('Simpson Cauchy [-1,1] err: %.1e\n', abs(my_simp(cauchy_pdf, -1, 1, 2000) - 0.5)) > mean_integrand = @(x) x.*cauchy_pdf(x); > sym = my_trap(mean_integrand, -100, 100, 10000); > same('odd integrand on a symmetric window cancels', sym, 0) > for R = [10 100 1000 10000] % open the window asymmetrically > q = my_simp(mean_integrand, -R, 2*R, 400000); > fprintf('R=%6d window integral %.6f hand value %.6f\n', R, q, log((1+4*R^2)/(1+R^2))/(2*pi)); > end > fprintf('limit ln(2)/pi = %.6f\n', log(2)/pi) % 0.220636 > ``` > *Look for:* Simpson error at machine level, the symmetric cancellation at $\sim10^{-16}$, windows marching to $0.220636$. > > **Step 4: watch the moments on growing windows, then ask the library.** `integral` is MATLAB's adaptive integrator; give it the divergent integrals and record everything it prints. > ```matlab > for R = [1 2 3 5 10] % normal: converges to 1 > q = my_simp(@(x) x.^2.*normal_pdf(x), -R, R, 20000); > fprintf('normal R=%4.1f %.6f exact %.6f\n', R, q, erf(R/sqrt(2)) - 2*R*normal_pdf(R)); > end > for R = [5 50 500 5000] % Cauchy: grows like 2R/pi forever > q = my_simp(@(x) x.^2.*cauchy_pdf(x), -R, R, 400000); > fprintf('Cauchy R=%6d %10.4f exact %10.4f 2R/pi %10.4f\n', R, q, (2/pi)*(R - atan(R)), 2*R/pi); > end > warning('on', 'all') > total = integral(cauchy_pdf, -Inf, Inf) % total probability: fine > try, m1 = integral(mean_integrand, -Inf, Inf), catch ME, disp(ME.message), end > try, mh = integral(mean_integrand, 0, Inf), catch ME, disp(ME.message), end > try, m2 = integral(@(x) x.^2.*cauchy_pdf(x), -Inf, Inf), catch ME, disp(ME.message), end > ``` > *Look for:* the same window numbers as the Python reference ($0.198748 \to 1.000000$; $2.3088, 30.8437, 317.3112, 3182.0990$). MATLAB's `integral` behaves differently from SciPy on the divergent cases (expect warnings about maximum intervals or non-convergence, possibly a symmetric-cancellation $0$); the *values and messages are library-specific and the meaninglessness is the point*. Record verbatim what your version prints. > > **Step 5: price Black Monday under both models.** > ```matlab > p_norm = 0.5*erfc(((0.20 + mu)/sd)/sqrt(2)); > fprintf('normal P = %.3e -> once per %.2e years\n', p_norm, 1/p_norm/252) > norm_market = @(x) exp(-((x - mu)/sd).^2/2)/(sd*sqrt(2*pi)); > fprintf('Simpson on [-1, -0.2]: %.3e\n', my_simp(norm_market, -1.0, -0.20, 20000)) > gam = 0.012; > p_cau = atan(gam/(0.20 + mu))/pi; > fprintf('Cauchy P = %.6f -> once per %.1f trading days\n', p_cau, 1/p_cau) > cauchy_market = @(x) (gam/pi)./((x - mu).^2 + gam^2); > band = (atan((-0.20 - mu)/gam) - atan((-1.0 - mu)/gam))/pi; > same('Cauchy band [-100%, -20%] quadrature vs arctan', ... > my_simp(cauchy_market, -1.0, -0.20, 20000), band) > fprintf('model disagreement: %.1e\n', p_cau/p_norm) > ``` > *Look for:* $5.699\times10^{-63}$, $0.019028$, band check at $\sim10^{-17}$, disagreement $3.3\times10^{60}$. > > **Step 6: meet reality (the shipped returns).** `readtable` handles the header; the quantiles are computed by interpolation on the sorted sample (the type-7 convention, matching numpy and R defaults), so all four languages report the same $\hat\gamma$. > ```matlab > T = readtable('eustocks_returns.csv'); r = T.DAX; n = numel(r); > m = mean(r); s = std(r); med = median(r); > qs = sort(r); qtl = @(p) interp1(0:n-1, qs, p*(n-1)); % type-7 quantile > gam_hat = (qtl(0.75) - qtl(0.25))/2; > fprintf('n=%d mean %.4f%% sd %.4f%% min %.3f%% max %.3f%%\n', n, 100*m, 100*s, 100*min(r), 100*max(r)) > m2 = mean((r - m).^2); m4 = mean((r - m).^4); > fprintf('median %.4f%% IQR/2 %.4f%% excess kurtosis %.2f\n', 100*med, 100*gam_hat, m4/m2^2 - 3) > z_worst = (min(r) - m)/s; > fprintf('worst day z = %.2f, normal P = %.2e\n', z_worst, 0.5*erfc(-z_worst/sqrt(2))) > for t = 2:5 > obs = sum(abs(r - m) > t*s); > fprintf('t=%d: observed %3d normal %7.2f Cauchy %7.2f\n', t, obs, ... > n*erfc(t/sqrt(2)), n*(2/pi)*atan(gam_hat/(t*s))); > end > fw = mean(abs(r - m) <= 3*s); fb = mean(abs(r - m) > 3*s); > same('within + beyond fractions sum to 1', fw + fb, 1) > edges = linspace(min(r), max(r), 80); grid = linspace(min(r), max(r), 800); > histogram(r, edges, 'Normalization', 'pdf', 'FaceAlpha', 0.45); hold on > plot(grid, exp(-((grid-m)/s).^2/2)/(s*sqrt(2*pi))) > plot(grid, (gam_hat/pi)./((grid-med).^2 + gam_hat^2)) > set(gca, 'YScale', 'log'); ylim([1e-3 200]); xlabel('daily log return') > legend('DAX daily log returns', 'normal (mean, sd)', 'Cauchy (median, IQR/2)') > ``` > *Look for:* the same summary numbers as the Python reference ($n=1859$, mean $0.0652\%$, sd $1.0301\%$, kurtosis $6.28$, counts $90/24/6/2$); the fraction check exact at $0$. > [!example]- R > > **Step 1: define the helper, then warm up on the hand numbers.** The helper reports the largest difference; every check runs through it. Base R has no `erf`; use the identity $\operatorname{erf}(x)=2\Phi(x\sqrt2)-1$ with `pnorm` (that is what the code below does). > ```r > same <- function(name, X, Y) cat(sprintf("%s: max difference %.1e\n", name, max(abs(X - Y)))) > normal_pdf <- function(x) exp(-x^2/2)/sqrt(2*pi) > cauchy_pdf <- function(x) 1/(pi*(1 + x^2)) > c_within <- function(t) 2*atan(t)/pi # Cauchy P(|X| <= t), exact > n_within <- function(t) 2*pnorm(t) - 1 # normal P(|X| <= t) = erf(t/sqrt(2)) > same("Cauchy P(|X| <= 1) equals 1/2 exactly", c_within(1), 0.5) > cat(sprintf("normal P(|X| <= 1) = %.6f\n", n_within(1))) # 0.682689 > cat(sprintf("within 3: %.6f %.6f\n", c_within(3), n_within(3))) # 0.795167 0.997300 > cat(sprintf("beyond-3 ratio: %.1f\n", (1 - c_within(3))/(1 - n_within(3)))) # 75.9 > mu <- 0.0005; sd_ <- 0.012 # the 1987 market statistics > cat(sprintf("crash z-score: %.4f\n", (-0.20 - mu)/sd_)) # -16.7083 > ``` > *Look for:* the exact-$\tfrac12$ check at $0.0$, then $0.682689$, $0.795167$ vs $0.997300$, ratio $75.9$, $z=-16.7083$. (`sd` is a base function name; the underscore keeps it intact.) > > **Step 2: build the quadrature rules and verify the 68-95-99.7 rule.** 1-based odd/even index slices, unlike Python. > ```r > my_trap <- function(f, a, b, n) { > x <- seq(a, b, length.out = n+1); y <- f(x); h <- (b - a)/n > h*(sum(y) - 0.5*(y[1] + y[n+1])) > } > my_simp <- function(f, a, b, n) { # n must be even > x <- seq(a, b, length.out = n+1); y <- f(x); h <- (b - a)/n > h/3*(y[1] + 4*sum(y[seq(2, n, 2)]) + 2*sum(y[seq(3, n-1, 2)]) + y[n+1]) > } > for (t in 1:3) { > exact <- n_within(t) > cat(sprintf("[-%d,%d] exact %.6f trap err %.1e simp err %.1e\n", t, t, exact, > abs(my_trap(normal_pdf, -t, t, 100) - exact), abs(my_simp(normal_pdf, -t, t, 100) - exact))) > } > ns <- c(4, 8, 16, 32, 64) > et <- sapply(ns, function(n) abs(my_trap(normal_pdf, -1, 1, n) - n_within(1))) > es <- sapply(ns, function(n) abs(my_simp(normal_pdf, -1, 1, n) - n_within(1))) > cat(sprintf("doubling ratios: trap %.2f simp %.2f\n", et[4]/et[5], es[4]/es[5])) > ``` > *Look for:* trapezoid error $1.6\times10^{-5}$ and Simpson error $8.6\times10^{-10}$ at $n=100$ on $[-1,1]$; final doubling ratios $4.00$ and $16.02$. > > **Step 3: spring the Cauchy trap.** > ```r > cat(sprintf("Simpson Cauchy [-1,1] err: %.1e\n", abs(my_simp(cauchy_pdf, -1, 1, 2000) - 0.5))) > mean_integrand <- function(x) x*cauchy_pdf(x) > sym <- my_trap(mean_integrand, -100, 100, 10000) > same("odd integrand on a symmetric window cancels", sym, 0) > for (R in c(10, 100, 1000, 10000)) { # open the window asymmetrically > q <- my_simp(mean_integrand, -R, 2*R, 400000) > cat(sprintf("R=%6d window integral %.6f hand value %.6f\n", R, q, log((1+4*R^2)/(1+R^2))/(2*pi))) > } > cat(sprintf("limit ln(2)/pi = %.6f\n", log(2)/pi)) # 0.220636 > ``` > *Look for:* Simpson error at machine level, the symmetric cancellation at $\sim10^{-16}$, windows marching to $0.220636$. > > **Step 4: watch the moments on growing windows, then ask the library.** `integrate` is R's adaptive integrator; wrap the divergent calls in `tryCatch` because R sometimes *stops with an error* instead of warning. > ```r > for (R in c(1, 2, 3, 5, 10)) { # normal: converges to 1 > q <- my_simp(function(x) x^2*normal_pdf(x), -R, R, 20000) > cat(sprintf("normal R=%4.1f %.6f exact %.6f\n", R, q, (2*pnorm(R) - 1) - 2*R*normal_pdf(R))) > } > for (R in c(5, 50, 500, 5000)) { # Cauchy: grows like 2R/pi forever > q <- my_simp(function(x) x^2*cauchy_pdf(x), -R, R, 400000) > cat(sprintf("Cauchy R=%6d %10.4f exact %10.4f 2R/pi %10.4f\n", R, q, (2/pi)*(R - atan(R)), 2*R/pi)) > } > print(integrate(cauchy_pdf, -Inf, Inf)) # total probability: fine > tryCatch(print(integrate(mean_integrand, -Inf, Inf)), error = function(e) print(e$message)) > tryCatch(print(integrate(mean_integrand, 0, Inf)), error = function(e) print(e$message)) > tryCatch(print(integrate(function(x) x^2*cauchy_pdf(x), -Inf, Inf)), error = function(e) print(e$message)) > ``` > *Look for:* the same window numbers as the Python reference. On the divergent calls, R typically stops with "the integral is probably divergent" (an error, not a value); SciPy silently returned $0.0$ for the same mean. Same mathematics, four different failure modes; record verbatim what your version says. > > **Step 5: price Black Monday under both models.** The normal tail is `pnorm(z)` directly (R's own erf-style special function). > ```r > z <- (-0.20 - mu)/sd_ > p_norm <- pnorm(z) > cat(sprintf("normal P = %.3e -> once per %.2e years\n", p_norm, 1/p_norm/252)) > norm_market <- function(x) exp(-((x - mu)/sd_)^2/2)/(sd_*sqrt(2*pi)) > cat(sprintf("Simpson on [-1, -0.2]: %.3e\n", my_simp(norm_market, -1.0, -0.20, 20000))) > gam <- 0.012 > p_cau <- atan(gam/(0.20 + mu))/pi > cat(sprintf("Cauchy P = %.6f -> once per %.1f trading days\n", p_cau, 1/p_cau)) > cauchy_market <- function(x) (gam/pi)/((x - mu)^2 + gam^2) > band <- (atan((-0.20 - mu)/gam) - atan((-1.0 - mu)/gam))/pi > same("Cauchy band [-100%, -20%] quadrature vs arctan", > my_simp(cauchy_market, -1.0, -0.20, 20000), band) > cat(sprintf("model disagreement: %.1e\n", p_cau/p_norm)) > ``` > *Look for:* $5.699\times10^{-63}$, $0.019028$, band check at $\sim10^{-17}$, disagreement $3.3\times10^{60}$. > > **Step 6: meet reality (the shipped returns).** `quantile`'s default (type 7) matches numpy, so $\hat\gamma$ agrees across languages. If you save the figure under headless `Rscript`, wrap the plotting in a `draw()` closure and use `png(...); draw(); dev.off()` (`dev.copy` alone writes nothing headless). > ```r > d <- read.csv("eustocks_returns.csv"); r <- d$DAX; n <- length(r) > m <- mean(r); s <- sd(r); med <- median(r) > qs <- quantile(r, c(0.25, 0.75)); gam_hat <- (qs[2] - qs[1])/2 > cat(sprintf("n=%d mean %.4f%% sd %.4f%% min %.3f%% max %.3f%%\n", > n, 100*m, 100*s, 100*min(r), 100*max(r))) > m2 <- mean((r - m)^2); m4 <- mean((r - m)^4) > cat(sprintf("median %.4f%% IQR/2 %.4f%% excess kurtosis %.2f\n", 100*med, 100*gam_hat, m4/m2^2 - 3)) > z_worst <- (min(r) - m)/s > cat(sprintf("worst day z = %.2f, normal P = %.2e\n", z_worst, pnorm(z_worst))) > for (t in 2:5) { > obs <- sum(abs(r - m) > t*s) > cat(sprintf("t=%d: observed %3d normal %7.2f Cauchy %7.2f\n", t, obs, > n*2*pnorm(-t), n*(2/pi)*atan(gam_hat/(t*s)))) > } > fw <- mean(abs(r - m) <= 3*s); fb <- mean(abs(r - m) > 3*s) > same("within + beyond fractions sum to 1", fw + fb, 1) > grid <- seq(min(r), max(r), length.out = 800) > hist(r, breaks = 80, freq = FALSE, col = rgb(0.4, 0.6, 0.8, 0.45), xlab = "daily log return", main = "") > lines(grid, exp(-((grid - m)/s)^2/2)/(s*sqrt(2*pi)), col = "darkorange") > lines(grid, (gam_hat/pi)/((grid - med)^2 + gam_hat^2), col = "forestgreen") > legend("topleft", c("normal (mean, sd)", "Cauchy (median, IQR/2)"), > col = c("darkorange", "forestgreen"), lty = 1) > ``` > *Look for:* the same summary numbers as the Python reference ($n=1859$, mean $0.0652\%$, sd $1.0301\%$, kurtosis $6.28$, counts $90/24/6/2$); the fraction check exact at $0$. (Base `hist` draws a linear axis; the log-scale tail view is in the shipped script if you want it.) > [!example]- Mathematica (numerical only) > > **Step 1: define the helper, then warm up on the hand numbers.** The helper reports the largest difference; every check runs through it. (`E`, `N`, `D`, `C` are protected names; use your own. Keep everything in machine numbers.) > ```wolfram > same[name_, x_, y_] := Print[name, ": max difference ", Max@Abs[x - y]]; > normalPdf[x_] := Exp[-x^2/2]/Sqrt[2 Pi]; > cauchyPdf[x_] := 1/(Pi (1 + x^2)); > cWithin[t_] := 2 ArcTan[t]/Pi; (* Cauchy P(|X| <= t), exact *) > nWithin[t_] := Erf[t/Sqrt[2.]]; (* normal P(|X| <= t), exact *) > same["Cauchy P(|X| <= 1) equals 1/2 exactly", cWithin[1.], 0.5] > Print["normal P(|X| <= 1) = ", nWithin[1.]] (* 0.682689 *) > Print["within 3: ", cWithin[3.], " ", nWithin[3.]] (* 0.795167 0.997300 *) > Print["beyond-3 ratio: ", (1 - cWithin[3.])/(1 - nWithin[3.])] (* 75.9 *) > mu = 0.0005; sd = 0.012; (* the 1987 market statistics *) > Print["crash z-score: ", (-0.20 - mu)/sd] (* -16.7083 *) > ``` > *Look for:* the exact-$\tfrac12$ check at $0.$, then $0.682689$, $0.795167$ vs $0.997300$, ratio $75.9$, $z=-16.7083$. > > **Step 2: build the quadrature rules and verify the 68-95-99.7 rule.** `Subdivide` builds the grid; the arithmetic is listable, so the whole grid evaluates at once. > ```wolfram > myTrap[f_, a_, b_, n_] := Module[{y = f[Subdivide[N[a], N[b], n]], h = (b - a)/n}, > h (Total[y] - 0.5 (First[y] + Last[y]))]; > mySimp[f_, a_, b_, n_] := Module[{y = f[Subdivide[N[a], N[b], n]], h = (b - a)/n}, > h/3 (First[y] + 4 Total[y[[2 ;; -2 ;; 2]]] + 2 Total[y[[3 ;; -3 ;; 2]]] + Last[y])]; > Do[Module[{exact = nWithin[N[t]]}, > Print["[-", t, ",", t, "] exact ", exact, > " trap err ", Abs[myTrap[normalPdf, -t, t, 100] - exact], > " simp err ", Abs[mySimp[normalPdf, -t, t, 100] - exact]]], {t, 1, 3}] > ns = {4, 8, 16, 32, 64}; > et = Table[Abs[myTrap[normalPdf, -1, 1, n] - nWithin[1.]], {n, ns}]; > es = Table[Abs[mySimp[normalPdf, -1, 1, n] - nWithin[1.]], {n, ns}]; > Print["doubling ratios: trap ", et[[4]]/et[[5]], " simp ", es[[4]]/es[[5]]] > ``` > *Look for:* trapezoid error $1.6\times10^{-5}$ and Simpson error $8.6\times10^{-10}$ at $n=100$ on $[-1,1]$; final doubling ratios $4.00$ and $16.02$. > > **Step 3: spring the Cauchy trap.** > ```wolfram > Print["Simpson Cauchy [-1,1] err: ", Abs[mySimp[cauchyPdf, -1, 1, 2000] - 0.5]] > meanIntegrand[x_] := x cauchyPdf[x]; > sym = myTrap[meanIntegrand, -100, 100, 10000]; > same["odd integrand on a symmetric window cancels", sym, 0.] > Do[Module[{q = mySimp[meanIntegrand, -rr, 2 rr, 400000]}, > Print["R=", rr, " window integral ", q, " hand value ", > Log[(1 + 4 rr^2)/(1 + rr^2)]/(2 Pi) // N]], {rr, {10., 100., 1000., 10000.}}] > Print["limit ln(2)/pi = ", Log[2.]/Pi] (* 0.220636 *) > ``` > *Look for:* Simpson error at machine level, the symmetric cancellation at $\sim10^{-16}$, windows marching to $0.220636$. > > **Step 4: watch the moments on growing windows, then ask the library.** `NIntegrate` is the numerical integrator (no symbolic `Integrate` anywhere); the divergent calls will emit `ncvb`/`slwcon` messages, which are the data here. > ```wolfram > Do[Module[{q = mySimp[Function[x, x^2 normalPdf[x]], -rr, rr, 20000]}, > Print["normal R=", rr, " ", q, " exact ", Erf[rr/Sqrt[2.]] - 2 rr normalPdf[rr]]], > {rr, {1., 2., 3., 5., 10.}}] > Do[Module[{q = mySimp[Function[x, x^2 cauchyPdf[x]], -rr, rr, 400000]}, > Print["Cauchy R=", rr, " ", q, " exact ", (2/Pi) (rr - ArcTan[rr]) // N, > " 2R/pi ", 2 rr/Pi // N]], {rr, {5., 50., 500., 5000.}}] > Print[NIntegrate[cauchyPdf[x], {x, -Infinity, Infinity}]] (* total probability: fine *) > Print[NIntegrate[meanIntegrand[x], {x, -Infinity, Infinity}]] (* watch the messages *) > Print[NIntegrate[meanIntegrand[x], {x, 0, Infinity}]] > Print[NIntegrate[x^2 cauchyPdf[x], {x, -Infinity, Infinity}]] > ``` > *Look for:* the same window numbers as the Python reference. On the divergent calls Mathematica typically prints `NIntegrate::ncvb` or `NIntegrate::slwcon` and returns some finite number anyway; SciPy silently returned $0.0$ for the mean and $-1.64$ (negative!) for the second moment. Record verbatim what your version prints. > > **Step 5: price Black Monday under both models.** > ```wolfram > pNorm = 0.5 Erfc[((0.20 + mu)/sd)/Sqrt[2.]]; > Print["normal P = ", pNorm, " -> once per ", 1/pNorm/252., " years"] > normMarket[x_] := Exp[-((x - mu)/sd)^2/2]/(sd Sqrt[2 Pi]); > Print["Simpson on [-1, -0.2]: ", mySimp[normMarket, -1., -0.20, 20000]] > gam = 0.012; > pCau = ArcTan[gam/(0.20 + mu)]/Pi; > Print["Cauchy P = ", pCau, " -> once per ", 1/pCau, " trading days"] > cauchyMarket[x_] := (gam/Pi)/((x - mu)^2 + gam^2); > band = (ArcTan[(-0.20 - mu)/gam] - ArcTan[(-1.0 - mu)/gam])/Pi; > same["Cauchy band [-100%, -20%] quadrature vs arctan", > mySimp[cauchyMarket, -1., -0.20, 20000], band] > Print["model disagreement: ", pCau/pNorm] > ``` > *Look for:* $5.699\times10^{-63}$, $0.019028$, band check at $\sim10^{-17}$, disagreement $3.3\times10^{60}$. > > **Step 6: meet reality (the shipped returns).** `Import` reads the CSV (header row first); `Quantile` with the parameter list below is the type-7 convention, matching the other three languages. > ```wolfram > raw = Import["eustocks_returns.csv", "CSV"]; > r = N[Rest[raw][[All, 2]]]; n = Length[r]; (* column 2 = DAX *) > m = Mean[r]; s = StandardDeviation[r]; med = Median[r]; > {q1, q3} = Quantile[r, {1/4, 3/4}, {{1, -1}, {0, 1}}]; gamHat = (q3 - q1)/2; > Print["n=", n, " mean ", 100 m, "% sd ", 100 s, "% min ", 100 Min[r], "% max ", 100 Max[r], "%"] > m2 = Mean[(r - m)^2]; m4 = Mean[(r - m)^4]; > Print["median ", 100 med, "% IQR/2 ", 100 gamHat, "% excess kurtosis ", m4/m2^2 - 3] > zWorst = (Min[r] - m)/s; > Print["worst day z = ", zWorst, ", normal P = ", 0.5 Erfc[-zWorst/Sqrt[2.]]] > Do[Module[{obs = Count[Abs[r - m], u_ /; u > t s]}, > Print["t=", t, ": observed ", obs, " normal ", n Erfc[t/Sqrt[2.]], > " Cauchy ", n (2/Pi) ArcTan[gamHat/(t s)]]], {t, 2, 5}] > fw = Count[Abs[r - m], u_ /; u <= 3 s]/n; fb = Count[Abs[r - m], u_ /; u > 3 s]/n; > same["within + beyond fractions sum to 1", N[fw + fb], 1.] > grid = Subdivide[Min[r], Max[r], 800]; > fitN = Exp[-((grid - m)/s)^2/2]/(s Sqrt[2 Pi]); > fitC = (gamHat/Pi)/((grid - med)^2 + gamHat^2); > Show[Histogram[r, 80, "PDF", ChartStyle -> Opacity[0.45]], > ListLinePlot[{Transpose[{grid, fitN}], Transpose[{grid, fitC}]}, > PlotLegends -> {"normal (mean, sd)", "Cauchy (median, IQR/2)"}], > AxesLabel -> {"daily log return", "density"}] > ``` > *Look for:* the same summary numbers as the Python reference ($n=1859$, mean $0.0652\%$, sd $1.0301\%$, kurtosis $6.28$, counts $90/24/6/2$); the fraction check exact at $0$. What the guided run should produce: ![DAX daily log returns (1991 to 1998) with the moment-fitted normal and the middle-half-fitted Cauchy, on linear and log density scales: the center looks normal, the tails do not](Media/fattails_returns.png) ![Truncated second moment versus window half-width: the normal curve locks onto 1 by R = 5 while the Cauchy curve grows like 2R/pi without bound](Media/fattails_moments.png) > [!warning] Verify against ground truth (required) > Record these checked numbers; "it ran" is not verification. > - Warm-up: `same("Cauchy P(|X| <= 1) equals 1/2 exactly", ...)` at $0.0$ (exact, the arctangent identity); normal $0.682689$; within-3 values $0.795167$ vs $0.997300$; beyond-3 ratio $75.9$; crash $z=-16.7083$. > - Rules: at $n=100$ on $[-1,1]$, trapezoid error $1.6\times10^{-5}$, Simpson error $8.6\times10^{-10}$; per-doubling error ratios settle at $4.00$ and $16.02$ (orders $2.00$ and $4.00$). > - Machine-precision exactness checks: Simpson on the Cauchy over $[-1,1]$ at $n=2000$ reproduces $\tfrac12$ with error $0.0$; the symmetric-window mean integrand cancels to $\sim3\times10^{-16}$; the within-plus-beyond fraction check in Step 6 is exactly $0$. > - The moving window: $[-R,2R]$ integrals $0.219449,\ 0.220624,\ 0.220635,\ 0.220636$ at $R=10,10^2,10^3,10^4$, matching the hand value $\tfrac{1}{2\pi}\ln\tfrac{1+4R^2}{1+R^2}$ each time and converging to $\tfrac{\ln2}{\pi}=0.220636$ (not $0$: the "mean" moved). > - Truncated second moments: normal $0.198748,\ 0.738536,\ 0.970709,\ 0.999985,\ 1.000000$ at $R=1,2,3,5,10$; Cauchy $2.3088,\ 30.8437,\ 317.3112,\ 3182.0990$ at $R=5,50,500,5000$ (tracking $2R/\pi$, never converging). > - Library integrator (SciPy 1.17 reference run): Cauchy mean on $(-\infty,\infty)$ returned $0.0$ with error estimate $0.0$ and *no warning*; half-line mean returned $13.05\pm2.6$ with a max-subdivisions warning; second moment returned $-1.63662$ (negative, for a positive integrand) with a "probably divergent" warning. Other languages fail differently; record yours verbatim. > - Black Monday: normal $P=5.699\times10^{-63}$ (once per $6.96\times10^{59}$ years); Cauchy $P=0.019028$ (once per $52.6$ trading days); band check through the helper at $\sim3.6\times10^{-17}$; model disagreement $3.3\times10^{60}$. > - DAX ($n=1859$): mean $0.0652\%$/day, sd $1.0301\%$/day, min $-9.628\%$, max $5.076\%$, median $0.0473\%$, IQR/2 $0.5520\%$, excess kurtosis $6.28$; worst-day $z=-9.41$ (normal $P=2.48\times10^{-21}$); counts beyond $t\cdot$sd: observed $90, 24, 6, 2$ versus normal $84.58, 5.02, 0.12, 0.001$ and Cauchy $309.84, 209.21, 157.62, 126.37$ at $t=2,3,4,5$. > [!note] Fat tails cluster (look at your histogram's dates) > The DAX's worst day in the sample ($-9.63\%$, business-time position 35, mid-August 1991) is the [Soviet coup attempt](https://en.wikipedia.org/wiki/1991_Soviet_coup_attempt); its best day ($+5.08\%$) comes *two days later* when the coup collapsed; the second-worst day ($-6.01\%$, position 1651, late October 1997) is the [Asian financial crisis](https://en.wikipedia.org/wiki/1997_Asian_financial_crisis) mini-crash. Extreme days arrive in bunches around events, violating the independence that moment-based models quietly assume. You will quantify tail weight here; clustering is a second, separate failure mode worth naming in your write-up. ### Generalization: tail probability versus threshold (Step 7) The DAX table compared the models at four thresholds. Now sweep the threshold and watch the disagreement grow without bound. For the standard versions of the two distributions, both two-sided tails are exact special functions (no quadrature needed, which is itself the lesson: by this point you have *verified* these functions with your own rules): $P_{\mathcal N}(|X|>t)=\operatorname{erfc}\!\Bigl(\tfrac{t}{\sqrt2}\Bigr), \qquad P_{C}(|X|>t)=1-\tfrac{2}{\pi}\arctan t=\tfrac{2}{\pi}\arctan\tfrac1t,$ where the last equality is the reflection identity $\arctan t+\arctan\tfrac1t=\tfrac{\pi}{2}$ (check it at machine precision; it is this section's built-in exactness check). For large $t$ the asymptotics explain everything you will see: $\operatorname{erfc}(t/\sqrt2)\sim\sqrt{\tfrac{2}{\pi}}\,\tfrac{e^{-t^{2}/2}}{t}$ dies exponentially while $\tfrac{2}{\pi}\arctan\tfrac1t\sim\tfrac{2}{\pi t}$ dies only like $1/t$, so the ratio grows like $e^{t^{2}/2}$: *faster than exponentially*. Every extra standard deviation of "extreme" multiplies the models' disagreement by a bigger factor than the last one. *Prediction (before running):* from the $t=3$ ratio of $75.9$ in the warm-up, estimate the $t=6$ ratio, then check how far off you were. > [!example]- Step 7 code, Python (reference) > > **Step 7: sweep the threshold.** Two exact tails, their ratio, and the reflection identity through the helper. > ```python > same("arctan(t) + arctan(1/t) equals pi/2", atan(3.0) + atan(1/3.0), pi/2, tol=1e-15) > for t in range(1, 7): > tn, tc = erfc(t/sqrt(2)), (2/pi)*atan(1.0/t) > print(f"t={t} normal {tn:.3e} Cauchy {tc:.4f} ratio {tc/tn:.2e}") > ts = np.linspace(0.5, 6, 400) > tn = np.array([erfc(t/sqrt(2)) for t in ts]); tc = (2/np.pi)*np.arctan(1/ts) > fig, ax = plt.subplots(1, 2, figsize=(11, 4)) > ax[0].semilogy(ts, tn, label="normal"); ax[0].semilogy(ts, tc, label="Cauchy") > ax[0].set_xlabel("threshold t (sd units)"); ax[0].set_ylabel("two-sided tail"); ax[0].legend() > ax[1].semilogy(ts, tc/tn); ax[1].set_xlabel("threshold t (sd units)") > ax[1].set_ylabel("Cauchy tail / normal tail"); plt.tight_layout(); plt.show() > ``` > *Look for:* the identity check at $0.0$; ratios $1.58,\ 6.49,\ 75.9,\ 2.46\times10^{3},\ 2.19\times10^{5},\ 5.33\times10^{7}$ at $t=1..6$; on the log axis the normal tail is a parabola opening downward while the Cauchy tail is nearly flat. > [!example]- Step 7 code, MATLAB > > **Step 7: sweep the threshold.** > ```matlab > same('arctan(t) + arctan(1/t) equals pi/2', atan(3) + atan(1/3), pi/2) > for t = 1:6 > tn = erfc(t/sqrt(2)); tc = (2/pi)*atan(1/t); > fprintf('t=%d normal %.3e Cauchy %.4f ratio %.2e\n', t, tn, tc, tc/tn); > end > ts = linspace(0.5, 6, 400); > tn = erfc(ts/sqrt(2)); tc = (2/pi)*atan(1./ts); > subplot(1,2,1); semilogy(ts, tn, ts, tc); xlabel('threshold t (sd units)') > ylabel('two-sided tail'); legend('normal', 'Cauchy') > subplot(1,2,2); semilogy(ts, tc./tn); xlabel('threshold t (sd units)') > ylabel('Cauchy tail / normal tail') > ``` > *Look for:* the same six ratios as the Python reference; the identity check at machine level. > [!example]- Step 7 code, R > > **Step 7: sweep the threshold.** The normal two-sided tail is `2*pnorm(-t)` (equal to $\operatorname{erfc}(t/\sqrt2)$). > ```r > same("arctan(t) + arctan(1/t) equals pi/2", atan(3) + atan(1/3), pi/2) > for (t in 1:6) { > tn <- 2*pnorm(-t); tc <- (2/pi)*atan(1/t) > cat(sprintf("t=%d normal %.3e Cauchy %.4f ratio %.2e\n", t, tn, tc, tc/tn)) > } > ts <- seq(0.5, 6, length.out = 400) > tn <- 2*pnorm(-ts); tc <- (2/pi)*atan(1/ts) > par(mfrow = c(1, 2)) > matplot(ts, cbind(tn, tc), type = "l", log = "y", lty = 1, col = c("navy", "darkorange"), > xlab = "threshold t (sd units)", ylab = "two-sided tail") > legend("bottomleft", c("normal", "Cauchy"), col = c("navy", "darkorange"), lty = 1) > plot(ts, tc/tn, type = "l", log = "y", xlab = "threshold t (sd units)", > ylab = "Cauchy tail / normal tail") > par(mfrow = c(1, 1)) > ``` > *Look for:* the same six ratios as the Python reference; the identity check at machine level. > [!example]- Step 7 code, Mathematica (numerical only) > > **Step 7: sweep the threshold.** > ```wolfram > same["arctan(t) + arctan(1/t) equals pi/2", ArcTan[3.] + ArcTan[1/3.], N[Pi/2]] > Do[Module[{tn = Erfc[t/Sqrt[2.]], tc = (2/Pi) ArcTan[1./t]}, > Print["t=", t, " normal ", tn, " Cauchy ", tc, " ratio ", tc/tn]], {t, 1, 6}] > LogPlot[{Erfc[t/Sqrt[2.]], (2/Pi) ArcTan[1/t]}, {t, 0.5, 6}, > AxesLabel -> {"threshold t (sd units)", "two-sided tail"}, > PlotLegends -> {"normal", "Cauchy"}] > LogPlot[((2/Pi) ArcTan[1/t])/Erfc[t/Sqrt[2.]], {t, 0.5, 6}, > AxesLabel -> {"threshold t (sd units)", "Cauchy tail / normal tail"}] > ``` > *Look for:* the same six ratios as the Python reference; the identity check at machine level. What the generalization should produce: ![Two-sided tail probability versus threshold on a log axis (exponential versus polynomial decay), and the Cauchy-to-normal tail ratio exploding past 10^7 by t = 6](Media/fattails_tails.png) > [!warning] Verify the generalization (required) > - The reflection identity $\arctan t+\arctan\tfrac1t=\tfrac{\pi}{2}$ through the helper at $0.0$ (machine precision, the built-in exactness check). > - Normal two-sided tails: $3.173\times10^{-1},\ 4.550\times10^{-2},\ 2.700\times10^{-3},\ 6.334\times10^{-5},\ 5.733\times10^{-7},\ 1.973\times10^{-9}$ at $t=1,\dots,6$. > - Cauchy two-sided tails: $0.5000,\ 0.2952,\ 0.2048,\ 0.1560,\ 0.1257,\ 0.1051$ at $t=1,\dots,6$. > - Ratios: $1.58,\ 6.49,\ 75.9,\ 2.46\times10^{3},\ 2.19\times10^{5},\ 5.33\times10^{7}$; each unit of $t$ multiplies the disagreement by more than the last (the $e^{t^2/2}$ signature). ### Your tails: three contamination dials Real market tails do not come from exotic textbook densities; they come from *regime mixing*: most days are calm, a few days are crises. The finale has you manufacture fat tails yourself with the simplest mixing model there is, the contaminated normal (a two-component [mixture distribution](https://en.wikipedia.org/wiki/Mixture_distribution)): with probability $1-p$ the day is calm and the return is drawn from $N(0,\sigma^2)$; with probability $p$ it is a crisis day drawn from $N(0,(k\sigma)^2)$. Two dials: how often ($p$) and how hard ($k$). The payoff of building your object instead of choosing one is that every diagnostic has an exact target you can reach by hand. Conditioning on which component fired (a weighted average of two normal moments, nothing more), with $\sigma=1$: $\mathbb{E}[X^2] = 1-p+pk^2, \qquad \mathbb{E}[X^4] = 3\bigl(1-p+pk^4\bigr),$ so the [excess kurtosis](https://en.wikipedia.org/wiki/Kurtosis) is $\boxed{\;\kappa(p,k)=\frac{3\,(1-p+pk^4)}{(1-p+pk^2)^2}-3,\;}$ and the tail probability is a mixture of two normal tails your erf machinery already computes, with $\sigma_{\text{mix}}=\sqrt{1-p+pk^2}$: $P\bigl(|X|>t\,\sigma_{\text{mix}}\bigr)=(1-p)\cdot 2\bar{\Phi}(t\,\sigma_{\text{mix}})+p\cdot 2\bar{\Phi}(t\,\sigma_{\text{mix}}/k).$ Run the study at three dial settings (the personalities): 1. **Setting A, pure ($p=0$).** The plain normal, as control: $\kappa=0$, and every tail count must match Step 6's normal predictions. 2. **Setting B, mild ($p=0.05$, $k=3$).** One day in twenty lives in a 3-sigma regime. Exact $\kappa=4.6531$; per $1859$ days (the DAX sample size) the mixture predicts $71.7,\ 22.7,\ 10.7,\ 4.5$ exceedances at $t=2,3,4,5$. Compare the middle of that table to the DAX's observed $90/24/6/2$: two humble dials already put you in reality's neighborhood. 3. **Setting C, wild ($p=0.01$, $k=10$).** Rare 10-sigma catastrophes. Exact $\kappa=73.5$, and a paradox to test: *fewer* 2-sigma days than the pure normal ($23.3$ versus $84.6$ per $1859$) but thousands of times more 5-sigma days ($8.9$ versus $0.001$). Fat tails are not "more volatility everywhere"; they are quiet stretches punctured by catastrophes. For each setting: (i) compute $\kappa$ and the tail table by hand from the formulas, before any code; (ii) verify the two moment formulas mechanically by running your own Simpson rule on $x^2 f(x)$ and $x^4 f(x)$ over a wide window; this check is deterministic, so the helper must report agreement at $\sim10^{-10}$, your machine-precision anchor for the finale; (iii) simulate $N=10{,}000$ days (one uniform decides calm-or-crisis, one normal draw supplies the return; no new samplers) and hold the empirical tail counts and sample kurtosis against your exact targets, expecting sampling scatter rather than machine agreement. Watch the wild setting especially: its sample kurtosis lands well below the exact $73.5$ at $N=10{,}000$, and you should explain why (the kurtosis estimator itself lives in the tail it is trying to measure). Close by tuning the dial to reality: find a $(p,k)$ whose $\kappa$ matches the DAX's $6.28$; there is a whole curve of solutions, and $p=5\%$, $k=3.3$ is one of them ($\kappa=6.24$). Check your tuned setting's predicted tail table against the observed $90/24/6/2$ and close with a verdict: which features of eight years of returns does a two-dial world capture, and which does it miss? #### Student task loop for implementation, analysis and reflection 1. **Predict** $\kappa$ and the tail table for each setting by hand from the dials, before running anything. 2. **Implement** the two-route check (exact formula versus your Simpson moments) and the simulation. 3. **Compare** the three layers honestly: formula versus quadrature at machine precision, simulation versus exact targets within sampling scatter, tuned dial versus the DAX. 4. **Interpret** what $p$ and $k$ mean for anyone averaging, optimizing, or pricing risk: how often the crisis regime strikes, how hard it hits, and why the calm stretches in between are the trap. ## Reflection Framework Address these (2-3 focused questions per category, no more): ### Convergence and breakdown - Your Simpson rule integrated the Cauchy density over $[-1,1]$ to machine precision, yet the Cauchy mean does not exist. Reconcile those two facts: what exactly diverges, and where does the quadrature rule stop being the issue? - The library integrators gave four different behaviors on the same divergent integrals (silent zero, warned garbage, negative "moment", or a hard error). Which failure mode is the most dangerous in production risk software, and why? ### Models and markets - The DAX table shows the normal model nearly exact at $t=2$, off by a factor of $50$ at $t=4$, and off by a factor of nearly $2000$ at $t=5$, while the Cauchy overshoots at every threshold. What does "the right model" even mean when one is calibrated to the center and the other to the tails? - Institutions kept normal-based risk models long after 1987 (and after 1998, and 2008). Given what you now know about which computations *require* finite moments, what would a bank actually have to give up to abandon them? ### Mathematical insights - **[Cauchy principal value](https://en.wikipedia.org/wiki/Cauchy_principal_value) versus absolute convergence:** your symmetric windows gave $0$ at machine precision and your $[-R,2R]$ windows gave $\tfrac{\ln2}{\pi}$, both "correct" limits of truncations. Why does probability theory demand absolute convergence for a mean, rather than accepting the principal value? - **The [central limit theorem](https://en.wikipedia.org/wiki/Central_limit_theorem)'s fine print:** the CLT needs a finite variance. In fact, the average of $N$ independent Cauchy draws is Cauchy with the *same* scale (averaging buys nothing, ever). What does that do to the idea of "estimating the mean return from more data," and how did your Object C's running mean display it? - **[Kurtosis](https://en.wikipedia.org/wiki/Kurtosis) as a diagnostic with a blind spot:** excess kurtosis measures tail weight through the fourth moment, but for $\nu\le4$ that moment does not exist, so the estimate never settles. When is a diverging diagnostic itself the most informative reading? ## (Optional) Mathematical Extensions These go beyond the project's stated level; attempt them if interested. - **The [Student's t](https://en.wikipedia.org/wiki/Student%27s_t-distribution) bridge:** the t family interpolates between your two poles ($\nu=1$ is Cauchy, $\nu\to\infty$ is normal), with moments existing up to order $\nu-1$. Redo the moment-window study across $\nu$ and watch moments switch on one at a time. - **[Stable distributions](https://en.wikipedia.org/wiki/Stable_distribution) and the generalized CLT:** sums of infinite-variance variables converge to stable laws (the Cauchy is the $\alpha=1$ case), which is why fat tails are not washed out by aggregation. - **[Extreme value theory](https://en.wikipedia.org/wiki/Extreme_value_theory):** estimate a tail index from your Object C's sample (the [Hill estimator](https://en.wikipedia.org/wiki/Heavy-tailed_distribution#Hill's_tail-index_estimator) is a dozen lines) and compare it to the true decay exponent. - **Prove the averaging claim:** using characteristic functions, show the mean of $N$ standard Cauchy draws is standard Cauchy, then verify it by simulation. ## (Optional) Real-World Context ### Applications - **[Value at Risk](https://en.wikipedia.org/wiki/Value_at_risk):** the industry-standard risk number is a tail quantile; your threshold sweep is exactly the computation whose model-dependence VaR critics point to. - **[Long-Term Capital Management](https://en.wikipedia.org/wiki/Long-Term_Capital_Management):** a Nobel-laureate-staffed fund lost $4.6 billion in 1998 when correlated fat-tail moves arrived that its models priced as many-sigma impossibilities. - **[Catastrophe modeling](https://en.wikipedia.org/wiki/Catastrophe_modeling) and reinsurance:** hurricane and earthquake losses are the textbook fat-tail business; extreme value theory (not moment fitting) is the working tool. - **[Black swan theory](https://en.wikipedia.org/wiki/Black_swan_theory):** Taleb's popularization of exactly the phenomenon you computed: the events that matter most live where moment-based models see nothing. ### Technical challenges - **[Tail-index estimation](https://en.wikipedia.org/wiki/Heavy-tailed_distribution):** the tail exponent controls everything, but by definition it must be estimated from the rare observations, so its own error bars are wide. - **[Volatility clustering](https://en.wikipedia.org/wiki/Volatility_clustering):** your DAX extremes arrived in bunches (August 1991, October 1997), violating the independence assumptions that both of your models share. - **[Numerical integration of heavy tails](https://en.wikipedia.org/wiki/Numerical_integration):** adaptive integrators assume the mass is where they look; slowly decaying integrands defeat their stopping rules, as your Step 4 printouts documented. ### Why it matters Fed the same two parameters, your two models priced the same one-day crash at $5.7\times10^{-63}$ versus $0.019$: sixty orders of magnitude apart, with reality (about one such day per century of trading) between them. Every risk report, capital requirement, and insurance premium built on a distributional assumption inherits a version of that gap, and the eight years of DAX data you analyzed produced two 5-sigma days where the normal model budgeted $0.001$ of one.