Regime detection: hidden Markov models and structural breaks.

Regimes: the market does not have a single behavior

Lessons 1.6.1 and 1.6.2 measured that volatility comes in streaks. A step further is to assume that the market alternates between a few regimes (calm uptrend, nervous sideways market, panic selloff), each with its own mean and volatility, and that it moves from one to another with some persistence. If you knew which regime you were in, you could adjust exposure, position sizing or the strategy (lesson 2.2.2). The problem is that the regime is not visible: only the returns are.

There are two ways to think about a change in behavior, and this lesson covers both:

• Recurring regimes. The market enters and leaves the same states many times. The natural model is a hidden Markov model (HMM), popularized in econometrics by Hamilton (1989).
• Structural break. Something changes once and does not come back: a regulation, a change of monetary regime, a company that changes its business. The natural tool is a structural change test that searches for the date.

The data are the two synthetic series in regimenes2.parquet (NYSE sessions from 2000 to 2025). In markov, the returns are normal with a mean and volatility decided by a three-state Markov chain; in ruptura, there is a single change on 2013-05-22. In both, the true state of each day is in the estado_real column: the known answer against which each method is measured. They are made-up data: the parameters come from the generator, not from any market.

The hidden Markov chain

An HMM has two pieces. The hidden chain $s_t \in \{0, 1, 2\}$ moves from one state to another with a transition matrix $A$, with $A_{ij} = P(s_{t+1} = j \mid s_t = i)$. The emissions are the returns: $r_t \mid s_t = k \sim N(\mu_k, \sigma_k^2)$. The true matrix of the markov series is

with a diagonal close to 1: the regimes are persistent. The average length of a stay in state $k$ is $1 / (1 - A_{kk})$: 125 sessions in the bull state, 50 in the sideways state and 33 in the bear state.

stateDiagram-v2
  direction LR
  Bull --> Sideways: 0.006
  Sideways --> Bull: 0.012
  Sideways --> Bear: 0.008
  Bear --> Sideways: 0.020
  Bull --> Bear: 0.002
  Bear --> Bull: 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 = ["bull", "sideways", "bear"]
resumen = mk.groupby("estado_real")["rendimiento_log"].agg(
    ["size", "mean", "std"])
resumen.index = nombres
resumen.round(5)

          size     mean      std
bull      3479  0.00086  0.00700
sideways  2407 -0.00018  0.01417
bear       652 -0.00115  0.02670

The figure at the top shows the series with its regimes: the price drops coincide with the bear stretches, but not all turbulent stretches are drops, and between one regime and another there is no signal visible to the naked eye. The volatilities measured on the days of each state match those of the generator (0.007, 0.014 and 0.028 daily); the means, only roughly (0.00086, −0.00018 and −0.00115, versus 0.0006, 0 and −0.0018): with a few hundred days per state, the sample mean strays far from the true one, and the volatility hardly at all.

Fitting the HMM: several starts

hmmlearn estimates a Gaussian HMM with the EM (expectation-maximization) algorithm: it alternates between computing the probability of each state on each day with the current parameters and re-estimating the parameters with those probabilities. EM raises the likelihood at every step, but it can get stuck in a local maximum, and which one it reaches depends on the starting point, which hmmlearn draws with random_state. The correct practice is to try several starts:

from hmmlearn.hmm import GaussianHMM

# hmmlearn warns through logging when an EM iteration lowers the
# likelihood in the sixth digit ("Model is not converging") and
# uses an operation that numpy 2.5 declares deprecated
# (DeprecationWarning). Neither changes the result.
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):
    """States sorted from lowest to highest volatility."""
    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

Three of the eight starts (seeds 0, 1 and 6) end up in bad local maxima, with log-likelihoods of about 20,146 and accuracies of 37% to 59%: the model has split the states another way (for example, it has split the calm regime in two and mixed the sideways one with the bear one). Another (seed 3) stops halfway, with 20,152 and 86%. The other four reach 20,282-20,283 and are right 94% of the time. The practical rule does not need the truth: keep the start with the highest likelihood, which here is also one of the most accurate.

The best model and label switching

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
bull      0.00082     0.00086  0.00726   0.00700
sideways -0.00022    -0.00018  0.01438   0.01417
bear     -0.00088    -0.00115  0.02646   0.02670
           bull  sideways   bear
bull      0.990     0.008  0.002
sideways  0.012     0.977  0.011
bear      0.009     0.037  0.955

The _real columns are the means and volatilities of the days that were actually in each state. The volatilities are recovered almost exactly; the means, less well (the bear one comes out at −0.00088, versus −0.00115 on its days and −0.0018 in the generator), because a daily mean is estimated with much more error than a volatility: it is the same problem as the mean break further down. The transition matrix has the right structure, with the persistence of the bear state somewhat underestimated.

The por_volatilidad step is not cosmetic. The labels 0, 1, 2 that hmmlearn returns are arbitrary: if the states are permuted, the likelihood is the same, so nothing forces state 0 to be the calm one. This is label switching. With the four good starts, the raw label that most days of each true regime receive:

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

          seed 2  seed 4  seed 5  seed 7
verdad
bull           1       0       0       1
sideways       0       2       1       0
bear           2       1       2       2

The four models are practically the same, but they number the regimes differently: in one the bull regime is 0 and in another it is 1. Comparing the raw state with the truth would give absurd accuracies, and two fits of the same model (or yesterday's and today's in a system that re-estimates) can change the number of each regime. Sorting them by volatility, or by mean, is what gives them meaning.

Filtering or smoothing: what was known each day

With the estimated parameters, the question that matters remains: which regime is the market in today? There are two different answers:

• The filtered probability, $P(s_t = k \mid r_1, \dots, r_t)$, uses the returns up to today. It is what a system knows at the moment of deciding.
• The smoothed probability, $P(s_t = k \mid r_1, \dots, r_T)$, uses the whole sample, including what happened after $t$. It is what hmmlearn's predict_proba returns (forward-backward algorithm), and the Viterbi state from predict also uses the whole sample.

flowchart TB
  P[Returns 1..T] --> F[Forward]
  F --> PF["Filtered: up to t"]
  F --> B[Backward]
  B --> PS["Smoothed: up to T"]
  PF --> U[Usable for trading]
  PS --> N[Only for study]

The filtered one is coded in a few lines with the forward algorithm: yesterday's probability is projected with the transition matrix and corrected with the density of today's return in each state.

from scipy.stats import norm

def filtrar(x, pi0, A, medias, vols):
    """P(state_t | returns up to t): forward algorithm."""
    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]          # correction with today's data
        alfa[t] = a / a.sum()
        prior = alfa[t] @ A          # prediction for tomorrow
    return alfa

filtrada = filtrar(r, pi0, A, medias, vols)
suavizada = mejor.predict_proba(X)[:, orden]
for nombre, p in [("filtered", filtrada), ("smoothed", suavizada)]:
    print(f"{nombre:9}  accuracy {(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)

filtered   accuracy 0.889
smoothed   accuracy 0.953
filtrada  bull  sideways  bear
verdad
bull      3281       185    13
sideways   257      2045   105
bear         5       162   485

The filtered probability is right on 88.9% of days and the smoothed one on 95.3%. The confusion matrix of the filtered one shows where the errors are: mostly between neighbors (bull with sideways, sideways with bear) and almost never between the extremes. The seven-point difference between the two is information from the future. It shows better around each entry into the bear regime:

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):
    """Sessions until P(bear) > 0.5 at each entry."""
    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), "entries; median delay:",
      {"filtrada": retraso(filtrada).median(),
       "suavizada": retraso(suavizada).median()})
evento.loc[[-10, -5, -2, -1, 0, 1, 2, 3, 5, 10, 20]].round(2)

28 entries; median delay: {'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

The smoothed probability of the bear regime starts to rise days before it begins: it is 0.36 the day before, when the filtered one is 0.07. Nothing in the returns up to that day announced the change; the smoothed one “knows” because it has seen the turbulent days that come afterward. The filtered one reacts when it sees the first bad days and crosses 0.5 with a median delay of 3 sessions; the smoothed one, on the same day.

What look-ahead does to a backtest

A simple rule: stay invested except when the probability of the bear regime exceeds 0.5, deciding with the previous day's probability. With the filtered and with the smoothed probability, versus buy and hold and versus an oracle that knows the true state of each day:

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

With the filtered probability, the regime filter raises the Sharpe from 0.33 to 0.44: in this data there is a real bear regime and getting out of it helps. With the smoothed one, the same code gives 0.58, almost the oracle's (0.62). The difference between 0.44 and 0.58 is not an edge: it is the future leaking into the past. And even the filtered version has a more subtle look-ahead: its parameters were estimated with the 26 years, including the ones the filter “walks through”.

Walk-forward: parameters without the future

The honest version estimates the model on an initial stretch and applies it, without re-estimating it, to the rest. With 2000-2012 to estimate and 2013-2025 to test:

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("vols estimated with 2000-2012:",
      np.sqrt(m_ent.covars_.ravel()[o]).round(4))
print(f"accuracy 2013-2025, filtered: "
      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)

vols estimated with 2000-2012: [0.0076 0.0142 0.026 ]
accuracy 2013-2025, filtered: 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

The parameters estimated with the first thirteen years resemble the true ones, and out of sample the filtered probability is right on 91.6% of days; the regime filter raises the 2013-2025 Sharpe from 0.67 to 0.82, with 95% of the time invested. The accuracy is higher than over the whole sample, and not on merit: the generator uses the same parameters throughout the 26 years, so what was learned in 2000-2012 still holds afterward. In a real market regimes are not stable (the volatility of the “calm regime” of 2005 is not that of 2017) and the model degrades out of sample. This is the comparison that always has to be made, and the one that tells how much of an edge remains when the future is removed.

A structural break: mean and variance

The ruptura series has independent normal returns with a single change on 2013-05-22: the daily mean goes from 0.0004 to 0 and the volatility, from 0.008 to 0.013. Two changes on the same date, one small and one large relative to the noise.

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

The volatility came out as planted (0.0078 and 0.0129). The mean did not: the later stretch, which should have a mean of 0, came out at 0.00032, almost the same as the earlier one (0.00038). With a volatility of 1% a day, the chance of 3,000 sessions moves the sample mean by ±0.00015 or more, of the same order as the planted change.

The classic test for a known date is the Chow test (1960): fit the model (here, just a mean) on the whole sample and on each stretch, and compare the residuals with an F statistic. Applied to the returns, it tests the mean; applied to their squares, the variance.

def chow_media(x, k):
    """Chow F for a change in the mean at position 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("mean:     F = {:.3f}, p = {:.3f}".format(*chow_media(rp, corte)))
print("variance: F = {:.1f}, p = {:.1e}".format(
    *chow_media(rp ** 2, corte)))
# sessions per stretch to detect 0.0004 with 80% power
n = 2 * (stats.norm.ppf(0.975) + stats.norm.ppf(0.8)) ** 2
print(f"sessions needed per stretch: "
      f"{n * rp.var() / 0.0004 ** 2:,.0f}")

mean:     F = 0.057, p = 0.811
variance: F = 598.2, p = 1.7e-126
sessions needed per stretch: 11,048

The change in variance is unambiguous (F = 598) and the change in mean, invisible (p = 0.81). It is not a failure of the test: detecting a change of 0.04% a day with that volatility would take about 11,000 sessions on each side, more than forty years. It is the same arithmetic that makes it so hard to prove that a strategy has an edge (lesson 2.1.8): the means of returns are estimated very poorly.

Searching for the date without knowing it

In practice the date is not known. The natural thing is to compute the Chow F at each candidate date (leaving out 15% at each end, so that both stretches have data) and keep the maximum. That maximum is not compared with an F(1, n): it is the best of thousands of tests, and its distribution is different (Andrews, 1993). Its critical value can be obtained by simulating series of the same size without a break:

def barrido_chow(x, recorte=0.15):
    """Chow F for a change in the mean at each candidate position."""
    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 = {"mean": [], "variance": []}
for _ in range(300):
    z = rng.standard_normal(len(rp))
    nulo["mean"].append(barrido_chow(z)[1].max())
    nulo["variance"].append(barrido_chow(z ** 2)[1].max())

for nombre, serie in [("mean", rp), ("variance", rp ** 2)]:
    k, F = barrido_chow(serie)
    print(f"{nombre:8}  max F {F.max():6.1f} on "
          f"{rp.index[k[F.argmax()]].date()}  5% critical "
          f"{np.percentile(nulo[nombre], 95):.1f}")
print(f"single-date critical value, F(1, n): "
      f"{stats.f.ppf(0.95, 1, len(rp) - 2):.2f}")

mean      max F   10.1 on 2021-05-13  5% critical 8.8
variance  max F  598.5 on 2013-05-20  5% critical 7.9
single-date critical value, F(1, n): 3.84

In the variance, the sweep finds 2013-05-20, two sessions from the truth, with an F 75 times the critical value. In the mean, it finds a “break” on 2021-05-13 with F = 10.1, above the simulated critical value (8.8) and far above the 3.84 of a single date. It is false: the only planted mean break is the 2013 one, and it is invisible. Two things explain it. The F assumes a constant variance, and the variance changed in 2013 (standardizing each stretch by its volatility, the maximum falls to 8.9); the rest is chance: a 5% test finds what is not there one time in twenty, and a sweep is a test. With 3.84 as the reference, it would have been “found” much more easily.

A test specific to the variance is the Inclán and Tiao test (1994), based on the cumulative sum of squares $C_k = \sum_{t \le k} r_t^2$. If the variance is constant, $C_k / C_T$ grows in a straight line with $k/T$; the statistic measures the largest deviation,

which without a break is asymptotically distributed as the maximum of a Brownian bridge (the Kolmogorov distribution, with a 5% critical value of 1.358). The two stretches taken separately serve as break-free controls:

from scipy.stats import kstwobign

def inclan_tiao(x):
    """Inclán and Tiao statistic and the position of its maximum."""
    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"5% critical value: {kstwobign.ppf(0.95):.3f}")
for nombre, x in [("full series", rp), ("before only", antes),
                  ("after only", despues)]:
    est, k = inclan_tiao(x)
    print(f"{nombre:15} IT {est:6.2f}  max on {x.index[k].date()}")

5% critical value: 1.358
full series     IT  13.55  max on 2013-05-17
before only     IT   0.74  max on 2002-04-25
after only      IT   0.56  max on 2024-02-14

On the full series the statistic is 13.5, ten times the critical value, with the maximum on 2013-05-17, three sessions from the truth. On the two break-free stretches it stays below the critical value (0.74 and 0.56): the controls give no false positives. The test assumes independent, normal returns; with fat tails and volatility clustering (lesson 1.6.2) it over-rejects, and corrected versions are used on real series.

Common mistakes with regimes and breaks

• Smoothing and trading. hmmlearn's predict_proba and predict use the whole sample. For trading, use filtered probabilities with parameters estimated beforehand.
• A single start. EM converges to local maxima: in this data, three out of eight starts give a bad model.
• Reading the state numbers. Labels are arbitrary; sort them by volatility or by mean.
• Choosing the number of states by eye. Compare with information criteria (exercise 1) and, above all, out of sample.
• Sweeping dates with single-date critical values. The maximum of many Fs needs its own critical value.
• Expecting to see mean breaks. A mean that is small relative to the volatility is not detected even with decades of data; a variance break, in days.
• Refitting the model every day without care. A re-estimated HMM can change the order of its states from one day to the next; if the system does not sort them, it changes “regime” without anything having changed.

Exercises

1. Fit HMMs with 2, 3 and 4 states to the markov series (four starts each) and compare their BIC (hmmlearn's bic method). How many states does it choose? Keep in mind that some start with 4 states may fail.

# with 4 states, some start leaves a state with no days: hmmlearn
# divides by zero (RuntimeWarning) and then fails (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:      # an empty state: invalid covariance
            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. Fit a two-state HMM to the ruptura series. Does it find the break? What does its transition matrix say?

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"accuracy {(estado == real_r).mean():.4f}; state changes: "
      f"{(np.diff(estado) != 0).sum()}; first turbulent day: "
      f"{rp.index[np.argmax(estado == 1)].date()}")
pd.DataFrame(m2.transmat_[np.ix_(o2, o2)]).round(4)

accuracy 0.9997; state changes: 1; first turbulent day: 2013-05-20
        0       1
0  0.9997  0.0003
1  0.0000  1.0000

3. A blog post presents a strategy that reduces exposure when an HMM detects the bear regime, with a Sharpe of 1.1 versus 0.6 for the index. The code uses model.fit(rendimientos) over 20 years and then model.predict(rendimientos). What would you ask before believing it?