CNN 1D, LSTM y GRU: ventanas, escalado y fugas.

Secuencias, no filas

En las lecciones anteriores cada fila era independiente de las demás: un vector de features con su etiqueta. Una serie temporal se puede tratar así si las features ya resumen la historia (el HAR resume 22 días en tres medias). Las redes de secuencias se saltan ese paso: reciben la ventana de los últimos días tal cual, una matriz de $L$ pasos por $k$ variables, y aprenden su propio resumen.

• CNN 1D: filtros de unos pocos pasos que recorren la ventana y detectan patrones locales (un salto, un cambio de pendiente), los mismos en cualquier posición.
• LSTM (Hochreiter y Schmidhuber, 1997) y GRU (Cho y otros, 2014): leen la ventana paso a paso y mantienen un estado interno; unas puertas aprendidas deciden qué olvidar y qué añadir en cada paso. La GRU tiene dos puertas y la LSTM tres, más un estado de celda aparte: la GRU tiene menos parámetros con el mismo tamaño de estado.

El banco de pruebas es el de volatilidad del nivel 5: seis activos sintéticos (índice, acción, bono, oro, divisa y cripto) de 2008 a 2025 con su varianza realizada diaria. El proceso verdadero es un HAR en el logaritmo de la varianza integrada, con efecto apalancamiento (la parte negativa del rendimiento de ayer sube la volatilidad), saltos en la acción y la cripto, y cambios de nivel en fechas conocidas (README de los datos). Un HAR bien estimado es casi el óptimo: la pregunta es si una red llega ahí sin que se lo digan, y si lo pasa.

El corte es cronológico: entrenamiento de 2008 a 2015, validación (para la parada temprana) 2016-2017 y prueba 2018-2025, la misma que la tabla de previsiones del README. Primero, el nivel medio del log RV en cada tramo:

import os
os.environ["KERAS_BACKEND"] = "jax"      # antes del primer import keras
import numpy as np
import pandas as pd
import statsmodels.api as sm
import keras
from keras import layers

v = pd.read_parquet("vol5_diario.parquet")
verdad = pd.read_parquet("vol5_verdad.parquet")
v = v.merge(verdad[["fecha", "activo", "iv_prevista"]],
            on=["fecha", "activo"]).sort_values(["activo", "fecha"])
v["lrv"] = np.log(v["rv_5min"])
F_VAL = np.datetime64("2016-01-01")
F_PRUEBA = np.datetime64("2018-01-01")
tramo = np.select([v["fecha"] < F_VAL, v["fecha"] < F_PRUEBA],
                  ["entreno", "validación"], "prueba")
niveles = v.groupby(["activo", tramo])["lrv"].mean().unstack()
niveles = niveles[["entreno", "validación", "prueba"]]
niveles["prueba - entreno"] = niveles["prueba"] - niveles["entreno"]
niveles.round(2)

        entreno  validación  prueba  prueba - entreno
activo
accion    -8.99       -9.31   -9.24             -0.25
bono     -11.51      -11.61  -11.34              0.17
cripto    -6.95       -7.22   -7.32             -0.37
divisa   -10.60      -10.13  -10.41              0.19
indice    -9.85       -9.54   -9.72              0.13
oro       -9.56       -9.92   -9.75             -0.18

La acción y la cripto bajan de nivel entre el entrenamiento y la prueba (−0,25 y −0,37 en logaritmos; en la acción el cambio de 2015 ya está dentro del entrenamiento, que mezcla los dos regímenes); el bono sube (+0,17, por el cambio de 2022) y la divisa también (+0,19: su cambio de nivel de enero de 2015 deja casi todo el entrenamiento en el régimen anterior). Una serie no estacionaria en nivel es lo normal en finanzas, y condiciona cómo se escala.

Ventanas sin mirar al futuro

La función que corta las ventanas es la pieza donde más fácil es equivocarse en un paso. La regla: la ventana que termina en el día $t$ contiene los días $t-L+1$ a $t$, y su objetivo es el día $t+1$. Con sliding_window_view de NumPy cada ventana es una vista, sin copiar datos. Antes de usarla con datos reales, se prueba con una serie cuyos valores son su propia posición:

def ventanas(z, L):
    """Ventanas de L filas que terminan en t (t = L-1 … n-2) y el
    valor de la columna 0 en t + 1."""
    w = np.lib.stride_tricks.sliding_window_view(
        z, (L, z.shape[1]))[:, 0]
    return w[:-1], z[L:, 0]

juguete = np.arange(8, dtype="float32")[:, None]
Xj, yj = ventanas(juguete, 3)
for xv, yv in zip(Xj[:3], yj[:3]):
    print(xv.ravel().tolist(), "->", float(yv))
print("ventanas:", len(Xj), "| último objetivo:", float(yj[-1]))
assert all(xv[-1, 0] + 1 == yv for xv, yv in zip(Xj, yj))

[0.0, 1.0, 2.0] -> 3.0
[1.0, 2.0, 3.0] -> 4.0
[2.0, 3.0, 4.0] -> 5.0
ventanas: 5 | último objetivo: 7.0

Con 8 valores y $L = 3$ salen 5 ventanas: la primera, [0, 1, 2], predice el 3, y la última predice el 7, el último valor de la serie. La aserción comprueba en todas que el objetivo es el día siguiente al último de la ventana. Si alguien cambia la función, la aserción falla antes de entrenar nada.

El escalado también puede mirar al futuro

Las redes entrenan mejor con entradas de escala parecida, así que se tipifica cada serie: log RV menos su media, entre su desviación. La pregunta es qué media. La de toda la muestra incluye 2018-2025, y si la serie cambia de nivel, esa media le dice al modelo dónde va a estar. Para medir cuánto vale esa información sin entrenar ninguna red, la previsión más simple posible: predecir que el log RV de cada día de prueba será la media del activo. Con la media del entrenamiento es una previsión legítima; con la de toda la muestra, no:

def stats(s, hasta):
    """Media y desviación de log RV y desviación del rendimiento,
    con las filas anteriores a `hasta`."""
    m = s["fecha"] < hasta
    return (s.loc[m, "lrv"].mean(), s.loc[m, "lrv"].std(),
            s.loc[m, "rendimiento"].std())

filas = []
for a, s in v.groupby("activo"):
    prueba = s["fecha"] >= F_PRUEBA
    real = s.loc[prueba, "lrv"]
    m_ent = stats(s, F_VAL)[0]
    m_todo = s["lrv"].mean()                 # usa 2018-2025
    filas.append({"activo": a,
                  "media de entreno": np.mean((real - m_ent) ** 2),
                  "media de toda la muestra":
                      np.mean((real - m_todo) ** 2)})
fuga = pd.DataFrame(filas).set_index("activo")
fuga.loc["media"] = fuga.mean()
fuga.round(3)

        media de entreno  media de toda la muestra
activo
accion             0.544                     0.491
bono               0.292                     0.275
cripto             0.795                     0.689
divisa             0.254                     0.221
indice             0.418                     0.402
oro                0.371                     0.341
media              0.446                     0.403

Calcular la media con toda la muestra baja el MSE de prueba de esa previsión de 0,446 a 0,403 de media, un 9,6 % menos, sin ningún modelo: más de un 13 % en la cripto y la divisa. Una red que recibe entradas y objetivos tipificados así empieza con esa ventaja regalada, porque su sesgo inicial (predecir 0) ya es la media del futuro. Es la misma fuga que el StandardScaler ajustado antes de partir los datos (lección 5.2.4), y en series con cambios de nivel pesa más que en un panel transversal.

La regla: media y desviación de cada activo con su tramo de entrenamiento, y esas mismas cifras para validación y prueba. Cada ventana lleva dos variables por día: el log RV tipificado y el rendimiento dividido por su desviación de entrenamiento. Se guarda también la desviación de cada activo, para devolver los errores a unidades de log RV, y el oráculo (iv_prevista, la esperanza exacta de la varianza integrada con la información hasta el día anterior, que solo existe porque los datos son sintéticos):

L = 22
ACTIVOS = sorted(v["activo"].unique())
partes = {k: [] for k in ("X", "y", "f", "a", "sd", "orac")}
for i, a in enumerate(ACTIVOS):
    s = v[v["activo"] == a]
    m, sd, sr = stats(s, F_VAL)              # solo entrenamiento
    z = np.c_[(s["lrv"] - m) / sd,
              s["rendimiento"] / sr].astype("float32")
    Xa, ya = ventanas(z, L)
    partes["X"].append(Xa)
    partes["y"].append(ya)
    partes["f"].append(s["fecha"].to_numpy()[L:])
    partes["a"].append(np.full(len(ya), i))
    partes["sd"].append(np.full(len(ya), sd))
    lor = (np.log(s["iv_prevista"].to_numpy()) - m) / sd
    partes["orac"].append(lor[L:])
X, y, fch, act, SD, orac = (np.concatenate(partes[k]) for k in
                            ("X", "y", "f", "a", "sd", "orac"))
ent = fch < F_VAL
val = (fch >= F_VAL) & (fch < F_PRUEBA)
te = fch >= F_PRUEBA
print("X:", X.shape, "| entreno, validación, prueba:",
      int(ent.sum()), int(val.sum()), int(te.sum()))

def mse(p, m):
    """MSE en unidades de log RV (deshace la escala de cada activo)."""
    return float(np.mean(((p[m] - y[m]) * SD[m]) ** 2))

X: (27042, 22, 2) | entreno, validación, prueba: 11958 3018 12066

Son 27.042 ventanas de 22 días por 2 variables: 11.958 de entrenamiento, 3.018 de validación y 12.066 de prueba. Las ventanas del final del entrenamiento y las del principio de la validación no se solapan en el objetivo (cada una predice un solo día), pero sí comparten días de entrada; con una sola frontera y la parada temprana como único uso de la validación, no hace falta purgar más.

Las baselines

Cuatro previsiones sin red, todas por activo y estimadas con el entrenamiento:

• Ingenua: el log RV de mañana es el de hoy.
• HAR: OLS sobre las medias de 1, 5 y 22 días del log RV, calculadas con la propia ventana.
• HAR con apalancamiento: añade $\ell_t = \min(r_t, 0)$, la parte negativa del rendimiento de hoy (el término $\gamma \ell$ del proceso verdadero).
• Oráculo calibrado: una recta sobre el logaritmo de iv_prevista. La recta corrige el sesgo entre la varianza integrada y la realizada medida (la RV no ve la parte nocturna y lleva error de medida y saltos). Es el suelo: ninguna previsión legítima puede bajar de él salvo por suerte.

def har(X):
    l = X[:, :, 0]
    return np.c_[l[:, -1], l[:, -5:].mean(1), l.mean(1)]

def por_activo(Z):
    """OLS de cada activo con su entrenamiento; predice todo."""
    p = np.zeros_like(y)
    for i in range(len(ACTIVOS)):
        mi = act == i
        A_ = np.c_[np.ones(mi.sum()), Z[mi]]
        b = np.linalg.lstsq(A_[ent[mi]], y[mi & ent], rcond=None)[0]
        p[mi] = A_ @ b
    return p

H = har(X)
HL = np.c_[H, np.minimum(X[:, -1, 1], 0)]   # parte negativa de r_t
pred = {"ingenua (RV de hoy)": X[:, -1, 0],
        "HAR": por_activo(H),
        "HAR + apalancamiento": por_activo(HL),
        "oráculo calibrado": por_activo(orac[:, None])}
resumen = pd.Series({k: mse(p, te) for k, p in pred.items()},
                    name="MSE prueba")
resumen.round(4)

ingenua (RV de hoy)     0.3184
HAR                     0.2137
HAR + apalancamiento    0.2114
oráculo calibrado       0.1962
Name: MSE prueba, dtype: float64

La ingenua tiene un MSE de 0,318; el HAR baja a 0,2137, el apalancamiento a 0,2114 y el oráculo marca 0,1962. Entre el HAR con apalancamiento y el oráculo queda un 7 %: es lo máximo que una red puede ganar aquí. Cualquier resultado por debajo de 0,1962 en la prueba es, casi seguro, una fuga.

Una forma frecuente de cometerla: suavizar con una media móvil centrada, que en pandas es un argumento (center=True) y en muchas librerías de indicadores es la opción por defecto de algún filtro. La celda cambia la media semanal del HAR por la centrada en $t$, que usa de $t-2$ a $t+2$:

c5 = []
for a in ACTIVOS:
    s = v[v["activo"] == a]
    m, sd, _ = stats(s, F_VAL)
    zc = ((s["lrv"] - m) / sd).rolling(5, center=True).mean()
    c5.append(zc.to_numpy()[L - 1:-1])    # alineada con el final t
c5 = np.nan_to_num(np.concatenate(c5))
H_mal = np.c_[H[:, 0], c5, H[:, 2]]   # t-2 … t+2 en vez de t-4 … t
print(f"HAR con la media semanal centrada: MSE "
      f"{mse(por_activo(H_mal), te):.4f}")

HAR con la media semanal centrada: MSE 0.1328

Un MSE de 0,1328, un 32 % por debajo del oráculo. El modelo no ha aprendido nada sobre la volatilidad: la media centrada contiene el día que se quiere prever ($t+1$) y el siguiente. Con una red el error es igual de fácil de cometer y más difícil de ver, porque no hay coeficientes que mirar: una capa Bidirectional sobre una secuencia que llega hasta el objetivo, un filtro aplicado a la serie entera antes de cortar ventanas o un relleno simétrico en una convolución que predice en cada paso. La comparación con el oráculo, cuando existe, y con una baseline honesta, cuando no, es la defensa.

Tres arquitecturas en Keras

Las tres reciben la misma entrada, (22, 2), y terminan en una Dense(1). Con el backend JAX, como en 5.6.2:

flowchart LR
  V[Ventana 22 x 2] --> C[Conv1D causal]
  V --> LS[LSTM 16]
  V --> G[GRU 16]
  C --> F[Flatten]
  F --> D[Dense 1]
  LS --> D
  G --> D

• La CNN 1D aplica 16 filtros de 5 días con relleno causal: la salida en la posición $t$ solo depende de $t-4$ a $t$. Después aplana las 22 posiciones y las combina en la capa densa; así conserva dónde está cada patrón (un GlobalAveragePooling1D lo perdería, y en volatilidad lo reciente importa más).
• La LSTM y la GRU leen los 22 días en orden y pasan su estado final, de 16 números, a la capa densa.

def red(tipo, n=16):
    x_in = keras.Input((L, 2))
    if tipo == "CNN 1D":
        h = layers.Conv1D(n, 5, activation="relu",
                          padding="causal")(x_in)
        h = layers.Flatten()(h)
    elif tipo == "LSTM":
        h = layers.LSTM(n)(x_in)
    else:
        h = layers.GRU(n)(x_in)
    m = keras.Model(x_in, layers.Dense(1)(h),
                    name=tipo.replace(" ", ""))
    m.compile(keras.optimizers.Adam(1e-3), "mse")
    return m

for tipo in ("CNN 1D", "LSTM", "GRU"):
    print(f"{tipo:<7}{red(tipo).count_params():>5} parámetros")

CNN 1D   529 parámetros
LSTM    1233 parámetros
GRU      977 parámetros

Las tres tienen entre 500 y 1.250 parámetros, frente a los 5 del HAR con apalancamiento por activo (30 en total). La LSTM tiene más que la GRU con el mismo estado: $4 \times (16 \times (2 + 16) + 16) = 1.216$ en la capa frente a 960.

El entrenamiento usa Adam (0,001), lotes de 512 y como mucho 30 épocas, con parada temprana a las 3 épocas sin mejora en la validación y los pesos de la mejor época. Los seis activos se entrenan juntos (un modelo de panel con las series tipificadas): las redes necesitan más datos que un HAR, y la tipificación por activo hace las series comparables.

historia = {}
for tipo in ("CNN 1D", "LSTM", "GRU"):
    keras.utils.set_random_seed(0)
    m = red(tipo)
    parada = keras.callbacks.EarlyStopping(
        patience=3, restore_best_weights=True)
    h = m.fit(X[ent], y[ent], epochs=30, batch_size=512, verbose=0,
              validation_data=(X[val], y[val]), callbacks=[parada])
    historia[tipo] = h.history["val_loss"]
    pred[tipo] = m.predict(X, batch_size=4096, verbose=0).ravel()
    mejor = int(np.argmin(h.history["val_loss"])) + 1
    print(f"{tipo:<7} épocas {len(h.history['loss']):>2}, mejor "
          f"validación en la {mejor}")
tabla = pd.Series({k: mse(p, te) for k, p in pred.items()},
                  name="MSE prueba").sort_values()
tabla.round(4)

CNN 1D  épocas 27, mejor validación en la 24
LSTM    épocas 19, mejor validación en la 16
GRU     épocas 19, mejor validación en la 16
oráculo calibrado       0.1962
GRU                     0.2106
LSTM                    0.2107
HAR + apalancamiento    0.2114
HAR                     0.2137
CNN 1D                  0.2137
ingenua (RV de hoy)     0.3184
Name: MSE prueba, dtype: float64

La parada temprana actuó en las tres: la LSTM y la GRU se pararon en la época 19 con su mejor validación en la 16, la CNN en la 27. En prueba (la figura de la cabecera pone las siete previsiones en una escala):

• Las dos recurrentes quedan casi empatadas con el HAR con apalancamiento: 0,2106 la GRU, 0,2107 la LSTM, 0,2114 el HAR con apalancamiento.
• La CNN 1D, con 0,2137, iguala al HAR sin apalancamiento: no aprovecha la información de los rendimientos tan bien.
• Todas quedan lejos del oráculo, 0,1962. La distancia es la misma para la red y para el modelo lineal.

¿Es la diferencia mayor que el ruido?

Una diferencia de 0,0008 en el MSE no dice nada sin su error típico. La prueba de Diebold y Mariano (1995) compara dos previsiones con la serie de diferencias de sus pérdidas: si la media de $d_t = e_{a,t}^2 - e_{b,t}^2$ no se distingue de cero, las dos previsiones son igual de buenas. Las pérdidas diarias están autocorrelacionadas (la volatilidad tiene memoria y los errores también), así que el error típico de la media se calcula con Newey-West. Aquí, $d_t$ es la media de los seis activos de cada día:

def diebold_mariano(a, b, lags=10):
    """Media de la diferencia diaria de pérdidas (a - b), promediada
    entre activos, con error típico de Newey-West."""
    e = pd.DataFrame({"f": fch[te],
                      "d": ((pred[a][te] - y[te]) * SD[te]) ** 2
                      - ((pred[b][te] - y[te]) * SD[te]) ** 2})
    d = e.groupby("f")["d"].mean().to_numpy()
    r = sm.OLS(d, np.ones_like(d)).fit(
        cov_type="HAC", cov_kwds={"maxlags": lags})
    return r.params[0], r.bse[0], r.tvalues[0]

for a, b in [("GRU", "HAR + apalancamiento"),
             ("LSTM", "HAR + apalancamiento"),
             ("HAR + apalancamiento", "HAR"),
             ("oráculo calibrado", "HAR + apalancamiento")]:
    dif, et, t = diebold_mariano(a, b)
    print(f"{a:<21}- {b:<21}{dif:+.4f} ± {et:.4f}  t = {t:+.1f}")

GRU                  - HAR + apalancamiento -0.0008 ± 0.0010  t = -0.8
LSTM                 - HAR + apalancamiento -0.0007 ± 0.0009  t = -0.8
HAR + apalancamiento - HAR                  -0.0023 ± 0.0005  t = -4.6
oráculo calibrado    - HAR + apalancamiento -0.0152 ± 0.0013  t = -11.6

La GRU y la LSTM mejoran al HAR con apalancamiento en 0,0008 y 0,0007, con un error típico de 0,0010: $t = -0{,}8$. No se distinguen. La prueba tiene potencia de sobra para ver diferencias reales: el término de apalancamiento mejora al HAR con $t = -4{,}6$, y el oráculo al HAR con apalancamiento con $t = -11{,}6$. Lo que pasa es que la red y el modelo lineal aprenden lo mismo. En estos datos era previsible: el proceso verdadero es un HAR con apalancamiento, y una red solo puede ganar en lo que el modelo lineal no captura (los saltos, la no linealidad del apalancamiento, las interacciones entre activos).

Por activo, la media esconde comportamientos opuestos:

por = pd.DataFrame({k: [mse(pred[k], te & (act == i))
                        for i in range(len(ACTIVOS))]
                    for k in ("HAR", "HAR + apalancamiento", "GRU",
                              "oráculo calibrado")},
                   index=ACTIVOS)
por["GRU / HAR+ap."] = por["GRU"] / por["HAR + apalancamiento"]
por.round(3)

          HAR  HAR + apalancamiento    GRU  oráculo calibrado  GRU / HAR+ap.
accion  0.309                 0.308  0.296              0.283          0.962
bono    0.134                 0.133  0.136              0.130          1.018
cripto  0.398                 0.397  0.385              0.354          0.969
divisa  0.121                 0.121  0.121              0.115          1.000
indice  0.174                 0.162  0.171              0.157          1.056
oro     0.147                 0.147  0.155              0.139          1.052

La GRU gana en la acción (cociente 0,962) y en la cripto (0,969), los dos activos con saltos en el proceso verdadero: tras un salto, la RV del día se dispara sin que la volatilidad de fondo cambie tanto, y una red que ve la forma de la ventana puede descontarlo mejor que un HAR lineal, que trata un salto como volatilidad persistente. Pierde en el índice (1,056) y el oro (1,052): en el índice el apalancamiento es el más fuerte de los seis ($\gamma = 0{,}18$), y un modelo de panel con una sola red para los seis activos lo recoge peor que un HAR con su propio coeficiente.

Esto sugiere el siguiente paso razonable, que no es una red más grande: un HAR por activo con un término para los saltos (la variación bipotencia, bv_5min, separa la parte continua de la de saltos) o una red pequeña que corrija los residuos del HAR. Y recuerda la regla del módulo: la red se compara con la mejor baseline lineal que se pueda construir, no con la más débil.

Secuencias más allá de la volatilidad

Las mismas piezas sirven para otras secuencias del curso, con las mismas trampas:

Y para rendimientos, no para volatilidad, recuerda el control del README: la media condicional de rendimiento es constante en estos datos, así que ninguna red puede predecirla. Una LSTM que «acierta» la dirección de mañana en vol5_diario está midiendo ruido o tiene una fuga.

Ejercicios

1. Entrena la GRU solo con la primera columna de la ventana (el log RV, sin los rendimientos). ¿Cuánto empeora en prueba? ¿Qué parte del modelo verdadero explica la diferencia?

keras.utils.set_random_seed(0)
x_in = keras.Input((L, 1))
m1 = keras.Model(x_in, layers.Dense(1)(layers.GRU(16)(x_in)))
m1.compile(keras.optimizers.Adam(1e-3), "mse")
m1.fit(X[ent][:, :, :1], y[ent], epochs=30, batch_size=512, verbose=0,
       validation_data=(X[val][:, :, :1], y[val]),
       callbacks=[keras.callbacks.EarlyStopping(
           patience=3, restore_best_weights=True)])
pred["GRU sin rendimientos"] = m1.predict(
    X[:, :, :1], batch_size=4096, verbose=0).ravel()
for k in ("HAR", "GRU sin rendimientos", "GRU"):
    print(f"{k:<21}{mse(pred[k], te):.4f}")

HAR                  0.2137
GRU sin rendimientos 0.2117
GRU                  0.2106

2. Calcula el R² de la regresión del log RV de prueba sobre la previsión (el cuadrado de la correlación) para el HAR con apalancamiento, la GRU y el oráculo, activo a activo, y compáralo con la tabla del README de los datos.

def r2(k, i):
    m = te & (act == i)
    return np.corrcoef(pred[k][m], y[m])[0, 1] ** 2

pd.DataFrame({k: [r2(k, i) for i in range(len(ACTIVOS))]
              for k in ("HAR + apalancamiento", "GRU",
                        "oráculo calibrado")},
             index=ACTIVOS).round(2)

        HAR + apalancamiento   GRU  oráculo calibrado
accion                  0.36  0.39               0.41
bono                    0.50  0.49               0.51
cripto                  0.40  0.42               0.46
divisa                  0.45  0.45               0.48
indice                  0.60  0.57               0.61
oro                     0.57  0.56               0.60

3. Un compañero mejora la GRU cambiando layers.GRU(16) por layers.Bidirectional(layers.GRU(16)) y presenta una bajada del MSE de prueba. Con las ventanas de esta lección, ¿puede ser una fuga? ¿Y si el modelo predijera la volatilidad de cada día de la ventana (secuencia a secuencia)?