Regime Detection (`wraquant.regimes`)¶

fit_gaussian_hmm(returns, n_states=2, covariance_type='full', n_init=10, n_iter=200, tol=0.0001, random_state=42)[source]¶

Fit a Gaussian Hidden Markov Model to financial return data.

A Gaussian HMM assumes that observed returns are generated by a latent (hidden) Markov chain with K states. Each state k has its own Gaussian emission distribution N(mu_k, sigma_k). The model learns:

Transition matrix A[i,j] = P(s_t = j | s_{t-1} = i)
Emission means mu_k for each state
Emission covariances Sigma_k for each state
Initial state distribution pi_k = P(s_0 = k)

The EM algorithm (Baum-Welch) is used for parameter estimation. Because EM can converge to local optima, multiple random restarts (n_init) are performed and the model with the highest log-likelihood is retained.

Typical interpretation for 2-state model on equity returns:

Low-volatility state (“risk-on” / bull): Positive mean return, low standard deviation, high Sharpe ratio. Markets are calm.
High-volatility state (“risk-off” / crisis): Near-zero or negative mean return, high standard deviation. Markets are stressed.

States are automatically re-ordered so that state 0 has the lowest volatility and state K-1 has the highest volatility. This ensures consistent labeling across runs.

Parameters:

returns (Series | ndarray) – Return series (simple or log returns). If a pandas Series, NaN values are dropped and the index is preserved.
n_states (int, default: 2) – Number of hidden states. Common choices are 2 (bull/bear) or 3 (bull/neutral/bear). Values above 5 risk overfitting.
covariance_type (str, default: 'full') – Type of covariance matrix. One of 'full', 'diag', 'spherical', 'tied'.
n_init (int, default: 10) – Number of random restarts for EM. The model with the highest log-likelihood is selected. More restarts reduce the chance of a poor local optimum.
n_iter (int, default: 200) – Maximum number of EM iterations per restart.
tol (float, default: 0.0001) – Convergence threshold for EM.
random_state (int, default: 42) – Base random seed for reproducibility.

Returns:

states (np.ndarray): Most likely state sequence (Viterbi), shape (T,), integers in [0, K-1].
state_probs (np.ndarray): Smoothed state probabilities, shape (T, K). Each row sums to 1.
transition_matrix (np.ndarray): Transition probability matrix, shape (K, K). Row i sums to 1.
means (np.ndarray): Emission means, shape (K,).
covariances (np.ndarray): Emission covariances, shape (K,) (variance per state for univariate data).
startprob (np.ndarray): Initial state distribution, shape (K,).
log_likelihood (float): Total log-likelihood of the data.
aic (float): Akaike Information Criterion.
bic (float): Bayesian Information Criterion.
n_states (int): Number of states.
steady_state (np.ndarray): Ergodic (stationary) distribution of the Markov chain, shape (K,).
avg_duration (np.ndarray): Expected duration in each state, 1 / (1 - A[i,i]), shape (K,).
model (GaussianHMM): The fitted hmmlearn model object.
index (pd.Index | None): Original index if input was a pandas Series.

Return type:

Example

>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(42)
>>> # Simulate 2-regime returns
>>> low_vol = rng.normal(0.001, 0.008, 300)
>>> high_vol = rng.normal(-0.002, 0.025, 200)
>>> returns = pd.Series(np.concatenate([low_vol, high_vol]))
>>> result = fit_gaussian_hmm(returns, n_states=2)
>>> print(result['means'])        # Per-state means
>>> print(result['transition_matrix'])  # Transition probabilities
>>> print(result['steady_state'])  # Long-run regime probabilities

Notes

The model uses hmmlearn’s GaussianHMM. For very long series (>10,000 observations), consider downsampling or using the variational Bayes variant.

See also

fit_ms_regression: Markov-switching regression with statsmodels. gaussian_mixture_regimes: GMM-based regime classification. regime_statistics: Compute per-regime financial statistics. regime_transition_analysis: Analyze regime transition dynamics. wraquant.regimes.base.detect_regimes: Unified high-level interface. wraquant.backtest.position.regime_signal_filter: Gate signals by regime. wraquant.viz.dashboard.regime_dashboard: Visualize regime analysis.

fit_hmm(returns, n_states=2)[source]¶

Fit a Gaussian HMM to return data (legacy interface).

This is a thin compatibility wrapper around fit_gaussian_hmm. New code should use fit_gaussian_hmm directly for richer output.

Parameters:

returns (Series) – Simple return series.
n_states (int, default: 2) – Number of hidden states.

Return type:

Any

Returns:

Fitted hmmlearn.hmm.GaussianHMM model.

See also

fit_gaussian_hmm: Full-featured HMM fitting with rich output.

fit_ms_autoregression(endog, k_regimes=2, order=1, switching_ar=True, switching_variance=True, exog_tvtp=None)[source]¶

Fit a Markov-switching autoregression (MS-AR) model.

This wraps statsmodels.tsa.regime_switching.MarkovAutoregression, which extends the standard Markov-switching regression to include autoregressive lags whose coefficients can also switch across regimes.

Hamilton (1989) framework:

The seminal Hamilton (1989) paper introduced the Markov-switching autoregressive model to study business-cycle dynamics. The key insight is that the data-generating process can shift between distinct regimes (e.g., expansion vs. recession) governed by a latent Markov chain, and within each regime the series follows an AR(p) process with regime-specific parameters:

\[y_t = \mu_{s_t} + \phi_{1,s_t}\,(y_{t-1} - \mu_{s_{t-1}}) + \cdots + \phi_{p,s_t}\,(y_{t-p} - \mu_{s_{t-p}}) + \sigma_{s_t}\,\epsilon_t\]

When to use MS-AR vs. plain HMM:

Use MS-AR when the series exhibits serial correlation within each regime (e.g., GDP growth, interest rates, momentum returns). The autoregressive terms capture local dynamics that a plain HMM ignores.
Use a plain Gaussian HMM (fit_gaussian_hmm) when returns are approximately i.i.d. within each regime and only the mean and variance shift. This is typical for daily equity returns.
Use MS-Regression (fit_ms_regression) when you want regime-dependent intercepts/exogenous coefficients but no AR lags.

Parameters:

endog (Series | ndarray) – Endogenous (dependent) variable. Typically a return or growth-rate series.
k_regimes (int, default: 2) – Number of regimes (hidden states). 2 is the most common choice (e.g., expansion/recession).
order (int, default: 1) – Autoregressive lag order p. Higher orders capture longer memory within each regime.
switching_ar (bool, default: True) – If True, the AR coefficients switch across regimes. If False, only the intercept and/or variance switch.
switching_variance (bool, default: True) – If True, the innovation variance switches across regimes. Almost always desirable for financial data.
exog_tvtp (ndarray | DataFrame | None, default: None) – Optional exogenous variables for time-varying transition probabilities (TVTP). Shape (T, q). When provided, the transition probabilities are modelled as logistic functions of these variables, allowing macro indicators to influence regime-switching dynamics.

Returns:

smoothed_probs (np.ndarray): Smoothed regime probabilities, shape (T, K).
filtered_probs (np.ndarray): Filtered regime probabilities, shape (T, K).
states (np.ndarray): Most likely state sequence (argmax of smoothed probabilities), shape (T,).
transition_matrix (np.ndarray): Transition matrix, shape (K, K).
expected_durations (np.ndarray): Expected duration in each regime, 1 / (1 - P[k,k]), shape (K,).
regime_params (dict): Per-regime parameters including 'mean_k', 'sigma2_k', and 'ar_k_j' (AR coefficient j for regime k).
aic (float): Akaike Information Criterion.
bic (float): Bayesian Information Criterion.
summary (str): Text summary from statsmodels.
model_result: The fitted statsmodels result object.

Return type:

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> # Simulate AR(1) with switching mean and persistence
>>> returns = pd.Series(np.concatenate([
...     rng.normal(0.002, 0.008, 250),
...     rng.normal(-0.001, 0.020, 250),
... ]))
>>> result = fit_ms_autoregression(returns, k_regimes=2, order=1)
>>> print(result['transition_matrix'])
>>> print(result['regime_params'])

Notes

Requires statsmodels >= 0.13. Convergence can be slow for higher-order AR models or many regimes. The Hamilton filter is used for inference.

References

Hamilton, J. D. (1989). “A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle.” Econometrica, 57(2), 357-384.

See also

fit_ms_regression: Markov-switching regression (no AR terms). fit_gaussian_hmm: HMM without autoregressive dynamics. select_n_states: BIC-based state selection.

fit_ms_regression(endog, k_regimes=2, trend='c', switching_variance=True, exog=None, n_iter=200)[source]¶

Fit a Markov-switching regression model using statsmodels.

This wraps statsmodels.tsa.regime_switching.MarkovRegression, which estimates a linear regression where the intercept (and optionally variance) can switch between regimes according to a hidden Markov chain.

When to use: When you want to model how the data-generating process changes across regimes, including regime-dependent means and volatilities. More interpretable than a raw HMM for economic modeling because you get regression coefficients per regime.

The model:

\[y_t = \mu_{s_t} + x_t' \beta_{s_t} + \sigma_{s_t} \epsilon_t\]

where \(s_t \in \{0, 1, \ldots, K-1\}\) follows a first-order Markov chain with transition matrix \(P\).

Parameters:

endog (Series | ndarray) – Endogenous (dependent) variable. Typically a return series.
k_regimes (int, default: 2) – Number of regimes (hidden states).
trend (str, default: 'c') – Trend specification for the mean equation. 'c' (constant), 't' (linear trend), 'ct' (both), 'n' (no trend).
switching_variance (bool, default: True) – If True, the variance switches across regimes. Almost always desirable for financial data.
exog (DataFrame | ndarray | None, default: None) – Optional exogenous regressors, shape (T, p). The coefficients on these regressors are not switching by default (the intercept and variance switch).
n_iter (int, default: 200) – Maximum number of EM iterations.

Returns:

smoothed_probs (np.ndarray): Smoothed regime probabilities, shape (T, K). smoothed_probs[:, k] is the probability of being in regime k at each time step given all data.
filtered_probs (np.ndarray): Filtered regime probabilities (using only past data), shape (T, K).
states (np.ndarray): Most likely state sequence (argmax of smoothed probabilities), shape (T,).
transition_matrix (np.ndarray): Transition matrix, shape (K, K).
regime_params (dict): Regime-specific parameters. Keys are 'mean_0', 'mean_1', …, 'sigma2_0', etc.
log_likelihood (float): Log-likelihood of the fitted model.
aic (float): Akaike Information Criterion.
bic (float): Bayesian Information Criterion.
summary (str): Text summary from statsmodels.
model_result (MarkovRegressionResults): The full statsmodels result object.

Return type:

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(np.concatenate([
...     rng.normal(0.001, 0.01, 200),
...     rng.normal(-0.002, 0.03, 200),
... ]))
>>> result = fit_ms_regression(returns, k_regimes=2)
>>> print(result['transition_matrix'])
>>> print(result['regime_params'])

Notes

Requires statsmodels >= 0.13. The model uses Hamilton’s (1989) filter for inference. Convergence can be sensitive to starting values.

See also

fit_gaussian_hmm: HMM using hmmlearn (more flexible emission: distributions).

regime_statistics: Compute financial statistics per regime.

fit_multivariate_hmm(returns, n_states=2, covariance_type='full', n_init=10, n_iter=200, random_state=42)[source]¶

Fit a multivariate Gaussian HMM for multi-asset regime detection.

Extends the univariate HMM to a vector of asset returns, enabling detection of regimes that manifest as shifts in cross-asset correlations, not just individual-asset volatility.

Why multivariate regime detection matters:

Single-asset HMMs detect volatility regimes but miss a crucial phenomenon: correlation regime shifts. During crises, asset correlations spike toward 1, destroying diversification precisely when it is needed most. A multivariate HMM captures this by estimating a full covariance matrix per regime.

Typical findings on equity/bond data:

Low-vol regime: Low equity vol, low bond vol, moderate negative equity-bond correlation (diversification works).
High-vol regime: High equity vol, elevated bond vol, correlations shift (equities become more correlated with each other, equity-bond correlation may flip sign).

The model estimates per-regime:

Mean return vector \(\mu_k \in \mathbb{R}^d\)
Full covariance matrix \(\Sigma_k \in \mathbb{R}^{d \times d}\)
From which per-regime correlation matrices are derived.

Parameters:

returns (DataFrame) – DataFrame of multi-asset returns, shape (T, d) where d is the number of assets. Column names are preserved for labeling. NaN rows are dropped.
n_states (int, default: 2) – Number of hidden states (regimes).
covariance_type (str, default: 'full') – Covariance parameterisation. 'full' is strongly recommended for capturing correlation shifts.
n_init (int, default: 10) – Number of random EM restarts.
n_iter (int, default: 200) – Maximum EM iterations per restart.
random_state (int, default: 42) – Random seed.

Returns:

states (np.ndarray): Most likely state sequence (Viterbi), shape (T,).
state_probs (np.ndarray): Smoothed state probabilities, shape (T, K).
means (np.ndarray): Per-regime mean vectors, shape (K, d).
covariances (np.ndarray): Per-regime covariance matrices, shape (K, d, d).
per_regime_correlations (dict[int, np.ndarray]): Per-regime correlation matrices, keyed by regime index.
transition_matrix (np.ndarray): Transition matrix (K, K).
log_likelihood (float): Total log-likelihood.
aic (float): AIC.
bic (float): BIC.
model: The fitted hmmlearn model.
index (pd.Index | None): Preserved DatetimeIndex.
columns (list[str]): Asset names.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> # Low-vol regime: low correlation
>>> cov_low = np.array([[0.0001, 0.00002], [0.00002, 0.00015]])
>>> low = rng.multivariate_normal([0.001, 0.0008], cov_low, 300)
>>> # High-vol regime: high correlation
>>> cov_high = np.array([[0.0009, 0.0006], [0.0006, 0.0008]])
>>> high = rng.multivariate_normal([-0.001, -0.0005], cov_high, 200)
>>> df = pd.DataFrame(np.vstack([low, high]), columns=['SPY', 'QQQ'])
>>> result = fit_multivariate_hmm(df, n_states=2)
>>> print(result['per_regime_correlations'][0])
>>> print(result['per_regime_correlations'][1])

See also

fit_gaussian_hmm: Univariate HMM. regime_conditional_moments: Per-regime moments from states.

gaussian_mixture_regimes(returns, n_components=2, covariance_type='full', n_init=10, random_state=42)[source]¶

Classify returns into regimes using a Gaussian Mixture Model.

Unlike HMMs, GMMs do not model temporal dependence between states. Each observation is independently assigned to the most likely Gaussian component. This makes GMMs:

Faster to fit than HMMs.
Simpler to interpret (no transition dynamics).
Less powerful for capturing regime persistence.

Good for quick regime labeling or as a baseline before fitting an HMM.

The model:

\[p(r_t) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(r_t; \mu_k, \Sigma_k)\]

Parameters:

returns (Series | ndarray) – Return series. NaN values are dropped.
n_components (int, default: 2) – Number of mixture components (regimes).
covariance_type (str, default: 'full') – Covariance type ('full', 'tied', 'diag', 'spherical').
n_init (int, default: 10) – Number of random restarts.
random_state (int, default: 42) – Random seed.

Returns:

states (np.ndarray): Component assignments, shape (T,).
state_probs (np.ndarray): Posterior probabilities, shape (T, K).
means (np.ndarray): Component means, shape (K,).
covariances (np.ndarray): Component variances, shape (K,).
weights (np.ndarray): Mixing proportions, shape (K,).
aic (float): AIC.
bic (float): BIC.
model (GaussianMixture): Fitted sklearn model.
index (pd.Index | None): Preserved index.

Return type:

Example

>>> result = gaussian_mixture_regimes(returns, n_components=2)
>>> print(result['means'])   # Component means
>>> print(result['weights']) # Mixing proportions
>>> # Combine with regime_statistics for financial interpretation
>>> stats = regime_statistics(returns, result['states'])

See also

fit_gaussian_hmm: HMM with temporal regime dynamics. regime_statistics: Per-regime financial statistics.

predict_regime(model, returns)[source]¶

Predict regime states and probabilities for new data.

Given a fitted HMM (from fit_gaussian_hmm or hmmlearn directly), decode the most likely hidden state sequence and compute posterior state probabilities for new observations.

Parameters:

model (Any) – A fitted HMM model with predict and predict_proba methods (e.g., from hmmlearn), or the dict returned by fit_gaussian_hmm (in which case the 'model' key is used).
returns (Series | ndarray) – Return series to classify. NaN values are dropped.

Returns:

states (np.ndarray): Most likely state at each time step.
state_probs (np.ndarray): Posterior probabilities, shape (T, K).
index (pd.Index | None): Preserved index if input was a pandas Series.

Return type:

Example

>>> result = fit_gaussian_hmm(train_returns, n_states=2)
>>> pred = predict_regime(result, test_returns)
>>> print(pred['states'])
>>> print(pred['state_probs'][:5])

regime_aware_portfolio(regime_probs, regime_weights)[source]¶

Compute blended portfolio weights using regime probabilities.

Given the current probability of being in each regime and the optimal portfolio weights for each regime, compute a probability-weighted blend. This bridges regime detection (regimes/) with portfolio optimisation (opt/).

The blended weight for asset j is:

\[w_j = \sum_{k=1}^{K} P(s_t = k) \cdot w_j^{(k)}\]

Parameters:

regime_probs (ndarray) – Current regime probabilities, shape (K,). Must sum to 1.
regime_weights (ndarray) – Optimal weights for each regime, shape (K, n_assets). Row k contains the weights for regime k.

Returns:

Blended portfolio weights, shape (n_assets,). Sums to 1 if each row of regime_weights sums to 1.

Return type:

ndarray

Example

>>> # 60% probability of bull, 40% probability of bear
>>> probs = np.array([0.6, 0.4])
>>> # Bull: 80% equity / 20% bonds; Bear: 30% equity / 70% bonds
>>> weights = np.array([[0.8, 0.2], [0.3, 0.7]])
>>> blended = regime_aware_portfolio(probs, weights)
>>> print(blended)  # [0.6*0.8 + 0.4*0.3, 0.6*0.2 + 0.4*0.7]
[0.6  0.4]

See also

fit_gaussian_hmm: Obtain regime probabilities. rolling_regime_probability: Time-varying regime probabilities.

regime_conditional_moments(returns, states)[source]¶

Compute per-regime mean vector and covariance/correlation matrices.

After fitting a multivariate HMM (or assigning regimes by any method), this function extracts the sample moments for each regime. These moments are the building blocks for regime-aware portfolio construction: the mean vector feeds into the expected-return input and the covariance matrix replaces the unconditional covariance in a mean-variance optimiser.

Parameters:

returns (DataFrame) – DataFrame of multi-asset returns, shape (T, d).
states (ndarray) – Integer regime labels, shape (T,), aligned with returns.

Return type:

dict[int, dict[str, ndarray]]

Returns:

Dictionary keyed by regime index (int), where each value is a dict with:

mean (np.ndarray): Sample mean vector, shape (d,).
cov (np.ndarray): Sample covariance matrix, shape (d, d).
corr (np.ndarray): Sample correlation matrix, shape (d, d).

Example

>>> result = fit_multivariate_hmm(returns_df, n_states=2)
>>> moments = regime_conditional_moments(
...     returns_df, result['states']
... )
>>> # Use regime-0 covariance for portfolio optimisation
>>> cov_bull = moments[0]['cov']
>>> mu_bull = moments[0]['mean']

See also

fit_multivariate_hmm: Multi-asset regime detection. regime_aware_portfolio: Blend weights across regimes.

regime_statistics(returns, states)[source]¶

Compute descriptive statistics for each detected regime.

After fitting an HMM or GMM, this function summarizes the financial characteristics of each regime. Typical interpretation for a 2-state model on equity returns:

Regime 0 (low vol): “risk-on” or “bull market” – positive mean return, low volatility, high Sharpe ratio.
Regime 1 (high vol): “risk-off” or “crisis” – negative or zero mean, high volatility, negative Sharpe.

Parameters:

returns (Series | ndarray) – Return series aligned with states. Can be a pandas Series or a numpy array.
states (ndarray) – Integer regime labels (0, 1, …, K-1), same length as returns.

Returns:

mean, std, sharpe, sortino_ratio, min, max, skewness, kurtosis, max_drawdown, VaR_95, CVaR_95, n_observations, pct_time, avg_duration.

sharpe: Annualized Sharpe ratio assuming 252 trading days.
sortino_ratio: Annualized Sortino ratio (downside deviation only, 252 trading days).
max_drawdown: Maximum peak-to-trough drawdown within the regime’s return segments.
VaR_95: 5th percentile of returns (historical Value at Risk at the 95 % confidence level).
CVaR_95: Mean of returns below the VaR_95 threshold (Conditional Value at Risk / Expected Shortfall).
avg_duration: Mean consecutive time steps spent in the regime.

Return type:

Example

>>> result = fit_gaussian_hmm(returns, n_states=2)
>>> stats = regime_statistics(returns, result['states'])
>>> print(stats)
            mean       std    sharpe  ...  avg_duration
regime
0       0.000523  0.008123  1.0223   ...          15.3
1      -0.001842  0.024561 -1.1917   ...           7.8

See also

fit_gaussian_hmm: Fit HMM to obtain state assignments. wraquant.regimes.base.regime_report: Comprehensive regime analysis. wraquant.backtest.position.regime_conditional_sizing: Adjust positions by regime. regime_transition_analysis: Analyze transition dynamics.

regime_transition_analysis(states, transition_matrix=None)[source]¶

Analyze regime transition dynamics.

Computes empirical and model-based statistics about how regimes switch. Useful for understanding regime persistence and expected timing of regime changes.

Parameters:

states (ndarray) – Integer regime labels, shape (T,).
transition_matrix (ndarray | None, default: None) – Model-estimated transition matrix, shape (K, K). If None, the empirical transition matrix is computed from states.

Returns:

empirical_transition_matrix (np.ndarray): Transition counts normalised to probabilities, shape (K, K).
transition_matrix (np.ndarray): Model or empirical transition matrix used for analysis.
avg_duration (np.ndarray): Expected duration in each regime (time steps), 1 / (1 - P[i,i]).
steady_state (np.ndarray): Ergodic distribution from the transition matrix.
regime_counts (dict): Number of distinct “visits” to each regime.
regime_durations (dict): List of consecutive durations for each regime.
expected_return_time (np.ndarray): Expected number of steps to return to each regime from stationarity.

Return type:

Example

>>> result = fit_gaussian_hmm(returns, n_states=2)
>>> analysis = regime_transition_analysis(
...     result['states'], result['transition_matrix']
... )
>>> print(f"Bull market lasts ~{analysis['avg_duration'][0]:.0f} days")
>>> print(f"Steady state: {analysis['steady_state']}")

See also

fit_gaussian_hmm: Fit HMM to obtain states and transition matrix. regime_statistics: Financial statistics per regime.

rolling_regime_probability(returns, n_states=2, window=None, min_window=60, n_init=5)[source]¶

Compute time-varying regime probabilities using rolling/expanding HMM.

Fits a Gaussian HMM at each time step using either a rolling or expanding window, producing a time series of regime probabilities that adapts to new data. Useful for real-time regime monitoring.

Parameters:

returns (Series) – Return series.
n_states (int, default: 2) – Number of hidden states.
window (int | None, default: None) – Rolling window size. If None, uses an expanding window (all data up to time t).
min_window (int, default: 60) – Minimum number of observations required before fitting the first model.
n_init (int, default: 5) – Number of EM restarts per window.

Return type:

Returns:

pd.DataFrame with columns prob_0, prob_1, …, prob_{K-1} indexed by the return series index. Rows before min_window are NaN.

Example

>>> probs = rolling_regime_probability(returns, n_states=2, window=252)
>>> probs.plot(title="Regime Probabilities Over Time")

Notes

This is computationally expensive. For 1000 observations with window=252, approximately 750 HMM fits are performed. Consider using a large window or expanding mode for efficiency.

See also

fit_gaussian_hmm: Single-shot HMM fitting.

select_n_states(returns, max_states=5, n_init=10, n_iter=200, random_state=42)[source]¶

Select the optimal number of HMM states using the BIC criterion.

Fits Gaussian HMMs with n_states from 2 to max_states and selects the model that minimises the Bayesian Information Criterion (BIC).

Why BIC over AIC for HMM state selection:

The BIC imposes a stronger penalty on model complexity than the AIC, which is critical for HMM selection because:

HMM parameter count grows quadratically with the number of states (transition matrix has K*(K-1) free parameters), so without a strong penalty the AIC tends to overfit by selecting too many states.
BIC is consistent: as sample size grows, BIC selects the true model (if it is in the candidate set). AIC tends to overfit even asymptotically.
In practice, financial regime models with 2-3 states are usually adequate; BIC correctly penalises 4+ state models that capture noise rather than true regimes.

The BIC is defined as:

\[\text{BIC} = k \ln(T) - 2 \ln(\hat{L})\]

where k is the number of free parameters, T is the sample size, and \(\hat{L}\) is the maximised likelihood.

Parameters:

returns (Series | ndarray) – Return series. NaN values are dropped.
max_states (int, default: 5) – Maximum number of states to consider (inclusive). The search covers [2, 3, ..., max_states].
n_init (int, default: 10) – Number of random restarts per model fit.
n_iter (int, default: 200) – Maximum EM iterations per restart.
random_state (int, default: 42) – Base random seed.

Returns:

optimal_n_states (int): The number of states with the lowest BIC.
scores (dict[int, float]): Mapping from n_states to its BIC value.
best_model (dict): The full result dict from fit_gaussian_hmm for the optimal model.

Return type:

Example

>>> result = select_n_states(returns, max_states=5)
>>> print(f"Optimal states: {result['optimal_n_states']}")
>>> print(result['scores'])
{2: -2345.6, 3: -2310.1, 4: -2290.3, 5: -2275.0}

See also

fit_gaussian_hmm: Fit a single HMM with a fixed number of states.

kalman_filter(observations, F, H, Q, R, x0, P0)[source]¶

Run a linear Kalman filter over a sequence of observations.

The Kalman filter is the optimal recursive estimator for a linear Gaussian state-space model:

\[ \begin{align}\begin{aligned}x_t = F \, x_{t-1} + w_t, \quad w_t \sim N(0, Q)\\y_t = H \, x_t + v_t, \quad v_t \sim N(0, R)\end{aligned}\end{align} \]

At each time step, the filter produces:

Prediction: \(\hat{x}_{t|t-1}\) and \(P_{t|t-1}\)
Update: \(\hat{x}_{t|t}\) and \(P_{t|t}\) after incorporating observation \(y_t\)

Financial applications:

Estimating latent factors or unobserved states
Noise reduction (e.g., estimating “true” volatility from noisy proxies)
Preprocessing for regime detection

Parameters:

observations (ndarray) – Observation matrix of shape (T, m) where T is the number of time steps and m is the observation dimension. For univariate observations, shape (T, 1).
F (ndarray) – State transition matrix (n, n).
H (ndarray) – Observation matrix (m, n).
Q (ndarray) – Process noise covariance (n, n).
R (ndarray) – Observation noise covariance (m, m).
x0 (ndarray) – Initial state estimate (n,) or (n, 1).
P0 (ndarray) – Initial state covariance (n, n).

Returns:

filtered_states (np.ndarray): Filtered state estimates, shape (T, n).
filtered_covs (np.ndarray): Filtered covariance matrices, shape (T, n, n).
predicted_states (np.ndarray): One-step-ahead predicted states, shape (T, n).
predicted_covs (np.ndarray): One-step-ahead predicted covariances, shape (T, n, n).
innovations (np.ndarray): Innovation (residual) sequence, shape (T, m).
innovation_covs (np.ndarray): Innovation covariances, shape (T, m, m).
kalman_gains (np.ndarray): Kalman gain matrices, shape (T, n, m).
log_likelihood (float): Total log-likelihood computed from the innovation sequence.

Return type:

Example

>>> # Track a random walk with noisy observations
>>> F = np.array([[1.0]])
>>> H = np.array([[1.0]])
>>> Q = np.array([[0.01]])   # Process noise
>>> R = np.array([[0.1]])    # Observation noise
>>> x0 = np.array([0.0])
>>> P0 = np.array([[1.0]])
>>> obs = np.cumsum(np.random.randn(100, 1) * 0.1) + \
...       np.random.randn(100, 1) * 0.3
>>> result = kalman_filter(obs, F, H, Q, R, x0, P0)
>>> print(result['log_likelihood'])
>>> print(result['filtered_states'][-5:])

See also

kalman_smoother: Backward pass for optimal smoothed estimates. kalman_regression: Time-varying coefficient regression. unscented_kalman: Nonlinear state estimation.

kalman_smoother(observations, F, H, Q, R, x0, P0)[source]¶

Run the Kalman filter followed by the RTS backward smoother.

The Rauch-Tung-Striebel (RTS) smoother uses all observations (past and future) to produce the optimal state estimate at each time step. The smoothed estimates have lower variance than the filtered estimates.

When to use: Offline analysis where you have the complete time series and want the best possible state estimates.

Parameters:

observations (ndarray) – Observations, shape (T, m).
F (ndarray) – State transition matrix (n, n).
H (ndarray) – Observation matrix (m, n).
Q (ndarray) – Process noise covariance (n, n).
R (ndarray) – Observation noise covariance (m, m).
x0 (ndarray) – Initial state (n,).
P0 (ndarray) – Initial covariance (n, n).

Returns:

smoothed_states (np.ndarray): Smoothed state estimates, shape (T, n).
smoothed_covs (np.ndarray): Smoothed covariance matrices, shape (T, n, n).

Return type:

Example

>>> result = kalman_smoother(obs, F, H, Q, R, x0, P0)
>>> # Smoothed states are more accurate than filtered
>>> print(result['smoothed_states'][-5:])
>>> print(result['filtered_states'][-5:])

See also

kalman_filter: Forward pass only. kalman_regression: Time-varying betas using DLM.

kalman_regression(y, X, state_noise=0.0001, obs_noise=None, x0=None, P0_scale=1.0, smooth=True)[source]¶

Estimate time-varying regression coefficients using a Kalman filter.

Models the regression coefficients as a random walk:

\[ \begin{align}\begin{aligned}\beta_t = \beta_{t-1} + w_t, \quad w_t \sim N(0, Q)\\y_t = X_t' \beta_t + v_t, \quad v_t \sim N(0, R)\end{aligned}\end{align} \]

This is a Dynamic Linear Model (DLM) where the state vector is the coefficient vector. Common financial applications:

Time-varying betas: Estimate CAPM beta that changes over time
Dynamic hedge ratios: For pairs trading or hedging
Factor exposure monitoring: Track how factor loadings evolve

Parameters:

y (Series | ndarray) – Dependent variable, shape (T,).
X (DataFrame | ndarray) – Regressors, shape (T, p). Include a column of ones for an intercept.
state_noise (float, default: 0.0001) – Process noise variance for each coefficient. Controls how quickly coefficients can change. Larger values allow faster adaptation.
obs_noise (float | None, default: None) – Observation noise variance. If None, estimated as the variance of OLS residuals.
x0 (ndarray | None, default: None) – Initial coefficient estimate, shape (p,). If None, OLS estimates are used.
P0_scale (float, default: 1.0) – Scale for initial covariance (P0 = P0_scale * I).
smooth (bool, default: True) – If True, apply the RTS smoother for better estimates.

Returns:

coefficients (np.ndarray): Time-varying coefficients, shape (T, p). If smooth=True, these are smoothed.
coefficient_covs (np.ndarray): Covariance of coefficients, shape (T, p, p).
residuals (np.ndarray): Observation residuals, shape (T,).
log_likelihood (float): Log-likelihood.
filtered_coefficients (np.ndarray): Filtered (causal) coefficient estimates, shape (T, p).

Return type:

Example

>>> import numpy as np, pandas as pd
>>> # Time-varying beta estimation
>>> market = np.random.randn(500) * 0.01
>>> # True beta changes from 1.0 to 1.5 halfway
>>> true_beta = np.where(np.arange(500) < 250, 1.0, 1.5)
>>> stock = true_beta * market + np.random.randn(500) * 0.005
>>> X = np.column_stack([np.ones(500), market])
>>> result = kalman_regression(stock, X, state_noise=1e-5)
>>> betas = result['coefficients'][:, 1]  # Time-varying beta
>>> print(f"Beta at t=100: {betas[100]:.2f}")  # ~1.0
>>> print(f"Beta at t=400: {betas[400]:.2f}")  # ~1.5

Notes

The state_noise parameter is critical. Too large and the estimates are noisy; too small and they adapt too slowly. A good starting point is 1e-4 to 1e-6 for daily returns.

See also

kalman_filter: General-purpose Kalman filter. kalman_smoother: RTS backward smoother.

unscented_kalman(observations, state_dim, f_func, h_func, Q, R, x0, P0, alpha=0.001, beta=2.0, kappa=0.0)[source]¶

Unscented Kalman Filter for nonlinear state estimation.

The UKF uses sigma points to propagate the state distribution through nonlinear dynamics and observation functions, avoiding the need for Jacobian computation (unlike the Extended Kalman Filter).

The state-space model:

\[ \begin{align}\begin{aligned}x_t = f(x_{t-1}) + w_t, \quad w_t \sim N(0, Q)\\y_t = h(x_t) + v_t, \quad v_t \sim N(0, R)\end{aligned}\end{align} \]

When to use: When f or h are nonlinear. Examples include stochastic volatility models, nonlinear factor models, or target tracking with angular observations.

Parameters:

observations (ndarray) – Observations, shape (T, m).
state_dim (int) – Dimension of the state vector (n).
f_func (Callable[[ndarray], ndarray]) – State transition function. Takes x of shape (n,) and returns predicted state of shape (n,).
h_func (Callable[[ndarray], ndarray]) – Observation function. Takes x of shape (n,) and returns predicted observation of shape (m,).
Q (ndarray) – Process noise covariance (n, n).
R (ndarray) – Observation noise covariance (m, m).
x0 (ndarray) – Initial state estimate (n,).
P0 (ndarray) – Initial covariance (n, n).
alpha (float, default: 0.001) – Spread of sigma points (typically 1e-4 to 1).
beta (float, default: 2.0) – Prior knowledge of distribution (2 is optimal for Gaussian).
kappa (float, default: 0.0) – Secondary scaling parameter (0 or 3 - n).

Returns:

filtered_states (np.ndarray): Filtered state estimates, shape (T, n).
filtered_covs (np.ndarray): Filtered covariances, shape (T, n, n).
innovations (np.ndarray): Innovation sequence, shape (T, m).
log_likelihood (float): Approximate log-likelihood.

Return type:

Example

>>> # Nonlinear observation model: observe angle to a target
>>> def f(x):
...     return x  # Random walk dynamics
>>> def h(x):
...     return np.array([np.arctan2(x[1], x[0])])  # Angle
>>> result = unscented_kalman(
...     observations, state_dim=2, f_func=f, h_func=h,
...     Q=np.eye(2)*0.01, R=np.array([[0.1]]),
...     x0=np.zeros(2), P0=np.eye(2),
... )

Notes

The UKF with the standard unscented transform uses 2n + 1 sigma points. Computational cost scales as O(n^2) per step (vs O(n^3) for EKF with Jacobian).

See also

kalman_filter: Linear Kalman filter. kalman_smoother: Linear Kalman smoother.

online_changepoint(data, hazard=0.005)[source]¶

Bayesian online change-point detection.

Implements the algorithm from Adams & MacKay (2007). Uses a normal model with unknown mean and known variance (estimated from data).

Parameters:

data (Series) – Time series to monitor.
hazard (float, default: 0.005) – Prior probability that a change-point occurs at each time step (constant hazard rate 1/mean_run_length).

Return type:

Returns:

Series of estimated run lengths (the most probable run length at each time step). A drop to zero indicates a detected change-point.

pelt_changepoint(data, *, penalty='bic', model='l2', min_size=5, jump=1, n_bkps=None)[source]¶

Pruned Exact Linear Time (PELT) optimal changepoint detection.

PELT finds the exact global optimum of the penalised cost function in O(n) expected time (under mild assumptions). It is the recommended method when you want optimal changepoints rather than an approximate or sequential solution.

The algorithm minimises:

\[\sum_{i=0}^{K} C(y_{\tau_i : \tau_{i+1}}) + \beta \, K\]

where C is the segment cost, K is the number of changepoints, and beta is the penalty.

Interpretation guidance:

A changepoint at index t means the statistical properties of the series changed between observation t-1 and t.
segment_stats gives the mean, variance, and length of each segment so you can characterise the regimes between breaks.
confidence is a heuristic based on the cost reduction at each changepoint relative to the global cost. Higher values indicate more confident breaks.

Parameters:

data (Series | ndarray) – Time series to analyse. If a pd.Series, the index is preserved in the output.
penalty (str | float, default: 'bic') –
Penalty method or explicit float value. Supported strings:
- "bic" – Bayesian Information Criterion (default). Works well when the number of changepoints is unknown.
- "aic" – Akaike Information Criterion. Tends to over-segment; prefer BIC unless you want sensitivity.
- "mbic" – Modified BIC (Liu et al., 2016).
If a float is given, it is used directly as the penalty value beta.
model (str, default: 'l2') –
Cost function for each segment. Options:
- "l2" – least-squares cost (Gaussian change in mean).
- "l1" – absolute deviation cost (robust to outliers).
- "rbf" – kernel cost (nonparametric, captures complex distribution changes).
- "linear" – linear regression cost (for trend breaks).
- "normal" – full Gaussian cost (mean and variance).
- "ar" – autoregressive model cost (order-1 AR).
min_size (int, default: 5) – Minimum segment length. Shorter segments are prohibited. Increase to suppress spurious breaks.
jump (int, default: 1) – Subsample factor for candidate changepoints. jump=1 considers every index; larger values trade accuracy for speed on very long series.
n_bkps (int | None, default: None) – If given, ignore penalty and find exactly this many changepoints (uses the Pelt predict with n_bkps).

Returns:

changepoints (list[int]): Detected changepoint indices (0-based). The last element is always the series length.
n_changepoints (int): Number of detected changepoints (excludes the trailing length marker).
segment_stats (list[dict]): Per-segment statistics with keys start, end, mean, var, length.
confidence (np.ndarray): Heuristic confidence score for each changepoint (between 0 and 1).
penalty_value (float): Effective penalty value used.
model (str): Cost model used.
cost (float): Total cost of the segmentation.

Return type:

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(5, 1, 200),
...     rng.normal(0, 2, 200),
... ])
>>> result = pelt_changepoint(pd.Series(data))
>>> print(result["changepoints"])  # ~ [200, 400, 600]
>>> print(result["segment_stats"])

References

Killick, R., Fearnhead, P. & Eckley, I. A. (2012). “Optimal Detection of Changepoints with a Linear Computational Cost.” J. Amer. Statist. Assoc., 107(500).

See also

binary_segmentation: Hierarchical top-down changepoint detection. window_changepoint: Sliding-window change score time series. online_changepoint: Bayesian online detection (sequential).

binary_segmentation(data, *, n_bkps=None, penalty='bic', model='l2', min_size=5, jump=1)[source]¶

Binary segmentation with hierarchical changepoint ordering.

Binary segmentation is a fast, greedy, top-down algorithm that recursively splits the series at the point of maximum cost improvement. Although sub-optimal compared to PELT, it is computationally faster and produces a natural hierarchy of changepoints ordered by significance.

Interpretation guidance:

hierarchical_order[0] is the single most important changepoint in the series. If you can only act on one break, act on this one.
The ordering degrades gracefully: if you want K changepoints, take the first K entries of hierarchical_order.

Parameters:

data (Series | ndarray) – Time series to analyse.
n_bkps (int | None, default: None) – Number of changepoints to detect. If None, the penalty is used to determine the number automatically.
penalty (str | float, default: 'bic') – Penalty method or explicit float. Same options as pelt_changepoint(). Ignored when n_bkps is set.
model (str, default: 'l2') – Cost function ("l1", "l2", "rbf", "linear", "normal").
min_size (int, default: 5) – Minimum segment length.
jump (int, default: 1) – Candidate subsample factor.

Returns:

changepoints (list[int]): Detected changepoints (last element is always series length).
n_changepoints (int): Number of changepoints.
hierarchical_order (list[int]): Changepoints sorted from most significant (largest cost improvement) to least.
segment_stats (list[dict]): Per-segment start, end, mean, var, length.
model (str): Cost model used.

Return type:

Example

>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(5, 2, 100),
...     rng.normal(0, 1, 200),
... ])
>>> result = binary_segmentation(pd.Series(data), n_bkps=2)
>>> # First entry is the dominant changepoint
>>> print(result["hierarchical_order"])

References

Scott, A. J. & Knott, M. (1974). “A Cluster Analysis Method for Grouping Means in the Analysis of Variance.” Biometrics.

See also

pelt_changepoint: Globally optimal changepoint detection. window_changepoint: Sliding-window approach.

window_changepoint(data, *, width=50, model='l2', min_size=5, jump=1, n_bkps=None, penalty='bic')[source]¶

Sliding-window changepoint detection with change-score time series.

Compares the distributions on either side of a sliding window to produce a continuous change score at every time step. Peaks in this score indicate likely changepoints. This approach is well suited for online/streaming applications where you want a continuously updated change intensity signal.

Interpretation guidance:

The change_score series is analogous to a “structural break intensity” indicator. Threshold it or combine it with other signals to trigger regime-switch alerts.
The score is non-negative and unitless; higher values indicate stronger evidence for a distributional shift.

Parameters:

data (Series | ndarray) – Time series to analyse.
width (int, default: 50) – Half-width of the sliding window. Total window is 2 * width. Larger windows are more robust but less responsive.
model (str, default: 'l2') – Cost function ("l2", "l1", "rbf", etc.).
min_size (int, default: 5) – Minimum segment length.
jump (int, default: 1) – Candidate subsample factor.
n_bkps (int | None, default: None) – If given, return exactly this many changepoints.
penalty (str | float, default: 'bic') – Penalty method or float. Used when n_bkps is None.

Returns:

change_score (pd.Series or np.ndarray): Change score at each time step. Same type as input.
changepoints (list[int]): Detected changepoint indices.
n_changepoints (int): Number of changepoints.
width (int): Window half-width used.

Return type:

Example

>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(3, 1, 200),
... ])
>>> result = window_changepoint(pd.Series(data), width=30)
>>> score = result["change_score"]
>>> # Plot score to see the peak near index 200

References

Truong, C., Oudre, L. & Vayer, N. (2020). “Selective review of offline change point detection methods.” Signal Processing.

See also

pelt_changepoint: Offline optimal detection. online_changepoint: Bayesian online detection.

cusum_control_chart(data, *, target=None, std_est=None, k=0.5, h=5.0)[source]¶

CUSUM (Cumulative Sum) control chart for process monitoring.

The CUSUM chart monitors the cumulative sum of deviations from a target value. It is one of the most powerful sequential tests for detecting small, persistent shifts in the mean of a process.

The upper and lower CUSUM statistics are:

\[ \begin{align}\begin{aligned}S^+_t = \max(0,\; S^+_{t-1} + z_t - k)\\S^-_t = \max(0,\; S^-_{t-1} - z_t - k)\end{aligned}\end{align} \]

where \(z_t = (x_t - \mu_0) / \sigma\) is the standardised observation, k is the allowance (slack) parameter, and h is the decision threshold. An alarm is triggered when either \(S^+_t > h\) or \(S^-_t > h\).

Interpretation guidance:

upper_cusum rising indicates the process mean is drifting above target. lower_cusum rising indicates drift below target.
alarm_points are the times at which the CUSUM crossed the control limit. These are your structural break candidates.
arl (average run length) is the expected number of observations between false alarms under the null hypothesis. Higher ARL means fewer false positives. With k = 0.5 and h = 5, the in-control ARL is approximately 465 (for a standard normal process).

Parameters:

data (Series | ndarray) – Time series to monitor.
target (float | None, default: None) – Target (in-control) mean. If None, the overall sample mean is used.
std_est (float | None, default: None) – Standard deviation estimate. If None, the sample standard deviation is used.
k (float, default: 0.5) – Allowance (slack) parameter, in units of standard deviation. Controls the size of the shift you want to detect. Default 0.5 detects a 1-sigma shift optimally.
h (float, default: 5.0) – Decision interval (control limit), in units of standard deviation. Larger h reduces false alarms but increases detection delay.

Returns:

upper_cusum (np.ndarray): Upper CUSUM statistic at each time step.
lower_cusum (np.ndarray): Lower CUSUM statistic.
control_limit (float): The h threshold used.
alarm_points (list[int]): Indices where an alarm was triggered (CUSUM exceeded h).
alarm_direction (list[str]): Direction of each alarm ("upper" or "lower").
arl (float): Estimated average run length (between alarms) based on the observed alarm sequence. inf if no alarms were triggered.
target (float): Target mean used.
std (float): Standard deviation estimate used.

Return type:

Example

>>> rng = np.random.default_rng(42)
>>> # In-control for 200 obs, then shift mean by +1 sigma
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(1, 1, 200),
... ])
>>> result = cusum_control_chart(pd.Series(data), target=0, std_est=1)
>>> print(result["alarm_points"])  # first alarm near index 200
>>> print(f"ARL: {result['arl']:.0f}")

References

Page, E. S. (1954). “Continuous Inspection Schemes.” Biometrika, 41(1/2).

Montgomery, D. C. (2009). Introduction to Statistical Quality Control. 6th ed. Wiley.

See also

online_changepoint: Bayesian approach to sequential detection. window_changepoint: Sliding-window approach.

label_regimes(states, returns)[source]¶

Assign descriptive labels to numeric regime states.

States are sorted by mean return: the state with the highest mean return is labeled "bull", the lowest "bear", and any intermediate states "neutral_1", "neutral_2", etc.

Parameters:

states (Series) – Integer regime state series.
returns (Series) – Corresponding return series (same index).

Return type:

Returns:

Series of string regime labels.

volatility_regime_labels(returns, *, window=21, n_levels=3, quantiles=None)[source]¶

Label regimes based on realised volatility quantiles.

A simple, model-free approach that classifies each period by where its rolling volatility falls within the historical distribution. No fitting, no hidden states – just raw vol percentiles.

Interpretation guidance:

"low_vol" periods typically correspond to trending or complacent markets. Strategy-wise, favour momentum and carry.
"high_vol" periods correspond to stressed or mean-reverting markets. Favour defensive positioning or mean-reversion.
"medium_vol" is the transition zone.

Parameters:

returns (Series | ndarray) – Return series.
window (int, default: 21) – Rolling window for realised volatility estimation. Default 21 (roughly one trading month).
n_levels (int, default: 3) – Number of volatility levels. Default 3 produces low_vol / medium_vol / high_vol. Use 2 for a binary split or 4+ for finer granularity.
quantiles (list[float] | None, default: None) – Explicit quantile boundaries. If provided, overrides n_levels. Must have n_levels - 1 elements, each in (0, 1).

Return type:

Returns:

pd.Series of string labels (e.g., "low_vol", "medium_vol", "high_vol"). NaN-filled for the warm-up period where rolling volatility is unavailable.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(rng.normal(0, 0.01, 500))
>>> labels = volatility_regime_labels(returns, n_levels=3)
>>> print(labels.value_counts())

See also

trend_regime_labels: Label by trend direction. composite_regime_labels: Combine vol + trend labels.

trend_regime_labels(returns, *, fast_window=10, slow_window=50, hysteresis=0.0005)[source]¶

Label regimes based on moving average slope with hysteresis.

Uses a dual moving-average crossover system with a hysteresis band to avoid whipsaw signals. The result is a clean, three-state classification: uptrend, downtrend, or sideways.

Interpretation guidance:

"uptrend": Fast MA is above slow MA by more than the hysteresis threshold. Bullish bias.
"downtrend": Fast MA is below slow MA by more than the hysteresis threshold. Bearish bias.
"sideways": The two MAs are within the hysteresis band. No directional conviction – favour range-bound strategies.

Parameters:

returns (Series | ndarray) – Return series.
fast_window (int, default: 10) – Fast moving average window (periods).
slow_window (int, default: 50) – Slow moving average window (periods).
hysteresis (float, default: 0.0005) – Minimum difference between fast and slow MA (in return units) required to declare a trend. Larger values suppress whipsaws but delay signals.

Return type:

Returns:

pd.Series of string labels ("uptrend", "downtrend", "sideways"). NaN-filled during warm-up.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(rng.normal(0.001, 0.01, 500))
>>> labels = trend_regime_labels(returns)
>>> print(labels.value_counts())

See also

volatility_regime_labels: Label by vol level. composite_regime_labels: Combine vol + trend labels.

composite_regime_labels(returns, *, vol_window=21, fast_window=10, slow_window=50, hysteresis=0.0005, n_vol_levels=2)[source]¶

Combine volatility and trend regimes into composite states.

Creates 4-6 composite labels by crossing trend direction (uptrend / downtrend / sideways) with volatility level (low / high or low / medium / high). Common composite states:

bull_calm: Uptrend + low vol. The best environment for passive equity holding.
bull_volatile: Uptrend + high vol. Often late-cycle or recovery rallies.
bear_calm: Downtrend + low vol. Grinding bear markets.
bear_volatile: Downtrend + high vol. Crisis periods (2008, March 2020).
sideways_calm: Range-bound, quiet.
sideways_volatile: Choppy, difficult to trade.

Interpretation guidance:

The composite label captures both direction and turbulence, which together determine the optimal strategy. For instance, momentum strategies work in bull_calm but fail in bear_volatile.

Parameters:

returns (Series | ndarray) – Return series.
vol_window (int, default: 21) – Window for rolling volatility.
fast_window (int, default: 10) – Fast MA window for trend.
slow_window (int, default: 50) – Slow MA window for trend.
hysteresis (float, default: 0.0005) – Trend hysteresis threshold.
n_vol_levels (int, default: 2) – 2 or 3 volatility levels.

Return type:

Returns:

pd.Series of string composite labels. NaN-filled during warm-up.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(rng.normal(0.001, 0.01, 500))
>>> labels = composite_regime_labels(returns)
>>> print(labels.value_counts())

See also

volatility_regime_labels: Volatility-only labeling. trend_regime_labels: Trend-only labeling. regime_duration_analysis: Analyse how long each composite

state typically lasts.

regime_duration_analysis(states)[source]¶

Analyse how long each regime typically lasts.

Computes the survival function, hazard rate, and expected remaining duration for each regime. This helps answer questions like “we’ve been in a bull regime for 60 days – how much longer can we expect it to last?”

Interpretation guidance:

survival_curve[k]: Probability that a regime-k spell lasts at least d periods. A slowly-decaying curve means the regime tends to persist.
hazard_rate[k]: Instantaneous probability of exiting regime k after having been in it for d periods. If the hazard rate is approximately constant, regime duration is memoryless (geometric distribution, consistent with Markov). An increasing hazard rate means longer spells are more likely to end soon.
expected_remaining[k]: Given that we are currently in regime k and have been for d periods, how many more periods should we expect? Computed from the empirical survival function.

Parameters:

states (Series | ndarray) – Integer regime labels, shape (T,).

Returns:

durations (dict[int, list[int]]): List of spell durations for each regime.
survival_curve (dict[int, pd.Series]): Kaplan-Meier style survival curve for each regime, indexed by duration.
hazard_rate (dict[int, pd.Series]): Empirical hazard rate for each regime, indexed by duration.
expected_remaining (dict[int, pd.Series]): Expected remaining duration conditional on having survived d periods, indexed by duration.
summary (pd.DataFrame): Per-regime summary with mean_duration, median_duration, max_duration, n_spells.

Return type:

Example

>>> states = np.array([0]*50 + [1]*30 + [0]*80 + [1]*40)
>>> result = regime_duration_analysis(states)
>>> print(result["summary"])
>>> # Survival curve for regime 0
>>> print(result["survival_curve"][0])

See also

regime_stability_score: Composite stability metric. composite_regime_labels: Generate regime labels to analyse.

regime_stability_score(states, transition_matrix=None)[source]¶

Measure how stable the detected regimes are.

Stability quantifies the temporal persistence of regimes. Stable regimes last many periods and transitions are infrequent. A stability score close to 1 indicates regimes that are highly persistent; a score near 0 indicates rapid, noisy switching that is likely untradable.

How it works:

Three sub-metrics are combined (equal-weighted average):

Average duration (normalised): longer average stays = more stable. Normalised as 1 - 1/avg_duration so it maps to [0, 1].
Transition frequency (inverse): fewer transitions per unit time = more stable. Computed as 1 - n_transitions / (T-1).
Probability entropy (inverse): when the transition matrix rows have low entropy (i.e., dominated by the diagonal), the chain is persistent. Normalised by log(K) where K is the number of regimes.

Interpretation guidance:

> 0.8: Highly stable regimes, suitable for discrete portfolio switches.
0.5 – 0.8: Moderate stability, consider probability-weighted blending rather than hard switches.
< 0.5: Unstable, likely overfitting or too many regimes.

Parameters:

states (ndarray | Series) – Integer regime labels, shape (T,).
transition_matrix (ndarray | None, default: None) – Optional (K, K) transition matrix. If None, an empirical matrix is computed from states.

Returns:

stability_score (float): Composite stability score in [0, 1].
avg_duration (float): Average regime duration in periods.
transition_frequency (float): Fraction of time steps with a regime change.
probability_entropy (float): Mean row entropy of the transition matrix (nats).
n_transitions (int): Total number of regime transitions.
per_regime_duration (dict[int, float]): Average duration per regime.

Return type:

Example

>>> states = np.array([0]*100 + [1]*100 + [0]*100)
>>> result = regime_stability_score(states)
>>> print(f"Stability: {result['stability_score']:.2f}")
0.99

See also

regime_separation_score: Distribution distance between regimes. compare_regime_methods: Multi-method comparison.

regime_separation_score(returns, states)[source]¶

Measure how well-separated the regime distributions are.

Separation quantifies whether the regimes have meaningfully distinct statistical properties. High separation means the regimes carry actionable information; low separation means the classification is arbitrary.

How it works:

For each pair of regimes (i, j):

Bhattacharyya distance: measures divergence between two Gaussian distributions. Larger = more separated.
Overlap coefficient: fraction of probability mass that overlaps when the two regime distributions are superimposed. Smaller = more separated.

The composite score is the average pairwise Bhattacharyya distance, normalised to [0, 1] via a sigmoid transform.

Interpretation guidance:

> 0.7: Well-separated regimes with distinct return distributions.
0.3 – 0.7: Moderate separation. Regimes differ but there is meaningful overlap.
< 0.3: Poor separation. Consider reducing the number of regimes.

Parameters:

returns (ndarray | Series) – Return series, shape (T,).
states (ndarray | Series) – Integer regime labels, shape (T,).

Returns:

separation_score (float): Composite separation in [0, 1].
pairwise_bhattacharyya (dict[tuple, float]): Bhattacharyya distance for each regime pair (i, j).
pairwise_overlap (dict[tuple, float]): Overlap coefficient for each pair.
per_regime_stats (dict[int, dict]): Mean and std per regime.

Return type:

Example

>>> rng = np.random.default_rng(0)
>>> returns = np.concatenate([
...     rng.normal(0.01, 0.01, 200),
...     rng.normal(-0.02, 0.03, 200),
... ])
>>> states = np.array([0]*200 + [1]*200)
>>> result = regime_separation_score(returns, states)
>>> print(f"Separation: {result['separation_score']:.2f}")

See also

regime_stability_score: Temporal persistence of regimes. compare_regime_methods: Multi-method comparison.

regime_predictability(states, transition_matrix=None, test_fraction=0.3)[source]¶

Assess whether regime transitions can be predicted.

A regime model is useful for trading only if you can anticipate transitions before they happen (or at least classify the current regime in real time). This function evaluates predictability using transition-probability-based forecasting.

How it works:

Split the state sequence into train / test at (1 - test_fraction) * T.
Estimate the transition matrix from the training set.
For each test observation, predict the next state as the most probable successor under the estimated transition matrix.
Compute out-of-sample accuracy and compare to a naive baseline (always predict the most common state).

Interpretation guidance:

accuracy > baseline + 0.1: The transition structure adds predictive value. Use regime probabilities to adjust portfolio weights.
accuracy ~ baseline: No predictive value from transitions. Consider whether the regime labeling is too noisy or the number of states is wrong.
economic_value: A rough estimate of the annualised Sharpe improvement from perfect regime prediction vs buy-and-hold. This uses per-regime Sharpe ratios and is indicative only.

Parameters:

states (ndarray | Series) – Integer regime labels, shape (T,).
transition_matrix (ndarray | None, default: None) – If provided, used instead of the training-set estimate.
test_fraction (float, default: 0.3) – Fraction of data reserved for testing.

Returns:

predictability_score (float): Composite predictability in [0, 1]. Combines accuracy lift and transition regularity.
accuracy (float): Out-of-sample one-step prediction accuracy.
baseline_accuracy (float): Accuracy of always predicting the most common state.
accuracy_lift (float): accuracy - baseline_accuracy.
transition_matrix_train (np.ndarray): Transition matrix estimated from training data.

Return type:

Example

>>> # Highly predictable two-state chain
>>> rng = np.random.default_rng(0)
>>> states = np.zeros(500, dtype=int)
>>> for t in range(1, 500):
...     if states[t-1] == 0:
...         states[t] = 0 if rng.random() < 0.95 else 1
...     else:
...         states[t] = 1 if rng.random() < 0.90 else 0
>>> result = regime_predictability(states)
>>> print(f"Accuracy: {result['accuracy']:.2f}")

See also

regime_stability_score: Regime persistence. regime_separation_score: Distribution distinctness. compare_regime_methods: Multi-method comparison.

compare_regime_methods(returns, methods=None, method_kwargs=None)[source]¶

Compare multiple regime detection methods on the same data.

Runs each method on the provided return series and evaluates stability, separation, and predictability. Returns a DataFrame for easy comparison.

Interpretation guidance:

Use this to pick the best regime model for your asset or strategy. The “best” model depends on your use-case: high stability for discrete allocation switches, high separation for regime-conditional signals, high predictability for transition timing.

Parameters:

returns (Series | ndarray) – Return series, shape (T,).
methods (dict[str, Callable[..., Any]] | None, default: None) –
Dictionary mapping method names to callables. Each callable must accept returns and return a tuple (states, transition_matrix) where states is an integer array and transition_matrix is (K, K) or None.

If None, a default set is used (requires detect_regimes from wraquant.regimes.base).
method_kwargs (dict[str, dict[str, Any]] | None, default: None) – Per-method keyword arguments. Keys should match methods keys.

Returns:

stability, separation, predictability
avg_duration, n_transitions, n_regimes
composite (equal-weighted average of the three scores)

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = np.concatenate([
...     rng.normal(0.01, 0.01, 200),
...     rng.normal(-0.02, 0.03, 200),
... ])
>>> df = compare_regime_methods(pd.Series(returns))
>>> print(df.sort_values("composite", ascending=False))

See also

regime_stability_score: Stability metric details. regime_separation_score: Separation metric details. regime_predictability: Predictability metric details.

pomegranate_hmm(data, n_states=2)[source]¶

Fit a Hidden Markov Model using pomegranate.

Parameters:

data (Series | ndarray) – Univariate observations (e.g. returns).
n_states (int, default: 2) – Number of hidden states.

Returns:

Dictionary containing:

states – predicted hidden state sequence (1-D int array).
means – estimated mean of each state’s emission distribution.
model – the fitted pomegranate DenseHMM object.
n_states – number of hidden states.

Return type:

filterpy_kalman(observations, F, H, Q, R, x0=None, P0=None)[source]¶

Run a Kalman filter using filterpy.

Parameters:

observations (ndarray | Series) – Observed measurements. Shape (T,) for univariate or (T, dim_z) for multivariate observations.
F (ndarray) – State transition matrix. Shape (dim_x, dim_x).
H (ndarray) – Observation (measurement) matrix. Shape (dim_z, dim_x).
Q (ndarray) – Process noise covariance. Shape (dim_x, dim_x).
R (ndarray) – Measurement noise covariance. Shape (dim_z, dim_z).
x0 (ndarray | None, default: None) – Initial state estimate. Shape (dim_x,) or (dim_x, 1). Defaults to zeros.
P0 (ndarray | None, default: None) – Initial covariance estimate. Shape (dim_x, dim_x). Defaults to identity.

Returns:

Dictionary containing:

filtered_states – filtered state estimates, shape (T, dim_x).
filtered_covariances – filtered covariance matrices, shape (T, dim_x, dim_x).
log_likelihood – total log-likelihood.
residuals – measurement residuals, shape (T, dim_z).

Return type:

river_drift_detector(stream, method='adwin')[source]¶

Detect concept drift in a data stream using river.

Processes each observation sequentially and records indices where a drift is detected.

Parameters:

stream (ndarray | Series | list[float]) – Sequential stream of numeric observations.
method (str, default: 'adwin') –
Drift detection method:
- 'adwin' – Adaptive Windowing (ADWIN)
- 'ddm' – Drift Detection Method
- 'eddm' – Early Drift Detection Method
- 'page_hinkley' – Page-Hinkley test

Returns:

Dictionary containing:

drift_indices – list of indices where drift was detected.
n_drifts – total number of drifts detected.
method – drift detection method used.

Return type:

dynamax_lgssm(observations, state_dim=2, emission_dim=None, n_iters=100)[source]¶

Fit a Linear Gaussian State Space Model using dynamax (JAX-based).

Parameters:

observations (ndarray | Series) – Observations of shape (T, obs_dim) or (T,) for univariate.
state_dim (int, default: 2) – Dimensionality of the latent state.
emission_dim (int | None, default: None) – Observation dimension. Inferred from observations if None.
n_iters (int, default: 100) – Number of EM iterations.

Returns:

Dictionary containing:

filtered_means – filtered state means, shape (T, state_dim).
filtered_covs – filtered state covariances.
smoothed_means – smoothed state means.
params – learned model parameters.
log_likelihoods – EM log-likelihood trace.

Return type:

pykalman_filter(observations, transition_matrices=None, observation_matrices=None, n_em_iter=10)[source]¶

Kalman filter with EM parameter learning via pykalman.

Unlike the pure-numpy Kalman in regimes/kalman.py, pykalman’s KalmanFilter can learn the transition/observation matrices and noise covariances from data via Expectation-Maximization. Use this when you don’t know the system parameters.

Parameters:

observations (ndarray | Series) – Observations of shape (T, obs_dim) or (T,) for univariate data.
transition_matrices (ndarray | None, default: None) – Initial guess for the state transition matrix. When None, pykalman learns it from data via EM.
observation_matrices (ndarray | None, default: None) – Initial guess for the observation matrix. When None, pykalman learns it from data via EM.
n_em_iter (int, default: 10) – Number of EM iterations for parameter learning.

Returns:

Dictionary containing:

filtered_means – filtered state estimates, shape (T, state_dim).
filtered_covs – filtered state covariance matrices, shape (T, state_dim, state_dim).
smoothed_means – smoothed (RTS) state estimates, shape (T, state_dim).
smoothed_covs – smoothed state covariance matrices, shape (T, state_dim, state_dim).
learned_params – dict with learned transition, observation, transition_cov, and observation_cov matrices.
log_likelihood – total log-likelihood of the data under the learned model.

Return type:

Example

>>> import numpy as np
>>> from wraquant.regimes.integrations import pykalman_filter
>>> obs = np.cumsum(np.random.default_rng(0).normal(0, 1, 100))
>>> result = pykalman_filter(obs, n_em_iter=5)
>>> result["smoothed_means"].shape
(100, 1)

Notes

Reference: Shumway & Stoffer (1982). “An Approach to Time Series Smoothing and Forecasting Using the EM Algorithm.” Journal of Time Series Analysis, 3(4), 253-264.

See also

filterpy_kalman: When system parameters are known a priori.
dynamax_lgssm: JAX-based alternative with GPU support.

Base Detection¶

Regime detection abstractions and unified interfaces.

Provides a standardized framework for regime detection so that all methods (HMM, GMM, Markov-switching, changepoint) return consistent results and can be swapped without changing downstream code.

class RegimeResult[source]¶

Bases: object

Standardized container for regime detection results.

All regime detection functions return this object, ensuring consistent access to states, probabilities, and statistics regardless of the detection method used.

states¶: Integer regime labels (T,). State 0 is always the lowest-volatility regime.

probabilities¶: Posterior regime probabilities (T, K). Row t gives P(regime=k | data) at time t.

transition_matrix¶: (K, K) matrix where entry (i,j) is P(next_regime=j | current_regime=i). Rows sum to 1.

n_regimes¶: Number of detected regimes.

means¶: Per-regime mean returns. Shape (K,) for univariate, (K, n_assets) for multivariate.

covariances¶: Per-regime variance/covariance. Shape (K,) for univariate, (K, n_assets, n_assets) for multivariate.

statistics¶: pd.DataFrame with per-regime summary stats (mean, std, sharpe, max_drawdown, duration, pct_time).

method¶: String identifying the detection method used.

model¶: The underlying fitted model object (for advanced use).

metadata¶: Additional method-specific results.

Parameters:

states (ndarray)
probabilities (ndarray)
transition_matrix (ndarray)
n_regimes (int)
means (ndarray)
covariances (ndarray)
statistics (DataFrame)
method (str)
model (Any, default: None)
metadata (dict, default: <factory>)

states: ndarray¶

probabilities: ndarray¶

transition_matrix: ndarray¶

n_regimes: int¶

means: ndarray¶

covariances: ndarray¶

statistics: DataFrame¶

method: str¶

model: Any = None¶

metadata: dict¶

property current_regime: int¶: Most likely current regime (last observation).

property current_probabilities: ndarray¶: Regime probabilities for the most recent observation.

property steady_state: ndarray¶: Long-run (ergodic) regime distribution.

property expected_durations: ndarray¶

1 / (1 - p_ii).

Type:: Expected duration in each regime

filter_signals(signals)[source]¶

Filter trading signals based on regime probabilities.

Delegates to wraquant.backtest.position.regime_signal_filter to zero out signals in unfavourable regimes.

Parameters:: signals (Series | ndarray) – Raw trading signal series (same length as states).
Return type:: Series | ndarray
Returns:: Regime-filtered signals.

plot()[source]¶

Plot regime detection results using viz dashboard.

Return type:: Any
Returns:: Plotly figure object from viz.dashboard.regime_dashboard.

summary()[source]¶

Human-readable summary of regime detection results.

Return type:: str
Returns:: Multi-line string with regime statistics.

__init__(states, probabilities, transition_matrix, n_regimes, means, covariances, statistics, method, model=None, metadata=<factory>)¶

Parameters:

states (ndarray)
probabilities (ndarray)
transition_matrix (ndarray)
n_regimes (int)
means (ndarray)
covariances (ndarray)
statistics (DataFrame)
method (str)
model (Any, default: None)
metadata (dict, default: <factory>)

Return type:

None

class RegimeDetector[source]¶

Bases: ABC

Base class for regime detection methods.

Subclass this to create custom regime detectors that integrate seamlessly with wraquant’s regime framework.

abstractmethod fit(returns)[source]¶

Fit the detector to return data.

Parameters:: returns (Series | DataFrame | ndarray) – Return series or multi-asset return DataFrame.
Return type:: RegimeDetector
Returns:: self, for method chaining.

abstractmethod predict(returns)[source]¶

Assign regime labels to each observation.

Parameters:: returns (Series | DataFrame | ndarray) – Return series.
Return type:: ndarray
Returns:: Integer state labels, shape (T,).

abstractmethod predict_proba(returns)[source]¶

Return posterior regime probabilities.

Parameters:: returns (Series | DataFrame | ndarray) – Return series.
Return type:: ndarray
Returns:: Probabilities array, shape (T, K).

abstractmethod to_result()[source]¶

Package fitted results into a RegimeResult.

Return type:: RegimeResult
Returns:: RegimeResult with all standard fields populated.

detect_regimes(returns, method='hmm', n_regimes=2, **kwargs)[source]¶

Detect market regimes using the specified method.

This is the primary entry point for regime detection. It dispatches to the appropriate specialized function and returns a standardized RegimeResult regardless of method.

Parameters:

returns (Series | DataFrame | ndarray) – Return series or multi-asset return DataFrame.
method (str, default: 'hmm') –
Detection method. Options: - "hmm" – Gaussian Hidden Markov Model (default).

Best for: capturing temporal regime persistence.
- "gmm" – Gaussian Mixture Model. Best for: quick classification without temporal structure.
- "ms_regression" – Markov-switching regression. Best for: regime-dependent mean and variance modeling.
- "changepoint" – Bayesian online changepoint detection. Best for: detecting abrupt structural breaks.
- "kmeans" – K-means on rolling statistics. Best for: simple, fast, interpretable regime labels.
n_regimes (int, default: 2) – Number of regimes to detect.
**kwargs (Any) – Additional arguments passed to the underlying method.

Return type:

RegimeResult

Returns:

RegimeResult with standardized fields.

Example

>>> result = detect_regimes(returns, method="hmm", n_regimes=2)
>>> print(f"Current regime: {result.current_regime}")
>>> print(f"Regime probabilities: {result.current_probabilities}")
>>> print(result.statistics)

regime_report(returns, result, risk_free_rate=0.0, annualization=252)[source]¶

Generate a comprehensive regime analysis report.

Parameters:

returns (Series | DataFrame) – Return series used for detection.
result (RegimeResult) – RegimeResult from detect_regimes or any detection function.
risk_free_rate (float, default: 0.0) – Annual risk-free rate for Sharpe calculations.
annualization (int, default: 252) – Trading days per year.

Returns:

summary – pd.DataFrame with per-regime statistics
current_regime – Current regime info (label, probability, duration so far)
transition_analysis – Expected durations, steady state, visit counts
regime_history – pd.DataFrame with date, regime, probability columns
allocation_suggestion – Suggested portfolio adjustment based on regime
risk_assessment – Per-regime VaR, CVaR, max drawdown

Return type:

Hidden Markov Models¶

Hidden Markov Model and mixture-model regime detection.

This module provides production-quality regime detection for financial time series using Gaussian Hidden Markov Models (HMMs), Markov-switching models, and Gaussian mixture models (GMMs).

When to use each model:

Gaussian HMM (fit_gaussian_hmm): The workhorse. Models returns as emissions from a latent Markov chain. Captures temporal persistence of regimes (markets stay in bull/bear for extended periods).
Markov-Switching Regression (fit_ms_regression): When you want regime-dependent regression coefficients or intercepts. Useful for modeling time-varying factor exposures.
Gaussian Mixture (gaussian_mixture_regimes): Quick classification without temporal structure. Good for labeling return distributions when you do not care about transition dynamics.

References

Hamilton, J. D. (1989). “A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle.” Econometrica, 57(2).

Ang, A. & Bekaert, G. (2002). “Regime Switches in Interest Rates.” Journal of Business & Economic Statistics, 20(2).

fit_gaussian_hmm(returns, n_states=2, covariance_type='full', n_init=10, n_iter=200, tol=0.0001, random_state=42)[source]¶

Fit a Gaussian Hidden Markov Model to financial return data.

A Gaussian HMM assumes that observed returns are generated by a latent (hidden) Markov chain with K states. Each state k has its own Gaussian emission distribution N(mu_k, sigma_k). The model learns:

Transition matrix A[i,j] = P(s_t = j | s_{t-1} = i)
Emission means mu_k for each state
Emission covariances Sigma_k for each state
Initial state distribution pi_k = P(s_0 = k)

Typical interpretation for 2-state model on equity returns:

Low-volatility state (“risk-on” / bull): Positive mean return, low standard deviation, high Sharpe ratio. Markets are calm.
High-volatility state (“risk-off” / crisis): Near-zero or negative mean return, high standard deviation. Markets are stressed.

States are automatically re-ordered so that state 0 has the lowest volatility and state K-1 has the highest volatility. This ensures consistent labeling across runs.

Parameters:

returns (Series | ndarray) – Return series (simple or log returns). If a pandas Series, NaN values are dropped and the index is preserved.
n_states (int, default: 2) – Number of hidden states. Common choices are 2 (bull/bear) or 3 (bull/neutral/bear). Values above 5 risk overfitting.
covariance_type (str, default: 'full') – Type of covariance matrix. One of 'full', 'diag', 'spherical', 'tied'.
n_init (int, default: 10) – Number of random restarts for EM. The model with the highest log-likelihood is selected. More restarts reduce the chance of a poor local optimum.
n_iter (int, default: 200) – Maximum number of EM iterations per restart.
tol (float, default: 0.0001) – Convergence threshold for EM.
random_state (int, default: 42) – Base random seed for reproducibility.

Returns:

states (np.ndarray): Most likely state sequence (Viterbi), shape (T,), integers in [0, K-1].
state_probs (np.ndarray): Smoothed state probabilities, shape (T, K). Each row sums to 1.
transition_matrix (np.ndarray): Transition probability matrix, shape (K, K). Row i sums to 1.
means (np.ndarray): Emission means, shape (K,).
covariances (np.ndarray): Emission covariances, shape (K,) (variance per state for univariate data).
startprob (np.ndarray): Initial state distribution, shape (K,).
log_likelihood (float): Total log-likelihood of the data.
aic (float): Akaike Information Criterion.
bic (float): Bayesian Information Criterion.
n_states (int): Number of states.
steady_state (np.ndarray): Ergodic (stationary) distribution of the Markov chain, shape (K,).
avg_duration (np.ndarray): Expected duration in each state, 1 / (1 - A[i,i]), shape (K,).
model (GaussianHMM): The fitted hmmlearn model object.
index (pd.Index | None): Original index if input was a pandas Series.

Return type:

Example

>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(42)
>>> # Simulate 2-regime returns
>>> low_vol = rng.normal(0.001, 0.008, 300)
>>> high_vol = rng.normal(-0.002, 0.025, 200)
>>> returns = pd.Series(np.concatenate([low_vol, high_vol]))
>>> result = fit_gaussian_hmm(returns, n_states=2)
>>> print(result['means'])        # Per-state means
>>> print(result['transition_matrix'])  # Transition probabilities
>>> print(result['steady_state'])  # Long-run regime probabilities

Notes

The model uses hmmlearn’s GaussianHMM. For very long series (>10,000 observations), consider downsampling or using the variational Bayes variant.

See also

fit_hmm(returns, n_states=2)[source]¶

Fit a Gaussian HMM to return data (legacy interface).

This is a thin compatibility wrapper around fit_gaussian_hmm. New code should use fit_gaussian_hmm directly for richer output.

Parameters:

returns (Series) – Simple return series.
n_states (int, default: 2) – Number of hidden states.

Return type:

Any

Returns:

Fitted hmmlearn.hmm.GaussianHMM model.

See also

fit_gaussian_hmm: Full-featured HMM fitting with rich output.

predict_regime(model, returns)[source]¶

Predict regime states and probabilities for new data.

Given a fitted HMM (from fit_gaussian_hmm or hmmlearn directly), decode the most likely hidden state sequence and compute posterior state probabilities for new observations.

Parameters:

model (Any) – A fitted HMM model with predict and predict_proba methods (e.g., from hmmlearn), or the dict returned by fit_gaussian_hmm (in which case the 'model' key is used).
returns (Series | ndarray) – Return series to classify. NaN values are dropped.

Returns:

states (np.ndarray): Most likely state at each time step.
state_probs (np.ndarray): Posterior probabilities, shape (T, K).
index (pd.Index | None): Preserved index if input was a pandas Series.

Return type:

Example

>>> result = fit_gaussian_hmm(train_returns, n_states=2)
>>> pred = predict_regime(result, test_returns)
>>> print(pred['states'])
>>> print(pred['state_probs'][:5])

fit_ms_regression(endog, k_regimes=2, trend='c', switching_variance=True, exog=None, n_iter=200)[source]¶

Fit a Markov-switching regression model using statsmodels.

The model:

\[y_t = \mu_{s_t} + x_t' \beta_{s_t} + \sigma_{s_t} \epsilon_t\]

where \(s_t \in \{0, 1, \ldots, K-1\}\) follows a first-order Markov chain with transition matrix \(P\).

Parameters:

endog (Series | ndarray) – Endogenous (dependent) variable. Typically a return series.
k_regimes (int, default: 2) – Number of regimes (hidden states).
trend (str, default: 'c') – Trend specification for the mean equation. 'c' (constant), 't' (linear trend), 'ct' (both), 'n' (no trend).
switching_variance (bool, default: True) – If True, the variance switches across regimes. Almost always desirable for financial data.
exog (DataFrame | ndarray | None, default: None) – Optional exogenous regressors, shape (T, p). The coefficients on these regressors are not switching by default (the intercept and variance switch).
n_iter (int, default: 200) – Maximum number of EM iterations.

Returns:

smoothed_probs (np.ndarray): Smoothed regime probabilities, shape (T, K). smoothed_probs[:, k] is the probability of being in regime k at each time step given all data.
filtered_probs (np.ndarray): Filtered regime probabilities (using only past data), shape (T, K).
states (np.ndarray): Most likely state sequence (argmax of smoothed probabilities), shape (T,).
transition_matrix (np.ndarray): Transition matrix, shape (K, K).
regime_params (dict): Regime-specific parameters. Keys are 'mean_0', 'mean_1', …, 'sigma2_0', etc.
log_likelihood (float): Log-likelihood of the fitted model.
aic (float): Akaike Information Criterion.
bic (float): Bayesian Information Criterion.
summary (str): Text summary from statsmodels.
model_result (MarkovRegressionResults): The full statsmodels result object.

Return type:

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(np.concatenate([
...     rng.normal(0.001, 0.01, 200),
...     rng.normal(-0.002, 0.03, 200),
... ]))
>>> result = fit_ms_regression(returns, k_regimes=2)
>>> print(result['transition_matrix'])
>>> print(result['regime_params'])

Notes

Requires statsmodels >= 0.13. The model uses Hamilton’s (1989) filter for inference. Convergence can be sensitive to starting values.

See also

fit_gaussian_hmm: HMM using hmmlearn (more flexible emission: distributions).

regime_statistics: Compute financial statistics per regime.

regime_statistics(returns, states)[source]¶

Compute descriptive statistics for each detected regime.

After fitting an HMM or GMM, this function summarizes the financial characteristics of each regime. Typical interpretation for a 2-state model on equity returns:

Regime 0 (low vol): “risk-on” or “bull market” – positive mean return, low volatility, high Sharpe ratio.
Regime 1 (high vol): “risk-off” or “crisis” – negative or zero mean, high volatility, negative Sharpe.

Parameters:

returns (Series | ndarray) – Return series aligned with states. Can be a pandas Series or a numpy array.
states (ndarray) – Integer regime labels (0, 1, …, K-1), same length as returns.

Returns:

mean, std, sharpe, sortino_ratio, min, max, skewness, kurtosis, max_drawdown, VaR_95, CVaR_95, n_observations, pct_time, avg_duration.

sharpe: Annualized Sharpe ratio assuming 252 trading days.
sortino_ratio: Annualized Sortino ratio (downside deviation only, 252 trading days).
max_drawdown: Maximum peak-to-trough drawdown within the regime’s return segments.
VaR_95: 5th percentile of returns (historical Value at Risk at the 95 % confidence level).
CVaR_95: Mean of returns below the VaR_95 threshold (Conditional Value at Risk / Expected Shortfall).
avg_duration: Mean consecutive time steps spent in the regime.

Return type:

Example

>>> result = fit_gaussian_hmm(returns, n_states=2)
>>> stats = regime_statistics(returns, result['states'])
>>> print(stats)
            mean       std    sharpe  ...  avg_duration
regime
0       0.000523  0.008123  1.0223   ...          15.3
1      -0.001842  0.024561 -1.1917   ...           7.8

See also

regime_transition_analysis(states, transition_matrix=None)[source]¶

Analyze regime transition dynamics.

Computes empirical and model-based statistics about how regimes switch. Useful for understanding regime persistence and expected timing of regime changes.

Parameters:

states (ndarray) – Integer regime labels, shape (T,).
transition_matrix (ndarray | None, default: None) – Model-estimated transition matrix, shape (K, K). If None, the empirical transition matrix is computed from states.

Returns:

empirical_transition_matrix (np.ndarray): Transition counts normalised to probabilities, shape (K, K).
transition_matrix (np.ndarray): Model or empirical transition matrix used for analysis.
avg_duration (np.ndarray): Expected duration in each regime (time steps), 1 / (1 - P[i,i]).
steady_state (np.ndarray): Ergodic distribution from the transition matrix.
regime_counts (dict): Number of distinct “visits” to each regime.
regime_durations (dict): List of consecutive durations for each regime.
expected_return_time (np.ndarray): Expected number of steps to return to each regime from stationarity.

Return type:

Example

>>> result = fit_gaussian_hmm(returns, n_states=2)
>>> analysis = regime_transition_analysis(
...     result['states'], result['transition_matrix']
... )
>>> print(f"Bull market lasts ~{analysis['avg_duration'][0]:.0f} days")
>>> print(f"Steady state: {analysis['steady_state']}")

See also

fit_gaussian_hmm: Fit HMM to obtain states and transition matrix. regime_statistics: Financial statistics per regime.

gaussian_mixture_regimes(returns, n_components=2, covariance_type='full', n_init=10, random_state=42)[source]¶

Classify returns into regimes using a Gaussian Mixture Model.

Unlike HMMs, GMMs do not model temporal dependence between states. Each observation is independently assigned to the most likely Gaussian component. This makes GMMs:

Faster to fit than HMMs.
Simpler to interpret (no transition dynamics).
Less powerful for capturing regime persistence.

Good for quick regime labeling or as a baseline before fitting an HMM.

The model:

\[p(r_t) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(r_t; \mu_k, \Sigma_k)\]

Parameters:

returns (Series | ndarray) – Return series. NaN values are dropped.
n_components (int, default: 2) – Number of mixture components (regimes).
covariance_type (str, default: 'full') – Covariance type ('full', 'tied', 'diag', 'spherical').
n_init (int, default: 10) – Number of random restarts.
random_state (int, default: 42) – Random seed.

Returns:

states (np.ndarray): Component assignments, shape (T,).
state_probs (np.ndarray): Posterior probabilities, shape (T, K).
means (np.ndarray): Component means, shape (K,).
covariances (np.ndarray): Component variances, shape (K,).
weights (np.ndarray): Mixing proportions, shape (K,).
aic (float): AIC.
bic (float): BIC.
model (GaussianMixture): Fitted sklearn model.
index (pd.Index | None): Preserved index.

Return type:

Example

>>> result = gaussian_mixture_regimes(returns, n_components=2)
>>> print(result['means'])   # Component means
>>> print(result['weights']) # Mixing proportions
>>> # Combine with regime_statistics for financial interpretation
>>> stats = regime_statistics(returns, result['states'])

See also

fit_gaussian_hmm: HMM with temporal regime dynamics. regime_statistics: Per-regime financial statistics.

rolling_regime_probability(returns, n_states=2, window=None, min_window=60, n_init=5)[source]¶

Compute time-varying regime probabilities using rolling/expanding HMM.

Fits a Gaussian HMM at each time step using either a rolling or expanding window, producing a time series of regime probabilities that adapts to new data. Useful for real-time regime monitoring.

Parameters:

returns (Series) – Return series.
n_states (int, default: 2) – Number of hidden states.
window (int | None, default: None) – Rolling window size. If None, uses an expanding window (all data up to time t).
min_window (int, default: 60) – Minimum number of observations required before fitting the first model.
n_init (int, default: 5) – Number of EM restarts per window.

Return type:

Returns:

pd.DataFrame with columns prob_0, prob_1, …, prob_{K-1} indexed by the return series index. Rows before min_window are NaN.

Example

>>> probs = rolling_regime_probability(returns, n_states=2, window=252)
>>> probs.plot(title="Regime Probabilities Over Time")

Notes

This is computationally expensive. For 1000 observations with window=252, approximately 750 HMM fits are performed. Consider using a large window or expanding mode for efficiency.

See also

fit_gaussian_hmm: Single-shot HMM fitting.

regime_aware_portfolio(regime_probs, regime_weights)[source]¶

Compute blended portfolio weights using regime probabilities.

The blended weight for asset j is:

\[w_j = \sum_{k=1}^{K} P(s_t = k) \cdot w_j^{(k)}\]

Parameters:

regime_probs (ndarray) – Current regime probabilities, shape (K,). Must sum to 1.
regime_weights (ndarray) – Optimal weights for each regime, shape (K, n_assets). Row k contains the weights for regime k.

Returns:

Blended portfolio weights, shape (n_assets,). Sums to 1 if each row of regime_weights sums to 1.

Return type:

ndarray

Example

>>> # 60% probability of bull, 40% probability of bear
>>> probs = np.array([0.6, 0.4])
>>> # Bull: 80% equity / 20% bonds; Bear: 30% equity / 70% bonds
>>> weights = np.array([[0.8, 0.2], [0.3, 0.7]])
>>> blended = regime_aware_portfolio(probs, weights)
>>> print(blended)  # [0.6*0.8 + 0.4*0.3, 0.6*0.2 + 0.4*0.7]
[0.6  0.4]

See also

fit_gaussian_hmm: Obtain regime probabilities. rolling_regime_probability: Time-varying regime probabilities.

fit_ms_autoregression(endog, k_regimes=2, order=1, switching_ar=True, switching_variance=True, exog_tvtp=None)[source]¶

Fit a Markov-switching autoregression (MS-AR) model.

Hamilton (1989) framework:

\[y_t = \mu_{s_t} + \phi_{1,s_t}\,(y_{t-1} - \mu_{s_{t-1}}) + \cdots + \phi_{p,s_t}\,(y_{t-p} - \mu_{s_{t-p}}) + \sigma_{s_t}\,\epsilon_t\]

When to use MS-AR vs. plain HMM:

Use MS-AR when the series exhibits serial correlation within each regime (e.g., GDP growth, interest rates, momentum returns). The autoregressive terms capture local dynamics that a plain HMM ignores.
Use a plain Gaussian HMM (fit_gaussian_hmm) when returns are approximately i.i.d. within each regime and only the mean and variance shift. This is typical for daily equity returns.
Use MS-Regression (fit_ms_regression) when you want regime-dependent intercepts/exogenous coefficients but no AR lags.

Parameters:

endog (Series | ndarray) – Endogenous (dependent) variable. Typically a return or growth-rate series.
k_regimes (int, default: 2) – Number of regimes (hidden states). 2 is the most common choice (e.g., expansion/recession).
order (int, default: 1) – Autoregressive lag order p. Higher orders capture longer memory within each regime.
switching_ar (bool, default: True) – If True, the AR coefficients switch across regimes. If False, only the intercept and/or variance switch.
switching_variance (bool, default: True) – If True, the innovation variance switches across regimes. Almost always desirable for financial data.
exog_tvtp (ndarray | DataFrame | None, default: None) – Optional exogenous variables for time-varying transition probabilities (TVTP). Shape (T, q). When provided, the transition probabilities are modelled as logistic functions of these variables, allowing macro indicators to influence regime-switching dynamics.

Returns:

smoothed_probs (np.ndarray): Smoothed regime probabilities, shape (T, K).
filtered_probs (np.ndarray): Filtered regime probabilities, shape (T, K).
states (np.ndarray): Most likely state sequence (argmax of smoothed probabilities), shape (T,).
transition_matrix (np.ndarray): Transition matrix, shape (K, K).
expected_durations (np.ndarray): Expected duration in each regime, 1 / (1 - P[k,k]), shape (K,).
regime_params (dict): Per-regime parameters including 'mean_k', 'sigma2_k', and 'ar_k_j' (AR coefficient j for regime k).
aic (float): Akaike Information Criterion.
bic (float): Bayesian Information Criterion.
summary (str): Text summary from statsmodels.
model_result: The fitted statsmodels result object.

Return type:

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> # Simulate AR(1) with switching mean and persistence
>>> returns = pd.Series(np.concatenate([
...     rng.normal(0.002, 0.008, 250),
...     rng.normal(-0.001, 0.020, 250),
... ]))
>>> result = fit_ms_autoregression(returns, k_regimes=2, order=1)
>>> print(result['transition_matrix'])
>>> print(result['regime_params'])

Notes

Requires statsmodels >= 0.13. Convergence can be slow for higher-order AR models or many regimes. The Hamilton filter is used for inference.

References

Hamilton, J. D. (1989). “A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle.” Econometrica, 57(2), 357-384.

See also

fit_ms_regression: Markov-switching regression (no AR terms). fit_gaussian_hmm: HMM without autoregressive dynamics. select_n_states: BIC-based state selection.

select_n_states(returns, max_states=5, n_init=10, n_iter=200, random_state=42)[source]¶

Select the optimal number of HMM states using the BIC criterion.

Fits Gaussian HMMs with n_states from 2 to max_states and selects the model that minimises the Bayesian Information Criterion (BIC).

Why BIC over AIC for HMM state selection:

The BIC imposes a stronger penalty on model complexity than the AIC, which is critical for HMM selection because:

HMM parameter count grows quadratically with the number of states (transition matrix has K*(K-1) free parameters), so without a strong penalty the AIC tends to overfit by selecting too many states.
BIC is consistent: as sample size grows, BIC selects the true model (if it is in the candidate set). AIC tends to overfit even asymptotically.
In practice, financial regime models with 2-3 states are usually adequate; BIC correctly penalises 4+ state models that capture noise rather than true regimes.

The BIC is defined as:

\[\text{BIC} = k \ln(T) - 2 \ln(\hat{L})\]

where k is the number of free parameters, T is the sample size, and \(\hat{L}\) is the maximised likelihood.

Parameters:

returns (Series | ndarray) – Return series. NaN values are dropped.
max_states (int, default: 5) – Maximum number of states to consider (inclusive). The search covers [2, 3, ..., max_states].
n_init (int, default: 10) – Number of random restarts per model fit.
n_iter (int, default: 200) – Maximum EM iterations per restart.
random_state (int, default: 42) – Base random seed.

Returns:

optimal_n_states (int): The number of states with the lowest BIC.
scores (dict[int, float]): Mapping from n_states to its BIC value.
best_model (dict): The full result dict from fit_gaussian_hmm for the optimal model.

Return type:

Example

>>> result = select_n_states(returns, max_states=5)
>>> print(f"Optimal states: {result['optimal_n_states']}")
>>> print(result['scores'])
{2: -2345.6, 3: -2310.1, 4: -2290.3, 5: -2275.0}

See also

fit_gaussian_hmm: Fit a single HMM with a fixed number of states.

fit_multivariate_hmm(returns, n_states=2, covariance_type='full', n_init=10, n_iter=200, random_state=42)[source]¶

Fit a multivariate Gaussian HMM for multi-asset regime detection.

Extends the univariate HMM to a vector of asset returns, enabling detection of regimes that manifest as shifts in cross-asset correlations, not just individual-asset volatility.

Why multivariate regime detection matters:

Typical findings on equity/bond data:

Low-vol regime: Low equity vol, low bond vol, moderate negative equity-bond correlation (diversification works).
High-vol regime: High equity vol, elevated bond vol, correlations shift (equities become more correlated with each other, equity-bond correlation may flip sign).

The model estimates per-regime:

Mean return vector \(\mu_k \in \mathbb{R}^d\)
Full covariance matrix \(\Sigma_k \in \mathbb{R}^{d \times d}\)
From which per-regime correlation matrices are derived.

Parameters:

returns (DataFrame) – DataFrame of multi-asset returns, shape (T, d) where d is the number of assets. Column names are preserved for labeling. NaN rows are dropped.
n_states (int, default: 2) – Number of hidden states (regimes).
covariance_type (str, default: 'full') – Covariance parameterisation. 'full' is strongly recommended for capturing correlation shifts.
n_init (int, default: 10) – Number of random EM restarts.
n_iter (int, default: 200) – Maximum EM iterations per restart.
random_state (int, default: 42) – Random seed.

Returns:

states (np.ndarray): Most likely state sequence (Viterbi), shape (T,).
state_probs (np.ndarray): Smoothed state probabilities, shape (T, K).
means (np.ndarray): Per-regime mean vectors, shape (K, d).
covariances (np.ndarray): Per-regime covariance matrices, shape (K, d, d).
per_regime_correlations (dict[int, np.ndarray]): Per-regime correlation matrices, keyed by regime index.
transition_matrix (np.ndarray): Transition matrix (K, K).
log_likelihood (float): Total log-likelihood.
aic (float): AIC.
bic (float): BIC.
model: The fitted hmmlearn model.
index (pd.Index | None): Preserved DatetimeIndex.
columns (list[str]): Asset names.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> # Low-vol regime: low correlation
>>> cov_low = np.array([[0.0001, 0.00002], [0.00002, 0.00015]])
>>> low = rng.multivariate_normal([0.001, 0.0008], cov_low, 300)
>>> # High-vol regime: high correlation
>>> cov_high = np.array([[0.0009, 0.0006], [0.0006, 0.0008]])
>>> high = rng.multivariate_normal([-0.001, -0.0005], cov_high, 200)
>>> df = pd.DataFrame(np.vstack([low, high]), columns=['SPY', 'QQQ'])
>>> result = fit_multivariate_hmm(df, n_states=2)
>>> print(result['per_regime_correlations'][0])
>>> print(result['per_regime_correlations'][1])

See also

fit_gaussian_hmm: Univariate HMM. regime_conditional_moments: Per-regime moments from states.

regime_conditional_moments(returns, states)[source]¶

Compute per-regime mean vector and covariance/correlation matrices.

Parameters:

returns (DataFrame) – DataFrame of multi-asset returns, shape (T, d).
states (ndarray) – Integer regime labels, shape (T,), aligned with returns.

Return type:

dict[int, dict[str, ndarray]]

Returns:

Dictionary keyed by regime index (int), where each value is a dict with:

mean (np.ndarray): Sample mean vector, shape (d,).
cov (np.ndarray): Sample covariance matrix, shape (d, d).
corr (np.ndarray): Sample correlation matrix, shape (d, d).

Example

>>> result = fit_multivariate_hmm(returns_df, n_states=2)
>>> moments = regime_conditional_moments(
...     returns_df, result['states']
... )
>>> # Use regime-0 covariance for portfolio optimisation
>>> cov_bull = moments[0]['cov']
>>> mu_bull = moments[0]['mean']

See also

fit_multivariate_hmm: Multi-asset regime detection. regime_aware_portfolio: Blend weights across regimes.

Kalman Filtering¶

Kalman filter, smoother, and time-varying regression.

This module provides pure-numpy implementations of the Kalman filter, Rauch-Tung-Striebel smoother, time-varying coefficient regression (dynamic linear model), and the Unscented Kalman Filter for nonlinear state estimation.

When to use each:

kalman_filter: Linear state-space models where you want filtered (causal) state estimates. Core building block.
kalman_smoother: When you have all data and want the best possible state estimates (non-causal, uses future data).
kalman_regression: Estimate time-varying betas (e.g., hedge ratios, factor exposures) that evolve as random walks.
unscented_kalman: Nonlinear dynamics or observation models where the standard Kalman filter’s linearity assumption breaks down.

References

Durbin, J. & Koopman, S. J. (2012). Time Series Analysis by State Space Methods. Oxford University Press.

Haykin, S. (2001). Kalman Filtering and Neural Networks. Wiley.

Julier, S. J. & Uhlmann, J. K. (2004). “Unscented Filtering and Nonlinear Estimation.” Proceedings of the IEEE, 92(3).

kalman_filter(observations, F, H, Q, R, x0, P0)[source]¶

Run a linear Kalman filter over a sequence of observations.

The Kalman filter is the optimal recursive estimator for a linear Gaussian state-space model:

\[ \begin{align}\begin{aligned}x_t = F \, x_{t-1} + w_t, \quad w_t \sim N(0, Q)\\y_t = H \, x_t + v_t, \quad v_t \sim N(0, R)\end{aligned}\end{align} \]

At each time step, the filter produces:

Prediction: \(\hat{x}_{t|t-1}\) and \(P_{t|t-1}\)
Update: \(\hat{x}_{t|t}\) and \(P_{t|t}\) after incorporating observation \(y_t\)

Financial applications:

Estimating latent factors or unobserved states
Noise reduction (e.g., estimating “true” volatility from noisy proxies)
Preprocessing for regime detection

Parameters:

observations (ndarray) – Observation matrix of shape (T, m) where T is the number of time steps and m is the observation dimension. For univariate observations, shape (T, 1).
F (ndarray) – State transition matrix (n, n).
H (ndarray) – Observation matrix (m, n).
Q (ndarray) – Process noise covariance (n, n).
R (ndarray) – Observation noise covariance (m, m).
x0 (ndarray) – Initial state estimate (n,) or (n, 1).
P0 (ndarray) – Initial state covariance (n, n).

Returns:

filtered_states (np.ndarray): Filtered state estimates, shape (T, n).
filtered_covs (np.ndarray): Filtered covariance matrices, shape (T, n, n).
predicted_states (np.ndarray): One-step-ahead predicted states, shape (T, n).
predicted_covs (np.ndarray): One-step-ahead predicted covariances, shape (T, n, n).
innovations (np.ndarray): Innovation (residual) sequence, shape (T, m).
innovation_covs (np.ndarray): Innovation covariances, shape (T, m, m).
kalman_gains (np.ndarray): Kalman gain matrices, shape (T, n, m).
log_likelihood (float): Total log-likelihood computed from the innovation sequence.

Return type:

Example

>>> # Track a random walk with noisy observations
>>> F = np.array([[1.0]])
>>> H = np.array([[1.0]])
>>> Q = np.array([[0.01]])   # Process noise
>>> R = np.array([[0.1]])    # Observation noise
>>> x0 = np.array([0.0])
>>> P0 = np.array([[1.0]])
>>> obs = np.cumsum(np.random.randn(100, 1) * 0.1) + \
...       np.random.randn(100, 1) * 0.3
>>> result = kalman_filter(obs, F, H, Q, R, x0, P0)
>>> print(result['log_likelihood'])
>>> print(result['filtered_states'][-5:])

See also

kalman_smoother: Backward pass for optimal smoothed estimates. kalman_regression: Time-varying coefficient regression. unscented_kalman: Nonlinear state estimation.

kalman_smoother(observations, F, H, Q, R, x0, P0)[source]¶

Run the Kalman filter followed by the RTS backward smoother.

When to use: Offline analysis where you have the complete time series and want the best possible state estimates.

Parameters:

observations (ndarray) – Observations, shape (T, m).
F (ndarray) – State transition matrix (n, n).
H (ndarray) – Observation matrix (m, n).
Q (ndarray) – Process noise covariance (n, n).
R (ndarray) – Observation noise covariance (m, m).
x0 (ndarray) – Initial state (n,).
P0 (ndarray) – Initial covariance (n, n).

Returns:

smoothed_states (np.ndarray): Smoothed state estimates, shape (T, n).
smoothed_covs (np.ndarray): Smoothed covariance matrices, shape (T, n, n).

Return type:

Example

>>> result = kalman_smoother(obs, F, H, Q, R, x0, P0)
>>> # Smoothed states are more accurate than filtered
>>> print(result['smoothed_states'][-5:])
>>> print(result['filtered_states'][-5:])

See also

kalman_filter: Forward pass only. kalman_regression: Time-varying betas using DLM.

kalman_regression(y, X, state_noise=0.0001, obs_noise=None, x0=None, P0_scale=1.0, smooth=True)[source]¶

Estimate time-varying regression coefficients using a Kalman filter.

Models the regression coefficients as a random walk:

\[ \begin{align}\begin{aligned}\beta_t = \beta_{t-1} + w_t, \quad w_t \sim N(0, Q)\\y_t = X_t' \beta_t + v_t, \quad v_t \sim N(0, R)\end{aligned}\end{align} \]

This is a Dynamic Linear Model (DLM) where the state vector is the coefficient vector. Common financial applications:

Time-varying betas: Estimate CAPM beta that changes over time
Dynamic hedge ratios: For pairs trading or hedging
Factor exposure monitoring: Track how factor loadings evolve

Parameters:

y (Series | ndarray) – Dependent variable, shape (T,).
X (DataFrame | ndarray) – Regressors, shape (T, p). Include a column of ones for an intercept.
state_noise (float, default: 0.0001) – Process noise variance for each coefficient. Controls how quickly coefficients can change. Larger values allow faster adaptation.
obs_noise (float | None, default: None) – Observation noise variance. If None, estimated as the variance of OLS residuals.
x0 (ndarray | None, default: None) – Initial coefficient estimate, shape (p,). If None, OLS estimates are used.
P0_scale (float, default: 1.0) – Scale for initial covariance (P0 = P0_scale * I).
smooth (bool, default: True) – If True, apply the RTS smoother for better estimates.

Returns:

coefficients (np.ndarray): Time-varying coefficients, shape (T, p). If smooth=True, these are smoothed.
coefficient_covs (np.ndarray): Covariance of coefficients, shape (T, p, p).
residuals (np.ndarray): Observation residuals, shape (T,).
log_likelihood (float): Log-likelihood.
filtered_coefficients (np.ndarray): Filtered (causal) coefficient estimates, shape (T, p).

Return type:

Example

>>> import numpy as np, pandas as pd
>>> # Time-varying beta estimation
>>> market = np.random.randn(500) * 0.01
>>> # True beta changes from 1.0 to 1.5 halfway
>>> true_beta = np.where(np.arange(500) < 250, 1.0, 1.5)
>>> stock = true_beta * market + np.random.randn(500) * 0.005
>>> X = np.column_stack([np.ones(500), market])
>>> result = kalman_regression(stock, X, state_noise=1e-5)
>>> betas = result['coefficients'][:, 1]  # Time-varying beta
>>> print(f"Beta at t=100: {betas[100]:.2f}")  # ~1.0
>>> print(f"Beta at t=400: {betas[400]:.2f}")  # ~1.5

Notes

The state_noise parameter is critical. Too large and the estimates are noisy; too small and they adapt too slowly. A good starting point is 1e-4 to 1e-6 for daily returns.

See also

kalman_filter: General-purpose Kalman filter. kalman_smoother: RTS backward smoother.

unscented_kalman(observations, state_dim, f_func, h_func, Q, R, x0, P0, alpha=0.001, beta=2.0, kappa=0.0)[source]¶

Unscented Kalman Filter for nonlinear state estimation.

The UKF uses sigma points to propagate the state distribution through nonlinear dynamics and observation functions, avoiding the need for Jacobian computation (unlike the Extended Kalman Filter).

The state-space model:

\[ \begin{align}\begin{aligned}x_t = f(x_{t-1}) + w_t, \quad w_t \sim N(0, Q)\\y_t = h(x_t) + v_t, \quad v_t \sim N(0, R)\end{aligned}\end{align} \]

When to use: When f or h are nonlinear. Examples include stochastic volatility models, nonlinear factor models, or target tracking with angular observations.

Parameters:

observations (ndarray) – Observations, shape (T, m).
state_dim (int) – Dimension of the state vector (n).
f_func (Callable[[ndarray], ndarray]) – State transition function. Takes x of shape (n,) and returns predicted state of shape (n,).
h_func (Callable[[ndarray], ndarray]) – Observation function. Takes x of shape (n,) and returns predicted observation of shape (m,).
Q (ndarray) – Process noise covariance (n, n).
R (ndarray) – Observation noise covariance (m, m).
x0 (ndarray) – Initial state estimate (n,).
P0 (ndarray) – Initial covariance (n, n).
alpha (float, default: 0.001) – Spread of sigma points (typically 1e-4 to 1).
beta (float, default: 2.0) – Prior knowledge of distribution (2 is optimal for Gaussian).
kappa (float, default: 0.0) – Secondary scaling parameter (0 or 3 - n).

Returns:

filtered_states (np.ndarray): Filtered state estimates, shape (T, n).
filtered_covs (np.ndarray): Filtered covariances, shape (T, n, n).
innovations (np.ndarray): Innovation sequence, shape (T, m).
log_likelihood (float): Approximate log-likelihood.

Return type:

Example

>>> # Nonlinear observation model: observe angle to a target
>>> def f(x):
...     return x  # Random walk dynamics
>>> def h(x):
...     return np.array([np.arctan2(x[1], x[0])])  # Angle
>>> result = unscented_kalman(
...     observations, state_dim=2, f_func=f, h_func=h,
...     Q=np.eye(2)*0.01, R=np.array([[0.1]]),
...     x0=np.zeros(2), P0=np.eye(2),
... )

Notes

The UKF with the standard unscented transform uses 2n + 1 sigma points. Computational cost scales as O(n^2) per step (vs O(n^3) for EKF with Jacobian).

See also

kalman_filter: Linear Kalman filter. kalman_smoother: Linear Kalman smoother.

Changepoint Detection¶

Change-point detection methods for structural break analysis.

Provides Bayesian online change-point detection (Adams & MacKay 2007), PELT (Pruned Exact Linear Time), binary segmentation, sliding-window change-point detection, and CUSUM control charts. These tools identify abrupt shifts in the statistical properties of a time series – a critical prerequisite for regime-aware investing.

online_changepoint(data, hazard=0.005)[source]¶

Bayesian online change-point detection.

Implements the algorithm from Adams & MacKay (2007). Uses a normal model with unknown mean and known variance (estimated from data).

Parameters:

data (Series) – Time series to monitor.
hazard (float, default: 0.005) – Prior probability that a change-point occurs at each time step (constant hazard rate 1/mean_run_length).

Return type:

Returns:

Series of estimated run lengths (the most probable run length at each time step). A drop to zero indicates a detected change-point.

pelt_changepoint(data, *, penalty='bic', model='l2', min_size=5, jump=1, n_bkps=None)[source]¶

Pruned Exact Linear Time (PELT) optimal changepoint detection.

The algorithm minimises:

\[\sum_{i=0}^{K} C(y_{\tau_i : \tau_{i+1}}) + \beta \, K\]

where C is the segment cost, K is the number of changepoints, and beta is the penalty.

Interpretation guidance:

A changepoint at index t means the statistical properties of the series changed between observation t-1 and t.
segment_stats gives the mean, variance, and length of each segment so you can characterise the regimes between breaks.
confidence is a heuristic based on the cost reduction at each changepoint relative to the global cost. Higher values indicate more confident breaks.

Parameters:

data (Series | ndarray) – Time series to analyse. If a pd.Series, the index is preserved in the output.
penalty (str | float, default: 'bic') –
Penalty method or explicit float value. Supported strings:
- "bic" – Bayesian Information Criterion (default). Works well when the number of changepoints is unknown.
- "aic" – Akaike Information Criterion. Tends to over-segment; prefer BIC unless you want sensitivity.
- "mbic" – Modified BIC (Liu et al., 2016).
If a float is given, it is used directly as the penalty value beta.
model (str, default: 'l2') –
Cost function for each segment. Options:
- "l2" – least-squares cost (Gaussian change in mean).
- "l1" – absolute deviation cost (robust to outliers).
- "rbf" – kernel cost (nonparametric, captures complex distribution changes).
- "linear" – linear regression cost (for trend breaks).
- "normal" – full Gaussian cost (mean and variance).
- "ar" – autoregressive model cost (order-1 AR).
min_size (int, default: 5) – Minimum segment length. Shorter segments are prohibited. Increase to suppress spurious breaks.
jump (int, default: 1) – Subsample factor for candidate changepoints. jump=1 considers every index; larger values trade accuracy for speed on very long series.
n_bkps (int | None, default: None) – If given, ignore penalty and find exactly this many changepoints (uses the Pelt predict with n_bkps).

Returns:

changepoints (list[int]): Detected changepoint indices (0-based). The last element is always the series length.
n_changepoints (int): Number of detected changepoints (excludes the trailing length marker).
segment_stats (list[dict]): Per-segment statistics with keys start, end, mean, var, length.
confidence (np.ndarray): Heuristic confidence score for each changepoint (between 0 and 1).
penalty_value (float): Effective penalty value used.
model (str): Cost model used.
cost (float): Total cost of the segmentation.

Return type:

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(5, 1, 200),
...     rng.normal(0, 2, 200),
... ])
>>> result = pelt_changepoint(pd.Series(data))
>>> print(result["changepoints"])  # ~ [200, 400, 600]
>>> print(result["segment_stats"])

References

Killick, R., Fearnhead, P. & Eckley, I. A. (2012). “Optimal Detection of Changepoints with a Linear Computational Cost.” J. Amer. Statist. Assoc., 107(500).

See also

binary_segmentation: Hierarchical top-down changepoint detection. window_changepoint: Sliding-window change score time series. online_changepoint: Bayesian online detection (sequential).

binary_segmentation(data, *, n_bkps=None, penalty='bic', model='l2', min_size=5, jump=1)[source]¶

Binary segmentation with hierarchical changepoint ordering.

Interpretation guidance:

hierarchical_order[0] is the single most important changepoint in the series. If you can only act on one break, act on this one.
The ordering degrades gracefully: if you want K changepoints, take the first K entries of hierarchical_order.

Parameters:

data (Series | ndarray) – Time series to analyse.
n_bkps (int | None, default: None) – Number of changepoints to detect. If None, the penalty is used to determine the number automatically.
penalty (str | float, default: 'bic') – Penalty method or explicit float. Same options as pelt_changepoint(). Ignored when n_bkps is set.
model (str, default: 'l2') – Cost function ("l1", "l2", "rbf", "linear", "normal").
min_size (int, default: 5) – Minimum segment length.
jump (int, default: 1) – Candidate subsample factor.

Returns:

changepoints (list[int]): Detected changepoints (last element is always series length).
n_changepoints (int): Number of changepoints.
hierarchical_order (list[int]): Changepoints sorted from most significant (largest cost improvement) to least.
segment_stats (list[dict]): Per-segment start, end, mean, var, length.
model (str): Cost model used.

Return type:

Example

>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(5, 2, 100),
...     rng.normal(0, 1, 200),
... ])
>>> result = binary_segmentation(pd.Series(data), n_bkps=2)
>>> # First entry is the dominant changepoint
>>> print(result["hierarchical_order"])

References

Scott, A. J. & Knott, M. (1974). “A Cluster Analysis Method for Grouping Means in the Analysis of Variance.” Biometrics.

See also

pelt_changepoint: Globally optimal changepoint detection. window_changepoint: Sliding-window approach.

window_changepoint(data, *, width=50, model='l2', min_size=5, jump=1, n_bkps=None, penalty='bic')[source]¶

Sliding-window changepoint detection with change-score time series.

Interpretation guidance:

The change_score series is analogous to a “structural break intensity” indicator. Threshold it or combine it with other signals to trigger regime-switch alerts.
The score is non-negative and unitless; higher values indicate stronger evidence for a distributional shift.

Parameters:

data (Series | ndarray) – Time series to analyse.
width (int, default: 50) – Half-width of the sliding window. Total window is 2 * width. Larger windows are more robust but less responsive.
model (str, default: 'l2') – Cost function ("l2", "l1", "rbf", etc.).
min_size (int, default: 5) – Minimum segment length.
jump (int, default: 1) – Candidate subsample factor.
n_bkps (int | None, default: None) – If given, return exactly this many changepoints.
penalty (str | float, default: 'bic') – Penalty method or float. Used when n_bkps is None.

Returns:

change_score (pd.Series or np.ndarray): Change score at each time step. Same type as input.
changepoints (list[int]): Detected changepoint indices.
n_changepoints (int): Number of changepoints.
width (int): Window half-width used.

Return type:

Example

>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(3, 1, 200),
... ])
>>> result = window_changepoint(pd.Series(data), width=30)
>>> score = result["change_score"]
>>> # Plot score to see the peak near index 200

References

Truong, C., Oudre, L. & Vayer, N. (2020). “Selective review of offline change point detection methods.” Signal Processing.

See also

pelt_changepoint: Offline optimal detection. online_changepoint: Bayesian online detection.

cusum_control_chart(data, *, target=None, std_est=None, k=0.5, h=5.0)[source]¶

CUSUM (Cumulative Sum) control chart for process monitoring.

The CUSUM chart monitors the cumulative sum of deviations from a target value. It is one of the most powerful sequential tests for detecting small, persistent shifts in the mean of a process.

The upper and lower CUSUM statistics are:

\[ \begin{align}\begin{aligned}S^+_t = \max(0,\; S^+_{t-1} + z_t - k)\\S^-_t = \max(0,\; S^-_{t-1} - z_t - k)\end{aligned}\end{align} \]

Interpretation guidance:

upper_cusum rising indicates the process mean is drifting above target. lower_cusum rising indicates drift below target.
alarm_points are the times at which the CUSUM crossed the control limit. These are your structural break candidates.
arl (average run length) is the expected number of observations between false alarms under the null hypothesis. Higher ARL means fewer false positives. With k = 0.5 and h = 5, the in-control ARL is approximately 465 (for a standard normal process).

Parameters:

data (Series | ndarray) – Time series to monitor.
target (float | None, default: None) – Target (in-control) mean. If None, the overall sample mean is used.
std_est (float | None, default: None) – Standard deviation estimate. If None, the sample standard deviation is used.
k (float, default: 0.5) – Allowance (slack) parameter, in units of standard deviation. Controls the size of the shift you want to detect. Default 0.5 detects a 1-sigma shift optimally.
h (float, default: 5.0) – Decision interval (control limit), in units of standard deviation. Larger h reduces false alarms but increases detection delay.

Returns:

upper_cusum (np.ndarray): Upper CUSUM statistic at each time step.
lower_cusum (np.ndarray): Lower CUSUM statistic.
control_limit (float): The h threshold used.
alarm_points (list[int]): Indices where an alarm was triggered (CUSUM exceeded h).
alarm_direction (list[str]): Direction of each alarm ("upper" or "lower").
arl (float): Estimated average run length (between alarms) based on the observed alarm sequence. inf if no alarms were triggered.
target (float): Target mean used.
std (float): Standard deviation estimate used.

Return type:

Example

>>> rng = np.random.default_rng(42)
>>> # In-control for 200 obs, then shift mean by +1 sigma
>>> data = np.concatenate([
...     rng.normal(0, 1, 200),
...     rng.normal(1, 1, 200),
... ])
>>> result = cusum_control_chart(pd.Series(data), target=0, std_est=1)
>>> print(result["alarm_points"])  # first alarm near index 200
>>> print(f"ARL: {result['arl']:.0f}")

References

Page, E. S. (1954). “Continuous Inspection Schemes.” Biometrika, 41(1/2).

Montgomery, D. C. (2009). Introduction to Statistical Quality Control. 6th ed. Wiley.

See also

online_changepoint: Bayesian approach to sequential detection. window_changepoint: Sliding-window approach.

Regime Scoring¶

Regime quality assessment and method comparison.

Evaluates how good a set of detected regimes actually is by measuring stability, separability, and predictability. These scores help you decide whether your regime model is capturing meaningful market structure or just fitting noise.

Three orthogonal quality dimensions¶

Stability (regime_stability_score()): Are the regimes persistent? Rapidly flickering regimes are hard to trade.
Separation (regime_separation_score()): Are the regimes statistically distinct? Overlapping distributions offer no actionable information.
Predictability (regime_predictability()): Can regime transitions be forecast? A model is useful only if it gives you a head-start on the next shift.

Use compare_regime_methods() to run all three metrics across multiple detection algorithms on the same data.

regime_stability_score(states, transition_matrix=None)[source]¶

Measure how stable the detected regimes are.

How it works:

Three sub-metrics are combined (equal-weighted average):

Average duration (normalised): longer average stays = more stable. Normalised as 1 - 1/avg_duration so it maps to [0, 1].
Transition frequency (inverse): fewer transitions per unit time = more stable. Computed as 1 - n_transitions / (T-1).
Probability entropy (inverse): when the transition matrix rows have low entropy (i.e., dominated by the diagonal), the chain is persistent. Normalised by log(K) where K is the number of regimes.

Interpretation guidance:

> 0.8: Highly stable regimes, suitable for discrete portfolio switches.
0.5 – 0.8: Moderate stability, consider probability-weighted blending rather than hard switches.
< 0.5: Unstable, likely overfitting or too many regimes.

Parameters:

states (ndarray | Series) – Integer regime labels, shape (T,).
transition_matrix (ndarray | None, default: None) – Optional (K, K) transition matrix. If None, an empirical matrix is computed from states.

Returns:

stability_score (float): Composite stability score in [0, 1].
avg_duration (float): Average regime duration in periods.
transition_frequency (float): Fraction of time steps with a regime change.
probability_entropy (float): Mean row entropy of the transition matrix (nats).
n_transitions (int): Total number of regime transitions.
per_regime_duration (dict[int, float]): Average duration per regime.

Return type:

Example

>>> states = np.array([0]*100 + [1]*100 + [0]*100)
>>> result = regime_stability_score(states)
>>> print(f"Stability: {result['stability_score']:.2f}")
0.99

See also

regime_separation_score: Distribution distance between regimes. compare_regime_methods: Multi-method comparison.

regime_separation_score(returns, states)[source]¶

Measure how well-separated the regime distributions are.

How it works:

For each pair of regimes (i, j):

Bhattacharyya distance: measures divergence between two Gaussian distributions. Larger = more separated.
Overlap coefficient: fraction of probability mass that overlaps when the two regime distributions are superimposed. Smaller = more separated.

The composite score is the average pairwise Bhattacharyya distance, normalised to [0, 1] via a sigmoid transform.

Interpretation guidance:

> 0.7: Well-separated regimes with distinct return distributions.
0.3 – 0.7: Moderate separation. Regimes differ but there is meaningful overlap.
< 0.3: Poor separation. Consider reducing the number of regimes.

Parameters:

returns (ndarray | Series) – Return series, shape (T,).
states (ndarray | Series) – Integer regime labels, shape (T,).

Returns:

separation_score (float): Composite separation in [0, 1].
pairwise_bhattacharyya (dict[tuple, float]): Bhattacharyya distance for each regime pair (i, j).
pairwise_overlap (dict[tuple, float]): Overlap coefficient for each pair.
per_regime_stats (dict[int, dict]): Mean and std per regime.

Return type:

Example

>>> rng = np.random.default_rng(0)
>>> returns = np.concatenate([
...     rng.normal(0.01, 0.01, 200),
...     rng.normal(-0.02, 0.03, 200),
... ])
>>> states = np.array([0]*200 + [1]*200)
>>> result = regime_separation_score(returns, states)
>>> print(f"Separation: {result['separation_score']:.2f}")

See also

regime_stability_score: Temporal persistence of regimes. compare_regime_methods: Multi-method comparison.

regime_predictability(states, transition_matrix=None, test_fraction=0.3)[source]¶

Assess whether regime transitions can be predicted.

How it works:

Split the state sequence into train / test at (1 - test_fraction) * T.
Estimate the transition matrix from the training set.
For each test observation, predict the next state as the most probable successor under the estimated transition matrix.
Compute out-of-sample accuracy and compare to a naive baseline (always predict the most common state).

Interpretation guidance:

accuracy > baseline + 0.1: The transition structure adds predictive value. Use regime probabilities to adjust portfolio weights.
accuracy ~ baseline: No predictive value from transitions. Consider whether the regime labeling is too noisy or the number of states is wrong.
economic_value: A rough estimate of the annualised Sharpe improvement from perfect regime prediction vs buy-and-hold. This uses per-regime Sharpe ratios and is indicative only.

Parameters:

states (ndarray | Series) – Integer regime labels, shape (T,).
transition_matrix (ndarray | None, default: None) – If provided, used instead of the training-set estimate.
test_fraction (float, default: 0.3) – Fraction of data reserved for testing.

Returns:

predictability_score (float): Composite predictability in [0, 1]. Combines accuracy lift and transition regularity.
accuracy (float): Out-of-sample one-step prediction accuracy.
baseline_accuracy (float): Accuracy of always predicting the most common state.
accuracy_lift (float): accuracy - baseline_accuracy.
transition_matrix_train (np.ndarray): Transition matrix estimated from training data.

Return type:

Example

>>> # Highly predictable two-state chain
>>> rng = np.random.default_rng(0)
>>> states = np.zeros(500, dtype=int)
>>> for t in range(1, 500):
...     if states[t-1] == 0:
...         states[t] = 0 if rng.random() < 0.95 else 1
...     else:
...         states[t] = 1 if rng.random() < 0.90 else 0
>>> result = regime_predictability(states)
>>> print(f"Accuracy: {result['accuracy']:.2f}")

See also

regime_stability_score: Regime persistence. regime_separation_score: Distribution distinctness. compare_regime_methods: Multi-method comparison.

compare_regime_methods(returns, methods=None, method_kwargs=None)[source]¶

Compare multiple regime detection methods on the same data.

Runs each method on the provided return series and evaluates stability, separation, and predictability. Returns a DataFrame for easy comparison.

Interpretation guidance:

Use this to pick the best regime model for your asset or strategy. The “best” model depends on your use-case: high stability for discrete allocation switches, high separation for regime-conditional signals, high predictability for transition timing.

Parameters:

returns (Series | ndarray) – Return series, shape (T,).
methods (dict[str, Callable[..., Any]] | None, default: None) –
Dictionary mapping method names to callables. Each callable must accept returns and return a tuple (states, transition_matrix) where states is an integer array and transition_matrix is (K, K) or None.

If None, a default set is used (requires detect_regimes from wraquant.regimes.base).
method_kwargs (dict[str, dict[str, Any]] | None, default: None) – Per-method keyword arguments. Keys should match methods keys.

Returns:

stability, separation, predictability
avg_duration, n_transitions, n_regimes
composite (equal-weighted average of the three scores)

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = np.concatenate([
...     rng.normal(0.01, 0.01, 200),
...     rng.normal(-0.02, 0.03, 200),
... ])
>>> df = compare_regime_methods(pd.Series(returns))
>>> print(df.sort_values("composite", ascending=False))

See also

regime_stability_score: Stability metric details. regime_separation_score: Separation metric details. regime_predictability: Predictability metric details.

Regime Labels¶

Regime labeling, classification, and duration analysis.

Provides rule-based and statistical approaches to labeling market regimes without requiring a fitted model. These functions are useful for backtesting, for creating training labels for supervised regime classifiers, and for generating interpretable regime descriptions.

label_regimes(states, returns)[source]¶

Assign descriptive labels to numeric regime states.

States are sorted by mean return: the state with the highest mean return is labeled "bull", the lowest "bear", and any intermediate states "neutral_1", "neutral_2", etc.

Parameters:

states (Series) – Integer regime state series.
returns (Series) – Corresponding return series (same index).

Return type:

Returns:

Series of string regime labels.

regime_statistics(returns, states)[source]¶

Compute descriptive statistics for each regime.

Parameters:

returns (Series) – Return series.
states (Series) – Integer regime state series (same index).

Return type:

Returns:

DataFrame indexed by regime state with columns for mean, std, skew, count, and fraction of total observations.

volatility_regime_labels(returns, *, window=21, n_levels=3, quantiles=None)[source]¶

Label regimes based on realised volatility quantiles.

A simple, model-free approach that classifies each period by where its rolling volatility falls within the historical distribution. No fitting, no hidden states – just raw vol percentiles.

Interpretation guidance:

"low_vol" periods typically correspond to trending or complacent markets. Strategy-wise, favour momentum and carry.
"high_vol" periods correspond to stressed or mean-reverting markets. Favour defensive positioning or mean-reversion.
"medium_vol" is the transition zone.

Parameters:

returns (Series | ndarray) – Return series.
window (int, default: 21) – Rolling window for realised volatility estimation. Default 21 (roughly one trading month).
n_levels (int, default: 3) – Number of volatility levels. Default 3 produces low_vol / medium_vol / high_vol. Use 2 for a binary split or 4+ for finer granularity.
quantiles (list[float] | None, default: None) – Explicit quantile boundaries. If provided, overrides n_levels. Must have n_levels - 1 elements, each in (0, 1).

Return type:

Returns:

pd.Series of string labels (e.g., "low_vol", "medium_vol", "high_vol"). NaN-filled for the warm-up period where rolling volatility is unavailable.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(rng.normal(0, 0.01, 500))
>>> labels = volatility_regime_labels(returns, n_levels=3)
>>> print(labels.value_counts())

See also

trend_regime_labels: Label by trend direction. composite_regime_labels: Combine vol + trend labels.

trend_regime_labels(returns, *, fast_window=10, slow_window=50, hysteresis=0.0005)[source]¶

Label regimes based on moving average slope with hysteresis.

Uses a dual moving-average crossover system with a hysteresis band to avoid whipsaw signals. The result is a clean, three-state classification: uptrend, downtrend, or sideways.

Interpretation guidance:

"uptrend": Fast MA is above slow MA by more than the hysteresis threshold. Bullish bias.
"downtrend": Fast MA is below slow MA by more than the hysteresis threshold. Bearish bias.
"sideways": The two MAs are within the hysteresis band. No directional conviction – favour range-bound strategies.

Parameters:

returns (Series | ndarray) – Return series.
fast_window (int, default: 10) – Fast moving average window (periods).
slow_window (int, default: 50) – Slow moving average window (periods).
hysteresis (float, default: 0.0005) – Minimum difference between fast and slow MA (in return units) required to declare a trend. Larger values suppress whipsaws but delay signals.

Return type:

Returns:

pd.Series of string labels ("uptrend", "downtrend", "sideways"). NaN-filled during warm-up.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(rng.normal(0.001, 0.01, 500))
>>> labels = trend_regime_labels(returns)
>>> print(labels.value_counts())

See also

volatility_regime_labels: Label by vol level. composite_regime_labels: Combine vol + trend labels.

composite_regime_labels(returns, *, vol_window=21, fast_window=10, slow_window=50, hysteresis=0.0005, n_vol_levels=2)[source]¶

Combine volatility and trend regimes into composite states.

Creates 4-6 composite labels by crossing trend direction (uptrend / downtrend / sideways) with volatility level (low / high or low / medium / high). Common composite states:

bull_calm: Uptrend + low vol. The best environment for passive equity holding.
bull_volatile: Uptrend + high vol. Often late-cycle or recovery rallies.
bear_calm: Downtrend + low vol. Grinding bear markets.
bear_volatile: Downtrend + high vol. Crisis periods (2008, March 2020).
sideways_calm: Range-bound, quiet.
sideways_volatile: Choppy, difficult to trade.

Interpretation guidance:

The composite label captures both direction and turbulence, which together determine the optimal strategy. For instance, momentum strategies work in bull_calm but fail in bear_volatile.

Parameters:

returns (Series | ndarray) – Return series.
vol_window (int, default: 21) – Window for rolling volatility.
fast_window (int, default: 10) – Fast MA window for trend.
slow_window (int, default: 50) – Slow MA window for trend.
hysteresis (float, default: 0.0005) – Trend hysteresis threshold.
n_vol_levels (int, default: 2) – 2 or 3 volatility levels.

Return type:

Returns:

pd.Series of string composite labels. NaN-filled during warm-up.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = pd.Series(rng.normal(0.001, 0.01, 500))
>>> labels = composite_regime_labels(returns)
>>> print(labels.value_counts())

See also

volatility_regime_labels: Volatility-only labeling. trend_regime_labels: Trend-only labeling. regime_duration_analysis: Analyse how long each composite

state typically lasts.

regime_duration_analysis(states)[source]¶

Analyse how long each regime typically lasts.

Interpretation guidance:

survival_curve[k]: Probability that a regime-k spell lasts at least d periods. A slowly-decaying curve means the regime tends to persist.
hazard_rate[k]: Instantaneous probability of exiting regime k after having been in it for d periods. If the hazard rate is approximately constant, regime duration is memoryless (geometric distribution, consistent with Markov). An increasing hazard rate means longer spells are more likely to end soon.
expected_remaining[k]: Given that we are currently in regime k and have been for d periods, how many more periods should we expect? Computed from the empirical survival function.

Parameters:

states (Series | ndarray) – Integer regime labels, shape (T,).

Returns:

durations (dict[int, list[int]]): List of spell durations for each regime.
survival_curve (dict[int, pd.Series]): Kaplan-Meier style survival curve for each regime, indexed by duration.
hazard_rate (dict[int, pd.Series]): Empirical hazard rate for each regime, indexed by duration.
expected_remaining (dict[int, pd.Series]): Expected remaining duration conditional on having survived d periods, indexed by duration.
summary (pd.DataFrame): Per-regime summary with mean_duration, median_duration, max_duration, n_spells.

Return type:

Example

>>> states = np.array([0]*50 + [1]*30 + [0]*80 + [1]*40)
>>> result = regime_duration_analysis(states)
>>> print(result["summary"])
>>> # Survival curve for regime 0
>>> print(result["survival_curve"][0])

See also

regime_stability_score: Composite stability metric. composite_regime_labels: Generate regime labels to analyse.

Integrations¶

Advanced regime detection integrations using optional packages.

Provides wrappers around pomegranate, filterpy, river, and dynamax for Hidden Markov Models, Kalman filtering, online drift detection, and Linear Gaussian State Space Models.

pomegranate_hmm(data, n_states=2)[source]¶

Fit a Hidden Markov Model using pomegranate.

Parameters:

data (Series | ndarray) – Univariate observations (e.g. returns).
n_states (int, default: 2) – Number of hidden states.

Returns:

Dictionary containing:

states – predicted hidden state sequence (1-D int array).
means – estimated mean of each state’s emission distribution.
model – the fitted pomegranate DenseHMM object.
n_states – number of hidden states.

Return type:

filterpy_kalman(observations, F, H, Q, R, x0=None, P0=None)[source]¶

Run a Kalman filter using filterpy.

Parameters:

observations (ndarray | Series) – Observed measurements. Shape (T,) for univariate or (T, dim_z) for multivariate observations.
F (ndarray) – State transition matrix. Shape (dim_x, dim_x).
H (ndarray) – Observation (measurement) matrix. Shape (dim_z, dim_x).
Q (ndarray) – Process noise covariance. Shape (dim_x, dim_x).
R (ndarray) – Measurement noise covariance. Shape (dim_z, dim_z).
x0 (ndarray | None, default: None) – Initial state estimate. Shape (dim_x,) or (dim_x, 1). Defaults to zeros.
P0 (ndarray | None, default: None) – Initial covariance estimate. Shape (dim_x, dim_x). Defaults to identity.

Returns:

Dictionary containing:

filtered_states – filtered state estimates, shape (T, dim_x).
filtered_covariances – filtered covariance matrices, shape (T, dim_x, dim_x).
log_likelihood – total log-likelihood.
residuals – measurement residuals, shape (T, dim_z).

Return type:

river_drift_detector(stream, method='adwin')[source]¶

Detect concept drift in a data stream using river.

Processes each observation sequentially and records indices where a drift is detected.

Parameters:

stream (ndarray | Series | list[float]) – Sequential stream of numeric observations.
method (str, default: 'adwin') –
Drift detection method:
- 'adwin' – Adaptive Windowing (ADWIN)
- 'ddm' – Drift Detection Method
- 'eddm' – Early Drift Detection Method
- 'page_hinkley' – Page-Hinkley test

Returns:

Dictionary containing:

drift_indices – list of indices where drift was detected.
n_drifts – total number of drifts detected.
method – drift detection method used.

Return type:

dynamax_lgssm(observations, state_dim=2, emission_dim=None, n_iters=100)[source]¶

Fit a Linear Gaussian State Space Model using dynamax (JAX-based).

Parameters:

observations (ndarray | Series) – Observations of shape (T, obs_dim) or (T,) for univariate.
state_dim (int, default: 2) – Dimensionality of the latent state.
emission_dim (int | None, default: None) – Observation dimension. Inferred from observations if None.
n_iters (int, default: 100) – Number of EM iterations.

Returns:

Dictionary containing:

filtered_means – filtered state means, shape (T, state_dim).
filtered_covs – filtered state covariances.
smoothed_means – smoothed state means.
params – learned model parameters.
log_likelihoods – EM log-likelihood trace.

Return type:

pykalman_filter(observations, transition_matrices=None, observation_matrices=None, n_em_iter=10)[source]¶

Kalman filter with EM parameter learning via pykalman.

Parameters:

observations (ndarray | Series) – Observations of shape (T, obs_dim) or (T,) for univariate data.
transition_matrices (ndarray | None, default: None) – Initial guess for the state transition matrix. When None, pykalman learns it from data via EM.
observation_matrices (ndarray | None, default: None) – Initial guess for the observation matrix. When None, pykalman learns it from data via EM.
n_em_iter (int, default: 10) – Number of EM iterations for parameter learning.

Returns:

Dictionary containing:

filtered_means – filtered state estimates, shape (T, state_dim).
filtered_covs – filtered state covariance matrices, shape (T, state_dim, state_dim).
smoothed_means – smoothed (RTS) state estimates, shape (T, state_dim).
smoothed_covs – smoothed state covariance matrices, shape (T, state_dim, state_dim).
learned_params – dict with learned transition, observation, transition_cov, and observation_cov matrices.
log_likelihood – total log-likelihood of the data under the learned model.

Return type: