“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.
What If David Ricardo Was Wrong? A New Econometric Challenge to Comparative Advantage
Based on: Gómez Julián, J. M. (2025). “Teorías del comercio internacional versus resultados de los tratados comerciales.” Revista Cubana de Economía Internacional, 12(1), 36–57. Read the original paper (Spanish)
Most people who have taken an introductory economics course have encountered a deceptively simple idea: countries should specialise in what they do relatively best, even if another country is better at producing everything. This is the doctrine of comparative advantage, largely attributed to David Ricardo’s early-nineteenth-century work on trade between England and Portugal. It has become one of the most cited justifications for free trade and for the architecture of modern trade agreements.
A 2025 paper by Juan Manuel Gómez Julián, published in the Revista Cubana de Economía Internacional, asks a provocative question: does the actual data from trade agreements support comparative advantage — or does it point back to the older, simpler idea of absolute advantage? His answer, reached through a combination of historical analysis, mathematical reasoning, and modern econometric modelling, is likely to unsettle a good deal of conventional trade-policy thinking.
The Two Competing Ideas, in Plain Language
Before diving into the paper’s contribution, it helps to be absolutely clear about what is at stake. Imagine two countries:
Country A can produce both wheat and steel more efficiently (faster, cheaper, with fewer resources) than Country B.
Country B is less efficient at producing both goods.
Absolute advantage (Adam Smith, 1776) says: Country A is simply better at both. Country B has no obvious reason to compete head-to-head, and trade between them will be shaped by the sheer gap in productive capability.
Comparative advantage (David Ricardo, 1817) says: hold on — even though Country A is better at both, it is proportionally better at steel than at wheat. Country B, while worse at everything, is relatively less terrible at wheat. So if Country A focuses on steel and Country B focuses on wheat, and they trade, both end up better off. Absolute superiority does not matter; what matters is the ratio of efficiencies within each country.
This idea is elegant. It is also, as Gómez Julián argues, surprisingly fragile when tested against real-world data.
What the Paper Actually Does
Gómez Julián approaches the question from three complementary angles, which gives the paper unusual methodological breadth.
1. Mathematical Generalisation
First, he examines how well each theory holds up when you push it mathematically — that is, when you ask whether the logic remains sound under more general and realistic assumptions than the original two-country, two-good textbook models. Comparative advantage, he finds, depends on a narrow set of assumptions (identical technologies in certain respects, constant costs, no transport costs, full employment) that tend to collapse when the model is made more realistic. Absolute advantage, by contrast, remains coherent under a wider range of conditions.
2. Historical Context
Second, the paper traces the intellectual history. Ricardo developed comparative advantage in a world where the nature of production was fundamentally different from today’s globalised, technology-intensive economy. The author argues that the theory was a product of its time — useful as a thought experiment, but not a reliable guide for modern trade policy, especially when the technological gap between trading partners is vast.
3. Econometric Evidence
This is where the paper makes its most distinctive contribution. Gómez Julián uses two families of statistical models to test which theory better explains the actual outcomes of trade agreements:
Computable General Equilibrium (CGE) models — large-scale simulation models that attempt to represent the entire economy, sector by sector, and then simulate what happens when a trade agreement changes tariffs, quotas, or market access. These are widely used by institutions like the World Bank and the WTO.
Objective Bayesian Generalised Linear Models (GLMs) — a modern statistical approach that uses Bayesian inference (updating beliefs with data) with minimal subjective assumptions (“objective” priors). This allows the researcher to let the data speak more freely, without imposing strong preconceptions about what the answer “should” be.
The combined results point in the same direction: trade outcomes between countries with significant technological asymmetries are better explained by absolute advantage than by comparative advantage.
What Does This Mean in Practice?
The practical implications are significant, and they run against the grain of mainstream trade-policy advice for the past several decades.
If comparative advantage is the correct lens, then free trade between any two countries — rich or poor, technologically advanced or not — is mutually beneficial almost by definition. The policy prescription is straightforward: liberalise, sign agreements, reduce barriers.
But if absolute advantage is the better model, then the structure of the agreement matters enormously. A trade deal between a highly industrialised country and a predominantly agricultural one is not inherently win-win. It may lock the less-developed country into low-value-added exports while flooding its markets with manufactured goods that undercut local industry. The technological and wage asymmetries between the signatories become the central concern, not an afterthought.
In other words, Gómez Julián’s findings suggest that trade agreements should be designed with deliberate attention to the power imbalances and productive capacities of the parties involved — not simply signed on the assumption that any trade is good trade.
Why This Matters Beyond Economics
If you are a political scientist, a policy analyst, or simply someone who follows geopolitics, this debate is far from academic. Trade agreements are among the most consequential instruments of foreign policy and domestic economic strategy. They shape industrial policy, labour markets, migration patterns, and even geopolitical alliances.
The question of whether a trade deal is “fair” or “beneficial” depends on which economic theory you use to evaluate it. If the dominant theory is wrong — or at least incomplete — then decades of trade policy advice may have systematically underestimated the risks of liberalisation between unequal partners.
This does not mean protectionism is the answer. But it does mean that the terms of engagement matter. A trade agreement that accounts for technological gaps, includes provisions for technology transfer, and builds in adjustment mechanisms is a very different instrument from one that simply eliminates tariffs between unequal economies and calls it a day.
Beyond the Phillips Curve — A Marxist Reinterpretation of Inflation
Political EconomyJuly 2025 · 8 min read
Beyond the Phillips Curve
A new study argues that inflation isn’t just about too much money chasing too few goods — it’s about how the capitalist class converts technological advantage into permanent profit.
Based on: Gómez Julián, J. M. (2025). “Más allá de la curva de Phillips: inflación, cambio tecnológico y acumulación de plusvalía en las economías capitalistas contemporáneas.” Realidad Económica, IADE.
Most of us were taught a tidy story: when unemployment falls, inflation rises, and vice versa. This trade-off — called the Phillips Curve — has anchored central bank policy for decades. But what if that story is not just incomplete, but fundamentally misleading?
A recent paper published in Realidad Económica by José Mauricio Gómez Julián argues exactly that. Using over fifty years of U.S. data (1968–2021), the study finds no significant long-run relationship between inflation and unemployment. Instead, it identifies a surprising positive link between technological change and inflation — and uses that finding to build a Marxist reinterpretation of what inflation actually does inside a capitalist economy.
It’s a paper that challenges both mainstream economics and the popular imagination. Let me walk you through it.
The Phillips Curve: A Love Story with Complications
In 1958, New Zealand economist A.W. Phillips noticed an elegant regularity in British data: wages tended to rise faster when unemployment was low. Later economists generalized this into a policy menu: want less unemployment? Accept a bit more inflation. Want to tame prices? Brace for a recession.
This trade-off became gospel in the 1960s. Central bankers thought they could fine-tune the economy like a thermostat — dial inflation up or down by adjusting demand. But the 1970s shattered that confidence. The U.S. experienced stagflation: high inflation and high unemployment at the same time, something the Phillips Curve said shouldn’t happen.
Since then, economists have debated whether the Phillips Curve is dead, dormant, or merely sleeping. Gómez Julián sides with a more radical verdict: the long-run Phillips Curve doesn’t just flatten — it was never there to begin with.
What does “long run” mean here?
Mainstream economists already accept that the long-run Phillips Curve is vertical (meaning no permanent trade-off). But Gómez Julián goes further: he finds that even in shorter cycles, the supposed inverse relationship is statistically fragile — easily dissolved once you account for other variables, especially technological change.
The Data, the Tools, and What They Found
The study uses three complementary statistical approaches — each chosen for a reason:
Bayesian Correlations
Unlike classical statistics, which gives you a yes-or-no answer (“significant at 5%”), Bayesian analysis lets you say something more nuanced: “Given the data, here is the probability that this relationship is positive, negative, or nonexistent.” Applied to U.S. inflation and unemployment, the Bayesian results show no consistent inverse relationship. The data simply doesn’t support the Phillips Curve story with any confidence.
Granger Causality
This is a standard econometric test that asks: does knowing today’s unemployment help you predict tomorrow’s inflation (or vice versa)? If the Phillips Curve were real, the answer should be yes. Gómez Julián finds that the answer is generally no. Unemployment does not Granger-cause inflation in the U.S. data. What does show predictive power? Research and development spending.
Error Correction Models (ECM)
These models examine whether variables that drift apart over time eventually pull back together — like two dancers who briefly separate but remain on the same floor. The ECM results confirm that inflation and unemployment do not share a stable long-run equilibrium. They are, statistically speaking, dancing to different music.
· · ·
The Surprising Link: Technology Drives Inflation
Here is the paper’s most provocative finding: R&D expenditure and inflation move together positively. When firms invest more in technology, inflation tends to rise — not fall, as you might expect from a productivity-enhancement standpoint.
Why would better technology lead to higher prices? To answer this, Gómez Julián turns to Marx — specifically, to the distinction between two types of surplus value.
Capitalist innovates (new machinery, process)
→
Extraordinary surplus value (temporary advantage)
→
Rivals adopt technology
→
Inflation absorbs the gap
→
Relative surplus value (permanent for the class)
Fig. 1 — The mechanism proposed by Gómez Julián, simplified.
Two Kinds of Surplus Value: A Quick Primer
If you’re not steeped in Marxist theory, don’t worry — the distinction is intuitive.
Absolute surplus value is what a capitalist gets by making workers work longer or harder for the same pay. It’s the old-fashioned squeeze. Relative surplus value, by contrast, comes from making production cheaper — through technology, efficiency, better organization — so that the value of labor-power (i.e., the cost of maintaining a worker) falls, even if wages don’t.
Now imagine a single firm introduces a breakthrough technology. It can produce goods faster and cheaper than its competitors. For a while, it earns extraordinary surplus value — a premium profit that exists only because it’s ahead of the pack. But here’s the catch: once competitors adopt the same technology, that advantage vanishes. The extraordinary surplus value disappears.
Gómez Julián’s argument is that inflation is the mechanism through which this temporary advantage gets converted into a permanent one. How? As the innovating firm’s higher productivity drives down unit costs, prices don’t fall proportionally — instead, the general price level adjusts upward. The gap between the old cost structure and the new one gets absorbed by inflation. What was a one-time windfall for the innovator becomes a structural shift in profitability for the entire capitalist class.
Inflation, in this reading, is not a policy error or a monetary accident. It is a functional mechanism of capitalist accumulation — one that converts technological advantage into lasting class-wide profit.
What This Means for Policy
If the paper is right, the implications are significant:
For central bankers: If inflation isn’t primarily a monetary phenomenon — if it’s rooted in the structural dynamics of production and profit — then raising interest rates to fight inflation is treating the symptom, not the disease. You might cool the economy, but you’re not addressing the engine that generates inflation in the first place.
For mainstream economists: The Phillips Curve may be less a stable empirical law and more a historical coincidence — a relationship that appeared to hold in a particular postwar context and has been propped up by theoretical convenience ever since. The paper adds to a growing body of evidence that the curve has become unreliable as a guide to policy.
For non-economists: This paper reframes inflation as a political question, not just a technical one. If inflation systematically benefits capital at the expense of labor — by preserving the gains of innovation for the capitalist class while workers’ purchasing power erodes — then debates about inflation are, at their core, debates about distribution and power.
A note of caution
The paper uses R&D spending as a proxy for technological change. This is standard in the literature, but it’s not a direct measure of innovation. R&D spending can reflect many things — tax incentives, defense contracts, speculative bubbles in tech. The correlation Gómez Julián finds is suggestive and theoretically grounded, but it warrants further investigation with additional proxies and across different economies.
· · ·
A Challenge to Orthodoxy
What makes this paper worth reading — whether you agree with it or not — is that it does something many economists avoid: it takes a heterodox theoretical framework seriously and tests it empirically. This isn’t armchair Marxism. It’s Bayesian statistics, Granger causality, and error correction models applied to five decades of data. The methodology is conventional; the interpretation is not.
The mainstream view treats inflation as essentially a monetary phenomenon — too much money, not enough stuff. Milton Friedman’s famous dictum that “inflation is always and everywhere a monetary phenomenon” still echoes through central banks worldwide. Gómez Julián doesn’t deny that money supply matters. But he argues it’s not the whole story — and may not even be the most important part.
In his framework, the relationship between technology, surplus value, and prices is structural. It doesn’t depend on whether a central bank is dovish or hawkish. It’s embedded in the logic of capitalist production itself.
So, Is the Phillips Curve Dead?
Probably not entirely. There are short-run contexts where demand pressures do push prices up, and the Phillips Curve captures something real about those moments. But the paper pushes us to ask harder questions: What determines the baseline around which those fluctuations occur? Why has inflation behaved the way it has over half a century, regardless of the unemployment rate?
Gómez Julián offers a provocative answer: inflation is the economy’s way of metabolizing technological progress into profit. It’s not a bug in the system. It’s a feature.
Whether you find that convincing depends, in part, on your theoretical priors. But the data doesn’t lie about what it doesn’t show: a reliable Phillips Curve. And that, at minimum, should give everyone — mainstream, heterodox, and curious layperson alike — something to think about.
A new research uses probability theory — and sixty years of U.S. economic data — to test one of the most consequential (and most overlooked) assumptions in political economy.
The Assumption Hiding in Plain Sight
If you’ve ever read a Marxist analysis of how profits equalize across industries, you’ve probably encountered something called the average rate of profit. The idea is straightforward: competition between capitalists drives different rates of profit in different sectors toward a common, system-wide average. This is one of the pillars of Marx’s theory of value in Capital, Volume III.
But there’s a quieter assumption underneath this one — so quiet that most discussions never mention it explicitly. To arrive at a uniform profit rate, Marx first assumes a uniform rate of surplus value across all productive sectors. In plainer terms: he assumes that the degree to which workers are exploited — the ratio of unpaid labor to paid labor — is roughly the same everywhere, whether you work in steel manufacturing, food processing, or textiles.
Adam Smith proposed this idea before Marx. Smith argued that if one job were obviously more exploitative (in the sense of yielding far more unpaid surplus per dollar of wages paid), workers and capital would flow toward or away from it until the differences vanished. Marx adopted this observation and, as scholar Jonathan Cogliano notes, elevated it to “the status of a central economic law” within his framework.
Yet the assumption has been challenged from multiple directions — Marxist and non-Marxist alike. Is it actually justified? Or is it a convenient simplification that distorts our understanding of how capitalism works?
José Mauricio Gómez Julián, of the Universidad Latina de Costa Rica, decided to approach the question from an unexpected angle: probability theory. His paper, published in Ciencia Económica (2022), asks whether the mathematical law that would need to hold for this assumption to be valid actually does hold — and then checks the answer against six decades of real-world data from the United States.
The Mathematical Backbone: The Law of Large Numbers
If you’ve taken any statistics course, you’ve likely met the Law of Large Numbers (LLN). It tells us that as you observe more and more instances of something random — coin flips, dice rolls, stock returns — the average of those observations settles down toward the true expected value.
There are two versions:
The Weak Law (WLLN): With enough observations, the sample average is probably close to the expected value.
The Strong Law (SLLN): With enough observations, the sample average is almost certainly equal to the expected value — a much stronger guarantee.
Gómez Julián’s insight is this: if you think of each productive sector of the economy as a random variable representing that sector’s rate of surplus value, then the LLN tells you what happens to the average across sectors as the number of sectors grows large. In mathematical language:
Strong Law: The probability that the average surplus-value rate across sectors equals the global expected value, in the limit, is exactly 1.
Weak Law: The probability that the average deviates from the global expected value by more than any tiny amount shrinks to zero.
If either version holds, you get the result Marx needs: across a sufficiently large number of sectors, the rates of surplus value converge to a common value — uniformity, or at least a powerful tendency toward it.
The Catch: Independence and Identical Distributions
Here’s where things get interesting — and where the classical LLN hits a wall.
The textbook version of the LLN requires two conditions:
Independence: The random variables (sectoral surplus-value rates) must be statistically independent of each other.
Identical distribution: Each variable must follow the same probability distribution.
Neither condition holds for the real economy. And Gómez Julián is admirably upfront about this. Sectors are deeply intertwined — the steel industry depends on mining, manufacturing depends on steel, services depend on consumer spending powered by manufacturing wages. The idea that one sector’s surplus-value rate has no relationship to another’s is economically unrealistic. Furthermore, different industries have different cost structures, different labor intensities, and different technologies. There is no reason their surplus-value rates should follow the same statistical distribution.
So does this kill the argument? Not at all. In fact, it’s the most intellectually interesting part of the paper.
Non-Classical Varieties: When the Rules Relax
Over the past several decades, mathematicians and econometricians have developed non-classical versions of the LLN that weaken or entirely drop the independence and identical-distribution requirements. Gómez Julián surveys several of these:
Li, Rao, and Wang (1995) showed the LLN holds for random variables arranged on a lattice structure under certain conditions — a structure that, as it happens, economic data naturally exhibits.
Adler and Rosalsky (1987) proved the law for weighted sums of independent, identically distributed random variables belonging to a normalized sum, generalizing the classical case.
Chen and Sung (2016) extended those results further: the variables no longer need to be identically distributed. They only need to be “stochastically dominated” by a single random variable, with certain weighting conditions.
Sung (2011) showed that the strong law can hold even when variables are dependent on each other, provided their probability moments (roughly, their averages and variability) satisfy certain finiteness conditions.
The crucial point: these results collectively tell us that the LLN’s convergence conclusion can survive even when the classical assumptions are substantially violated — which is exactly the situation with sectoral surplus-value rates.
Gómez Julián argues that the economic dynamics described by Smith — workers and capital moving between sectors in response to unequal advantages — are precisely the kind of compensatory dependence mechanism that these non-classical versions accommodate. The variables aren’t independent, but their dependence is structured in a way that still drives convergence.
What the Data Actually Shows
The theoretical argument is compelling, but Gómez Julián doesn’t stop there. He turns to sixty years of U.S. data (1960–2020), sourced from the Bureau of Economic Analysis (BEA), to see what the empirical evidence says.
He calculates sectoral surplus-value rates using macroeconomic data on gross operating surplus (representing surplus labor time) and employee compensation (representing necessary labor time), following a standard operationalization of Marx’s categories. After carefully determining which sectors qualify as “productive” in the Marxist sense — a nontrivial task, since the service sector includes activities with very different relationships to surplus-value production — he arrives at 36 productive sectors.
Here’s what the statistical analysis found:
Finding 1: No Identical Distributions
A probability distribution fitting exercise (using the Bayesian Information Criterion) revealed that the 36 sectors’ surplus-value rates follow a patchwork of different distributions — Log-Normal, Cauchy, Uniform, Weibull, and Logistic — with none following a normal distribution. The identical-distribution requirement of the classical LLN is not met.
Finding 2: No Statistical Independence
A Pearson correlation analysis across all 630 possible sector pairs yielded a mean correlation of about 0.08 and a median of about 0.14. While these may look small, a deeper cut reveals that roughly 40% of sector pairs have correlations of 0.3 or above — a level that’s practically meaningful. The sectors are not independent. This makes intuitive sense: industries are connected through supply chains, labor markets, and shared macroeconomic conditions.
Finding 3: Differences Tend Toward Zero
This is the key finding. When Gómez Julián computed the differences between each sector’s surplus-value rate and the global average (both the mean and the median), he found that these differences exhibit a strong tendency toward reciprocal nullification — positive differences roughly cancel out negative ones. The sum of all differences relative to the global mean was essentially zero (on the order of 10⁻¹⁴). The mean of differences relative to the global median was 0.0012 — vanishingly small.
Distributional fitting on these differences revealed they follow a Cauchy distribution (when measured against the global mean) or a uniform distribution (against the global median), with the medians of these distributions sitting very close to zero.
In plain language: sectors deviate from the average in different directions, and those deviations largely cancel each other out.
Why This Matters
Gómez Julián’s paper makes two types of contributions that are worth distinguishing:
For Marxist political economy: If the uniformity assumption holds — even approximately, even as a tendency rather than an iron law — then a large body of research on the long-run behavior of the average rate of profit, both within countries and across the global economy, is on sounder footing than critics have suggested. Researchers studying the tendency of the rate of profit to fall (or not) can continue to work without needing to explicitly model sector-by-sector differences in exploitation rates, at least for aggregate, long-run analyses.
For probability theory and economics: The paper demonstrates a productive intersection between a specific question in political economy and the deep mathematics of convergence theorems. It shows that the non-classical LLN theorems aren’t just abstract curiosities — they have direct relevance to understanding real economic phenomena. The structured dependence between economic sectors isn’t a bug that invalidates the mathematical framework; it’s a feature that the right version of the framework already accounts for.
A Few Honest Caveats — And Why They No Longer Apply
The original 2022 paper was unusually transparent about its limitations, and that transparency is one of its strengths. Rather than forcing the data into inappropriate statistical procedures, it openly acknowledged where the available inferential tools broke down.
At the time, three important caveats remained.
First, formal hypothesis testing had to be abandoned.
The reason was purely statistical rather than economic. Classical inferential procedures—Student’s t tests, Wilcoxon tests, and most conventional non-parametric alternatives—are built on assumptions that the data simply did not satisfy. Sectoral surplus-value rates are neither independent nor identically distributed. They are linked through supply chains, technological change, capital mobility, and macroeconomic shocks. Even bootstrap procedures could not fully solve the problem because ordinary resampling may weaken dependence between resamples while leaving the internal dependence structure fundamentally unchanged. Consequently, the 2022 paper relied primarily on descriptive statistics together with probability-theoretic arguments instead of formal significance testing.
Second, the classification of productive sectors inevitably involved theoretical judgment.
Although the paper carefully justified the inclusion and exclusion of economic activities using Marxian categories and modern national accounting, reasonable scholars could still debate where certain services belong within the circuit of capital.
Third, the empirical evidence came exclusively from the United States.
The descriptive regularities were remarkably strong, but demonstrating that the same convergence mechanism operates under different institutional settings naturally remained an empirical question.
Those were genuine limitations in 2022.
Today, however, the first—and arguably the most important—of them has largely been overcome.
A much more comprehensive methodological paper (Gómez Julián, 2026; SSRN 5172185) develops an entirely new inferential framework specifically designed for exactly the type of data that made the original analysis difficult: dependent, heterogeneous, and unbalanced observations. Instead of trying to force classical statistical tests to work outside the assumptions under which they were derived, the newer paper constructs hypothesis testing from the ground up for this class of problems.
The key innovation is recognizing that the convergence of sectoral surplus-value rates is fundamentally a law-of-large-numbers problem under dependence, not an independent-samples problem. The framework therefore combines three complementary asymptotic structures—triangular arrays (TAC), correlation-weighted sums (WSC), and mixingale processes (MPC)—which respectively model hierarchical dependence, contemporaneous intersectoral dependence, and temporal dependence. Rather than treating these as competing approaches, the paper proves conditions under which they become metrically equivalent and therefore support the same inferential conclusions.
The inferential consequences are substantial.
Instead of abandoning significance testing because dependence invalidates classical procedures, the new framework explicitly extends the Neyman-Pearson paradigm to dependent observations, derives dependence-aware confidence regions, establishes rigorous Type I error control under strong-mixing assumptions, and integrates Bayesian and frequentist inference into a single coherent architecture. Robust procedures—including fixed-b heteroskedasticity-and-autocorrelation-robust inference, block bootstrap techniques that preserve dependence, adaptive conformal inference, composite and Whittle likelihoods, and hierarchical Bayesian estimation—serve as mutually reinforcing validation mechanisms rather than isolated alternatives.
In other words, what had been acknowledged as a methodological limitation in the 2022 paper became the central research question of the later work.
Rather than concluding that inference was impossible under dependence, the subsequent research asks a more fundamental question: what should hypothesis testing look like when dependence is the normal state of the data rather than an exception? The result is a unified inferential framework specifically intended for datasets that violate the assumptions of classical statistics—precisely the situation encountered with sectoral surplus-value rates.
The other caveats remain, although they are considerably less problematic than before. The classification of productive sectors continues to depend on theoretical interpretation, because that issue belongs to political economy rather than statistics. Likewise, expanding the empirical analysis to additional countries remains a desirable avenue for future research. Yet the principal statistical objection—that no valid inferential procedure existed for dependent sectoral data—has now been directly addressed through a purpose-built mathematical framework.
Looking back, the 2022 paper can therefore be read as identifying an important statistical obstacle, while the later work attempts to remove it. Together, the two papers form a coherent research program: first demonstrating that the convergence hypothesis is theoretically plausible and descriptively supported, and then developing the inferential machinery required to test that hypothesis rigorously without relying on unrealistic assumptions of independence or identical distributions.
Gómez Julián, J.M. (2022). Sobre la validez del supuesto de uniformidad en las tasas de plusvalía sectorial desde la teoría de las probabilidades. Ciencia Económica, 11(17).DOI: 10.22201/fe.24484962e.2022.11.17.2
Gómez Julián, J.M. (2026). Hypothesis Testing for Dependent Variables with Unbalanced Data: A Unified Framework: Theory, Robustness, and Software. SSRN Electronic Journal. DOI: 10.2139/ssrn.5172185.
You can also find this libraryat CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
An introduction to the R package that brings together bivariate hurdle regression, horseshoe regularization, and multi-method causal inference under one roof.
The Problem Nobody Teaches You in Grad School
You have two count variables measured over time. Maybe they are insurgent attacks and counterinsurgent operations in a conflict zone. Maybe they are weekly disease case counts and mortality figures. Maybe they are criminal incidents and arrests in a city.
You know the standard toolkit: Poisson regression, negative binomial if there is overdispersion, maybe a zero-inflated model if zeros pile up. But then the real world intervenes:
Both series have far too many zeros. Not the kind of “a few extra zeros” that zero-inflated models handle — the kind where 60–80% of observations are zero, and the rest are counts. This is the hallmark of hurdle data: a process that first decides whether anything happens at all, and then, if it does, decides how much.
The two series influence each other over time. Do counterinsurgent operations today predict insurgent attacks next month? Does disease incidence Granger-cause mortality with a lag? You need cross-lagged dynamics, not two separate regressions.
You have many potential predictors — economic indicators, population density, climate variables, regime dummies — and your sample size is modest. Overfitting is a real danger.
You want causal claims, not just correlations, and a single method is not enough. You need converging evidence.
If this sounds familiar, the R package bivarhr was built for exactly your situation.
What bivarhr Does
bivarhr is an open-source R package (MIT-licensed, available on GitHub) that provides a unified workflow for:
Bivariate hurdle negative binomial regression — jointly modeling two zero-heavy count series with separate zero-generating and count-generating processes.
Horseshoe priors — automatic Bayesian regularization that shrinks irrelevant coefficients toward zero while letting strong signals survive.
Bayesian Model Averaging (BMA) via stacking — combining predictions across many candidate models (different lag orders, different regularization strengths) instead of betting everything on one.
Multi-method causal inference — transfer entropy, Dynamic Bayesian Networks, Hidden Markov Models, VARX models, synthetic control, and sensitivity analysis.
Rigorous validation — temporal placebo tests, rolling out-of-sample evaluation, extreme bounds analysis, and counterfactual average treatment effects.
All of this runs on top of Stan via the cmdstanr interface, which means full Hamiltonian Monte Carlo sampling with the No-U-Turn Sampler (NUTS) — the gold standard for Bayesian computation.
The Core Model: Bivariate Hurdle Negative Binomial
Let us unpack what the model actually does, without unnecessary jargon.
The Hurdle Idea
A hurdle model splits the data-generating process into two questions:
Did anything happen? (Zero vs. non-zero.) This is modeled with a Bernoulli distribution and a logit link — essentially a logistic regression.
If something happened, how much? (Positive counts only.) This is modeled with a truncated negative binomial distribution and a log link.
Mathematically, for each time point and series :
where is the probability of crossing the hurdle (a non-zero count), is the conditional mean of the negative binomial, and is the dispersion parameter.
The key distinction from zero-inflated models: a hurdle model does not distinguish between “structural zeros” and “sampling zeros.” It simply asks whether the count is zero or not, and if not, how large. This is conceptually cleaner and often more appropriate for event count data.
Bivariate Means Cross-Lags
The “bivariate” part means that the model is estimated jointly for two series, and . What makes this powerful is that the design matrices for can include lagged values of , and vice versa. Four specifications are available:
Specification
C → I
I → C
Interpretation
A
Yes
No
Granger-causes
B
No
Yes
Granger-causes
C
Yes
Yes
Bidirectional causality
D
No
No
No cross-series effects
By comparing the predictive performance of these specifications, you can formally test whether one series has leading information about the other — a Granger-causality test conducted within a fully Bayesian framework.
Horseshoe Priors: Letting the Data Decide
When you have many candidate predictors and a modest sample size, ordinary maximum likelihood will happily overfit. The horseshoe prior solves this elegantly.
Each coefficient gets a prior that is a mixture of a tight “spike” near zero and a broad “slab” away from zero. The balance between spike and slab is governed by:
A global shrinkage parameter that controls overall sparsity — how many coefficients the model expects to be non-zero.
A local shrinkage parameter for each coefficient — allowing individual signals to escape shrinkage if the data support it.
The regularized version used by bivarhr (following Piironen & Vehtari, 2017) adds a slab regularization term that prevents the unbounded coefficient estimates that can afflict the original horseshoe. The result: noise is suppressed, genuine signals are preserved, and you do not need to manually select variables.
In practice, you set a prior guess for the fraction of coefficients you expect to be non-zero (e.g., 0.1 for strong sparsity, 0.5 for moderate). The model does the rest.
Bayesian Model Averaging: Stop Picking One Model
One of the most consequential analytical choices is the lag order: how many past time steps of should enter the equation for ? Different lag orders can tell different stories.
Rather than picking one, bivarhr uses stacking (Yao et al., 2018) to combine predictions across many models. The algorithm finds optimal weights such that the combined predictive distribution:
maximizes the expected log predictive density (ELPD), estimated via Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO).
This means the final inference is not tied to a single model. Models with good out-of-sample predictive performance get high weight; poor models are down-weighted automatically. The stacking weights themselves are informative: if specification A () consistently gets higher weight than specification D (no cross-lags), that is evidence of a genuine cross-series effect.
The Causal Inference Toolkit
A bivariate regression with cross-lags is a necessary but not sufficient condition for causal claims. bivarhr therefore provides six complementary causal inference methods:
Transfer Entropy
Transfer entropy measures how much knowing the past of series reduces your uncertainty about the current value of series , beyond what ‘s own past already tells you. It is an information-theoretic generalization of Granger causality that makes no linearity assumptions.
bivarhr computes transfer entropy for three data transformations — raw counts, rates (counts per unit exposure), and binary presence/absence — each pre-whitened via an appropriate GLM to remove confounding trends. Permutation-based significance tests (with Benjamini-Hochberg correction) control the false discovery rate.
Dynamic Bayesian Networks (DBN)
A DBN learns a directed acyclic graph over discretized versions of , , and any regime variable, with edges restricted to flow from to . The learned structure reveals which variables are direct parents of which — does yesterday’s regime predict today’s attacks, or only through yesterday’s operations?
Hidden Markov Models (HMM)
An HMM with multivariate Poisson emissions assumes that both and are driven by a latent (hidden) state that evolves over time. The inferred state sequence can reveal regimes — a “calm” state with low counts for both series, an “escalation” state with high counts, and so on — without you having to define them in advance.
VARX Models
A Vector Autoregression with Exogenous variables provides a classical time-series check. bivarhr fits bivariate VAR models and reports stability, serial correlation, normality, and ARCH diagnostics — a complementary perspective to the Bayesian hurdle model.
Sensitivity Analysis
Using the framework of Cinelli and Hazlett (2020), bivarhr quantifies how robust an OLS regression result is to an unobserved confounder. The output is a “robustness value”: how strong would an omitted confounder need to be (in terms of partial ) to explain away the estimated effect? This gives a concrete, interpretable measure of fragility.
Synthetic Control
A Bayesian Structural Time Series (BSTS) model constructs a synthetic counterfactual: what would series have looked like in the absence of a treatment (e.g., a policy intervention, a conflict escalation)? The difference between observed and counterfactual is the estimated causal effect, with full posterior uncertainty intervals.
Validation: Proving Your Results Are Not Artifacts
Temporal Placebo Test
Shuffle the time indices of your data, refit the model, and compare the ELPD to the original. If the original model’s ELPD is substantially higher, the temporal structure is genuine — not an artifact of overfitting.
Rolling Out-of-Sample Evaluation
Split the data at multiple cut points (60%, 70%, 80%, 90% of the sample), fit the model on the training portion, and forecast the remaining observations. The RMSE of these forecasts gives an honest assessment of predictive accuracy.
Extreme Bounds Analysis (EBA)
Re-fit the model with every combination of control variables and check whether the key coefficients (e.g., cross-lag effects) remain stable. If a coefficient flips sign or loses significance when you add or remove a single control, it is fragile.
G-Computation for Average Treatment Effects
Counterfactual ATEs answer the question: “What would happen to if we set the cross-lags and transition variables to zero?” This is the causal estimand that many applied researchers actually care about, computed via posterior predictive simulation.
Getting Started
Installation
# Install the package from GitHub
devtools::install_github("isadorenabi/bivarhr")
# Install cmdstanr (required for model fitting)
install.packages("cmdstanr",
repos=c("https://stan-dev.r-universe.dev",
getOption("repos")))
cmdstanr::install_cmdstan()
# Optional: install packages for full causal inference functionality
The applications extend well beyond conflict analysis. Any domain with paired zero-heavy count time series is a candidate:
Public health: Disease incidence and mortality; hospital admissions and readmissions; vaccination rates and outbreak counts.
Criminal justice: Criminal incidents and law enforcement responses; drug seizures and overdose deaths.
Economics: Firm entry and exit counts; patent applications and citations; trade flow counts between countries.
Ecology: Predator and prey counts; species occurrence and extinction events; colonization and local disappearance.
Social science: Protest events and government responses; legislative proposals and vetoes; online misinformation posts and fact-checking activity.
Under the Hood: Technical Details
For those who want to look inside:
Component
Implementation
Sampler
Stan NUTS (No-U-Turn Sampler) via cmdstanr
Shrinkage
Regularized horseshoe (Piironen & Vehtari, 2017)
Model comparison
with stacking weights (Yao et al., 2018)
Convergence
, , zero divergences
LOO reliability
Pareto for all observations
Dispersion
Truncated log-normal prior on
The Stan model code is fully transparent — you can inspect it via get_hurdle_model() — and the generated quantities block produces posterior predictive checks, pointwise log-likelihoods, and fitted values out of the box.
Diagnostic Checklist
Before trusting your results, bivarhr encourages you to verify:
Metric
Target
Warning sign
< 1.01
> 1.05 means chains did not converge
(bulk)
> 400
< 400 means inefficient sampling
Divergences
0
Any divergence suggests geometric pathologies
Pareto
< 0.7
> 0.7 means LOO approximation is unreliable
difference
> 5 for “strong” preference
< 2 means models are practically equivalent
The Bigger Picture
What makes bivarhr distinctive is not any single component — hurdle models exist, horseshoe priors exist, transfer entropy exists — but the integration. The package recognizes that in applied work, the bottleneck is not any one method’s sophistication. It is the workflow: getting from raw data to a validated causal claim, with uncertainty quantification, in a reproducible pipeline.
By combining a carefully specified bivariate hurdle model with automatic regularization, principled model averaging, and six complementary causal inference methods — all validated through placebo tests and out-of-sample evaluation — bivarhr gives applied researchers a framework that is both rigorous and practical.
If your work involves paired count time series with excess zeros and you care about causal inference, this package deserves a serious look.
You can also find this libraryat CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
An introduction to an R package that brings Bayesian inference, panel data econometrics, and rigorous validation to one of political economy’s most enduring empirical debates.
Why This Package Exists
Here is a question that has occupied economists for over two centuries: when you pay for something, does the price you pay bear any systematic relationship to the labor required to make it?
Adam Smith thought so. David Ricardo refined the idea. Karl Marx built an entire theory of exploitation on it. And since the mid-twentieth century, empirical researchers have been trying to measure the strength of this correspondence with real-world data.
The challenge has always been methodological. The datasets are panel data — prices observed across many economic sectors over many time periods — and they demand techniques that respect both the cross-sectional structure (different industries behave differently) and the temporal dimension (relationships can shift over time). A simple scatterplot of values against prices, however illustrative, will not settle the question.
valueprhr is an R package built to close this methodological gap. It provides a complete, reproducible pipeline: from raw price matrices to model estimation, from Bayesian inference to out-of-sample validation, from structural break detection to side-by-side model comparison. It was designed for political economy, but as we will see, its toolkit applies to any panel data problem where you need to assess the correspondence between two variables across entities and time.
The Core Idea (in Plain Language)
In the classical and Marxian tradition, the value of a commodity is determined by the total labor time — direct and indirect — required to produce it. If a table requires 10 hours of socially necessary labor and a chair requires 5, the table’s value is twice the chair’s.
This gives rise to what economists call direct prices (denoted pd): prices that are strictly proportional to the labor embodied in each commodity. They represent what prices would be if they perfectly mirrored labor content.
But capitalism does not work that way. Capital flows between sectors seeking the highest return, and competition tends to equalize the rate of profit across industries. The prices that emerge from this process are called prices of production (denoted pπ). They redistribute surplus value: sectors with higher organic composition of capital (more machinery relative to labor) tend to have prices of production above their direct prices, and vice versa.
The central empirical question is: despite this redistribution, how closely do direct prices and prices of production correspond?
The standard test is a log-linear regression:
ln(pπit) = α + β · ln(pdit) + uit
where i indexes sectors and t indexes time periods.
Three hypotheses are at stake:
β ≈ 1: a one-percent increase in direct prices is associated with roughly a one-percent increase in production prices (proportionality).
R² ≈ 1: direct prices explain the vast majority of the variation in production prices.
Stability: the relationship holds consistently across time periods.
If all three hold, the labor theory of value has strong empirical support. valueprhr gives you the tools to test each one rigorously.
What’s Inside the Package
valueprhr organizes its functionality into six modules. Here is what each does and why it matters.
1. Data Preparation
Real-world data rarely arrives in the format econometric methods require. The package accepts two data frames in wide format (rows = years, columns = sectors) — one for direct prices, one for production prices — and converts them into the long-format panel structure that econometric models expect.
#> Year Sector direct production log_direct log_production
#> 1 1960 Agriculture 45.2 48.1 3.81 3.87
#> 2 1961 Agriculture 46.0 49.0 3.83 3.89
#> ...
The function prepare_log_matrices() does the same job but returns matrix format, which is what the multivariate methods (PLS, CCA) need.
2. Panel Data Models
This is where the core econometrics happens. The package implements two complementary specifications:
Two-Way Fixed Effects (FE) controls for both sector-specific and time-specific unobserved heterogeneity:
Yit = αi + γt + β · Xit + εit
In plain terms: every sector has its own baseline (some sectors are systematically more expensive), every year has its own macroeconomic conditions (inflation, crises), and the model isolates the within variation to estimate the core relationship.
#> Estimate = 0.9754, SE = 0.0123, t = 79.30, p = 0.0000
The cluster_type = "group" option computes cluster-robust standard errors at the sector level, which accounts for serial correlation within each sector’s time series.
Mundlak Correlated Random Effects (CRE) takes a different route. Instead of dummy variables for every sector, it decomposes the predictor into a within-sector component (how Xit deviates from sector i‘s average) and a between-sector component (the sector average itself):
Yit = α + βW · (Xit − X̄i) + βB · X̄i + ui + εit
In data science language: this is a way to control for group-level confounders without the computational cost of N dummy variables. If βW = βB, the within and between effects are the same, and a simpler Random Effects model suffices. If they differ, the relationship between values and prices operates differently within a sector over time than across sectors.
#> -> Fail to reject H0: RE/CRE specification is consistent
The function test_mundlak_specification() formalizes this check. A low p-value means you should stick with Fixed Effects; a high p-value means the simpler model is adequate.
The package also includes a Panel Granger Causality test (the Dumitrescu-Hurlin procedure), which tests whether past values of direct prices help predict current production prices — and vice versa.
panel_granger_test(panel, lags=c(1, 2))
#> direction lag W_stat Z_stat p_value significant
#> 1 direct -> production 1 8.432 3.126 0.0018 TRUE
#> 2 direct -> production 2 6.215 2.441 0.0146 TRUE
#> 3 production -> direct 1 5.890 2.103 0.0354 TRUE
#> 4 production -> direct 2 4.012 1.332 0.1828 FALSE
3. Bayesian Models
Classical (frequentist) estimation gives you a single point estimate for β. Bayesian methods give you a full probability distribution over possible values, incorporating your prior beliefs and updating them with the data.
In econometric language: instead of β̂ = 0.975 ± 0.012, you get a posterior distribution showing that β lies between 0.95 and 1.00 with 95% probability.
The package offers two Bayesian approaches:
Sector-by-Sector Bayesian GLM fits an independent Bayesian linear model for each sector, using weakly informative priors (the rstanarm package handles the MCMC sampling via Stan). Each sector gets its own slope and intercept, along with Leave-One-Out Cross-Validation (LOO-CV) scores.
In data science language: LOO-CV is a principled way to assess out-of-sample predictive performance without holding out data. The LOOIC (LOO Information Criterion) is the Bayesian analogue of AIC — lower is better.
Bayesian Hierarchical Model goes further by pooling information across sectors. Instead of treating each sector in isolation, it assumes that sector-specific slopes are drawn from a common population distribution:
βi ~ N(μβ, σβ2)
Sectors with less data “borrow strength” from the population mean. This is especially valuable when some sectors have short time series.
When the number of sectors (N) is large relative to the number of time periods (T), standard regression becomes unstable. This is the “small T, large N” problem common in panel data. The package offers three multivariate techniques to handle it:
Partial Least Squares (PLS) extracts latent components that explain covariance between direct prices and production prices. It handles multicollinearity gracefully and is widely used in chemometrics, genomics, and now in value-price analysis.
Canonical Correlation Analysis (CCA) finds linear combinations of direct prices and production prices that are maximally correlated. In econometric language: CCA extracts the “shared economic signal” — the common factor driving both sets of prices.
Standard k-fold cross-validation violates temporal ordering. If you train on 1960–1990 and test on 1985–1990, future information leaks into the training set. The package implements two time-aware approaches:
Rolling Window CV trains on t₀ … tW, tests on tW+1 … tW+H, then rolls the window forward.
cv<-rolling_window_cv(
panel,
window_sizes=c(20, 30),
step_size=2,
test_horizon=3
)
print(cv$summary)
Leave-One-Sector-Out (LOSO) trains on all sectors except one and predicts the held-out sector. This tests cross-sectional generalization: does the value-price relationship estimated from other sectors hold for agriculture? For mining? For finance?
loso<-leave_one_sector_out(panel)
print(loso$summary)
#> metric mean sd
#> 1 RMSE 0.04521 0.01832
#> 2 MAE 0.03587 0.01456
#> 3 R_squared 0.96120 0.02340
An average R² above 0.96 in LOSO-CV means the relationship generalizes robustly across sectors.
6. Structural Break Tests
Has the value-price relationship been stable over time? Or did it shift at some point — due to globalization, a methodological change in data construction, a technological revolution, or a regime shift in profit rate equalization?
The package aggregates the panel to a time series and applies a battery of tests:
A non-significant result is actually good news here: it means the value-price correspondence has been structurally stable across the entire sample period.
The Full Pipeline in One Command
If you want to run everything at once — data preparation, FE and CRE models, cross-validation, structural break tests, and model comparison — the package offers a single entry point:
results<-run_full_analysis(
direct,
production,
run_bayesian=FALSE, # Set TRUE if you have rstanarm installed
run_cv=TRUE,
run_breaks=TRUE,
verbose=TRUE
)
# Access everything
print(results$comparison)
print(results$cv_summary)
cat(format_break_results(results$breaks))
You can then export the comparison table and CV results to CSV:
export_results_csv(
results$comparison,
results$cv_summary,
output_dir="results/"
)
Beyond Political Economy: General Panel Data Applications
Although valueprhr was built for the specific question of value-price correspondence, its methods are general-purpose panel data tools. Any research problem involving the relationship between two variables observed across entities and time can benefit from the package:
Health economics: Does out-of-pocket spending track underlying treatment costs across regions over time?
Environmental economics: Do carbon prices reflect the embodied emissions of goods across industries?
Education: Do standardized test scores correspond to instructional expenditure across school districts over decades?
Finance: Do book values predict market valuations across sectors?
Any two-variable panel regression where you need fixed effects, Mundlak decomposition, robust standard errors, time-aware cross-validation, or structural break detection.
The key requirement is that your data has a panel structure (entities × time) and that the distributional assumptions of the models are reasonable for your context. The methods — two-way FE, Mundlak CRE, Bayesian hierarchical models, PLS, CCA, rolling-window CV, structural break tests — are econometric staples that transcend any particular application domain.
Installation and Dependencies
The package requires R ≥ 4.1.0. Core functionality depends only on base R and the Metrics package. Extended features (Bayesian models, panel data infrastructure, structural break tests) are handled through soft dependencies that are loaded on demand:
# Install from GitHub
install.packages("devtools")
devtools::install_github("isadorenabi/valueprhr")
# Optional: install all suggested packages at once
Note for Bayesian models:rstanarm requires a working C++ toolchain — Rtools on Windows, Xcode Command Line Tools on macOS, or build-essential on Linux.
What Makes This Package Methodologically Different
Three features distinguish valueprhr from a hand-rolled analysis:
Time-aware validation. Most applied work reports in-sample R² as evidence of fit. valueprhr pairs every model with rolling-window and leave-one-sector-out cross-validation, giving you out-of-sample performance that is honest about temporal dependence and cross-sectional generalization.
The Mundlak decomposition. By splitting effects into within-sector and between-sector components, the package lets you test whether the value-price relationship operates at the sector level (structural), at the temporal level (cyclical), or both. This is a nuance that most empirical studies in this literature overlook.
Bayesian hierarchical pooling. Sectors with short time series are a common headache. The hierarchical model lets small sectors borrow statistical strength from the population, producing more stable estimates than independent sector-by-sector regressions.
A Note on the Underlying Data
The wiki documentation mentions that market price indices used in this framework are constructed by temporally disaggregating the aggregate Consumer Price Index using the Input-Output matrix as a structural indicator, relying on closed-form Bayesian solutions from the BayesianDisaggregation library. This is a methodological detail worth understanding: the sectoral prices are not raw market quotes but statistically consistent decompositions of the macroeconomic aggregate. This ensures that the estimated sectoral price movements add up to the observed CPI, a property that many ad hoc sectoral price datasets lack.
Citation
If you use valueprhr in your research:
@software{gomezjulian2025valueprhr,
author = {Gómez Julián, José Mauricio},
title = {valueprhr: Value-Price Analysis with Bayesian and Panel Data Methods},
year = {2025},
url = {https://github.com/isadorenabi/valueprhr},
note = {R package version 0.1.0}
}
Author: José Mauricio Gómez Julián — ORCID — isadore.nabi@pm.me
The labor theory of value is either one of the most important ideas in the history of economics or one of the most contested. Either way, it deserves better tools than a spreadsheet and a prayer. valueprhr brings the full machinery of modern econometrics to the question — and lets the data speak for itself.
You can also find this libraryat CRAN and download it directly from R and RStudio.
LISTEN TO THIS POST AS A PODCAST:
You found a relationship in your data. It’s statistically significant. The coefficients are clean. But here’s the question that keeps rigorous researchers up at night: will this relationship still hold six months from now?
That question — whether a relationship between economic variables is real and durable, not just a fluke of the particular sample you happened to look at — is exactly what the EconCausal R package was built to answer.
What Problem Does EconCausal Solve?
Causal inference in time series data is notoriously tricky. Standard regression can tell you that two variables move together. It cannot tell you why, or whether that co-movement will persist when the economic regime shifts — which it inevitably will.
EconCausal addresses this challenge by combining three well-established econometric methodologies into a single, reproducible pipeline. What makes it distinctive is not any single technique — each component has solid academic pedigree — but rather the orchestration: a standardized protocol that forces every candidate relationship through rigorous pre-testing, proper temporal validation, and explicit decision rules before declaring a finding.
Think of it as a quality-control assembly line for causal claims in time series.
The Three Approaches
EconCausal implements three methodological frameworks, each capturing a different aspect of what “causality” can mean in temporal data. Using all three gives you complementary perspectives rather than relying on a single method’s assumptions.
1. Error Correction Models with MARS (ECM-MARS)
The idea: Many economic variables share a long-run equilibrium relationship but deviate from it in the short run. An Error Correction Model (ECM) captures both dimensions — the long-run cointegrating relationship and the short-run adjustment dynamics.
EconCausal enhances the classical ECM by replacing its usual linear regression engine with Multivariate Adaptive Regression Splines (MARS), a flexible non-parametric method developed by Jerome Friedman in the early 1990s. This matters because real-world adjustment mechanisms are often nonlinear: variables may correct slowly when deviations are small but snap back aggressively when deviations cross certain thresholds. MARS can detect these threshold effects automatically without you having to specify them in advance.
What the protocol does in practice:
Tests whether the variables are integrated of order one (I(1)) — a prerequisite for cointegration analysis
Runs both the Engle-Granger and Johansen cointegration tests, applying an “either” rule: if either test finds evidence of cointegration, the analysis proceeds
Controls for serial correlation and heteroskedasticity using HAC-consistent standard errors
Fits the ECM with MARS as the regression engine
Validates the relationship using rolling-origin cross-validation
Evaluates the model using a dual criterion (more on this below)
Why two cointegration tests? Because Engle-Granger and Johansen have complementary strengths. Engle-Granger is intuitive and straightforward but can miss cointegrating relationships in multivariate settings. Johansen’s approach is more powerful in higher dimensions but relies on stricter distributional assumptions. Requiring only one of them to flag significance is a pragmatic compromise that reduces false negatives without being reckless about false positives.
2. Bayesian Structural Time Series (BSTS)
The idea: Decompose a time series into interpretable structural components — trend, seasonality, regression effects — and do so within a fully Bayesian framework that quantifies uncertainty at every level.
BSTS, originally popularized by Google’s CausalImpact package and the work of Scott and Varian (2014), models the data-generating process as a state-space system. The key innovation in EconCausal’s implementation is the use of spike-and-slab priors for automatic variable selection. This is a Bayesian regularization technique that effectively asks: “Which candidate predictors actually belong in this model?” — and it does so probabilistically, rather than relying on stepwise procedures or arbitrary p-value cutoffs.
What makes it useful for causal inference:
The structural components absorb confounding patterns (trends, seasonality) that might otherwise masquerade as causal relationships
Spike-and-slab priors guard against overfitting by aggressively shrinking irrelevant predictors toward zero
The Bayesian framework produces probabilistic forecasts — not just point predictions, but full posterior distributions — which naturally give you calibrated prediction intervals
It handles the “what would have happened without the intervention” question, which is central to impact evaluation
3. Bayesian GLM with AR(1) Errors (BGLM-AR1)
The idea: Fit a Bayesian generalized linear model but explicitly account for temporal dependence in the residuals through a first-order autoregressive error structure.
This might sound like a niche technical detail, but it addresses one of the most common mistakes in applied econometrics: ignoring autocorrelation in regression residuals. When residuals are correlated over time, standard errors are wrong, confidence intervals are too narrow, and you end up being more confident in your results than you should be.
By modeling the residual autocorrelation directly with an AR(1) structure and using Hamiltonian Monte Carlo (HMC) sampling via Stan, the BGLM-AR1 approach produces valid inference even when the data have strong temporal dependence. The priors are weakly informative and calibrated for standardized variables, which means they provide regularization without dominating the data.
The dual evaluation criterion: This is where EconCausal’s protocol really distinguishes itself. A relationship is only accepted if it improves both:
Predictive density (measured by Expected Log Predictive Density, or ELPD, estimated via Pareto-Smoothed Importance Sampling Leave-One-Out cross-validation — PSIS-LOO)
Point forecast accuracy (measured by RMSE)
Requiring improvement on both dimensions simultaneously is stricter than most practitioners apply. A model might improve point forecasts while worsening the calibration of its uncertainty estimates, or vice versa. EconCausal demands both.
What Makes EconCausal Different: The Protocol
You might be thinking: “Most of these individual techniques already exist in other R packages — brms and Stan for Bayesian modeling, the ecm package for error correction models, bsts and CausalImpact for structural time series.” And you’d be right.
The novelty of EconCausal is not in inventing new statistical methods. It lies in integrating them into a single, standardized protocol with explicit decision rules, thresholds, and validation procedures that run automatically.
Here’s what that means concretely:
Temporal Stability Validation
Rather than evaluating a model on a single train-test split, EconCausal uses Leave-Future-Out (LFO) cross-validation with configurable windows and horizons. The system repeatedly:
Trains on data up to a cutoff point
Forecasts into the future
Evaluates forecast accuracy
Slides the cutoff forward
Repeats
A relationship is only declared “stable” if it performs consistently across multiple temporal folds — not just one lucky window. The package reports a support metric: the proportion of temporal folds where the relationship holds up.
Dual Decision Rules with Explicit Thresholds
Every candidate relationship must clear both the Bayesian bar (ELPD improvement) and the frequentist bar (RMSE reduction). These are not vague guidelines — the protocol applies explicit numerical thresholds and includes tie-breaking rules for ambiguous cases.
A Single Reproducible Pipeline
From data loading through pre-testing, model fitting, temporal validation, and final go/no-go decisions — everything runs in one pipeline. This eliminates the researcher degrees of freedom that plague so many applied studies, where hundreds of analytical choices are made invisibly along the way.
An independent novelty assessment included in the package’s documentation compared this protocol against both the academic literature and existing software frameworks. The finding: while all the building blocks are state-of-the-art and well-documented in the literature, no existing academic paper or software package was found that implements the same end-to-end protocol with the same combination of pre-tests, dual decision rules, support thresholds, and temporal validation.
Approach
Components = State of the Art?
Identical Protocol in Literature?
Identical Protocol in Existing Software?
BGLM-AR(1)
Yes
No
No
ECM-MARS
Yes
No
No (though ecm + earth approximates the model-fitting step)
BSTS
Yes
No
No
Who Should Use EconCausal?
The package was originally designed for economic research — specifically, examining production-circulation relationships — but its applicability extends to any domain where you need to answer questions like:
Does X cause Y, or is it the other way around? The temporal structure of the methods naturally distinguishes between “X leads Y” and “Y leads X.”
Is this relationship robust across time? The rolling-origin validation directly tests temporal stability.
Are there nonlinear dynamics at play? The MARS component in ECM-MARS can detect threshold effects and nonlinear adjustment speeds.
How certain can we be? The Bayesian components (BSTS and BGLM-AR1) provide full uncertainty quantification, not just point estimates.
This makes it useful for:
Academic researchers conducting empirical macroeconomics, financial economics, or applied econometrics
Policy analysts evaluating whether proposed policy relationships will hold across different economic conditions
Data scientists in finance who need to move beyond correlation to directionality and stability
Graduate students learning causal inference methods who want a pedagogical tool that enforces methodological discipline
The package depends on a modern Bayesian computing stack — cmdstanr for Hamiltonian Monte Carlo sampling, loo for PSIS-LOO computation, and bsts for structural time series components — so make sure your R environment is up to date.
A Note on What EconCausal Is Not
It’s worth being explicit about the boundaries:
It is not a magic wand. If your data are fundamentally unsuitable for causal inference (too few observations, structural breaks that dominate the signal, measurement error in the variables), EconCausal will not rescue you. Its value lies in applying rigorous methods correctly and transparently — not in overcoming the limitations of your data.
It does not replace domain knowledge. Statistical evidence for a stable, directional relationship is necessary but not sufficient for a causal claim. You still need a theoretical mechanism, an understanding of potential confounders, and judgment about the plausibility of the identifying assumptions.
It does not invent new mathematics. Its contribution is the protocol — the disciplined, reproducible integration of established techniques into a pipeline with explicit, auditable decision rules.
The Bottom Line
EconCausal occupies a valuable niche in the R ecosystem. While individual components of its methodology are available elsewhere — Bayesian GLMs in brms, BSTS in the bsts package, ECM in the ecm package — no other tool packages them into a single protocol that systematically tests temporal stability, applies dual evaluation criteria, and produces reproducible, auditable results.
If your work involves making causal claims from time series data — and you care about whether those claims will survive contact with the future — EconCausal deserves a serious look.
References for the underlying methodologies: Engle & Granger (1987), Johansen (1991), Friedman (1991), Scott & Varian (2014), Vehtari et al. (2017), Hyndman & Athanasopoulos (2021). See the package documentation for the complete bibliography.
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.
