---
title: "2.1.7 Detección de regímenes: HMM y cambios estructurales"
description: "HMM con hmmlearn y contrastes de ruptura sobre series sintéticas con regímenes conocidos: precisión, anticipación del suavizado y fecha del cambio."
url: "https://pythonparatrading.com/curso-python-trading/nivel-2/regimenes-hmm-y-cambios-estructurales/"
lang: "es"
image: "https://pythonparatrading.com/og/curso-python-trading/nivel-2/regimenes-hmm-y-cambios-estructurales.png"
---

# 2.1.7 · Detección de regímenes: modelos de Markov ocultos y cambios estructurales

> Los mercados alternan épocas tranquilas y agitadas, y a veces cambian de una vez para siempre. Esta lección ajusta un modelo de Markov oculto a una serie sintética cuyos regímenes se conocen día a día, mide cuánto acierta y cuánto acierta de más cuando mira al futuro, y localiza una ruptura estructural con contrastes que no saben dónde está.

- Curso: [Nivel 2: Estrategias sistemáticas y backtesting riguroso](https://pythonparatrading.com/curso-python-trading/nivel-2/)
- Módulo: 2.1 · Estadística y econometría para trading
- Duración: 75 min
- Requisitos: [0.7.2 · Probabilidad y distribuciones, correlación frente a causalidad y regresión simple](https://pythonparatrading.com/curso-python-trading/nivel-0/probabilidad-distribuciones-y-regresion/), [1.6.2 · Hechos estilizados: colas gruesas, agrupamiento, autocorrelación y apalancamiento](https://pythonparatrading.com/curso-python-trading/nivel-1/hechos-estilizados-de-los-mercados/)

## Al terminar esta lección

- Ajustar un modelo de Markov oculto gaussiano con `hmmlearn`, elegir entre varios arranques por verosimilitud y resolver el intercambio de etiquetas.
- Medir la precisión del modelo frente a los estados verdaderos con una matriz de confusión.
- Calcular probabilidades filtradas desde cero y distinguirlas de las suavizadas, que usan información futura.
- Detectar una ruptura estructural en la media o en la varianza (Chow, barrido con valor crítico simulado, Inclán-Tiao) y saber qué rupturas son difíciles de ver.

## 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

$$
A = \begin{pmatrix} 0{,}992 & 0{,}006 & 0{,}002 \\ 0{,}012 & 0{,}980 & 0{,}008 \\ 0{,}010 & 0{,}020 & 0{,}970 \end{pmatrix},
$$

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.

```mermaid
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
```

```python
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)
```
```output
         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:

```python
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)
```
```output
         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
```

*Fig. — Ocho arranques del mismo HMM de tres estados: log-verosimilitud frente al acierto sobre los estados verdaderos.* Dispersión de ocho puntos. Tres arranques (semillas 0, 1 y 6) quedan en log-verosimilitudes de unas 20.146 con aciertos del 37 % al 59 %; la semilla 3, en 20.152 con un 86 %; los cuatro restantes, destacados, en 20.282 a 20.283 con un 94 %. La mayor verosimilitud coincide con el mayor acierto.

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

```python
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)
```
```output
           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:

```python
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)
```
```output
         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.

```mermaid
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.

```python
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)
```
```output
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:

```python
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)
```
```output
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
```

*Fig. — Probabilidad del régimen bajista alrededor de sus 28 entradas: la suavizada sube antes de que empiece; la filtrada, después.* Dos líneas de la sesión −10 a la +20 respecto al primer día del régimen bajista, con la media de 28 entradas. La filtrada se mantiene en torno a 0,1 hasta el día 0 (0,07 el día anterior), sube a 0,47 el día 3 y a 0,63 el día 10. La suavizada, destacada, empieza a subir una semana antes: 0,17 cinco días antes, 0,36 el día anterior y 0,51 el día 0.

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.

<aside class="callout" data-kind="riesgo">

Usar probabilidades suavizadas (o el estado de Viterbi, o un modelo ajustado con toda la muestra) para decidir la
posición de cada día es sesgo de anticipación (lección 2.5.1), y de los difíciles de ver: el código parece correcto,
la probabilidad «es del día t» y el resultado mejora. Un filtro de régimen solo vale con probabilidades filtradas y
con parámetros estimados antes de cada fecha.

</aside>

## 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:

```python
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)
```
```output
                    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
```

*Fig. — Sharpe de un filtro de régimen (fuera del mercado si P(bajista) > 0,5 el día anterior) según la probabilidad que use, en la serie markov.* Cuatro barras horizontales: Comprar y mantener, 0,33; Filtro con la filtrada, 0,44; Filtro con la suavizada: anticipación, 0,58; Oráculo: conoce el estado de hoy, 0,62. La suavizada, destacada, queda casi a la altura del oráculo sin ninguna ventaja real: la diferencia con la filtrada es anticipación.

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:

```python
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)
```
```output
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.

```python
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)
```
```output
                 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.

```python
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}")
```
```output
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:

```python
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}")
```
```output
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
```

*Fig. — Serie con ruptura, 2000–2025: volatilidad de 63 sesiones (arriba) y estadístico de Inclán-Tiao de cada fecha (abajo).* Arriba, la volatilidad anualizada de 63 sesiones oscila alrededor del 12,4 % hasta mayo de 2013 y alrededor del 20,6 % después, con una línea discontinua en la ruptura real del 2013-05-22. Abajo, el estadístico de Inclán-Tiao sube desde 2000 hasta un máximo de 13,5 el 2013-05-17, destacado, y baja hasta 0 al final; la línea del valor crítico del 5 %, 1,358, queda muy por debajo.

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,

$$
IT = \sqrt{T/2}\, \max_k \left| \frac{C_k}{C_T} - \frac{k}{T} \right| ,
$$

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:

```python
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()}")
```
```output
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.

**Solución.**

```python
# 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)
```
```output
   log_verosimilitud      bic
2            20154.6 -40247.7
3            20281.6 -40440.3
4            20284.0 -40365.9
```

El BIC más bajo es el de 3 estados, el número verdadero. Con 4, la verosimilitud apenas mejora y el BIC lo penaliza
por los parámetros de más; con 2, el modelo no puede separar el lateral del bajista y la verosimilitud cae mucho. En
datos reales el BIC suele pedir más estados de los que son útiles (las colas gruesas se «explican» con estados
extra); la decisión final se toma fuera de la muestra.

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

**Solución.**

```python
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)
```
```output
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
```

El HMM encuentra la ruptura casi exacta y asigna un solo cambio de estado. Su matriz de transición lo delata: la
probabilidad de volver del estado agitado al tranquilo es 0, y la de pasar del tranquilo al agitado, minúscula. Un HMM
ajustado a una ruptura única describe un estado «absorbente», no un régimen que va y viene; si un HMM ajustado a datos
reales sale así, lo que hay es un cambio estructural y no regímenes recurrentes.

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?

**Solución.**

Tres cosas. Primero, si `predict` (Viterbi) o `predict_proba` se aplican a la misma muestra con la que se ajustó:
las dos usan rendimientos posteriores a cada día, como la probabilidad suavizada de esta lección, que subía el Sharpe
de 0,44 a 0,58 sin ventaja real. Segundo, si los parámetros se estimaron con los 20 años: aunque se filtrara bien,
el modelo conoce la volatilidad de las crisis que luego «detecta». Tercero, si la posición de cada día usa la
probabilidad del día anterior, cuántos arranques, estados y umbrales se probaron (lección 2.1.8) y qué pasa con
costes. La prueba que lo resolvería es un walk-forward con probabilidades filtradas y reestimación solo con datos
pasados.

## Autoevaluación

1. Ajustas el mismo HMM de tres estados con ocho semillas distintas y obtienes log-verosimilitudes entre 20.146 y 20.283 y aciertos entre el 37 % y el 94 %. ¿Qué haces?

   - A. Promediar los ocho modelos
   - B. Quedarte con el de mayor log-verosimilitud: el algoritmo EM converge a máximos locales y la semilla decide a cuál
   - C. Quedarte con el de mayor acierto
   - D. Usar la semilla 0, que es la convención

   *Respuesta B: Quedarte con el de mayor log-verosimilitud: el algoritmo EM converge a máximos locales y la semilla decide a cuál.* En datos reales no hay estados verdaderos con los que medir el acierto; la verosimilitud sí se puede comparar. En la serie del curso, el máximo de log-verosimilitud coincide con los modelos que aciertan el 94 %.

2. El estado 0 del HMM resulta ser el de mayor volatilidad, aunque en la verdad el 0 es el alcista. ¿Es un error?

   - A. Sí: hay que cambiar de semilla hasta que coincidan
   - B. No: las etiquetas de un HMM son arbitrarias (intercambio de etiquetas). Hay que ordenarlas por algo con significado, como la volatilidad
   - C. Sí: el modelo no ha convergido
   - D. No, pero significa que el mercado estaba al revés

   *Respuesta B: No: las etiquetas de un HMM son arbitrarias (intercambio de etiquetas). Hay que ordenarlas por algo con significado, como la volatilidad.* La verosimilitud no cambia si se permutan los estados, así que el número de cada estado no significa nada. Ordenar por volatilidad (o por media) hace comparables modelos y ajustes.

3. ¿Por qué no se pueden usar las probabilidades de `predict_proba` de hmmlearn para decidir la posición de cada día en un backtest?

   - A. Porque no suman 1
   - B. Porque son probabilidades suavizadas: la de cada día se calcula con todos los rendimientos de la muestra, también los posteriores
   - C. Porque hmmlearn es lento
   - D. Porque están en logaritmos

   *Respuesta B: Porque son probabilidades suavizadas: la de cada día se calcula con todos los rendimientos de la muestra, también los posteriores.* El algoritmo hacia delante y hacia atrás usa toda la muestra. En la serie del curso, la probabilidad suavizada del régimen bajista ya vale 0,36 el día antes de que empiece, cuando la filtrada vale 0,07; un filtro de salida con la suavizada sube el Sharpe de 0,44 a 0,58 sin ninguna ventaja real.

4. La media diaria de la serie con ruptura pasa de 0,0004 a 0 el 2013-05-22, y el contraste de Chow en esa fecha da p = 0,81. ¿Qué concluyes?

   - A. Que no hubo ruptura
   - B. Que el contraste tiene muy poca potencia para cambios de media tan pequeños frente a la volatilidad: harían falta décadas de datos a cada lado
   - C. Que el contraste está mal programado
   - D. Que la ruptura fue en la varianza y no puede haberla en la media

   *Respuesta B: Que el contraste tiene muy poca potencia para cambios de media tan pequeños frente a la volatilidad: harían falta décadas de datos a cada lado.* Con una volatilidad diaria del 1 %, distinguir una media de 0,04 % de una de 0 con un 80 % de potencia exige unas 11.000 sesiones a cada lado. Además, en la muestra la media posterior salió en 0,032 %: el azar se comió casi todo el cambio.

5. Barres todas las fechas posibles con un Chow y el máximo de F vale 10,1. El valor crítico de una F(1, n) al 5 % es 3,84. ¿Hay ruptura?

   - A. Sí, con mucha holgura
   - B. No se puede saber con ese valor crítico: el máximo de muchos F no se distribuye como un F. Con el crítico simulado del máximo (8,8 en esta muestra) el rechazo es justo, y en la serie del curso es un falso positivo
   - C. No, porque 10,1 es poco
   - D. Sí, porque cualquier F mayor que 1 indica ruptura

   *Respuesta B: No se puede saber con ese valor crítico: el máximo de muchos F no se distribuye como un F. Con el crítico simulado del máximo (8,8 en esta muestra) el rechazo es justo, y en la serie del curso es un falso positivo.* Elegir la fecha que maximiza F es hacer miles de contrastes y quedarse con el mejor: es un problema de pruebas múltiples (lección 2.1.8). El valor crítico del máximo hay que obtenerlo por simulación o de las tablas de Andrews (1993).

## Ideas clave

- Un HMM supone que un estado oculto con dinámica de Markov decide la media y la volatilidad de cada día; EM lo estima con máximos locales: varios arranques y el de mayor verosimilitud.
- Las etiquetas son arbitrarias: ordénalas por volatilidad antes de interpretar o comparar.
- En la serie del curso, el mejor modelo acierta el estado del 94 % de los días (Viterbi); la probabilidad filtrada, la única utilizable en tiempo real, el 89 %.
- Las probabilidades suavizadas usan el futuro: en un backtest inflan el Sharpe de un filtro de régimen de 0,44 a 0,58.
- Una ruptura de varianza se localiza con precisión de días; una de media pequeña es casi invisible. Barrer fechas exige valores críticos del máximo, no del F de una fecha.

## Recursos

### Research

- [Volatility Targeting en Python](https://pythonparatrading.com/volatility-targeting/) · Research-0004 · 8 min

### Documentación

- [hmmlearn: tutorial (GaussianHMM, varios arranques, predict_proba)](https://hmmlearn.readthedocs.io/en/stable/tutorial.html) · hmmlearn.readthedocs.io

### Libros

- [Hamilton, J. D. (1989): A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle. Econometrica 57(2), 357-384](https://doi.org/10.2307/1912559) · doi.org
- [Inclán, C. y Tiao, G. C. (1994): Use of Cumulative Sums of Squares for Retrospective Detection of Changes of Variance. JASA 89(427), 913-923](https://doi.org/10.1080/01621459.1994.10476824) · doi.org
- [Andrews, D. W. K. (1993): Tests for Parameter Instability and Structural Change With Unknown Change Point. Econometrica 61(4), 821-856](https://doi.org/10.2307/2951764) · doi.org
- [Chow, G. C. (1960): Tests of Equality Between Sets of Coefficients in Two Linear Regressions. Econometrica 28(3), 591-605](https://doi.org/10.2307/1910133) · doi.org

### Datos

- [Datos del nivel 2: series con regímenes de Markov y con una ruptura, 2000-2025 (Parquet)](https://pythonparatrading.com/curso/datos/nivel-2/regimenes2.parquet) · Descarga · regimenes2.parquet

Lección anterior: [2.1.6 · Filtro de Kalman para hedge ratios dinámicos](https://pythonparatrading.com/curso-python-trading/nivel-2/filtro-de-kalman-hedge-ratio-dinamico/)

Lección siguiente: [2.1.8 · El problema de las pruebas múltiples: Bonferroni, FDR y el «factor zoo»](https://pythonparatrading.com/curso-python-trading/nivel-2/pruebas-multiples-bonferroni-fdr-y-factor-zoo/)