Deriva de datos y de concepto: detectarla y qué hacer.

Dos derivas y solo una duele

Un modelo supervisado aprende la relación entre unas features $X$ y un objetivo $y$ a partir de una muestra. Lo que aprende es una aproximación a $P(y \mid X)$, y lo que usa para estimarla son ejemplos sacados de $P(X)$. Cualquiera de las dos puede cambiar cuando el modelo ya está en producción, y la literatura las separa (Gama et al., 2014; Quiñonero-Candela et al., 2009):

La deriva de concepto, a su vez, puede ser brusca (un cambio de normativa, un competidor que arbitra la señal), gradual (una anomalía que se publica y se va diluyendo) o recurrente (los regímenes de calma y estrés que vuelven).

flowchart TB
  D[Deriva] --> X["Datos: P(X)"]
  D --> C["Concepto: P(y|X)"]
  C --> B[Brusca]
  C --> G[Gradual]
  C --> R[Recurrente]

La deriva de datos es la que se ve: basta comparar la distribución de las features de hoy con la del entrenamiento, y eso se puede hacer cada día, sin esperar a ninguna etiqueta. La deriva de concepto es la que duele, y no se ve mirando las features: solo aparece en los resultados, cuando ya han llegado las etiquetas que el modelo predijo mal. En datos reales nunca se sabe cuál de las dos ha ocurrido ni cuándo. Los datos deriva5_* del nivel sí lo dicen.

El laboratorio: seis segmentos con fecha conocida

deriva5_muestras.parquet tiene 30 activos ficticios en 3.018 sesiones de la NYSE (2014-2025), seis features tipificadas y el rendimiento del día siguiente ya alineado. El objetivo es lineal en las features con una relación señal/ruido de rendimientos diarios: un IC verdadero de 0,10, un $R^2$ del 1 %. El README del nivel publica los segmentos:

El protagonista es un modelo fijo: un OLS estimado con 2014-2016 y que nunca se reentrena. Su calidad se mide como en todo el nivel, con el IC diario: la correlación transversal entre su predicción y el rendimiento del día siguiente de los 30 activos (lección 5.1.3).

import numpy as np
import pandas as pd

d = pd.read_parquet("deriva5_muestras.parquet")
seg = pd.read_csv("deriva5_verdad.csv", parse_dates=["desde"])
F = [f"f{k}" for k in range(1, 7)]
fechas = pd.DatetimeIndex(np.sort(d["fecha"].unique()))
n, NA = len(fechas), d["activo"].nunique()
X = d[F].to_numpy(float).reshape(n, NA, 6)  # sesión × activo × f
Y = d["objetivo"].to_numpy(float).reshape(n, NA)
s = np.searchsorted(seg["desde"], fechas, side="right") - 1

base = s == 0                                # 2014-2016
X1 = np.column_stack([np.ones(base.sum() * NA),
                      X[base].reshape(-1, 6)])
beta = np.linalg.lstsq(X1, Y[base].ravel(), rcond=None)[0]
print("coeficientes del modelo fijo (pb):",
      np.round(beta[1:] * 1e4, 1))

def ic_diario(P, Y):
    """Correlación transversal de cada sesión entre P e Y."""
    P = P - P.mean(1, keepdims=True)
    Yc = Y - Y.mean(1, keepdims=True)
    return ((P * Yc).sum(1)
            / np.sqrt((P**2).sum(1) * (Yc**2).sum(1)))

ic = pd.Series(ic_diario(X @ beta[1:], Y), index=fechas)
tabla = seg[["segmento", "desde", "tipo"]].copy()
tabla["desde"] = tabla["desde"].dt.date
tabla["ic_fijo"] = ic.groupby(s).mean().round(3).to_numpy()
tabla

coeficientes del modelo fijo (pb): [ 8.9  6.1 -4.5 -0.  -0.   1. ]
   segmento       desde               tipo  ic_fijo
0         0  2014-01-02               base    0.120
1         1  2017-01-03              datos    0.117
2         2  2019-07-01           concepto   -0.019
3         3  2021-01-04  datos_irrelevante   -0.026
4         4  2022-07-01   concepto_gradual   -0.035
5         5  2023-07-03            estable   -0.051

El OLS recupera la relación base con su error de muestreo (8,9, 6,1 y −4,5 pb frente a 8, 5 y −4) y ceros para las features sin efecto. Su IC medio no se mueve con la deriva de datos de 2017 (0,117 frente a 0,120) y se hunde con la de concepto de 2019 (−0,019). Desde entonces no se recupera: en 2021 la deriva de f4 es irrelevante para él, y la gradual de 2022 le quita además la parte de f2 (−0,051 al final).

Deriva de datos: PSI y KS

El índice de estabilidad de la población (PSI, Population Stability Index) viene del riesgo de crédito. Se cortan los datos de referencia en $B$ tramos (aquí, deciles), se mide qué fracción $p_i$ de la referencia y $q_i$ del periodo actual cae en cada uno y se suma:

Es una divergencia de Kullback-Leibler simetrizada sobre el histograma: cero si los dos histogramas coinciden, y crece con cualquier diferencia, de media, de dispersión o de forma. La regla práctica habitual lee un PSI por debajo de 0,1 como estable y por encima de 0,25 como deriva fuerte.

def psi(ref, cur, nb=10):
    """Population Stability Index con deciles de la referencia."""
    cortes = np.quantile(ref, np.linspace(0, 1, nb + 1)[1:-1])
    p = np.bincount(np.searchsorted(cortes, ref), minlength=nb)
    q = np.bincount(np.searchsorted(cortes, cur), minlength=nb)
    p = np.maximum(p / len(ref), 1e-6)
    q = np.maximum(q / len(cur), 1e-6)
    return float(np.sum((q - p) * np.log(q / p)))

ref = X[base]
filas = {}
for k in range(1, 6):
    cur = X[s == k]
    filas[f"{k} {seg['tipo'][k]}"] = [
        psi(ref[..., j].ravel(), cur[..., j].ravel())
        for j in range(6)]
pd.DataFrame(filas, index=F).T.round(2)

                       f1   f2   f3    f4   f5   f6
1 datos              0.97  0.0  0.5  0.00  0.0  0.0
2 concepto           0.00  0.0  0.0  0.00  0.0  0.0
3 datos_irrelevante  0.00  0.0  0.0  1.23  0.0  0.0
4 concepto_gradual   0.00  0.0  0.0  1.27  0.0  0.0
5 estable            0.00  0.0  0.0  1.28  0.0  0.0

El PSI ve exactamente las derivas de datos plantadas: 0,97 en f1 y 0,50 en f3 en el segmento 1 (la dispersión duplicada de f3 también cuenta), 1,23 en f4 desde 2021, y ceros en todo lo demás. Sobre todo, ceros en el segmento 2, el de la deriva de concepto, y en el gradual. La figura de la cabecera pone las dos cosas juntas: el PSI salta en 2017 y en 2021 sin que el IC se entere, y el IC se hunde en julio de 2019 sin que ningún PSI se mueva.

La alternativa clásica es la prueba de Kolmogorov-Smirnov de dos muestras (scipy.stats.ks_2samp), que compara las funciones de distribución empíricas y da un p-valor. Antes de usar cualquiera de las dos como alarma hay que saber qué dan cuando no pasa nada: se compara cada mes de 2016 con 2014-2015, sin ningún cambio plantado.

from scipy import stats

ref = X[fechas.year <= 2015].reshape(-1, 6)   # 2014-2015
meses = pd.period_range("2016-01", "2016-12", freq="M")
psis, pvals = [], []
for m in meses:
    cur = X[fechas.to_period("M") == m].reshape(-1, 6)
    for j in range(6):
        psis.append(psi(ref[:, j], cur[:, j]))
        pvals.append(stats.ks_2samp(ref[:, j], cur[:, j]).pvalue)
psis, pvals = np.array(psis), np.array(pvals)
print(f"72 pruebas (12 meses × 6 features) sin deriva")
print(f"PSI: mediana {np.median(psis):.3f}, máximo {psis.max():.3f}")
print(f"KS con p < 0,05: {np.mean(pvals < 0.05):.0%} "
      f"(lo esperado: 5 %)")
ks_ano = [stats.ks_2samp(ref[:, j],
                         X[fechas.year == 2016][..., j].ravel())
          for j in range(6)]
print("2016 entero, p del KS:",
      [f"{r.pvalue:.3f}" for r in ks_ano])
print("estadístico D:", [round(float(r.statistic), 3)
                         for r in ks_ano])

72 pruebas (12 meses × 6 features) sin deriva
PSI: mediana 0.018, máximo 0.055
KS con p < 0,05: 33% (lo esperado: 5 %)
2016 entero, p del KS: ['0.000', '0.146', '0.009', '0.000', '0.315', '0.018']
estadístico D: [0.029, 0.016, 0.023, 0.04, 0.013, 0.022]

Sin deriva, el PSI mensual se queda en valores de ruido (mediana 0,018, máximo 0,055), muy por debajo del 0,25. El KS no: el 33 % de las 72 pruebas sale con $p < 0{,}05$, seis veces lo esperado, y con 2016 entero frente a 2014-2015 cuatro de las seis features «cambian» con $p < 0{,}05$ (tres con $p \le 0{,}009$). Hay dos razones. Las features son AR(1) con $\rho = 0{,}5$ dentro de cada activo: las filas no son independientes, el tamaño efectivo de la muestra es varias veces menor que el número de filas y el p-valor sale demasiado pequeño. Y con miles de filas, diferencias minúsculas son significativas: los estadísticos D están entre 0,013 y 0,040: diferencias de unas centésimas en la función de distribución.

Deriva multivariante: validación adversaria

El PSI y el KS miran cada feature por separado. Una deriva en la relación entre features (dos que dejan de estar correlacionadas) no se ve así. La validación adversaria mira todas a la vez: se entrena un clasificador que distinga los datos de referencia de los actuales. Si no puede (AUC ≈ 0,5), las dos muestras son indistinguibles; si puede, la importancia del clasificador dice qué features han cambiado. Rabanser et al. (2019) comparan varios detectores de deriva en varios conjuntos de datos: el mejor en su estudio es un contraste de dos muestras sobre una representación reducida, y los clasificadores que separan dominios resultan útiles para describir la deriva y ver si es dañina.

Dos detalles deciden si el AUC significa algo. Los pliegues deben ser de meses enteros (filas de días vecinos son casi iguales, como en la lección 5.4.1) y el AUC se calcula en cada pliegue y se promedia. Si se juntan las predicciones de todos los pliegues y se calcula un solo AUC, sale por debajo de 0,5 sin deriva: cada pliegue tiene una proporción de clases distinta y el clasificador, sin nada que aprender, predice esa proporción, al revés en el pliegue de prueba.

import lightgbm as lgb
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import StratifiedGroupKFold

mref = fechas.year == 2016                 # el último año visto
mes = np.repeat((fechas.year * 100 + fechas.month).to_numpy(), NA)
mes = mes.reshape(n, NA)

def adversaria(m):
    """AUC de un clasificador que separa 2016 de las sesiones m."""
    A = np.vstack([X[mref].reshape(-1, 6), X[m].reshape(-1, 6)])
    c = np.r_[np.zeros(mref.sum() * NA), np.ones(m.sum() * NA)]
    g = np.r_[mes[mref].ravel(), mes[m].ravel()]
    aucs, imp = [], np.zeros(6)
    cv = StratifiedGroupKFold(4, shuffle=True, random_state=0)
    for tr, te in cv.split(A, c, g):          # meses enteros
        clf = lgb.LGBMClassifier(
            n_estimators=60, num_leaves=8, min_child_samples=100,
            random_state=0, n_jobs=4, deterministic=True,
            force_row_wise=True, verbose=-1).fit(A[tr], c[tr])
        p = clf.predict_proba(A[te])[:, 1]
        aucs.append(roc_auc_score(c[te], p))  # AUC de cada pliegue
        imp += clf.booster_.feature_importance("gain")
    return np.mean(aucs), F[int(imp.argmax())], imp.max() / imp.sum()

semestre = fechas.year.astype(str) + np.where(
    fechas.month <= 6, "-S1", "-S2")
adv = {}
for sm in semestre.unique():
    if sm.startswith("2016"):
        continue
    auc, top, peso = adversaria(semestre == sm)
    adv[sm] = {"auc": auc, "feature": top, "peso": peso,
               "ic_fijo": ic[semestre == sm].mean()}
adv = pd.DataFrame(adv).T.astype({"auc": float, "peso": float,
                                  "ic_fijo": float})
adv.loc[["2014-S1", "2015-S2", "2017-S1", "2019-S1", "2019-S2",
         "2020-S2", "2021-S1", "2023-S1", "2025-S2"]].round(3)

           auc feature   peso  ic_fijo
2014-S1  0.493      f3  0.203    0.138
2015-S2  0.518      f4  0.198    0.103
2017-S1  0.826      f1  0.535    0.142
2019-S1  0.824      f1  0.533    0.123
2019-S2  0.519      f2  0.198   -0.031
2020-S2  0.496      f1  0.185   -0.025
2021-S1  0.782      f4  0.905   -0.029
2023-S1  0.780      f4  0.916   -0.026
2025-S2  0.790      f4  0.919   -0.027

En los semestres sin deriva el AUC queda entre 0,49 y 0,52, lo que dan dos muestras iguales. En 2017-2019 sube a 0,83, con f1 como feature más importante (un 53 % de la ganancia). Desde 2021 queda en 0,78-0,79, con f4 (92 %). Y en 2019-S2 y 2020-S2, tras el cambio de concepto, vuelve a 0,52 y 0,50 mientras el IC del modelo fijo es negativo (−0,031 y −0,025). La validación adversaria responde muy bien a su pregunta, si los datos de hoy se parecen a los del entrenamiento, y esa pregunta no es la que importa para la cuenta de resultados.

Para qué sirve, entonces, vigilar las features: para diagnosticar (cuando el rendimiento cae, saber si las features cambiaron ayuda a entender por qué), para avisar de errores de datos (un proveedor que cambia una unidad, una columna que llega vacía) y para vigilar mientras las etiquetas no llegan. Con un objetivo a 21 días, el rendimiento de hoy se conoce dentro de un mes; las features, hoy.

Vigilar el rendimiento: media móvil, CUSUM y Page-Hinkley

La deriva de concepto solo se ve en el rendimiento. El IC diario es una serie ruidosa: con 30 activos y un IC verdadero de 0,12, su desviación típica es 0,18 (casi $1/\sqrt{29}$), así que un día con IC negativo no dice nada. El detector más sencillo, el del README del nivel, es la media móvil de 63 sesiones del IC con una alarma cuando cae por debajo de cero.

cambio = pd.Timestamp("2019-07-01")
mm = ic.rolling(63).mean()
fecha_alarma = mm[(mm.index >= cambio) & (mm < 0)].index[0]
retraso = ((fechas >= cambio) & (fechas <= fecha_alarma)).sum()
print("primera media móvil < 0:", fecha_alarma.date(),
      f"({retraso} sesiones después del cambio)")
antes = mm["2017-01-01":"2019-06-28"]
print(f"mínimo antes del cambio: {antes.min():.3f}")

primera media móvil < 0: 2019-08-29 (43 sesiones después del cambio)
mínimo antes del cambio: 0.071

Avisa el 2019-08-29, 43 sesiones después del cambio, sin falsas alarmas antes (la media móvil no bajó de 0,071 en 2017-2019). Es fácil de explicar y lenta: una ventana de 63 días necesita que la mitad de su contenido sea posterior al cambio para que la media cruce el cero.

El CUSUM de Page (1954) es el detector secuencial clásico. Para una caída de la media desde $\mu_0$:

$S_t$ acumula cuánto queda el IC por debajo de $\mu_0 - k$ y se reinicia a cero cuando el modelo va bien. La holgura $k$ se suele fijar en la mitad del salto que se quiere detectar (aquí, de 0,12 a 0: $k = 0{,}06$) y el umbral $h$ decide el equilibrio entre falsas alarmas y retraso. Es la misma recursión que el filtro CUSUM de la lección 5.3.5, que allí muestreaba precios y aquí vigila un modelo. Page-Hinkley es su variante sin $\mu_0$: compara cada dato con la media desde el último reinicio.

x = ic.to_numpy()
t0 = np.searchsorted(fechas, pd.Timestamp("2017-01-03"))
mu0 = ic[base].mean()

def cusum(x, mu0, k, h):
    """CUSUM unilateral de Page para una caída de la media."""
    S, alarmas = 0.0, []
    for t, v in enumerate(x):
        S = max(0.0, S + (mu0 - k) - v)
        if S > h:
            alarmas.append(t)
            S = 0.0                       # reinicio tras la alarma
    return alarmas

def page_hinkley(x, delta, lam):
    """Page-Hinkley: caída frente a la media desde el reinicio."""
    m = M = suma = 0.0
    cnt, alarmas = 0, []
    for t, v in enumerate(x):
        cnt += 1
        suma += v
        m += suma / cnt - v - delta
        M = min(M, m)
        if m - M > lam:
            alarmas.append(t)
            m = M = suma = 0.0
            cnt = 0
    return alarmas

def resumen(alarmas):
    f = fechas[t0:][alarmas]
    falsas = int((f < cambio).sum())
    tras = f[f >= cambio]
    r = int(((fechas >= cambio) & (fechas <= tras[0])).sum())
    return falsas, tras[0].date(), r

print(f"IC medio de referencia (2014-2016): {mu0:.3f}")
for h in [0.5, 1.0, 2.0, 3.0]:
    fa, f, r = resumen(cusum(x[t0:], mu0, 0.06, h))
    print(f"CUSUM h={h}: {fa:2d} falsas alarmas 2017-2019, "
          f"primera tras el cambio {f} ({r} sesiones)")
for lam in [1.0, 2.0, 3.0]:
    fa, f, r = resumen(page_hinkley(x[t0:], 0.03, lam))
    print(f"Page-Hinkley lambda={lam}: {fa:2d} falsas alarmas, "
          f"{f} ({r} sesiones)")

IC medio de referencia (2014-2016): 0.120
CUSUM h=0.5: 19 falsas alarmas 2017-2019, primera tras el cambio 2019-07-08 (5 sesiones)
CUSUM h=1.0:  3 falsas alarmas 2017-2019, primera tras el cambio 2019-07-16 (11 sesiones)
CUSUM h=2.0:  0 falsas alarmas 2017-2019, primera tras el cambio 2019-07-19 (14 sesiones)
CUSUM h=3.0:  0 falsas alarmas 2017-2019, primera tras el cambio 2019-08-01 (23 sesiones)
Page-Hinkley lambda=1.0:  8 falsas alarmas, 2019-07-11 (8 sesiones)
Page-Hinkley lambda=2.0:  1 falsas alarmas, 2019-07-18 (13 sesiones)
Page-Hinkley lambda=3.0:  0 falsas alarmas, 2019-07-29 (20 sesiones)

El CUSUM con $h = 2$ avisa el 2019-07-19, a las 14 sesiones, sin ninguna falsa alarma en los dos años y medio anteriores, que tienen además la deriva de datos de 2017. Con $h = 0{,}5$ avisa a las 5 sesiones, a cambio de 19 falsas alarmas en ese periodo, unas siete al año. Page-Hinkley se comporta parecido: con $\lambda = 3$, 20 sesiones sin falsas alarmas. La mitad de arriba de la figura sigue el estadístico del CUSUM (sin reinicio) en 2019: plano hasta julio y, desde el cambio, una subida casi recta.

Falsas alarmas frente a retraso

Un detector no se juzga por un episodio. El cambio de 2019 es una sola realización, y 14 sesiones pueden ser suerte. Para medir el equilibrio se simulan muchas historias: 2.000 caminos de diez años con IC sacados al azar de los días sin cambio (2017-2019-S1) para contar falsas alarmas, y caminos de un año con IC de después del cambio (2019-S2 a 2020) para medir el retraso.

tc = np.searchsorted(fechas, cambio)
t21 = np.searchsorted(fechas, pd.Timestamp("2021-01-01"))
sano, roto = x[t0:tc], x[tc:t21]     # IC antes y después
rng = np.random.default_rng(0)
NS, T = 2000, 2520                   # 2.000 caminos de 10 años
a = rng.choice(sano, (NS, T))
b = rng.choice(roto, (NS, 252))

def sim_cusum(h, k=0.06):
    S, fa = np.zeros(NS), np.zeros(NS)
    for t in range(T):
        S = np.maximum(0, S + (mu0 - k) - a[:, t])
        fa += S > h
        S[S > h] = 0
    S, r = np.zeros(NS), np.full(NS, 252)
    for t in range(252):
        S = np.maximum(0, S + (mu0 - k) - b[:, t])
        r[(S > h) & (r == 252)] = t + 1
    return fa.mean() / 10, np.median(r)

def sim_media(w):
    c = np.cumsum(a, 1)
    bajo = (c[:, w:] - c[:, :-w]) / w < 0
    fa = (bajo[:, 1:] & ~bajo[:, :-1]).sum(1).mean() / 10
    c = np.cumsum(np.hstack([np.zeros((NS, 1)), a[:, -w:], b]), 1)
    tras = ((c[:, w:] - c[:, :-w]) / w)[:, 1:] < 0
    r = np.where(tras.any(1), tras.argmax(1) + 1, 252)
    return fa, np.median(r)

filas = [("CUSUM", h, *sim_cusum(h)) for h in [1, 1.5, 2, 2.5, 3]]
filas += [("media móvil", w, *sim_media(w)) for w in [21, 42, 63]]
pd.DataFrame(filas, columns=["detector", "parámetro",
                             "falsas/año", "retraso mediano"]
             ).round(3)

      detector  parámetro  falsas/año  retraso mediano
0        CUSUM        1.0       1.030             11.0
1        CUSUM        1.5       0.193             17.0
2        CUSUM        2.0       0.038             23.0
3        CUSUM        2.5       0.007             30.0
4        CUSUM        3.0       0.001             36.0
5  media móvil       21.0       0.204             18.0
6  media móvil       42.0       0.002             36.0
7  media móvil       63.0       0.000             54.0

El retraso típico del CUSUM con $h = 2$ es de 23 sesiones, no 14, con 0,04 falsas alarmas al año (una cada 25 años). La media móvil de 63 sesiones tarda 54 de mediana: los 43 de los datos también fueron algo afortunados. La figura pone las dos familias en el mismo plano. A igual tasa de falsas alarmas, el CUSUM llega unas sesiones antes: con 0,2 al año, 17 sesiones frente a 18; con unas 0,002-0,007, 30 frente a 36. La ventaja existe (el CUSUM es óptimo para este problema en el sentido de Lorden), pero es modesta. La decisión importante es dónde ponerse en la curva: cuánto cuesta una falsa alarma (apagar un modelo que funciona) frente a cuánto cuesta cada día de retraso (operar con un modelo roto).

Regímenes: el filtro de dos estados

La lección 2.1.7 detectaba regímenes con un modelo de Markov oculto y distinguía entre probabilidades filtradas (con la información hasta cada día) y suavizadas (con toda la muestra, también la futura). La misma idea vigila un modelo: dos estados, «funciona» (IC medio $\mu_0 = 0{,}12$) y «no funciona» (IC medio 0), con la misma desviación, una probabilidad de cambiar de estado de 0,002 al día y la recursión del filtro:

donde $A$ es la matriz de transición y $\phi(x_t)$ la densidad del IC de hoy en cada estado. Los parámetros se fijan de antemano, con 2014-2016, en vez de estimarlos con toda la muestra: estimarlos con datos de 2019-2025 sería mirar el futuro aunque se use el filtro.

from scipy.stats import norm

medias = np.array([mu0, 0.0])     # el modelo funciona / no funciona
sd = ic[base].std()
p = 0.998                         # duración esperada: 500 sesiones
A = np.array([[p, 1 - p], [1 - p, p]])
pi, prob = np.array([1.0, 0.0]), []
for v in x[t0:]:
    pi = (pi @ A) * norm.pdf(v, medias, sd)  # predecir y actualizar
    pi = pi / pi.sum()
    prob.append(pi[1])
prob = pd.Series(prob, index=fechas[t0:])
print(f"desv. del IC diario: {sd:.3f}")
print("P(roto) máxima antes del cambio:",
      round(prob[:"2019-06-28"].max(), 2))
print("primera P(roto) > 0,9:",
      prob[(prob.index >= cambio) & (prob > 0.9)].index[0].date())
print(prob.resample("QE").mean()["2019":"2020"].round(2)
      .to_string())

desv. del IC diario: 0.184
P(roto) máxima antes del cambio: 0.63
primera P(roto) > 0,9: 2019-07-19
2019-03-31    0.02
2019-06-30    0.01
2019-09-30    0.83
2019-12-31    0.98
2020-03-31    0.98
2020-06-30    0.98
2020-09-30    0.97
2020-12-31    0.99
Freq: QE-DEC

La probabilidad de que el modelo esté roto pasa de 0,9 el 2019-07-19, el mismo día que el CUSUM. No es casualidad: con dos estados gaussianos de medias conocidas, el logaritmo del cociente de verosimilitudes de cada día es lineal en $x_t$ con pendiente proporcional a la diferencia de medias, y el CUSUM con $k$ en el punto medio acumula exactamente ese cociente. El filtro añade la probabilidad previa de cambio y da una cifra que se lee directamente. Su falsa alarma más cercana antes del cambio llegó a 0,63, la mitad de abajo de la figura.

Los regímenes recurrentes piden otra respuesta. En el panel del nivel el momentum se invierte en cada rebote (lección 5.1.1): un detector que avise cuando el rebote ya ha empezado llega tarde a un régimen que dura unos meses, y reentrenar con datos del rebote enseña al modelo una relación que se va a invertir otra vez. Si los regímenes vuelven, es mejor que el modelo los vea como una feature observable y aprenda la interacción (lección 5.10.1); la detección de deriva sirve para lo que no vuelve.

Qué hacer: reentrenar, reponderar, apagar

Una alarma no dice qué hacer. Primero se diagnostica: si las features también cambiaron (PSI, validación adversaria), puede ser deriva de datos con efecto real en el modelo, un error de datos o las dos cosas; si no, es deriva de concepto.

flowchart TB
  M[Monitor del IC] --> A{¿Alarma?}
  A -- no --> M
  A -- sí --> D{¿PSI o AUC altos?}
  D -- no --> C[Concepto]
  D -- sí --> E[Datos o errores]
  C --> R[Apagar, reentrenar]
  E --> R

Las respuestas posibles, comparadas en un walk-forward que reentrena el primer día de cada mes con los objetivos ya conocidos:

• Reentrenar con ventana expansiva: todo el pasado.
• Reentrenar con ventana móvil: solo el último año (252 sesiones) o medio año (126).
• Reponderar: todo el pasado con pesos que caen a la mitad cada 126 sesiones.
• El oráculo: los coeficientes verdaderos de cada día, el techo.

def ols(i0, i1, semivida=None):
    """OLS con las sesiones i0..i1-1 (pesos exponenciales)."""
    A = np.column_stack([np.ones((i1 - i0) * NA),
                         X[i0:i1].reshape(-1, 6)])
    y = Y[i0:i1].ravel()
    if semivida:
        edad = np.repeat(np.arange(i1 - i0)[::-1], NA)
        w = np.sqrt(0.5 ** (edad / semivida))
        A, y = A * w[:, None], y * w
    return np.linalg.lstsq(A, y, rcond=None)[0][1:]

mes_n = (fechas.year * 12 + fechas.month).to_numpy()
reentreno = {t for t in range(t0, n) if mes_n[t] != mes_n[t - 1]}

def walk_forward(entrenar):
    """Reentrena el primer día de cada mes con lo que ya se sabe."""
    P, b = np.zeros((n, NA)), None
    for t in range(t0, n):
        if t in reentreno:
            b = entrenar(t - 1)      # objetivos conocidos hasta t
        P[t] = X[t] @ b
    return P

b_verdad = pd.read_csv("deriva5_coef_verdad.csv").iloc[:, 1:]
esquemas = {
    "fijo (2014-2016)": lambda t: beta[1:],
    "expansiva": lambda t: ols(0, t),
    "móvil 252": lambda t: ols(t - 252, t),
    "móvil 126": lambda t: ols(t - 126, t),
    "pesos, semivida 126": lambda t: ols(0, t, 126),
}
pred = {k: walk_forward(f) for k, f in esquemas.items()}
pred["oráculo (b verdaderos)"] = np.einsum(
    "tik,tk->ti", X, b_verdad.to_numpy())
res = {}
for k, P in pred.items():
    v = pd.Series(ic_diario(P[t0:], Y[t0:]), index=fechas[t0:])
    res[k] = {"2017-2025": v.mean(),
              "2019-S2 a 2020": v["2019-07":"2020-12"].mean(),
              "2022-S2 a 2023-S1": v["2022-07":"2023-06"].mean()}
pd.DataFrame(res).T.round(3)

                        2017-2025  2019-S2 a 2020  2022-S2 a 2023-S1
fijo (2014-2016)            0.007          -0.019             -0.035
expansiva                   0.052           0.001              0.027
móvil 252                   0.101           0.055              0.105
móvil 126                   0.101           0.066              0.098
pesos, semivida 126         0.101           0.054              0.104
oráculo (b verdaderos)      0.111           0.091              0.108

En todo el periodo 2017-2025 el modelo fijo apenas suma (IC 0,007), la expansiva llega a 0,052 y las tres variantes que olvidan el pasado, a 0,101, cerca del oráculo (0,111). La diferencia está en los 18 meses que siguen al cambio de concepto: la expansiva da 0,001, porque hasta 2019 acumula cinco años de un efecto de f1 positivo; la móvil de 126 sesiones da 0,066 y la de 252, 0,055. En el periodo gradual de 2022-2023 todas las adaptativas siguen al oráculo (0,098-0,105 frente a 0,108): un cambio repartido en doce meses es fácil de seguir con una ventana de uno.

Es lo contrario de lo que la lección 5.4.3 encontró en el panel semanal, donde ganaba la expansiva, y no hay contradicción. Allí los efectos son pequeños y oscilan sin tendencia; aquí el cambio es grande frente al ruido y permanente. La ventana óptima depende de eso (Pesaran y Timmermann, 2007, citados en la lección 5.4.3).

La figura sigue el coeficiente de f1 estimado cada mes. La móvil de 126 sesiones cruza el cero a finales de 2019 y la de 252, en enero de 2020; después oscilan alrededor de −8 con su error de estimación. La expansiva ya bajaba antes del cambio (de 9 a 6,5 en 2017-2019: el desplazamiento de f1 en 2017, sin cambio en la media del objetivo, atenúa la pendiente de un OLS con una sola constante para los dos periodos) y no se vuelve negativa hasta noviembre de 2024.

La combinación con el detector: el modelo fijo hasta la alarma del CUSUM; después, apagado unas sesiones (sin posición) y reentreno solo con los datos posteriores a la alarma.

alarma = t0 + cusum(x[t0:], mu0, 0.06, 2.0)[0]
print("alarma del CUSUM (h = 2):", fechas[alarma].date())

def tras_alarma(espera):
    """Modelo fijo hasta la alarma; apagado `espera` sesiones;
    después, OLS solo con datos desde la alarma (cada mes)."""
    P, b = np.zeros((n, NA)), beta[1:]
    for t in range(t0, n):
        if alarma <= t < alarma + espera:
            b = np.zeros(6)
        elif t >= alarma + espera and (t in reentreno
                                       or t == alarma + espera):
            b = ols(alarma, t - 1)
        P[t] = X[t] @ b
    with np.errstate(invalid="ignore"):   # apagado: P = 0, IC 0/0
        v = ic_diario(P[t0:], Y[t0:])
    return pd.Series(np.nan_to_num(v), index=fechas[t0:])

for espera in [21, 63, 126]:
    v = tras_alarma(espera)
    print(f"apagado {espera:3d} sesiones: IC 2017-2025 "
          f"{v.mean():.3f} · 2019-S2 a 2020 "
          f"{v['2019-07':'2020-12'].mean():.3f}")

alarma del CUSUM (h = 2): 2019-07-19
apagado  21 sesiones: IC 2017-2025 0.103 · 2019-S2 a 2020 0.075
apagado  63 sesiones: IC 2017-2025 0.101 · 2019-S2 a 2020 0.063
apagado 126 sesiones: IC 2017-2025 0.100 · 2019-S2 a 2020 0.054

Con 21 sesiones apagado, el IC de los 18 meses posteriores al cambio es 0,075, mejor que cualquier ventana fija y a 0,016 del oráculo. Esperar más cuesta: con 126 sesiones apagado baja a 0,054, porque apagar también apaga f2 y f3, que seguían funcionando. Apagar es una medida de emergencia (no operar con un modelo del que se sospecha), no una solución: el valor está en volver pronto con un modelo estimado con datos del régimen nuevo.

movil = {fechas[t].date(): ols(t - 252, t) * 1e4
         for t in sorted(reentreno)}
movil = pd.DataFrame(movil, index=F).T
verdad = b_verdad.set_axis(F, axis=1).set_index(fechas) * 100
fechas_ver = ["2019-07-01", "2019-10-01", "2020-01-02",
              "2022-07-01", "2023-01-03", "2023-07-03", "2024-07-01"]
comp = pd.DataFrame({
    f: {"b1 móvil": movil.loc[pd.Timestamp(f).date(), "f1"],
        "b1 verdad": verdad.loc[f, "f1"],
        "b2 móvil": movil.loc[pd.Timestamp(f).date(), "f2"],
        "b2 verdad": verdad.loc[f, "f2"],
        "b5 móvil": movil.loc[pd.Timestamp(f).date(), "f5"],
        "b5 verdad": verdad.loc[f, "f5"]} for f in fechas_ver}).T
comp.round(1)

            b1 móvil  b1 verdad  b2 móvil  b2 verdad  b5 móvil  b5 verdad
2019-07-01       7.9       -8.0       4.0        5.0       0.9        0.0
2019-10-01       2.4       -8.0       4.9        5.0       1.5        0.0
2020-01-02      -0.5       -8.0       4.7        5.0       0.9        0.0
2022-07-01      -7.2       -8.0       2.6        5.0       1.3        0.0
2023-01-03      -7.9       -8.0       2.2        2.4       3.0        3.1
2023-07-03      -8.2       -8.0       1.8        0.0       4.3        6.0
2024-07-01     -10.1       -8.0      -1.0        0.0       6.6        6.0

La tabla muestra los límites de la ventana móvil de 252 sesiones con los coeficientes verdaderos. En julio de 2019 el modelo aún cree en el efecto viejo (7,9 frente a −8,0), en octubre va por 2,4 y en enero de 2020, por −0,5. En la deriva gradual, el coeficiente de f2 baja (2,6 en julio de 2022, 1,8 en julio de 2023, cuando la verdad ya es 0) y el de f5 sube (de 1,3 a 4,3, frente a 0 y 6): la ventana de un año va con unos seis meses de retraso, que es lo que cuesta promediar un año de datos ruidosos.

Ejercicios

1. La deriva de datos de 2017 es un covariate shift de manual. La corrección clásica es reponderar el entrenamiento por el cociente de densidades $P_{\text{actual}}(x)/P_{\text{base}}(x)$, que se estima con el clasificador adversario como $p/(1-p)$. Estima el OLS de 2014-2016 con esos pesos (actual = 2017-2019-S1) y compara coeficientes e IC con el modelo fijo sin pesos.

actual = (fechas >= "2017-01-01") & (fechas < "2019-07-01")
A = np.vstack([X[base].reshape(-1, 6), X[actual].reshape(-1, 6)])
c = np.r_[np.zeros(base.sum() * NA), np.ones(actual.sum() * NA)]
clf = lgb.LGBMClassifier(
    n_estimators=60, num_leaves=8, min_child_samples=100,
    random_state=0, n_jobs=4, deterministic=True,
    force_row_wise=True, verbose=-1).fit(A, c)
pb = clf.predict_proba(X[base].reshape(-1, 6))[:, 1]
w = pb / (1 - pb)
w = w / w.mean()
sw = np.sqrt(w)
beta_w = np.linalg.lstsq(X1 * sw[:, None], Y[base].ravel() * sw,
                         rcond=None)[0]
print("sin pesos (pb):", np.round(beta[1:] * 1e4, 1))
print("con pesos (pb):", np.round(beta_w[1:] * 1e4, 1))
print(f"tamaño efectivo: {w.sum()**2 / (w**2).sum() / w.size:.0%} "
      "de las filas")
for nombre, b in [("sin pesos", beta), ("con pesos", beta_w)]:
    v = ic_diario(X[actual] @ b[1:], Y[actual]).mean()
    print(f"IC 2017-2019-S1 {nombre}: {v:.3f}")

sin pesos (pb): [ 8.9  6.1 -4.5 -0.  -0.   1. ]
con pesos (pb): [ 5.5  4.4 -4.3  0.9 -2.   1.2]
tamaño efectivo: 13% de las filas
IC 2017-2019-S1 sin pesos: 0.117
IC 2017-2019-S1 con pesos: 0.112

2. El CUSUM vigilaba el modelo fijo. Vigila ahora el modelo de ventana móvil de 252 sesiones con el mismo detector ($k = 0{,}06$, $h = 2$, $\mu_0$ = su IC medio en 2017-2018). ¿Cuándo avisa y qué pasa en la deriva gradual?

ic_movil = pd.Series(ic_diario(pred["móvil 252"][t0:], Y[t0:]),
                     index=fechas[t0:])
mu_m = ic_movil[:"2018-12-31"].mean()
al = cusum(ic_movil.to_numpy(), mu_m, 0.06, 2.0)
print(f"IC medio 2017-2018: {mu_m:.3f}")
print("alarmas:", [str(f.date()) for f in ic_movil.index[al]])
print(ic_movil.groupby(ic_movil.index.year).mean().round(3)
      .to_string())

IC medio 2017-2018: 0.112
alarmas: ['2019-07-19', '2019-08-23', '2019-11-26']
2017    0.112
2018    0.112
2019    0.058
2020    0.084
2021    0.097
2022    0.105
2023    0.097
2024    0.133
2025    0.115

3. Diseña la vigilancia de un modelo semanal de selección de acciones cuya etiqueta es el rendimiento a 21 sesiones (solapada, como ret_sig_21d del panel). ¿Qué vigilas cada día, qué cada semana y cómo calibras los umbrales?