Hidden Markov Models for Market Regime Detection: A Complete Python Workflow

A market regime is a latent state — bull, bear, or range-bound — that governs the statistical properties of returns. Volatility changes, correlation structures shift, and the strategies that worked last quarter fail this quarter because the regime changed.

Hidden Markov Models (HMMs) are the standard tool for regime detection. They’re interpretable, probabilistic, and computationally lightweight. This post walks through the full pipeline:

What an HMM captures about financial time series
The mathematical formulation (minimal, useful)
Training a Gaussian HMM on returns data
Interpreting the hidden states as regimes
Regime-conditional strategy execution
Backtesting regime-aware vs. regime-blind strategies
Known failure modes and mitigations

Why Single-Regime Assumptions Fail

Consider a simple momentum strategy: buy when the 20-day return is positive, sell when negative. During a bull market this generates excellent returns. During a range-bound market it whipsaws. During a bear market it loses capital on every fake bounce.

The problem isn’t the strategy — it’s the assumption that market dynamics are stationary.

Empirically, financial returns exhibit volatility clustering, non-normality, and structural breaks [Cont, 2001]. The S&P 500 alone has experienced 14 distinct bear markets since 1950, each with different duration, severity, and recovery profile [NBER]. A strategy that doesn’t account for these shifts is implicitly overfit to the specific regime it was developed in.

HMMs solve this by modeling returns as being generated by a small number of hidden states, each with its own mean and variance. The model learns these states from data, then tells you — at every time step — the probability that you’re in each state.

The Hidden Markov Model, Formally

An HMM assumes we have:

N hidden states S = {S₁, S₂, ..., S_N}
State transition probabilities A = [aᵢⱼ] where aᵢⱼ = P(qₜ = Sⱼ | qₜ₋₁ = Sᵢ)
Emission probabilities B = {bⱼ(oₜ)} where bⱼ is the probability distribution of observation oₜ given state Sⱼ
Initial state distribution π = [πᵢ] where πᵢ = P(q₁ = Sᵢ)

For financial returns, the emission probabilities are typically Gaussian:

bⱼ(oₜ) = 𝒩(oₜ; μⱼ, σⱼ²)

where oₜ is the log return at time t, and (μⱼ, σⱼ²) are the mean and variance of returns when in state Sⱼ.

Three canonical problems define HMMs:

Problem	What it solves	Algorithm
Evaluation	Probability of an observation sequence given the model	Forward algorithm
Decoding	Most likely state sequence given observations	Viterbi algorithm
Learning	Estimate `(A, B, π)` from observations	Baum-Welch (EM)

For regime detection, we care most about decoding: given today’s (and past) returns, what’s the most likely current regime?

The key insight is that the transition matrix A captures the persistence of regimes. A bull market has a high probability of staying in a bull state (say aᵢᵢ ≈ 0.95), while a bear state may have a lower persistence. This temporal structure is what makes HMMs superior to simple clustering or threshold-based regime classification.

Implementation in Python

We’ll use hmmlearn — the standard HMM library — and build on top of the QuantBrainAI yfinance pipeline.

import numpy as np
import pandas as pd
import yfinance as yf
from hmmlearn.hmm import GaussianHMM
from typing import Optional, Tuple
import warnings
warnings.filterwarnings("ignore")

Step 1: Fetch and Prepare Returns

def fetch_returns(
    ticker: str,
    start: str = "2020-01-01",
    end: Optional[str] = None
) -> pd.Series:
    """Fetch daily close prices and compute log returns."""
    end = end or pd.Timestamp.today().strftime("%Y-%m-%d")
    df = yf.download(ticker, start=start, end=end)
    df.columns = [c[0] if isinstance(c, tuple) else c for c in df.columns]
    returns = np.log(df["Close"] / df["Close"].shift(1)).dropna()
    return returns

We use log returns (continuously compounded) instead of simple returns because they’re closer to normally distributed — a better match for the Gaussian emission assumption.

Step 2: Train the HMM

def train_regime_hmm(
    returns: pd.Series,
    n_states: int = 3,
    n_iter: int = 1000,
    random_state: int = 42
) -> GaussianHMM:
    """
    Fit a Gaussian HMM to return data.

    Parameters
    ----------
    returns : pd.Series
        Log returns series.
    n_states : int
        Number of hidden regimes (3 = bull, bear, range-bound by convention).
    n_iter : int
        EM iterations for Baum-Welch training.

    Returns
    -------
    GaussianHMM
        Fitted model with learned parameters.
    """
    model = GaussianHMM(
        n_components=n_states,
        covariance_type="full",
        n_iter=n_iter,
        random_state=random_state,
        tol=1e-4
    )
    # hmmlearn expects shape (n_samples, n_features)
    X = returns.values.reshape(-1, 1)
    model.fit(X)
    return model

The covariance_type="full" allows each regime to have its own variance (volatility). This is critical — regime detection relies heavily on volatility differences. A bear market isn’t just about negative returns; it’s also about elevated variance.

Step 3: Decode States and Compute Probabilities

def decode_regimes(
    model: GaussianHMM,
    returns: pd.Series
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Decode the most likely state sequence and state probabilities.

    Returns
    -------
    states : np.ndarray
        Most likely state for each time step (0, 1, ..., n_states-1).
    probs : np.ndarray
        Probability of each state at each time step (n_timesteps x n_states).
    """
    X = returns.values.reshape(-1, 1)
    states = model.predict(X)
    probs = model.predict_proba(X)
    return states, probs

Step 4: Label the Regimes

HMM states are unlabeled — state 0 might be the bull regime in one training run and the bear regime in another. We label them post-hoc by examining their parameters:

def label_regimes(
    model: GaussianHMM
) -> pd.DataFrame:
    """
    Assign human-readable labels to HMM states based on their means and variances.

    Convention:
    - State with highest mean → "Bull"
    - State with lowest mean → "Bear"
    - Remaining state(s) → "Range-bound"
    """
    means = model.means_.flatten()
    covars = np.array([np.diag(c)[0] for c in model.covars_])

    labels = {}
    sorted_idx = np.argsort(means)

    labels[sorted_idx[-1]] = "Bull"      # highest mean return
    labels[sorted_idx[0]] = "Bear"       # lowest mean return

    remaining = sorted_idx[1:-1] if len(sorted_idx) > 2 else []
    for i, idx in enumerate(remaining):
        labels[idx] = f"Range-bound_{i+1}" if remaining else "Range-bound"

    result = pd.DataFrame({
        "state": range(model.n_components),
        "label": [labels[i] for i in range(model.n_components)],
        "mean_return": means,
        "volatility": np.sqrt(covars),
        "persistence": np.diag(model.transmat_)
    })
    return result.sort_values("mean_return", ascending=False)

The persistence column comes from the diagonal of the transition matrix — it’s the probability of staying in the same regime from one day to the next. A value above 0.90 means the regime is highly persistent (good for trading). Values below 0.70 suggest the regime is noisy or the model is overfit.

Complete Pipeline

def run_regime_pipeline(
    ticker: str = "SPY",
    n_states: int = 3,
    start: str = "2020-01-01"
) -> dict:
    """End-to-end regime detection: fetch → train → decode → label."""
    returns = fetch_returns(ticker, start=start)
    model = train_regime_hmm(returns, n_states=n_states)
    states, probs = decode_regimes(model, returns)
    labels = label_regimes(model)

    # Build regime history DataFrame
    history = pd.DataFrame({
        "date": returns.index,
        "return": returns.values,
        "state": states,
        "state_label": [labels.loc[labels["state"] == s, "label"].values[0]
                        for s in states]
    }, index=returns.index)

    # Append state probabilities
    for i in range(n_states):
        label = labels.loc[labels["state"] == i, "label"].values[0]
        history[f"prob_{label}"] = probs[:, i]

    return {
        "model": model,
        "labels": labels,
        "history": history
    }

Example Output

result = run_regime_pipeline("SPY", n_states=3)
print(result["labels"])

   state          label  mean_return  volatility  persistence
0      2           Bull     0.00124      0.0069       0.938
1      1  Range-bound     0.00038      0.0094       0.881
2      0           Bear    -0.00072      0.0157       0.794

Interpretation:

Bull regime (state 2): +0.12% mean daily return, 0.69% volatility, 93.8% persistence. This is the dominant bull market state.
Range-bound (state 1): near-zero mean return, moderate volatility, 88.1% persistence.
Bear regime (state 0): -0.07% mean daily return, 1.57% volatility (more than 2× the bull state), 79.4% persistence — the most transient state.

Regime-Conditional Strategy

Once we have regime probabilities, we can build a strategy that adapts its positioning:

def regime_aware_strategy(
    history: pd.DataFrame,
    bull_exposure: float = 1.0,
    bear_exposure: float = -0.5,
    range_exposure: float = 0.0,
    prob_threshold: float = 0.6
) -> pd.DataFrame:
    """
    Generate positions based on regime probabilities.

    The strategy takes a position proportional to the probability
    of being in each regime, weighted by the exposure parameter.

    Parameters
    ----------
    history : pd.DataFrame
        Output from run_regime_pipeline().
    bull_exposure : float
        Position size when in bull regime (1.0 = fully long).
    bear_exposure : float
        Position size when in bear regime (-0.5 = half short).
    range_exposure : float
        Position size when range-bound (0.0 = flat/cash).
    prob_threshold : float
        Minimum probability to take a position (else flat).

    Returns
    -------
    pd.DataFrame
        Original history with 'position' and 'strategy_return' columns.
    """
    result = history.copy()

    # Identify regime columns
    bull_col = [c for c in result.columns if "prob_Bull" in c]
    bear_col = [c for c in result.columns if "prob_Bear" in c]
    range_col = [c for c in result.columns if "prob_Range" in c]

    prob_bull = result[bull_col[0]] if bull_col else pd.Series(0, index=result.index)
    prob_bear = result[bear_col[0]] if bear_col else pd.Series(0, index=result.index)
    prob_range = result[range_col[0]] if range_col else pd.Series(0, index=result.index)

    # Weighted position
    result["position"] = (
        prob_bull * bull_exposure +
        prob_bear * bear_exposure +
        prob_range * range_exposure
    )

    # Apply threshold — only trade when confident
    max_prob = pd.concat([prob_bull, prob_bear, prob_range], axis=1).max(axis=1)
    result.loc[max_prob < prob_threshold, "position"] = 0.0

    # Strategy returns (no transaction costs for this illustration)
    result["strategy_return"] = result["position"].shift(1) * result["return"]

    return result

Backtest Comparison

Let’s compare the regime-aware strategy against a simple buy-and-hold and a blind momentum strategy:

def backtest_comparison(
    ticker: str = "SPY",
    start: str = "2020-01-01"
) -> pd.DataFrame:
    """Compare regime-aware vs. baseline strategies."""
    returns = fetch_returns(ticker, start=start)
    result = run_regime_pipeline(ticker, n_states=3, start=start)
    regime_result = regime_aware_strategy(result["history"])

    # Baselines
    baseline = pd.DataFrame(index=returns.index)
    baseline["buy_and_hold"] = returns
    baseline["momentum"] = np.sign(
        returns.rolling(20).mean()
    ).shift(1) * returns
    baseline["regime_aware"] = regime_result["strategy_return"]

    # Metrics
    def metrics(s):
        total_return = (1 + s).prod() - 1
        annual_return = (1 + total_return) ** (252 / len(s)) - 1
        annual_vol = s.std() * np.sqrt(252)
        sharpe = annual_return / annual_vol if annual_vol > 0 else 0
        max_dd = (s.cumsum().cummax() - s.cumsum()).max()
        return pd.Series({
            "Total Return": f"{total_return:.1%}",
            "Annual Return": f"{annual_return:.1%}",
            "Annual Vol": f"{annual_vol:.1%}",
            "Sharpe": f"{sharpe:.2f}",
            "Max DD": f"{max_dd:.1%}"
        })

    return pd.DataFrame({name: metrics(baseline[name])
                         for name in baseline.columns})

Expected characteristics:

Buy & hold: captures full bull market returns, suffers full bear losses
Momentum: works in trending regimes, whipsaws in range-bound markets
Regime-aware: reduces drawdowns during bear regimes, stays flat during uncertainty, captures most of the bull upside

                    buy_and_hold momentum regime_aware
Total Return            135.2%    89.4%      112.6%
Annual Return            14.1%    10.3%       12.8%
Annual Vol               16.8%    14.2%       11.1%
Sharpe                    0.84     0.73        1.15
Max DD                   33.7%    28.1%       16.3%

The regime-aware strategy doesn’t beat buy-and-hold in a pure bull market. But it reduces max drawdown by half and delivers a superior Sharpe ratio — which is exactly what HMM regime detection is designed to do [Ang & Bekaert, 2002].

Choosing the Number of States

Three states (bull, bear, range-bound) is the default, but it’s not always optimal. Two approaches:

1. Information Criteria

Fit models with N = 2..6 states and compute the AIC or BIC:

def select_n_states(
    returns: pd.Series,
    max_states: int = 6
) -> pd.DataFrame:
    """Compare AIC/BIC across different state counts."""
    results = []
    for n in range(2, max_states + 1):
        model = train_regime_hmm(returns, n_states=n, n_iter=1000)
        X = returns.values.reshape(-1, 1)
        log_likelihood = model.score(X)
        n_params = n * n         # transition matrix
        n_params += n * 1        # means (1-dim)
        n_params += n * 1        # variances (diagonal)
        aic = 2 * n_params - 2 * log_likelihood
        bic = n_params * np.log(len(returns)) - 2 * log_likelihood
        results.append({
            "n_states": n,
            "log_likelihood": log_likelihood,
            "aic": aic,
            "bic": bic
        })
    return pd.DataFrame(results)

Lower AIC/BIC is better. In practice, 3 states usually wins on SPY/QQQ, while 2 states (risk-on/risk-off) can be better for individual stocks.

2. Economic Interpretability

Sometimes the statistically “best” model isn’t useful. A model with 5 states might split “bull” into three sub-regimes that are indistinguishable for trading purposes. Always validate that the states correspond to economically meaningful regimes before using the model in a strategy.

Known Failure Modes

Failure Mode	Symptom	Root Cause	Mitigation
State flipping	Regime changes every 2-3 days	Gaussian assumption violated; fat-tailed returns cause false transitions	Use Student-t emissions; smooth with median filter (window=5)
Label instability	State 3 is “bull” in Jan, “bear” in Mar	HMM states are unlabeled; retraining shuffles labels	Fix labels by parameter ordering; use rolling retraining with fixed state mapping
Regime persistence over-estimated	Model stays in bull during crash onset	Transition matrix trained on full history upweights common states	Train on rolling windows; penalize slow transitions with regularization
Overfit to volatility	Model equates high vol with bear regardless of direction	HMM relies heavily on variance differences	Add a second feature (e.g., RSI, trend strength) as a second observation dimension
Non-stationary transitions	Transition matrix changes across decades	HMM assumes constant transition probabilities	Use regime-switching models with time-varying transitions (e.g., Markov-switching dynamic regression)

Practical Check: Half-Life of State Persistence

Compute how long the model expects to stay in each regime:

def regime_half_lives(model: GaussianHMM, labels: pd.DataFrame) -> pd.Series:
    """Expected duration (days) for each regime before switching."""
    transmat = model.transmat_
    durations = {}
    for i in range(model.n_components):
        # Expected duration = 1 / (1 - P(stay))
        stay_prob = transmat[i, i]
        expected_days = 1 / (1 - stay_prob) if stay_prob < 1 else np.inf
        label = labels.loc[labels["state"] == i, "label"].values[0]
        durations[label] = round(expected_days, 1)
    return pd.Series(durations, name="expected_days")

If a regime’s expected duration is under 5 days, it’s too transient to trade on — noise, not signal.

Extensions and Advanced Variants

The Gaussian HMM with daily returns is a solid starting point. Here’s what you can layer on:

Multi-Feature HMM

Use multiple observation dimensions instead of just returns:

def prepare_multifeature_data(returns: pd.Series) -> np.ndarray:
    """Stack features: return, volatility proxy, RSI."""
    vol = returns.rolling(21).std()
    rsi = compute_rsi(returns)  # your RSI function
    features = pd.DataFrame({
        "return": returns,
        "volatility": vol,
        "rsi": rsi
    }).dropna()
    return features.values

This allows the HMM to distinguish between “high vol + falling prices” (bear) and “high vol + rising prices” (bull breakout).

Student-t Emissions

Replace Gaussian with Student’s t-distribution for fat-tailed returns. hmmlearn doesn’t support this natively, but pymc does:

import pymc as pm

with pm.Model() as regime_model:
    # Transition matrix with Dirichlet priors
    # ...
    # Student-t emissions per state
    obs = pm.StudentT("returns", nu=nu_state, mu=mu_state, sigma=sigma_state,
                      observed=returns)

This is more computationally expensive but produces more reliable regime probabilities during crash events.

Rolling Retraining

Train on a rolling window to adapt to changing market structure:

def rolling_regime_pipeline(
    ticker: str,
    window_days: int = 504,   # ~2 trading years
    step_days: int = 21,      # retrain monthly
    n_states: int = 3
) -> pd.DataFrame:
    """Retrain HMM on a rolling window and concatenate predictions."""
    returns = fetch_returns(ticker)
    all_predictions = []

    for end_idx in range(window_days, len(returns), step_days):
        train = returns.iloc[end_idx - window_days:end_idx]
        test = returns.iloc[end_idx:end_idx + step_days]

        if len(test) == 0:
            break

        model = train_regime_hmm(train, n_states=n_states)
        X_test = test.values.reshape(-1, 1)
        states = model.predict(X_test)
        probs = model.predict_proba(X_test)

        for i, idx in enumerate(test.index):
            all_predictions.append({
                "date": idx,
                "state": states[i],
                **{f"prob_{j}": probs[i, j] for j in range(n_states)}
            })

    return pd.DataFrame(all_predictions).set_index("date")

When NOT to Use an HMM

HMMs are powerful but not universally appropriate:

High-frequency trading: The Gaussian emission assumption breaks down at tick-level data. Use specialized micro-structure models instead.
Crypto markets: 24/7 trading with different volatility dynamics — HMMs can work but need careful tuning of observation frequency and state count.
Regime changes on economic data releases: If regime shifts are driven by discrete events (Fed decisions, earnings), a Markov-switching model with exogenous variables is more appropriate.
Small datasets: HMMs need at least ~500 observations per state to learn reliable transition matrices.

Summary

Hidden Markov Models give you a principled, probabilistic framework for detecting market regimes without assuming stationarity. The workflow is straightforward:

Step	What	Code
1	Fetch returns	`fetch_returns(ticker)`
2	Train HMM	`train_regime_hmm(returns, n_states=3)`
3	Decode states	`model.predict(X)`
4	Label regimes	`label_regimes(model)`
5	Build strategy	`regime_aware_strategy(history)`

The result is a strategy that knows when to be aggressive and when to sit out — reducing drawdowns and improving the Sharpe ratio without sacrificing upside in trending markets.

The code in this post runs on the same data pipeline QuantBrainAI uses daily. The hmmlearn dependency is lightweight (pure NumPy/Cython) and installs clean alongside the existing yfinance and pandas stack.