causa

Causal discovery, identification, and estimation for Go. Pure standard library. Zero dependencies.

Status: early development — v0.16.0 released. Granger causality shipped in v0.1.0, PC-stable constraint-based discovery in v0.2.0, DirectLiNGAM directional discovery in v0.3.0, linear-SEM interventions + counterfactuals (the do-operator) in v0.4.0, FCI latent-confounder discovery (returning a PAG) in v0.5.0, the Shpitser–Pearl ID algorithm for causal-effect identification in v0.6.0, the IDC algorithm for conditional-effect identification in v0.7.0, the IDP algorithm for identification over a PAG (an equivalence class) in v0.8.0, the CIDP algorithm for conditional identification over a PAG in v0.9.0, numeric evaluation of the PAG estimand in v0.10.0, continuous (linear-Gaussian) estimand evaluation in v0.11.0, bootstrap uncertainty quantification of the evaluated effect in v0.12.0, and opt-in selection-bias discovery (FCI rules R5–R7) in v0.13.0, bounded discovery, bootstrap, and dense-state execution in v0.14.0, multivariate VAR/conditional Granger plus dependence-preserving bootstrap in v0.15.0, and temporal validity diagnostics plus multiplicity-safe Granger scans in v0.16.0. Pre-1.0, minor versions may break the API. Nothing below is claimed as shipped until it is implemented, tested against ground-truth datasets, and benchmarked. This README is kept honest by policy: capabilities are labeled exactly as they are.

What

causa is a general-purpose causal discovery, identification, and estimation library written in pure Go (stdlib only, CGO-free). Its algorithms cover different data regimes: pairwise Granger, reduced-form VAR, and conditional Granger are time-series methods; moving-block and stationary bootstrap preserve serial dependence. PC, FCI, DirectLiNGAM, Gaussian fitting, and the i.i.d. row bootstrap expect independent observational rows unless explicitly documented otherwise. The library is being built to power deterministic root-cause analysis in AIOpsFlow, but it has no dependency on or domain types from any host application.

Why

Correlation is not causation, and nowhere does that bite harder than in production telemetry: a network saturation event raises database CPU and API latency at the same time, and naive correlation blames the database. Distinguishing real causal structure from coincident movement requires actual causal mathematics — conditional independence, structural asymmetries, intervention modeling — not just pattern matching.

Today that mathematics lives almost entirely in Python (DoWhy, causal-learn, gCastle, Tigramite). A Go service that needs causal reasoning must ship a Python sidecar with a scientific stack — hundreds of megabytes of runtime and a supply chain to audit. causa exists to remove that hop: causal inference as a plain Go import, small enough to embed anywhere Go runs.

Data regimes and assumptions

The algorithms are composable, but their sampling assumptions are not interchangeable:

API family	Expected input	Load-bearing assumptions
GrangerTest	Two aligned chronological series	Appropriate lag order and approximately stationary innovation process
VAR fitting, stability/whiteness diagnostics, and VARGrangerTest / VARGrangerScan	Aligned multivariate chronological series	Linear reduced-form VAR, covariance stationarity, adequate lags, innovation residuals, multiplicity control for scans
MovingBlockBootstrap, StationaryBootstrap	Ordered observations from one approximately stationary process	Weak dependence and a scientifically defensible block length
PCStable, FCI	Variables observed over independent rows	Valid conditional-independence test, faithfulness; PC also requires causal sufficiency
DirectLiNGAM	Independent observational rows	Linear acyclic model, causal sufficiency, mutually independent non-Gaussian noise
Identify, IdentifyConditional, IDP/CIDP	A supplied causal graph or PAG	The graph and identification scope are valid; these APIs do not infer temporal structure
SampleGaussian, SampleDistribution, Bootstrap	Independent rows	Distribution fitting and the ordinary percentile bootstrap treat rows as i.i.d.

Raw chronological telemetry must not be passed to PC, FCI, DirectLiNGAM, or the i.i.d. bootstrap as if adjacent timestamps were independent. Detrending alone does not establish independence. For temporal uncertainty use a block bootstrap; for linear predictive dynamics use VAR/Granger. Neither choice repairs nonstationarity or turns predictive precedence into an intervention claim.

Bounded execution and memory safety

The zero-value APIs remain usable. Callers handling untrusted dimensions or latency-sensitive workloads can additionally bound the expensive paths:

• PCStableContext and FCIContext accept cancellation, expose DiscoveryDiagnostics, and honor MaxTests. A canceled or exhausted run returns no partial graph.
• NewDistributionWithOptions, SampleDistributionWithOptions, and Expr.EvaluateWithOptions reject integer-overflowing or oversized dense state spaces before allocation. The default maximum is DefaultMaxDenseCells; reviewed callers may override it.
• BootstrapContext, the moving/stationary block variants, and their Gaussian-effect helpers provide cooperative cancellation. BootstrapOptions.MaxResamples bounds callback work and replicate storage.
• FitVARWithOptions, SelectVARLagsWithOptions, VARGrangerTestWithOptions, and VARGrangerScan reject overflowing or oversized regression designs. The default cap is DefaultMaxVARDesignCells.
• VARModel.Stability bounds both companion-matrix cells and shifted-QR iterations. The default cell cap admits a 64-dimensional companion matrix; larger reviewed workloads must opt in.
• VARGrangerScanContext supports cooperative cancellation and rejects hypothesis families above DefaultMaxVARGrangerTests before fitting. Cancellation or failure never returns a partial family.

These are safety bounds, not statistical convergence guarantees. Increasing a budget can make a search more complete, but cannot repair violated causal assumptions or insufficient data.

Roadmap

Capability	Method	Status
Granger causality	Pairwise OLS autoregressions (QR-fitted) + F-test	Released in v0.1.0 — ground-truth-validated and benchmarked; flags confounders by design (see below)
Multivariate temporal prediction	Reduced-form VAR, AIC/BIC/HQIC lag selection, companion stability, Portmanteau whiteness, multiplicity-adjusted conditional Granger scans, moving/stationary bootstrap	Released through v0.16.0 — independently cross-checked against base R; unit-boundary fits are explicitly indeterminate; controls observed histories, not hidden causes; predictive rather than structural
Constraint-based discovery	PC-stable algorithm (conditional-independence tests) → CPDAG	Released in v0.2.0 — ground-truth-validated and benchmarked; recovers a Markov equivalence class, not a unique DAG (see below)
Directional discovery	DirectLiNGAM (deterministic, non-Gaussian noise) → causal order + weighted DAG	Released in v0.3.0 — ground-truth-validated and benchmarked; identifies a fully directed model when the noise is non-Gaussian (see below)
Interventions / counterfactuals	Linear SEM + do-operator (forward substitution; Pearl abduction–action–prediction)	Released in v0.4.0 — exact for a fully specified linear SEM; the general do-calculus identification problem (latent confounders) remains research (see below)
Latent-confounder discovery	FCI (Possible-D-SEP + Zhang's rules) → PAG	Released in v0.5.0 — ground-truth-validated and benchmarked; drops causal sufficiency, reporting latent common causes as bidirected (↔) edges; no selection bias by default, opt-in since v0.13.0 (see below)
Causal-effect identification	Shpitser–Pearl ID → symbolic estimand + discrete evaluator	Released in v0.6.0 — validated against brute-force truth on random latent SCMs; decides identifiability of P(y | do(x)) in a diagram with latent confounders, or proves non-identifiability with a hedge (see below)
Conditional-effect identification	Shpitser–Pearl IDC (do-calculus Rule 2 + m-separation + ID)	Released in v0.7.0 — validated against brute-force truth on random latent SCMs; identifies P(y | do(x), z), the effect of x on y within a context z (see below)
Identification over an equivalence class	Jaber–Zhang–Bareinboim IDP (pc-components + regions, Prop. 6/7 over induced PAGs)	Released in v0.8.0 — decides identifiability of P(y | do(x)) from a PAG (what FCI returns), not a single asserted diagram; decision cross-checked case-for-case against the reference PAGId implementation; symbolic (render-only) estimand; refuses a selection-bias PAG (see below)
Conditional identification over an equivalence class	Jaber–Zhang–Bareinboim CIDP (PAG do-calculus Rule 2 + definite-status m-separation + IDP)	Released in v0.9.0 — decides identifiability of P(y | do(x), z) from a PAG, the contextual effect; decision cross-checked against PAGId::CIDP; refuses a selection-bias PAG (see below)
Numeric evaluation of a PAG effect	Discrete evaluator over the identified IDP/CIDP estimand	Released in v0.10.0 — turns an identified PAG effect into the interventional table P(y | do(x)[, z]) from a discrete joint; validated against brute-force truth on random latent SCMs (PAG per SCM via an oracle FCI)
Continuous (linear-Gaussian) evaluation	Canonical-form Gaussian factor algebra over the identified estimand (Expr.EvaluateGaussian)	Released in v0.11.0 — evaluates an identified estimand on a normal observational joint, returning P(y | do(x)) as a Gaussian; exact for a linear-Gaussian model; validated against the closed-form structural effect (SEM.TotalEffect) on random latent SCMs
Uncertainty quantification	Nonparametric bootstrap over the evaluated effect (Bootstrap, Expr.BootstrapGaussianEffect) + distribution fitting (SampleGaussian, SampleDistribution)	Released in v0.12.0 — resamples the data to turn a point-estimated causal effect into a confidence interval; validated by its coverage on known linear-Gaussian SCMs (a nominal 95% interval covers the truth ≈95% of the time)
Selection-bias discovery	Zhang's rules R5–R7 in FCI (opt-in FCIOptions.SelectionBias) → PAG with undirected (—) edges	Released in v0.13.0 — admits selection bias into discovery, sound and complete for the class that also permits selection (Zhang 2008); off by default (PAG byte-identical to earlier versions)
Identification under selection bias	Recovering P(y | do(x)) from a selection-biased PAG (IDP/CIDP)	Research — a March 2026 preprint proposes complete identification from ancestral graphs under selection bias; it has not yet been independently reproduced here. The identifiers still refuse a PAG with an — edge (ErrSelectionBiasUnsupported) rather than return an unvalidated estimand

Granger tells you that series A helps predict series B — necessary but not sufficient for causation (confounders fool it). PC and LiNGAM answer different structural questions under stronger i.i.d./SCM assumptions; they are not an automatic upgrade to run directly on raw time series. PC recovers a Markov equivalence class, while DirectLiNGAM can identify a full direction only when its linear, causally sufficient, non-Gaussian-noise model is credible.

Constraint-based discovery (PC-stable)

PCStable(data, names, opts) recovers causal structure from a panel of continuous variables (one slice per variable). It is the order-independent PC-stable algorithm of Colombo & Maathuis (2014): it thins a complete graph down to a skeleton with conditional-independence tests, orients unshielded colliders, and closes under Meek's rules R1–R4. The default conditional-independence test (FisherZTest) is the linear-Gaussian partial correlation — computed by QR-residualization that reuses the same Householder solver as the Granger path — transformed by Fisher's z; the CITest extension point lets you supply another test for non-Gaussian or discrete data.

The rows supplied to PC must be treated as independent observations; serially autocorrelated samples invalidate the default Fisher-z calibration. The test of conditional independence $x_i \perp x_j \mid S$: residualize both variables on $[1, S]$, correlate the residuals, and refer the variance-stabilized statistic to the standard normal,

$$r_{ij\cdot S} = \mathrm{corr}\big(x_i - \hat{x}_i^{(S)},\; x_j - \hat{x}_j^{(S)}\big), \qquad z = \tfrac{1}{2}\,\ln\frac{1+r_{ij\cdot S}}{1-r_{ij\cdot S}},$$ $$T = \sqrt{n - |S| - 3}\;\lvert z\rvert \;\overset{H_0}{\sim}\; \mathcal{N}(0,1), \qquad p = 2\big(1 - \Phi(T)\big) = \mathrm{erfc}\big(T/\sqrt{2}\big);$$

the edge is deleted (independence accepted) when $p > \alpha$ (default $\alpha = 0.05$), and the levels stop growing once $n - |S| - 3 < 1$ (the honest small-sample cap described below).

Output is a CPDAG, not a DAG. Constraint-based discovery identifies structure only up to Markov equivalence. A directed edge A → C is compelled (every DAG consistent with the data agrees on it); an undirected edge A — B is reversible (the equivalence class contains DAGs orienting it either way). A chain A → B → C and a fork A ← B → C are indistinguishable from observational data, so their edges come back undirected; a collider A → C ← B is oriented, because it alone implies A ⫫ B. Read an undirected edge as "direction not identified", never as "no causal link".

Assumptions (standard for constraint-based discovery, and stated because they are easy to violate in practice): causal sufficiency (no unobserved common cause of two measured variables), faithfulness (every conditional independence in the data is entailed by the graph), and a correct CI test (the default assumes linear-Gaussian data). On small samples there is an honest cap: conditioning sets stop growing once n − |S| − 3 < 1, so independencies that need larger conditioning sets than the sample supports cannot be tested and some edges a larger sample would remove may remain.

Latent-confounder discovery (FCI)

FCI(data, names, opts) runs the Fast Causal Inference algorithm of Spirtes, Glymour & Scheines and returns a Partial Ancestral Graph (PAG). It is the constraint-based sibling of PCStable with one crucial difference: it drops causal sufficiency. PC assumes every common cause is measured; FCI does not, and where the data imply an unobserved common cause it says so — an edge A ↔ B — instead of inventing a spurious direct link. This is the difference between "the database causes the latency" and "something you are not measuring drives both".

A PAG is read endpoint by endpoint, not as whole arrows. Each edge end carries one of three marks: an arrowhead (>, a definite non-ancestor — compelled across the whole equivalence class), a tail (-, a definite ancestor), or a circle (o, undetermined). So the edge types mean:

• A → B — a definite direction: B is not a cause of A in any graph of the class.
• A ↔ B — a latent common cause: neither causes the other; something unobserved drives both. This is the structure causal sufficiency cannot express, and the reason FCI exists.
• A o→ B — the arrowhead at B is compelled, but A's end is undetermined.
• A o—o B — wholly undetermined direction (as an undirected CPDAG edge is).

Reading a circle as "no relationship" is the classic error — it means direction not identified.

Algorithm. FCI starts from the PC-stable skeleton (shared code, same order-independence), then does what PC cannot: the Possible-D-SEP refinement. The PC skeleton can leave an edge that is removable only by conditioning on variables not adjacent to the pair; FCI hunts for exactly those separating sets. A vertex $v$ is in $\mathrm{Possible\text{-}D\text{-}SEP}(a)$ when some path $a \dots v$ has, for every consecutive triple $\langle x, y, z\rangle$ on it, either $y$ a collider on the path ($x \mathbin{}!!\to y \leftarrow!!\mathbin{} z$) or ${x, y, z}$ a triangle:

$$v \in \mathrm{Possible\text{-}D\text{-}SEP}(a) \iff \exists\, \text{path } a \dots v \ \text{s.t. } \forall\, \langle x,y,z\rangle:\ y \text{ is a collider on the path} \ \lor\ x, z \text{ adjacent.}$$

For every still-adjacent pair FCI tests conditional independence against subsets of this set, deleting the edge (and recording the separating set) when a test passes. The Possible-D-SEP sets are snapshotted before any deletion, so — exactly like the PC-stable skeleton — the result does not depend on the order pairs are visited. Unshielded colliders are then oriented (from the refined separating sets), and the graph is completed by Zhang's orientation rules R1–R4 and R8–R10, applied to closure: R1–R3 propagate arrowheads that would otherwise create an unlicensed collider or cycle, R4 is the discriminating-path rule, and R8–R10 add the compelled tails that distinguish a definite → from an o→.

Selection bias — off by default, opt-in since v0.13.0. By default FCI assumes the sampling process induces no selection (no conditioning by how the data were collected): it runs R1–R4 and R8–R10, the sound and complete rule set for the latent-confounder-only class, and never produces a tail–tail (—) edge. Setting FCIOptions{SelectionBias: true} admits selection bias by additionally running Zhang's rules R5, R6 and R7 — the only rules that introduce — edges. An — edge records an ancestral relation to a latent selection variable (a variable the sampling implicitly conditioned on), so under this mode a tail at A on A — B means A is an ancestor of B or of a selection variable. The skeleton, collider and R1–R4/R8–R10 phases are identical in both modes; on data with no selection bias the two return the same PAG (R5–R7 have nothing to fire on), so the flag is safe to leave off and free to turn on. The result is again sound and complete (Zhang 2008) for the larger class that permits selection.

Note the downstream identifiers do not yet handle selection: IdentifyPAG / IdentifyConditionalPAG refuse a PAG carrying an — edge with ErrSelectionBiasUnsupported rather than return a silently wrong estimand (identification under selection bias is separate research; a March 2026 preprint now proposes an algorithm, but it has not yet been independently reproduced here — see the roadmap).

Assumptions (stated because, as everywhere here, violating them silently returns a plausible but wrong graph): independent observational rows, faithfulness, a correct CI test (the default FisherZTest assumes linear-Gaussian data), and — unless SelectionBias is set — no selection bias. Causal sufficiency is not required — that is the entire point. The honest small-sample cap carries over from PC: conditioning sets stop growing once n − |S| − 3 < 1.

Reference: Spirtes, Glymour & Scheines, Causation, Prediction, and Search (2nd ed., 2000), ch. 6 (FCI, Possible-D-SEP); Zhang, "On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias", Artificial Intelligence 172 (2008) (rules R1–R10); Colombo & Maathuis, JMLR 15 (2014) (order-independent skeleton).

Directional discovery (DirectLiNGAM)

DirectLiNGAM(data, names, opts) estimates a Linear Non-Gaussian Acyclic Model (LiNGAM) from a panel of continuous variables and returns a full causal order plus the weighted coefficient matrix B of the structural model x = B·x + e. It is the deterministic DirectLiNGAM method of Shimizu et al. (JMLR 2011) — not ICA-LiNGAM: there is no random ICA initialization, so, like the rest of causa, the same input always yields the same output. It iteratively finds the most exogenous variable by the paper's entropy-based independence measure (a maximum-entropy negentropy approximation built from E[log cosh] and E[x·exp(−x²/2)] moments), regresses it out of the rest, and recurses; the connection strengths are then least-squares estimates on the original data. Reused throughout is the same Householder-QR OLS solver as the Granger and PC paths.

The model and the mathematics behind the direction choice: for standardized data the structural model is

$$x = B\,x + e,$$

with $B$ strictly lower-triangular in the causal order and $e$ mutually independent, non-Gaussian disturbances. The differential entropy of a standardized variable $u$ is estimated by the maximum-entropy approximation (Hyvärinen 1998),

$$\hat{H}(u) \;\approx\; \tfrac{1}{2}\big(1+\ln 2\pi\big) \;-\; k_1\big(\mathbb{E}[\ln\cosh u] - \gamma\big)^2 \;-\; k_2\big(\mathbb{E}[\,u\,e^{-u^2/2}\,]\big)^2,$$

with $k_1 = 79.047$, $k_2 = 7.4129$, $\gamma = 0.37457$. For a candidate cause $x_i$ against $x_j$, the directional statistic is the log-likelihood ratio of the two directions,

$$T \;=\; \big(H(x_j) + H(r_i^{(j)})\big) \;-\; \big(H(x_i) + H(r_j^{(i)})\big),$$

where $r_i^{(j)}$ is the standardized OLS residual of $x_i$ on $x_j$ (and vice versa): a positive $T$ favors $x_i \to x_j$. Each candidate is scored by $\sum_j \min(0, T)^2$ — penalizing only the pairs that testify against its exogeneity — and the minimizer (a true root scores ≈ 0) is peeled off; ties break to the lowest index, keeping the whole procedure deterministic.

Where PC leaves an edge undirected, LiNGAM directs it. A chain A → B → C and a fork A ← B → C are one Markov equivalence class — indistinguishable to a constraint-based method — but DirectLiNGAM separates them, because non-Gaussian noise breaks the symmetry that made them equivalent. The output is a single fully directed DAG, not an equivalence class.

Assumptions (stated bluntly because violating them silently returns a plausible but wrong model): linearity, acyclicity, causal sufficiency, and — the load-bearing one — non-Gaussian, mutually independent noise (at most one disturbance may be Gaussian). On Gaussian noise the model is fundamentally unidentifiable: a linear-Gaussian SEM and its reverse fit the data equally well, so the recovered order is arbitrary and meaningless. DirectLiNGAM cannot detect this and returns a fully oriented model regardless; a LiNGAM result is only trustworthy when the non-Gaussianity assumption genuinely holds. This failure mode is pinned by an honest-failure test. One conscious divergence from the reference implementations: coefficient pruning is a simple absolute-magnitude threshold (LiNGAMOptions.PruneThreshold), not adaptive-lasso, which would require an L1 optimizer this stdlib-only library does not carry.

Interventions and counterfactuals (SEM + do-operator)

NewSEM(names, B, intercept) and FitSEM(data, names, order) build a linear structural equation model

$$x = c + B\,x + e,$$

with $B$ acyclic — $B_{ij}$ is the direct causal effect of $x_j$ on $x_i$, non-zero only when $j$ precedes $i$ in a causal order — an intercept vector $c$, and mutually independent, mean-zero disturbances $e$. NewSEM takes a hand-specified matrix (and topologically sorts it, rejecting a cycle); FitSEM estimates $B$ and $c$ from data by regressing each variable on its predecessors in a given causal order via the same Householder-QR solver as the rest of the library. The order a DirectLiNGAM run recovers is exactly such an input, so discovery and intervention compose: FitSEM(data, r.Nodes(), r.CausalOrder()).

The do-operator. Intervene(do) returns the interventional expectation $\mathbb{E}[x \mid \mathrm{do}(x_S = \xi_S)]$. Fixing the intervened variables severs their incoming edges; every other variable takes its structural mean given its parents. In causal order this is a single forward substitution,

$$\mathbb{E}[x_i \mid \mathrm{do}(\cdot)] = \begin{cases} \xi_i & i \in S,\\[2pt] c_i + \sum_j B_{ij}\,\mathbb{E}[x_j \mid \mathrm{do}(\cdot)] & i \notin S, \end{cases}$$

equivalently the reduced form $x = (I - B)^{-1}(c + e)$ with the intervened components clamped. TotalEffect(from, to) returns the total causal effect — the sum over every directed path of the product of its edge coefficients, the $(I-B)^{-1}$ entry — computed as $\mathbb{E}[\text{to} \mid \mathrm{do}(\text{from}=1)] - \mathbb{E}[\text{to} \mid \mathrm{do}(\text{from}=0)]$.

Counterfactuals. Counterfactual(observed, do) answers "given that we observed $x$, what would $x$ have been had we done $\mathrm{do}$?" by Pearl's abduction–action–prediction:

$$\underbrace{e = (I - B)\,x_{\text{obs}} - c}_{\text{abduction}}, \qquad \underbrace{\text{fix } \mathrm{do}\text{-vars, cut their in-edges}}_{\text{action}}, \qquad \underbrace{x'_i = c_i + \sum_j B_{ij}\,x'_j + e_i}_{\text{prediction}},$$

re-propagating the structural equations while holding the abducted disturbances fixed. With an empty intervention the counterfactual reproduces the observation exactly — the consistency the method satisfies by construction.

Scope (stated bluntly, as everywhere here). This is the identified linear-SEM case: it assumes you already hold the full structural model — every direct effect and the acyclic structure. It does not solve the general do-calculus identification problem — recovering an interventional distribution from an observational one plus a partially known graph with latent confounders — which is exactly what Identify (below) does. What ships here is exact for a fully specified linear SEM, where interventions and counterfactuals reduce to forward substitution through the structural equations. Linearity, acyclicity and causal sufficiency carry over from the model you supply; a mis-specified $B$ returns a plausible but wrong answer, silently.

Causal-effect identification (ID algorithm)

Identify(g, y, x) answers the question the linear-SEM path cannot: given only a causal diagram — not a fully specified model — can the interventional distribution P(y | do(x)) be computed from observational data at all, and if so, what is the formula? It is the Shpitser–Pearl ID algorithm, sound and complete for a single causal diagram: if it returns an estimand the effect is identifiable, and if it returns a hedge the effect is provably not identifiable — no algorithm could do better.

The input is a Diagram, an acyclic directed mixed graph (ADMG) with two edge kinds: directed X → Y (direct causation) and bidirected X ↔ Y (an unobserved common cause — a latent confounder). Causal sufficiency is not assumed; that is the whole point. This is the graph a v0.5.0 FCI run is evidence about — though Identify takes a single asserted diagram, not a PAG (identifying over an equivalence class is the separate IDP algorithm, v0.8.0, below).

Identifiability, and why it can fail. An effect is identifiable when every causal model consistent with the diagram that agrees on the observational P(V) also agrees on P(y | do(x)). When two models can match on all observables yet disagree on the interventional answer, no estimand exists — the bow arc X → Y with X ↔ Y is the canonical failure. The algorithm detects these via hedges (a C-forest structure); IDResult.Identifiable reports which case holds and IDResult.Hedge names the obstructing vertices.

The estimand and the discrete evaluator. When identifiable, the result is a symbolic Expr over the observational joint, built from marginalization, products and conditionals — the Tian–Pearl C-component factorization. For the classic back-door graph it reads

$$P(y \mid do(x)) \;=\; \sum_{z} P(y \mid x, z)\,P(z),$$

and for the front-door graph X → M → Y with X ↔ Y (identifiable despite the latent confounder)

$$P(y \mid do(x)) \;=\; \sum_{m} P(m \mid x)\sum_{x'} P(y \mid m, x')\,P(x').$$

The estimand is exact but not algebraically simplified, so it may not read character-for-character like the textbook form — its value is what is guaranteed. Expr.Evaluate turns it into numbers: give it a discrete observational joint (Distribution) and it returns the interventional table P(y | do(x)). This is how the implementation is validated — against brute-force truncated factorization on random discrete latent SCMs, the parameterization-independent proof of correctness — and it is a useful capability in its own right.

Scope. Sound and complete identification for one causal diagram without selection bias; non-identifiability is returned as a result, not an error. Out of scope (and stated so plainly): selection bias, identification over a PAG (IDP), and estimating the identified estimand from continuous samples — the evaluator is for discrete joints. Reference: Shpitser & Pearl, "Identification of Joint Interventional Distributions in Recursive Semi-Markovian Causal Models" (AAAI 2006); Tian & Pearl, "A General Identification Condition for Causal Effects" (AAAI 2002).

Conditional effects. IdentifyConditional(g, y, x, z) identifies P(y | do(x), z) — the effect of x on y within the subpopulation where the context z takes a value — the query that actually drives contextual decisions ("what does this remediation do to the SLO for this kind of service?"). It is the Shpitser–Pearl IDC algorithm: do-calculus Rule 2 moves each conditioning variable that behaves like an intervention out of z and into x — decided by an m-separation test (d-separation generalized to mixed graphs, X ⫫ Y | Z accounting for bidirected edges) in the graph with the edges into x and out of the candidate removed — until none remains movable; the residual is then

$$P(y \mid do(x), z) \;=\; \frac{P_x(y, z)}{P_x(z)} \;=\; \frac{P_x(y, z)}{\sum_{y} P_x(y, z)},$$

with the joint P_x(y, z) handed to the ID algorithm above. With z empty it is exactly Identify; if the joint is not identifiable, neither is the conditional effect. Reference: Shpitser & Pearl, "Identification of Conditional Interventional Distributions" (UAI 2006).

Identification over a PAG (IDP)

IdentifyPAG(g, y, x) (or the (*PAG).Identify method) closes the loop with discovery: it decides whether P(y | do(x)) is identifiable directly from a PAG — the Markov equivalence class that FCI returns — rather than from a single diagram you had to assert. This is the guarantee an autonomous system actually needs: when the graph itself was learned from data, the honest question is not "is the effect identifiable in one graph I picked?" but "is it identifiable in every graph the data leaves possible?" It is the Jaber–Zhang–Bareinboim IDP algorithm, sound and complete for this problem under the no-selection-bias scope. A PAG carrying an undirected (—) edge is outside that scope — such edges encode selection bias (which FCI produces only in its opt-in SelectionBias mode) — so IdentifyPAG and IdentifyConditionalPAG refuse it with ErrSelectionBiasUnsupported rather than return a silently wrong estimand; identification under selection bias is deferred (see the roadmap).

The difficulty over a single diagram is the circle mark. An edge X o→ Y leaves open both X → Y (where do(x) matters) and X ↔ Y (where it does not); the two disagree on P(y | do(x)), so the effect is not identifiable — a single undetermined endpoint flips the answer. IDP works by the same C-component logic as ID, lifted to PAGs: it reduces the query over the possible ancestors of Y, splitting the working set by buckets (o–o components), pc-components and regions, and peeling off marginalizable pieces via Propositions 6 and 7 — with edge visibility (whether a directed edge can be hiding confounding) judged on the full PAG throughout.

// Z ↔ X → Y: X is confounded with Z, but the edge X → Y is VISIBLE, so the
// effect is identifiable across the whole equivalence class.
g, _ := causa.NewPAG([]string{"Z", "X", "Y"}, []causa.PAGEdge{
    {A: 0, B: 1, MarkA: causa.Arrow, MarkB: causa.Arrow}, // Z ↔ X
    {A: 1, B: 2, MarkA: causa.Tail, MarkB: causa.Arrow},  // X → Y
})
r, _ := g.Identify([]int{2}, []int{1})
fmt.Println(r.Identifiable) // true

Scope and validation. IDP returns the identifiability decision — the sound, tested guarantee, cross-checked case-for-case against the authors' reference PAGId implementation on byte-identical adjacency matrices — together with the symbolic Prop. 6/7 estimand. To turn it into numbers, call (*PAGIDResult).Evaluate (see below); the raw Estimand is the exact formula and String()s faithfully. Out of scope for now: identification under selection bias — a selection-biased PAG (one with an — edge, as FCI's SelectionBias mode can produce) is refused, not identified. Reference: Jaber, Ribeiro, Zhang & Bareinboim, "Causal Identification under Markov Equivalence: Calculus, Algorithm, and Completeness" (NeurIPS 2022).

Conditional effects over a PAG. IdentifyConditionalPAG(g, y, x, z) (or (*PAG).IdentifyConditional) is to IdentifyPAG what IDC is to ID: it decides whether the contextual effect P(y | do(x), z) is identifiable from the equivalence class — "what does this remediation do to the SLO for this class of service?", now honest about latent confounding. It is the Jaber–Zhang–Bareinboim CIDP algorithm: the PAG do-calculus Rule 2 (a definite-status m-separation test in the manipulated PAG P_{\overline{W},\underline{X}} — edges into W removed, then visible edges out of X removed) moves each context variable that behaves like an intervention from z into x, after which the residual P_x(y, z) / \sum_y P_x(y, z) is handed to IDP. Same scope as IDP: the identifiability decision is the validated guarantee (cross-checked against PAGId::CIDP), and the estimand is evaluated with (*PAGIDResult).Evaluate. With z empty it is exactly IdentifyPAG.

Numeric evaluation of a PAG effect

(*PAGIDResult).Evaluate(joint) turns an identified IDP or CIDP effect into actual numbers: given a discrete observational joint P(V), it returns the interventional table P(y | do(x)) — or P(y | do(x), z) — read with ProbAt exactly as an Identify estimand's Evaluate output. The raw Prop. 6/7 estimand is a c-factor expression that, evaluated directly, carries "spectator" variables (non-ancestors of the outcome the c-factor form leaves in, which come out constant); Evaluate reduces them away and returns a clean table over Y ∪ X (∪ Z). So call Evaluate, not Estimand.Evaluate.

r, _ := g.Identify([]int{y}, []int{x})   // identifiable
tab, _ := r.Evaluate(joint)              // joint is a discrete *Distribution over P(V)
p, _ := tab.ProbAt(map[int]int{x: 1, y: 1}) // P(Y=1 | do(X=1))

Validation. Correctness is checked against brute-force truth on random discrete latent SCMs, the same parameterization-independent standard used for Identify: for each SCM the PAG is derived by an oracle FCI (a d-separation oracle fed to this package's own FCI, so the PAG is exactly the one the SCM's equivalence class implies), the observational joint and the true P(y | do(x)[, z]) are computed by marginalizing latents and intervening, and the evaluated estimand is required to match — across random parameterizations, including latent-confounded, Prop-6-heavy estimands. As with Expr.Evaluate, the evaluator is for discrete joints.

Continuous (linear-Gaussian) evaluation

Expr.EvaluateGaussian(joint) is the continuous companion to Expr.Evaluate: where the discrete evaluator reads a probability table, this reads a normal observational joint N(μ, Σ) (NewGaussian) and returns P(y | do(x)) as a Gaussian, exact for a linear-Gaussian model. It is the estimand backend most telemetry needs — continuous metrics rather than binned states.

The identified estimand is the same symbolic Expr the ID/IDC algorithms produce; only the numeric backend changes. A Gaussian factor in canonical (information) form — φ(x) = exp(−½·xᵀ K x + hᵀx + g) — closes exactly the three operations an estimand is built from: products add (K, h, g), a conditional's ratio subtracts them, and marginalizing a variable block out is a Schur complement. So EvaluateGaussian walks the same AST as Evaluate, swapping dense tables for Gaussian factors. The result is a GaussianFactor over the query's free variables Y ∪ X; fix the intervention with Condition to read the outcome's mean and covariance:

r, _ := causa.Identify(g, []int{y}, []int{x})    // identifiable
f, _ := r.Estimand.EvaluateGaussian(joint)       // joint is a *GaussianDistribution N(μ, Σ)
hi, _ := f.Condition(map[int]float64{x: 1})       // P(Y | do(X=1)) as a Gaussian
lo, _ := f.Condition(map[int]float64{x: 0})       // P(Y | do(X=0))
m1, _ := hi.MeanAt(y)
m0, _ := lo.MeanAt(y)                              // m1 − m0 is the causal slope dE[Y|do(X)]/dX

The interventional distribution P(y | do(x)) is a distribution over y for each fixed x, so the estimand factor is proper in the outcome but flat in the intervention (its x-block is degenerate by design) — never marginalize x out; Condition on it. A degenerate input joint, or an intermediate factor whose eliminated block is not positive definite, is reported as ErrNotPositiveDefinite rather than silently returning garbage — the Cholesky factorization that drives the algebra doubles as the positive-definiteness test.

Validation. Correctness is checked against closed-form structural truth on random linear-Gaussian latent SCMs, the continuous analogue of the discrete SCM harness: for each diagram a random linear-Gaussian model is drawn (one latent root per bidirected edge), the exact observed covariance Σ = (I−B)⁻¹ Ω (I−B)⁻ᵀ is formed, and the interventional slope read off the evaluated estimand is required to match the structural total effect the independent SEM.TotalEffect do-operator path computes — across random parameterizations, over back-door and (latent-confounded) front-door estimands. Scope. Exact for a linear-Gaussian model given the distribution; estimating Σ from finite samples and quantifying the resulting uncertainty is v0.12.0 (bootstrap).

Uncertainty quantification (bootstrap)

Everything above evaluates an estimand on a known distribution. On real data you hold a finite sample, so the effect is an estimate with sampling error. v0.12.0 closes that gap with a nonparametric bootstrap: resample the rows with replacement, recompute the effect on each resample, and read a confidence interval off the spread — turning "the number" into "the number ± robustness".

Two primitives bridge raw data to the evaluators (a gap until now — nothing fit a distribution from samples): SampleGaussian(data) fits a normal by its sample mean and covariance (for EvaluateGaussian), and SampleDistribution(data, card) builds an empirical joint (for Evaluate). Expr.BootstrapGaussianEffect composes the whole continuous path — fit → evaluate → unit contrast — into an interval for the causal slope dE[y | do(x)]/dx:

r, _ := causa.Identify(g, []int{y}, []int{x})   // identifiable back-door / front-door query
ci, _ := r.Estimand.BootstrapGaussianEffect(data, x, y,
    causa.BootstrapOptions{Resamples: 1000, Level: 0.95, Seed: 1})
fmt.Printf("effect %.3f, 95%% CI [%.3f, %.3f]\n", ci.Point, ci.Lower, ci.Upper)
significant := ci.Lower > 0 || ci.Upper < 0   // does the interval exclude no-effect?

The engine underneath, Bootstrap(n, stat, opts), is general: it resamples n row indices and calls your stat closure, so it bootstraps any scalar — a discrete P(y | do(x)), a regression slope, whatever stat returns — not just the Gaussian effect. A degenerate resample (say a collinear draw with no positive-definite covariance) is skipped; the percentile interval and a bootstrap standard error come back in a BootstrapResult.

Validation. A confidence interval is judged by its coverage: on a known linear-Gaussian SCM with a known true effect, we draw many fresh samples, build a 95% interval from each, and count how often it contains the truth — it lands near 95% (≈0.93 here, honestly slightly conservative for the percentile method), across back-door and latent-confounded front-door effects. Scope. The ordinary Bootstrap is the basic percentile bootstrap — no bias-correction/acceleration (BCa) — under an i.i.d.-rows assumption; it quantifies sampling variability, not model misspecification.

For ordered observations, MovingBlockBootstrap samples overlapping fixed-length contiguous blocks, while StationaryBootstrap uses Politis–Romano geometrically distributed circular blocks. Both expose context-aware variants and Gaussian-effect helpers. They retain serial dependence inside resampled runs, but still require an approximately stationary, weakly dependent process; block length is an explicit modeling choice, never guessed silently. The implementations follow Künsch (1989) and Politis–Romano (1994).

Granger causality

GrangerTest(cause, effect, lags) is available since v0.1.0 (import path github.com/jousudo/causa). It fits a restricted autoregression of effect on its own lags and an unrestricted one that adds cause's lags, both via a Householder-QR least-squares solver, and reports the F-statistic and its p-value:

$$\text{restricted:}\quad y_t = c + \sum_{i=1}^{p} a_i\,y_{t-i} + \varepsilon_t \qquad\quad \text{unrestricted:}\quad y_t = c + \sum_{i=1}^{p} a_i\,y_{t-i} + \sum_{i=1}^{p} b_i\,x_{t-i} + \varepsilon_t$$ $$F = \frac{(\mathrm{RSS}_r - \mathrm{RSS}_u)/p}{\mathrm{RSS}_u/(n - 2p - 1)} \;\overset{H_0}{\sim}\; F(p,\; n-2p-1), \qquad H_0:\ b_1 = \dots = b_p = 0,$$

with the p-value evaluated through the regularized incomplete beta function (continued-fraction form). "x Granger-causes y" is the rejection of $H_0$ — x's past improves the prediction of y beyond y's own past. Known limitation — confounding: if a hidden common cause drives both series, pairwise Granger can report causality even when no direct edge exists. This is inherent to the method, not a defect; the behavior is pinned by TestGrangerFlagsConfounder and documented on the function.

Multivariate VAR and conditional Granger

FitVAR(data, names, lags) estimates the reduced-form system

$$y_t = c + A_1 y_{t-1} + \cdots + A_p y_{t-p} + u_t,$$

where data[variable][time] is aligned at a common interval. SelectVARLags compares orders on the same response timestamps using the Lütkepohl forms of AIC, BIC, and Hannan–Quinn. VARModel.Stability computes the roots of the VAR(1) companion representation and reports one of three states: stable, unstable, or indeterminate. The indeterminate band around unit modulus is deliberate: numerical noise near a unit root must not become a high-confidence decision.

VARModel.WhitenessTest implements the multivariate Portmanteau test. Its small-sample-adjusted statistic is

$$Q_h^* = T^2 \sum_{\ell=1}^{h}\frac{1}{T-\ell} \operatorname{tr}\!\left(C_\ell^\top C_0^{-1} C_\ell C_0^{-1}\right) \;\overset{H_0}{\sim}\;\chi^2_{k^2(h-p)},$$

where the null is zero residual autocovariance through lag $h$. A small p-value rejects residual whiteness; a large value does not prove independence. ResidualAutocorrelation remains available for locating individual residual relationships.

The conservative workflow validates the fitted dynamics before using any direction to raise an application's confidence:

selection, _ := causa.SelectVARLags(data, []string{"cpu", "latency", "traffic"}, 12)
model, _ := causa.FitVAR(data, []string{"cpu", "latency", "traffic"}, selection.SelectedBIC)
stability, _ := model.Stability(nil)
whiteness, _ := model.WhitenessTest(selection.SelectedBIC+10, true)
rejectWhiteness, _ := whiteness.Reject(0.05)
if stability.Status != causa.VARStable || rejectWhiteness {
    // Fail closed: do not promote directions from a misspecified temporal model.
}
scan, _ := causa.VARGrangerScan(data, model.Nodes(), selection.SelectedBIC, nil)
for _, finding := range scan.Findings {
    if finding.Significant {
        fmt.Printf("%s -> %s adjusted-p=%.4g\n",
            finding.Test.CauseName, finding.Test.EffectName, finding.AdjustedPValue)
    }
}

VARGrangerTest tests one prespecified direction. VARGrangerScan tests every ordered pair inside the same full VAR and adjusts the complete family. Its zero-value default is Holm, which controls family-wise error without an independence assumption. Benjamini–Hochberg is an explicit, less-conservative FDR option under its usual independence or positive-dependence conditions; PValueNone is available only when the caller deliberately owns multiplicity elsewhere.

With two variables the conditional test is numerically equivalent to GrangerTest; with measured common drivers it can remove pairwise omitted-variable false positives. It cannot control a variable that was not measured. This is a predictive reduced-form VAR, not an SVAR: contemporaneous directions are not identified, and a significant lag test is not an intervention effect. Unit roots require a real unit-root/cointegration workflow; regime changes, nonlinear dynamics, irregular sampling, hidden causes, and post-selection reuse of the same data require additional analysis. The stability diagnostic is not a unit-root test, and Portmanteau failure to reject is not proof that the model is true.

The numeric implementation is locked to independent base-R oracles in scripts/var_oracle.R and scripts/var_validity_oracle.R. Method references: Hosking, “The Multivariate Portmanteau Statistic,” JASA 75 (1980), doi:10.1080/01621459.1980.10477520; Holm, “A Simple Sequentially Rejective Multiple Test Procedure,” Scandinavian Journal of Statistics 6 (1979), JSTOR 4615733; Benjamini and Hochberg, “Controlling the False Discovery Rate,” JRSS B 57 (1995), doi:10.1111/j.2517-6161.1995.tb02031.x.

Design principles

• Standard library only. No third-party imports, ever. CGO-free. The whole supply chain is the Go toolchain.
• Deterministic and reproducible. Same input, same output — causal claims must be auditable.
• Documented like it matters. Every exported symbol carries godoc; every feature ships with testable Example functions and benchmarks.
• Validated, not asserted. Algorithms are tested against synthetic ground-truth suites and cross-checked against reference implementations before any release claims them.
• Semantic versioning from v0.1.0. Pre-1.0, minor versions may break the API; this is stated here instead of discovered by users.

Installation

go get github.com/jousudo/causa

Contributing

Issues, discussion and code contributions are welcome; note that pre-1.0, minor versions may still adjust the API. See CONTRIBUTING.md for the stdlib-only constraint, godoc/Example/benchmark expectations, and PR etiquette.

License

Apache License 2.0