PCMCI+ Causal Discovery

• PCMCI+ is a constraint-based causal discovery algorithm that reconstructs causal graphs from time series by testing conditional independencies through PC-style condition-selection and MCI tests.
• It efficiently detects both lagged and contemporaneous causal relationships while addressing challenges such as autocorrelation, high dimensionality, and nonlinear dependencies.
• Applications in climate science, epidemiology, and psychology demonstrate PCMCI+'s robust performance and practical benefits in uncovering complex causal dynamics.

Searching arXiv for PCMCI+ and related time-series causal discovery papers to ground the article in the original literature. PCMCI+ is a constraint-based causal discovery algorithm for multivariate time series that extends PCMCI to infer both lagged and contemporaneous causal relations under autocorrelation, high dimensionality, and potentially nonlinear dependence. It operates on a time series graph whose nodes are variables indexed by time and whose edges represent direct cause–effect relations across lags, using a two-stage procedure that combines a PC-style condition-selection phase with Momentary Conditional Independence (MCI) tests, followed by orientation rules for contemporaneous links. In the original development, PCMCI addressed large nonlinear time series datasets and left contemporaneous links undirected (Runge et al., 2017); PCMCI+ added explicit contemporaneous skeleton discovery and orientation, and was formulated as order-independent and consistent in the oracle case under causal sufficiency and faithfulness assumptions (Runge, 2020).

1. Conceptual setting and formal problem

PCMCI+ addresses the problem of reconstructing a causal graph from an observed multivariate discrete-time stochastic process $X_t = (X_t^1,\dots,X_t^N)$, where causal effects may occur with delays $0 \le \tau \le \tau_{\max}$. In this representation, lagged edges $X^i_{t-\tau} \to X^j_t$ for $\tau>0$ encode directed causal influence forward in time, while contemporaneous edges at $\tau=0$ represent instantaneous dependencies at the time resolution of the data (Runge et al., 2017, Runge, 2020).

The method assumes a stationary, time-invariant structural causal process of the form

$X^j_t = f_j(P(X^j_t), \eta^j_t),$

where $P(X^j_t)$ denotes the set of causal parents of $X^j_t$, the $f_j$ are measurable functions, and the noise terms are mutually and serially independent (Runge, 2020). A link exists when the driver and target are not conditionally independent given the appropriate remainder of the graph. In the lagged setting used by PCMCI, this is expressed by conditioning on the rest of the past; FullCI directly tests this definition but is high-dimensional and can reduce effect sizes through overconditioning (Runge et al., 2017).

The assumptions underlying causal interpretation are standard for constraint-based causal discovery in time series: the Causal Markov condition, Faithfulness or Adjacency Faithfulness, causal sufficiency, stationarity or time-invariant mechanisms, and temporal order for lagged effects (Runge et al., 2017, Runge, 2020). Several application papers restate these assumptions with domain-specific wording. In wearable-sensor epidemiology, for example, causal sufficiency was operationalized by including likely confounders such as temperature, relative humidity, and activity level (Arvind et al., 2023). In psychological time series, stationarity, Markov, faithfulness, causal sufficiency, and time-ordering were explicitly adopted, with lagged links oriented by temporal order and contemporaneous links handled separately (Vitanza et al., 7 Jul 2025).

A recurring technical motivation is that strong autocorrelation degrades conventional CI-based methods. The 2020 PCMCI+ paper states that standard CI methods such as PC suffer from low recall and partially inflated false positives under strong autocorrelation, whereas PCMCI+ improves CI reliability by optimizing conditioning sets and can even benefit from autocorrelation for contemporaneous orientation through informative lagged–contemporaneous triples (Runge, 2020).

2. Algorithmic architecture

PCMCI+ inherits PCMCI’s two core phases and adds a contemporaneous conditioning and orientation stage. In PCMCI, the pipeline consists of condition selection followed by MCI testing (Runge et al., 2017). PCMCI+ extends this with a separate contemporaneous skeleton phase and graphical orientation rules for lag-0 edges (Runge, 2020).

The first phase is a PC-style condition-selection procedure. For each target variable $Y_t$, the algorithm initializes a broad candidate set of lagged parents up to $0 \le \tau \le \tau_{\max}$0; in PCMCI+ this may also include contemporaneous candidates depending on the phase (Runge et al., 2017, Debeire et al., 2023). Candidates are then pruned iteratively through CI tests with gradually increasing conditioning-set size. Rather than enumerating all subsets, the method restricts attention to subsets formed from the strongest currently retained candidates, ordered by a dependency statistic or minimum test statistic accumulated during pruning (Runge et al., 2017, Runge, 2020). This produces a compact, high-recall superset of the true parents, reducing the dimensionality of subsequent tests.

The second phase is MCI testing. For a candidate lagged link $0 \le \tau \le \tau_{\max}$1, PCMCI tests

$0 \le \tau \le \tau_{\max}$2

that is, dependence conditional on the estimated parents of the target excluding the driver and the estimated parents of the driver (Runge et al., 2017). The extra conditioning on the driver-side parent set is central: it accounts for autocorrelation in the source, produces innovation-like residuals under the null, and improves false-positive calibration (Runge et al., 2017). In the broader PCMCI+/MCI formulation, the same logic appears as conditioning on selected parents of the effect and an approximation to the parents of the cause (Debeire et al., 2023).

PCMCI+ modifies this basic design in two ways. First, it performs a contemporaneous skeleton phase in which lag-0 adjacencies are tested using the same MCI logic, but with conditioning subsets drawn from contemporaneous neighbors and lagged selected sets (Runge, 2020). Second, it orients the surviving contemporaneous edges using collider and propagation rules adapted from constraint-based graphical discovery (Runge, 2020, Vitanza et al., 7 Jul 2025).

A concise way to characterize PCMCI+ is that it separates lagged pruning, contemporaneous pruning, and final orientation while preserving the low-dimensional conditioning philosophy of MCI. This suggests that the method’s main innovation is not a new CI statistic but an optimization of which variables enter each CI test.

3. Conditional independence testing and statistical mechanics

PCMCI+ is test-agnostic. Its behavior depends both on the graph-search logic and on the CI test used to instantiate each independence query. The original PCMCI paper implemented three main classes of CI tests: ParCorr for linear-Gaussian relations, GPDC for nonlinear additive relations, and a fully nonparametric kNN conditional mutual information estimator (Runge et al., 2017).

For ParCorr, residuals are obtained by regressing $0 \le \tau \le \tau_{\max}$3 and $0 \le \tau \le \tau_{\max}$4 on the conditioning set $0 \le \tau \le \tau_{\max}$5, and the residual correlation is tested. The partial correlation is

$0 \le \tau \le \tau_{\max}$6

with test statistic

$0 \le \tau \le \tau_{\max}$7

and a $0 \le \tau \le \tau_{\max}$8-distribution null with $0 \le \tau \le \tau_{\max}$9 degrees of freedom (Runge et al., 2017). The same form is reiterated in later PCMCI+ discussions (Runge, 2020, Vitanza et al., 7 Jul 2025).

For nonlinear additive settings, GPDC first residualizes via Gaussian process regression and then tests dependence of residuals using distance correlation. The original paper emphasizes that the null distribution of the copula-transformed residual distance correlation can be precomputed by sample size, avoiding permutations (Runge et al., 2017). For fully nonparametric inference, CMI is estimated by a kNN estimator,

$X^i_{t-\tau} \to X^j_t$0

with significance assessed using a local permutation scheme in the original PCMCI formulation (Runge et al., 2017). Application papers frequently use the nearest-neighbor CMI estimator, sometimes referred to as CMIknn, to capture nonlinear dependencies in small or noisy time series (Arvind et al., 2023, Vitanza et al., 7 Jul 2025).

The statistical rationale for MCI is central. The original PCMCI paper states that MCI typically has much lower dimensionality than FullCI and “provably” larger or equal effect size, summarized by

$X^i_{t-\tau} \to X^j_t$1

which improves power (Runge et al., 2017). The PCMCI+ paper extends this argument to contemporaneous links, showing that the PCMCI+ MCI construction can have strictly larger effect size than a variant that does not condition on both sides’ lagged parent sets (Runge, 2020).

Autocorrelation handling is therefore not incidental. Conditioning only on the target’s past, or using unconditional dependence measures such as Corr, dCor, MI, or TE-like tests, tends to inflate false positives in autocorrelated data (Runge et al., 2017). PCMCI+ was explicitly designed to counter this by conditioning on parents of both variables and by structuring the graph search to keep those conditioning sets compact (Runge, 2020).

4. Contemporaneous discovery and orientation

The distinction between PCMCI and PCMCI+ is sharpest at $X^i_{t-\tau} \to X^j_t$2. The 2017 PCMCI paper leaves contemporaneous dependencies undirected and states that orientation rules for $X^i_{t-\tau} \to X^j_t$3 are not specified there (Runge et al., 2017). PCMCI+ was introduced specifically to discover and orient contemporaneous links in addition to lagged ones (Runge, 2020).

In PCMCI+, contemporaneous edges are first treated as adjacencies in the skeleton. During the MCI-based contemporaneous skeleton phase, each pair $X^i_{t-\tau} \to X^j_t$4 is tested for conditional independence using conditioning sets that combine small subsets of contemporaneous neighbors with selected lagged adjacencies of the target and source (Runge, 2020). If independence is found at any admissible conditioning set, the adjacency is removed and the corresponding separating set is stored.

Orientation then proceeds by rules adapted to time series. The PCMCI+ paper describes a collider orientation phase for unshielded triples and a subsequent orientation phase with rules $X^i_{t-\tau} \to X^j_t$5–$X^i_{t-\tau} \to X^j_t$6 (Runge, 2020). The rules operate under time-order and stationarity restrictions, so that only a limited family of triples is relevant. The collider phase can use “None,” conservative, or majority rules depending on how separating subsets involving the middle node are aggregated (Runge, 2020). The later psychological study summarizes this extension as a contemporaneous conditioning phase followed by standard PC orientation rules—collider, propagation, and common-child rules—after a refined skeleton is obtained (Vitanza et al., 7 Jul 2025).

One explicit orientation condition widely cited for a collider is: $X^i_{t-\tau} \to X^j_t$7 if $X^i_{t-\tau} \to X^j_t$8 and $X^i_{t-\tau} \to X^j_t$9 are nonadjacent, both adjacent to $\tau>0$0, and $\tau>0$1 is not in the separating set for $\tau>0$2 (Vitanza et al., 7 Jul 2025). In the fuller PCMCI+ formalism, propagation and cycle-avoidance rules then orient further contemporaneous edges when implied by already oriented structures (Runge, 2020).

The method is designed to be order-independent. PCMCI+ adopts PC-stable updates, removing adjacencies only after completing each conditioning-cardinality level, and uses consistent conflict marking with conservative or majority collider handling (Runge, 2020). This matters because contemporaneous orientation is otherwise especially sensitive to variable ordering in constraint-based algorithms.

A common misconception is to treat all same-time edges as uninterpretable. PCMCI+ instead interprets contemporaneous dependencies as effects at delays shorter than the sampling rate, a convention explicitly used in the sensor-data epidemiology study (Arvind et al., 2023). At the same time, the literature stresses that some contemporaneous edges may remain unoriented if orienting triples are absent or the rules are inconclusive (Runge, 2020, Vitanza et al., 7 Jul 2025).

5. Theoretical properties, complexity, and limitations

PCMCI and PCMCI+ are formulated with asymptotic guarantees under standard assumptions. For PCMCI, the original paper states that under causal sufficiency, Markov, faithfulness, stationarity, and a consistent CI test, the method recovers the true time series graph in the infinite-sample limit (Runge et al., 2017). For PCMCI+, the 2020 paper states soundness of the skeleton under the assumptions, completeness of the returned CPDAG under Adjacency Faithfulness with conservative collider handling, and completeness under standard Faithfulness with majority or standard rules (Runge, 2020).

The 2020 paper also proves a lemma characterizing the output of the lagged skeleton phase: the retained lagged set contains the true lagged parents and, at most, the lagged parents of contemporaneous ancestors (Runge, 2020). This explains why the first phase is intended as a permissive but efficient screening stage rather than a final edge-deciding step.

Complexity is reduced relative to naïve PC-style subset enumeration. The original PCMCI paper gives the worst-case number of CI calls in the condition-selection phase as $\tau>0$3, with MCI calls $\tau>0$4, hence overall worst-case polynomial in $\tau>0$5 and $\tau>0$6 (Runge et al., 2017). PCMCI+ preserves this efficiency profile by restricting subset enumeration in the lagged phase and by confining exponential worst-case behavior to contemporaneous subset search among same-time adjacencies rather than among all lagged embeddings (Runge, 2020).

The main limitations are also consistently stated across the literature. Violations of causal sufficiency due to hidden confounding can bias edge inference and especially complicate contemporaneous orientation (Runge et al., 2017, Runge, 2020). Small samples and highly nonlinear or high-dimensional conditioning can reduce power, particularly for GPDC and CMI-based tests (Runge et al., 2017). Nonstationarity can affect CI calibration and identifiability, a limitation emphasized in applied studies on psychology and industrial data (Vitanza et al., 7 Jul 2025, Ma et al., 17 Nov 2025).

A plausible implication is that PCMCI+ is most reliable when the graph is sparse, the sampling process is sufficiently regular, and domain knowledge helps choose $\tau>0$7 and relevant covariates. The literature strongly favors this interpretation, even where formal sensitivity analyses are incomplete.

6. Empirical performance and applications

The empirical evidence for PCMCI and PCMCI+ comes from both synthetic benchmarks and domain applications.

In the original PCMCI paper, a climate teleconnection case study examined the relation between the Nino 3.4 index and British Columbia land-surface air temperature using monthly data from 1979–2017 with $\tau>0$8. A known physical direction, Nino $\tau>0$9 BCT with lag around two months, was used as a reference. Correlation was significant in both lag directions and confounded, whereas FullCI detected $\tau=0$0 with partial correlation approximately $\tau=0$1 and $\tau=0$2. PCMCI yielded an MCI partial correlation approximately $\tau=0$3 with $\tau=0$4 and maintained high detection power $\tau=0$5 when adding further variables, because PC1 excluded irrelevant covariates and preserved effect size (Runge et al., 2017).

The same paper reported synthetic benchmarks for random linear and nonlinear networks with $\tau=0$6–$\tau=0$7, $\tau=0$8 cross-links, $\tau=0$9, and varying autocorrelation strengths. In linear models with ParCorr and $X^j_t = f_j(P(X^j_t), \eta^j_t),$0, Corr had high false positives, FullCI power dropped from about $X^j_t = f_j(P(X^j_t), \eta^j_t),$1 at $X^j_t = f_j(P(X^j_t), \eta^j_t),$2 to about $X^j_t = f_j(P(X^j_t), \eta^j_t),$3 at $X^j_t = f_j(P(X^j_t), \eta^j_t),$4, PC showed conservative p-values and false-positive inflation for strong autocorrelation, whereas PCMCI maintained high and robust power, correct false-positive control, and retained $X^j_t = f_j(P(X^j_t), \eta^j_t),$5 true parents up to $X^j_t = f_j(P(X^j_t), \eta^j_t),$6 with far fewer conditions than FullCI (Runge et al., 2017). In nonlinear settings, PCMCI-GPDC and PCMCI-CMI outperformed alternatives in power or false-positive control depending on the test, while high-dimensional FullCI variants often collapsed (Runge et al., 2017).

PCMCI+’s dedicated evaluation focused on lagged and contemporaneous discovery in synthetic additive models with linear and nonlinear functions, strong autocorrelation up to $X^j_t = f_j(P(X^j_t), \eta^j_t),$7, $X^j_t = f_j(P(X^j_t), \eta^j_t),$8 contemporaneous links, and average cross in-degree $X^j_t = f_j(P(X^j_t), \eta^j_t),$9 (Runge, 2020). Compared with PC, GCresPC, and VAR-LiNGAM, PCMCI+ had robust lagged and contemporaneous true positive rates, controlled false positive rates well under strong autocorrelation, and achieved markedly higher recall for contemporaneous orientation, often more than twice that of comparators (Runge, 2020). The same study reported that recall for contemporaneous orientation increases with autocorrelation because lagged–contemporaneous triples become informative (Runge, 2020).

Application work illustrates the method’s cross-domain use:

Domain	Configuration	Reported findings
Wearable-sensor epidemiology (Arvind et al., 2023)	113 asthmatic adolescents, 48-hour series, $P(X^j_t)$0 hours, partial correlation and nearest-neighbor CMI	Personalized PM2.5 $P(X^j_t)$1 respiratory-rate relations with peaks within the first hour and around 5–7 hours
Psychology (Vitanza et al., 7 Jul 2025)	45 participants, 22 symptoms, $P(X^j_t)$2, CMIknn with $P(X^j_t)$3, $P(X^j_t)$4	Idiographic mixed graphs with lag-1 and contemporaneous edges; aggregated networks reveal disorder-specific structures
Industrial foundry process (Ma et al., 17 Nov 2025)	Cluster-specific causal analysis of furnace cycles using Tigramite and hybrid CI integration	Recurrent voltage $P(X^j_t)$5 cooling-water-temperature and energy–temperature–weight relations across operational modes

The epidemiology study is notable for using PCMCI+ specifically because it can detect both lagged and contemporaneous dependencies and explicitly addresses strong autocorrelation via MCI testing (Arvind et al., 2023). It reported that at 15-minute resolution the cohort-wide probability of causal links from PM2.5 to respiratory rate showed three clear linear-test peaks at roughly 1 hour, 5 hours, and 7 hours, while nonlinear dependencies showed a sharp peak within the first hour, a decline to around 5 hours, and then gradual increase toward 8 hours (Arvind et al., 2023).

The psychology study used PCMCI+ with CMIknn to infer nonlinear symptom networks for 45 individuals and found pronounced idiographic heterogeneity, while aggregation by diagnosis revealed disorder-specific mechanisms consistent with prior psychopathological literature (Vitanza et al., 7 Jul 2025). It also included a synthetic validation in which PCMCI+ with CMIknn recovered linear and nonlinear lagged causality whereas VAR, TE, and PCMCI+ with partial correlation failed for nonlinear mechanisms (Vitanza et al., 7 Jul 2025).

The foundry study applied PCMCI+ within clusters of melting cycles to approximate stationarity and combined multiple CI tests into a hybrid graph. Across clusters, it found a stable core involving total energy consumption, furnace temperature, and material weight, and a recurrent delayed voltage $P(X^j_t)$6 cooling-water-temperature relation with lag of at least three 10-second steps (Ma et al., 17 Nov 2025).

7. Extensions, uncertainty quantification, and practical use

Later work has focused on improving uncertainty quantification and stability. A prominent extension is Bagged-PCMCI+, which wraps PCMCI+ in a block-bootstrap procedure that preserves temporal dependencies and lag structure, then aggregates the resulting graphs by majority voting (Debeire et al., 2023). For bootstrap replicate $P(X^j_t)$7, PCMCI+ is run on a resampled series, and a link-wise confidence is computed as

$P(X^j_t)$8

These confidences are interpreted as selection frequencies under block-bootstrap resampling rather than posterior probabilities (Debeire et al., 2023).

The bagging study reports that Bagged-PCMCI+ improves precision and recall relative to base PCMCI+, with especially pronounced gains in challenging regimes involving short samples, large numbers of variables, and high autocorrelation (Debeire et al., 2023). Using $P(X^j_t)$9 replicates, the mean absolute error between bootstrap confidence and estimated true link frequencies was reported to be on the order of $X^j_t$0–$X^j_t$1 for absent and existing links overall, with somewhat larger errors for contemporaneous links (Debeire et al., 2023). The method is therefore presented as both an uncertainty estimator and a variance-reduction device.

Practical guidance across the literature is comparatively consistent. $X^j_t$2 should cover plausible maximum causal delays and can be chosen using domain knowledge and unconditional lag functions or cross-correlations (Runge et al., 2017, Arvind et al., 2023). The CI test should match the expected mechanism class: ParCorr for linear-Gaussian settings, GPDC for nonlinear additive-noise relations, and CMI or CMIknn for highly flexible nonlinear dependence when sufficient sample size and runtime are available (Runge et al., 2017, Runge, 2020). High-resolution data can exacerbate multiple testing and false positives, so coarser resampling may be useful, particularly for nonlinear tests (Arvind et al., 2023). Application studies also emphasize including likely confounders in the multivariate system and checking information retention across resolutions or resampling schemes (Arvind et al., 2023).

Software implementation is centered on the Tigramite Python package, which provides PCMCI and PCMCI+ together with multiple CI tests and plotting utilities (Runge et al., 2017, Arvind et al., 2023, Vitanza et al., 7 Jul 2025). The method’s continued use across climate science, epidemiology, psychology, and industrial monitoring suggests that its main practical value lies in combining explicit time-order constraints with low-dimensional CI testing in settings where naïve dependence measures are overwhelmed by autocorrelation and dimensionality.