Hechos estilizados: colas gruesas, agrupamiento, autocorrelación y apalancamiento.

Cuatro regularidades

Un hecho estilizado es una regularidad estadística que aparece en las rentabilidades de muchos activos, mercados y épocas, aunque los detalles cambien. Cont (2001) los recoge en un catálogo que es la referencia clásica11Nota 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.. Esta lección mide los cuatro que más condicionan un sistema de trading:

flowchart TB
  H[Hechos estilizados] --> C[Colas gruesas]
  H --> A[Agrupamiento]
  H --> D[Autocorr. débil]
  H --> L[Apalancamiento]
  C --> R1[VaR y stops]
  A --> R2[Sizing por vol]
  D --> R3[Señales frágiles]
  L --> R4[Riesgo asimétrico]

La lección 0.7.2 ya comparó las colas de un ETF ficticio con la normal y con una t de Student. Aquí se usan los dos índices del nivel 1, ISX y FIX, con siete años de sesiones de la NYSE, y se mide cada hecho con la herramienta que le corresponde.

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 (cripto, diario)"] = resumen(r_xsn)
tabla["Universo (mediana de 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 (cripto, diario)       730.0      -0.01             2.43
Universo (mediana de 57)  2765.0       0.02             0.59

ISX y FIX tienen una curtosis en exceso de 6,7 y 8,5: muy lejos del 0 de la normal. Las otras dos filas son la advertencia del aviso: el token XSN y los 57 valores del universo salen de otros modelos del generador del curso, con colas más finas (2,4 en XSN) o casi normales (0,6 de mediana en el universo). Los hechos estilizados de unos datos sintéticos son los que su autor decidió poner. En el resto de la lección se usan ISX y FIX.

Colas gruesas

La curtosis resume las colas en un número, y es un número inestable (se verá más abajo). Lo que importa para el riesgo es más directo: cuántos días extremos hay frente a los que predice una normal con la misma media y desviación típica.

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

En ISX hay 35 días a más de tres desviaciones típicas de la media donde la normal espera 4,7; 14 a más de cuatro donde espera 0,11, y 8 a más de cinco donde espera una milésima. Cuanto más lejos en la cola, mayor el error de la normal: siete veces en ±3σ, más de cien en ±4σ, miles en ±5σ.

Los cuantiles lo enseñan de otra forma. Si las rentabilidades fueran normales, el cuantil del 1 % de las estandarizadas sería −2,33; el del 0,1 %, −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

La figura de la cabecera dibuja esta comparación completa (un gráfico cuantil-cuantil): si ISX fuera normal, los puntos estarían sobre la diagonal. En el centro casi lo están; en los extremos se separan. Y hay un detalle menos visible: en el centro, los cuantiles del 25 y el 75 % de ISX están más cerca de cero que los de la normal (−0,44 y 0,41 frente a ±0,67), mientras que el del 0,1 % está en −5,56 en lugar de −3,09. Una distribución de colas gruesas con la misma desviación típica tiene más días tranquilos y más días extremos; le faltan los intermedios.

El índice de cola

Las colas de las rentabilidades se parecen más a una ley de potencia que a la cola exponencial de la normal: la probabilidad de una pérdida mayor que x cae como x^(−α). El estimador de Hill estima α con las k pérdidas más grandes:

donde x₍₁₎ ≥ x₍₂₎ ≥ … son las pérdidas ordenadas de mayor a menor.

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}

Con estas muestras, α sale entre 1,7 y 3,1, y depende de cuántas pérdidas se usen (k): el estimador es sensible a esa elección y conviene mirar varios valores. La interpretación importa más que la cifra exacta: con una cola de potencia de índice α solo existen los momentos de orden menor que α. Con α por debajo de 4, la curtosis teórica no es finita y la muestral no se estabiliza nunca:

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

De 1,4 a 6,4 en ISX y de 2,2 a 6,9 en FIX según el año. Un único número de curtosis para describir un activo dice menos de lo que parece: un año sin días extremos tiene curtosis baja y el siguiente, con dos, altísima.

Agregar no siempre normaliza

Al sumar rentabilidades diarias en semanales o mensuales, el teorema central del límite empuja hacia la normal y la curtosis debería bajar. Ocurre, pero despacio y con ruido:

filas = {}
for h, nombre in [(1, "diaria"), (5, "semanal"), (21, "mensual")]:
    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
diaria   1759.0  6.65   8.54
semanal   351.0  2.89  12.69
mensual    83.0  9.95   5.84

En ISX la curtosis cae de 6,7 en diario a 2,9 en semanal, pero vuelve a subir a 9,9 en mensual: con 83 meses, uno solo (el desplome de marzo de 2020) basta para disparar la estimación. En FIX pasa al revés: sube a 12,7 en semanal y baja a 5,8 en mensual. La agregación reduce las colas, pero las muestras cortas de horizontes largos no permiten medirlo bien. Por eso los riesgos de cola se estudian con datos diarios o intradía y luego se agregan con un modelo.

Agrupamiento de volatilidad y autocorrelación débil

El segundo y el tercer hecho van juntos. La autocorrelación (ACF) de una serie en el retardo k es la correlación entre la serie y ella misma k sesiones antes. Para las rentabilidades es casi nula: la dirección de ayer apenas dice nada de la de hoy. Para su valor absoluto, que mide el tamaño del movimiento, es claramente positiva y decae despacio.

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"Banda del 95 %: ±{banda:.3f}")
pd.DataFrame(res).round(3)

Banda del 95 %: ±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

La autocorrelación de las rentabilidades de ISX en el retardo 1 es de 0,045, justo dentro de la banda de ruido (±0,047), y solo 3 de los 20 primeros retardos se salen de ella en ISX (4 en FIX). La de |r| está fuera de la banda en los 20 retardos de las dos series: 0,25 en el retardo 1 y todavía 0,18 a 20 sesiones en ISX. Es el agrupamiento de volatilidad: los días agitados vienen en rachas. Es lo que hace que la EWMA de la lección 1.6.1 funcione y la base de los modelos GARCH (lección 2.1.4).

El contraste de Ljung-Box resume los primeros retardos en un solo p-valor (hipótesis nula: ninguna autocorrelación hasta el retardo 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

Para |r| el p-valor es cero a efectos prácticos: el agrupamiento es inequívoco. Para r, el contraste rechaza la ausencia de autocorrelación al 5 % en las dos series. ¿Hay entonces algo que predecir en la dirección? Poco: con autocorrelaciones de 0,05-0,10, una regla que siga la dirección de ayer se comería esa ventaja en costes (lección 1.1.5), y en estos datos esa autocorrelación sale del propio modelo (la tendencia entre fechas fijas cambia la rentabilidad media de un tramo a otro). En datos reales pasa algo parecido: autocorrelaciones pequeñas, que cambian de signo entre periodos y que los costes suelen anular. Un contraste significativo no es una estrategia.

El agrupamiento también se ve sin ACF, en probabilidades:

x = r["ISX"]
extremo = x.abs() >= x.abs().quantile(0.9)    # decil superior
p_base = extremo.mean()
p_tras = extremo[extremo.shift(1, fill_value=False)].mean()
print(f"P(día extremo) = {p_base:.1%}")
print(f"P(día extremo | ayer extremo) = {p_tras:.1%}")

P(día extremo) = 10.0%
P(día extremo | ayer extremo) = 24.4%

Uno de cada diez días está en el decil de movimientos más grandes por definición; después de uno de ellos, la probabilidad de otro se multiplica por más de dos.

Una forma útil de pensarlo es como dos estados que persisten: el mercado pasa largos periodos en calma, un choque lo lleva a un estado agitado y de ahí sale poco a poco. La EWMA y los GARCH son maneras de estimar en qué punto de ese ciclo está hoy.

stateDiagram-v2
  state "Agitación" as Agitacion
  [*] --> Calma
  Calma --> Calma: día tranquilo
  Calma --> Agitacion: choque
  Agitacion --> Agitacion: más choques
  Agitacion --> Calma: se disipa

Efecto apalancamiento

El cuarto hecho es una asimetría: la volatilidad tiende a subir más después de una caída que después de una subida del mismo tamaño. Se llama efecto apalancamiento por una explicación clásica (cuando cae la cotización, sube el apalancamiento financiero de la empresa), aunque en índices y en otros activos se atribuye también a la reacción de los inversores ante las pérdidas. Se mide relacionando la rentabilidad de hoy con la volatilidad futura:

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

Los dos extremos tienen más volatilidad futura que el centro: es otra vez el agrupamiento (un día grande, de cualquier signo, anuncia más días grandes). La asimetría está en los extremos: tras el decil de las peores sesiones, la volatilidad de las diez siguientes es del 26,9 % anual; tras el de las mejores, del 23,8 %. El mínimo, 13,3 %, está en el decil 5, el de los días casi planos.

La forma más limpia de medirlo elimina el efecto del nivel de volatilidad: estandariza cada día por su volatilidad prevista (EWMA, lección 1.6.1) y compara lo que viene después de un día de z < −1 con lo que viene después de uno de 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[z² de mañana] = "
          f"{z2_siguiente[filtro].mean():.2f}  (n = {filtro.sum()})")
print(f"Curtosis de r: {x.kurt():.2f}; de z: {z.kurt():.2f}")

       z < -1: E[z² de mañana] = 1.84  (n = 194)
 -1 <= z <= 1: E[z² de mañana] = 1.05  (n = 1264)
        z > 1: E[z² de mañana] = 1.28  (n = 238)
Curtosis de r: 6.65; de z: 6.87

Si la volatilidad prevista recogiera todo, el cuadrado de mañana estandarizado valdría 1 de media en las tres filas. Después de un día tranquilo vale casi exactamente eso (1,05); después de una caída de más de una desviación típica, 1,84; después de una subida de la misma magnitud, 1,28. Las caídas mueven la volatilidad bastante más que las subidas, y más de lo que la EWMA, que es simétrica, ya ha incorporado.

La última línea responde a una pregunta pendiente de la lección 1.6.1: estandarizar por una volatilidad cambiante no quita las colas. La curtosis de z sigue en torno a 7. En estos datos es así por construcción: la distribución de cada día es una t de Student, con colas gruesas por sí misma, y la volatilidad cambiante añade más colas a la serie completa. En datos reales también suele quedar curtosis después de estandarizar, y por eso los modelos de riesgo usan distribuciones t o históricas en lugar de la normal.

Lo que cada hecho cambia en un sistema

Una consecuencia común a los cuatro: los datos sintéticos que se usen para probar un sistema (nivel 2, lección 2.5.7) tienen que reproducir estos hechos. Un paseo aleatorio con rentabilidades normales e independientes, el modelo más sencillo, no tiene ninguno de ellos, y un sistema probado solo con él no ha visto nunca un mercado de verdad.

Ejercicios

1. Repite la tabla de días extremos con las rentabilidades estandarizadas por la EWMA (z de la última celda). ¿Se acercan a la 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())} observados, "
          f"{esperados:.3f} esperados")

|z| > 3: 23 observados, 4.579 esperados
|z| > 4: 11 observados, 0.107 esperados
|z| > 5: 4 observados, 0.001 esperados

2. Calcula la correlación entre la rentabilidad de hoy y la volatilidad de 10 sesiones de los días siguientes, en ISX y en FIX. ¿Qué signo tiene y qué tamaño?

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. Un colega te enseña un backtest de una estrategia de reversión que «gana en 9 de cada 10 operaciones» probada con rentabilidades simuladas normales e independientes. ¿Qué hechos estilizados faltan en esa simulación y qué pueden cambiar?