Stylized facts: fat tails, clustering, autocorrelation and leverage.

Four regularities

A stylized fact is a statistical regularity that shows up in the returns of many assets, markets and periods, even if the details change. Cont (2001) collects them in a catalog that is the classic reference11Note 1 · R. Cont, “Empirical properties of asset returns: stylized facts and statistical issues”, Quantitative Finance 1(2), 2001, 223-236, https://doi.org/10.1080/713665670.. This lesson measures the four that most shape a trading system:

flowchart TB
  H[Stylized facts] --> C[Fat tails]
  H --> A[Clustering]
  H --> D[Weak autocorr.]
  H --> L[Leverage]
  C --> R1[VaR and stops]
  A --> R2[Vol-based sizing]
  D --> R3[Fragile signals]
  L --> R4[Asymmetric risk]

Lesson 0.7.2 already compared the tails of a fictional ETF with the normal and with a Student's t. Here the two level 1 indices, ISX and FIX, are used, with seven years of NYSE sessions, and each fact is measured with the tool that suits it.

import numpy as np
import pandas as pd
from scipy import stats

isx = pd.read_csv("etf_indice_diario.csv", parse_dates=["fecha"],
                  index_col="fecha")["indice"]
fix = pd.read_csv("fix_indice_diario.csv", parse_dates=["fecha"],
                  index_col="fecha")["indice"]
r = np.log(pd.DataFrame({"ISX": isx, "FIX": fix})).diff().dropna()

xsn = pd.read_parquet("xsn_1h.parquet").set_index("ts")["spot"]
r_xsn = np.log(xsn.resample("D").last()).diff().dropna()
uni = pd.read_parquet("universo_precios.parquet").pivot(
    index="fecha", columns="ticker", values="cierre")
ru = np.log(uni).diff()

def resumen(s):
    return {"n": s.count(), "asimetria": s.skew(),
            "curtosis_exceso": s.kurt()}

tabla = {c: resumen(r[c]) for c in r}
tabla["XSN (crypto, daily)"] = resumen(r_xsn)
tabla["Universe (median of 57)"] = {
    "n": ru.count().median(), "asimetria": ru.skew().median(),
    "curtosis_exceso": ru.kurt().median()}
pd.DataFrame(tabla).T.round(2)

                              n  asimetria  curtosis_exceso
ISX                      1759.0       0.19             6.65
FIX                      1759.0      -0.28             8.54
XSN (crypto, daily)       730.0      -0.01             2.43
Universe (median of 57)  2765.0       0.02             0.59

ISX and FIX have an excess kurtosis of 6.7 and 8.5: very far from the normal's 0. The other two rows are the point of the warning: the XSN token and the 57 stocks of the universe come from other models in the course's generator, with thinner tails (2.4 for XSN) or almost normal ones (a median of 0.6 in the universe). The stylized facts of synthetic data are the ones its author decided to put in. The rest of the lesson uses ISX and FIX.

Fat tails

Kurtosis sums up the tails in a single number, and it is an unstable number (as shown below). What matters for risk is more direct: how many extreme days there are versus those predicted by a normal with the same mean and standard deviation.

filas = {}
for c in r:
    z = (r[c] - r[c].mean()) / r[c].std()
    for k in (3, 4, 5):
        esperados = len(z) * 2 * stats.norm.sf(k)
        filas[(c, k)] = {"observados": int((z.abs() > k).sum()),
                         "normal": esperados,
                         "veces": (z.abs() > k).sum() / esperados}
pd.DataFrame(filas).T.round(3)

       observados  normal     veces
ISX 3        35.0   4.749     7.370
    4        14.0   0.111   125.651
    5         8.0   0.001  7933.043
FIX 3        24.0   4.749     5.054
    4        12.0   0.111   107.701
    5         9.0   0.001  8924.674

ISX has 35 days more than three standard deviations from the mean where the normal expects 4.7; 14 beyond four where it expects 0.11, and 8 beyond five where it expects one thousandth. The further out in the tail, the larger the normal's error: seven times at ±3σ, more than a hundred at ±4σ, thousands at ±5σ.

Quantiles show it another way. If returns were normal, the 1% quantile of the standardized returns would be −2.33; the 0.1% quantile, −3.09:

z = (r["ISX"] - r["ISX"].mean()) / r["ISX"].std()
p = [0.001, 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99, 0.999]
pd.DataFrame({"normal": stats.norm.ppf(p),
              "ISX": z.quantile(p).to_numpy()},
             index=[f"{100 * q:g}%" for q in p]).round(2)

       normal   ISX
0.1%    -3.09 -5.56
1%      -2.33 -2.62
5%      -1.64 -1.46
25%     -0.67 -0.44
50%      0.00 -0.01
75%      0.67  0.41
95%      1.64  1.50
99%      2.33  3.25
99.9%    3.09  5.44

The figure at the top draws this full comparison (a quantile-quantile plot): if ISX were normal, the points would lie on the diagonal. In the center they almost do; at the extremes they pull away. And there is a less visible detail: in the center, ISX's 25% and 75% quantiles are closer to zero than the normal's (−0.44 and 0.41 versus ±0.67), while the 0.1% quantile is at −5.56 instead of −3.09. A fat-tailed distribution with the same standard deviation has more quiet days and more extreme days; it is short on the ones in between.

The tail index

The tails of returns look more like a power law than like the normal's exponential tail: the probability of a loss larger than x falls like x^(−α). The Hill estimator estimates α with the k largest losses:

where x₍₁₎ ≥ x₍₂₎ ≥ … are the losses sorted from largest to smallest.

def hill(perdidas, k):
    x = np.sort(perdidas[perdidas > 0])[::-1]
    return 1 / np.mean(np.log(x[:k] / x[k]))

for c in r:
    est = {k: round(float(hill(-r[c].to_numpy(), k)), 2)
           for k in (25, 50, 100, 200)}
    print(c, est)

ISX {25: 2.99, 50: 2.69, 100: 2.25, 200: 1.97}
FIX {25: 2.47, 50: 3.05, 100: 2.45, 200: 1.73}

With these samples, α comes out between 1.7 and 3.1, and it depends on how many losses are used (k): the estimator is sensitive to that choice and it is worth looking at several values. The interpretation matters more than the exact figure: with a power-law tail of index α, only moments of order lower than α exist. With α below 4, theoretical kurtosis is not finite and the sample kurtosis never settles:

r.groupby(r.index.year).kurt().round(2).T

fecha  2019  2020  2021  2022  2023  2024  2025
ISX    1.40  3.07  3.24  5.31  2.88  5.39  6.39
FIX    5.07  2.84  6.87  5.46  4.34  2.24  6.00

From 1.4 to 6.4 in ISX and from 2.2 to 6.9 in FIX depending on the year. A single kurtosis number to describe an asset says less than it seems: a year without extreme days has low kurtosis and the next, with two of them, a very high one.

Aggregating does not always normalize

When daily returns are summed into weekly or monthly ones, the central limit theorem pushes toward the normal and kurtosis should fall. It does, but slowly and noisily:

filas = {}
for h, nombre in [(1, "daily"), (5, "weekly"), (21, "monthly")]:
    agg = np.log(pd.DataFrame({"ISX": isx, "FIX": fix})).diff(h)
    agg = agg.dropna().iloc[::h]
    filas[nombre] = {"n": len(agg), **agg.kurt().round(2)}
pd.DataFrame(filas).T

              n   ISX    FIX
daily    1759.0  6.65   8.54
weekly    351.0  2.89  12.69
monthly    83.0  9.95   5.84

In ISX kurtosis falls from 6.7 at daily frequency to 2.9 at weekly, but rises again to 9.9 at monthly: with 83 months, a single one (the March 2020 crash) is enough to blow up the estimate. In FIX the opposite happens: it rises to 12.7 at weekly and falls to 5.8 at monthly. Aggregation reduces the tails, but the short samples of long horizons do not allow it to be measured well. That is why tail risks are studied with daily or intraday data and then aggregated with a model.

Volatility clustering and weak autocorrelation

The second and third facts go together. The autocorrelation (ACF) of a series at lag k is the correlation between the series and itself k sessions earlier. For returns it is almost zero: yesterday's direction says next to nothing about today's. For their absolute value, which measures the size of the move, it is clearly positive and decays slowly.

from statsmodels.tsa.stattools import acf

n = len(r)
banda = 1.96 / np.sqrt(n)
res = {}
for c in r:
    a, b = acf(r[c], nlags=20), acf(r[c].abs(), nlags=20)
    res[c] = {"r_k1": a[1], "|r|_k1": b[1], "|r|_k5": b[5],
              "|r|_k10": b[10], "|r|_k20": b[20],
              "r_fuera_banda": int((np.abs(a[1:]) > banda).sum()),
              "|r|_fuera_banda": int((np.abs(b[1:]) > banda).sum())}
print(f"95% band: ±{banda:.3f}")
pd.DataFrame(res).round(3)

95% band: ±0.047
                    ISX     FIX
r_k1              0.045   0.020
|r|_k1            0.248   0.294
|r|_k5            0.253   0.242
|r|_k10           0.198   0.273
|r|_k20           0.180   0.148
r_fuera_banda     3.000   4.000
|r|_fuera_banda  20.000  20.000

The autocorrelation of ISX returns at lag 1 is 0.045, just inside the noise band (±0.047), and only 3 of the first 20 lags fall outside it in ISX (4 in FIX). That of |r| is outside the band at all 20 lags in both series: 0.25 at lag 1 and still 0.18 at 20 sessions in ISX. This is volatility clustering: turbulent days come in streaks. It is what makes the EWMA of lesson 1.6.1 work and the basis of GARCH models (lesson 2.1.4).

The Ljung-Box test sums up the first lags in a single p-value (null hypothesis: no autocorrelation up to lag 10):

from statsmodels.stats.diagnostic import acorr_ljungbox

lb = {f"{c} {serie}": acorr_ljungbox(x, lags=[10]).iloc[0]
      for c in r for serie, x in [("r", r[c]), ("|r|", r[c].abs())]}
pd.DataFrame(lb).T.round(4)

           lb_stat  lb_pvalue
ISX r      22.2468     0.0139
ISX |r|   808.6959     0.0000
FIX r      47.3878     0.0000
FIX |r|  1222.0450     0.0000

For |r| the p-value is zero for all practical purposes: the clustering is unmistakable. For r, the test rejects the absence of autocorrelation at 5% in both series. Is there, then, something to predict in the direction? Little: with autocorrelations of 0.05-0.10, a rule that follows yesterday's direction would eat that edge in costs (lesson 1.1.5), and in this data that autocorrelation comes from the model itself (the trend between fixed dates changes the mean return from one segment to the next). On real data something similar happens: small autocorrelations, which change sign between periods and which costs usually wipe out. A significant test is not a strategy.

Clustering can also be seen without the ACF, in probabilities:

x = r["ISX"]
extremo = x.abs() >= x.abs().quantile(0.9)    # top decile
p_base = extremo.mean()
p_tras = extremo[extremo.shift(1, fill_value=False)].mean()
print(f"P(extreme day) = {p_base:.1%}")
print(f"P(extreme day | extreme yesterday) = {p_tras:.1%}")

P(extreme day) = 10.0%
P(extreme day | extreme yesterday) = 24.4%

One day in ten is in the decile of the largest moves by definition; after one of them, the probability of another more than doubles.

A useful way to think about it is as two states that persist: the market spends long periods calm, a shock moves it into a turbulent state, and it comes out of it gradually. The EWMA and GARCH models are ways of estimating where in that cycle the market is today.

stateDiagram-v2
  state "Turbulence" as Turbulent
  [*] --> Calm
  Calm --> Calm: quiet day
  Calm --> Turbulent: shock
  Turbulent --> Turbulent: more shocks
  Turbulent --> Calm: dissipates

Leverage effect

The fourth fact is an asymmetry: volatility tends to rise more after a fall than after a rise of the same size. It is called the leverage effect after a classic explanation (when the share price falls, the company's financial leverage rises), although in indices and other assets it is also attributed to investors' reaction to losses. It is measured by relating today's return to future volatility:

fut = np.sqrt(252 * x.pow(2).rolling(10).mean().shift(-10))
decil = pd.qcut(x, 10, labels=range(1, 11))
tabla = pd.DataFrame({
    "r_media_%": 100 * x.groupby(decil, observed=True).mean(),
    "vol_futura_%": 100 * fut.groupby(decil, observed=True).mean()})
tabla.round(2)

     r_media_%  vol_futura_%
ISX
1        -2.26         26.90
2        -0.94         18.63
3        -0.52         15.86
4        -0.27         15.45
5        -0.06         13.27
6         0.14         14.12
7         0.33         15.14
8         0.58         14.79
9         1.00         17.88
10        2.51         23.82

Both extremes have more future volatility than the center: it is clustering again (a big day, of either sign, heralds more big days). The asymmetry is at the extremes: after the decile of the worst sessions, the volatility of the next ten is 26.9% annualized; after that of the best, 23.8%. The minimum, 13.3%, is in decile 5, that of almost flat days.

The cleanest way to measure it removes the effect of the volatility level: standardize each day by its forecast volatility (EWMA, lesson 1.6.1) and compare what follows a day with z < −1 with what follows one with z > +1.

sig = np.sqrt(x.pow(2).ewm(alpha=0.06, adjust=False).mean().shift(1))
z = (x / sig).iloc[63:]
z2_siguiente = z.shift(-1) ** 2
for nombre, filtro in [("z < -1", z < -1), ("-1 <= z <= 1",
                       z.abs() <= 1), ("z > 1", z > 1)]:
    print(f"{nombre:>13}: E[tomorrow's z²] = "
          f"{z2_siguiente[filtro].mean():.2f}  (n = {filtro.sum()})")
print(f"Kurtosis of r: {x.kurt():.2f}; of z: {z.kurt():.2f}")

       z < -1: E[tomorrow's z²] = 1.84  (n = 194)
 -1 <= z <= 1: E[tomorrow's z²] = 1.05  (n = 1264)
        z > 1: E[tomorrow's z²] = 1.28  (n = 238)
Kurtosis of r: 6.65; of z: 6.87

If the forecast volatility captured everything, tomorrow's standardized square would average 1 in all three rows. After a quiet day it is almost exactly that (1.05); after a fall of more than one standard deviation, 1.84; after a rise of the same magnitude, 1.28. Falls move volatility considerably more than rises, and more than the EWMA, which is symmetric, has already taken in.

The last line answers a question left open in lesson 1.6.1: standardizing by a changing volatility does not remove the tails. The kurtosis of z stays at around 7. In this data that is so by construction: each day's distribution is a Student's t, fat-tailed in itself, and the changing volatility adds more tail to the full series. On real data some kurtosis also tends to remain after standardizing, and that is why risk models use t or historical distributions instead of the normal.

What each fact changes in a system

A consequence common to all four: the synthetic data used to test a system (level 2, lesson 2.5.7) has to reproduce these facts. A random walk with normal, independent returns, the simplest model, has none of them, and a system tested only on it has never seen a real market.

Exercises

1. Repeat the extreme days table with the returns standardized by the EWMA (z from the last cell). Do they get closer to the normal?

ze = (z - z.mean()) / z.std()
for k in (3, 4, 5):
    esperados = len(ze) * 2 * stats.norm.sf(k)
    print(f"|z| > {k}: {int((ze.abs() > k).sum())} observed, "
          f"{esperados:.3f} expected")

|z| > 3: 23 observed, 4.579 expected
|z| > 4: 11 observed, 0.107 expected
|z| > 5: 4 observed, 0.001 expected

2. Compute the correlation between today's return and the 10-session volatility of the following days, in ISX and in FIX. What sign does it have and how large is it?

for c in r:
    f10 = r[c].pow(2).rolling(10).mean().shift(-10)
    print(f"{c}: {r[c].corr(np.sqrt(f10)):.3f}")

ISX: -0.066
FIX: -0.083

3. A colleague shows you a backtest of a mean-reversion strategy that “wins 9 out of 10 trades”, tested with simulated normal, independent returns. Which stylized facts are missing from that simulation and what can they change?