--- title: "2.1.7 Regime detection: HMMs and structural breaks" description: "HMMs with hmmlearn and break tests on synthetic series with known regimes: accuracy, look-ahead of smoothing and the date of the change." url: "https://pythonparatrading.com/en/python-trading-course/level-2/regime-detection-hmm-structural-breaks/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-2/regime-detection-hmm-structural-breaks.png" --- # 2.1.7 · Regime detection: hidden Markov models and structural breaks > Markets alternate between calm and turbulent periods, and sometimes they change once and for all. This lesson fits a hidden Markov model to a synthetic series whose regimes are known day by day, measures how often it is right and how much more often it is right when it looks into the future, and locates a structural break with tests that do not know where it is. - Course: [Level 2: Systematic strategies and rigorous backtesting](https://pythonparatrading.com/en/python-trading-course/level-2/) - Module: 2.1 · Statistics and econometrics for trading - Duration: 75 min - Prerequisites: [0.7.2 · Probability and distributions, correlation versus causation and simple regression](https://pythonparatrading.com/en/python-trading-course/level-0/probability-distributions-and-regression/), [1.6.2 · Stylized facts: fat tails, clustering, autocorrelation and leverage](https://pythonparatrading.com/en/python-trading-course/level-1/stylized-facts-of-financial-markets/) ## By the end of this lesson - Fit a Gaussian hidden Markov model with `hmmlearn`, choose among several starts by likelihood and resolve label switching. - Measure the accuracy of the model against the true states with a confusion matrix. - Compute filtered probabilities from scratch and tell them apart from smoothed ones, which use future information. - Detect a structural break in the mean or the variance (Chow, a sweep with a simulated critical value, Inclán-Tiao) and know which breaks are hard to see. ## 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 $$ A = \begin{pmatrix} 0.992 & 0.006 & 0.002 \\ 0.012 & 0.980 & 0.008 \\ 0.010 & 0.020 & 0.970 \end{pmatrix}, $$ 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. ```mermaid 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 ``` ```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 = ["bull", "sideways", "bear"] resumen = mk.groupby("estado_real")["rendimiento_log"].agg( ["size", "mean", "std"]) resumen.index = nombres resumen.round(5) ``` ```output 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: ```python 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) ``` ```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. — Eight starts of the same three-state HMM: log-likelihood against accuracy on the true states.* Scatter of eight points. Three starts (seeds 0, 1 and 6) sit at log-likelihoods of about 20,146 with accuracies of 37% to 59%; seed 3 at 20,152 with 86%; the other four, highlighted, at 20,282 to 20,283 with 94%. The highest likelihood matches the highest accuracy. 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 ```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 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: ```python 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) ``` ```output 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. ```mermaid 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. ```python 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) ``` ```output 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: ```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): """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) ``` ```output 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 ``` *Fig. — Probability of the bear regime around its 28 entries: the smoothed one rises before it starts; the filtered one, after.* Two lines from session −10 to +20 around the first day of the bear regime, averaged over 28 entries. The filtered probability stays near 0.1 until day 0 (0.07 the day before), rises to 0.47 on day 3 and 0.63 on day 10. The smoothed one, highlighted, starts rising a week earlier: 0.17 five days before, 0.36 the day before and 0.51 on day 0. 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: ```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 of a regime filter (out of the market if P(bear) > 0.5 the day before) depending on the probability it uses, on the markov series.* Four horizontal bars: Buy and hold, 0.33; Filter on the filtered, 0.44; Filter on the smoothed: look-ahead, 0.58; Oracle: knows today's state, 0.62. The smoothed one, highlighted, almost matches the oracle without any real edge: the gap to the filtered one is look-ahead. 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: ```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("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) ``` ```output 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. ```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 ``` 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. ```python 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}") ``` ```output 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: ```python 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}") ``` ```output 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 ``` *Fig. — Series with a break, 2000–2025: 63-session volatility (top) and the Inclán-Tiao statistic at each date (bottom).* Top, the annualised 63-session volatility moves around 12.4% until May 2013 and around 20.6% afterwards, with a dashed line at the true break on 2013-05-22. Bottom, the Inclán-Tiao statistic rises from 2000 to a maximum of 13.5 on 2013-05-17, highlighted, and falls to 0 at the end; the 5% critical value line, 1.358, sits far below. 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, $$ IT = \sqrt{T/2}\, \max_k \left| \frac{C_k}{C_T} - \frac{k}{T} \right| , $$ 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: ```python 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()}") ``` ```output 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. **Solution.** ```python # 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) ``` ```output log_verosimilitud bic 2 20154.6 -40247.7 3 20281.6 -40440.3 4 20284.0 -40365.9 ``` The lowest BIC is that of 3 states, the true number. With 4, the likelihood barely improves and the BIC penalizes it for the extra parameters; with 2, the model cannot separate the sideways regime from the bear one and the likelihood falls a lot. With real data the BIC tends to ask for more states than are useful (fat tails get “explained” with extra states); the final decision is made out of sample. 2. Fit a two-state HMM to the `ruptura` series. Does it find the break? What does its transition matrix say? **Solution.** ```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"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) ``` ```output accuracy 0.9997; state changes: 1; first turbulent day: 2013-05-20 0 1 0 0.9997 0.0003 1 0.0000 1.0000 ``` The HMM finds the break almost exactly and assigns a single change of state. Its transition matrix gives it away: the probability of going back from the turbulent state to the calm one is 0, and that of going from the calm one to the turbulent one, tiny. An HMM fitted to a single break describes an “absorbing” state, not a regime that comes and goes; if an HMM fitted to real data comes out like this, what there is is a structural change and not recurring regimes. 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? **Solution.** Three things. First, whether `predict` (Viterbi) or `predict_proba` are applied to the same sample the model was fitted on: both use returns after each day, like the smoothed probability in this lesson, which raised the Sharpe from 0.44 to 0.58 without any real edge. Second, whether the parameters were estimated with the 20 years: even if filtering were done properly, the model knows the volatility of the crises it later “detects”. Third, whether each day's position uses the previous day's probability, how many starts, states and thresholds were tried (lesson 2.1.8) and what happens with costs. The test that would settle it is a walk-forward with filtered probabilities and re-estimation only with past data. ## Self-assessment 1. You fit the same three-state HMM with eight different seeds and get log-likelihoods between 20,146 and 20,283 and accuracies between 37% and 94%. What do you do? - A. Average the eight models - B. Keep the one with the highest log-likelihood: the EM algorithm converges to local maxima and the seed decides which one - C. Keep the one with the highest accuracy - D. Use seed 0, which is the convention *Answer B: Keep the one with the highest log-likelihood: the EM algorithm converges to local maxima and the seed decides which one.* With real data there are no true states to measure accuracy against; likelihood, on the other hand, can be compared. In the course series, the maximum log-likelihood coincides with the models that are right 94% of the time. 2. State 0 of the HMM turns out to be the one with the highest volatility, even though in the truth state 0 is the bull one. Is it an error? - A. No: HMM labels are arbitrary (label switching). They have to be sorted by something meaningful, such as volatility - B. Yes: you have to change the seed until they match - C. Yes: the model has not converged - D. No, but it means the market was upside down *Answer A: No: HMM labels are arbitrary (label switching). They have to be sorted by something meaningful, such as volatility.* The likelihood does not change if the states are permuted, so the number of each state means nothing. Sorting by volatility (or by mean) makes models and fits comparable. 3. Why can't the probabilities from hmmlearn's `predict_proba` be used to decide each day's position in a backtest? - A. Because they do not sum to 1 - B. Because hmmlearn is slow - C. Because they are in logarithms - D. Because they are smoothed probabilities: each day's is computed with all the returns in the sample, including the later ones *Answer D: Because they are smoothed probabilities: each day's is computed with all the returns in the sample, including the later ones.* The forward-backward algorithm uses the whole sample. In the course series, the smoothed probability of the bear regime is already 0.36 the day before it begins, when the filtered one is 0.07; an exit filter using the smoothed one raises the Sharpe from 0.44 to 0.58 without any real edge. 4. The daily mean of the series with a break goes from 0.0004 to 0 on 2013-05-22, and the Chow test on that date gives p = 0.81. What do you conclude? - A. That there was no break - B. That the test is badly coded - C. That the test has very little power for changes in the mean this small relative to volatility: decades of data would be needed on each side - D. That the break was in the variance and there cannot be one in the mean *Answer C: That the test has very little power for changes in the mean this small relative to volatility: decades of data would be needed on each side.* With a daily volatility of 1%, telling a mean of 0.04% from one of 0 with 80% power requires about 11,000 sessions on each side. Moreover, in the sample the later mean came out at 0.032%: chance ate almost all of the change. 5. You sweep all possible dates with a Chow test and the maximum F is 10.1. The 5% critical value of an F(1, n) is 3.84. Is there a break? - A. It cannot be known with that critical value: the maximum of many Fs is not distributed as an F. With the simulated critical value of the maximum (8.8 in this sample) the rejection is marginal, and in the course series it is a false positive - B. Yes, by a wide margin - C. No, because 10.1 is small - D. Yes, because any F greater than 1 indicates a break *Answer A: It cannot be known with that critical value: the maximum of many Fs is not distributed as an F. With the simulated critical value of the maximum (8.8 in this sample) the rejection is marginal, and in the course series it is a false positive.* Choosing the date that maximizes F is running thousands of tests and keeping the best one: it is a multiple testing problem (lesson 2.1.8). The critical value of the maximum has to be obtained by simulation or from the Andrews (1993) tables. ## Key points - An HMM assumes that a hidden state with Markov dynamics decides each day's mean and volatility; EM estimates it with local maxima: several starts and keep the one with the highest likelihood. - Labels are arbitrary: sort them by volatility before interpreting or comparing. - In the course series, the best model gets the state right on 94% of days (Viterbi); the filtered probability, the only one usable in real time, on 89%. - Smoothed probabilities use the future: in a backtest they inflate the Sharpe of a regime filter from 0.44 to 0.58. - A variance break is located to within days; a small mean break is almost invisible. Sweeping dates requires critical values of the maximum, not of a single date's F. ## Resources ### Research - [Volatility Targeting in Python](https://pythonparatrading.com/en/volatility-targeting-2/) · Research-0004 · 7 min ### Documentation - [hmmlearn: tutorial (GaussianHMM, multiple starts, predict_proba)](https://hmmlearn.readthedocs.io/en/stable/tutorial.html) · hmmlearn.readthedocs.io ### Books - [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. and 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 ### Data - [Level 2 data: series with Markov regimes and with a break, 2000-2025 (Parquet)](https://pythonparatrading.com/curso/datos/nivel-2/regimenes2.parquet) · Download · regimenes2.parquet Previous lesson: [2.1.6 · Kalman filter for dynamic hedge ratios](https://pythonparatrading.com/en/python-trading-course/level-2/kalman-filter-dynamic-hedge-ratio/) Next lesson: [2.1.8 · The multiple testing problem: Bonferroni, FDR and the “factor zoo”](https://pythonparatrading.com/en/python-trading-course/level-2/multiple-testing-bonferroni-fdr-factor-zoo/)