Fat tails, or why the normal lies

import numpy as np
import qadata
from scipy import stats

x = qadata.fat_tailed_returns().to_numpy()   # planted t(4), 20% ann vol
n = len(x)

z = (x - x.mean()) / x.std()                 # numpy std: ddof=0, population
S = (z ** 3).mean()
K = (z ** 4).mean()

print(f"n = {n}, daily std = {x.std(ddof=1):.6%}, ann vol = {x.std(ddof=1) * np.sqrt(252):.2%}")
print(f"skewness : hand {S:.15f}  scipy {stats.skew(x):.15f}")
print(f"kurtosis : hand {K:.15f}  scipy {stats.kurtosis(x, fisher=False):.15f}")
print(f"excess kurtosis: {K - 3:.4f}   (a normal would be ~0)")

import numpy as np
import qadata
from scipy import stats

x = qadata.fat_tailed_returns().to_numpy()
n = len(x)
z = (x - x.mean()) / x.std(ddof=1)

for k in (3, 4):
    obs = int((np.abs(z) > k).sum())
    exp = 2 * stats.norm.sf(k) * n
    print(f"|z| > {k}: observed {obs:3d}   normal expects {exp:6.2f}")

zmax = np.abs(z).max()
p = 2 * stats.norm.sf(zmax)
print(f"biggest move: {x[np.abs(z).argmax()]:+.4%} = {zmax:.2f} sigma")
print(f"normal frequency: once per {1 / p / 252:,.0f} years")

import numpy as np
import qadata
from scipy import stats

x = qadata.fat_tailed_returns().to_numpy()
n = len(x)
rng = np.random.default_rng(99)
g = x.mean() + x.std(ddof=1) * rng.standard_normal(n)   # the matched twin

for name, v in [("t(4) series", x), ("matched normal", g)]:
    z = (v - v.mean()) / v.std()
    S, K = (z ** 3).mean(), (z ** 4).mean()
    jb = n / 6 * (S ** 2 + (K - 3) ** 2 / 4)
    res = stats.jarque_bera(v)
    print(f"{name:15s}: JB hand {jb:10.4f}  scipy {res.statistic:10.4f}  p = {res.pvalue:.4f}")
print(f"5% critical value (chi2, 2 dof): {stats.chi2.ppf(0.95, 2):.4f}")