Espartaco

Marxist Philosophy of Science

“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

EnglishEspañol

Tag: statistical software

  • Marx, Adam Smith, and the Law of Large Numbers

    Marx, Adam Smith, and the Law of Large Numbers

    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:

    1. Independence: The random variables (sectoral surplus-value rates) must be statistically independent of each other.
    2. 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.

  • When the Whole Is All You See: A Bayesian Way to Recover the Parts with BayesianDisaggregation

    When the Whole Is All You See: A Bayesian Way to Recover the Parts with BayesianDisaggregation

    You can also find this library at 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.

  • topologyR: Turning Time Series into Shapes to Test What Your Models Quietly Assume

    topologyR: Turning Time Series into Shapes to Test What Your Models Quietly Assume

    You can also find this library at CRAN and download it directly from R and RStudio.

    LISTEN TO THIS POST AS A PODCAST:

    There is a habit so embedded in quantitative work that most practitioners never think to question it. You have a time series — quarterly GDP, an EEG channel, a temperature record — and at some point you fit a smooth curve through it, interpolate a missing value, or estimate a “long-run trend.” All of these moves rest on a single, seldom-checked assumption: that the data form one continuous whole, that a single smooth function can legitimately pass through every point.

    But what if they don’t? What if your series is, structurally, two or three disjoint pieces glued together by the calendar — pieces between which no continuous function can travel? In that case, the spline you just fitted is not an approximation of reality; it is a mathematical fiction painted over a fracture.

    topologyR is an R package that lets you check this before you model. It takes a numeric time series, converts it into a graph, converts that graph into a topological space, and then asks the one question that determines whether global continuous methods are even valid: is this space one connected piece, or several?

    It sounds abstract. It is abstract — but the consequence is concrete. The package is the work of José Mauricio Gómez Julián, and it is open-source, with a GitHub repository, a detailed Wiki, and a companion research paper archived on Zenodo. What follows is a tour of what the package does, why it matters, and where it fits in the broader landscape of topological data analysis.

    The Hidden Assumption

    Think about what happens when you impute a missing value in a time series using a cubic spline. The spline assumes that the points on either side of the gap belong to the same continuous process — that the missing value lies somewhere along a smooth bridge between them. If the series has actually undergone a structural break, a regime change, or a discontinuity between those points, the spline will happily produce a number, and that number will be wrong in a way no confidence interval can capture.

    This is not a niche problem. It appears in econometrics (trend estimation across business cycles), in neuroscience (coherence across brain-state transitions), in climatology (warming trends across regime shifts). The methodological error is always the same: assuming continuity without first verifying that continuity is mathematically possible.

    topologyR’s contribution is to make that verification explicit, parameter-free, and exact.

    From Numbers to Shapes: The Pipeline in Three Steps

    The package’s workflow has an elegant, almost architectural logic. You feed it a series of numbers; it returns a topological verdict. Between input and output, three transformations occur.

    Step 1: The Series Becomes a Graph

    The first move is borrowed from network science: the visibility graph. Imagine your time series plotted as a mountain range — each observation is a peak or a valley at a given time. Two points are connected by an edge if you could stand on one and see the other, with no taller peak blocking the line of sight.

    topologyR implements two flavours. The Horizontal Visibility Graph (HVG) connects two points if every point between them is strictly lower than the shorter of the two — a horizontal line of sight. It runs in linear time and captures the skeleton of the series’ ups and downs. The Natural Visibility Graph (NVG) is more generous: it connects two points if every intermediate point lies below the straight line joining them, regardless of the heights of the endpoints. It is denser, richer, and runs in O(n log n) expected time. The NVG always contains the HVG as a subgraph.

    Both are parameter-free. There is no threshold to tune, no bandwidth to select, no ε to agonise over. The graph is determined entirely by the data’s own geometry. This matters enormously: it eliminates the single largest source of arbitrariness in the entire pipeline.

    Step 2: The Graph Becomes a Topology

    Here is where topologyR departs from ordinary network analysis. A graph tells you who is adjacent to whom. A topology tells you something deeper: what the neighbourhood structure of the entire space looks like — which collections of points form coherent open regions, and how those regions combine.

    The construction follows a method introduced by Nada, El Atik, and Atef in 2018. For each vertex v in the graph, you form its closed neighbourhood — the vertex itself plus all its direct neighbours. This family of closed neighbourhoods serves as a subbase. You then close it under finite intersections to obtain a base, and close the base under arbitrary unions to obtain the full topology.

    If those words feel heavy, think of it this way: the subbase is a rough draft of “who belongs with whom.” Intersecting neighbourhoods refines the draft — “the points that both neighbourhoods agree on.” Taking unions completes the picture — “every region that can be assembled from these building blocks.” The result is a genuine topological space, complete with open sets satisfying the standard axioms, sitting on top of your time series like a scaffolding you didn’t know was there.

    Step 3: The Topology Reveals Its Connectivity

    Now comes the decisive question. A topological space is connected if it cannot be split into two non-empty open pieces — if there is no clean fracture running through it. For finite spaces, there is a beautiful theorem, due to McCord (1966) and Stong (1966), that makes this check exact and tractable. The specialization preorder orders the points by how their neighbourhoods nest inside one another, and the connected components of the resulting structure are precisely the topological connected components.

    The crucial practical point: this computation works directly on the base — the refined building blocks — without ever needing to enumerate the full topology (which can be exponentially large). It runs in polynomial time, and the components it returns are exact, not approximate.

    The Decision Rule

    Everything so far converges on a single, actionable verdict. topologyR hands you a connectedness decision, and that decision has a direct methodological consequence:

    • If the induced topology is connected, then your data are consistent with a single continuous process. Global continuous methods — splines, kriging, polynomial interpolation, moving-average imputation, kernel methods — are mathematically supported. You may proceed.
    • If the induced topology is disconnected, then no single continuous function can cover the entire series. Global continuous methods are invalid by construction. You must segment the series along the connected components the package identifies, and model each piece independently — with regime-switching models, component-wise imputation, or finite mixtures.

    This is the package’s value proposition: a reproducible, topology-first workflow that decides, before you touch a model, whether global continuity is a justified assumption or a silent error.

    Global versus Local

    A subtlety worth flagging: the rule depends on what you are trying to learn. Global properties — a secular trend, a Hurst exponent, total neural synchronisation, a centennial warming signal — depend on relationships among all points and require topological connectivity to be valid. Local properties — instantaneous volatility in a small window, point-to-point rates of change, low-order autocorrelation — are defined on restricted neighbourhoods and remain valid within each connected component, regardless of whether the whole series is one piece or several. The package gives you the component structure to make that distinction operational.

    Time Has an Arrow: Directed Topologies and Irreversibility

    So far, the construction has treated the visibility graph as undirected — time flows, but the edges don’t care which way. That discards information. Time series are inherently directional: time runs from past to future, and many real systems are irreversible — they behave differently forwards and backwards. Economic expansions creep upward over years; recessions collapse in quarters. Neurons fire and recover on different timescales. The undirected graph cannot see this asymmetry.

    topologyR’s directed mode fixes this. With directed = TRUE, each visibility edge is oriented from the earlier time point to the later one, producing a directed acyclic graph (a DAG) in which the time index is a natural topological order. From this directed graph, the package extracts two neighbourhood structures: the forward neighbourhood (who can I see ahead of me?) and the backward neighbourhood (who behind me can see me?).

    Applying the Nada construction to each yields two topologies: a forward topology τ⁺ and a backward topology τ⁻. The pair (X, τ⁺, τ⁻) forms what Kelly (1963) called a bitopological space — a set equipped with two topologies rather than one. The divergence between them is a direct, topological measurement of temporal irreversibility.

    Irreversibility Indices

    In a perfectly reversible process — symmetric dynamics, no privileged direction — the two topologies coincide: τ⁺ ≅ τ⁻. They have the same number of connected components, the same base size, the same connectivity. In an irreversible process, they pull apart.

    topologyR quantifies this with several indices. The component irreversibility measures the normalised difference in the number of connected components between the forward and backward topologies: zero means symmetric, one means maximally asymmetric. The base irreversibility does the same for the sizes of the topological bases. The asymmetry direction — the signed difference in component counts — tells you which way the arrow points: a positive value means the forward topology is more connected (fewer components) than the backward one.

    That last point has a concrete physical interpretation. Consider a time series with gradual expansions and abrupt contractions — the classic shape of a business cycle, where GDP creeps up over years and drops in a quarter. During a gradual rise, forward visibility is relatively unobstructed: looking ahead from a point on the upslope, you can see far. After an abrupt drop, backward visibility is blocked: looking back from the trough, the cliff face hides earlier points. This asymmetry means the forward topology should be more connected than the backward topology — fewer forward components, more backward fragmentation. The package predicts, and the data confirm, a positive asymmetry direction for such series.

    The Alexandrov Layer and the Resolution Hierarchy

    There is a third topology lurking in the directed graph, and it is older than the Nada construction by several decades. The Alexandrov topology τ_A, introduced by Alexandrov in 1937, is the topology whose open sets are the upsets of the reachability relation — the sets that, once you enter them, contain everything reachable downstream. For each vertex, its minimal open set is the collection of all vertices reachable from it via directed paths.

    topologyR computes this efficiently: a reverse-order bitset propagation that processes vertices from last to first, OR-ing reachability sets together in O(nm/64) time, reusing the same high-performance bitset infrastructure as the Nada engine.

    The relationship between the Alexandrov and Nada topologies is precise and informative: τ_A is always a subset of the forward Nada topology. The Alexandrov base captures pure order structure — “who can reach whom” — while the Nada intersection closure generates additional sets that are not upsets, catching finer-grained structure. The difference in base sizes, |B_Nada| − |B_A|, tells you exactly how much extra topological information the Nada pipeline extracts beyond the raw ordering. A large gap means the closure operations are doing real work; a small gap means the order structure already tells most of the story.

    Under the Hood: Performance Without Compromise

    Topological enumeration is, in the worst case, exponential — the number of open sets can in principle double with every additional element. This is an inherent mathematical fact, not a software limitation. But topologyR is engineered so that the decision you actually need — connectedness — never requires that enumeration.

    The connectivity computation works on the base alone, via the specialization preorder, in O(n² · ⌈B/64⌉) time. The C++ backend (via Rcpp) represents every subset as a packed array of 64-bit words, so set operations reduce to machine-level bitwise instructions. A compile-time template dispatch selects single-word operations for series up to 64 points, two-word for up to 128, three-word for up to 192 — zero loop overhead, branch-free. Beyond that, a runtime fallback handles arbitrary sample sizes, and OpenMP parallelisation is available where the build supports it.

    The practical upshot: you can run the connectivity decision on series with thousands of points without ever touching the exponential regime. Safety limits (max_base_sets, max_open_sets) cap the intersection and union closures with informative termination flags, so if a computation does hit resource limits, you know exactly where and why — and the connectivity result remains valid as long as the base closure completes.

    A Real-World Test: Reading the Business Cycle

    The paper accompanying the package applies the framework to quarterly U.S. real GDP growth from 1992 to 2024 — 129 observations spanning over three decades. The bitopological analysis recovers a positive asymmetry direction: the forward topology is more connected than the backward one, exactly as predicted for a series with gradual expansions and abrupt contractions.

    The undirected topology partitions the series into six connected components, each corresponding to a distinct macroeconomic regime. Strikingly, the COVID-19 contraction and its rebound — the deepest and fastest swing in the sample — are classified as a single topological episode: one connected component spanning the collapse and recovery, reflecting the fact that the visibility structure treats the V-shaped episode as one structural unit rather than two separate events.

    This is the kind of insight the package is designed to produce: not a forecast, not a parameter estimate, but a structural classification that tells you where the legitimate boundaries in your data lie — and, critically, whether a global model is appropriate at all.

    Where It Sits: Complementary, Not Competing

    It is important to be clear about what topologyR is not. It is not a general-purpose topological data analysis (TDA) engine. Packages like GUDHI, Ripser, TDAstats, and scikit-TDA compute persistent homology — multi-scale features across all dimensions, capturing higher-order structures (loops, voids) via Betti numbers β₁, β₂ and their persistence across scales. That is a richer and harder enterprise.

    topologyR has a narrower and more focused aim: it zeroes in on β₀ — connectedness — for one-dimensional series, using graph-induced topologies, and it turns that single invariant into an actionable decision rule for method selection. Think of it as a pre-model governance tool: a rigorous gatekeeper that runs before you choose your modelling strategy, telling you whether the continuity assumptions your favourite methods require are actually justified by the data’s structure.

    The two approaches are complementary. For early-warning detection, precursor signals, or multi-channel structure, persistent homology is the right tool. For the binary question “can I legitimately fit a global continuous model to this series?”, topologyR gives a direct, interpretable, and mathematically exact answer. A natural hybrid workflow uses topologyR as a pre-test and persistent homology for deeper multi-scale analysis.

    Honest Limitations

    No tool is universal, and topologyR is transparent about its boundaries:

    1. Graph choice matters. HVG and NVG produce different graphs, and therefore potentially different topologies. The NVG, being denser, tends to produce fewer connected components. The package encourages comparing both and interpreting the difference — the gap itself is diagnostic.
    2. Sampling and noise. Sparse sampling can mimic disconnection; minor overlaps can mimic connection. The connectedness verdict should be treated as prima facie evidence, not absolute truth — especially near the boundary.
    3. β₀ only. The approach focuses on connectedness. It will not capture loops, voids, or higher-order patterns that persistent homology can detect. If your question is about cycles or multi-scale structure rather than fragmentation, you need the heavier machinery.
    4. Enumeration is exponential; connectivity is not. This is handled honestly: the connectivity decision is polynomial and scalable; full topology enumeration (needed for pairwise connectedness in the bitopological sense) is capped by safety limits with transparent reporting.

    The Bigger Picture

    What makes topologyR more than a clever technical exercise is its epistemological stance. It transforms a step that is normally a tacit habit — assuming continuity — into an explicit, testable, mathematical procedure. In doing so, it removes arbitrariness from one of the most consequential decisions in applied quantitative work: the choice between global and segmented methods.

    The package’s central theorem — that the Nada construction extends to directed graphs and yields a bitopological space whose asymmetry quantifies irreversibility — is formalised in Lean 4 against Mathlib, so the mathematical foundation is not merely asserted but machine-checked. The implementation is CRAN-compliant, passes R CMD check --as-cran cleanly, and ships with 68 unit tests covering visibility graphs, topology generation, connectivity, directed topology, Alexandrov topology, and bitopological analysis.

    For anyone who works with time series and has ever fitted a spline, run a kriging model, or estimated a trend — which is to say, for most of applied quantitative science — topologyR offers something rare: a way to check, before you model, whether the smoothness you are about to assume is a property of your data or a story you are telling yourself.

    Links and Credits

    The package is authored by José Mauricio Gómez Julián and released under the MIT licence. It requires R ≥ 4.0.0 with Rcpp and ggplot2. The companion paper, “Bitopological Spaces from Directed Graphs: Extending the Nada Construction to Capture Temporal Irreversibility,” develops the full mathematical theory, including the central theorem, the Alexandrov sublayer, specialization preorder, pairwise connectedness, polynomial-time algorithms, and the Lean 4 formalisation.

    If you use topologyR in your research, please cite the repository release.