“Is that to say we are against Free Trade? No, we are for Free Trade, because by Free Trade all economical laws, with their most astounding contradictions, will act upon a larger scale, upon the territory of the whole earth; and because from the uniting of all these contradictions in a single group, where they will stand face to face, will result the struggle which will itself eventuate in the emancipation of the proletariat.”
Karl Heinrich Marx · Marx-Engels Collected Works, Vol. VI, p. 290
You can also find this library at CRAN and download it directly from R and RStudio.
Also, we recommend viewing the mind map summary at the end of the article to better understand the relationship between the functions of the package.
LISTEN TO THIS POST AS A PODCAST:
ESCUCHA ESTA PUBLICACIÓN COMO PODCAST
R Library Review
Meet gdpar
General Dynamic Parameter Models via Reference Anchoring
In the fleeting calculus of a two-second decision—overtaking a car on a narrow road—the human brain performs a remarkable statistical trick. It does not build a model of the approaching driver from scratch. Instead, it retrieves a baseline: the average driver, representing typical reaction times and modal aggression. In a split second, it reads the specific signals of the actual driver—relative speed, vehicle type, micro-movements—and estimates how this specific driver deviates from the baseline. The decision to overtake emerges from that synthesis.
This cognitive recipe—population reference + individual deviation—is the philosophical bedrock of the R package gdpar (General Dynamic Parameter models via Reference Anchoring) by José Mauricio Gómez Julián. The package takes this intuition, formalizes it as a rigorous statistical decomposition, proves the conditions under which it is mathematically identifiable, ships a Stan-based Bayesian engine to estimate it, and layers on causal inference, geometry-adaptive sampling, and dependence-robust inference.
❦
The Anatomy of Deviation
Every layer of gdpar is an elaboration of a single, elegant equation. For each observation $i$ with covariates $x_i$:
Read it as: the parameter of individual $i$ equals a population reference, plus a deviation that is itself a function of the individual’s covariates and of the reference itself.
That final clause is where the architecture pivots from classical statistics. The deviation $\Delta$ does not merely depend on who you are (your covariates $x_i$); it depends on what the reference is. If you transplant the model to a new population, the deviation function behaves differently because $\theta_{\text{ref}}$ is one of its arguments. This structural dependence is the defining feature of “reference anchoring.” It distinguishes gdpar from random-effects or varying-coefficient models, where the deviation is structurally separate from the reference.
So, what is the shape of $\Delta$? The package singles out a specific functional form called the Additive–Multiplicative–Modulated (AMM) decomposition:
Three mechanisms, cleanly separated and independently interpretable:
$a(x)$ — A pure additive shift. Think of this as a traditional fixed-effect driven by covariates.
$b(x)\odot\theta_{\text{ref}}$ — A covariate-dependent scaling of the reference (using the Hadamard/elementwise product). This is where “the deviation depends on the reference” enters multiplicatively.
$W(\theta_{\text{ref}})\,x$ — Covariates are mixed through a matrix $W$ that is, itself, tuned by the reference. This is the explicit, structural reference-dependent channel.
Standard models drop out as special cases. Set $\Delta \equiv 0$ and you have fixed-effects regression. Set $W \equiv 0$ and you have a hierarchical model with multiplicative interaction. Set $b \equiv 0$ and you have a varying-coefficient model. The AMM is the smallest natural family that contains all three and elevates the reference to an active argument of the deviation.
❦
The Three Estimation Engines
gdpar defines three complementary engines for estimating $\Delta$. Crucially, only one is executable in the current release—a deliberate choice to promise a mathematical scope that exceeds the executable surface, and to say so honestly.
Path
Engine
Representation
Status
Path 1
Hierarchical Bayesian (Stan)
Parametric AMM
✅ Operational
Path 2
Varying-coefficient (splines)
Smooth $\beta(z)$
🚧 Conceptual
Path 3
Hypernetwork / Neural Net
Net generates $\theta_i$
🚧 Conceptual
Paths 2 and 3 are documented to “reference grade”—full asymptotic theory (contraction rates, Bernstein–von Mises) is developed in the Wiki—but they abort with gdpar_unsupported_feature_error if invoked. Path 1 places priors on every component ($\theta_{\text{ref}}, a, b, W$) and samples the joint posterior with HMC, yielding native, full-posterior uncertainty.
A Tale of Two Posteriors: EB vs. FB
Within Path 1, gdpar offers two inferential regimes. Full Bayes (FB) via gdpar() samples the joint posterior, remaining most faithful to the cognitive analogy. Empirical Bayes (EB) via gdpar_eb() estimates the hyperparameters by maximizing a marginal likelihood via a Laplace approximation, then samples the remaining parameters conditionally.
The EB vs FB Comparator
Rather than forcing a choice, gdpar treats them as parallel routes. It ships a dedicated comparator, gdpar_compare_eb_fb(), which quantifies agreement on $\theta_{\text{ref}}$ and the reduced parameter vector $\xi$. The Wiki develops the theory to first-class depth: EB and FB lower-level posteriors agree asymptotically (Theorem 7A), while EB intervals under-cover by $O(n^{-1})$ (Proposition 7B). If you have ever wondered if EB is “good enough” for your data, gdpar lets you answer that empirically.
Distributional Regression: Every Parameter is a Slot
gdpar is not constrained to modeling the mean. A probability distribution has multiple parameters—location, scale, shape, tail index, zero-inflation probability—and each one can carry its own AMM decomposition. The package indexes these by $k = 1, \dots, K$:
$$ \theta_i^{(k)} = \theta_{\text{ref}}^{(k)} + \Delta^{(k)}(x_i, \theta_{\text{ref}}^{(k)}), \qquad k = 1, \dots, K $$
The built-in roster covers Gaussian, Poisson, negative binomial, Bernoulli, Beta, Gamma, Student-$t$, Tweedie, ZIP, ZINB, and hurdle families. Zero-inflated and hurdle models receive an especially elegant treatment: both the zero-inflation probability $\pi_i$ and the count parameter $\theta_i$ are anchored to their respective references—a dual deviation design.
❦
The Causal Bridge
Because the AMM form produces individual parameters, individual treatment effects emerge naturally. gdpar_causal_bridge() implements a T-learner: fit the anchored model separately under treatment and control, then read the conditional average treatment effect (CATE) at $x_i$ as the difference of the anchored individual predictions:
A second layer, gdpar_compare_meta_learners(), benchmarks the AMM-based learner against external meta-learners via pluggable adapters: grf::causal_forest on the R side and EconML’s CausalForestDML on the Python side (via reticulate). The framework’s causal claims are benchmarked, not asserted.
Mechanics & Clockwork
Several engineering decisions elevate gdpar from a theoretical exercise to a serious computational environment:
Stan Code Generator: Composes programs from canonical pieces—AMM blocks for $p=1$ and $p \geq 1$, EB marginal/conditional blocks, distributional-$K$ blocks—selected by the resolved $(K, p, \text{family}, W, \text{parametrization}, \text{group})$. The $W$ basis supports B-splines with Stan-side Cox–de Boor evaluation, ensuring differentiability inside HMC.
Identifiability Pre-flight: Before any sampling, gdpar_check_identifiability() runs a Gram-matrix check (Proposition 1C), a per-coordinate cross-component check (C4-bis) for $p > 1$, and a per-group anti-aliasing check (C7). If your design is non-identifiable, you find out before the sampler burns your CPU, accompanied by a structured gdpar_identifiability_error naming the dependent directions.
Data-Driven Reparametrization: Treats the parametrization of $b(x) \odot \theta_{\text{ref}}$ as a pre-fit decision. A short pilot computes an information ratio, dispatching to CP, NCP, or—gdpar‘s root-cause resolution—a linear reparametrization that samples the product $\theta_{\text{ref}} \cdot b$ directly, sidestepping bilinear funnels altogether.
Opt-in Power Tools
Two advanced capabilities are switched off by default, documented as thoroughly as the core path.
1. Geometry-Adaptive Sampling
Hierarchical AMM posteriors can be geometrically hostile—funnels, near-determinism, heavy tails. The opt-in geometry engine climbs a ladder of Riemannian metrics: Euclidean → Fisher/SoftAbs → sub-Riemannian → relativistic/Finsler. A certifying orchestrator diagnoses the pathology, selects a metric, tunes the integrator, and emits a certificate. If full sampling is certified infeasible, a Laplace fallback provides a plug-in posterior with ELPD on par with mgcv-REML or INLA-Laplace.
2. Dependence-Robust Inference
gdpar does not model temporal or spatial dependence in its point structure; instead, it makes the inference robust to dependence (a working-independence + sandwich-variance stance in the spirit of Liang & Zeger, 1986). You receive diagnostics (Durbin–Watson, Ljung–Box, Moran’s $I$) and robust SEs via block bootstrap—moving or circular blocks in time (with the Politis–White flat-top automatic block length), tiled randomized-origin blocks in space. Point estimates remain pristine; only the uncertainty is made honest.
⚠️ Honest Limitations
The Wiki is admirably forthright about scope. Only Path 1 is executable in 0.1.0. Dependence is not modelled—only the inference is made robust. The package’s mathematical scope exceeds its executable surface by design. Read the “Implementation status” notes carefully before relying on a feature.
❦
TL;DR
gdpar takes one of the most natural ideas in human prediction—predict an individual as a deviation from a population reference, where the deviation itself depends on the reference—and transforms it into a fully specified, identifiability-checked, Stan-powered Bayesian regression framework. It is theoretically rigorous, computationally serious, and unusually honest about what it does and does not yet do. If your work involves individual heterogeneity, distributional regression, or causal effect estimation with principled uncertainty, gdpar demands a careful look.
You can also find this libraryat CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
The Problem Every Quantitative Researcher Knows
Imagine you are staring at a spreadsheet with 50 columns and a few hundred rows of macroeconomic indicators. Somewhere inside this matrix, a handful of variables carry the signal you care about. The rest is noise — or worse, confounding variation that masquerades as signal. You need to answer three questions: Which variables actually matter? What is the latent structure driving them? And how persistent are the components you extract?
If you have worked with panel data for any length of time, you know that these questions are rarely addressed by a single tool. You run a LASSO in one environment, a PCA in another, an ADF test in a third. Each method lives in a different package with different input formats, different assumptions, and different output structures. Stitching the results together is left to you.
SignalY is an R package that was built to solve exactly this fragmentation problem. Developed by José Mauricio Gómez Julián and released under the MIT license, it provides a unified framework for signal extraction from panel data through multivariate time series analysis. Its design rests on three analytical pillars — column selection, series decomposition, and persistence analysis — that can be used independently or chained together through a single orchestrating function.
This post walks through the problem SignalY solves, the methods it implements, and why the combination matters for applied econometrics and data science.
Pillar 1: Which Variables Matter?
The first challenge in any high-dimensional analysis is selection. When you have dozens of potential predictors, you need a principled way to determine which ones carry structural information and which ones are along for the ride.
The Horseshoe Prior
SignalY approaches this problem through Bayesian sparse regression with the Horseshoe prior. The Horseshoe, introduced by Carvalho, Polson, and Scott (2010) and refined for practical variable selection by Piironen and Vehtari (2017), is a global-local shrinkage prior with a distinctive property: it is aggressive around zero (shrinking noise variables strongly toward zero) while maintaining heavy tails (allowing true signals to escape shrinkage). This dual behavior makes it particularly well-suited for sparse problems where you expect only a few variables to matter, but you do not know which ones.
The math behind this is elegant. Each coefficient βⱼ is given a prior with two layers of shrinkage:
A local parameter λⱼ that controls how much each individual coefficient is shrunk.
A global parameter τ that governs overall sparsity.
The result is a shrinkage profile where most coefficients collapse toward zero — the global pull — while a small number of coefficients stand apart, barely affected — the local escape. This is what gives the Horseshoe its name: the prior density looks like the shape of a horseshoe, with a sharp spike at zero and long, flat arms extending outward.
SignalY’s fit_horseshoe() function estimates this model and provides built-in shrinkage profile diagnostics, so you can visually inspect which variables survived shrinkage and by how much.
Four Ways to Select Variables
Fitting a model is one thing; converting the posterior into a concrete variable selection is another. SignalY offers four distinct selection strategies, each with different strengths:
Projection predictive selection (select_by_projection()): This is the most theoretically robust approach. It works by projecting the full posterior onto candidate submodels and selecting the smallest submodel whose predictive distribution is close enough to the full model. The reference is Piironen and Vehtari (2017), and the implementation respects the posterior geometry rather than relying on ad-hoc thresholds.
Credible interval exclusion (select_by_credible_interval()): Selects variables whose posterior credible intervals do not include zero. Intuitive and easy to interpret, though it can be conservative in high-dimensional settings.
Shrinkage-based selection (select_by_shrinkage()): Uses the kappa (shrinkage fraction) parameters to identify variables that escaped shrinkage. This is useful when you want to understand the degree of shrinkage, not just the binary in-or-out question.
Magnitude-based screening (select_by_magnitude()): A straightforward effect-size filter. Useful as a first pass or when you need to combine Bayesian inference with a frequentist-style screening step.
The fact that SignalY provides all four in a coherent pipeline — not as separate, unrelated functions — is the key design decision. You can run all four and cross-validate the results, or choose the one that best matches your inferential philosophy.
Beyond Regression: Factor Discovery
Sometimes the question is not “which of my 50 variables matter?” but rather “what are the few latent factors driving all 50?” SignalY addresses this through two complementary methods:
PCA with block bootstrap (pca_bootstrap()): Standard principal component analysis, but with block bootstrap confidence intervals that account for temporal dependence in time series data. It also includes entropy-based topology analysis, which measures the informational content of each component.
Dynamic Factor Models (estimate_dfm()): Implements the Bai and Ng (2002) information criteria for automatic determination of the number of static factors, then fits a VAR (Vector Autoregression) on the factor dynamics. This captures not just what the latent factors are, but how they evolve over time.
Pillar 2: What Is the Underlying Structure?
Once you know which variables matter (or have constructed a composite signal), the next question is decomposition: what are the trend, cycle, and residual components of your series?
This is where signal processing meets econometrics, and SignalY implements three methodologically distinct approaches, each with its own strengths.
Wavelet Decomposition
filter_wavelet() implements the Maximal Overlap Discrete Wavelet Transform (MODWT) using Daubechies wavelets, following the framework of Percival and Walden (2000).
Unlike a Fourier transform, which decomposes a signal into infinite sinusoids (losing all time information), a wavelet decomposes a signal into localized, finite-length oscillations at different scales. The MODWT variant is particularly useful for time series because it does not decimate the data (no downsampling), meaning the output length matches the input length at every scale.
In practice, the wavelet decomposition separates a series into:
Detail coefficients (D1, D2, D3, …): capturing oscillations at progressively coarser time scales — high-frequency noise in D1-D2, business-cycle frequencies in D3-D4, and longer cycles in higher levels.
Smooth coefficients (S): the low-frequency approximation that captures the trend.
SignalY includes multi-resolution variance analysis, which tells you how much of the total variance is explained at each scale. This is invaluable for understanding whether your series is dominated by high-frequency noise, medium-term fluctuations, or long-run trends.
Empirical Mode Decomposition
filter_emd() implements Empirical Mode Decomposition (Huang et al., 1998), a fundamentally different approach. Where wavelets impose a predetermined basis (Daubechies, Haar, Symmlet, etc.), EMD is fully data-adaptive. It works by iteratively sifting the signal — identifying local extrema, fitting envelopes, and subtracting the mean — until it extracts Intrinsic Mode Functions (IMFs) that satisfy specific oscillatory conditions.
The key advantage of EMD is that it makes no assumptions about stationarity or linearity. The IMFs are defined by the data itself, not by a mathematical basis. This makes EMD particularly powerful for:
Non-stationary signals whose frequency content changes over time.
Non-linear oscillations that cannot be captured by fixed-basis decompositions.
Signals where the “natural” decomposition is not known a priori.
The trade-off is that EMD can be sensitive to end effects and mode mixing, though SignalY’s implementation includes standard mitigations.
HP-GC Bayesian Filter
filter_hpgc() implements the Grant and Chan (2017) unobserved-components Hodrick-Prescott filter, estimated via MCMC. This is a significant upgrade over the traditional HP filter, which requires you to manually set the smoothing parameter λ (the famous λ = 1600 for quarterly data, or λ = 6.25 for annual data, or any of the other arbitrary rules of thumb floating around the literature).
The HP-GC approach formulates the decomposition as a Bayesian unobserved-components model:
A trend component whose second differences are penalized (this is the smoothness prior, equivalent to the HP penalty).
A cycle component modeled as an AR(2) process.
The smoothing parameter λ is estimated from the data via MCMC, not fixed by the user.
This removes one of the most criticized aspects of the classical HP filter — its sensitivity to the arbitrary choice of λ — while preserving its interpretability. The output includes the estimated trend, cycle, and residual, each with full posterior distributions.
filter_all(): Compare All Three
A particularly useful design choice is the filter_all() function, which runs all three decomposition methods on the same series and returns the results in a comparable format. This is not just a convenience function; it is an epistemological statement. No single decomposition method is universally correct. By running all three and comparing, you can identify components that are robust across methods (strong signal) versus components that depend on the specific decomposition assumptions (potentially method artifact).
Pillar 3: How Persistent Is the Signal?
The third question — what is the persistence regime of your series or its components — is critical for downstream modeling. If your extracted trend is a random walk, that has very different implications than if it is a stationary AR process. If your cycle is near-unit-root, standard mean-reversion models will fail.
A Comprehensive Unit Root Battery
test_unit_root() runs four classical tests with complementary null hypotheses:
Test
Null Hypothesis
Key Feature
Augmented Dickey-Fuller (ADF)
Unit root exists
Most widely used; sensitive to lag selection
Phillips-Perron (PP)
Unit root exists
Non-parametric correction for serial correlation
KPSS
Series is stationary
Reversed null; useful as cross-check against ADF
Elliott-Rothenberg-Stock (ERS)
Unit root exists
Point optimal test with higher power near unity
The critical insight is that no single test is definitive. The ADF and PP tests can fail to reject a false unit root (low power near unity). The KPSS test has the opposite null hypothesis, so it can detect stationarity that ADF misses. By running all four and synthesizing the results, SignalY provides a more robust classification than any individual test.
The automated synthesis follows a standard decision logic:
If ADF/PP/ERS reject unit root and KPSS fails to reject stationarity → stationary.
If ADF/PP/ERS fail to reject and KPSS rejects → unit root.
Mixed results → borderline / near-unit-root, flagged for careful interpretation.
This automated synthesis is not a black box; the individual test statistics and p-values are all available for inspection. But the synthesis gives you a quick, defensible classification without manually cross-referencing four separate test outputs.
The Orchestrator: One Call, Full Pipeline
The signal_analysis() function is the centerpiece of SignalY’s design philosophy. A single call can run the complete analysis pipeline:
…and returns a unified result object with print(), summary(), and plot() methods. The plot() method generates interactive plotly dashboards with filter trends, coefficient profiles, PCA loadings, and DFM factor panels.
The formula interface (Y ~ X1 + X2 + X3) makes it feel like a standard R regression call, while the methods argument lets you mix and match analytical layers as needed.
How Well Does It Work?
The Wiki includes recovery benchmarks on synthetic data with known ground truth. These are worth highlighting because they address the most important question: does this actually work?
Task
Method
Recovery Metric
Factor structure (3 latent factors)
PCA / DFM
r > 0.95, exact factor count
Sparse variable selection (5 of 50)
Horseshoe
F1 > 0.85, Precision > 0.90
Logarithmic trend recovery
EMD
r > 0.95 with true trend
Multi-scale cycle extraction
Wavelet (D3+D4)
r > 0.70 with true cycle
Stochastic trend + AR(2) cycle
HP-GC Bayesian
Trend r > 0.90, cycle r > 0.50
Stationarity classification
Unit Root Battery
4/4 correct on synthetic data
A few things stand out:
The Horseshoe achieves over 90% precision in a 5-of-50 sparse selection problem. This means that when it says a variable matters, it is almost always right. The F1 score above 0.85 indicates a good balance between precision and recall.
Factor recovery is near-perfect (r > 0.95), and the DFM correctly identifies the exact number of latent factors.
Wavelet cycle extraction at r > 0.70 and HP-GC cycle extraction at r > 0.50 reflect the inherent difficulty of extracting cyclical components from noisy data. These are realistic numbers, not inflated claims.
Unit root classification achieves 100% accuracy on synthetic data with clear-cut cases. Real-world data is messier, but this validates the synthesis logic.
Who Should Use SignalY?
SignalY is built for three overlapping communities:
Economists and econometricians working with panel or multivariate time series data who need to move from raw data to structural inference — identifying relevant variables, extracting latent factors, decomposing signals, and characterizing persistence — without stitching together five different packages.
Quantitative researchers in finance, macro, or political economy who face high-dimensional predictor sets and need principled Bayesian variable selection rather than stepwise regression or arbitrary LASSO tuning.
Data scientists working on signal processing problems where the signals are non-stationary, non-linear, or embedded in high-dimensional panels, and where the standard Python signal processing toolkit does not provide the statistical rigor needed for publication-quality inference.
Getting Started
Installation is straightforward:
# From GitHub
remotes::install_github("IsadoreNabi/SignalY")
library(SignalY)
# Minimal workflow
data<-data.frame(Y=as.vector(Y), X)
result<-signal_analysis(data=data, y_formula="Y",
methods=c("pca", "wavelet", "unitroot"))
plot(result)
The package is MIT-licensed, actively maintained (current version 1.1.2), and designed to work with standard R data frames.
The Bigger Picture
What makes SignalY interesting is not any single method — the Horseshoe prior, MODWT, EMD, and unit root tests all exist in other packages. The value is in the integration. By placing Bayesian sparse regression, spectral decomposition, and persistence analysis inside a single coherent framework with a unified interface, SignalY enables workflows that are difficult to replicate otherwise:
Run a Horseshoe regression to select variables, then decompose the fitted signal with wavelets, then test the stationarity of the extracted components — all without changing packages, data formats, or mental models.
Compare wavelet, EMD, and HP-GC decompositions of the same series to identify robust components versus method-dependent artifacts.
Use the DFM to discover latent factors, then test each factor’s persistence regime to inform your downstream modeling choices.
In applied econometrics, the quality of your inference depends on the coherence of your pipeline. SignalY makes that coherence a feature rather than a chore.
SignalY is developed by José Mauricio Gómez Julián. The source code, documentation, and wiki are available at github.com/IsadoreNabi/SignalY under the MIT License.
You can also find this library at CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
Imagine you’re monitoring a manufacturing line, and the sensor readings suddenly look… different. The numbers aren’t obviously wrong — they’re within range — but something has shifted. Or picture an economist staring at decades of GDP data, trying to pinpoint exactly when a country’s growth model fundamentally changed. Or a doctor watching a patient’s vitals, waiting for the moment “stable” tips into “critical.”
All of these scenarios share the same underlying question: when did the system switch from one regime to another?
This is the problem of regime change detection (also called changepoint detection), and it’s one of the most practically important — and theoretically rich — problems in time series analysis. A new open-source R package called RegimeChange tackles this problem with unusual breadth, combining classical statistics, Bayesian inference, and even deep learning under a single, coherent interface. In this post, I’ll walk you through what it does, how it works, and why it matters — in plain language, without cutting corners on the technical details.
The Problem: Signal vs. Noise
Every observable system fluctuates. The central question of changepoint detection is deceptively simple: is this fluctuation just noise (random variation within the same regime), or is it signal (evidence that the system has transitioned to a qualitatively different state)?
This is harder than it sounds. There’s an irreducible tension between two types of error:
False alarm (Type I): You detect a change where none exists. In manufacturing, this means shutting down a line for no reason. In medicine, it’s a false positive that triggers unnecessary intervention.
Missed detection (Type II): You fail to detect a real change. In finance, this means holding a position through a regime shift. In epidemiology, it means missing the start of an outbreak.
Every detection algorithm negotiates this trade-off through its threshold: how much evidence do you demand before declaring that the world has changed? Set it too low and you cry wolf; set it too high and you react too late. There is no universal answer — the right threshold depends on the relative cost of each type of error in your specific context.
What Makes RegimeChange Different
Several R packages already address changepoint detection — the venerable changepoint package, wbs, not, ecp, and the Python ruptures library, to name a few. Each has strengths and limitations. RegimeChange distinguishes itself by integrating approaches that are usually siloed:
Frequentist and Bayesian methods together. Most packages pick one paradigm. RegimeChange gives you both, plus deep learning, all callable through the same detect_regimes() function.
Offline and online modes. You can analyze a complete dataset retrospectively (“when did the changes occur?”) or process streaming data in real time (“is a change happening right now?”) using the same library.
Native uncertainty quantification. Every changepoint estimate comes with confidence intervals (via bootstrap for frequentist methods) or full posterior distributions (for Bayesian methods). This isn’t just a point estimate — it’s a statement about how confident you should be.
Robustness to messy data. Real-world data has outliers, heavy tails, and autocorrelation. RegimeChange includes configurable robust estimation and AR(1) pre-whitening — features that, as the benchmarks show, make a dramatic difference on contaminated data.
An optional Julia backend. For large datasets, the package can transparently dispatch to high-performance Julia implementations while keeping the R interface you’re used to.
Let’s look at each of these in more detail.
The Three Method Families
RegimeChange organizes its detection algorithms into three families. Here’s what each one does, explained simply.
Frequentist Methods
These are the classical, well-established workhorses of changepoint detection.
CUSUM (Cumulative Sum) is the granddaddy of them all, dating to E. S. Page’s 1954 paper. The idea is beautifully intuitive: you accumulate deviations from a target (or “baseline”) value over time. If the system is in the same regime, positive and negative deviations cancel out, and the cumulative sum hovers near zero. When a change occurs, deviations start piling up in one direction, and the statistic grows. When it crosses a threshold, you sound the alarm. CUSUM runs in O(n) time and is ideal for detecting a single change in mean, especially in online monitoring scenarios.
PELT (Pruned Exact Linear Time), introduced by Killick, Fearnhead, and Eckley in 2012, is the modern gold standard for detecting multiple changepoints. It uses dynamic programming to find the globally optimal segmentation of your data — the set of changepoints that minimizes a cost function (essentially: how well do the segments fit the data, penalized by the number of segments). The “pruning” part is clever: it discards candidate changepoint locations that provably can’t be optimal, which brings the average-case complexity down from O(n²) to O(n). RegimeChange supports several penalty criteria — BIC, AIC, MBIC, MDL — or you can specify a manual numeric penalty.
Binary Segmentation is a simpler, greedy alternative: find the single best changepoint, split the data there, then recurse on each segment. It’s fast (O(n log n)) but doesn’t guarantee the global optimum. Wild Binary Segmentation (WBS), from Fryzlewicz (2014), improves on this by searching over many random sub-intervals rather than the full dataset, making it much more robust when changepoints are close together.
The package also includes FPOP (Functional Pruning Optimal Partitioning), which maintains piecewise-quadratic cost functions for even better theoretical guarantees; E-Divisive, a nonparametric method using energy statistics that can detect changes in any aspect of the distribution; Kernel-based CPD, which maps data into a reproducing kernel Hilbert space to capture complex distributional shifts; and NOT (Narrowest-Over-Threshold), which excels at precise localization.
Bayesian Methods
Where frequentist methods give you a point estimate and (if you ask) a confidence interval, Bayesian methods give you something richer: a full probability distribution over when the change might have occurred.
BOCPD (Bayesian Online Changepoint Detection), from Adams and MacKay (2007), is the flagship. At each time step, it maintains a posterior distribution over the “run length” — the number of observations since the last changepoint. When a new data point arrives, it updates this distribution: either the run length grows by one (no change), or it resets to zero (change detected). The probability of a reset at each step is governed by a “hazard function,” typically a geometric distribution representing a constant probability of change per time step.
What makes BOCPD powerful is that it naturally produces, at every time point, the probability that a changepoint just occurred. This is a fundamentally different kind of output than a binary “change/no-change” flag — it’s a calibrated measure of confidence that updates with every new observation. RegimeChange implements BOCPD with several conjugate prior models: Normal-Gamma (for unknown mean and variance), Normal with known variance, Gamma-Poisson (for count data), and Normal-Wishart (for multivariate data).
Shiryaev-Roberts is another Bayesian sequential method, based on work by Shiryaev from the 1960s. It accumulates likelihood ratios from all possible past changepoint times and is asymptotically optimal for minimizing detection delay — that is, for detecting changes as quickly as possible once they occur.
Deep Learning Methods
For complex, nonlinear patterns that classical methods struggle with, RegimeChange offers optional deep learning detectors (requiring keras and tensorflow):
Autoencoders are trained to reconstruct “normal” patterns. When the data changes regime, reconstruction error spikes — the model can’t rebuild what it hasn’t seen before.
Temporal Convolutional Networks (TCNs) use dilated causal convolutions to model long-range temporal dependencies, operating either in supervised mode (if you have labeled changepoints for training) or unsupervised mode (predicting the next value and flagging prediction errors).
Transformers bring self-attention to the problem, capturing long-range relationships in the time series.
Contrastive Predictive Coding (CPC) learns representations through self-supervised contrastive learning, detecting changes where learned embeddings shift significantly.
An ensemble mode combines multiple deep learning methods, requiring a minimum number of them to agree before declaring a changepoint — a strategy that trades some sensitivity for robustness.
One Interface to Rule Them All
The genius of RegimeChange is that you don’t need to learn a different API for each method. Everything flows through a single function:
library(RegimeChange)
# Generate data with a changepoint at t=200
set.seed(42)
data<-c(rnorm(200, 0, 1), rnorm(200, 2, 1))
# Detect using PELT (the default offline method)
result<-detect_regimes(data, method="pelt")
print(result)
plot(result)
# Switch to Bayesian online detection
result<-detect_regimes(data, method="bocpd",
prior=normal_gamma(mu0=0, kappa0=1,
alpha0=1, beta0=1))
# Use CUSUM for online monitoring
result<-detect_regimes(data, method="cusum",
mode="online", threshold=5)
You specify what kind of change you’re looking for (type = "mean", "variance", "both", "trend", or "distribution"), how many changes you expect ("single", "multiple", or a specific integer), and the penalty criterion for model complexity. The function returns a regime_result object containing the detected changepoints, their confidence intervals, segment-by-segment statistics (mean, variance, etc.), and — for Bayesian methods — the full posterior distribution.
Online Detection
For real-time monitoring, you create a persistent detector object and feed it data one observation at a time:
detector<-regime_detector(method="bocpd",
prior=normal_gamma(),
threshold=0.5)
for (xindata_stream) {
detector<-update(detector, x)
if (detector$last_result$alarm) {
message("Changepoint detected at time ", detector$last_result$t)
detector<-reset(detector)
}
}
This is the pattern you’d use for industrial monitoring, fraud detection, or epidemiological surveillance — anywhere data arrives sequentially and you need to react in real time.
Robustness: Where RegimeChange Really Shines
Real data is messy. It has outliers, heavy tails, and autocorrelation. Standard changepoint methods, which typically assume independent, normally distributed observations, can fail dramatically under these conditions.
RegimeChange addresses this with a configurable robustness layer in its PELT implementation. When you set robust = TRUE, the package:
Winsorizes the data (clips extreme values beyond a specified quantile) to limit the influence of outliers.
Uses Huber M-estimation or Tukey biweight loss functions instead of standard squared-error loss. These functions grow linearly (Huber) or flatten out entirely (Tukey) for large residuals, so a single outlier can’t dominate the cost calculation.
Employs the Qn scale estimator (from Rousseeuw and Croux, 1993), which is far more efficient than the median absolute deviation (MAD) — 82% efficiency vs. 37% — while maintaining the same 50% breakdown point (meaning up to half the data can be contaminated before the estimator breaks).
You can also set robust = "auto", which analyzes the data’s contamination level — comparing the MAD to the standard deviation and counting outlier proportions — and automatically selects the appropriate robustness level ("none", "mild", "moderate", or "aggressive").
For autocorrelated data, the correct_ar = TRUE option estimates the AR(1) coefficient and applies a pre-whitening transformation (subtracting the lagged value scaled by the autocorrelation coefficient) before running detection. This prevents autocorrelation from inflating false positive rates.
The Benchmark Numbers
The package’s wiki reports extensive benchmarks against established R packages (changepoint, wbs, not, ecp) across 17 scenarios with 50 replications each. The results are striking:
Overall: RegimeChange’s automatic selector achieved the highest mean F1 score (0.804) across all methods tested.
Contaminated data: The robust mode achieved F1 = 0.998 across all contamination scenarios, versus 0.479 for the standard changepoint PELT — a 108% improvement on messy data.
Subtle variance changes (2:1 ratio): RegimeChange scored 0.667 vs. 0.400 for reference packages — a 67% improvement.
Strong autocorrelation (AR coefficient 0.7): RegimeChange scored 0.900 vs. 0.720 — a 25% improvement.
In 16 of 17 scenarios, RegimeChange either won or tied. The single loss was marginal (0.017 in one heavy-tailed scenario).
Uncertainty Quantification: Beyond Point Estimates
A changepoint location of “observation 200” is only half the story. The other half is: how confident are we? RegimeChange takes this seriously, providing two types of uncertainty:
Location uncertainty — how precise is the estimate? For frequentist methods, the package uses block bootstrap resampling: it repeatedly resamples blocks of the data (preserving local dependence structure), reruns detection, and computes 95% confidence intervals from the distribution of bootstrap estimates. For Bayesian methods, uncertainty comes naturally from the posterior distribution over run lengths.
Existence uncertainty — is there even a change at all? BOCPD produces, at each time point, the posterior probability that a changepoint just occurred. This is arguably more informative than a p-value: it’s a direct probabilistic statement that you can use for decision-making.
The Julia Backend
R is wonderful for data analysis, but it can be slow for computationally intensive tasks. RegimeChange ships with an optional Julia backend — a complete reimplementation of the core algorithms (PELT, FPOP, BOCPD, CUSUM, Kernel CPD, WBS, and multivariate PELT) in Julia, a language designed for numerical performance.
The integration is seamless. You initialize Julia once:
init_julia()
julia_available() # Returns TRUE if ready
After that, the package automatically dispatches to Julia for large datasets (n > 1000 by default) while using R for smaller ones. You can benchmark the difference with benchmark_backends(). The Julia implementation includes thoughtful numerical engineering: Welford’s algorithm for numerically stable running mean/variance, Kahan compensated summation for cumulative sums, log-domain computations to avoid underflow in BOCPD, and careful handling of ill-conditioned covariance matrices in the multivariate case.
If Julia isn’t installed, the package simply falls back to the R implementation — no errors, no fuss.
Evaluation and Comparison
RegimeChange includes a comprehensive evaluation framework. If you know the true changepoints (e.g., in simulated data or a well-studied real dataset), you can compute a full battery of metrics:
This gives you precision, recall, and F1 score (with a tolerance window for matching), Hausdorff distance (the worst-case localization error), Rand Index and Adjusted Rand Index (segmentation agreement), and the covering metric (a weighted IoU measure). You can also run compare_methods() to test multiple algorithms on the same data and see the results side by side.
The package ships with several built-in datasets: well_log (geophysical measurements), industrial_sensor data, economic_cycles, and simulated_changepoints — giving you realistic examples to experiment with.
When to Use What
Here’s a practical guide based on the library’s design and benchmark results:
Your Situation
Recommendation
Clean data, clear mean shifts
PELT with default settings
Outliers or heavy tails
PELT with robust = TRUE (or "auto")
Autocorrelated time series
PELT with correct_ar = TRUE
Need real-time detection
BOCPD or CUSUM via regime_detector()
Need probability of change (not just yes/no)
BOCPD
Nonparametric, no distributional assumptions
E-Divisive or Kernel CPD
Complex nonlinear patterns (and you have data to train on)
Deep learning methods
High-dimensional multivariate data
Sparse Projection CPD
Not sure which method to use
Ensemble mode, or let method = "auto" choose
The Bigger Picture
What I find compelling about RegimeChange is its philosophical stance. The wiki frames changepoint detection not just as a technical problem but as an epistemological one: how do we distinguish real change from noise? The likelihood ratio — the ratio of how probable an observation is under the “changed” hypothesis vs. the “unchanged” hypothesis — is literally a quantification of surprise. Each observation casts a vote, evidence accumulates, and the decision threshold represents how much surprise we demand before acting.
This framing connects the technical machinery to the real-world stakes. Whether you’re monitoring a factory, tracking a disease outbreak, or studying climate data, the fundamental structure is the same: the world was in one state, it transitioned to another, and you need to detect that transition as reliably and quickly as possible.
RegimeChange gives you a unified toolkit for that task — one that spans the classical-to-modern spectrum, handles messy data gracefully, tells you not just where the change is but how confident you should be, and scales from your laptop to a Julia-accelerated backend when you need speed.
Getting Started
RegimeChange is available on GitHub and installs like any R package from source:
It requires R ≥ 4.0.0 and depends on ggplot2, rlang, cli, and magrittr. The Julia backend is optional (requires Julia ≥ 1.6 and the JuliaCall R package), as are the deep learning methods (require keras and tensorflow). The package ships with three vignettes — an introduction, an offline detection guide, and a Bayesian methods tutorial — and the wiki is a comprehensive reference covering everything from mathematical foundations to detailed benchmarks.
For researchers and practitioners who need rigorous, flexible, and robust changepoint detection in R, RegimeChange is well worth a close look.
You can also find this libraryat CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
You have a national Consumer Price Index. It is a single number per year. What you actually want is the price index for each sector of the economy — manufacturing, services, agriculture, construction — because those individual paths are what your model, your policy analysis, or your investment thesis really needs. You know how the sectors combine into the whole: the weights are public, or at least knowable. What you do not know is the sectoral numbers themselves. You only ever see their weighted sum.
This is the disaggregation problem. It sounds like bookkeeping. It is actually a quietly profound statistical question, and the R package BayesianDisaggregation — built by José Mauricio Gómez Julián — tackles it with a degree of intellectual honesty that is rare in software. This post is a deep, plain-language tour of what the package does, why it works the way it does, and what it teaches about building statistical tools that tell the truth.
If you want to follow along with the code, the project’s wiki has installation instructions, function references, and worked examples. This post deliberately stays code-free so the ideas come through.
The problem, precisely
Let’s make it concrete. You observe an aggregate index — call it the CPI — over T years. You also have a matrix of weights W: for each year and each of K sectors, W tells you that sector’s share of the total. The weights sum to one within each year. The aggregate is, up to measurement noise, the weighted sum of the latent sectoral indices:
CPI at year t ≈ the weighted sum of the K sectoral indices, using that year’s weights.
The goal: recover the K sectoral indices — call them φ — from the single aggregate series and the known weights.
Here is the catch, and it is the entire heart of the matter. At every single year, the aggregate pins down one linear combination of the K sectors. The remaining K−1 directions are completely unconstrained by the data. With four sectors, you have one equation and four unknowns per year. The system is, in a precise mathematical sense, under-determined.
This is not a numerical annoyance you can engineer away. It is structural. Any method that hands you a single, sharp sectoral path from an aggregate alone is, whether it admits it or not, smuggling in assumptions to fill the gap — and most methods do not tell you how much of the answer is data and how much is assumption.
The question is not can you disaggregate. You always can. The question is: can you do it honestly, showing your work, and carrying the right amount of uncertainty forward?
The wrong way: a cautionary tale
The first version of the package — the 0.1.x line — advertised “MCMC-free Bayesian disaggregation.” The pitch was appealing: no slow Markov chain Monte Carlo sampling, just clean closed-form math. The implementation used a family of deterministic update rules — weighted, multiplicative, Dirichlet, adaptive — over the prior weight matrix.
It did not work. Not in the sense of crashing. In the much more dangerous sense of appearing to work while silently not doing the one thing it claimed to do.
The package’s own author, in a moment of radical honesty that is worth pausing on, audited the 0.1.x family and catalogued six foundational defects, labeled F1 through F6:
F1 — the aggregate never entered the computation. The “posterior” was derived entirely from the prior weight matrix. The actual observed CPI — the one piece of real evidence — was never used. The method was not conditioning on data; it was rearranging priors.
F2 — the Dirichlet concentration cancelled on renormalization. A parameter that was supposed to control how concentrated the sectoral estimates were simply vanished in the algebra when weights were normalized to sum to one.
F3 — the temporal pattern cancelled too. A component meant to encode smoothness over time also disappeared in the renormalization.
F4 — the “efficiency” term was a fixed constant. It looked like a data-dependent quality score; it was actually invariant.
F5 — there were no recovery tests. No one had ever generated synthetic data with a known truth and checked whether the method got it back.
F6 — a correlation helper cheated. It computed both Pearson and Spearman correlation and reported whichever was larger, a form of silent data-snooping.
The most damning of these is F1. A Bayesian method that does not condition on data is not a Bayesian method. It is a deterministic transformation dressed in Bayesian vocabulary. And the worst part is that it looked reasonable — it returned numbers, they moved in plausible directions, and nothing crashed.
The author’s response is, I think, the most important thing about this whole project. Rather than patching the defects one by one — adding the CPI here, fixing the concentration there — the author recognized that the foundational problem (not using the data) cannot be fixed within a deterministic re-weighting framework. The fix is not a patch. The fix is a fundamentally different model: one in which the aggregate enters as genuine evidence.
So the entire 0.1.x family was deleted. Every function — bayesian_disaggregate(), compute_L_from_P(), spread_likelihood(), the four update rules, the grid search, the save function, the cheating correlation helper — all of it, gone. In its place: two Bayesian engines that actually condition on the data. The package version jumped to 0.2.0, the DESCRIPTION was edited to remove any claim of novelty, and the documentation was rewritten to be explicit about what was removed and why.
This is what intellectual honesty in software looks like. It is not common. We should notice it when it happens.
The right idea: the aggregate as evidence
The conceptual move is simple to state and deep in consequence. Instead of treating the relationship between the aggregate and the sectors as a renormalization identity — a bit of algebra you apply after the fact — treat it as an observation density. The aggregate CPI is data. It is evidence about the latent sectors. The model should condition on it.
In the package’s canonical engine, this means: the latent sectoral indices φ evolve over time as a random walk with drift, and the observed CPI is generated from the weighted sum of the sectors with some observation noise. The aggregate is not a constraint imposed after estimation; it is the likelihood. The sectors come out the other side as a posterior distribution — not a single number, but a full cloud of plausible values with credible intervals.
This is the difference between solving an equation and updating beliefs in light of evidence. The first gives you a point. The second gives you a distribution. And as we will see, the distribution is the whole point.
Two engines, one trade-off
The package offers two ways to do the Bayesian inference, and the choice between them is a clean, well-explained trade-off between richness and exactness.
The state-space MCMC engine
This is the canonical, full-featured model. The sectoral indices live in log space, which guarantees they stay positive — a natural constraint for price indices that a model in raw levels could violate. Each sector’s log-index follows a random walk with its own drift and its own innovation scale (the amount of jitter per period).
Two layers of hierarchical structure make this more than K independent random walks:
Partial pooling on the drift. Each sector has its own drift, but the drifts are drawn from a common distribution. This means sectors share information about their average growth rate without being forced to be identical — the classic shrinkage trade-off.
Partial pooling on the innovation scale. Similarly, each sector’s volatility is drawn from a shared distribution. Sectors borrow strength in estimating how jittery they are.
The initial cross-section — the starting levels of the sectors at the first period — is anchored at the aggregate level with an estimable dispersion parameter. This is a subtle but important point. In the old, broken 0.1.x family, the “concentration” parameter was supposed to control how spread out the sectors were, but it cancelled in the algebra and had no effect. In the new model, the dispersion is a genuine parameter that the data and priors can actually estimate. It does not cancel. It does real work.
Finally, the observation: the CPI is modeled as coming from a Student-t distribution centered at the weighted sum of the sectors, with an estimable scale. The Student-t (rather than a Gaussian) makes the model robust to outliers in the aggregate — a heavy-tailed observation can be accommodated without wrecking the fit. If you prefer, you can switch to a plain Gaussian observation.
Because this model is not conjugate — the log transform, the Student-t, and the hierarchical structure break the neat algebra that would allow a closed-form solution — it is fit by Hamiltonian Monte Carlo via Stan (using either the cmdstanr or rstan backend). HMC is the gold standard for this kind of model: it handles the correlated, high-dimensional parameter space efficiently and comes with reliable diagnostics. The package runs four chains by default, checks the R-hat convergence statistic and the number of divergent transitions, and returns posterior draws of every sectoral index at every period.
The closed-form conjugate engine
The second engine is the linear-Gaussian counterpart. The sectoral indices evolve as a random walk in levels (not logs), and the aggregate is observed with Gaussian noise. This model is conjugate, which means its exact posterior can be computed in closed form — no MCMC, no sampling, no convergence diagnostics. The tool is the Kalman filter combined with the Rauch-Tung-Striebel smoother: the filter passes forward through time, updating beliefs about the sectors given each new CPI observation, and the smoother passes backward, refining those beliefs using future information.
If you want joint posterior draws — not just the smoothed means and variances, but actual correlated samples from the full posterior — the package uses the Durbin-Koopman simulation smoother, a elegant technique that produces draws with the correct cross-time and cross-sector covariance structure. These draws are not marginal approximations; they are genuine samples from the joint posterior.
This engine is the “correct realization of the original MCMC-free posterior idea.” The 0.1.x family wanted a closed-form Bayesian answer; the problem was not that closed-form is un-Bayesian (conjugacy is perfectly Bayesian), but that the old method did not use the data. This engine uses the data — the aggregate enters as the observation equation — and it does so in exact, closed form.
The trade-off is explicit and documented. The closed-form engine buys you speed and mathematical exactness. It costs you three things: positivity (levels can drift slightly negative, which logs prevent), robustness (Gaussian observations are sensitive to outliers, which the Student-t handles), and the cross-sector hierarchy (there is no partial pooling in the linear-Gaussian model). The MCMC engine buys you all three of those, at the cost of sampling.
Both engines return the same thing: an array of posterior draws of the sectoral indices, with dimensions [periods × sectors × draws], plus summary tables of medians and 95% credible bands.
The honesty at the center: weak identification
Here is where the package earns its deepest respect. It would be tempting, after building a real Bayesian model that conditions on the data, to declare victory and hand users sharp sectoral estimates. The package does not do this. Instead, it is explicit about a fact that most disaggregation methods gloss over: the sectoral split is only weakly identified.
Remember the under-determination: at each period, the aggregate pins down one linear combination of the K sectors. The remaining K−1 directions are governed by the prior — the cross-sectional dispersion and the temporal smoothness — not by the data. This means the posterior intervals for individual sectors are wide, and they are influenced by the prior. This is not a bug. It is not a limitation to be engineered away. It is the correct representation of what the data can and cannot tell you.
The package’s own recovery tests — which generate synthetic data from the model’s own data-generating process, where the true sectoral paths are known — confirm this directly. The aggregate is recovered essentially perfectly: the correlation between the fitted aggregate and the true aggregate is above 0.95, often essentially 1.0. The aggregate is strongly identified. But the sectoral coverage — the fraction of times the true sectoral path falls within the 95% credible band — is around 0.84, with the bands being deliberately wide. The test asserts that coverage exceeds 0.70, a conservative threshold, because the package refuses to claim sectoral precision that the data cannot deliver.
This is “rigour by layers”: assert tightly what is identified, assert conservatively what is not. It would be easy to tune the prior to produce narrower, more impressive-looking bands. The package deliberately does not.
Why the full posterior matters: propagating uncertainty
If the sectoral estimates are uncertain — wide bands, prior-influenced — then what good are they? This is where the design of the package reveals its purpose. The sectoral indices are not the final product. They are input to a downstream model. In the author’s research program, they feed a nested Ornstein-Uhlenbeck model of price gravitation. But the principle is general: any time a disaggregated estimate feeds a second-stage analysis, the uncertainty in the first stage should flow into the second.
The package handles this by multiple imputation, following Rubin’s rules. Each posterior draw of the sectoral indices is treated as one imputation — one plausible version of the truth. The downstream model is fit once per imputation, and the results are combined. The effect is that the weak per-sector identification — the wide bands, the prior influence — is carried forward into the downstream uncertainty intervals rather than being discarded. You do not plug in a point estimate and pretend it is the truth. You plug in the whole cloud and let the cloud’s shape propagate.
The package’s documentation is candid about a consequence: because disaggregation is under-determined, the random-walk smoother prior dilutes the reversion signal, biasing the downstream reversion speed toward slowness by a modest, quantified amount (roughly 13–26%). Crucially, the direction is conservative — the true gravitation is at least as fast as reported — and the fraction of missing information that is propagated is about 0.4. This is not hidden. It is measured, reported, and flagged as a property of honest under-determination.
A sharp contrast: an ad-hoc method that simply added noise to a point estimate did not produce proper imputations and gave sub-nominal coverage when routed through Rubin’s rules. The coherent posterior from the conjugate engine — the one that actually conditions on the data — did. The mathematics of multiple imputation demands proper posterior draws; garbage in, garbage out.
How it is validated
The package’s validation strategy is worth studying because it embodies a philosophy: test the identified quantity tightly, test the unidentified quantity conservatively, and test the computation itself exactly.
Three layers, kept deliberately separate:
Smoke tests run always, on every check. They confirm that both engines compile, sample, and return the correct [periods × sectors × draws] array structure on synthetic data. They catch breakage.
Recovery tests are gated behind an environment flag because they require actually compiling and sampling the Stan model, which is slow. They generate data from the model’s own data-generating process — the same random walk with drift, the same partial pooling, the same aggregate observation — so the true sectoral paths are known. Then they check: does the aggregate come back essentially perfectly? (Yes, correlation above 0.95.) Do the sectoral bands cover the truth at a reasonable rate? (Yes, above 0.70, honestly wide.) The recovery test is well-posed because the simulator uses the same process as the model. If the model cannot recover its own data, something is wrong with the sampler or the implementation. If it can, you have a meaningful baseline.
Golden tests run always and are the most stringent. They use Stan’s generate_quantities function — which deterministically recomputes derived quantities from frozen parameter draws, with no random number generation involved — and demand a bit-for-bit match against a frozen reference output. This catches any change to the model’s computed quantities: if someone edits the Stan code and the log-likelihood values shift by even one bit, the test fails. The reference fixture is generated by the same code path, isolating the CSV serialization so the comparison is exact.
This is not the “does it run?” school of testing. It is the “does it compute the right thing, and does it compute exactly the same thing tomorrow?” school.
Where it sits among existing methods
The package is careful — almost unusually careful — about situating itself relative to the existing literature. It makes no claim of being the first or only solution to this problem. The documentation uses the phrase “we did not find” rather than “we are the first,” and the DESCRIPTION file was explicitly edited to remove any “novel/original” claim.
The adjacent traditions, and what each misses:
Biproportional balancing (RAS, IPF) iteratively scales a matrix to match new margins. It is deterministic: no posterior, no credible intervals, no treatment of the aggregate as evidence. It is a useful accounting tool, not an inference method.
Temporal disaggregation (Denton, Chow-Lin, Fernández) distributes a low-frequency aggregate to higher frequency using an indicator series. This is a temporal problem — splitting annual into quarterly — not a cross-sectional one. It assumes you already have the sectoral decomposition and just need finer time resolution.
Forecast reconciliation (MinT and related methods) projects inconsistent hierarchical forecasts onto a coherent subspace. It is forecast-centric and linear-algebraic: it corrects forecasts that do not add up, rather than recovering latent components from an aggregate by Bayesian updating.
Compositional or Dirichlet state-space models evolve simplex weights over time. They model how shares move, not how the components themselves are recovered conditioned on their weighted sum.
Each tradition addresses a real problem. None, as far as the package’s author could find, does exactly this: recover latent cross-sectional components from a single aggregate by conditioning on it as a genuine observation density and returning a posterior that can be propagated downstream. The claim is narrow and checkable, not sweeping.
The data pipeline
The package is not just a model; it is a usable tool. It includes hardened readers for the real inputs. A CPI reader pattern-matches on column headers (in English or Spanish — it recognizes “date,” “fecha,” “year,” “año” for the time column and “cpi,” “indice,” “price” for the value column), parses localized number formats (European-style decimals and thousands separators), collapses duplicate years by averaging, and returns a clean, sorted data frame. A weights reader loads a sector-by-year table, normalizes weights to the simplex within each year, and handles missing entries gracefully.
An alignment function intersects the years covered by the CPI and the weights, ensuring both cover the same periods before either engine runs. A convenience wrapper reads both files and runs the disaggregation in one call. And a simulator generates synthetic data from the model’s own data-generating process — the same random walk with drift, the same partial pooling, the same aggregate observation — so that recovery tests, examples, and exploratory analysis are always well-posed.
One data note worth flagging, because it is a common error: the model works in index levels, not rates of change. Feeding a percent-change series (inflation rate) instead of a level series (the CPI itself) is a category error — the aggregate would not be on the same scale as the weighted sum of the sectors. The CPI must be a level series, re-indexed to the same base as whatever the sectors will be compared against.
The bigger lesson
You could read this package as a technical contribution: a Bayesian state-space model for disaggregation with two engines, honest uncertainty, and a propagation contract. That reading is correct but incomplete.
The deeper lesson is about how to build statistical software that tells the truth. The 0.1.x family did not fail by crashing. It failed by producing plausible-looking output that did not depend on the data. That is the most dangerous failure mode in statistics, because there is no error message. The numbers look reasonable. The plots look smooth. Nothing warns you that the entire computation is a rearrangement of priors.
The author caught it — caught it in their own work, which is harder than catching it in someone else’s — by doing the unglamorous thing: generating data with a known truth and checking whether the method recovered it. When it did not, they did not patch. They deleted and rebuilt. And then they documented the deletion, in public, with the defects labeled and explained, so that anyone reading the history would understand not just what changed but why.
The resulting package has a quality that is hard to name but easy to feel when you read the source: every design choice has a reason, every reason is documented, and the documentation is honest about what the method can and cannot do. The aggregate is strongly identified; the sectors are weakly identified; the uncertainty is wide and prior-influenced; and all of that is surfaced, not hidden, because the whole point is to carry that uncertainty forward rather than fake it away.
In a field where it is tempting to claim sharp results from sparse data, this is a quiet act of integrity. The package does not solve the under-determination. Nothing can. It does something better: it honors it, by returning the wide, honest, propagatable posterior that the data actually supports.
The package, its source code, installation instructions, and full function reference are on GitHub, with extended documentation in the wiki. It is MIT-licensed and written in R, with the MCMC engine powered by Stan.
You can also find this libraryat CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
When Market Prices Gravitate: A Bayesian Look at an Old Question in Economics
An old question, asked again — properly
There is a question in economics that is older than most of the academic disciplines that border it. Do market prices — the noisy, day-to-day, here-and-now prices at which goods actually change hands — tend to settle toward some underlying center of gravity? And if they do, how fast, how violently, and through what mechanism?
Classical political economy, from Smith and Ricardo through Marx, thought they do. The idea was that behind the churning surface of market prices there sit “prices of production”: long-run, cost-anchored prices toward which actual prices are pulled, the way a spring pulls a weight back toward its rest position. In the Marxian version, there is one more layer underneath: those prices of production themselves are supposed to gravitate around “values,” the labour embodied in commodities. Whether any of this is true is an empirical question, and for a long time the empirical tools to answer it were not really up to the job.
A small R package called bayesianOU, written by José Mauricio Gómez Julián and hosted on GitHub, takes a serious swing at that question. It is not the first attempt to test price gravitation statistically, but it is one of the most technically careful I have seen, and it is built in a way that is instructive far beyond the Marxian debate that motivates it. What follows is a walkthrough of what the package does, why it is interesting, and — just as importantly — where it honestly admits its own limits.
The tool that makes it possible: the Ornstein-Uhlenbeck process
Strip the economics away for a moment and the statistical core of the package is a workhorse object from physics: the Ornstein-Uhlenbeck (OU) process. Imagine a particle moving in a fluid, attached to a spring. Brownian motion jiggles it randomly; the spring pulls it back toward a fixed point. The further it drifts away, the harder the pull. The result is a wiggly series that never settles but always tends to settle — a mean-reverting random walk.
The OU process is exactly the mathematical object you want when you suspect a variable is noisy but anchored. It has a “speed of reversion” (how hard the spring pulls) and an “equilibrium level” (where the spring’s rest point is). Estimate those, and you can say something quantitative about gravitation: not just “yes, prices come back,” but “they come back with a half-life of about nine years.”
That number — the half-life — is the prize. It is the difference between “market prices eventually settle” (which could mean anything) and “market prices settle on a timescale comparable to a business cycle” (which is a falsifiable, interpretable claim).
What the package actually builds
The package fits, by Bayesian inference, a family of models built on the OU process but considerably richer than the textbook version. There are two first-class models, sharing one inference engine.
The single-level model
The first model asks: do market prices revert toward an equilibrium that is a function of the prices of production, and what does that reversion look like once we let it be nonlinear, volatile, heavy-tailed, and structurally heterogeneous across sectors?
Each of those adjectives is doing real work, and each corresponds to a feature that simpler approaches handle poorly or not at all:
Nonlinear drift. A plain OU process pulls back with a force proportional to the deviation. The package allows a cubic correction, so the restoring force can strengthen super-linearly when prices are far from equilibrium. This matters: real markets may behave gently near the center and violently at the extremes, and a linear model cannot represent that.
Stochastic volatility. Financial data, and economic data generally, go through quiet stretches and turbulent ones. The package does not assume a single noise level; it lets the volatility itself wander over time, following its own mean-reverting process on the log-variance. This is the same idea that powers modern stochastic-volatility models in finance, and it is essential for not fooling yourself about the precision of your estimates.
Heavy tails. Economic shocks are not Gaussian. Crashes, booms, and policy shocks produce outliers that a normal distribution would call essentially impossible. The package uses Student-t innovations and estimates the degrees of freedom from the data, so the model can discover for itself just how fat-tailed the world is.
Hierarchical structure across sectors. An economy has dozens of sectors, and each one presumably has its own reversion speed, its own equilibrium, its own noise. Estimating each sector in isolation throws away the information that they are all part of the same economy. Estimating them all with one set of parameters pretends they are identical. The package takes the middle path — hierarchical, or “partial pooling,” priors — where each sector’s parameters are drawn from a shared distribution whose properties the model also estimates. Sectors borrow strength from one another without being forced into lockstep.
A time-varying coupling. This is the most economically loaded feature. The strength with which market prices track prices of production is allowed to depend on the aggregate profit rate (what the package calls TMG). When the general rate of profit is high, the pull of production prices on market prices may be one thing; when it is low, another. Whether that modulation exists, and in which direction, is a hypothesis the model can test rather than assume.
All of this is estimated jointly, with full Bayesian uncertainty, using Stan’s Hamiltonian Monte Carlo sampler. You do not get a point estimate of the reversion speed; you get a posterior distribution, and from it a credible interval and a probability statement like “there is a 95% chance the half-life is between six and eighteen years.”
The nested cascade
The second model is the more ambitious one, and it is where the package earns its “nested” branding. Instead of market prices reverting to a fixed equilibrium, they revert to a latent production price — a hidden, unobserved series that itself evolves over time according to its own OU process, driven by the general profit rate. And, if you turn on the third level, that latent production price in turn gravitates around an observed “value” index built directly from labour-content accounting.
So the full structure is a cascade: market price → latent production price → value. Each arrow is an OU reversion, each with its own speed, and the speeds are constrained so that the outer (market) layer reverts faster than the inner (production) layer — an economically natural separation of timescales, enforced softly so the data can push back.
The reason this matters is that it converts a slogan — “prices of production gravitate around values” — into a literal statistical hierarchy that can be fit to data and compared against alternatives. The headline empirical result, from a fit to 37 US sectors over 1960–2020, is a value-coupling coefficient essentially equal to one, with the posterior probability of it being positive effectively equal to one. In plain terms: in standardized units, prices of production track labour values almost one-for-one. That is a found result, not an assumed one — the prior on the coupling was centred at zero, deliberately neutral.
The inference engine, and why it is not a footnote
It would be easy to glance at the model description, nod, and move on. But how these quantities are estimated is half of what makes the package serious, and it is worth a paragraph for readers who do not think about MCMC every day.
Bayesian inference works by exploring the space of all parameter values consistent with both the data and the prior, and characterizing that space as a probability distribution. For models this complex — with latent volatility paths, hierarchical structure, and hundreds of parameters — you cannot do that with pencil and paper. You use a Markov chain Monte Carlo sampler, specifically Hamiltonian Monte Carlo, which borrows an idea from physics: give the parameter space a “potential energy” (the log-posterior) and a “kinetic energy” (a randomly chosen momentum), and let the system glide around the posterior like a ball rolling over a landscape.
Stan’s NUTS sampler automates this about as well as it can be automated, and the package uses it with within-chain parallelism (via Stan’s reduce_sum) to handle the fact that the likelihood must be summed over many timepoints and sectors. The diagnostics — R-hat for chain agreement, effective sample size, divergence counts — are surfaced through a validate_ou_fit function, and the package is explicit that you should look at them before believing anything.
Model comparison is done with PSIS-LOO, a clever technique that approximates leave-one-out cross-validation without refitting the model dozens of times, by reweighting the posterior draws using importance sampling. It is the modern standard, and the package is appropriately cautious about it: because the model has a latent volatility state at every observation, plain LOO is known to be optimistic, and the documentation says so plainly.
The honesty that makes it credible
Here is where the package surprised me, and here is why I think it deserves a wider audience than the Marxian-economics niche it lives in.
A naïve reading of the results would be triumphant: the value coupling is one-to-one, the reversion exists, the half-life is about nine years. But the package’s own validation section does something rare. It runs the model against legitimate rivals on genuinely held-out data — a full decade, 2011 to 2020 — and reports, without spin, that a random walk beats the OU model at forecasting, that a no-gravitation restriction ties or beats it, and that the value term adds no detectable predictive density.
That sounds like a refutation. The package argues, carefully, that it is nothing of the sort — and the argument is the most intellectually interesting thing here.
The key move is to distinguish two different questions. One is structural: does a reversion mechanism exist, and how fast is it? The other is predictive: can you forecast next year’s price better than a naïve benchmark? These are related but not identical, and for a slow process they come apart in a specific, predictable way.
If gravitation is real but slow — a half-life of nine years on a dataset whose test window is a decade — then over the forecast horizon the process looks, to first order, like a random walk. The reversion is there, but it is too weak to show up in a one-step or few-step prediction. The random walk, which assumes no reversion, will forecast almost as well, because over short horizons a barely-reverting process and a non-reverting one are nearly indistinguishable. So the random walk winning the forecasting horse race is not evidence against gravitation; it is evidence consistent with gravitation being slow.
This is not special pleading. It is a logical point about what different functionals of a model can and cannot tell you. The structural parameters — estimated from the joint likelihood over the whole panel, borrowing strength across 37 sectors and 61 years — use far more information than any single-series forecast. They can pin down a central tendency that a univariate test cannot. And the package shows, through simulation-based calibration and adversarial negative controls, that the estimation pipeline does not manufacture gravitation when none is present: feed it a true random walk and it reports a half-life of about fifty years; feed it a null value-coupling and the posterior honestly covers zero.
The low-kappa trap, and why it matters to everyone
The package names a difficulty it calls the low-kappa trap, and it is worth understanding because it is a trap that catches far more than Marxian price theory.
Kappa is the reversion speed. As kappa shrinks toward zero, the OU process approaches a pure random walk. The trouble is that there is no bright line separating “slow mean reversion” from “no mean reversion.” It is a continuum, and three distinct problems stack up exactly there:
Algebraically, reversion speed and discrete-time persistence are two sides of the same coin; kappa going to zero is the same as the autocorrelation going to one. There is no internal frontier.
Statistically, the power of a unit-root test — the standard tool for asking “is this a random walk?” — collapses exactly as the truth approaches the random walk boundary. With a finite sample and a half-life comparable to the sample length, the test simply cannot tell. This is a well-known result in econometrics, and it is why decades of “is the real exchange rate stationary?” papers argued past one another.
Numerically, if the reversion speed is parameterized to be strictly positive (as it must be, for the sampler to behave), then “the probability that kappa is greater than zero” is trivially one — it tells you nothing. The informative quantity is the half-life, and the probability that the half-life exceeds some sensible horizon.
The package’s response to the trap is instructive. It does not pretend the trap is not there. It states all three layers explicitly, reports the slow tail honestly (one sector has a non-trivial posterior probability of a half-life beyond forty years), and argues that the joint hierarchical posterior — which pools information across the whole panel — is a more powerful discriminator than any univariate test. That is a defensible position, and it is stated with the caveat attached rather than buried in a footnote.
This is the broader lesson. Anyone working with time series that might be slowly mean-reverting — interest rates, real exchange rates, commodity prices, climate variables, pollutant concentrations — runs into exactly this trap. The package’s framing of it, in three layers, is one of the clearest expositions I have read, and it would travel well into any of those domains.
What I appreciate, and what I would watch for
A few things stand out as genuinely good practice, and they are worth naming because they are rarer than they should be.
The separation of economic and sampler convergence. The package is scrupulous about not confusing two senses of “convergence.” Economic convergence — does the price revert? — is a statement about kappa and the half-life. Sampler convergence — did the MCMC chains agree? — is a statement about R-hat and divergences. These share a word and nothing else, and conflating them is a classic source of muddled reasoning. The documentation keeps them lexically distinct throughout.
Neutral priors on the load-bearing hypotheses. The prior on the profit-rate modulation is centred at zero. The prior on the value coupling is centred at zero. The package does not bake the answer into the question. When the posterior then moves clearly away from zero, that means something.
Out-of-sample integrity by construction. A subtle and common error in time-series work is “leakage”: accidentally letting future information contaminate the training procedure, so that out-of-sample results are secretly in-sample. The package offers a fit_window switch that keeps the two designs genuinely separate, and it computes the common-factor loadings from the training window only. This is the kind of plumbing detail that separates trustworthy work from work that just looks trustworthy.
The negative results are reported. Many packages, and most blog posts about them, would quietly omit the fact that a random walk out-forecasts the model. This one leads with it and then reasons about it. That is how a field accumulates reliable knowledge rather than just encouraging headlines.
What should a careful reader watch for? The half-life estimate of about nine years is, by the package’s own account, probably conservatively slow — a controlled study of the disaggregation step suggests the true figure may be closer to seven or eight. The cubic nonlinearity is a minor refinement on this data (its coefficient sits near its prior). The Student-t degrees of freedom and the stochastic-volatility scale are weakly identified when both are present, a known tension the documentation flags but does not resolve. And the headline value-coupling result, while striking, is measured on standardized levels that share a cost-price component by construction; the package defends this with a “wedge” argument — subtracting the shared component and testing the residual — but a sceptical reader should follow that argument itself rather than take it on trust.
None of these caveats undermine the project. They are the project. A statistical framework that cannot articulate its own soft spots is not a framework you should believe.
Why it is worth your time
You do not need to be a Marxian economist, or any kind of economist, to get something out of this package. If you work with time series that exhibit slow, noisy reversion toward a moving target — and a great deal of the physical and social world does — the modelling ideas here are directly portable: the nonlinear OU drift, the stochastic volatility, the hierarchical pooling across groups, the careful separation of structural estimation from forecasting, and the three-layer diagnosis of the low-reversion trap.
And if you are interested in the classical question of whether prices gravitate toward values, this is about as good a statistical treatment as you will find: modern machinery, honest reporting, and a willingness to let the data argue back against the theory that motivated the exercise in the first place.
The repository, the full mathematical specification, the validation blocks, and a frank discussion of every methodological decision live at github.com/IsadoreNabi/bayesianOU, with the wiki carrying the complete technical detail. Read the methodology notes before you quote a number; that is what they are there for.
You can also find this library at CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
There is a habit so embedded in quantitative work that most practitioners never think to question it. You have a time series — quarterly GDP, an EEG channel, a temperature record — and at some point you fit a smooth curve through it, interpolate a missing value, or estimate a “long-run trend.” All of these moves rest on a single, seldom-checked assumption: that the data form one continuous whole, that a single smooth function can legitimately pass through every point.
But what if they don’t? What if your series is, structurally, two or three disjoint pieces glued together by the calendar — pieces between which no continuous function can travel? In that case, the spline you just fitted is not an approximation of reality; it is a mathematical fiction painted over a fracture.
topologyR is an R package that lets you check this before you model. It takes a numeric time series, converts it into a graph, converts that graph into a topological space, and then asks the one question that determines whether global continuous methods are even valid: is this space one connected piece, or several?
It sounds abstract. It is abstract — but the consequence is concrete. The package is the work of José Mauricio Gómez Julián, and it is open-source, with a GitHub repository, a detailed Wiki, and a companion research paper archived on Zenodo. What follows is a tour of what the package does, why it matters, and where it fits in the broader landscape of topological data analysis.
The Hidden Assumption
Think about what happens when you impute a missing value in a time series using a cubic spline. The spline assumes that the points on either side of the gap belong to the same continuous process — that the missing value lies somewhere along a smooth bridge between them. If the series has actually undergone a structural break, a regime change, or a discontinuity between those points, the spline will happily produce a number, and that number will be wrong in a way no confidence interval can capture.
This is not a niche problem. It appears in econometrics (trend estimation across business cycles), in neuroscience (coherence across brain-state transitions), in climatology (warming trends across regime shifts). The methodological error is always the same: assuming continuity without first verifying that continuity is mathematically possible.
topologyR’s contribution is to make that verification explicit, parameter-free, and exact.
From Numbers to Shapes: The Pipeline in Three Steps
The package’s workflow has an elegant, almost architectural logic. You feed it a series of numbers; it returns a topological verdict. Between input and output, three transformations occur.
Step 1: The Series Becomes a Graph
The first move is borrowed from network science: the visibility graph. Imagine your time series plotted as a mountain range — each observation is a peak or a valley at a given time. Two points are connected by an edge if you could stand on one and see the other, with no taller peak blocking the line of sight.
topologyR implements two flavours. The Horizontal Visibility Graph (HVG) connects two points if every point between them is strictly lower than the shorter of the two — a horizontal line of sight. It runs in linear time and captures the skeleton of the series’ ups and downs. The Natural Visibility Graph (NVG) is more generous: it connects two points if every intermediate point lies below the straight line joining them, regardless of the heights of the endpoints. It is denser, richer, and runs in O(n log n) expected time. The NVG always contains the HVG as a subgraph.
Both are parameter-free. There is no threshold to tune, no bandwidth to select, no ε to agonise over. The graph is determined entirely by the data’s own geometry. This matters enormously: it eliminates the single largest source of arbitrariness in the entire pipeline.
Step 2: The Graph Becomes a Topology
Here is where topologyR departs from ordinary network analysis. A graph tells you who is adjacent to whom. A topology tells you something deeper: what the neighbourhood structure of the entire space looks like — which collections of points form coherent open regions, and how those regions combine.
The construction follows a method introduced by Nada, El Atik, and Atef in 2018. For each vertex v in the graph, you form its closed neighbourhood — the vertex itself plus all its direct neighbours. This family of closed neighbourhoods serves as a subbase. You then close it under finite intersections to obtain a base, and close the base under arbitrary unions to obtain the full topology.
If those words feel heavy, think of it this way: the subbase is a rough draft of “who belongs with whom.” Intersecting neighbourhoods refines the draft — “the points that both neighbourhoods agree on.” Taking unions completes the picture — “every region that can be assembled from these building blocks.” The result is a genuine topological space, complete with open sets satisfying the standard axioms, sitting on top of your time series like a scaffolding you didn’t know was there.
Step 3: The Topology Reveals Its Connectivity
Now comes the decisive question. A topological space is connected if it cannot be split into two non-empty open pieces — if there is no clean fracture running through it. For finite spaces, there is a beautiful theorem, due to McCord (1966) and Stong (1966), that makes this check exact and tractable. The specialization preorder orders the points by how their neighbourhoods nest inside one another, and the connected components of the resulting structure are precisely the topological connected components.
The crucial practical point: this computation works directly on the base — the refined building blocks — without ever needing to enumerate the full topology (which can be exponentially large). It runs in polynomial time, and the components it returns are exact, not approximate.
The Decision Rule
Everything so far converges on a single, actionable verdict. topologyR hands you a connectedness decision, and that decision has a direct methodological consequence:
If the induced topology is connected, then your data are consistent with a single continuous process. Global continuous methods — splines, kriging, polynomial interpolation, moving-average imputation, kernel methods — are mathematically supported. You may proceed.
If the induced topology is disconnected, then no single continuous function can cover the entire series. Global continuous methods are invalid by construction. You must segment the series along the connected components the package identifies, and model each piece independently — with regime-switching models, component-wise imputation, or finite mixtures.
This is the package’s value proposition: a reproducible, topology-first workflow that decides, before you touch a model, whether global continuity is a justified assumption or a silent error.
Global versus Local
A subtlety worth flagging: the rule depends on what you are trying to learn. Global properties — a secular trend, a Hurst exponent, total neural synchronisation, a centennial warming signal — depend on relationships among all points and require topological connectivity to be valid. Local properties — instantaneous volatility in a small window, point-to-point rates of change, low-order autocorrelation — are defined on restricted neighbourhoods and remain valid within each connected component, regardless of whether the whole series is one piece or several. The package gives you the component structure to make that distinction operational.
Time Has an Arrow: Directed Topologies and Irreversibility
So far, the construction has treated the visibility graph as undirected — time flows, but the edges don’t care which way. That discards information. Time series are inherently directional: time runs from past to future, and many real systems are irreversible — they behave differently forwards and backwards. Economic expansions creep upward over years; recessions collapse in quarters. Neurons fire and recover on different timescales. The undirected graph cannot see this asymmetry.
topologyR’s directed mode fixes this. With directed = TRUE, each visibility edge is oriented from the earlier time point to the later one, producing a directed acyclic graph (a DAG) in which the time index is a natural topological order. From this directed graph, the package extracts two neighbourhood structures: the forward neighbourhood (who can I see ahead of me?) and the backward neighbourhood (who behind me can see me?).
Applying the Nada construction to each yields two topologies: a forward topology τ⁺ and a backward topology τ⁻. The pair (X, τ⁺, τ⁻) forms what Kelly (1963) called a bitopological space — a set equipped with two topologies rather than one. The divergence between them is a direct, topological measurement of temporal irreversibility.
Irreversibility Indices
In a perfectly reversible process — symmetric dynamics, no privileged direction — the two topologies coincide: τ⁺ ≅ τ⁻. They have the same number of connected components, the same base size, the same connectivity. In an irreversible process, they pull apart.
topologyR quantifies this with several indices. The component irreversibility measures the normalised difference in the number of connected components between the forward and backward topologies: zero means symmetric, one means maximally asymmetric. The base irreversibility does the same for the sizes of the topological bases. The asymmetry direction — the signed difference in component counts — tells you which way the arrow points: a positive value means the forward topology is more connected (fewer components) than the backward one.
That last point has a concrete physical interpretation. Consider a time series with gradual expansions and abrupt contractions — the classic shape of a business cycle, where GDP creeps up over years and drops in a quarter. During a gradual rise, forward visibility is relatively unobstructed: looking ahead from a point on the upslope, you can see far. After an abrupt drop, backward visibility is blocked: looking back from the trough, the cliff face hides earlier points. This asymmetry means the forward topology should be more connected than the backward topology — fewer forward components, more backward fragmentation. The package predicts, and the data confirm, a positive asymmetry direction for such series.
The Alexandrov Layer and the Resolution Hierarchy
There is a third topology lurking in the directed graph, and it is older than the Nada construction by several decades. The Alexandrov topology τ_A, introduced by Alexandrov in 1937, is the topology whose open sets are the upsets of the reachability relation — the sets that, once you enter them, contain everything reachable downstream. For each vertex, its minimal open set is the collection of all vertices reachable from it via directed paths.
topologyR computes this efficiently: a reverse-order bitset propagation that processes vertices from last to first, OR-ing reachability sets together in O(nm/64) time, reusing the same high-performance bitset infrastructure as the Nada engine.
The relationship between the Alexandrov and Nada topologies is precise and informative: τ_A is always a subset of the forward Nada topology. The Alexandrov base captures pure order structure — “who can reach whom” — while the Nada intersection closure generates additional sets that are not upsets, catching finer-grained structure. The difference in base sizes, |B_Nada| − |B_A|, tells you exactly how much extra topological information the Nada pipeline extracts beyond the raw ordering. A large gap means the closure operations are doing real work; a small gap means the order structure already tells most of the story.
Under the Hood: Performance Without Compromise
Topological enumeration is, in the worst case, exponential — the number of open sets can in principle double with every additional element. This is an inherent mathematical fact, not a software limitation. But topologyR is engineered so that the decision you actually need — connectedness — never requires that enumeration.
The connectivity computation works on the base alone, via the specialization preorder, in O(n² · ⌈B/64⌉) time. The C++ backend (via Rcpp) represents every subset as a packed array of 64-bit words, so set operations reduce to machine-level bitwise instructions. A compile-time template dispatch selects single-word operations for series up to 64 points, two-word for up to 128, three-word for up to 192 — zero loop overhead, branch-free. Beyond that, a runtime fallback handles arbitrary sample sizes, and OpenMP parallelisation is available where the build supports it.
The practical upshot: you can run the connectivity decision on series with thousands of points without ever touching the exponential regime. Safety limits (max_base_sets, max_open_sets) cap the intersection and union closures with informative termination flags, so if a computation does hit resource limits, you know exactly where and why — and the connectivity result remains valid as long as the base closure completes.
A Real-World Test: Reading the Business Cycle
The paper accompanying the package applies the framework to quarterly U.S. real GDP growth from 1992 to 2024 — 129 observations spanning over three decades. The bitopological analysis recovers a positive asymmetry direction: the forward topology is more connected than the backward one, exactly as predicted for a series with gradual expansions and abrupt contractions.
The undirected topology partitions the series into six connected components, each corresponding to a distinct macroeconomic regime. Strikingly, the COVID-19 contraction and its rebound — the deepest and fastest swing in the sample — are classified as a single topological episode: one connected component spanning the collapse and recovery, reflecting the fact that the visibility structure treats the V-shaped episode as one structural unit rather than two separate events.
This is the kind of insight the package is designed to produce: not a forecast, not a parameter estimate, but a structural classification that tells you where the legitimate boundaries in your data lie — and, critically, whether a global model is appropriate at all.
Where It Sits: Complementary, Not Competing
It is important to be clear about what topologyR is not. It is not a general-purpose topological data analysis (TDA) engine. Packages like GUDHI, Ripser, TDAstats, and scikit-TDA compute persistent homology — multi-scale features across all dimensions, capturing higher-order structures (loops, voids) via Betti numbers β₁, β₂ and their persistence across scales. That is a richer and harder enterprise.
topologyR has a narrower and more focused aim: it zeroes in on β₀ — connectedness — for one-dimensional series, using graph-induced topologies, and it turns that single invariant into an actionable decision rule for method selection. Think of it as a pre-model governance tool: a rigorous gatekeeper that runs before you choose your modelling strategy, telling you whether the continuity assumptions your favourite methods require are actually justified by the data’s structure.
The two approaches are complementary. For early-warning detection, precursor signals, or multi-channel structure, persistent homology is the right tool. For the binary question “can I legitimately fit a global continuous model to this series?”, topologyR gives a direct, interpretable, and mathematically exact answer. A natural hybrid workflow uses topologyR as a pre-test and persistent homology for deeper multi-scale analysis.
Honest Limitations
No tool is universal, and topologyR is transparent about its boundaries:
Graph choice matters. HVG and NVG produce different graphs, and therefore potentially different topologies. The NVG, being denser, tends to produce fewer connected components. The package encourages comparing both and interpreting the difference — the gap itself is diagnostic.
Sampling and noise. Sparse sampling can mimic disconnection; minor overlaps can mimic connection. The connectedness verdict should be treated as prima facie evidence, not absolute truth — especially near the boundary.
β₀ only. The approach focuses on connectedness. It will not capture loops, voids, or higher-order patterns that persistent homology can detect. If your question is about cycles or multi-scale structure rather than fragmentation, you need the heavier machinery.
Enumeration is exponential; connectivity is not. This is handled honestly: the connectivity decision is polynomial and scalable; full topology enumeration (needed for pairwise connectedness in the bitopological sense) is capped by safety limits with transparent reporting.
The Bigger Picture
What makes topologyR more than a clever technical exercise is its epistemological stance. It transforms a step that is normally a tacit habit — assuming continuity — into an explicit, testable, mathematical procedure. In doing so, it removes arbitrariness from one of the most consequential decisions in applied quantitative work: the choice between global and segmented methods.
The package’s central theorem — that the Nada construction extends to directed graphs and yields a bitopological space whose asymmetry quantifies irreversibility — is formalised in Lean 4 against Mathlib, so the mathematical foundation is not merely asserted but machine-checked. The implementation is CRAN-compliant, passes R CMD check --as-cran cleanly, and ships with 68 unit tests covering visibility graphs, topology generation, connectivity, directed topology, Alexandrov topology, and bitopological analysis.
For anyone who works with time series and has ever fitted a spline, run a kriging model, or estimated a trend — which is to say, for most of applied quantitative science — topologyR offers something rare: a way to check, before you model, whether the smoothness you are about to assume is a property of your data or a story you are telling yourself.
The package is authored by José Mauricio Gómez Julián and released under the MIT licence. It requires R ≥ 4.0.0 with Rcpp and ggplot2. The companion paper, “Bitopological Spaces from Directed Graphs: Extending the Nada Construction to Capture Temporal Irreversibility,” develops the full mathematical theory, including the central theorem, the Alexandrov sublayer, specialization preorder, pairwise connectedness, polynomial-time algorithms, and the Lean 4 formalisation.
If you use topologyR in your research, please cite the repository release.
