Espartaco

“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: Análisis

  • 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.

    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.

  • Discovering the Equations Behind Your Data: A Look at EmpiricalDynamics library in R

    Discovering the Equations Behind Your Data: A Look at EmpiricalDynamics library in R

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

    You have a time series — quarterly GDP, a population count, a sensor reading, the price of something. You suspect a differential equation is governing it, but you don’t know which one. You could guess a form and fit it. Or you could hand the data to an algorithm that searches the space of possible equations and hands you back the law, written in symbols, with a measure of how well it fits and how complex it is.

    That second option is what EmpiricalDynamics does. It’s an R package — with a high-performance Julia backend — for discovering differential and stochastic differential equations directly from empirical time series. In this post I’ll walk through what it is, why it exists, and how it works under the hood. For full mathematical detail and the recovery-test suite, the project Wiki is the authoritative reference; this post is the map.


    The problem it solves

    There’s a field called equation discovery (or “symbolic regression for dynamics”). The idea: instead of assuming dy/dt = α + βy and estimating α, β, you let a genetic algorithm breed candidate expressions — combining variables with +, , ×, ÷, sin, exp, and so on — and you keep the ones that fit well per unit of complexity. The output is a Pareto front of equations trading accuracy against simplicity, from which you pick the one that represents a real law rather than memorised noise.

    EmpiricalDynamics, written by José Mauricio Gómez Julián (v0.1.5, MIT licence), wraps this idea into a complete, rigorous workflow aimed squarely at researchers working with observational data — economists, physicists, epidemiologists, anyone studying dynamical systems they didn’t generate in a lab. The motivation is honest: real-world data is noisy, gappy, and rarely matches a textbook form, so the toolkit has to be robust to all three.


    The architecture: a six-step pipeline

    The package is structured as a pipeline you can read top to bottom. Each step is one or more exported R functions; the heavy symbolic search optionally hands off to Julia.

    1. Preprocessing — numerical differentiation. To discover dZ/dt = f(Z, X), you first need dZ/dt. Differentiating noisy data is famously treacherous, so the package offers several methods. The flagship is Total Variation Regularization (TVR), which solves a convex optimisation balancing reconstruction fidelity against a penalty on jumps in the derivative — promoting piecewise-smooth estimates that tolerate trends and structural breaks. It’s backed by a cascading solver chain (CLARABEL → SCS → OSQP) with internal rescaling so tolerances behave regardless of data scale. Cheaper alternatives — Savitzky-Golay, smoothing splines, finite differences, spectral — are there for cleaner or periodic data, and suggest_differentiation_method() will recommend one based on your data’s characteristics.
    2. Exploration. Before fitting anything, explore_dynamics() and a family of phase-portrait, bivariate, and 3-D surface plots let you eyeball the functional structure — does dZ/dt look linear in Z, quadratic, oscillatory?
    3. Equation discovery. This is the heart. symbolic_search() runs the genetic algorithm and returns a Pareto front of candidate equations. You pick one with select_equation() using AIC, BIC, MDL, or a “knee” heuristic. If you already have a theoretical form (say, a Solow growth model or a logistic equation), fit_specified_equation() is more reliable — it fits your form with Levenberg-Marquardt nonlinear least squares rather than searching blindly. The search has three backends: "r_genetic" (pure R), "r_exhaustive" (small problems), and "julia" (the industrial path).
    4. Residual analysis and SDE construction. Here’s where the package earns its keep for serious work. Once you have a drift function , the residuals tell you whether there’s leftover stochastic structure. A diagnostic suite runs Ljung-Box (autocorrelation), ARCH-LM (conditional heteroscedasticity), Breusch-Pagan, Jarque-Bera, and a runs test. If the residuals carry structure, you’re not done — you have a stochastic differential equation dZ = f dt + g dW, and now you need to recover the diffusion g.
    5. Validation. Cross-validation — critically, block CV for time series, not random CV, which would destroy temporal dependence. Plus trajectory simulation from the discovered SDE and qualitative-behaviour checks (fixed points, stability, bifurcations) to confirm the equation reproduces the dynamics you actually observe.
    6. Output. LaTeX equations, coefficient tables, model-comparison tables, and full markdown/HTML reports — publication-ready.

    The Julia backend

    The R side handles statistics, diagnostics, and orchestration. The Julia side — inst/julia/symbolic_backend.jl, calling SymbolicRegression.jl — handles the expensive evolutionary search, multi-threaded across CPU cores. The Julia code defines a ScientificSearchConfig struct with the usual knobs (population size, iterations, parsimony penalty) and a notable extra: automatic detection of physical constants — π, e, φ, g, c, h, k_B — so that when your data is generated by π·sin(t), the search can return π·sin(t) rather than 3.14159·sin(t). The README reports recovering π and e to 10⁻⁸ precision from noisy data.

    The R↔Julia bridge is JuliaCall (with JuliaConnectoR as an alternative). Recent releases (see NEWS.md) fixed real bugs in this glue — a wrong hard-coded UUID that made setup_julia_backend() falsely report the backend as unavailable, a single-predictor TypeError, and a Hall-of-Fame extraction ParseError. The lesson is mundane but important: the engine was fine; the R-side marshalling had drifted out of sync with current SymbolicRegression.jl. It now works end-to-end.


    The one idea worth understanding deeply

    If you take away a single technical point, take this: the better TVR estimates the drift, the more it destroys the information needed to recover the diffusion.

    TVR works by smoothing the derivative — penalising total variation, which suppresses high-frequency components. But in an SDE, the diffusion coefficient g lives in those high-frequency residuals. So as you crank up TVR’s regularisation to get a clean drift, the residual-based diffusion estimate collapses. The Wiki documents this starkly: across solver generations, residual-based diffusion R² fell from 0.591 to 0.005 as the solver improved.

    EmpiricalDynamics resolves the tension by estimating diffusion from quadratic variation(ΔZ)²/Δt ≈ g²(X) — using the raw increments directly, bypassing TVR entirely. With this method the diffusion R² jumps to 0.985, regardless of how aggressively TVR smoothed the drift. That’s the kind of design decision that separates a toy from a toolkit: it acknowledges a fundamental statistical trade-off and routes around it.


    Does it actually work?

    The Wiki’s recovery test suite is the evidence, and it’s worth reading in full. The epistemological framing is unusually careful: a recovery test generates synthetic data from a known ground-truth equation, hands it to the algorithm “blind,” and checks whether the algorithm rediscovers the law. The Wiki is explicit that this is a necessary but not sufficient condition for genuine causal inference — spurious correlations and unobserved confounders remain possible, and there’s an irreducible simulation-to-reality gap. That honesty matters.

    The headline numbers:

    • Lorenz attractor (chaotic ODE): average R² ≈ 0.937 across the three Lorenz equations, with the TVR solver reporting 30/30 optimal convergences.
    • SDE recovery with drift 10·sin(X) − 2.5·Z³ and diffusion 0.10 + 0.06·|X|: drift R² = 0.841, diffusion R² = 0.985 (via quadratic variation), under a challenging signal-to-noise ratio of 0.29.

    These are strong results, and the Wiki is transparent about where the limits are — diffusion R² above 0.5 is often the practical ceiling when SNR < 1, and the symbolic search is stochastic, so run-to-run variability is real.


    Should you use it?

    If you’re a researcher with time series data and a hypothesis that a differential equation is at work — economic growth, epidemiological spread, predator-prey dynamics, interest-rate feedback — EmpiricalDynamics gives you a principled, end-to-end pipeline rather than a bag of tricks. It won’t hand you causality on a plate (nothing will), but it will tell you, defensibly, what functional form the data is consistent with, how confident you should be, and what the residuals imply about stochastic structure.

    The package is young (v0.1.5, first released late 2025) and the roadmap is ambitious — GPU acceleration, neural-guided search, Bayesian SDEs via Stan. But the core is solid, the methodology is rigorous, and the documentation — both the README and the Wiki — is unusually thorough. If equation discovery is relevant to your work, this is a project worth watching and, more to the point, worth using.

    Links: GitHub repository · Wiki (full theoretical foundations, recovery tests, and API reference)