Detección de regímenes: modelos de Markov ocultos y cambios estructurales.

Regímenes: el mercado no tiene un solo comportamiento

Las lecciones 1.6.1 y 1.6.2 midieron que la volatilidad viene en rachas. Un paso más es suponer que el mercado alterna unos pocos regímenes (calma alcista, lateral nervioso, caída con pánico), cada uno con su media y su volatilidad, y que pasa de uno a otro con cierta persistencia. Si se supiera en qué régimen se está, se podría ajustar la exposición, el tamaño de las posiciones o la estrategia (lección 2.2.2). El problema es que el régimen no se ve: solo se ven los rendimientos.

Hay dos maneras de pensar en un cambio de comportamiento, y esta lección trata las dos:

• Regímenes recurrentes. El mercado entra y sale de los mismos estados muchas veces. El modelo natural es un modelo de Markov oculto (HMM, hidden Markov model), que en econometría popularizó Hamilton (1989).
• Ruptura estructural. Algo cambia una vez y no vuelve: una regulación, un cambio de régimen monetario, una empresa que cambia de negocio. Lo natural es un contraste de cambio estructural que busca la fecha.

Los datos son las dos series sintéticas de regimenes2.parquet (sesiones de la NYSE de 2000 a 2025). En markov, los rendimientos son normales con media y volatilidad decididas por una cadena de Markov de tres estados; en ruptura, hay un único cambio el 2013-05-22. En las dos, el estado verdadero de cada día está en la columna estado_real: la respuesta conocida contra la que se mide cada método. Son datos inventados: los parámetros salen del generador, no de ningún mercado.

La cadena de Markov oculta

Un HMM tiene dos piezas. La cadena oculta $s_t \in \{0, 1, 2\}$ pasa de un estado a otro con una matriz de transición $A$, con $A_{ij} = P(s_{t+1} = j \mid s_t = i)$. Las emisiones son los rendimientos: $r_t \mid s_t = k \sim N(\mu_k, \sigma_k^2)$. La matriz verdadera de la serie markov es

con una diagonal cercana a 1: los regímenes son persistentes. La duración media de una estancia en el estado $k$ es $1 / (1 - A_{kk})$: 125 sesiones en el alcista, 50 en el lateral y 33 en el bajista.

stateDiagram-v2
  direction LR
  Alcista --> Lateral: 0,006
  Lateral --> Alcista: 0,012
  Lateral --> Bajista: 0,008
  Bajista --> Lateral: 0,020
  Alcista --> Bajista: 0,002
  Bajista --> Alcista: 0,010

import logging
import warnings

import numpy as np
import pandas as pd

reg = pd.read_parquet("regimenes2.parquet")
mk = (reg[reg["serie"] == "markov"]
      .dropna(subset=["rendimiento_log"]).set_index("fecha"))
r = mk["rendimiento_log"]
verdad = mk["estado_real"].to_numpy()
nombres = ["alcista", "lateral", "bajista"]
resumen = mk.groupby("estado_real")["rendimiento_log"].agg(
    ["size", "mean", "std"])
resumen.index = nombres
resumen.round(5)

         size     mean      std
alcista  3479  0.00086  0.00700
lateral  2407 -0.00018  0.01417
bajista   652 -0.00115  0.02670

La figura de la cabecera enseña la serie con sus regímenes: las caídas del precio coinciden con los tramos bajistas, pero no todos los tramos agitados son caídas, y entre un régimen y otro no hay ninguna señal visible a simple vista. Las volatilidades medidas en los días de cada estado coinciden con las del generador (0,007, 0,014 y 0,028 diarias); las medias, solo de lejos (0,00086, −0,00018 y −0,00115, frente a 0,0006, 0 y −0,0018): con unos cientos de días por estado, la media muestral se aleja mucho de la verdadera, y la volatilidad casi nada.

Ajustar el HMM: varios arranques

hmmlearn estima un HMM gaussiano con el algoritmo EM (expectation-maximization): alterna calcular la probabilidad de cada estado en cada día con los parámetros actuales y reestimar los parámetros con esas probabilidades. EM sube la verosimilitud en cada paso, pero puede quedarse en un máximo local, y a cuál llega depende del punto de partida, que hmmlearn sortea con random_state. La práctica correcta es probar varios arranques:

from hmmlearn.hmm import GaussianHMM

# hmmlearn avisa por logging cuando una iteración de EM baja la
# verosimilitud en la sexta cifra («Model is not converging») y
# usa una operación que numpy 2.5 declara obsoleta
# (DeprecationWarning). Ninguno de los dos cambia el resultado.
logging.getLogger("hmmlearn").setLevel(logging.ERROR)
warnings.filterwarnings("ignore", category=DeprecationWarning,
                        module="hmmlearn")

X = r.to_numpy().reshape(-1, 1)

def ajustar(X, semilla, k=3):
    return GaussianHMM(n_components=k, covariance_type="full",
                       n_iter=200, random_state=semilla).fit(X)

def por_volatilidad(modelo):
    """Estados ordenados de menor a mayor volatilidad."""
    return np.argsort(modelo.covars_.ravel())

def acierto(modelo, X, verdad):
    mapa = np.empty(modelo.n_components, int)
    mapa[por_volatilidad(modelo)] = np.arange(modelo.n_components)
    return (mapa[modelo.predict(X)] == verdad).mean()

modelos = {s: ajustar(X, s) for s in range(8)}
semillas = pd.DataFrame(
    {"log_verosimilitud": {s: m.score(X) for s, m in modelos.items()},
     "acierto": {s: acierto(m, X, verdad) for s, m in modelos.items()}})
semillas.rename_axis("semilla").round(3)

         log_verosimilitud  acierto
semilla
0                20146.518    0.370
1                20145.782    0.381
2                20281.648    0.944
3                20151.671    0.863
4                20281.611    0.946
5                20283.232    0.938
6                20146.817    0.586
7                20283.230    0.938

Tres de los ocho arranques (semillas 0, 1 y 6) acaban en máximos locales malos, con log-verosimilitudes de unas 20.146 y aciertos del 37 % al 59 %: el modelo ha repartido los estados de otra manera (por ejemplo, ha partido el régimen tranquilo en dos y ha mezclado el lateral con el bajista). Otro (semilla 3) se queda a medio camino, con 20.152 y un 86 %. Los otros cuatro llegan a 20.282-20.283 y aciertan el 94 %. La regla práctica no necesita la verdad: quedarse con el arranque de mayor verosimilitud, que aquí es también uno de los que más aciertan.

El mejor modelo y el intercambio de etiquetas

mejor = modelos[semillas["log_verosimilitud"].idxmax()]
orden = por_volatilidad(mejor)
medias = mejor.means_.ravel()[orden]
vols = np.sqrt(mejor.covars_.ravel()[orden])
A = mejor.transmat_[np.ix_(orden, orden)]
pi0 = mejor.startprob_[orden]
print(pd.DataFrame({"media": medias, "media_real": resumen["mean"],
                    "vol": vols, "vol_real": resumen["std"]},
                   index=nombres).round(5))
pd.DataFrame(A, index=nombres, columns=nombres).round(3)

           media  media_real      vol  vol_real
alcista  0.00082     0.00086  0.00726   0.00700
lateral -0.00022    -0.00018  0.01438   0.01417
bajista -0.00088    -0.00115  0.02646   0.02670
         alcista  lateral  bajista
alcista    0.990    0.008    0.002
lateral    0.012    0.977    0.011
bajista    0.009    0.037    0.955

Las columnas _real son las medias y volatilidades de los días que de verdad estuvieron en cada estado. Las volatilidades se recuperan casi exactas; las medias, peor (la del bajista sale en −0,00088, frente a −0,00115 en sus días y −0,0018 en el generador), porque una media diaria se estima con mucho más error que una volatilidad: es el mismo problema que la ruptura de media de más abajo. La matriz de transición tiene la estructura correcta, con la persistencia del bajista algo infravalorada.

El paso por_volatilidad no es cosmético. Las etiquetas 0, 1, 2 que devuelve hmmlearn son arbitrarias: si se permutan los estados, la verosimilitud es la misma, así que nada obliga al estado 0 a ser el tranquilo. Es el intercambio de etiquetas (label switching). Con los cuatro arranques buenos, la etiqueta en bruto que recibe la mayoría de los días de cada régimen verdadero:

real = pd.Series(verdad, name="verdad").map(dict(enumerate(nombres)))
buenas = semillas.index[semillas["acierto"] > 0.9]
etiquetas = {
    f"semilla {s}": pd.crosstab(real, modelos[s].predict(X))
    .idxmax(axis=1) for s in buenas}
pd.DataFrame(etiquetas).reindex(nombres)

         semilla 2  semilla 4  semilla 5  semilla 7
verdad
alcista          1          0          0          1
lateral          0          2          1          0
bajista          2          1          2          2

Los cuatro modelos son prácticamente el mismo, pero numeran los regímenes de forma distinta: en uno el alcista es el 0 y en otro el 1. Comparar el estado en bruto con la verdad daría aciertos absurdos, y dos ajustes del mismo modelo (o el de ayer y el de hoy en un sistema que reestima) pueden cambiar el número de cada régimen. Ordenarlos por volatilidad, o por media, es lo que les da significado.

Filtrar o suavizar: lo que se sabía cada día

Con los parámetros estimados queda la pregunta que importa: ¿en qué régimen está el mercado hoy? Hay dos respuestas distintas:

• La probabilidad filtrada, $P(s_t = k \mid r_1, \dots, r_t)$, usa los rendimientos hasta hoy. Es lo que un sistema sabe en el momento de decidir.
• La probabilidad suavizada, $P(s_t = k \mid r_1, \dots, r_T)$, usa toda la muestra, también lo que pasó después de $t$. Es lo que devuelve predict_proba de hmmlearn (algoritmo hacia delante y hacia atrás), y el estado de Viterbi de predict también usa toda la muestra.

flowchart TB
  P[Rendimientos 1..T] --> F[Hacia delante]
  F --> PF["Filtrada: hasta t"]
  F --> B[Hacia atrás]
  B --> PS["Suavizada: hasta T"]
  PF --> U[Sirve para operar]
  PS --> N[Solo para estudiar]

La filtrada se programa en unas líneas con el algoritmo hacia delante: la probabilidad de ayer se proyecta con la matriz de transición y se corrige con la densidad del rendimiento de hoy en cada estado.

from scipy.stats import norm

def filtrar(x, pi0, A, medias, vols):
    """P(estado_t | rendimientos hasta t): algoritmo hacia delante."""
    dens = norm.pdf(np.asarray(x)[:, None], medias, vols)
    alfa = np.zeros_like(dens)
    prior = pi0
    for t in range(len(dens)):
        a = prior * dens[t]          # corrección con el dato de hoy
        alfa[t] = a / a.sum()
        prior = alfa[t] @ A          # predicción para mañana
    return alfa

filtrada = filtrar(r, pi0, A, medias, vols)
suavizada = mejor.predict_proba(X)[:, orden]
for nombre, p in [("filtrada", filtrada), ("suavizada", suavizada)]:
    print(f"{nombre:9}  acierto {(p.argmax(1) == verdad).mean():.3f}")
etiqueta = dict(enumerate(nombres))
pd.crosstab(real, pd.Series(filtrada.argmax(1), name="filtrada")
            .map(etiqueta)).reindex(index=nombres, columns=nombres)

filtrada   acierto 0.889
suavizada  acierto 0.953
filtrada  alcista  lateral  bajista
verdad
alcista      3281      185       13
lateral       257     2045      105
bajista         5      162      485

La filtrada acierta el 88,9 % de los días y la suavizada el 95,3 %. La matriz de confusión de la filtrada muestra dónde están los errores: sobre todo entre vecinos (alcista con lateral, lateral con bajista) y casi nunca entre los extremos. La diferencia de siete puntos entre las dos es información del futuro. Se ve mejor alrededor de cada entrada en el régimen bajista:

entradas = np.flatnonzero((verdad[1:] == 2) & (verdad[:-1] != 2)) + 1
entradas = entradas[(entradas >= 10) & (entradas + 20 < len(r))]
dias = np.arange(-10, 21)
evento = pd.DataFrame(
    {"filtrada": [filtrada[entradas + d, 2].mean() for d in dias],
     "suavizada": [suavizada[entradas + d, 2].mean() for d in dias]},
    index=pd.Index(dias, name="dia"))

def retraso(p):
    """Sesiones hasta que P(bajista) > 0,5 en cada entrada."""
    salida = []
    for e in entradas:
        fin = e
        while fin < len(verdad) and verdad[fin] == 2:
            fin += 1
        dentro = np.flatnonzero(p[e:fin, 2] > 0.5)
        salida.append(dentro[0] if len(dentro) else np.nan)
    return pd.Series(salida)

print(len(entradas), "entradas; retraso mediano:",
      {"filtrada": retraso(filtrada).median(),
       "suavizada": retraso(suavizada).median()})
evento.loc[[-10, -5, -2, -1, 0, 1, 2, 3, 5, 10, 20]].round(2)

28 entradas; retraso mediano: {'filtrada': 3.0, 'suavizada': 0.0}
     filtrada  suavizada
dia
-10      0.11       0.11
-5       0.10       0.17
-2       0.08       0.27
-1       0.07       0.36
 0       0.20       0.51
 1       0.24       0.58
 2       0.39       0.67
 3       0.47       0.69
 5       0.55       0.71
 10      0.63       0.67
 20      0.46       0.47

La probabilidad suavizada del régimen bajista empieza a subir días antes de que empiece: vale 0,36 el día anterior, cuando la filtrada vale 0,07. Nada en los rendimientos hasta ese día anunciaba el cambio; la suavizada lo «sabe» porque ha visto los días agitados que vienen después. La filtrada reacciona al ver los primeros días malos y cruza el 0,5 con un retraso mediano de 3 sesiones; la suavizada, el mismo día.

Lo que la anticipación hace a un backtest

Una regla sencilla: estar invertido salvo cuando la probabilidad del régimen bajista supera 0,5, decidiendo con la probabilidad del día anterior. Con la filtrada y con la suavizada, frente a comprar y mantener y frente a un oráculo que conoce el estado verdadero de cada día:

def resultado(pos):
    pnl = r.loc[pos.index] * pos
    return pd.Series({"sharpe": pnl.mean() / pnl.std() * np.sqrt(252),
                      "rent_total": np.expm1(pnl.sum()),
                      "en_mercado": pos.mean()})

def fuera_si_bajista(p, index):
    return pd.Series(p[:, 2] < 0.5, index=index).astype(float)

siempre = pd.Series(1.0, index=r.index)
pd.DataFrame({
    "comprar_y_mantener": resultado(siempre),
    "filtrada": resultado(
        fuera_si_bajista(filtrada, r.index).shift(1, fill_value=1.0)),
    "suavizada": resultado(
        fuera_si_bajista(suavizada, r.index).shift(1, fill_value=1.0)),
    "oraculo": resultado(
        pd.Series(verdad != 2, index=r.index).astype(float)),
}).T.round(2)

                    sharpe  rent_total  en_mercado
comprar_y_mantener    0.33        5.01        1.00
filtrada              0.44        6.35        0.91
suavizada             0.58       10.05        0.90
oraculo               0.62       11.69        0.90

Con la filtrada, el filtro de régimen sube el Sharpe de 0,33 a 0,44: en estos datos hay un régimen bajista de verdad y salir de él ayuda. Con la suavizada, el mismo código da 0,58, casi el del oráculo (0,62). La diferencia entre 0,44 y 0,58 no es una ventaja: es el futuro filtrándose en el pasado. Y aun la versión filtrada tiene una anticipación más discreta: sus parámetros se estimaron con los 26 años, incluidos los que el filtro «recorre».

Walk-forward: parámetros sin futuro

La versión honesta estima el modelo con un tramo inicial y lo aplica, sin reestimarlo, al resto. Con 2000-2012 para estimar y 2013-2025 para probar:

ent = (r.index < "2013-01-01")
candidatos = [ajustar(X[ent], s) for s in range(5)]
m_ent = max(candidatos, key=lambda m: m.score(X[ent]))
o = por_volatilidad(m_ent)
f_oos = filtrar(r, m_ent.startprob_[o], m_ent.transmat_[np.ix_(o, o)],
                m_ent.means_.ravel()[o],
                np.sqrt(m_ent.covars_.ravel()[o]))
prueba = ~ent
print("vol. estimadas con 2000-2012:",
      np.sqrt(m_ent.covars_.ravel()[o]).round(4))
print(f"acierto 2013-2025, filtrada: "
      f"{(f_oos[prueba].argmax(1) == verdad[prueba]).mean():.3f}")
pos = fuera_si_bajista(f_oos, r.index).shift(1, fill_value=1.0)
pd.DataFrame({"comprar_y_mantener": resultado(siempre[prueba]),
              "hmm_fuera_de_muestra": resultado(pos[prueba])}
             ).T.round(2)

vol. estimadas con 2000-2012: [0.0076 0.0142 0.026 ]
acierto 2013-2025, filtrada: 0.916
                      sharpe  rent_total  en_mercado
comprar_y_mantener      0.67        3.93        1.00
hmm_fuera_de_muestra    0.82        4.72        0.95

Los parámetros estimados con los primeros trece años se parecen a los verdaderos, y fuera de la muestra la probabilidad filtrada acierta el 91,6 % de los días; el filtro de régimen sube el Sharpe de 2013-2025 de 0,67 a 0,82, con un 95 % del tiempo invertido. El acierto es más alto que en toda la muestra, y no por mérito: el generador usa los mismos parámetros durante los 26 años, así que lo aprendido en 2000-2012 sigue valiendo después. En un mercado real los regímenes no son estables (la volatilidad del «régimen tranquilo» de 2005 no es la de 2017) y el modelo se degrada fuera de la muestra. Esta es la comparación que hay que hacer siempre, y la que dice cuánto queda de una ventaja cuando se quita el futuro.

Una ruptura estructural: media y varianza

La serie ruptura tiene rendimientos normales independientes con un solo cambio el 2013-05-22: la media diaria pasa de 0,0004 a 0 y la volatilidad, de 0,008 a 0,013. Dos cambios en la misma fecha, uno pequeño y otro grande en relación con el ruido.

from scipy import stats

rp = (reg[reg["serie"] == "ruptura"]
      .dropna(subset=["rendimiento_log"])
      .set_index("fecha")["rendimiento_log"])
corte = rp.index.get_loc(pd.Timestamp("2013-05-22"))
antes, despues = rp.iloc[:corte], rp.iloc[corte:]
pd.DataFrame({"antes": [len(antes), antes.mean(), antes.std()],
              "despues": [len(despues), despues.mean(),
                          despues.std()]},
             index=["sesiones", "media", "desviacion"]).round(5)

                 antes     despues
sesiones    3365.00000  3173.00000
media          0.00038     0.00032
desviacion     0.00780     0.01295

La volatilidad salió como se plantó (0,0078 y 0,0129). La media no: el tramo posterior, que debería tener media 0, salió en 0,00032, casi lo mismo que el anterior (0,00038). Con una volatilidad del 1 % diario, el azar de 3.000 sesiones mueve la media muestral en ±0,00015 o más, del mismo orden que el cambio plantado.

El contraste clásico para una fecha conocida es el de Chow (1960): ajustar el modelo (aquí, solo una media) en toda la muestra y en cada tramo, y comparar los residuos con un estadístico F. Aplicado a los rendimientos, contrasta la media; aplicado a sus cuadrados, la varianza.

def chow_media(x, k):
    """F de Chow para un cambio de media en la posición k."""
    x = np.asarray(x, float)
    sse = ((x - x.mean()) ** 2).sum()
    a, b = x[:k], x[k:]
    sse2 = ((a - a.mean()) ** 2).sum() + ((b - b.mean()) ** 2).sum()
    F = (sse - sse2) / (sse2 / (len(x) - 2))
    return F, stats.f.sf(F, 1, len(x) - 2)

print("media:    F = {:.3f}, p = {:.3f}".format(*chow_media(rp, corte)))
print("varianza: F = {:.1f}, p = {:.1e}".format(
    *chow_media(rp ** 2, corte)))
# sesiones por tramo para detectar 0,0004 con un 80 % de potencia
n = 2 * (stats.norm.ppf(0.975) + stats.norm.ppf(0.8)) ** 2
print(f"sesiones necesarias por tramo: "
      f"{n * rp.var() / 0.0004 ** 2:,.0f}")

media:    F = 0.057, p = 0.811
varianza: F = 598.2, p = 1.7e-126
sesiones necesarias por tramo: 11,048

El cambio de varianza es inequívoco (F = 598) y el de media, invisible (p = 0,81). No es un fallo del contraste: para detectar un cambio de 0,04 % diario con esa volatilidad harían falta unas 11.000 sesiones a cada lado, más de cuarenta años. Es la misma aritmética que hace tan difícil demostrar que una estrategia tiene ventaja (lección 2.1.8): las medias de los rendimientos se estiman muy mal.

Buscar la fecha sin saberla

En la práctica la fecha no se conoce. Lo natural es calcular el F de Chow en cada fecha candidata (dejando fuera el 15 % de cada extremo, para que los dos tramos tengan datos) y quedarse con el máximo. Ese máximo no se compara con una F(1, n): es el mejor de miles de contrastes, y su distribución es otra (Andrews, 1993). Su valor crítico se puede obtener simulando series sin ruptura del mismo tamaño:

def barrido_chow(x, recorte=0.15):
    """F de Chow de cambio de media en cada posición candidata."""
    x = np.asarray(x, float)
    n = len(x)
    k = np.arange(int(recorte * n), int((1 - recorte) * n))
    c1, c2 = np.cumsum(x), np.cumsum(x ** 2)
    sse = c2[-1] - c1[-1] ** 2 / n
    sse_a = c2[k - 1] - c1[k - 1] ** 2 / k
    sse_b = (c2[-1] - c2[k - 1]) - (c1[-1] - c1[k - 1]) ** 2 / (n - k)
    sse2 = sse_a + sse_b
    return k, (sse - sse2) / (sse2 / (n - 2))

rng = np.random.default_rng(2028)
nulo = {"media": [], "varianza": []}
for _ in range(300):
    z = rng.standard_normal(len(rp))
    nulo["media"].append(barrido_chow(z)[1].max())
    nulo["varianza"].append(barrido_chow(z ** 2)[1].max())

for nombre, serie in [("media", rp), ("varianza", rp ** 2)]:
    k, F = barrido_chow(serie)
    print(f"{nombre:8}  máximo F {F.max():6.1f} el "
          f"{rp.index[k[F.argmax()]].date()}  crítico 5 % "
          f"{np.percentile(nulo[nombre], 95):.1f}")
print(f"crítico de una sola fecha, F(1, n): "
      f"{stats.f.ppf(0.95, 1, len(rp) - 2):.2f}")

media     máximo F   10.1 el 2021-05-13  crítico 5 % 8.8
varianza  máximo F  598.5 el 2013-05-20  crítico 5 % 7.9
crítico de una sola fecha, F(1, n): 3.84

En la varianza, el barrido encuentra el 2013-05-20, a dos sesiones de la verdad, con un F que multiplica por 75 el valor crítico. En la media, encuentra una «ruptura» el 2021-05-13 con F = 10,1, por encima del crítico simulado (8,8) y muy por encima del 3,84 de una sola fecha. Es falsa: la única ruptura de media plantada es la de 2013, y es invisible. Dos cosas la explican. El F supone una varianza constante, y la varianza cambió en 2013 (tipificando cada tramo por su volatilidad, el máximo baja a 8,9); el resto es azar: un contraste al 5 % encuentra lo que no hay una de cada veinte veces, y un barrido es un contraste. Con 3,84 como referencia, se habría «encontrado» con mucha más facilidad.

Un contraste específico para la varianza es el de Inclán y Tiao (1994), basado en la suma acumulada de los cuadrados $C_k = \sum_{t \le k} r_t^2$. Si la varianza es constante, $C_k / C_T$ crece en línea recta con $k/T$; el estadístico mide la mayor desviación,

que sin ruptura se distribuye asintóticamente como el máximo de un puente browniano (la distribución de Kolmogorov, con valor crítico al 5 % de 1,358). Como controles sin ruptura sirven los dos tramos por separado:

from scipy.stats import kstwobign

def inclan_tiao(x):
    """Estadístico de Inclán y Tiao y posición de su máximo."""
    x = np.asarray(x, float) - np.mean(x)
    c = np.cumsum(x ** 2)
    n = len(x)
    d = c / c[-1] - np.arange(1, n + 1) / n
    k = int(np.abs(d).argmax())
    return np.sqrt(n / 2) * abs(d[k]), k

print(f"crítico al 5 %: {kstwobign.ppf(0.95):.3f}")
for nombre, x in [("serie completa", rp), ("solo antes", antes),
                  ("solo después", despues)]:
    est, k = inclan_tiao(x)
    print(f"{nombre:15} IT {est:6.2f}  máximo el {x.index[k].date()}")

crítico al 5 %: 1.358
serie completa  IT  13.55  máximo el 2013-05-17
solo antes      IT   0.74  máximo el 2002-04-25
solo después    IT   0.56  máximo el 2024-02-14

En la serie completa el estadístico vale 13,5, diez veces el crítico, con el máximo el 2013-05-17, a tres sesiones de la verdad. En los dos tramos sin ruptura queda por debajo del crítico (0,74 y 0,56): los controles no dan falsos positivos. El contraste supone rendimientos independientes y normales; con colas gruesas y agrupamiento de volatilidad (lección 1.6.2) rechaza de más, y en series reales se usan versiones corregidas.

Errores típicos con regímenes y rupturas

• Suavizar y operar. predict_proba y predict de hmmlearn usan toda la muestra. Para operar, probabilidades filtradas con parámetros estimados antes.
• Un solo arranque. EM converge a máximos locales: en estos datos, tres de cada ocho arranques dan un modelo malo.
• Leer los números de estado. Las etiquetas son arbitrarias; ordénalas por volatilidad o por media.
• Elegir el número de estados a ojo. Se compara con criterios de información (ejercicio 1) y, sobre todo, fuera de la muestra.
• Barrer fechas con valores críticos de una fecha. El máximo de muchos F necesita su propio valor crítico.
• Esperar ver rupturas de media. Una media pequeña frente a la volatilidad no se detecta ni con décadas de datos; una ruptura de varianza, en días.
• Reajustar el modelo cada día sin cuidado. Un HMM reestimado puede cambiar el orden de sus estados de un día a otro; si el sistema no los ordena, cambia de «régimen» sin que haya cambiado nada.

Ejercicios

1. Ajusta HMM de 2, 3 y 4 estados a la serie markov (cuatro arranques cada uno) y compara su BIC (método bic de hmmlearn). ¿Cuántos estados elige? Ten en cuenta que algún arranque con 4 estados puede fallar.

# con 4 estados, algún arranque deja un estado sin días: hmmlearn
# divide por cero (RuntimeWarning) y después falla (ValueError)
warnings.filterwarnings("ignore", category=RuntimeWarning,
                        module="hmmlearn")
bic = {}
for k in [2, 3, 4]:
    mejor_k = None
    for s in range(4):
        try:
            m = ajustar(X, s, k=k)
        except ValueError:      # un estado vacío: covarianza no válida
            continue
        if mejor_k is None or m.score(X) > mejor_k.score(X):
            mejor_k = m
    bic[k] = {"log_verosimilitud": mejor_k.score(X),
              "bic": mejor_k.bic(X)}
pd.DataFrame(bic).T.round(1)

   log_verosimilitud      bic
2            20154.6 -40247.7
3            20281.6 -40440.3
4            20284.0 -40365.9

2. Ajusta un HMM de dos estados a la serie ruptura. ¿Encuentra la ruptura? ¿Qué dice su matriz de transición?

Xr = rp.to_numpy().reshape(-1, 1)
m2 = max((ajustar(Xr, s, k=2) for s in range(3)),
         key=lambda m: m.score(Xr))
o2 = por_volatilidad(m2)
mapa = np.empty(2, int)
mapa[o2] = np.arange(2)
estado = mapa[m2.predict(Xr)]
real_r = (rp.index >= "2013-05-22").astype(int)
print(f"acierto {(estado == real_r).mean():.4f}; cambios de estado: "
      f"{(np.diff(estado) != 0).sum()}; primer día agitado: "
      f"{rp.index[np.argmax(estado == 1)].date()}")
pd.DataFrame(m2.transmat_[np.ix_(o2, o2)]).round(4)

acierto 0.9997; cambios de estado: 1; primer día agitado: 2013-05-20
        0       1
0  0.9997  0.0003
1  0.0000  1.0000

3. Un artículo de un blog presenta una estrategia que reduce exposición cuando un HMM detecta el régimen bajista, con un Sharpe de 1,1 frente a 0,6 del índice. El código usa model.fit(rendimientos) sobre 20 años y después model.predict(rendimientos). ¿Qué le preguntarías antes de creértelo?