Finding Causal Signal in Complex Systems

Part III : A practical workflow for markets and economies: what predictive Shapley importance actually measures, why SHAP is not causality, and how to build a defensible pipeline anyway

Brian James Curry | Vector1 Research

TL;DR: Parts I and II established the distinction and the decision. This is the workflow.

Markets and economies violate many of the assumptions that make causal discovery clean. They contain latent common factors, feedback, regime change, temporal dependence, and large groups of correlated variables. The standard predictive toolkit random forests, permutation importance, and SHAP-style attribution is still useful. But it does not answer a causal question.

The central result of this article comes from a synthetic market system with 12 observed variables, one unobserved risk factor, and known ground truth. Across six independent replications:

• The Spearman correlation between predictive Shapley importance and the absolute true total causal effect was +0.21, with a standard deviation of 0.08.
• The Spearman correlation between predictive Shapley importance and membership in the target’s Markov blanket was +0.70, with a standard deviation of 0.14.
• The top-ranked feature was a child of the target in all six runs. In the structural model, that feature had exactly zero causal effect on the target.

The precise conclusion is not that “SHAP estimates the Markov blanket.” It does not, at least not as a general theorem. The result is narrower and more defensible:

In this simulation, a global Shapley allocation of predictive R² behaved much more like a noisy ranking of the target’s predictive neighborhood than a ranking of causal effects.

That distinction matters. A target’s predictive neighborhood can contain direct causes, downstream consequences, and variables that share a child with the target. It can also omit distant ancestors with large total effects once their mediators are observed.

This is not a reason to discard feature importance. It is a reason to give it the right job: conservative dimension reduction, not causal ordering and never automatic control selection.

1. Why markets are the hard case

Markets do not cooperate with textbook assumptions.

• Latent common factors: Asset returns, credit spreads, liquidity, and sentiment load on risk factors that are imperfectly measured or not measured at all. Causal sufficiency — the assumption that every relevant common cause has been observed — is often implausible.
• Feedback: Sentiment can move prices, while prices can move sentiment. A static directed acyclic graph cannot represent both directions at the same instant.
• Non-stationarity: The data-generating process changes. A graph learned in one monetary, volatility, or liquidity regime may not survive the next.
• Temporal dependence: Random train-test splits can leak information across time, and conventional independent-and-identically-distributed standard errors can understate uncertainty.
• Dimensionality with collinearity: Analysts may begin with dozens or hundreds of related variables. Conditional independence tests and variable-attribution methods become unstable precisely where they are most tempting to use.
• Limited intervention: Monetary policy, credit conditions, and market structure cannot usually be randomized by the analyst. Natural experiments may exist, but they are scarce and often local to one policy, population, or regime.

The workflow below is designed around those failures. It does not make them disappear.

2. The synthetic market system

The main simulation contains 12 observed variables including the target plus one unobserved risk factor. Innovations are non-Gaussian, and the graph is known by construction.

Python

EDGES = {
    ('policy_rate', 'liquidity'): -0.80,
    ('latent_risk', 'liquidity'): -0.30,
    ('liquidity', 'sentiment'): 0.60,
    ('latent_risk', 'sentiment'): -0.70,
    ('liquidity', 'credit_spread'): -0.50,
    ('sentiment', 'credit_spread'): -0.60,
    ('liquidity', 'capex'): 0.55,
    ('sentiment', 'capex'): 0.45,
    ('sentiment', 'demand_index'): 0.50,
    ('capex', 'demand_index'): 0.65,
    ('demand_index', 'equity_return'): 0.70,
    ('credit_spread', 'equity_return'): -0.55,
    ('liquidity', 'equity_return'): 0.40,
    ('equity_return', 'analyst_revisions'): 0.85,
    ('equity_return', 'media_volume'): 0.60,
    ('sentiment', 'media_volume'): 0.50,
    ('latent_risk', 'commodity_index'): 0.75,
    ('policy_rate', 'fx_rate'): 0.70,
}

The target is equity_return. Its local graph structure is:

• Direct parents: credit_spread, demand_index, and liquidity.
• Direct children: analyst_revisions and media_volume.
• Another parent of a child: sentiment, which is also called a spouse because it shares the child media_volume with the target.
• Markov blanket: the six variables above.
• Outside the blanket: capex, commodity_index, fx_rate, noise_var, and policy_rate.

The role labels are not mutually exclusive. In this graph, sentiment is a spouse relation through media_volume, but it is also an ancestor of equity_return through other paths. A spouse does not necessarily have zero causal effect. Markov-blanket membership is a statement about conditional prediction, not a complete classification of causal influence.[1][2]

The true total causal effects, computed by summing products over all directed paths into equity_return, are:

• liquidity: +1.456
• policy_rate: −1.165, despite being outside the blanket
• sentiment: +0.885
• demand_index: +0.700
• credit_spread: −0.550
• capex: +0.455, despite being outside the blanket
• analyst_revisions: 0.000
• media_volume: 0.000
• commodity_index: 0.000
• fx_rate: 0.000
• noise_var: 0.000

Before fitting a model, the core problem is already visible: policy_rate has the second-largest total causal effect in the system and is not in the target’s Markov blanket. Its effect is fully mediated by liquidity.

3. Stage 1 — Screening: the right and wrong jobs for feature importance

First, a terminology correction. SHAP, SAGE, permutation importance, and the measure used in this simulation are related, but they are not interchangeable.

• SHAP assigns Shapley-based attributions to a trained model’s output, originally at the level of an individual prediction.[3]
• SAGE allocates global predictive performance or loss across features while accounting for interactions.[4]

The simulation below refits a model for each sampled feature subset and allocates out-of-sample R² using Shapley values. I call it predictive Shapley R². It is closer to a Shapley allocation of universal predictive power than to ordinary TreeSHAP, but it is not a literal implementation of either standard SHAP or SAGE.

All three answer predictive attribution questions. None becomes causal simply because it uses Shapley values.

Python

def predictive_shapley_r2(
    model_factory,
    X_train,
    y_train,
    X_test,
    y_test,
    features,
    n_permutations=40,
    seed=0,
):
    """Monte Carlo Shapley allocation of out-of-sample R².

    A fresh model is fitted for every unique feature subset.
    model_factory should return a deterministically seeded estimator.
    """
    rng = np.random.default_rng(seed)
    contributions = {feature: [] for feature in features}
    cache = {}

    def value(subset):
        key = tuple(sorted(subset))
        if key not in cache:
            if not key:
                numerator = ((y_test - y_train.mean()) ** 2).sum()
                denominator = ((y_test - y_test.mean()) ** 2).sum()
                cache[key] = 1 - numerator / denominator
            else:
                model = model_factory()
                model.fit(X_train[list(key)], y_train)
                cache[key] = model.score(X_test[list(key)], y_test)
        return cache[key]

    for _ in range(n_permutations):
        subset = []
        previous = value(subset)

        for feature in rng.permutation(features):
            subset.append(feature)
            current = value(subset)
            contributions[feature].append(current - previous)
            previous = current

    return {
        feature: float(np.mean(values))
        for feature, values in contributions.items()
    }

The random forest predicts well: out-of-sample R² = 0.926.

In one representative run, the ranking was:

• analyst_revisions: predictive Shapley R² 0.2705 | permutation importance 0.9590 | abs. true effect 0.000 | direct child of the target
• capex: predictive Shapley R² 0.1842 | permutation importance 0.0014 | abs. true effect 0.455 | outside the blanket
• demand_index: predictive Shapley R² 0.1633 | permutation importance 0.0358 | abs. true effect 0.700 | direct parent
• liquidity: predictive Shapley R² 0.1081 | permutation importance 0.0082 | abs. true effect 1.456 | direct parent
• media_volume: predictive Shapley R² 0.1060 | permutation importance 0.1139 | abs. true effect 0.000 | direct child
• sentiment: predictive Shapley R² 0.0788 | permutation importance 0.0002 | abs. true effect 0.885 | parent of a target child and causal ancestor
• credit_spread: predictive Shapley R² 0.0297 | permutation importance 0.0185 | abs. true effect 0.550 | direct parent
• policy_rate: predictive Shapley R² 0.0045 | permutation importance 0.0004 | abs. true effect 1.165 | outside the blanket
• fx_rate: predictive Shapley R² 0.0008 | permutation importance −0.0002 | abs. true effect 0.000 | outside the blanket
• noise_var: predictive Shapley R² −0.0097 | permutation importance 0.0003 | abs. true effect 0.000 | outside the blanket
• commodity_index: predictive Shapley R² −0.0101 | permutation importance −0.0001 | abs. true effect 0.000 | outside the blanket

Read the top and bottom together.

The top-ranked feature, analyst_revisions, has a true causal effect of exactly zero. It is downstream of the target: revisions follow returns in this model; they do not produce them. policy_rate, with the second-largest causal effect in the system, sits near the bottom because its effect is mediated through liquidity. liquidity, the largest total cause, does not reach the top three.

Across six replications:

• Spearman correlation between predictive Shapley R² and absolute true causal effect: +0.207, with standard deviation 0.083.
• Spearman correlation between predictive Shapley R² and the binary Markov-blanket indicator: +0.702, with standard deviation 0.135.
• Top-ranked feature was a child of the target: 6 of 6 runs.
• Markov-blanket precision among the top six features: 89%.
• Mean rank of liquidity (total effect 1.456): 5.0 of 11.
• Mean rank of policy_rate (total effect 1.165): 7.7 of 11.

Six replications are enough to demonstrate a failure mode, not to establish a universal benchmark. The result should be read as a controlled counterexample to causal interpretation, not as a claim that every SHAP implementation will produce this exact ordering.

3.1 What the Markov blanket does — and does not — mean

In a faithful directed acyclic graph, the target’s Markov blanket contains its parents, its children, and the other parents of its children. Conditional on that set, the target is independent of the remaining variables. That makes a Markov boundary a natural object for prediction and feature selection.[1][5]

It does not make every blanket member a cause. The simulation shows two distinct errors that arise when predictive importance is read causally.

• False causal promotion: A downstream variable may be highly predictive because it records information produced by the target. Here, analyst_revisions is the strongest predictor and has no causal path back to returns.
• False causal omission: A distant ancestor may contribute little incremental prediction once its mediators are present. policy_rate acts on returns through liquidity, so its conditional predictive contribution is small even though its total effect is large.

The second failure is especially dangerous. A variable discarded during screening cannot be rescued by a later discovery or estimation stage.

3.2 Why screening still matters

The raw number of labeled directed acyclic graphs grows super-exponentially with the number of variables. If each labeled DAG receives equal prior mass, the structural search space has approximately:

• 6 variables: 3.78 million DAGs, or 21.9 bits.
• 8 variables: 783.7 billion DAGs, or 39.5 bits.
• 11 variables: 3.16 × 1⁰²² DAGs, or 74.7 bits.
• 20 variables: 2.34 × 1⁰⁷² DAGs, or 240.4 bits.
• 40 variables: 1.12 × 1⁰²⁷⁶ DAGs, or 917.0 bits.

Reducing 40 candidate variables to eight reduces this combinatorial description length from 917 bits to about 40 bits a 95.7% reduction.

That number is conceptual. It does not mean a discovery algorithm enumerates every DAG, guarantees a 23-fold runtime improvement, or proves that the screen retained the true graph. It shows why dimension reduction is economically valuable and statistically dangerous at the same time.

3.3 A conservative screening rule

Feature screening should optimize recall, not pretend to identify causes.

• Take a union, not an intersection: Combine several broad screens such as predictive Shapley R², mutual information, and bootstrapped sparse-model support — and retain the union. A false positive costs computation; a false negative can invalidate the entire analysis.
• Force-include plausible exogenous drivers and confounders: Policy changes, regulatory dates, commodity shocks, calendar structure, and known common causes should not be dropped simply because a flexible predictor assigns them little incremental importance.
• Use temporal and domain restrictions before statistical screening: Remove impossible edges and separate variable tiers using knowledge of measurement timing and system mechanics.
• Consider local screening around multiple variables, not only the final target: Target-only screening is designed to find a predictive neighborhood and can miss upstream structure. A union of local neighborhoods improves recall, although it may return enough variables to reduce the computational savings.
• Run sensitivity to the retained set: If the graph or estimand changes materially when borderline variables are restored, the screen is part of the uncertainty and must be reported as such.

4. Stage 2 — Discovery on the reduced system

Once the candidate set is manageable, method choice should follow the data-generating assumptions rather than habit.

PC or GES

• Appropriate as a baseline for acyclic, causally sufficient systems under Markov and faithfulness assumptions.
• Typically returns a Markov-equivalence class represented by a CPDAG, not a uniquely oriented graph.
• Fragile when conditional independence tests are weak, variables are highly correlated, or important common causes are unmeasured.[2]

FCI or RFCI

• Designed for settings that may contain latent confounders and selection variables.
• Returns a partial ancestral graph, or PAG, whose unresolved endpoints are part of the result rather than an inconvenience to erase.[6]
• A PAG may not support a unique adjustment set or point-identified causal effect. In that case, the honest output is ambiguity, a set of compatible effects, or failure to identify.

LiNGAM-family methods

• Can orient a linear acyclic system when structural errors are mutually independent and non-Gaussian.[7]
• Fat-tailed returns may satisfy one useful ingredient non-Gaussianity but they do not provide “orientation for free.” Volatility clustering, latent confounding, nonlinear dynamics, and dependent errors can violate the rest of the model.
• Use only after checking whether linearity, acyclicity, error independence, and the treatment of latent variables are plausible.

NOTEARS and its extensions

• Converts acyclicity into a smooth optimization constraint, making score-based structure learning computationally attractive.[8]
• The original NOTEARS model is linear. Nonlinear behavior requires a nonlinear extension; the optimization framework alone does not supply it.[9]

PCMCI-family methods

• Designed for high-dimensional time series with lagged dependencies and autocorrelation.[10]
• Plain PCMCI, PCMCI+, and latent-variable extensions do not share identical assumptions. Choose the version that matches whether contemporaneous effects and hidden common causes are allowed.
• Sampling frequency, maximum lag, stationarity, and conditional-independence testing remain substantive assumptions.

No method removes the need for structural judgment. Discovery output is a hypothesis set, not a finding. Bootstrap edges, vary tuning parameters and conditioning tests, report stability, and retain the equivalence class or PAG rather than forcing a single graph when the data do not identify one.

5. Stages 3–5 — Role assignment, identification, and estimation

Discovery and estimation solve different problems.

• Role assignment asks which variables are plausibly upstream, downstream, confounding, mediating, or colliding. When edge directions remain ambiguous, those roles remain ambiguous too.
• Identification asks whether the causal effect can be expressed from the observed distribution under the assumed graph through a backdoor set, front-door structure, instrument, or another identification strategy.[1]
• Estimation computes the identified quantity. Flexible estimators can reduce functional-form bias, but they cannot repair an unidentified estimand or a bad adjustment set.

In applied work, DoWhy is useful for making graph assumptions and identification steps explicit, while EconML provides estimators such as LinearDML and CausalForestDML for heterogeneous effects.[11][12] The mechanism below implements a simple cross-fitted, partially linear DML estimator directly.[13]

Python

def dml(df, treatment, outcome, controls, n_splits=5, seed=0):
    """Illustrative cross-fitted partially linear DML for i.i.d. data."""
    if controls:
        W = df[list(controls)].to_numpy()
    else:
        W = np.zeros((len(df), 1))

    T = df[treatment].to_numpy()
    Y = df[outcome].to_numpy()
    residual_t = np.zeros_like(T)
    residual_y = np.zeros_like(Y)

    def nuisance_model():
        return RandomForestRegressor(
            n_estimators=200,
            min_samples_leaf=10,
            random_state=seed,
            n_jobs=-1,
        )

    folds = KFold(
        n_splits=n_splits,
        shuffle=True,
        random_state=seed,
    )

    for train_idx, test_idx in folds.split(W):
        if controls:
            model_t = nuisance_model().fit(W[train_idx], T[train_idx])
            model_y = nuisance_model().fit(W[train_idx], Y[train_idx])

            residual_t[test_idx] = (
                T[test_idx] - model_t.predict(W[test_idx])
            )
            residual_y[test_idx] = (
                Y[test_idx] - model_y.predict(W[test_idx])
            )
        else:
            residual_t[test_idx] = T[test_idx] - T[train_idx].mean()
            residual_y[test_idx] = Y[test_idx] - Y[train_idx].mean()

    theta = (residual_t @ residual_y) / (residual_t @ residual_t)
    epsilon = residual_y - theta * residual_t
    se = np.sqrt(
        ((epsilon ** 2) @ (residual_t ** 2))
        / (residual_t @ residual_t) ** 2
    )

    return theta, se

This code is appropriate for the article’s independent synthetic draws. It is not a production template for market time series. Real temporal data require time-respecting or blocked cross-fitting and dependence-robust inference, such as HAC standard errors or a suitable block bootstrap.

5.1 The cost of choosing controls predictively

Now estimate the total effect of demand_index on equity_return. The true effect is +0.700. The DML estimator is held constant; only the adjustment set changes.

• No controls: estimate 1.173 | standard error 0.0092 | bias +0.473
• All other observed variables: estimate 0.284 | standard error 0.0098 | bias −0.416
• Top six predictive Shapley variables: estimate 0.261 | standard error 0.0102 | bias −0.439
• Graph-derived backdoor set {capex, sentiment}: estimate 0.681 | standard error 0.0168 | bias −0.019

Using the top six predictive features as the adjustment set produces a 63% relative error. The reason is structural. The selected controls include analyst_revisions and media_volume, both descendants of the outcome. That creates post-outcome adjustment. For media_volume, which is caused by both equity_return and sentiment, conditioning also opens a collider path.

The most dangerous part is the apparent precision:

• Predictively selected estimate: 0.261 ± 0.010, with an approximate 95% confidence interval of [0.241, 0.281].
• True effect: 0.700.

The wrong specification has a smaller standard error than the graph-derived specification, and its interval does not come close to the truth. Cross-fitting and flexible nuisance models worked as designed. The adjustment set was wrong.

Feature importance cannot select an identification set. The variables that best predict an outcome and the variables that identify a causal effect are different mathematical objects.

6. Stage 6 — When observed adjustment is not enough

The main simulation already contains an unobserved risk driver, but the chosen demand_index example has an observed backdoor set that blocks the relevant paths. That is not guaranteed in practice.

To isolate the harder failure, consider a second system in which an unobserved factor U confounds both treatment and outcome. The true treatment effect is +0.500.

• Adjusting only for observed covariates {W}: produces +0.078, a bias of −0.422.
• Adjusting for {W, U}, as if U were observed: produces +0.489, a bias of −0.011.

The observed-only estimate is attenuated by approximately 84%. The regression output alone cannot tell us that an unmeasured factor caused the attenuation.

Sensitivity analysis does not discover the missing variable, but it can quantify how strong omitted confounding would need to be to change the conclusion. Cinelli and Hazlett’s robustness value expresses that strength through partial R² measures.[14]

Python

f_squared = (theta / se) ** 2 / degrees_of_freedom
robustness_value = 0.5 * (
    np.sqrt(f_squared ** 2 + 4 * f_squared) - f_squared
)

In this example, the Robustness value RV(q = 1) is 0.092.

Under the equal-strength interpretation, omitted confounding associated with both the treatment and outcome residuals at partial R² values of about 9.2% could reduce the estimate to zero.

In many market systems, that is not an implausibly large hidden factor. The defensible conclusion is therefore not that the effect has been established. It is that the estimate is not robust to modest omitted confounding.

Practical responses include:

• Use FCI or an appropriate latent-variable time-series method when its assumptions fit, and retain the resulting PAG rather than converting uncertainty into a fully directed graph.
• Add theoretically justified observed proxies or factor controls, while recognizing that “de-factoring” can also remove causal variation if performed mechanically.
• Use negative controls, natural experiments, or instruments where their assumptions can be defended. A merely correlated instrument is not an instrument.
• Report method-appropriate sensitivity analysis with the point estimate. Robustness values are especially useful for linear omitted-variable-bias analysis; nonlinear, dynamic, or machine-learning estimands may require other sensitivity methods.
• Treat non-identification as a legitimate result.

7. Stage 7 — When acyclicity fails: represent feedback through time

Suppose sentiment moves returns and returns move sentiment at different lags:

Python

for t in range(1, n):
    sentiment[t] = (
        0.30 * sentiment[t - 1]
        + 0.55 * returns[t - 1]
        + noise()
    )

    returns[t] = (
        0.10 * returns[t - 1]
        - 0.35 * sentiment[t - 1]
        + noise()
    )

The simulated contemporaneous correlation is −0.074: small, symmetric, and directionally uninformative. The lagged model recovers both directions:

• return_(t−1) → sentiment_t: estimated +0.545; true coefficient +0.55.
• sentiment_(t−1) → return_t: estimated −0.357; true coefficient −0.35.

The contemporaneous relationship is near zero because the lagged effects have opposite signs and partially cancel. A static analysis could conclude that little is happening in a system with strong bidirectional coupling.

When the causal delays are correctly represented and there is no instantaneous cycle, indexing variables by time converts the feedback loop into an acyclic time-unrolled graph. That is the structure exploited by lag-aware discovery methods.

Temporal order is a powerful constraint: an observed event cannot cause an event already in its past. But it is not infallible. Coarse sampling can compress faster processes into the same period, asynchronous measurement can misorder events, and an omitted common driver can produce apparent lagged causation. A static graph forced to choose one direction in this system would be wrong. The correct answer is that both directions operate, at different lags, with different signs.

8. The full workflow

Stage 0 — Define the causal question and domain constraints

• State the treatment, outcome, population, horizon, and whether the target is a total, direct, or heterogeneous effect.
• Encode temporal order, variable tiers, forbidden edges, known exogenous drivers, and plausible common causes.
• Separate variables measured before treatment from mediators and post-outcome variables.

Stage 1 : Screen conservatively

• Use the union of multiple broad screens rather than a single importance ranking.
• Force-include domain-mandated drivers and confounders.
• Screen around multiple nodes if ancestral structure matters.
• Deliverable: A smaller candidate system with documented inclusion rules.
• Not delivered: Causal direction, an adjustment set, or a causal ranking.

Stage 2 : Discover structure under explicit assumptions

• Choose PC/GES, FCI/RFCI, LiNGAM, NOTEARS, or a PCMCI-family method based on acyclicity, latent variables, temporal structure, functional form, and error assumptions.
• Bootstrap edges and vary key specifications.
• Deliverable: A stable CPDAG, PAG, or set of candidate graphs — not necessarily one DAG.

Stage 3 : Assign roles without erasing ambiguity

• Distinguish plausible causes, effects, mediators, confounders, and colliders.
• Where directions remain unresolved, carry the uncertainty forward.
• Deliverable: Admissible causal roles across the surviving graph set.

Stage 4 : Identify the estimand

• Derive a valid backdoor set, front-door strategy, instrument, or other identifying expression.
• If compatible graphs imply different adjustment sets or effects, report partial identification or non-identification.
• Deliverable: An estimand justified by assumptions, not by feature importance.

Stage 5 : Estimate with a method matched to the data

• Use DML or causal forests when flexible nuisance estimation or heterogeneity is required.
• Use cross-fitting as required by the estimator.
• For time series and panels, use time-respecting folds and dependence-appropriate inference.
• Deliverable: An estimate and uncertainty conditional on the identification assumptions.

Stage 6 : Challenge the result

• Test alternative graph and adjustment specifications.
• Use placebo outcomes or treatments, negative controls, subset stability, regime splits, and sensitivity analysis where appropriate.
• Treat these as diagnostics, not proofs that the causal assumptions are true.
• Deliverable: A map of how the conclusion fails.

Stage 7 : Simulate decisions, not just coefficients

• Propagate surviving structural models under explicit interventions.
• Rank interventions by expected value and information value.
• Stress-test across regimes, graph uncertainty, and parameter uncertainty.
• Deliverable: A decision distribution, not a single deterministic forecast.

The load-bearing distinction is this: Screening asks what helps predict the target. Role assignment asks what is upstream of the target. Identification asks which variables must be conditioned on to estimate a particular effect. These are three different operations.

Conflating them is how analyst_revisions ends up in a trading thesis — or in an adjustment set.

9. Failure modes specific to markets

• Predictive R² treated as causal validation. The random forest achieved R² = 0.926 while ranking a zero-effect child first. Predictive accuracy validates prediction on the evaluation distribution, not causal interpretation.
• Ordinary SHAP, SAGE, and refit-based Shapley values treated as interchangeable. They allocate different value functions and can produce different rankings. Their common use of Shapley values does not make any of them causal.[3][4]
• Predictive importance used to select controls. In the simulation, this produced a 63% relative error with a narrow confidence interval.
• The ancestor problem ignored. policy_rate had the second-largest total effect and an average predictive-importance rank of 7.7 out of 11.
• Static graphs fitted to feedback systems. A contemporaneous correlation of −0.074 concealed two strong lagged effects with opposite signs.
• Random cross-validation used on time series. Leakage can inflate predictive performance and contaminate nuisance-model estimates. Time-respecting validation is part of the causal design, not a cosmetic modeling choice.
• Causal sufficiency assumed by silence. Failing to discuss latent confounding does not remove it.
• One graph treated as permanent. Structure should be tested across regimes and rolling windows. A changing graph may be the substantive finding.
• A point estimate reported when the graph identifies only a set. FCI can expose uncertainty about orientation; it does not authorize the analyst to resolve it by preference.

10. Limitations

This is a synthetic counterexample, not a market benchmark. The system is linear, the ground truth is known, and the six replications use 6,000 observations each. Real markets contain nonlinear thresholds, volatility clustering, strategic adaptation, measurement error, and regime change. The coefficients do not generalize.

The importance measure is predictive Shapley R², not ordinary SHAP. Local TreeSHAP, interventional SHAP, conditional SHAP, SAGE, permutation importance, and refit-based Shapley values allocate different objects. The broader causal warning applies to all predictive attributions, but the reported rank correlations belong only to the measure implemented here.

The Markov-blanket association is empirical, not proven. A Markov boundary is theoretically relevant to optimal prediction under suitable conditions, but finite samples, model misspecification, correlated variables, redundant features, and non-unique boundaries can distribute importance outside the ideal blanket.

The structural-bit calculation uses a uniform count over labeled DAGs. It measures raw combinatorial uncertainty, not an algorithm’s actual runtime or a realistic prior over market structures.

The DML demonstration is partially linear and i.i.d. DML protects inference against regularization bias in nuisance estimation under its conditions. It does not protect against unmeasured confounding, invalid controls, interference, weak overlap, non-stationarity, or dependent sampling.

The robustness value is not a confounding detector. It quantifies sensitivity within an omitted-variable-bias framework. It does not reveal whether a confounder exists or whether the causal graph is correct.[14]

None of this replaces credible intervention. Where a defensible natural experiment, discontinuity, staggered rollout, randomized test, or instrument exists, use it. Discovery can prioritize hypotheses and expose assumptions; it cannot manufacture exogenous variation.

11. Closing

The trilogy in one line each:

• Part I: Causal discovery and causal inference are different problems. The graph is an output of one and an input to the other.
• Part II: Observational structure and intervention are complementary sources of information. Neither automatically completes the other.
• Part III: Predictive attribution solves a third problem. It can be valuable without being causal.

In this synthetic system, predictive Shapley R² aligned far more closely with the target’s local predictive neighborhood than with total causal effects. It ranked a pure consequence first in all six runs while pushing a major upstream driver toward the bottom.

The lesson is not that prediction is inferior. It is that prediction answers a different question.

Feature importance tells you what a model used to predict. Discovery proposes which arrows may exist. Identification determines whether an effect is recoverable. Estimation tells you its magnitude under those assumptions. Collapse those steps, and a precise, well-validated model can still produce a confidently wrong causal number and in markets, someone may trade on it.

All reported figures come from a synthetic structural causal model with 12 observed variables, one unobserved risk factor, and known ground truth. Rank correlations are averaged across six independent replications with n = 6,000 per run. The companion script uses NumPy, pandas, scikit-learn, and SciPy and reproduces the reported simulation results. Production libraries discussed in the workflow include EconML, DoWhy, causal-learn, Tigramite, and lingam.

• Part I: “Causal Discovery Is Not Causal Inference.”
• Part II: “Causal Modeling: When to Discover, When to Infer.”

References

[1] Pearl, J. (2009). Causality: Models, Reasoning, and Inference (2nd ed.). Cambridge University Press.

[2] Spirtes, P., Glymour, C., & Scheines, R. (2000). Causation, Prediction, and Search (2nd ed.). MIT Press.

[3] Lundberg, S. M., & Lee, S.-I. (2017). “A Unified Approach to Interpreting Model Predictions.” Advances in Neural Information Processing Systems, 30.

[4] Covert, I., Lundberg, S. M., & Lee, S.-I. (2020). “Understanding Global Feature Contributions with Additive Importance Measures.” Advances in Neural Information Processing Systems, 33, 17212–17223.

[5] Koller, D., & Sahami, M. (1996). “Toward Optimal Feature Selection.” Proceedings of the 13th International Conference on Machine Learning.

[6] Colombo, D., Maathuis, M. H., Kalisch, M., & Richardson, T. S. (2012). “Learning High-Dimensional Directed Acyclic Graphs with Latent and Selection Variables.” The Annals of Statistics, 40(1), 294–321.

[7] Shimizu, S., Hoyer, P. O., Hyvärinen, A., & Kerminen, A. (2006). “A Linear Non-Gaussian Acyclic Model for Causal Discovery.” Journal of Machine Learning Research, 7, 2003–2030.

[8] Zheng, X., Aragam, B., Ravikumar, P., & Xing, E. P. (2018). “DAGs with NO TEARS: Continuous Optimization for Structure Learning.” Advances in Neural Information Processing Systems, 31.

[9] Zheng, X., Dan, C., Aragam, B., Ravikumar, P., & Xing, E. P. (2020). “Learning Sparse Nonparametric DAGs.” Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics.

[10] Runge, J., Nowack, P., Kretschmer, M., Flaxman, S., & Sejdinovic, D. (2019). “Detecting and Quantifying Causal Associations in Large Nonlinear Time Series Datasets.” Science Advances, 5(11), eaau4996.

[11] Sharma, A., & Kiciman, E. (2020). “DoWhy: An End-to-End Library for Causal Inference.” arXiv:2011.04216.

[12] Battocchi, K., Dillon, E., Hei, M., Lewis, G., Oka, P., Oprescu, M., & Syrgkanis, V. (2019). “EconML: A Python Package for ML-Based Heterogeneous Treatment Effects Estimation.” Version 0.x.

[13] Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). “Double/Debiased Machine Learning for Treatment and Structural Parameters.” The Econometrics Journal, 21(1), C1–C68.

[14] Cinelli, C., & Hazlett, C. (2020). “Making Sense of Sensitivity: Extending Omitted Variable Bias.” Journal of the Royal Statistical Society: Series B, 82(1), 39–67.

About the Author

Brian Curry is a Kansas City–based machine learning and AI researcher, data scientist, and the creator of Vector1 Research, a personal open-source research initiative working at the intersection of causal inference, complex systems, cognitive AI architecture, and economic systems.

He has more than 20 years of experience across analytics, digital product, marketing, and technology, including applied data science and machine learning work at Koch Industries, Tractor Supply Company, Hallmark, Garmin, and AT&T. His practitioner work focuses on production intelligence systems: causal marketing analytics, simulation, agentic systems, and the cloud infrastructure required to deploy them.

His research and open-source work includes:

• PyCausalSim: a Python framework for causal discovery and inference through simulation.
• Papilon: a framework for causal media modeling, complex systems analysis, optimization, and simulation.
• Daedalus: a platform for causal and economic analysis.
• Memory-Node Encapsulation and NeoCortex-M: architectures for memory-driven artificial intelligence.
• The Autonomy Economy: a research framework for the macroeconomics of advanced automation.