AutoCause: A Python framework that automates expert decisions in environmental time-series causal discovery

# AutoCause: A Python framework that automates expert decisions in environmental time-series causal discovery
Source: [https://arxiv.org/html/2608.00198](https://arxiv.org/html/2608.00198)
\[orcid=0000\-0002\-5163\-413X\]\\cormark\[1\]\\creditConceptualization, Methodology, Software, Validation, Formal analysis, Writing – original draft

1\]organization=ISR\-Lisbon, Instituto Superior Técnico, addressline=Av\. Rovisco Pais 1, city=Lisbon, postcode=1049\-001, country=Portugal

\[orcid=0000\-0003\-4277\-0292\]\\creditConceptualization, Supervision, Writing – review and editing

2\]organization=Digital Scholarship at Oxford, University of Oxford, city=Oxford, postcode=OX1 3BG, country=United Kingdom

\[orcid=0000\-0002\-2564\-8116\]\\creditSupervision, Writing – review and editing

3\]organization=Jet Propulsion Lab\., Caltech, city=Pasadena, postcode=CA 91109, state=California, country=USA

\[orcid=0000\-0002\-5655\-9562\]\\creditSupervision, Writing – review and editing

\\cortext

\[1\]Corresponding author

Miguel Arana\-Catania\[David R\. Ardila\[Rodrigo Ventura

###### Abstract

Environmental time\-series causal discovery requires expert decisions about method choice, conditional\-independence tests, lag horizons, sample\-size adequacy, multiple\-testing control, and evidence interpretation\. Applied inconsistently across datasets, these choices yield graphs that cannot be compared, reproduced, or audited\. We present AutoCause, an open\-source Python workflow that records each decision, derives defaults from an extended causal\-audit module, and admits domain\-informed overrides\. The workflow wraps four established causal\-discovery methods from three families, adds non\-causal reference models, and grades links by method\-count support\. On 145 datasets from DGP\-Atlas, TimeGraph, and a topology\-derived CausalRivers reference, the methods recover complementary parts of the reference graphs\. Majority\-supported links are more precise than single\-method links on the synthetic benchmarks but not against river topology\. AutoCause converts inconsistent expert practice into an auditable, repeatable analysis; causal interpretation remains with the analyst\. Available at[https://github\.com/marcoruizrueda/autocause](https://github.com/marcoruizrueda/autocause)\.

###### keywords:

causal discovery\\septime series\\sepenvironmental applications\\sepmulti\-method consensus\\sepassumption diagnostics\\sepopen\-source framework

## 1Introduction

Environmental models are commonly developed to answer questions that are causal in nature: which variables drive changes in an environmental response, which interactions are direct rather than mediated by other processes, which predictors are merely proxies for shared forcing, and which parts of a model are robust enough to support interpretation or management decisions\.

In most of these systems, controlled experiments are impractical and the analyst must work with whatever multivariate record the monitoring infrastructure provides\. A 10\-year discharge series from a river network, a 3\-year half\-hourly flux record from a peatland, or a decadal satellite atmospheric column all share the same interpretive problem: the variables are autocorrelated, seasonally forced, nonlinearly coupled, and influenced by drivers that no sensor measured\. Under these conditions, simple association \(Pearson or Spearman correlation\) can be informative but insufficient: two variables may be highly correlated because both respond to the annual cycle, because one mediates the effect of another, or because a third unobserved process drives both\.

Several families of causal\-discovery methods have been developed to address this conflation, reviewed recently bybrouillard2025grounding,niu2024,runge2023causal, andrunge2019inferring\. Each family rests on different mathematical assumptions about functional form, temporal memory, stationarity, and causal sufficiency\(assaad2022survey;moraffah2021causal;glymour2019review\)\. Constraint\-based methods such as PCMCI\+ require that a partial\-correlation or mutual\-information test can detect the conditional dependence, regression methods such as VARLiNGAM need non\-Gaussian residuals for full identifiability, and information\-theoretic approaches need large effective samples for their nearest\-neighbour estimators\. Violating any of these conditions produces a degraded or empty graph, and published comparisons find no single method that performs best across all data regimes\(assaad2022survey;nowack2020causal\)\. To support the method decision,runge2023causalorganize method choice through a questionnaire based on observable data characteristics, including linearity, stationarity, lag structure, and possible latent confounding\. The questionnaire structures the decision, but the practitioner must still evaluate each characteristic and translate the result into a software configuration\.

The questionnaire still leaves several steps to the analyst\. A Tigramite workflow\(runge2019detecting\), for example, requires data assessment, conditional\-independence \(CI\) test selection, a lag horizon, multiple\-testing handling, and interpretation of method\-specific output\. Applying a fixed PCMCI\+ and ParCorr configuration to seasonal or nonlinear benchmark data can produce low precision or recall\(ferdous2025timegraph\)\.brouillard2025groundingidentify the gap between methodological requirements and application settings as a barrier to using causal discovery on real data\.

AutoCause was developed to close this gap by formalising the configuration and evaluation steps surrounding established causal\-discovery algorithms\. The framework invokes Tigramite for constraint\-based methods and the LiNGAM library, but its contribution lies in the surrounding decision logic rather than in a new discovery algorithm\. A single function call ingests a pandas DataFrame and records the six decisions defined above\. Pre\-discovery diagnostics from an extended causal\-audit module\(ruiz2026causalaudit\)inform method, CI\-test, lag, and sample\-size decisions together with any required preprocessing\.

Once each method has completed its analysis and the correspondingpp\-values \(where available\) have undergone false\-discovery\-rate \(FDR\) correction, AutoCause counts how many of the four evaluated methods identify each candidate link and assigns the appropriate consensus\-support tier\. Dataset\-level surrogate diagnostics are reported separately and are not incorporated into the tier assignment\. The workflow generates one CSV file for each method, containing both the raw and FDR\-corrected results, together with the consensus\-tier edge list and cross\-method comparison figures\. When a reference graph is provided, the corresponding graph\-recovery metrics are also reported\. The execution configuration is stored alongside the generated outputs, allowing the analysis to be reproduced or repeated with modified settings\. Apart from requiring a regularly sampled time axis, the interface does not assume a particular application domain\. The same workflow can therefore be applied to river discharge measured every 6 hours, eddy\-covariance fluxes sampled every 30 minutes, or satellite retrievals acquired at daily intervals\.

The benchmarks used for the present evaluation comprise 145 datasets from three collections\. DGP\-Atlas\(ruiz2026causalaudit\), an atlas of data\-generating processes \(DGPs\), contributes 97 synthetic series across 10 families of controlled assumption violations\. TimeGraph\(ferdous2025timegraph\)contributes 18 categories that probe linear, nonlinear, trend\-seasonal, and missing\-data regimes\. CausalRivers\(stein2025causalrivers\)contributes 30 five\-station subgraphs from a real Bavarian discharge network evaluated against its physical topology\. Other environmental domains, from atmospheric chemistry and marine biogeochemistry to ecological sensor arrays, remain outside the present scope and will need their own evaluation before transfer is claimed\.

##### Contributions\.

This paper makes three contributions at different levels of the causal\-discovery workflow\.

At the software level, AutoCause provides the open\-source pipeline that implements the workflow described above\. Tigramite and the LiNGAM library supply the discovery algorithms; AutoCause adds the diagnostic, configuration, and evaluation layers that those libraries leave to the user\. Every recorded default can be replaced with a domain\-informed setting, and each method’s raw output remains available for independent inspection\. Among the releases compared in Table[2](https://arxiv.org/html/2608.00198#S2.T2), none combines pre\-discovery diagnostics, adaptive conditional\-independence\-test selection, method\-count support grading, surrogate falsification, and non\-causal reference models in one workflow\. This integration is the software\-level novelty of the paper\.

At the pre\-discovery level, the paper contributes three extensions to the original causal\-audit module\(ruiz2026causalaudit\)\. VARLiNGAM enters the recommendation space for datasets whose risk scores indicate linear, non\-Gaussian, and causally sufficient conditions\. A diagnostic\-extraction defect that silenced nonlinearity\-aware and seasonality\-aware routing has been corrected\. The abstention rule has been revised so that seasonal data receive a deseasonalization recommendation rather than an unconditional rejection\.

Finally, at the evaluation level, the benchmark experiments reported here cover graph recovery, method complementarity, consensus\-support precision, and dataset\-level surrogate behavior\. DGP\-Atlas is notable because, although its 97 synthetic graphs have been publicly available since 2025, no prior publication reports causal\-discovery F1 scores on them\. A single consensus\-support analysis also spans all three collections, whose reference graphs come from a generative model, a structured benchmark, and a physical river topology\. This common scale is what exposes the dependence of majority\-support precision on the type of reference graph\. The TimeGraph analysis compares the PCMCI\+ results produced by AutoCause with those reported byferdous2025timegraph\.

## 2Background and related work

AutoCause wraps methods whose intellectual roots span vector\-autoregressive \(VAR\) econometrics, information theory, and constraint\-based graphical modeling\. To situate these methods, we adopt the four\-family taxonomy ofassaad2022surveyand draw on the identifiability analysis ofrunge2023causal, who map each family’s assumptions to the environmental conditions under which it can or cannot recover causal structure\.

### 2\.1Why causal discovery matters for environmental modeling

Environmental monitoring produces multivariate time series from stream gauges, flux towers, atmospheric instruments, and ecological sensors\(runge2019inferring;delwiche2021fluxnetch4\)\. In these settings, researchers may ask whether drought response can be attributed to soil moisture, which measured variables precede methane flux changes in a peatland\(krich2020estimating\), or whether an upstream change is followed by a downstream water\-quality response\. Causal\-discovery methods test which temporal dependencies remain after conditioning on other measured variables\(runge2023causal;runge2019inferring\)\.

Once estimated, such graphs support process identification by pairing candidate driver\-response links with estimated lags\. They can inform variable selection when direct\-driver candidates are separated from seasonal or indirect associations, and they can motivate targeted field measurements or mechanistic\-model refinement\. Environmental applications range from climate\-model evaluation and constrained projections\(nowack2020causal\)to biosphere\-atmosphere interactions\(krich2020estimating\)and ecosystem methane dynamics\(delwiche2021fluxnetch4\)\.

These benefits do not imply that observational causal discovery can prove causality on its own\. Its conclusions remain conditional on assumptions about sampling, temporal resolution, causal sufficiency, stationarity, and the adequacy of the chosen conditional independence or structural model\. The practical value in environmental modeling therefore depends less on a single algorithm and more on a workflow whose assumptions are explicit, whose configuration is recorded, and whose outputs allow comparison across methodological paradigms\. AutoCause fills this role by wrapping methods from the paradigms surveyed below within a reproducible pipeline that includes pre\-discovery diagnostics, multi\-method analysis, dataset\-level surrogate diagnostics, and consensus\-support tiering\.

### 2\.2Methods for time\-series causal discovery

assaad2022surveydistinguish four families of time\-series causal\-discovery methods, a taxonomy adopted here\. VAR\-Granger\(granger1969investigating;toda1995statistical\)and VARLiNGAM\(hyvarinen2010estimation\)belong to the regression\-based family\. Transfer entropy and conditional mutual information form the information\-theoretic family; nearest\-neighbour conditional\-mutual\-information tests become more data\-demanding as the conditioning dimension increases\(runge2018conditional\)\. Constraint\-based methods such as PC, fast causal inference \(FCI\)\(spirtes2000causation\), PCMCI\+\(runge2020discovering\), and LPCMCI\(gerhardus2020high\)prune links that lack conditional\-independence support\. Continuous\-optimization methods such as DYNOTEARS encode graph constraints in an optimization objective\(pamfil2020dynotears\)\. Beyond these four families, learned\-prior approaches provide a complementary option\(thumm2026causaltime\_prior\)\. AutoCause includes methods from three of these families because their assumptions and finite\-sample behavior differ\.

The paradigms also differ in how they orient links, a distinction that governs the orientation comparison in Appendix[A\.3](https://arxiv.org/html/2608.00198#A1.SS3)\. Regression\-based and information\-theoretic methods orient every link they report: VAR\-Granger and transfer entropy by the time precedence of the lagged predictor, and VARLiNGAM additionally orients contemporaneous links through the non\-Gaussianity of the residuals\(shimizu2006lingam\)\. Constraint\-based methods orient lagged links by time order but identify contemporaneous links only up to the Markov equivalence class, returning them undirected when no unshielded collider or orientation rule fixes the direction\(runge2020discovering\)\. A symmetric association screen such as cross\-correlation assigns no direction\. These differences are intrinsic to the paradigms rather than to any implementation, so a comparison across methods reports adjacency first and treats orientation at the resolution each paradigm supports\.

### 2\.3Existing software frameworks

Given this diversity of methods and orientation conventions, a practitioner needs software that exposes not only the algorithms but also the configuration, preprocessing, and validation steps\. Tigramite\(runge2019detecting\)implements PCMCI\+, LPCMCI, and several conditional\-independence tests\. The LiNGAM package implements discovery based on independent component analysis \(ICA\) for non\-Gaussian models\(ikeuchi2023lingam\), and TETRAD\(tetrad2018\)provides a Java graphical interface\. Several frameworks expose multiple algorithms through a common application programming interface \(API\): Salesforce CausalAI\(arpit2023causalai\)includes time\-series and tabular methods, Causal\-learn\(zheng2024causallearn\)provides Python implementations of graphical\-model algorithms, and gCastle\(zhang2021gcastle\)emphasizes optimization\-based structure learning\.

Specialized extensions address narrower workflow elements\. CausalFlow\(castri2023fpcmci\)combines feature selection with PCMCI, VCDF\(vcdf2026\)evaluates cross\-temporal\-fold consistency, and CausalNex is an open\-source Python framework for Bayesian\-network modeling that implements DYNOTEARS for structure learning\(causalnexsoftware;zheng2018dags;pamfil2020dynotears\)\. Bootstrap aggregation measures stability across repeated fits of one method\(debeire2024bagged\), and a related robust cross\-validation procedure targets stability for validating agent\-based\-models\(rcv2026\)\. CausalMGM extends mixed graphical models to time\-series constraint\-based discovery\(andrews2019causalMGM\)\. Information\-theoretic alternatives include the tEDM R package\(LYU2026102435\)and the CausationEntropy software archive\(slote2025causationentropy\)\.

Tables[1](https://arxiv.org/html/2608.00198#S2.T1)and[2](https://arxiv.org/html/2608.00198#S2.T2)report the documented language, license, activity, and workflow capabilities of the compared releases\. The feature comparison distinguishes multi\-family execution, adaptive CI\-test selection, method\-count support grading, pre\-discovery diagnostics, lag estimation, and non\-causal reference models\. Surrogate diagnostics are listed as a dataset\-level falsification capability and are not treated as part of edge\-tier assignment\. AutoCause integrates all of these workflow steps around the discovery algorithms in a single reproducible call, which is why the upper block of Table[2](https://arxiv.org/html/2608.00198#S2.T2)shows checkmarks only in its column\.

Table 1:Software availability of compared frameworks, based on cited papers and repository documentation checked on 11 July 2026\.FrameworkLanguageLicenseSourceAutoCausePythonAGPLv3\+This workTigramitePythonGPL\-3\(runge2019detecting\)CausalAIPythonBSD\-3\-Clause\(arpit2023causalai\)CausalFlowPythonGPL\-3\(castri2023fpcmci\)gCastlePythonApache 2\.0\(zhang2021gcastle\)causal\-learnPythonMIT\(zheng2024causallearn\)VCDFPythonApache 2\.0\(vcdf2026\)CausalNexPythonApache 2\.0\(causalnexsoftware\)LiNGAMPythonMIT\(ikeuchi2023lingam\)TETRADJavaLGPL\(tetrad2018\)Table 2:Capability comparison of time\-series causal\-discovery frameworks\.✓\\checkmarkmeans the feature is documented as built in,∼\\simmeans partial or manual support, and – means the feature was not found in the cited release\. The lower block lists capabilities that AutoCause does not currently provide, included so that readers can judge what the framework omits\.CapabilityAutoCause Tigramite CausalAI CausalFlow gCastle VCDF CausalNex LiNGAM TETRAD Capabilities integrated by AutoCauseMulti\-paradigm✓\\checkmark–∼\\sim–∼\\sim–––∼\\simAdaptive CI\-test✓\\checkmark––––––––Evidence tiering✓\\checkmark––––––––Falsification✓\\checkmark∼\\sim–––∼\\sim–––Pre\-discovery diag\.✓\\checkmark––––––––Correlation baseline✓\\checkmark––––––––Predictive baseline✓\\checkmark––––––––Autoτmax\\tau\_\{\\max\}✓\\checkmark––––––––Paradigm diversity✓\\checkmark––––––––TS\-native✓\\checkmark✓\\checkmark✓\\checkmark✓\\checkmark∼\\sim✓\\checkmark∼\\sim✓\\checkmark∼\\simLatent confounders✓\\checkmark✓\\checkmark––∼\\sim–––✓\\checkmarkWhere AutoCause is partial or absentGPU acceleration––✓\\checkmark–✓\\checkmark–✓\\checkmark––Neural causal methods––✓\\checkmark–✓\\checkmark––––Continuous opt\. \(NOTEARS\)––✓\\checkmark–✓\\checkmark–✓\\checkmark––Interventional data––∼\\sim–––––✓\\checkmarkCyclic / feedback––––––––∼\\simMulti\-realization \(J\-PCMCI\+\)–✓\\checkmark–––––––Bootstrap aggregation–✓\\checkmark–––––––Graphical / no\-code interface––✓\\checkmark–––––✓\\checkmark

### 2\.4Benchmarks and datasets

Table 3:Benchmarks available for time\-series causal discovery, grouped by how the reference graph is obtained\.⋆\\starmarks the three used here\.‡indicates that the reference comes from a physical or engineered constraint external to the statistical model\. Vars andTTreport the number of variables and the series length; ODE and SDE denote ordinary and stochastic differential equations\.DatasetDomainVarsTChallengeRef\.SyntheticCauseMe/C4CClimate5–40150–500Autocorrelation\(runge2020causality4climate\)TimeGraph⋆\\starGeneric4–201 000Contemp\. edges\(ferdous2025timegraph\)CausalDynamicsODE/SDEVar\.Cont\.Stochastic\(causaldynamics2025\)DGP\-Atlas⋆\\starSynth\. VAR5–8500–1kViolations\(ruiz2026causalaudit\)Semi\-syntheticCausalTimeMedicine, Air Quality,
traffic20–36600–800Nonlinear\(cheng2024causaltime\)DREAM3/4Gene net\.10–10021Short series\(marbach2010dream\)NetSimfMRI5–100200HRF filter\(smith2011netsim\)Krebs cycleBiochem\.16Var\.Reaction net\(krebs\_benchmark2025\)Real\-worldCausalRivers⋆\\star‡Hydrology4945 yrSubgraphs\(stein2025causalrivers\)CIPCaD \(TEP\)‡Industrial331\.5kFeedback\(cipcad2022\)ESS Cryogenics‡Accelerator233HoursHigh\-dim\(mogensen2024industrial\)Sachs \(single\-cell\)Biology11Var\.Mixed interv\.\(sachs2005causal\)Available benchmarks span synthetic, semi\-synthetic, and real\-world domains \(Table[3](https://arxiv.org/html/2608.00198#S2.T3)\)\. Synthetic collections provide a known generating graph\. Semi\-synthetic collections fit generators to observations and impose a known structure, so their interpretation remains conditional on the generator\. Real\-world benchmarks provide externally constructed reference structures, such as river topology for CausalRivers, a fault\-propagation graph for CIPCaD, or a documented control architecture for ESS Cryogenics\.

Each benchmark type has a different limitation\. Synthetic data lack the measurement processes of field observations\. Semi\-synthetic data inherit assumptions from their generators\. Real\-world evaluations require part of the graph to be specified independently, and that reference may omit relevant processes\. Many environmental applications lack such an external reference\(brouillard2025grounding\)\.

We evaluate on the three benchmarks marked⋆\\starin Table[3](https://arxiv.org/html/2608.00198#S2.T3)\. They were chosen to cover three distinct evaluation regimes: a fully controlled synthetic setting where every link is known \(DGP\-Atlas\), a structured synthetic benchmark with both lagged and contemporaneous ground truth that allows comparison against published outputs \(TimeGraph\), and a real environmental setting where the reference comes from river topology rather than from a generative model \(CausalRivers\)\. Below, each collection is described\.

DGP\-Atlasprovides 500 synthetic VAR\(1\) processes organized in 10 families, each isolating a specific assumption violation \(nonstationarity, structural breaks, irregular sampling, high persistence, latent confounders, seasonality, polynomial nonlinearity, non\-Gaussian noise, mixed violations, and extreme cases\)\. Ground truth is the VAR\(1\) coefficient matrix and all edges are lagged \(lag 1\), with no contemporaneous structure\.

TheTimeGraphdataset comprises 18 groups using four variables, having a lag of two time steps, and containing 1000 samples for each group\. The dataset consists of four structural families\.AAgroups have a linear structure, whileBBgroups are linear structures but with nonlinearity of polynomials \(x2x^\{2\},x3x^\{3\}\)\.CCgroups have trends and seasons on top of the linear structure, andDDgroups have missing data blocks\. The confounded variant is marked by aCsuffix \(A1C, B1C, C1C, D1C\), in which one driving variable is removed from the observed set so that its influence appears as latent confounding\.

Finally, the fullCausalRiversbenchmark provides two real\-world river\-network datasets, and we use RiversBavaria: 494 gauging stations recording water level at 15\-minute intervals over 2019–2023, converted to discharge through rating curves maintained by the Bavarian Environment Agency\. The reference graph is the directed river topology itself, where each edge represents a direct upstream\-to\-downstream station pair established from digital elevation models and field surveys\. Because this reference encodes only channel connectivity, it does not capture shared rainfall, snowmelt routing, reservoir operations, tributary inflows, or other catchment processes that can generate statistical dependence between stations\. A link that a discovery method reports but the topology omits is therefore scored as a false discovery, even though it may reflect a genuine hydrological coupling that the reference simply does not represent\.

In order to investigate different types of methodological problems within this single river network, the benchmark uses three distinct sampling strategies for subgraphs\. Sampling of randomly connected subgraphs of five stations is used to check how well the algorithms perform in complex topologies where no particular structure can be assumed\. In root\-cause chains, the longest directed path contains all five stations sequentially to check how well the algorithm distinguishes between the direct and indirect links by ignoring the latter through conditioning\. For confounder subgraphs, in which there is at least one station with two or more downstream stations, the problem tested is how the methods respond to the presence of a common forcing variable\. In all these cases, sampling is done without overlap between the rivers to ensure diversity in the sampling subgraphs\.

## 3The AutoCause framework

The range of methods discussed above still requires practitioners to make six decisions before interpreting the results: which discovery method to run, which conditional\-independence test to use, how to set the lag window, whether the available sample is adequate, how to control the false discovery rate, and how to grade the resulting evidence\. AutoCause records a default or user\-supplied value for each decision and logs the preprocessing actions that support it\.

### 3\.1Architecture and inputs/outputs

The AutoCause package is fully written in Python with a main experiment script around the single methodrun\_causal\_discovery\_workflow\. There, the user provides a pandas DataFrame object with a time\-based index and column names for all observable variables, an output directory, and, optionally, a dictionary specifying what stages/methods are to be executed\. In the absence of the dictionary, all stages will operate with their default settings, while in the case of a partial dictionary, only the selected stages/methods will be executed\. Within the function, there are seven stages of processing \(Fig\.[1](https://arxiv.org/html/2608.00198#S3.F1)\) that use inputs from the preceding stages and generate outputs in corresponding directories\.

The pre\-discovery phase \(Stage 0\) estimates a search horizon and screens sample\-size adequacy\. The lag heuristic combines an autocorrelation zero\-crossing, a user\-supplied domain bound, and a data\-availability constraint\. The causal\-audit module then maps diagnostic summaries for nonstationarity, irregular sampling, persistence, possible causal insufficiency, nonlinearity, and seasonality to risk scores in\[0,1\]\[0,1\]\. Its decision tree records a method recommendation and any indicated preprocessing\.

Stage 1 selects the conditional\-independence test\. The Ramsey regression\-equation specification error test \(RESET\)\(ramsey1969tests\)and a distance\-correlation comparison\(szekely2007measuring\)assess nonlinearity, while a distributional diagnostic distinguishes standard and robust partial\-correlation routes\. The available routes are CMIknn, ParCorr, and RobustParCorr\. Stage 2 provides a non\-causal association reference by evaluating Pearson and Spearman correlations over lags1,…,τmax1,\\ldots,\\tau\_\{\\max\}and recording the minimumpp\-value and its lag\.

Stage 3 runs four evaluated causal methods from three families: VAR\-Granger and VARLiNGAM \(regression\), transfer entropy\(schreiber2000measuring;kraskov2004estimating\)\(information\-theoretic\), and PCMCI\+ \(constraint\-based\)\. LPCMCI is available but excluded from the benchmark loops under the stated wall\-time budget\. Where compatiblepp\-values are available, the stage applies method\-specific Benjamini–Hochberg correction and records the testing family\. Stage 4 fits a random\-forest predictive reference without assigning causal meaning to its importances\.

Stage 5 records method agreement and assigns the method\-count support tiers in Table[6](https://arxiv.org/html/2608.00198#S3.T6)\. Stage 6 computes dataset\-level iterative amplitude\-adjusted Fourier transform \(IAAFT\) surrogate diagnostics\(schreiber2000iaaft\)and temporal\-stability summaries\. These diagnostics do not define the edge tiers\. Stage 7 writes configurations, method outputs, metrics, and figures\. A supplied reference graph enables graph\-recovery scoring\. Otherwise, the workflow retains the discovered links and diagnostic evidence\.

Optional stages can be selected subject to their input dependencies\. Each dataset directory separates raw and corrected method outputs, consensus\-support labels, diagnostic figures, and the executed configuration\. Table[4](https://arxiv.org/html/2608.00198#S3.T4)maps the principal artifacts to their interpretive purpose\. A complete file manifest is included in the versioned paper release, and Section[4\.1](https://arxiv.org/html/2608.00198#S4.SS1)identifies the files used for the reported results\.

![Refer to caption](https://arxiv.org/html/2608.00198v1/x1.png)Figure 1:AutoCause pipeline\. A complete call runs Stages 0–7\. Optional stages can be selected subject to their input dependencies\.
### 3\.2How to use it in practice

A practitioner does not need to invoke every pipeline stage at once\. The interpretation sequence described here separates data assessment from causal\-discovery output\. Theaudit\-onlyoption returns the assumption\-risk assessment and method recommendation before discovery begins\. The recorded diagnostics cover stationarity, sample\-size adequacy, functional form, seasonality, possible causal insufficiency, and the temporal dependence horizon\. Association and predictive references can then show what structure is present before causal methods test which links remain after conditioning and method\-specific multiple\-testing control\.

The lagged\-correlation reference in Stage 2 provides the link\-level results used in the association\-to\-causation comparison in Appendix[A\.2](https://arxiv.org/html/2608.00198#A1.SS2)\. That analysis examines how much of the detected structure remains as stronger conditioning requirements are introduced\. Because lagged correlation tests marginal association, it may recover direct links, indirect paths, and dependencies induced by shared drivers\.

The random\-forest reference in Stage 4 fits oneRandomForestRegressorfor each target variable using lagged values of the remaining variables as predictors\. Permutation importance identifies predictors whose shuffled lagged values reduce out\-of\-sample predictive performance\. This model serves as a flexible predictive reference but it is neither a causal\-discovery method nor a formal upper bound on predictability\. Agreement with a causal method indicates that the corresponding lagged predictor contains predictive information, nevertheless, it does not establish causal direction\. Both non\-causal references are excluded from consensus tiering\.

Table 4:Artifacts written per dataset byrun\_causal\_discovery\_workflowand the practitioner question each answers, grouped by pipeline phase\.Non\-causal baselinescomprise lagged correlation \(Stage 2\) and the predictive baseline \(Stage 4\)\. InCausal discovery,<m\>spans all causal methods \(Granger, transfer entropy, PCMCI\+, LPCMCI, VARLiNGAM\)\.ArtifactStagePractitioner question answeredPre\-discovery diagnostics[causal\_audit/figures/prediscovery\_summary](https://arxiv.org/html/2608.00198v1/causal_audit/figures/prediscovery_summary)0cData risks, method and CI\-test recommendation[causal\_audit/figures/stationarity\_diagnostic](https://arxiv.org/html/2608.00198v1/causal_audit/figures/stationarity_diagnostic)0cStationarity over the record[causal\_audit/figures/sample\_size\_adequacy](https://arxiv.org/html/2608.00198v1/causal_audit/figures/sample_size_adequacy)0bTeffT\_\{\\text\{eff\}\}adequacy for the selected method[causal\_audit/figures/spectral\_density](https://arxiv.org/html/2608.00198v1/causal_audit/figures/spectral_density)0cSeasonality detection and deseasonalization flag[causal\_audit/figures/correlation\_and\_lagged](https://arxiv.org/html/2608.00198v1/causal_audit/figures/correlation_and_lagged)0cPairwise association at each lag[causal\_audit/figures/assumption\_deep\_dive](https://arxiv.org/html/2608.00198v1/causal_audit/figures/assumption_deep_dive)0cViolation magnitude per assumption[causal\_audit/figures/dependency\_network](https://arxiv.org/html/2608.00198v1/causal_audit/figures/dependency_network)0cMarginal dependency structure[causal\_audit/audit\_evidence\.json](https://arxiv.org/html/2608.00198v1/causal_audit/audit_evidence.json)0cRaw diagnostic test values \(machine\-readable\)[causal\_audit/diagnostic\_tables/\*\.csv](https://arxiv.org/html/2608.00198v1/causal_audit/diagnostic_tables/*.csv)0cStationarity, persistence, irregularity, and confounding results per pair[causal\_audit/risk\_profile\.json](https://arxiv.org/html/2608.00198v1/causal_audit/risk_profile.json)0cSix calibrated risk scores \(machine\-readable\)[causal\_audit/recommendation\_policy\.json](https://arxiv.org/html/2608.00198v1/causal_audit/recommendation_policy.json)0cRecommended method, CI test, and preprocessing[sample\_size\_adequacy\.json](https://arxiv.org/html/2608.00198v1/sample_size_adequacy.json)0bTeffT\_\{\\text\{eff\}\}per method, viable methods[power\_analysis\.json](https://arxiv.org/html/2608.00198v1/power_analysis.json)0bMinimum detectable effect sizeLag estimation[tau\_max\_estimation\.json](https://arxiv.org/html/2608.00198v1/tau_max_estimation.json)0aEstimatedτmax\\tau\_\{\\max\}and bounding constraintNon\-causal baselines[method/correlation/1\-raw/results\_correlation\.csv](https://arxiv.org/html/2608.00198v1/method/correlation/1-raw/results_correlation.csv)2Lagged Pearson and Spearman with best\-lag selection[figures/correlation/correlation\_heatmap\.png](https://arxiv.org/html/2608.00198v1/figures/correlation/correlation_heatmap.png)2Comparison of Pearson, Spearman, Kendall, and distance correlation[figures/correlation/method\_comparison\.png](https://arxiv.org/html/2608.00198v1/figures/correlation/method_comparison.png)2Agreement and disagreement across four correlation measures[figures/correlation/partial\_correlation\_network\.png](https://arxiv.org/html/2608.00198v1/figures/correlation/partial_correlation_network.png)2Associations surviving control for other variables[method/predictive\_baseline/1\-raw/results\_predictive\_baseline\.csv](https://arxiv.org/html/2608.00198v1/method/predictive_baseline/1-raw/results_predictive_baseline.csv)4Predictive importances with permutation p\-values[figures/per\_method/predictive\_baseline\_graph\.svg](https://arxiv.org/html/2608.00198v1/figures/per_method/predictive_baseline_graph.svg)4Predictive importance network[figures/per\_method/predictive\_baseline\_lags\.svg](https://arxiv.org/html/2608.00198v1/figures/per_method/predictive_baseline_lags.svg)4Lag with strongest predictive importanceCausal discovery[method/<m\>/1\-raw/results\_<m\>\.csv](https://arxiv.org/html/2608.00198v1/method/%3Cm%3E/1-raw/results_%3Cm%3E.csv)3Edges detected per method, with lags[figures/diagnostics/<m\>\_fdr](https://arxiv.org/html/2608.00198v1/figures/diagnostics/%3Cm%3E_fdr)3Edges surviving FDR correction per method[figures/diagnostics/<m\>\_dag](https://arxiv.org/html/2608.00198v1/figures/diagnostics/%3Cm%3E_dag)3Consistency of directed edges per method[figures/per\_method/<m\>\_graph\.svg](https://arxiv.org/html/2608.00198v1/figures/per_method/%3Cm%3E_graph.svg)3Recovered graph per method[figures/per\_method/<m\>\_pvalues\.svg](https://arxiv.org/html/2608.00198v1/figures/per_method/%3Cm%3E_pvalues.svg)3Significance level per edge per method[figures/per\_method/<m\>\_lags\.svg](https://arxiv.org/html/2608.00198v1/figures/per_method/%3Cm%3E_lags.svg)3Detected effect lags per method[figures/comparison/method\_comparison\.svg](https://arxiv.org/html/2608.00198v1/figures/comparison/method_comparison.svg)3Edges found by multiple methods[figures/comparison/pvalue\_comparison\.svg](https://arxiv.org/html/2608.00198v1/figures/comparison/pvalue_comparison.svg)3Significance comparison across methods[figures/diagnostics/lag\_analysis](https://arxiv.org/html/2608.00198v1/figures/diagnostics/lag_analysis)3Cross\-method lag analysis and agreementConsensus[ensemble\_edges\.csv](https://arxiv.org/html/2608.00198v1/ensemble_edges.csv)5Edges agreed by multiple methods[consensus/2\-core/consensus\.csv](https://arxiv.org/html/2608.00198v1/consensus/2-core/consensus.csv)5≥2\\geq\\\!2\-method agreement before tieringSurrogate diagnostics and consensus support[consensus/5\-tiers/consensus\_with\_tiers\.csv](https://arxiv.org/html/2608.00198v1/consensus/5-tiers/consensus_with_tiers.csv)5Method\-count consensus with support\-tier labels[falsification\_results\.csv](https://arxiv.org/html/2608.00198v1/falsification_results.csv)6Dataset\-level IAAFT surrogate edge\-rate resultsReproducibility and evaluation[experiment\_log\.json](https://arxiv.org/html/2608.00198v1/experiment_log.json)7Executed configuration \(reproducibility record\)[graph\_recovery\_metrics\.csv](https://arxiv.org/html/2608.00198v1/graph_recovery_metrics.csv)7Method performance against ground truth##### Usage\.

Listings[1](https://arxiv.org/html/2608.00198#LST1)–[3](https://arxiv.org/html/2608.00198#LST3)illustrate three levels of interaction with the software\. The first executes the complete workflow using the recorded diagnostic and implementation defaults\. The second restricts the run to pre\-discovery diagnostics and non\-causal reference models, which supports an initial inspection of the data before causal discovery\. The third supplies an explicit configuration in which the selected methods, conditional\-independence test, andτmax\\tau\_\{\\max\}are informed by both the diagnostic output and domain knowledge\.

Listing 1:Minimum invocation of AutoCause\. The complete workflow runs using the recorded diagnostic and implementation defaults, and the resulting configuration is stored in the output log\.1importpandasaspd

2fromframework\.core\.run\_workflowimportrun\_causal\_discovery\_workflow

3

4df=pd\.read\_csv\("my\_data\.csv",index\_col=0,parse\_dates=True\)

5result=run\_causal\_discovery\_workflow\(

6data\_df=df,

7output\_dir="results/my\_experiment",

8\)

9print\(result\["consensus"\]\["tier1\_edges"\]\)

Listing 2:Exploratory invocation\. The workflow first reports the pre\-discovery diagnostics and then runs the lagged\-correlation and predictive reference models\. These outputs describe the associative and predictive structure of the data before any causal\-discovery method is selected\.1result=run\_causal\_discovery\_workflow\(

2data\_df=df,

3output\_dir="results/exploration",

4enable\_causal\_audit=True,

5method\_config=\{

6"granger":\{"enabled":False\},

7"transfer\_entropy":\{"enabled":False\},

8"pcmci":\{"enabled":False\},

9"varlingam":\{"enabled":False\},

10"lpcmci":\{"enabled":False\},

11"correlation":\{"enabled":True\},

12"predictive\_baseline":\{"enabled":True\},

13\},

14\)

15

Listing 3:Full discovery run configured after the exploratory step in Listing[2](https://arxiv.org/html/2608.00198#LST2)\. VARLiNGAM is enabled because the diagnostics indicate non\-Gaussianity, whereas LPCMCI is excluded because of its computational cost\. The practitioner suppliesτmax\\tau\_\{\\max\}using domain knowledge1result=run\_causal\_discovery\_workflow\(

2data\_df=df,

3output\_dir="results/discovery",

4tau\_max=5,

5alpha=0\.05,

6sampling\_days=0\.25,

7method\_config=\{

8

9"granger":\{"enabled":True\},

10"transfer\_entropy":\{"enabled":True\},

11"pcmci":\{"enabled":True,

12"test\_method":"parcorr",

13"allow\_contemporaneous":True\},

14"varlingam":\{"enabled":True\},

15"lpcmci":\{"enabled":False\},

16"predictive\_baseline":\{"enabled":True\},

17"correlation":\{"enabled":True\},

18\},

19enable\_consensus=True,

20enable\_causal\_audit=True,

21true\_edges=\{\("X1","Y"\),\("X2","Y"\)\},

22undirected\_eval=True,

23\)

The framework includes an examples folder containing the three benchmark scripts, a high\-performance computing \(HPC\) submission template, and the scripts to reproduce the experiments in this paper\. The release corresponding to this paper is version0\.2\.0\. Every numerical claim corresponds to a file in the released experiment results\.

### 3\.3Methodological requirements and how AutoCause addresses them

The interpretation of causal\-discovery output depends on identifying assumptions and finite\-sample requirements that cannot be verified completely from observational data\. AutoCause records diagnostic evidence for six decision areas, summarised in Table[5](https://arxiv.org/html/2608.00198#S3.T5), and reports unresolved risks rather than certifying validity\.

Table 5:Six methodological requirements that AutoCause screens before and during discovery\. None of the diagnostics listed in the right column can verify causal validity on its own\. When domain knowledge of the system is available it should override any data\-driven estimate\.RequirementTheoretical basisAutoCause assessmentCausal search horizon \(τmax\\tau\_\{\\max\}\)Bounds the tested delays and should reflect the sampling design and plausible process timescales\(runge2023causal\)Estimated from autocorrelation\-function \(ACF\) decay when no physical value is supplied, then bounded by a user limit and an effective\-sample\-size constraint \(tau\_max\)Conditional independence testTest power depends on functional form; linear tests have low power when a dependency has little linear component\(runge2018conditional\)Assessed via RESET, distance correlation, and Shapiro–Wilk diagnostics; selects ParCorr, RobustParCorr, or CMIknn accordinglyMethod identifiabilityEach method assumes specific data properties \(linearity, non\-Gaussianity, causal sufficiency\)Six assumption\-violation risks quantified by causal\-audit; maps to the method whose conditions are best satisfiedSample size adequacyFinite\-sample power decreases as the conditioning problem grows\(runge2018conditional\)Per\-methodTminT\_\{\\min\}computed from the variable countNN,τmax\\tau\_\{\\max\}, and the missing fraction; warns or falls back when insufficientMultiple\-testing correctionTesting many candidate links raises the expected false\-discovery proportion\(benjamini1995controlling\)Benjamini–Hochberg correction applied when a method returns comparable per\-link or per\-lagpp\-values; uncorrected outputs are labelled \(fdr\_method\)Multi\-method corroborationRepeated detection can prioritize links but does not establish causality\(assaad2022survey\)Method\-count consensus assigns a support tier; dataset\-level surrogate diagnostics are reported separately##### Method selection\.

Causal\-audit maps the six risk scores to a recommended method and records the rule path that produced the recommendation\. The recommendation is advisory because observational diagnostics cannot verify every identifying assumption, and users can replace it with a domain\-informed method set\. The full benchmark loops run four causal methods to measure complementarity\. These results should therefore not be interpreted as an independent validation of the deployed selector\.

##### Lag\-window selection\.

Each causal\-discovery method requires a maximum lag,τmax\\tau\_\{\\max\}, that limits how far into the past the search for candidate causes extends\. A window that is too short excludes delayed effects occurring beyond the selected horizon\. A window that is too long increases both the number of tested links and the size of the conditioning sets, which can reduce statistical power\. AutoCause therefore allows the user to defineτmax\\tau\_\{\\max\}from process knowledge whenever such information is available\. For river discharge sampled every 6 hours, for example, an analyst may chooseτmax=5\\tau\_\{\\max\}=5when a 30\-hour horizon exceeds the expected flood\-wave travel time between connected gauges\. A study of soil\-moisture responses at daily resolution may require a longer lag window\. When no physically supported value is available, AutoCause derives two fallback estimates: the first zero crossing of the autocorrelation function\(box2015time\), and the largest lag that preserves at least five effective observations per estimated coefficient\(lutkepohl2005var\)\. The smaller estimate is retained, subject to any upper bound supplied by the user\. The run log stores the candidate values, the active constraint, and the final lag, allowing the analysis to be repeated with an alternative setting\. In the experiments reported here,τmax=5\\tau\_\{\\max\}=5is fixed for all methods so that lag selection does not contribute additional variation to the benchmark comparison\.

##### Adaptive conditional\-independence test\.

The conditional\-independence test used by PCMCI\+ directly affects which forms of dependence can be detected\. ParCorr is computationally efficient and is appropriate for approximately linear conditional relationships\. It can, however, fail when the association is strongly nonlinear, as in symmetric quadratic relationships for which linear correlation may vanish\. CMIknn estimates conditional mutual information with akk\-nearest\-neighbour procedure and can represent a broader range of functional relationships, but its computational cost and sample\-size requirements are substantially higher\(runge2018conditional\)\. AutoCause selects among the available tests using two pairwise diagnostics applied after preprocessing\. The first is the Ramsey RESET test\(ramsey1969tests\), which fits a linear model of each variable from its own lags and tests whether polynomial functions of the fitted values explain additional residual variation\. A significant result indicates that the linear specification may be inadequate\. The second diagnostic compares distance correlation\(szekely2007measuring\)with Pearson correlation for each variable pair\. A comparatively large distance\-correlation value indicates dependence that is not well represented by a linear measure\. If either diagnostic is triggered, PCMCI\+ is configured with CMIknn\. Otherwise, a Shapiro–Wilk test on the residual distributions determines whether the workflow uses ParCorr or RobustParCorr, a rank\-based alternative intended for departures from Gaussianity\. The run log records the selected conditional\-independence test, the diagnostic statistics, the correspondingpp\-values, and the rule that determined the selection\. Because RESET and distance correlation assess marginal or pairwise relationships, they may not detect nonlinear dependence that appears only after conditioning on other variables\. The user may therefore override the automatic selector and configure PCMCI\+ to use CMIknn directly\. The CI\-sensitivity analysis in Table[7](https://arxiv.org/html/2608.00198#S3.T7)complements this option by rerunning PCMCI\+ with ParCorr, RobustParCorr, and CMIknn, allowing the analyst to identify links whose detection changes with the selected test\.

##### Sample\-size screening and fallback\.

A nonparametric conditional\-independence test may still return a graph when the available sample is too small for stable estimation\. AutoCause therefore assesses sample\-size adequacy before running causal discovery\. It estimates the effective sample size,TeffT\_\{\\text\{eff\}\}, after accounting for the selected lag window and missing observations, and compares the result with method\-specific operational thresholds\. The default values are 50 for ParCorr, 200 for CMIknn, 75 for VAR\-Granger, and 120 for VARLiNGAM, with further adjustment for the number of variables andτmax\\tau\_\{\\max\}\. These thresholds serve as implementation safeguards informed by published guidance and observed estimator behavior\. They should not be interpreted as universal lower bounds or as guarantees of valid inference\. WhenTeffT\_\{\\text\{eff\}\}falls below the threshold for the selected test, the workflow issues a warning and may use a less sample\-intensive alternative, such as replacing CMIknn with ParCorr\. The run log stores the applicable threshold, the estimatedTeffT\_\{\\text\{eff\}\}, and any fallback action, allowing reviewers and users to assess whether the reported analysis was supported by an adequate sample\.

##### FDR handling\.

AutoCause applies Benjamini–Hochberg correction when a wrapped method returns comparablepp\-values over a clearly defined family of hypotheses\. For VAR\-Granger, the correction is applied across the candidate variable pairs evaluated at the selected lag\. PCMCI\+ uses the multiple\-testing procedure implemented within its discovery rule\. For transfer entropy, the correction covers the set of candidate lags tested for each variable pair\. VARLiNGAM and the random\-forest reference model do not return comparable per\-edgepp\-values, so their outputs are reported without Benjamini–Hochberg adjustment\. The corresponding result files explicitly indicate that no FDR correction was available\. This design does not impose a common statistical interpretation on methods that expose different inferential quantities\. Instead, AutoCause applies multiple\-testing correction only when a method returns statistics that define a coherent family of hypotheses\. The scope of each correction is stored in the run metadata\.

##### Multi\-method consensus and support tiering\.

The benchmark evaluation includes four causal\-discovery methods drawn from three methodological families: VAR\-Granger and VARLiNGAM, transfer entropy, and PCMCI\+\. LPCMCI is available in the software but is omitted from the full benchmark loops because of its computational cost\. Each discovered link is assigned a support tier according to the number of evaluated methods that detect it, using the fixed criteria in Table[6](https://arxiv.org/html/2608.00198#S3.T6)\. Under this four\-method composition, support from any three methods necessarily spans at least two methodological families so the family\-diversity condition therefore does not impose an additional restriction beyond the majority threshold\. The reported analysis consequently evaluates consensus by method count but does not determine whether agreement across different methodological families is more informative than agreement within the same family when the number of supporting methods is held constant\.

Table 6:Consensus\-support tiers used in the evaluation\. These labels summarize agreement among the four evaluated causal\-discovery methods and do not validate causality\.TierCriterionInterpretationTier\-1Detected by at least three of the four causal methods\. With the present method set, this necessarily spans at least two methodological families\.High consensus support\. The link is repeatedly detected under different modelling assumptions, but may still reflect indirect dependence or shared forcing\.Tier\-2Detected by exactly two causal methods\.Moderate consensus support\. The link has limited corroboration and requires method\-specific and domain\-based inspection\.Tier\-3Detected by one causal method only\.Single\-method support\. The link is sensitive to one method’s assumptions and should be treated as exploratory\.Evaluated causal methods: VAR\-Granger and VARLiNGAM \(regression\-based\), transfer entropy \(information\-theoretic\), and PCMCI\+ \(constraint\-based\)\. LPCMCI is available in the software but excluded from the full benchmark loops because of computational cost\. Lagged correlation and the random\-forest model are non\-causal references\. Dataset\-level surrogate diagnostics are reported separately in Appendix[A\.4](https://arxiv.org/html/2608.00198#A1.SS4)\.

##### Per\-edge diagnostic outputs\.

In addition to the inferred graph and consensus\-tier edge list, AutoCause generates the diagnostics listed in Table[7](https://arxiv.org/html/2608.00198#S3.T7), which support interpretation when no reference graph is available\. These outputs address whether an empty graph may reflect insufficient power, whether a reported link changes with the CI test, and whether its estimated strength varies across the record\.

Table 7:Per\-edge diagnostic outputs available to the practitioner after a discovery run\. Each addresses a specific interpretability question that arises when ground truth is unavailable\. MDES denotes the minimum detectable effect size\.DiagnosticWhat it reportsPractitioner question answeredPower analysis \(MDES\)Minimum detectable partial correlation givenTeffT\_\{\\text\{eff\}\},NN,τmax\\tau\_\{\\max\}“No edges found: insufficient power, or genuine independence?”CI\-test sensitivityEdges found by ParCorr, RobustParCorr, and CMIknn separately“Does this edge survive a change of conditional\-independence test?”Method\-specific link strengthNormalized statistic per edge \(partialR2R^\{2\}, normalized transfer entropy \(TE\), partial correlation\)“How does this method’s link statistic rank within the graph?”Temporal stabilityRolling\-window statistic across temporal blocks“Does the reported link remain similar across the record?”Lag confidenceBootstrap CI on the estimated lag per edge“What is the uncertainty range on the estimated lag?”

### 3\.4Pre\-discovery diagnostics: extending causal\-audit

Causal\-audit\(ruiz2026causalaudit\)is a pre\-discovery module that operationalizes the questionnaire\-based method selector ofrunge2023causalthrough statistical diagnostics\. It covers branches that can be assessed from one stochastic time series without interventional data: stationarity, possible causal insufficiency, contemporaneous dependence, and restricted structural\-model conditions\. A high causal\-insufficiency score can route the recommendation to LPCMCI, but it does not identify a hidden variable\. Methods that require multiple related datasets, such as J\-PCMCI\+ and seqICP, approaches designed for deterministic dynamics, such as convergent cross mapping \(CCM\), and continuous\-optimization methods, such as DYNOTEARS, are not included in the current implementation\. AutoCause summarizes six data characteristics as risk scores between 0 and 1: nonstationarity, irregular sampling, persistence, possible causal insufficiency, nonlinearity, and seasonality\. Each score is obtained through a logistic mapping and is accompanied by a bootstrap uncertainty interval\. A rule\-based decision tree then uses these scores to recommend a discovery method and conditional\-independence test\. When any prespecified risk threshold is exceeded, the workflow withholds an automatic recommendation and reports an abstention instead\.

This paper extends the original causal\-audit with three changes that emerged during DGP\-Atlas validation\. First, when all six risks for nonlinearity, confounding, nonstationarity, and irregularity score below their low\-risk thresholds and the Shapiro–Wilk test rejects Gaussianity, the decision tree now recommends VARLiNGAM, whose ICA step can orient edges that ParCorr alone leaves ambiguous\. Second, a defect in the diagnostic\-extraction code was causing the nonlinearity and seasonality risk scores to return their logistic prior intercept regardless of the actual spectral\-ratio and Spearman\-minus\-Pearson divergence observed in the data; correcting this defect restores the routing that directs seasonal series toward deseasonalization and nonlinear series toward CMIknn\. Third, the original abstention rule rejected every dataset whose seasonality risk exceeded 0\.60, even when deseasonalization followed by PCMCI\+ is well\-established practice\(box2015time;runge2020discovering\); the revised rule instead records a deseasonalization recommendation and proceeds\. Because all three changes were developed with DGP\-Atlas in view, Appendix[A\.3](https://arxiv.org/html/2608.00198#A1.SS3)treats the DGP\-Atlas recommendation accuracy as an internal consistency check rather than as independent external validation\.

Figure[2](https://arxiv.org/html/2608.00198#S3.F2)presents the pre\-discovery summary generated for each dataset\. It combines the estimated risk scores, autocorrelation profile, nonlinearity diagnostics, and recommended method in a single inspection panel\. This summary supports review of the automated configuration, but it does not replace domain knowledge about the sampling process, physically plausible delays, or potentially omitted drivers\.

![Refer to caption](https://arxiv.org/html/2608.00198v1/x2.png)Figure 2:Pre\-discovery summary for a CausalRivers subgraph comprising five Bavarian gauging stations on the Donau and Abens\. The panels show the standardized discharge series; six assumption\-violation risk scores with 95% bootstrap uncertainty intervals, grouped as low \(<0\.30<0\.30\), moderate \(0\.30–0\.60\), and high \(≥0\.60\\geq 0\.60\); variable\-specific autocorrelation functions with 95% significance bands and inset estimates ofTeffT\_\{\\text\{eff\}\}; pairwise nonlinearity diagnostics; and the resulting method and preprocessing recommendation\. In the final panel,Nonedenotes abstention because at least one prespecified high\-risk threshold was exceeded\.

## 4Experimental evaluation

The evaluation uses the three benchmark collections introduced above: DGP\-Atlas, TimeGraph, and CausalRivers, addressing these three questions:

- •How does the automated workflow perform relative to fixed or published configurations on the three benchmark collections?
- •Do the evaluated causal\-discovery methods recover complementary parts of the reference graph across datasets?
- •Does majority support identify links with higher precision, and under which environmental conditions does this ordering break down?

Reference\-graph information is not provided to AutoCause during an individual run\. It is introduced only after discovery to compute evaluation metrics\. DGP\-Atlas informed development of the pre\-discovery module, while TimeGraph influenced the design of the benchmark comparison\. The experiments therefore assess the software implementation in terms of graph recovery, complementarity among methods, precision across support tiers, and dataset\-level surrogate behavior\. They do not independently validate the automatic lag estimator, every sample\-size fallback rule, or the thresholds used by the causal\-audit module\. The reported findings should consequently be interpreted as benchmark evidence for the evaluated release, not as external validation of all six automated decisions\.

### 4\.1Setup, dataset selection, and evaluation metrics

The three benchmarks examine different aspects of the workflow\. DGP\-Atlas represents controlled violations of statistical assumptions, TimeGraph provides structured synthetic failure modes, and CausalRivers supplies environmental time series evaluated against a topology\-derived river reference\. For DGP\-Atlas and TimeGraph, the known generative graph is used as ground truth\. In CausalRivers, the directed river topology serves as an external reference, but it does not constitute a complete model of discharge dynamics\. Shared rainfall, snowmelt, reservoir regulation, tributary inflows, and other omitted processes may produce dependencies that are absent from the topology\. The primary evaluation is performed on graph skeletons so each reported link is treated as an undirected adjacency between its two endpoint variables, discarding its orientation mark and its estimated lag\. With TP, FP, and FN denoting the true\-positive, false\-positive, and false\-negative adjacency counts, skeleton precision isPrecision=TPTP\+FP\\mathrm\{Precision\}=\\frac\{\\mathrm\{TP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FP\}\}\. Recall, reported as the true\-positive rate isTPR=TPTP\+FN\\mathrm\{TPR\}=\\frac\{\\mathrm\{TP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FN\}\}\. The F1 score is the harmonic mean of the two,F1=2​Precision​TPRPrecision\+TPR\\mathrm\{F1\}=\\frac\{2\\,\\mathrm\{Precision\}\\,\\mathrm\{TPR\}\}\{\\mathrm\{Precision\}\+\\mathrm\{TPR\}\}\. Thus, the structural Hamming distance isSHD=FP\+FN\\mathrm\{SHD\}=\\mathrm\{FP\}\+\\mathrm\{FN\}, and the graph\-level false discovery rate isFDR=FPTP\+FP\\mathrm\{FDR\}=\\frac\{\\mathrm\{FP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FP\}\}\. This quantity is the realized proportion of false links in the predicted graph\. It differs from the expected false discovery rate targeted by the Benjamini–Hochberg multiple\-testing procedure\. When a method predicts no links, the released evaluation records graph FDR as zero\. This convention must be interpreted together with TPR, because the precision denominator is empty in such cases\.

On the synthetic benchmarks, FP and FN are measured against the known generating graph\. In contrast, on CausalRivers, they are measured against direct upstream\-to\-downstream adjacencies in the topology\-derived reference\. The area under the receiver operating characteristic \(AUROC\) curve and the area under the precision–recall \(AUPRC\) curve are reported only for methods that supply a continuous score comparable across candidate links\. This practically translates into using PCMCI\+ \(gives calibratedpp\-values for each individual edge\), VAR\-Granger \(F\-testpp\-values\), and transfer entropy \(permute\-test basedpp\-values\)\. In the case of VARLiNGAM, we have estimates of connection strength but no associatedpp\-values; in addition, random forest \(baseline\) gives us permutation importance values which do not make sense across pairs\.

Furthermore, there are distinct identifiability conditions for lagged and contemporaneous relationships, hence the orientation problem is considered independently of the skeleton construction\. While TimeGraph has both types of relationships and CausalRivers defines upstream\-to\-downstream oriented contemporaneous relationships, the algorithms have a different capability to orient the latter relationships\. Thus, the skeleton reconstruction procedure becomes the common objective measure for all three benchmark sets, and the orientation issue is considered separately in Appendix[A\.3](https://arxiv.org/html/2608.00198#A1.SS3)\.

##### Parameters and design choices\.

The main 145\-dataset evaluation usesα=0\.05\\alpha=0\.05andτmax=5\\tau\_\{\\max\}=5\. Benjamini–Hochberg correction is applied wherever per\-edgepp\-values are available: across candidate pairs for VAR\-Granger, within the PCMCI\+ procedure, and across candidate lags for each transfer\-entropy pair\. The lagged\-correlation reference instead selects the smallestpp\-value across lags and two correlation tests and compares it directly withα=0\.05\\alpha=0\.05, without an additional correction across pairs\. It is therefore the least conservative reference and is used only as an association screen\. VARLiNGAM and the random\-forest reference do not return per\-edgepp\-values and are reported without this correction\. No benchmark\-specific hyperparameter search was conducted; the remaining choices follow from the pre\-discovery diagnostics or stated operational constraints\. Five further design choices are fixed across experiments:

First,τmax=5\\tau\_\{\\max\}=5is held fixed across all three benchmarks rather than re\-estimated per dataset\. The framework’s automatic estimator \(Section[3\.3](https://arxiv.org/html/2608.00198#S3.SS3)\) is bypassed so that every method searches the same lag window\. The chosen value covers the true lag of two on TimeGraph with margin and matches the estimator’s recommendation on a representative DGP\-Atlas pre\-discovery audit\.

Second, the adaptive conditional\-independence test selector operates through the pre\-discovery module\. The evaluated release uses fixed operational thresholds for NonlinearityRisk, SeasonalityRisk, and IrregularityRisk to determine the conditional\-independence test and any preprocessing action\. Because DGP\-Atlas informed the development of this module, the corresponding thresholds are implementation defaults rather than universally validated cutoffs\. AutoCause stores their values in the run configuration, allowing users to replace them when domain knowledge or a sensitivity analysis justifies an alternative setting\.

Third, the DGP\-Atlas configuration includes ten dataset identifiers per family\. The pre\-discovery module flags the high\-persistence family for additional stationarity assessment, but those datasets remain in the stress\-test evaluation\. Some irregular\-sampling datasets also fall below the operational effective\-sample threshold; they are retained so that the reported results do not depend on the screening stage accepting the dataset\. Three nonlinear datasets are excluded because their explosive trajectories violate the bounded\-variance rule implemented in the release, leaving 97 DGP\-Atlas datasets for evaluation\.

Fourth, LPCMCI is excluded from the benchmark loops because the stated trial configuration exceeded the two\-hour wall\-time limit\. Fifth, we resample CausalRivers from 15\-minute to 6\-hour resolution and restrict the analysis to calendar year 2021, which keeps runtime tractable while providing approximately 1460 observations per station\. The lag horizon is fixed atτmax=5\\tau\_\{\\max\}=5, corresponding to a 30\-hour search window that exceeds the flood\-wave travel time between directly connected Alpine\-foreland stations\. This value is an operational benchmark setting chosen for cross\-method comparability and is not claimed to cover every hydrological delay in the network\. From each of the three topology classes described in Section[2\.4](https://arxiv.org/html/2608.00198#S2.SS4)we draw 10 subgraphs of five stations, giving 30 subgraphs total\.

### 4\.2Causal discovery results

Method performance appears in Figure[3](https://arxiv.org/html/2608.00198#S4.F3)and Table[8](https://arxiv.org/html/2608.00198#S4.T8)\. VARLiNGAM has the highest aggregate DGP\-Atlas F1\. On TimeGraph, the evaluated causal methods have lower recall in several nonlinear categories\. CausalRivers causal\-method point estimates lie between 0\.57 and 0\.60 and their intervals overlap\. No paired test is used to claim equality\. The random\-forest model measures lagged predictability and remains a non\-causal reference\.

![Refer to caption](https://arxiv.org/html/2608.00198v1/x3.png)Figure 3:Mean F1 score per method per benchmark with dataset\-level bootstrap 95% intervals from 2000 resamples\. Lagged correlation and the random\-forest predictive baseline are non\-causal references\. Source:make\_figures\.pyand the per\-dataset metric files specified in Section[4\.1](https://arxiv.org/html/2608.00198#S4.SS1)\.Table 8:Mean F1 score per method per benchmark using unordered skeleton pairs\. Bold marks the highest point estimate per benchmark, with ties bolded\. Non\-causal references are not ranked as causal methods\. LPCMCI was excluded from the benchmark because the trial configuration exceeded the predefined wall\-time limit\.MethodDGP\-AtlasTimeGraphCausalRiversCausal methodsVAR\-Granger0\.8640\.4820\.571VARLiNGAM0\.8890\.3960\.597Transfer entropy0\.7560\.4140\.602PCMCI\+ \(adaptive CI\)0\.8140\.4860\.600Non\-causal baselines referencesLagged correlation0\.8010\.6080\.602Random Forest0\.7880\.6960\.602Because the 10 DGP\-Atlas families represent different assumption violations, the aggregate F1 values in Table[8](https://arxiv.org/html/2608.00198#S4.T8)conceal substantial variation across regimes\. The disaggregated result can be seen in Figure[4](https://arxiv.org/html/2608.00198#S4.F4), which reports the family\-level results for the causal methods and the two non\-causal reference models\. VARLiNGAM obtains the highest mean F1 in six families: F1, F3, F4, F6, F8, and F10\. Its clearest advantage occurs on F8, which introduces non\-Gaussian noise, where it reachesF1=0\.93\\mathrm\{F1\}=0\.93\. This result is consistent with the ICA\-based identifiability assumptions of LiNGAM\(hyvarinen2010estimation;shimizu2006lingam\)\. PCMCI\+ reachesF1=0\.87\\mathrm\{F1\}=0\.87on the clean VAR family F1 andF1=0\.91\\mathrm\{F1\}=0\.91on F8\.

Preprocessing has a marked effect on the irregular\-sampling families\. On F3, linear interpolation before discovery raises PCMCI\+ fromF1=0\.00\\mathrm\{F1\}=0\.00toF1=0\.78\\mathrm\{F1\}=0\.78, while VARLiNGAM reachesF1=0\.91\\mathrm\{F1\}=0\.91under the same preprocessing\. On F9, which combines irregular sampling with additional violations, PCMCI\+ increases fromF1=0\.26\\mathrm\{F1\}=0\.26without interpolation toF1=0\.86\\mathrm\{F1\}=0\.86after the audit\-triggered step\. The lagged\-correlation baseline remains competitive, with F1 values between 0\.70 and 0\.93, but it does not lead any family\. F4 provides a separate stress test for persistence: although the pre\-discovery module flags the family for stationarity risk, both VARLiNGAM and PCMCI\+ retain mean F1 scores above 0\.90\.

![Refer to caption](https://arxiv.org/html/2608.00198v1/x4.png)Figure 4:Mean F1 score for each method across the 10 DGP\-Atlas families, including the lagged\-correlation and predictive reference models\. The preprocessing stage applies interpolation to F3 and F9 and removes seasonality from F6\. VARLiNGAM achieves the highest mean F1 in six families, while the lagged\-correlation baseline remains slightly below the best\-performing causal method in each family\.TimeGraph and CausalRivers break this pattern into per\-category and per\-topology detail \(Figures[6](https://arxiv.org/html/2608.00198#S4.F6)and[7](https://arxiv.org/html/2608.00198#S4.F7)\)\. On the linear TimeGraph categories A1 and A2, PCMCI\+ reaches F1 = 1\.00 \(TPR = 1\.00, FDR = 0\.00, SHD = 0\), matching the published baseline\. The polynomial\-nonlinear B categories tell a different story: ParCorr has no power against cubic and quadratic dependence, so PCMCI\+ and transfer entropy drop to near zero\. VAR\-Granger still detects the linear term embedded in the structural equations \(B1, F1≈\\approx0\.67\)\. The correlation and predictive baselines hold F1≈\\approx0\.67–0\.80 by capturing marginal association without conditioning\. Figure[5](https://arxiv.org/html/2608.00198#S4.F5)illustrates the preprocessing effect on one trend\-seasonal TimeGraph case and one irregularly sampled DGP\-Atlas case; Appendix[A\.3](https://arxiv.org/html/2608.00198#A1.SS3)provides the orientation analysis\.

![Refer to caption](https://arxiv.org/html/2608.00198v1/x5.png)Figure 5:Effect of audit\-driven preprocessing on two cases\. Both rows run PCMCI\+ \(ParCorr, Benjamini–Hochberg FDR\); the two right\-most columns differ only in whether the audit\-recommended preprocessing is applied\. Top row, TimeGraph C1 \(trend\-seasonal\): without deseasonalization three of four skeleton edges are recovered; with it all four are recovered, with one contemporaneous link oriented against the equivalence class\. Bottom row, DGP\-Atlas F3 \(irregular sampling, 27% missing\): without imputation no edges are recovered; with audit\-triggered interpolation the skeleton is exact \(TP = 6, FP = 0\) and all lagged links are oriented\. Green: correct skeleton and direction\. Blueo \-\- o: contemporaneous link left undirected within the Markov equivalence class\. Purple: recovered but oriented against the class\. Grey dotted: missed\. Orange: false positive\.The CausalRivers study evaluates a different river network, subgraph sample, and tuning protocol\(stein2025causalrivers\)\. Its published values are therefore contextual rather than direct baselines for the fixed\-configuration RiversBavaria evaluation reported here\.

The topology\-specific results \(Figure[7](https://arxiv.org/html/2608.00198#S4.F7)\) vary across methods and across the three subgraph classes defined in Section[2\.4](https://arxiv.org/html/2608.00198#S2.SS4): randomly connected subgraphs, root\-cause chains, and confounder subgraphs\. PCMCI\+ yields the lowest observed FDR on the evaluated root\-cause and confounder subgraphs\. This pattern is consistent with conditional\-independence testing removing some indirect associations, although the aggregate metrics do not isolate that mechanism\(runge2020discovering\)\. Transfer entropy attains the highest F1 point estimate on the confounder subgraphs\. Determining whether this advantage reflects nonlinear dependence induced by shared upstream drivers would require a path\-level analysis, which is outside the present evaluation\.

The connection between topologically adjacent stations creates strong lagged dependence in the discharge signal, explaining the high recall measured against the topology\-derived reference\. Predicted links outside that reference may arise from indirect propagation, unmeasured meteorological forcing, reservoir regulation, tributary inflows, or other catchment processes\. Because these mechanisms were not evaluated directly, they should be treated as possible explanations rather than identified causes of the observed error pattern\. Across the evaluated methods and topology classes, graph FDR ranges from 0\.49 to 0\.60\.

All CausalRivers results reported here are based on 30 subgraphs of five stations, discharge records from calendar year 2021, and a 6\-hour sampling interval\. The sample represents a limited rainfall–nival regime in Bavaria and differs from the flatter East German network examined bystein2025causalrivers\. Regulated rivers, snowmelt\-dominated catchments, or networks evaluated at different temporal resolutions may produce different method rankings\. Similar performance across the three selected topology classes therefore does not establish transfer to other river systems, climates, lag windows, or monitoring designs\.

The following analyses examine method complementarity \(Section[4\.3](https://arxiv.org/html/2608.00198#S4.SS3)\) and test whether majority support increases precision \(Section[4\.4](https://arxiv.org/html/2608.00198#S4.SS4)\)\. Figure[8](https://arxiv.org/html/2608.00198#S4.F8)illustrates the support\-tier rule on one DGP\-Atlas example\.

![Refer to caption](https://arxiv.org/html/2608.00198v1/x6.png)Figure 6:F1 score for each method across the 18 TimeGraph categories, evaluated on unordered skeleton pairs\. Lagged correlation and the predictive model are included as non\-causal reference methods\. The figure was generated withmake\_figures\.pyfrom the TimeGraph per\-family result files\.![Refer to caption](https://arxiv.org/html/2608.00198v1/x7.png)Figure 7:Skeleton TPR \(solid\) and graph FDR \(hatched\) for CausalRivers, grouped by subgraph topology, with bootstrap 95% intervals across subgraphs\. The three classes are the CausalRivers sampling strategies: randomly connected five\-station subgraphs, root\-cause chains in which the longest directed path contains all five stations, and confounder structures in which at least one station has two or more downstream stations\. PCMCI\+ has the lowest FDR point estimate on both the root\-cause and confounder subgraphs, with 0\.49 in each case\.![Refer to caption](https://arxiv.org/html/2608.00198v1/x8.png)Figure 8:Method\-count support for one DGP\-Atlas instance from the non\-Gaussian F8 family \(dgp\_002\), containing six variables and nine reference connections\. All results are evaluated on the unordered skeleton\. Panel \(a\) shows the reference graph\. Panel \(b\) shows Tier\-1 links detected by at least three methods, recovering 8 of 9 reference links with precision 0\.89\. Panel \(c\) shows Tier\-2 links detected by two methods, of which 1 of 2 matches the reference\. Panel \(d\) shows Tier\-3 links detected by one method, none of which matches the reference\. Orange edges are absent from the reference graph, and grey dashed edges denote missed reference connections\. The figure was generated withreproduce\_fig8\_tiers\.py\.
### 4\.3Multi\-method evaluation and comparison with reported results

The TimeGraph analysis compares the PCMCI\+ results produced by AutoCause with the published values reported byferdous2025timegraphon the same datasets\. The two evaluations differ in their preprocessing steps, graph\-matching conventions, and multiple\-testing procedures\. The comparison therefore characterizes the behavior of the two pipelines rather than isolating the effect of any single design choice\. It is a diagnostic benchmark, not a component ablation, and is separate from the 18\-category main evaluation\. Its ten rows run under a separate released configuration rather than the main\-evaluation defaults of Section[4\.1](https://arxiv.org/html/2608.00198#S4.SS1)\. Within this diagnostic run, the categories, sample length, variable count, lag setting, and ParCorr test are aligned; AutoCause retains its own input handling and multiple\-testing procedure\. SHD values are excluded because the two pipelines use incompatible graph conventions\. Attributing any performance difference to a single stage would require an isolated component experiment\.

Table 9:Published PCMCI\+ values\(ferdous2025timegraph\)and AutoCause PCMCI\+ values on ten reported TimeGraph rows\. Each column retains its pipeline’s graph\-matching convention; the comparison is diagnostic and does not estimate a component effect\. FDR is recorded as zero for an empty prediction under the convention defined in Section[4\.1](https://arxiv.org/html/2608.00198#S4.SS1)\.Published PCMCI\+AutoCause PCMCI\+CategoryNoiseTPRFDRTPRFDRA1 \(linear\)Gaussian1\.000\.001\.000\.00A1 \(linear\)Student\-t0\.670\.331\.000\.00B1 \(poly\. nonlinear\)Gaussian0\.001\.000\.000\.00B1 \(poly\. nonlinear\)Student\-t0\.330\.931\.000\.33C1 \(trend\-seasonal\)Gaussian0\.001\.001\.000\.33A1C \(linear \+ confounder\)Gaussian0\.670\.501\.000\.20A1C \(linear \+ confounder\)Student\-t0\.670\.501\.000\.00B1C \(poly\. nonlin\. \+ conf\.\)Gaussian0\.220\.000\.000\.00B1C \(poly\. nonlin\. \+ conf\.\)Student\-t0\.330\.941\.000\.33C1C \(trend\-seas\. \+ conf\.\)Gaussian0\.330\.790\.500\.33Mean0\.420\.600\.750\.15The pipelines agree on the clean A1 Gaussian row and differ on several seasonal, heavy\-tailed, and confounded rows\. Both have low TPR on at least one Gaussian polynomial\-nonlinear row\. Because graph conventions and processing pipelines differ, the table describes reported behavior rather than a matched effect estimate\. Attributing a difference to one decision rule would require a component ablation\.

Table[10](https://arxiv.org/html/2608.00198#S4.T10)compares PCMCI\+ with a per\-dataset oracle that selects the highest\-F1 result among the four evaluated causal methods and configurations\. The oracle is not deployable because it uses the reference graph\. The oracle achieves a higher F1 score than PCMCI\+ on each of the three benchmark collections, while the identity of the best\-performing method varies across datasets\. This establishes complementarity within the evaluated set, not the performance of a deployable selector\.

Table 10:PCMCI\+ versus a per\-dataset oracle that chooses the highest\-F1 result among the four evaluated causal methods and configurations using the reference graph\. The oracle bounds selectors restricted to this evaluated set; it is not deployable\. The final column reports the most frequent sole winner\. Skeleton F1 is used throughout\.BenchmarkPCMCI\+ aloneMulti\-method oracleΔ\\DeltaF1Most frequent best methodDGP\-Atlas0\.8140\.918\+0\.104\+0\.104VARLiNGAM \(27/97\)TimeGraph0\.4860\.639\+0\.153\+0\.153PCMCI\+ \(6/18\)CausalRivers0\.6000\.648\+0\.048\+0\.048PCMCI\+ \(14/30\)A separate question concerns the vote threshold used to turn those multi\-method outputs into a prioritized edge list\. Table[11](https://arxiv.org/html/2608.00198#S4.T11)compares the majority\-of\-three rule used for Tier\-1 with looser and stricter vote thresholds\. In the current four\-method set, the majority\-of\-three and the stated family\-aware rule are mathematically identical because any three methods necessarily span at least two methodological families\. The reported precision gain therefore comes from vote count, not from an independently demonstrated family\-diversity effect\. Requiring any two methods admits more links at lower precision, while unanimity admits fewer links at higher precision\. Whether cross\-family agreement carries information beyond the vote count alone requires an equal\-vote comparison between same\-family and cross\-family pairs, not attempted here\.

Table 11:Vote\-threshold precision on DGP\-Atlas and TimeGraph\. Rules admit a link by agreement among the four causal methods: majority\-3 \(the Tier\-1 rule\), any\-two, and unanimous\. The family\-aware and majority\-3 rules are identical for the present method set\.nnis the number of admitted edges, so stricter rules admit fewer at higher precision\. Bold marks the majority\-of\-three rule used for Tier\-1\.RuleDGP\-Atlas prec\.DGP\-AtlasnnTimeGraph prec\.TimeGraphnnFamily\-aware majority \(≥\\geq3 methods\)0\.84812910\.69449Majority\-3 \(≥\\geq3 methods\)0\.84812910\.69449Any two \(≥\\geq2 methods\)0\.77815000\.56280Unanimous \(all 4 causal methods\)0\.9058770\.70617
### 4\.4Consensus\-support precision evaluation

Figure[8](https://arxiv.org/html/2608.00198#S4.F8)illustrates method\-count consensus on one DGP\-Atlas instance\. Across DGP\-Atlas, Tier\-1 precision is 0\.848, compared with 0\.344 for Tier\-2 and 0\.126 for Tier\-3\. TimeGraph shows the same ordering: 0\.694, 0\.355, and 0\.300\. Majority support therefore prioritizes higher\-precision links on the two synthetic benchmarks\. These values are skeleton precision, not calibrated posterior probabilities or causal\-validity scores\.

The ordering does not hold on CausalRivers\. Relative to the topology\-derived reference, Tier\-1 precision is 0\.425, while Tier\-2 and Tier\-3 precision are 0\.667 and 0\.571\. The latter values use six and seven links, respectively\. The observed result is that majority agreement does not select the highest\-precision tier in this benchmark\. Omitted shared drivers are one possible explanation, but the evaluation does not identify their contribution\.

![Refer to caption](https://arxiv.org/html/2608.00198v1/x9.png)Figure 9:Consensus\-support precision per benchmark\. On DGP\-Atlas and TimeGraph, majority\-supported Tier\-1 exceeds Tier\-2 and Tier\-3\. On CausalRivers the ordering reverses; the Tier\-2 and Tier\-3 estimates use only six and seven edges\. The dashed line shows a common reference prevalence\. Sources:recompute\_tier\_metrics\.pyandverify\_causalrivers\_tiers\.py\.

## 5Discussion

AutoCause does not relax the assumptions required by time\-series causal discovery or turn observational records into causal proof\. It records assumption diagnostics, preprocessing, lag selection, conditional\-independence testing, FDR handling, method comparison, surrogate checks, and support tiers\. The analysis path is therefore inspectable even for environmental records with autocorrelation, seasonality, missing observations, nonlinear responses, or unmeasured common drivers\. None of the recorded steps guarantees identification\.

### 5\.1From method selection to auditable causal\-discovery workflow

No causal method is the sole best result on every evaluated dataset, consistent with broader benchmark evidence\(runge2023causal;assaad2022survey\)\. PCMCI\+ with ParCorr targets approximately linear conditional relationships, VARLiNGAM uses non\-Gaussian linear structure, information\-theoretic methods offer nonparametric tests at greater sample\-size cost, and LPCMCI represents latent\-confounding ambiguity with greater computational cost\. The workflow records which assumptions and constraints informed each route rather than treating a common API as method selection\.

Table[9](https://arxiv.org/html/2608.00198#S4.T9)describes two processing pipelines rather than an isolated component effect\. They agree on the clean linear row and differ elsewhere, but their graph conventions, preprocessing, input handling, and multiple\-testing procedures prevent attribution to one decision rule\. In nonlinear categories, an empty or unstable graph should be reported with the corresponding sample\-size and test diagnostics rather than interpreted as causal absence\.

AutoCause applies Benjamini–Hochberg correction when a method exposes a compatible testing family\(benjamini1995controlling\)\. The TimeGraph comparison does not isolate the contribution of Benjamini–Hochberg correction\. It does, however, show why the correction procedure, the associated family of hypotheses, and the convention used for empty predictions must be reported explicitly\.

### 5\.2Consensus support and its limits

The consensus\-support scheme converts agreement across methods into an ordering for further inspection\. On DGP\-Atlas and TimeGraph, links supported by three or four methods have higher precision than links reported by only one method\. These tiers summarize observed agreement across the evaluated methods; they are not posterior probabilities or guarantees of causal validity\.

CausalRivers shows that this ordering is not universal\. All methods report links outside the topology\-derived reference as well as reference adjacencies, so several methods can agree on a dependency that is absent from the direct topology, possibly because of indirect propagation or omitted common drivers\.

The tier\-precision inversion follows from this difference in reference structure\. On the synthetic benchmarks, known graphs and controlled violations make majority support a useful precision filter\. On river data, Tier\-1 marks a dependency that needs physical interpretation through topology, travel time, elevation, regulation, tributaries, catchment boundaries, and measured meteorological inputs\. Cross\-method agreement does not remove latent\-confounding or omitted\-variable limitations\(peters2017elements;spirtes2000causation\)\.

For environmental modelling, the support tier should be read together with pre\-discovery diagnostics and domain context\. Elevated seasonality, persistence, or causal\-insufficiency risk weakens the interpretation of any consensus\-supported link\. Tier\-1 identifies a priority for further analysis, not a direct causal conclusion\.

### 5\.3Falsification and the role of non\-causal baselines

IAAFT surrogates add a dataset\-level check beyond graph\-recovery metrics\. Each surrogate retains the marginal distribution and approximately preserves the power spectrum of every variable while disrupting dependence between variables\. Methods that continue to report many links under this null produce a high surrogate edge rate\. On the evaluated CausalRivers subset, PCMCI\+ has a lower surrogate edge rate than the regression\-based methods\. This result characterizes aggregate behavior under the surrogate null and does not establish the validity of any individual observed link\.

Lagged correlation and the random\-forest baseline address different questions from the causal\-discovery methods\. The former summarizes temporal association, while the latter captures predictive structure\. Neither provides evidence of causality by itself\. A gap between predictive and causal recovery can be consistent with nonlinearity, omitted drivers, insufficient sample size, or violated identifiability assumptions, but it does not identify which explanation applies\. The references therefore distinguish association and prediction from dependencies that survive a causal method’s conditioning and correction procedure\. The predictive baseline in particular requires careful reading\. It is excluded from the causal voting rule and its significance flags are not corrected in the same way as the causal methods\. Its advantage on some TimeGraph categories is a predictive reference point, not a directly comparable causal\-discovery result\. The asymmetry is diagnostically useful because it exposes structure that causal methods miss, but it should not be used to rank random forests as causal methods\.

### 5\.4Decisions that still require domain expertise

Automation does not remove the need for domain judgment\. Three decisions remain especially dependent on knowledge of the underlying system\.

The first concerns the maximum lag,τmax\\tau\_\{\\max\}\. An autocorrelation\-based estimate provides an operational starting point, but it does not establish a physically meaningful delay horizon\. Information about hydrological travel time, ecological response delay, sampling frequency, and known process memory should take precedence when available\. A lag window that is too short excludes delayed effects beyond the selected horizon\. A window that is too long increases the number of tested links, enlarges the conditioning sets, and can reduce statistical power\.

The second concerns nonlinear conditional\-independence testing\. CMIknn can detect a broader range of functional relationships than linear tests, but its performance depends strongly on effective sample size and conditioning dimension\. In the evaluated settings, some nonlinear regimes remain difficult to recover despite the use of a nonlinear test\. An empty or sparse result from CMIknn should therefore be interpreted together with the corresponding sample\-size and power diagnostics rather than as evidence that no causal dependence is present\.

The third concerns latent confounding\. No software default can establish that all relevant common drivers have been observed\. LPCMCI represents uncertainty caused by latent confounders, but its computational cost exceeded the wall\-time budget used in the benchmark evaluation\. When the pre\-discovery diagnostics indicate a high risk of causal insufficiency, stronger analysis may require additional variables, explicit domain constraints, independent validation, interventional evidence, or multiple related time series\. J\-PCMCI\+\(gunther2023jpcmci\)provides one possible extension for environmental datasets that share partially observed or latent contexts\.

### 5\.5Methodological, computational, and empirical limitations

AutoCause adds a decision and evaluation layer around established causal\-discovery libraries\. It does not replace the algorithms implemented by Tigramite or LiNGAM, and the present evaluation does not support claims of algorithmic superiority over those tools or other causal\-discovery frameworks\. The contribution lies in recording diagnostic results, applying configurable decision rules, comparing method\-specific outputs, and assigning qualified support tiers\. For this reason, the original output of each wrapped method remains available rather than being reduced to a single framework\-level score\.

The identifying assumptions of the wrapped methods remain a primary limitation\. The current implementation relies on forms of stationarity, adequate temporal resolution, and a lag window that contains the relevant delays without reducing the effective sample size excessively\. Interpolation and sample\-size screening provide operational responses to irregular sampling and missing observations, but interpolation can smooth short\-lived dynamics or introduce values not directly observed\. Likewise, causal\-sufficiency diagnostics can flag possible omitted drivers without identifying them\. False\-positive and false\-negative links therefore remain possible even when all recorded diagnostics and safeguards are satisfied\.

Edge orientation requires separate interpretation from skeleton recovery\. Temporal order determines the direction of lagged links, but some contemporaneous links cannot be uniquely oriented from observational conditional\-independence information\. Several directed graphs may therefore remain consistent with the same observed dependence structure\. Recovering the correct skeleton does not guarantee that every reported direction is correct, and orientation errors may alter the physical interpretation of the inferred process\. AutoCause preserves the orientation marks returned by each method, but stronger directional conclusions require domain\-specific timing constraints, interventional evidence, or other identifying information\.

Computational cost also restricts the current method set\. The runtime of CMIknn and LPCMCI increases with the number of variables, record length, lag horizon, and conditioning complexity\. The present consensus evaluation excludes interventional methods, feedback\-aware approaches, neural causal\-discovery models\(cheng2024cutsplus;chen2023cuts;cheng2025dycast;yang2026dycausal;zhang2025local\), and methods designed for multiple related datasets\. The automatic lag estimator, sample\-size thresholds, and recommendation rules also require further sensitivity analysis because the benchmark loops do not independently validate every decision rule\.

The empirical evaluation provides only limited evidence for transfer across environmental domains\. DGP\-Atlas and TimeGraph are synthetic benchmarks, while the CausalRivers analysis uses 30 Bavarian subgraphs of five stations, one calendar year, and one temporal aggregation\. The topology\-derived reference does not represent rainfall, snowmelt, reservoir regulation, tributary effects, or other processes that may influence discharge\. The observed rankings, error rates, and support\-tier behavior should therefore not be generalized to other catchments, environmental sensor networks, ecological records, or atmospheric datasets without additional evaluation\.

Method\-count consensus has a related limitation\. On DGP\-Atlas and TimeGraph, majority\-supported links have the highest observed precision\. On CausalRivers, that ordering reverses, showing that several methods can agree on a dependency that is absent from an incomplete reference graph or driven by a shared source of error\. The support tiers should therefore be used to prioritize links for further inspection, not as a substitute for a final causal graph\.

## 6Conclusions

AutoCause is an open\-source decision\-support workflow for configuring, running, and documenting causal\-discovery methods on multivariate environmental time series\. It records method selection, conditional\-independence\-test selection, lag\-window selection, sample\-size adequacy, FDR handling, and evidence grading\. Preprocessing actions are logged, and defaults can be overridden when physical knowledge or sensitivity analysis supports another setting\.

Across DGP\-Atlas, TimeGraph, and a topology\-derived CausalRivers reference, the evaluated methods show data\-regime\-dependent graph recovery\. Majority\-supported links are more precise than single\-method links on the synthetic benchmarks, but the ordering breaks on CausalRivers\. Consensus support must therefore be read as repeated detection rather than causal validation\.

These results support AutoCause as an auditable workflow layer, not a new discovery algorithm or an automatic source of causal truth\. Further work should validate decision thresholds independently, isolate component effects, expand environmental testing, and test whether cross\-family agreement adds information beyond vote count\.

## Appendix AAdditional analyses of causal evidence and identifiability

This appendix presents four supplementary analyses that extend the main evaluation without being necessary to interpret its central results\. They examine computational cost, the transition from marginal association to conditional causal testing, consistency between empirical performance and method\-specific identifiability conditions, and dataset\-level surrogate diagnostics\. Each analysis follows the metric definitions, graph conventions, and reference\-graph qualifications established in the main text\.

### A\.1Computational cost

Table[A\.1](https://arxiv.org/html/2608.00198#A1.T1)reports empirical wall\-clock time per DGP\-Atlas dataset on an Apple M3 Pro with 16 GB RAM\. The values combine differences in variable count, data properties, preprocessing, and selected CI test, so they should not be interpreted as a controlled scalability law\. The main operational result is that the nonlinear F7 family, which invokes CMIknn for part of the analysis, is substantially slower than the predominantly ParCorr\-based families\. A controlled scalability analysis that variesNN,TT, andτmax\\tau\_\{\\max\}independently is still required before making general complexity claims\.

Table A\.1:Mean wall\-clock time per dataset measured on an Apple M3 Pro with 16 GB of RAM\. Each family contains 10 datasets except F7, which retains 7 after excluding three explosive trajectories\. The DGP\-Atlas families contain 5–8 variables and series lengths ofT=500T=500–1000\. F7 includes runs in which the conditional\-independence test was selected automatically and CMIknn was used for part of the family; the remaining families rely predominantly on ParCorr\-based configurations\. Values were computed from the per\-dataset execution logs included in the released DGP\-Atlas results\.FamilyMean time \(s\)F1 \(clean VAR, 7 vars\)157F2 \(structural breaks, 6 vars\)183F3 \(irregular sampling, 5 vars\)140F4 \(high persistence, 5 vars\)161F5 \(latent confounders, 8 vars\)155F6 \(seasonality, 7 vars\)170F7 \(nonlinear, 8 vars\)904F8 \(non\-Gaussian, 5 vars\)161F9 \(mixed violations, 6 vars\)125F10 \(extreme cases, 6 vars\)231
### A\.2From marginal association to conditional causal testing

Non\-causal reference models are included in the software to show how much structure can be recovered without imposing causal assumptions\. In CausalRivers,stein2025causalriversreport that simple cross\-correlation and reverse\-physical baselines can attain scores comparable to those of several causal\-discovery methods\.assaad2022surveylikewise use naive approaches as reference points for quantifying the gain associated with stronger causal assumptions, whilerunge2018chaostreats lagged correlation as a minimally restrictive comparator when evaluating PCMCI\. The predictive reference extends this comparison to nonlinear structure\. Its role follows the distinction drawn byshmueli2010explainbetween explanation and prediction and the argument ofpeters2017elementsthat predictive accuracy alone does not establish causal direction\.

The edge counts expose the difference between marginal association and conditional discovery\. As shown in Table[8](https://arxiv.org/html/2608.00198#S4.T8), on DGP\-Atlas, lagged correlation has mean F1 = 0\.801, recall = 0\.988, and precision = 0\.684\. The high recall follows from testing marginal association across lags without FDR correction across pairs, which admits most true edges but also retains edges that conditional methods would remove\. PCMCI\+ has mean F1 = 0\.814 with higher precision and lower recall on the same skeleton endpoint\. PCMCI\+ has mean F1 = 0\.814 on the same skeleton endpoint\.

On TimeGraph, the lagged\-correlation reference attainsF1=0\.608\\mathrm\{F1\}=0\.608, recall=0\.722=0\.722, and precision=0\.555=0\.555, compared withF1=0\.486\\mathrm\{F1\}=0\.486for PCMCI\+ and0\.4820\.482for VAR\-Granger\. In several nonlinear and missing\-data categories, lagged correlation retains marginal associations that the conditional tests do not detect, while in the clean linear categories \(A1, A2\) PCMCI\+ reaches recall=1\.00=1\.00and precision=1\.00=1\.00and the correlation reference does not\. The aggregate advantage of lagged correlation therefore reflects category\-level variation rather than a uniform pattern across all data regimes\. On CausalRivers, the evaluated methods converge nearF1≈0\.60\\mathrm\{F1\}\\approx 0\.60, as high recall is offset by numerous links that are absent from the topology\-derived reference\.

These results should not be interpreted as evidence that the non\-causal reference models outperform the causal\-discovery methods\. Instead, they indicate whether substantial associative or predictive structure remains when a causal method reports few links\. Lagged correlation and the predictive reference are excluded from the consensus vote because neither provides conditional causal support\.

The patterns we observe are consistent with the assumptions and inferential reach of the methods under study\. When their conditional\-independence tests are well suited to the data, constraint\-based algorithms eliminate indirect links that are still detectable through marginal correlations\(runge2018chaos\)\. Partial correlation performs poorly when the underlying dependence is weakly linear, whereas CMIknn can capture more general functional dependencies, albeit with substantially higher demands on sample size and computation\(runge2018conditional\)\. Linear VAR approaches can still identify the linear component of a structural relationship that also includes nonlinear terms, which aligns with the partial Granger recovery seen in categories that mix linear and polynomial effects\(plagborgmoller2021local\)\. The similar performance scores attained by several methods on the small CausalRivers subgraphs call for a different type of caution, because the topology\-based reference is relatively dense given the total number of possible station pairs\(petersen2025negative\)and does not represent all hydrological mechanisms that generate dependence\. Consequently, across the three benchmarks, the difference between non\-causal and causal recovery depends on how well method assumptions match the data characteristics, rather than following a uniform performance ranking\.

### A\.3Consistency with theoretical identifiability conditions

The three extensions described in Section[3\.4](https://arxiv.org/html/2608.00198#S3.SS4)were developed against DGP\-Atlas, so Table[A\.2](https://arxiv.org/html/2608.00198#A1.T2)is an internal consistency check on those extensions rather than an independent validation\. It compares the modal causal\-audit recommendation on each family with the causal method that attains the highest mean F1 on that family\.

Table A\.2:Causal\-audit recommendation versus the highest\-mean\-F1 causal method on each DGP\-Atlas family\. Status distinguishes an exact method match from a partial match in which the recommended preprocessing is useful but another causal method has higher mean F1\. The random\-forest reference is excluded from the target because it is not a causal\-discovery method\.FamilyRecommendationHighest\-F1 causal methodF1ConfidenceStatusF1 \(clean\)VARLiNGAMVARLiNGAM0\.950\.85ExactF2 \(breaks\)PCMCI\+PCMCI\+0\.810\.85ExactF3 \(irregular\)PCMCI\+ \+ imputeVARLiNGAM0\.910\.83PartialF4 \(persistence\)VARLiNGAMVARLiNGAM0\.990\.84ExactF5 \(latent\)VARLiNGAMVAR\-Granger0\.770\.84NoF6 \(seasonality\)PCMCI\+VARLiNGAM0\.810\.73NoF7 \(nonlinear\)PCMCI\+VAR\-Granger0\.860\.85NoF8 \(non\-Gaussian\)VARLiNGAMVARLiNGAM0\.930\.85ExactF9 \(mixed\)PCMCI\+ \+ imputeVAR\-Granger0\.950\.84PartialF10 \(extreme\)VARLiNGAMVARLiNGAM0\.950\.84ExactA separate comparison asks whether the per\-method results follow the identifiability conditions described byrunge2023causalandassaad2022survey\. Table[A\.3](https://arxiv.org/html/2608.00198#A1.T3)maps selected benchmark cells to an assumption\-based expectation and the observed best method\.

Table A\.3:Assumption\-based expectation and observed mean\-F1 pattern on a representative subset of benchmark cells\. Cells where the random\-forest baseline leads \(RF\) correspond to regimes in which the evaluated conditional\-independence tests have low recall while the predictive reference maintains structure through association alone\.BenchmarkCellAssumption\-based expectationObserved patternTimeGraphA1, A2PCMCI\+ ParCorr \(clean linear, sufficient\)PCMCI\+ \(F1 = 1\.00\)TimeGraphB1, B2VARLiNGAM/TE \(linear test on polynomial fails\)RF \(F1 = 0\.80\)TimeGraphC1, C2VARLiNGAM/TE \(seasonality breaks ParCorr\)VARLiNGAM/TE \(F1 = 0\.80\)TimeGraphD2, D2CTE \(nonparametric on missing data\)RF \(sample size too low\)DGP\-AtlasF1 \(clean\)VARLiNGAM \(non\-Gaussian preferable to ParCorr\)VARLiNGAM \(F1 = 0\.95\)DGP\-AtlasF3 \(irregular\)VARLiNGAM/TE \(PCMCI\+ assumes regular grid\)VARLiNGAM \(F1 = 0\.91\); PCMCI\+ reaches 0\.78 after audit\-driven imputationDGP\-AtlasF5 \(latent\)LPCMCI \(admits latent confounders\)VAR\-Granger \(F1 = 0\.77, LPCMCI excluded\)DGP\-AtlasF7 \(nonlinear\)PCMCI\+ CMIknn \(theoretical\)VAR\-Granger \(F1 = 0\.86\); adaptive PCMCI\+ has lower recallDGP\-AtlasF8 \(non\-Gaussian\)VARLiNGAM \(ICA identifiability\)VARLiNGAM \(F1 = 0\.93\)#### A\.3\.1DGP\-Atlas diagnostic cases

Table[A\.2](https://arxiv.org/html/2608.00198#A1.T2)shows that the causal\-audit recommendation matches the method with highest F1, exactly on five families and partially on two others\. The remaining three families \(F5, F6, F7\) have aNostatus, meaning the recommended method is not the highest scoring one\. Each mismatch traces to a different source, ranging from evaluation\-target limitations through diagnostic ambiguity to finite\-sample behavior of nonlinear tests\.

The two partial matches arise from the irregular and mixed families \(F3, F9\), which exercise both interpolation and CI\-test routing\. Their released comparisons change more than one setting, so they show the behavior of the audit\-driven configuration path rather than an isolated interpolation effect\.

Of the three complete mismatches, the latent\-variable family \(F5\) highlights a shortcoming of the evaluation objective rather than of the recommendation mechanism\. Its scored graph includes only observed variables, and LPCMCI was not run due to the wall\-time constraint\. VARLiNGAM is recommended because the risks for nonlinearity, nonstationarity, and irregularity are all low, yet VAR\-Granger attains a higher F1 score on this family\. As the benchmark does not assess recovery of latent\-variable structure, this mismatch stems from the evaluation’s limited scope rather than from an error in the routing logic\.

The seasonality family \(F6\) illustrates a different kind of ambiguity\. The module suggests using PCMCI\+ with a deseasonalization preprocessing step, yet VARLiNGAM attains the best F1 score \(0\.813 versus 0\.731 for PCMCI\+\)\. This difference is relatively small, and the associated recommendation confidence \(0\.73\) is the lowest across all ten families, indicating that the diagnostic metrics do not strongly distinguish these two methods in this setting\.

The nonlinear family \(F7\) highlights a finite\-sample constraint\. Its structural equations combine linear and polynomial terms, and for the assessed values ofTT,NN, andτmax\\tau\_\{\\max\}, the RobustParCorr and CMIknn pipelines fail to recover all links\. VAR\-Granger performs best because it can identify the linear component of the structural equations\. This observation is tied to the particular configuration evaluated and should not be interpreted as establishing a general sample\-size threshold for CMIknn\.

#### A\.3\.2TimeGraph nonlinear and missing\-data cases

TimeGraph separates nonlinear\-test limitations from missing\-data limitations\. PCMCI\+ with ParCorr returns few or no links in polynomial categories where the linear component is weak\(runge2018chaos\)\. Missingness and lag truncation also reduce the effective sample, but the category\-level results in Figure[6](https://arxiv.org/html/2608.00198#S4.F6)show that incomplete series do not fail uniformly\.

Transfer entropy tests lags\{1,2,3,5\}\\\{1,2,3,5\\\}per pair and applies Benjamini–Hochberg correction across the fourpp\-values\. On B1, the ordered values do not satisfy the Benjamini–Hochberg step\-up criterion;α/4=0\.0125\\alpha/4=0\.0125is the first\-rank threshold, not a common threshold for every ordered value\. The correction therefore retains no tested lag for those weak signals\.

Category\-level differences are consistent with the inferential targets of the evaluated methods\. A linear VAR may recover the linear projection of a relationship that also contains polynomial terms\(plagborgmoller2021local\)\. Linear interpolation may also alter the residual distribution on which VARLiNGAM depends, although the present experiments do not isolate that mechanism\(shimizu2006lingam\)\. The random\-forest reference preserves predictive structure in some sparse categories, but it remains excluded from the causal consensus vote\.

The tested CMIknn configurations also return low or zero recall on the cubic\-quadratic B categories atT=1000T=1000andN=4N=4\. This result is specific to the evaluated signals, conditioning sets, and estimator settings\.

#### A\.3\.3Skeleton and orientation identifiability

Edge orientation requires a separate analysis from the skeleton because the endpoints have different identifiability conditions\. Under the benchmark’s sampling and graph conventions, a positive\-lag link is oriented from the variable att−τt\-\\tauto the variable attt\. A contemporaneous link may be identifiable only up to a Markov equivalence class when no collider or orientation rule fixes its direction\(spirtes2000causation;pearl2009causality\)\. PCMCI\+ returns an undirected mark \(o\-o\) when several orientations remain compatible with the tested independencies\(runge2020discovering\)\.

AutoCause reports each link at the orientation resolution its method identifies \(Section[2\.2](https://arxiv.org/html/2608.00198#S2.SS2)\)\. For the constraint\-based methods it preserves the Tigramite output convention, in which the array entry indexed by sourceiiatt−τt\-\\tauand targetjjatttorients lagged links by time order and retains the contemporaneous marker returned by the algorithm: an arrowhead where an unshielded collider or an orientation rule determines the direction, and an undirected mark \(o\-o\) where it does not\. Because the skeleton metric compares unordered pairs, it is insensitive to orientation and serves as the primary comparator \(Section[4\.1](https://arxiv.org/html/2608.00198#S4.SS1)\); orientation is examined separately below, distinguishing links oriented in the generative direction, links left undirected within the equivalence class, and links oriented against it\.

The released TimeGraph orientation audit distinguishes correct directions, reversals, and contemporaneous links left undirected\. Those outputs are provided underexperiments/timegraph\_validation/directed\_eval/\. They are not combined with skeleton F1 because a fixed\-direction score and an equivalence\-class\-aware score answer different questions\.

CausalRivers requires a separate interpretation because its reference specifies upstream\-to\-downstream links while methods may report lagged, contemporaneous, directed, or undirected outputs\. Skeleton recovery remains the primary cross\-method endpoint\. The directed scoring implementation and outputs are provided inexperiments/causalrivers\_validation/results\_directed/andexperiments/causalrivers\_validation/directed\_eval\.py\.

Table[A\.4](https://arxiv.org/html/2608.00198#A1.T4)reports three F1 score levels over the 30 subgraphs: skeleton adjacency, exact directed source\-target match, and the identifiable completed partially directed acyclic graph \(CPDAG\), which credits a correctly oriented link or a contemporaneous link correctly left undirected within the Markov equivalence class\. All methods recover the skeleton at comparable F1\. Exact\-direction F1 is lower for every method\. Under the identifiable CPDAG score, PCMCI\+ has the highest point estimate because its constraint\-based orientation rules leave unorientable contemporaneous links undirected rather than imposing a direction\. The random\-forest reference is excluded from causal\-method rankings\.

Table A\.4:Orientation performance on CausalRivers, reported as mean F1 across 30 subgraphs from three topology classes\. Skeleton F1 evaluates unordered adjacency, directed F1 requires an exact source–target match, and identifiable F1 credits either a correctly oriented link or a contemporaneous link that is correctly left undirected within its Markov equivalence class\. The random\-forest baseline is included only as a non\-causal predictive reference and is not ranked with the causal\-discovery methods\.MethodSkeletonDirectedIdentifiablePCMCI\+ \(ParCorr\)0\.600\.400\.57VARLiNGAM0\.600\.350\.42Transfer entropy0\.600\.290\.34VAR\-Granger0\.570\.350\.36Random forest0\.600\.350\.35

### A\.4Dataset\-level surrogate diagnostics

IAAFT surrogates provide a null model for cross\-variable dependence while retaining specified univariate properties\(schreiber2000iaaft\)\. For each variable, the procedure preserves the marginal distribution and approximately preserves the power spectrum\. Generating the variables independently disrupts their observed cross\-variable dependence\. The resulting null does not prove the absence of every causal mechanism\.

For each method and dataset, the surrogate edge rate \(SER\) is the mean fraction of possible directed links returned across 25 surrogate datasets\. The separation indicator equals one when the observed edge count exceeds the empirical surrogate edge\-count threshold defined in the released falsification script\. The reported separation rate is the mean of that indicator across datasets\. This dataset\-level diagnostic does not test individual links and is not used to assign the consensus\-support tiers in Section[4\.4](https://arxiv.org/html/2608.00198#S4.SS4)\. The finite\-sample rule is recorded inexperiments/falsification\_validation/results/falsification\_summary\.mdand the accompanying script\.

The analysis uses a stratified subset of 22 datasets: one DGP from each DGP\-Atlas family, six TimeGraph categories, and two CausalRivers subgraphs from each topology class\. PCMCI\+ uses ParCorr throughout this diagnostic so that the reported surrogate behavior is not mixed with adaptive CI\-test selection\. Table[A\.5](https://arxiv.org/html/2608.00198#A1.T5)reports the resulting SER and separation rates\.

Table A\.5:Dataset\-level IAAFT surrogate diagnostics on a 22\-dataset stratified subset with 25 surrogates per dataset\. SER is the mean fraction of possible directed links under the surrogate null\. Separation is the fraction of datasets that pass the finite\-sample edge\-count rule defined in the released falsification script\. Source:falsification\_summary\.mdandfalsification\_table\.tex\.DGP\-AtlasTimeGraphCausalRiversMethodSERSeparationSERSeparationSERSeparationVAR\-Granger0\.17690%0\.20733%0\.239100%VARLiNGAM0\.18770%0\.17733%0\.239100%PCMCI\+ \(ParCorr\)0\.25480%0\.06467%0\.042100%PCMCI\+ has the lowest SER on the evaluated TimeGraph and CausalRivers subsets\. The DGP\-Atlas subset includes a nonlinear F7 case for which fixed ParCorr is misspecified, but the analysis does not isolate that case’s contribution to aggregate SER\. A high separation rate means that a method reports more aggregate structure in the observed data than under the IAAFT null\. It does not validate individual links or classify links outside the CausalRivers topology reference as physically false\.

## Software and data availability

### Software availability

Name:AutoCause

Developer:Marco Ruiz, ISR\-Lisbon, Instituto Superior Técnico

Contact:marco\.rueda@tecnico\.ulisboa\.pt

Year first available:2026

Hardware:No specialized accelerator is required for the small benchmark cases\. Runtime and memory increase with the number of variables, record length, lag window, and selected method\. Tested on Apple M3 Pro, AMD Ryzen, and ARM\-based HPC nodes\.

Software:Python≥\\geq3\.10; tested on Python 3\.13\. A pinneduv\.lockfile records the software environment for dependency\-level reproducibility\.

Program language:Python

Program size:Approximately 27,000 Python source lines in the release, measured withcloc framework;pytest \-\-collect\-onlyreports 70 unit and integration tests for discovery methods, sample\-size diagnostics, consensus rules, and graph\-recovery metrics

Repository:License:AGPLv3\+

Cost:Free

Release:Version 0\.2\.0, tagged asPaper version\. The reviewer package contains the release metadata and software state used for the reported experiments\.

Error handling:Known insufficient\-data and assumption\-risk conditions are recorded as warnings; method failures are retained in the experiment log and are not converted into detected edges

Documentation:Usage examples, API reference, and HPC submission scripts are included in the repository

Reproducibility:The versioned paper release contains experiment scripts, figure\-generation scripts, configuration files, and output tables used for the reported analyses\. Every numerical value in the paper traces to a specific CSV or script in this release, and a file manifest maps each table and figure to its source\.

### Data availability statement

Benchmarks\.DGP\-Atlas is archived at[https://zenodo\.org/records/19409395](https://zenodo.org/records/19409395); TimeGraph is available at[https://github\.com/hferdous/TimeGraph](https://github.com/hferdous/TimeGraph); and CausalRivers is available at[https://github\.com/CausalRivers/causalrivers](https://github.com/CausalRivers/causalrivers)\.Software\.AutoCause is released under AGPLv3\+, and causal\-audit is available at[https://github\.com/marcoruizrueda/causal\-audit](https://github.com/marcoruizrueda/causal-audit)\. If the AutoCause repository is not public during review, the editor and reviewers receive a protected versioned package\.Reproducibility materials\.That package contains the software environment, experiment configurations, output CSV files, generated tables, figures, and reproduction instructions used for this article\.

## Funding

This work was supported by LARSyS through FCT funding \[grant numbers[10\.54499/LA/P/0083/2020](https://arxiv.org/html/2608.00198v1/10.54499/LA/P/0083/2020),[10\.54499/UIDP/50009/2020](https://arxiv.org/html/2608.00198v1/10.54499/UIDP/50009/2020), and[10\.54499/UIDB/50009/2020](https://arxiv.org/html/2608.00198v1/10.54499/UIDB/50009/2020)\]\. The funder had no involvement in study design, data collection, analysis, interpretation, manuscript preparation, or the decision to submit\.

## Declaration of competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper\.

\\printcredits

## Declaration of generative AI and AI\-assisted technologies in the manuscript preparation process

During the preparation of this work, the authors used ChatGPT \(OpenAI\) and Grammarly to reword and rephrase text, originally written by the authors\. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication\.

## Acknowledgments

The authors acknowledge institutional support from ISR\-Lisbon, Instituto Superior Técnico, and LARSyS\.

## References