ReferenceHidden regimes

Hidden regimes

Hidden-regime inference: recovering a switch nobody logged. Import from ergodic.latent.

One observed series; a hidden two-state process (call it the switch); a mechanism that changes how the series behaves while the switch is on. This pillar asks three questions and keeps them apart on purpose. recover() answers WHEN the switch was on and WHICH named mechanism block changed: it returns the posterior state path and a delta interval per block (level, spread, rhythm, amplitude, in the built-in vocabularies). certify() answers WHETHER a regime exists at all, vocabulary-free: it fits the likelihood's own one-state null, simulates control worlds from it, and places the observed two-state improvement in that measured null distribution. calibrate() measures a configured detector's false-fire rate on the same kind of control worlds, with a Wilson interval on the rate. Conflating the three is the standard failure of switching models: a fitted path always explains the data it was fitted to, so "which block moved" must never double as "a regime exists". consensus() lays several likelihoods' verdicts on one series side by side, absences (a vocabulary that cannot express a block) kept distinct from non-fires.

The model underneath is one hidden Markov chain whose path is summed out exactly by the forward algorithm, so NUTS samples only continuous parameters. What varies is the emission: a RegimeLikelihood says how the regime shows up in the data, declares its delta blocks with shrinkage priors, and everything else (chain priors, marginalization, sampling, decoding, label alignment, verdicts) is shared engine. Two likelihoods ship: spectral_whittle(), the fast flagship that reads level, spread, rhythm, and amplitude off framed local spectra, and switching_ar(), the autoregressive second implementation that proves the contract generalizes and is documented blind to amplitude (a sinusoid's AR coefficients do not depend on its amplitude).

Label switching is handled by alignment, never by ordering constraints: decoded state trajectories are correlated against an iterated reference and flipped draws have their deltas sign-corrected, so the posterior stays untruncated and a delta near zero is not pushed off it. Verdicts are equal-tailed credible intervals at a familywise-corrected level, gated by the likelihood's own sanity checks (the spectral wave gate). The measured sentinels behind those choices ride in the module docstrings: the spectral model fires on zero of twenty control worlds at the 99% familywise level and recovers a planted noise tripling as a log-spread delta near 1.1, and the switching AR fires six of six on level, spread, and rhythm arms while its amplitude blindness stands.

pymc is imported inside the functions that sample, so the bare install never loads it; pip install 'ergodic[latent]' adds it, and jax plus numpyro get the fast sampler path.

Also exported from ergodic: CalibrationReport, CertifyResult, Consensus, LatentResult, MarkovChain, RegimeLikelihood, calibrate, certify, consensus, latent_model, recover.

CalibrationReport#

class CalibrationReport(
    n: int,
    fires: int,
    rate: float,
    ci_lo: float,
    ci_hi: float,
    worlds: tuple[CalibrationWorld, ...],
) -> None

The measured false-fire rate of one configured detector on control worlds.

n control worlds were simulated from the likelihood's fitted one-state null on the observed series and each ran through the full recover pipeline; a world counts as fired when any block's verdict came back true. rate is fires / n and the interval is Wilson's, which stays honest at 0 of n and n of n where the textbook interval collapses to a point.

A detector that fires on control worlds is not slightly wrong, it is uninformative: every fire on real data then carries the same doubt. The spectral likelihood's own history says why this report exists: an intermediate version fired on 45% of control worlds and looked excellent on effect worlds while doing it.

CalibrationReport.summary#

def CalibrationReport.summary() -> str

One sentence, the number a decision reads.

CalibrationReport.to_frame#

def CalibrationReport.to_frame() -> pd.DataFrame

Per-world rows: world, seed, fired, fired_blocks.

CertifyResult#

class CertifyResult(observed: float, null: Array, p_value: float) -> None

Whether a regime exists at all: the observed improvement in its own null.

observed is twice the log-likelihood gain of the two-state model over the likelihood's fitted one-state null on the real series; null holds the same statistic recomputed on control worlds simulated from that null, one entry per world. p_value is the standard simulation p, (1 + #{null >= observed}) / (1 + B): with B control worlds the finest attainable evidence is 1 / (1 + B), and reporting zero would claim more resolution than the simulation bought.

This is deliberately a different object from verdicts(). A delta interval says which block moved given a regime; this record says whether the regime story beats a world that never had one, with the null distribution measured rather than assumed (the BIC alternative prices an extra state at a fixed ln n per parameter, and the reference lab measured that price wrong in both directions on these worlds).

CertifyResult.n_null#

property CertifyResult.n_null: int

How many control worlds the null distribution holds.

CertifyResult.to_frame#

def CertifyResult.to_frame() -> pd.DataFrame

The whole test as a frame: one row per statistic, observed first.

Columns statistic and source ("observed" or "null"), so the null histogram and the observed line plot straight off it.

Consensus#

class Consensus(
    target: str,
    models: tuple[str, ...],
    verdicts: Mapping[str, tuple[bool | None, ...]],
    level: float,
) -> None

Several likelihoods' verdicts on one series, laid side by side.

Each model speaks its own vocabulary, so a block can be missing from a model rather than quiet in it: verdicts maps each block to one entry per model, True/False where the model declares the block and None where it does not (the switching AR has no amplitude block, and that absence is structural, not a non-fire). agreed collects the blocks every declaring model settles the same way; contested the blocks where declaring models disagree, which is where a cross-model rule would adjudicate.

One rule is not available: the amp/spread rule (a spectral amplitude fire against a volatility model's spread silence separates a louder wave from a noisier floor) needs a vol_garch likelihood that ergodic does not have, so amp_vs_spread() raises by name.

Consensus.blocks#

property Consensus.blocks: tuple[str, ...]

Every block any model declares, in first-seen order.

Consensus.agreed#

property Consensus.agreed: dict[str, bool]

Blocks where every model that declares them says the same thing.

Consensus.contested#

property Consensus.contested: tuple[str, ...]

Blocks where declaring models disagree; a cross-model rule's domain.

Consensus.amp_vs_spread#

def Consensus.amp_vs_spread() -> Any

Reserved: the amplitude-versus-spread cross-model rule.

A spectral amplitude fire and a volatility model's spread silence together separate "the wave got louder" from "the noise floor rose", which no single built-in likelihood can do alone. The rule needs a vol_garch likelihood, which is not available, so this raises.

Consensus.to_frame#

def Consensus.to_frame() -> pd.DataFrame

One row per block: a verdict column per model, then the consensus.

Model verdicts are True/False/None (None: not in that model's vocabulary); consensus is "fires"/"quiet" when the declaring models agree and "contested" when they do not.

DeltaBlock#

class DeltaBlock(
    name: str,
    site: str,
    transform: Callable[[Array], Array] | None = None,
) -> None

One named mechanism block: a verdict-bearing site of the emission model.

name is the vocabulary word ("level", "spread"), site the full site name whose posterior draws measure the state-1-minus-state-0 difference. transform reduces a vector-valued site to the scalar the verdict reads; it receives the sign-aligned draws of the whole decode subsample at once (shape (S, ...)), because a sensible reduction (the dominant component of a lag-coefficient delta, say) is a property of the posterior, and returns shape (S,).

EmissionSpec#

class EmissionSpec(
    loglik: Any,
    deltas: tuple[DeltaBlock, ...],
    gates: tuple[Gate, ...] = (),
) -> None

What emissions() returns to the engine.

loglik is the emission log-density tensor: (T, K), or (T, K, K) for an emission that conditions on the previous state as well. deltas are the declared mechanism blocks, gates the vetoes.

Gate#

class Gate(
    name: str,
    blocks: tuple[str, ...],
    check: Callable[[Mapping[str, Array]], tuple[float, bool]],
) -> None

A likelihood-declared sanity condition that can veto named verdicts.

A delta interval can exclude zero for a parameter the data never pinned down, so a likelihood may require its own evidence before a verdict stands: the spectral model's wave gate demands that the fitted line actually clears the noise continuum before a frequency or amplitude verdict is allowed to fire. check receives the decode subsample's posterior draws under bare parameter names and returns the gate statistic and whether it passed. A gate should be label-invariant (built from max or symmetric functions over the two states), since it runs on raw, unaligned draws.

GateOutcome#

class GateOutcome(
    name: str,
    blocks: tuple[str, ...],
    value: float,
    passed: bool,
) -> None

One gate's measured statistic and whether it passed.

blocks names the delta blocks this gate vetoes when it fails; the statistic and its meaning are the declaring likelihood's own (the spectral wave gate's is the 3rd percentile of the line's log-power excess over the continuum under its peak).

LatentModel#

class LatentModel(
    model: RegimeLikelihood,
    options: Mapping[str, Any] = <factory>,
) -> None

A validated recover() configuration, held as a value.

Built by latent_model(); everything recover() accepts except the data and the target lives in options. run applies it. The calibration harness deep-copies one of these per control world, which is why the configuration is a frozen value rather than a closure.

LatentModel.run#

def LatentModel.run(data: TimeSeriesData, target: str) -> LatentResult

Recover with this configuration.

LatentResult#

class LatentResult(
    target: str,
    gammas: Array,
    flips: BoolArray,
    delta_draws: Mapping[str, Array],
    gates: tuple[GateOutcome, ...],
    tick_multiplier: int,
    offset: int,
    chain: MarkovChain,
    likelihood: Any,
    n_draws: int,
    divergences: int | None = None,
    rhat_max: float | None = None,
    rhat_max_unaligned: float | None = None,
) -> None

Posterior state path and mechanism deltas for one recovered switch.

gammas is (S, T_eff): the decode subsample's aligned posterior probability of state 1 per emission step. delta_draws maps each declared block to its (S,) sign-aligned draws. tick_multiplier and offset place the emission steps back on the original series' clock: step t covers ticks offset + t * tick_multiplier onward.

Arrays are copied in and frozen, so the result owns its draws.

LatentResult.n_decoded#

property LatentResult.n_decoded: int

How many posterior draws were decoded and aligned.

LatentResult.n_steps#

property LatentResult.n_steps: int

Emission steps in the decoded path.

LatentResult.blocks#

property LatentResult.blocks: tuple[str, ...]

The declared delta blocks, in declaration order.

LatentResult.gamma#

def LatentResult.gamma() -> Array

Posterior state marginals, (T_eff, K): the mean over decoded draws.

LatentResult.path#

def LatentResult.path() -> IntArray

The hard path at emission resolution: state 1 where its marginal clears one half.

LatentResult.tick_path#

def LatentResult.tick_path() -> IntArray

The hard path on the original clock, one entry per tick.

Each emission step's state is repeated over the ticks it covers. The array starts at tick offset of the original series and ends at offset + n_steps * tick_multiplier; ticks the preparation consumed or truncated carry no state and are simply absent.

LatentResult.deltas#

def LatentResult.deltas(
    level: float = 0.94,
) -> dict[str, tuple[float, float, float]]

Equal-tailed (lo, mean, hi) per block at level, no correction.

This is the descriptive read of the posterior; verdicts() is the corrected, gated decision.

LatentResult.verdicts#

def LatentResult.verdicts(level: float = 0.99) -> dict[str, bool]

Which blocks moved: zero-exclusion at familywise level, gates applied.

The per-block interval runs at the Bonferroni share 1 - (1 - level) / n_blocks, so with the default level and four blocks each interval is equal-tailed at 99.75%, and the chance any inert block fires stays at or under one percent. A block a failed gate names is vetoed whatever its interval says.

LatentResult.to_frame#

def LatentResult.to_frame() -> pd.DataFrame

The decoded path as a frame: one row per emission step.

Columns: step, tick (the step's first tick on the original clock), gamma (posterior probability of state 1), state (the hard path).

MarkovChain#

class MarkovChain(stay_alpha: float = 60.0, stay_beta: float = 1.0) -> None

The prior on the hidden chain: stay ~ Beta(stay_alpha, stay_beta) per state.

The stay probability is the whole personality of a two-state chain, and the prior on it is where dwell knowledge enters. The default, Beta(60, 1), concentrates near one: a regime with a mean dwell around 240 ticks stays put with probability about 0.996 per tick, and a prior that allows rapid flip-flopping invites the sampler to explain noise with a busy path. A likelihood whose emission step spans many ticks supplies its own default at its own resolution through default_chain().

The initial distribution is always Dirichlet(1, 1); there is nothing worth configuring in where a long chain happens to start.

PreparedData#

class PreparedData(
    obs: Mapping[str, Array],
    T_eff: int,
    tick_multiplier: int = 1,
    offset: int = 0,
) -> None

What prepare() hands the engine: numpy only.

obs holds the named arrays the emission reads (frame means and periodograms, or a lag matrix), already standardized. T_eff is the number of emission steps, which is also the length of the decoded path. tick_multiplier says how many ticks of the original series one emission step covers (96 for a spectral frame, 1 for an AR step), and offset how many ticks the preparation consumed at the start (the lag order for an autoregression). Together they place the decoded path back on the original clock.

RegimeLikelihood#

class RegimeLikelihood(k: int = 2) -> None

The unit of extension: how a hidden regime shows up in one series.

Subclass as a frozen dataclass and implement four methods. prepare turns the raw series into the numpy tensors the emission reads. emissions runs inside the active pm.Model(), registers every prior site as f"{name}_{param}", and returns the emission log-density tensor with the declared delta blocks and gates. emissions_np is the numpy mirror of the same density for one posterior draw, which is what the decode's forward-backward runs on; it has to be the same math, because a decode under a different density than the fit answers a different question. default_chain states the transition prior at this likelihood's own step resolution.

null_loglik and simulate_null are the one-state null model that certify() and calibrate() read: the log-likelihood of the best regime-free fit on the prepared series, and a fresh control world drawn from that same fit. The two must describe the same fitted null, or the p-values built on them answer a question nobody asked.

k is fixed at two states: one switch, on or off. More states change the delta vocabulary (which pair differs?) and are a design of their own, so anything else is refused here rather than half-supported.

RegimeLikelihood.default_chain#

def RegimeLikelihood.default_chain() -> MarkovChain

The transition prior at this likelihood's step resolution.

RegimeLikelihood.prepare#

def RegimeLikelihood.prepare(y: Array) -> PreparedData

Numpy-only: standardize, frame, or lag-slice y into emission inputs.

RegimeLikelihood.emissions#

def RegimeLikelihood.emissions(
    name: str,
    data: PreparedData,
) -> EmissionSpec

Register priors inside the active pm.Model() and return the emission spec.

RegimeLikelihood.emissions_np#

def RegimeLikelihood.emissions_np(
    data: PreparedData,
    draw: Mapping[str, Array],
) -> Array

The numpy mirror of emissions() for one draw: (T, K) log-densities.

RegimeLikelihood.null_loglik#

def RegimeLikelihood.null_loglik(data: PreparedData) -> float

The fitted one-state null's log-likelihood on the prepared series.

This is the regime-free half of certify()'s improvement statistic, so it must be conditioned exactly as the two-state marginal log-likelihood is (the autoregression conditions both on the first p ticks), or the two are not comparable.

RegimeLikelihood.simulate_null#

def RegimeLikelihood.simulate_null(
    data: PreparedData,
    rng: np.random.Generator,
) -> Array

A fresh control world drawn from the same fitted one-state null.

Returned in the prepared series' own (standardized) space, at the original series length, so it can go straight back through prepare(). certify() builds its null distribution from these worlds and calibrate() measures the detector's fire rate on them.

calibrate#

def calibrate(
    spec: LatentModel,
    data: TimeSeriesData,
    target: str,
    *,
    n_controls: int = 20,
    seed: int = 0,
) -> CalibrationReport

Is the configured detector honest: measured fire-rate on control worlds.

Fits the likelihood's one-state null on the observed series, simulates n_controls regime-free worlds from it, runs the full recover() pipeline (the spec's own budgets and options) on each, and reports how often any verdict fired, with a Wilson interval on the rate. World i runs on a deep copy of the spec with sampler seed seed + i, so the worlds are independent and the whole report reproduces from its seed.

This is the harness that turned the spectral likelihood from a detector that fired on 45% of control worlds into one that fired on 0 of 20 (the ENBW deflation and the per-frame detrend were the fixes, and this measurement is how they were found). Run it once per configuration you intend to trust, on a series from the domain you intend to trust it on: the control worlds inherit the observed series' own memory and scale, so the rate is about your detector on your data, not a generic promise.

certify#

def certify(
    data: TimeSeriesData,
    target: str,
    *,
    model: RegimeLikelihood,
    chain: MarkovChain | None = None,
    n_null: int = 99,
    seed: int = 0,
    refit_draws: int = 200,
    refit_warmup: int = 200,
) -> CertifyResult

WHETHER a regime exists at all: a simulated-null likelihood-ratio test.

The statistic is twice the log-likelihood gain of the two-state model over the likelihood's own fitted one-state null: the two-state marginal log-likelihood (hidden path summed out by the forward algorithm) at a cheap point estimate of the parameters, minus model.null_loglik. The null distribution is measured, not assumed: n_null control worlds are drawn from the fitted one-state model (model.simulate_null) and each is pushed through the identical computation. The p-value is (1 + #{null >= observed}) / (1 + n_null), so with the default 99 worlds the finest attainable evidence is 0.01, and reporting less would claim resolution the simulation did not buy.

The cheap refit is a short single-chain posterior sample (refit_draws after refit_warmup, seed fixed across worlds), with the statistic evaluated at the posterior-mean parameters. The alternative, a pm.find_MAP point, was measured on one world of each likelihood and rejected: on the spectral model the Laplace shrinkage pulls the MAP's deltas to zero and the "improvement" came out negative on a genuine level world (-7.3 against the short sample's +25.8) while also running about seven times slower (30 s against 4 s per world, the optimizer recompiling per world where the short sample amortizes); on the AR model both were adequate (91.2 against 90.2 on the same world) and only there was MAP faster. The calibration argument for the choice is that the observed series and every null world go through exactly the same procedure, budget and seed included, so under the null the observed statistic is exchangeable with the simulated ones and the p-value is honest by construction; a fancier fit for the observed side alone would inflate the observed statistic relative to its null and manufacture significance.

Budget honestly stated: one short fit per world, so the default n_null=99 costs about 100 short fits (minutes on the spectral flagship, longer on the tick-resolution AR). Trim n_null for a coarser but faster answer; the p-value's floor rises with it.

Returns a CertifyResult with the observed statistic, the null array, and the p-value.

consensus#

def consensus(
    results: Sequence[LatentResult],
    *,
    level: float = 0.99,
) -> Consensus

Several likelihoods' verdicts on one series, block by block.

Each model speaks its own vocabulary, so the table distinguishes three things a flat comparison would blur: a block both models fired, a block they disagree on, and a block one model cannot express at all (the switching AR's structural amplitude blindness is an absence, not a non-fire, and is shown as None). Verdicts are recomputed from each result's delta draws at the shared level, familywise-corrected and gated exactly as verdicts() reports them.

One cross-model rule is not available: adjudicating a spectral amplitude fire against a volatility model's spread verdict (a louder wave against a noisier floor) needs a vol_garch likelihood that ergodic does not have, so amp_vs_spread() raises by name.

latent_model#

def latent_model(*, model: RegimeLikelihood, **options: Any) -> LatentModel

Configuration now, data later: recover()'s keywords as a frozen spec.

Validated eagerly. A keyword recover() does not take is refused here, and so is a chain of the wrong type, because the consumer of this spec is calibrate(), where the first use of a bad option would otherwise be twenty NUTS runs in.

recover#

def recover(
    data: TimeSeriesData,
    target: str,
    *,
    model: RegimeLikelihood,
    chain: MarkovChain | None = None,
    graph: Any | None = None,
    draws: int = 1000,
    warmup: int = 1000,
    chains: int = 4,
    seed: int = 0,
    decode_draws: int = 150,
) -> LatentResult

Recover a hidden two-state regime from one observed series.

Fits the switching model model describes (its emission, its declared mechanism deltas, its gates) with the hidden path marginalized out, then decodes the path per posterior draw and aligns the labels. The result carries the posterior state marginals, the hard path, the sign-aligned delta draws, and the familywise-corrected verdicts.

Parameters data: A TimeSeriesData; rows are already in time order. Panel and hierarchical inputs are refused by name. target: The continuous column to read the regime from. model: A RegimeLikelihood: spectral_whittle() or switching_ar(), or your own subclass. chain: The transition prior as a MarkovChain. None means the likelihood's own default at its emission resolution, which is almost always what you want: dwell knowledge stated per tick is wrong by a factor of the frame length for a framed likelihood. graph: Not supported: one regime drives one series; pass nothing. draws, warmup, chains, seed: The sampler budget and seed. The defaults are the measured reference configuration (4 chains of 1000 after 1000 warmup). decode_draws: How many posterior draws to decode and align (evenly spaced over the posterior). The path marginals stabilize well before the full posterior is spent, and each decoded draw costs a forward-backward pass.

Returns LatentResult Frozen: state marginals, hard path, delta draws, gates, verdicts, diagnostics.

semi_markov#

def semi_markov(*, min_dwell: int) -> MarkovChain

Reserved: a chain with an explicit dwell-time distribution.

A geometric dwell is what a Markov chain can express, and the worlds this pillar was measured on have negative-binomial dwells, which the Beta stay prior only approximates. An explicit semi-Markov chain (dwell floor, or a full dwell distribution) needs its own forward pass and is not available, so this raises.