32 CDMs: The Unified Framework

Status: Draft

v0.4

32.1 Introduction

Part IV — Synthesis

Stratum. All three modelling roles, as one generative object (CDM) and as practice.

Levels. Revisit L1–L3 on a single CDM; case studies; then reporting and software as the practice of Reason (identify → estimate → display).

Deferred. Nothing architectural—competing-risks survival and other specialised estimators remain thin by design.

How to read. Chapter 28 is the spine made explicit; 28b applies it; reporting and software show how declared relations become computation without reinventing estimators. Formal structure conditions becoming; it is not an additional efficient agent. The glossary supplies the concise philosophical context.

Structure, dynamics, and observation meet in the book’s central modelling object: the Causal Dynamical Model (CDM). A CDM combines structural causal models with dynamical systems (ODEs, SDEs, networks) and state-space inference, with explicit SCM semantics for interventions and counterfactuals (Pearl 2009; Durbin and Koopman 2012). The framing applies to any complex dynamical system with structure, time, and partial observation; the executable examples use the Julia stack introduced in the preface.

Interpretive note. Edges record direct graphical relevance; structural assignments encode the causal mechanisms, including higher-order interactions. Shared exogenous draws are a declared cross-world coupling for counterfactual simulation, not an empirical definition of identity or a one-to-one translation of creative advance. Formal constraints may condition the admissible model class, but they are not additional event-like parents unless a mechanism and intervention semantics are supplied.

The CDM unifies:

All 3 modelling strata: Structural (Chapters 1–9), Dynamical (Chapters 10–18), Observable (observation and measurement; Chapters 19–27)
All 3 levels of Reason: Seeing (Chapters 1, 3, 10–12, 19–21), Doing (Chapters 4–6, 13–15, 22–24), Imagining (Chapters 7–9, 16–18, 25–27)
Practice of Reason (this Part): reporting standards and software patterns (Chapters 29–30) as identify → estimate → display

32.2 CDMs as a Unified Generative Model

At a high level, a CDM is a single generative model that combines:

a structural component (mechanisms and intervention semantics),
a dynamical component (state evolution through time),
an observation component (how data are generated from latent state),
and exogenous variables (unmodelled variation and noise).

\[ \begin{aligned} X_1 &\coloneqq f_1(C, U^x_1) \\ A_t &\coloneqq \pi(H_t, C, U^a_t) \\ X_{t+1} &\coloneqq f(X_t, A_t, C, U^x_{t+1}) \quad \text{for } t \geq 1 \\ Y_t &\coloneqq h(X_t, C, U^y_t) \end{aligned} \]

Here \(X_t\) is a realised (possibly latent) state, \(Y_t\) is a record, and \(A_t\) is an endogenous action when present. A measurement interaction and design may be made explicit where recording perturbs the process; the direct observation model is its non-perturbative special case. The temporal support, descriptive scale, identity criterion, and role of each object must be declared before interpreting the equations.

The three roles classify components of a generative model; Pearl’s levels classify queries. Their cross-tabulation is a coverage matrix, not an ontological product. Identification, estimation, design, and criticism operate on the model-query pair \((\mathcal M,Q)\) and on the observation assumptions that connect it to data.

32.2.1 Semantic declarations

For every CDM, record the referent, value type, temporal support, descriptive scale, identity criterion where endurance is claimed, and model role of each variable. For every relation, record whether it is causal influence, constitution, participation, measurement, or temporal order. For every intervention, record its target, replacement rule, timing, duration, information set, admissibility and solution conditions, and the mechanisms assumed unchanged. This metadata is part of the model specification, not optional philosophical commentary.

32.3 Synthesis of All 27 Chapters

32.3.1 Part I: Structural Stratum (Chapters 1-9)

Phase 1: Seeing in Structural (Chapters 1-3)

Chapter 1: The Causal Hierarchy and Three Strata, foundational framework
Chapter 2: The Primary Unit: The Dyad, fundamental unit of directed dependence
Chapter 3: Graph Theory and Causal Patterns, what structure looks like

Phase 2: Doing in Structural (Chapters 4-6)

Chapter 4: Structural Causal Models as Executable Mechanisms, how to modify structure
Chapter 5: Identification: When Can We Learn from Data?, what can be learned
Chapter 6: Do-Calculus: Rules for Interventions, symbolic rules for interventions

Phase 3: Imagining in Structural (Chapters 7-9)

Chapter 7: Counterfactuals: Unit-Level Alternatives at Structural Level, alternative structural configurations
Chapter 8: Transportability: Generalising Structural Claims, generalising across domains
Chapter 9: From Structure to Time: FEP and Attractors, transition to Dynamical

CDM Integration: CDMs incorporate Structural stratum concepts through graph structure \(G\) (from Chapters 1-3), SCM semantics (from Chapters 4-6), and counterfactual/transportability reasoning (from Chapters 7-8).

32.3.2 Part II: Dynamical Stratum (Chapters 10-18)

Phase 4: Seeing in Dynamical (Chapters 10-12)

Chapter 10: Deterministic Dynamics: ODEs as Causal Processes, observing deterministic dynamics
Chapter 11: Stochastic Dynamics: SDEs and Random Processes, observing stochastic dynamics
Chapter 12: State-Space Models: Inferring Structure from Observations, inferring dynamics from observations

Phase 5: Doing in Dynamical (Chapters 13-15)

Chapter 13: Intervening in Deterministic Systems, modifying deterministic dynamics
Chapter 14: Intervening in Stochastic Systems, modifying stochastic dynamics
Chapter 15: Advanced Dynamics: Regimes, Networks, and Resilience, advanced intervention topics

Phase 6: Imagining in Dynamical (Chapters 16-18)

Chapter 16: Counterfactual Dynamics: Alternative Trajectories, alternative dynamic trajectories
Chapter 17: Sensitivity Analysis and Robustness in Dynamics, sensitivity to assumptions
Chapter 18: From Dynamical to Observable: Records and Actualisation, transition to the Observable stratum

CDM Integration: CDMs incorporate Dynamical stratum concepts through latent process dynamics \(f(\cdot)\) (from Chapters 10-11, 13-15), state-space inference (from Chapter 12), and counterfactual dynamics (from Chapter 16).

32.3.3 Part III: Observable Stratum (Chapters 19–27)

Phase 7: Seeing in the Observable stratum (Chapters 19–21)

Chapter 19: Observational Methods: Learning from Data, G-methods and IPTW
Chapter 20: TMLE and Doubly Robust Estimation, robust estimation
Chapter 21: Model Validation with Observable Data, validating inferences

Phase 8: Doing in Measurement (Chapters 22-24)

Chapter 22: Interventional Reasoning: Forecasting Under Interventions, forecasting under interventions
Chapter 23: Policy Evaluation and Dynamic Treatment Strategies, evaluating policies
Chapter 24: Causal Decision-Making, making decisions with causal models

Phase 9: Imagining in Measurement (Chapters 25-27)

Chapter 25: Counterfactual Reasoning: Unit-Level Alternatives, alternative observable outcomes
Chapter 26: Hypothesis Generation from Counterfactuals, generating hypotheses
Chapter 27: Experimental Design: Optimal Measurements, designing optimal studies

CDM Integration: CDMs incorporate Measurement-role concepts through observation model \(h(\cdot)\) (from Chapters 19-21), interventional reasoning (from Chapters 22-24), and counterfactual reasoning leading to study design (from Chapters 25-27).

32.3.4 How CDMs Unify All Phases

CDMs provide a single framework that supports:

Seeing (Level 1) across all strata:

Structural: Graph structure, identification (Chapters 1-3, 5)
Dynamical: ODEs, SDEs, state-space inference (Chapters 10-12)
Observable: Observational methods, validation (Chapters 19-21)

Doing (Level 2) across all strata:

Structural: SCMs, do-calculus (Chapters 4-6)
Dynamical: Interventions in ODEs/SDEs (Chapters 13-15)
Observable: Interventional forecasting, policy evaluation, decision-making (Chapters 22-24)

Imagining (Level 3) across all strata:

Structural: Counterfactuals, transportability (Chapters 7-8)
Dynamical: Counterfactual dynamics, sensitivity (Chapters 16-17)
Observable: Counterfactual reasoning, hypothesis generation, study design (Chapters 25-27)

Exogenous variables \(U_t\) represent unmodelled variation and noise. Interventions \(do(A_t = a)\) modify the generating process, and counterfactual reasoning compares alternative outcomes for the same unit by holding fixed its exogenous realisation \(\mathbf{u}\).

32.4 What Is a CDM?

32.4.1 Definition

A CDM is a causal model with dynamical systems structure that includes:

Graph structure (when applicable): \(G\) encoding directed dependencies between variables
Latent process: \(X_t\) with dynamics \(f(\cdot)\)
Observation model: \(Y_t\) with measurement \(h(\cdot)\)
Explicit noise: Exogenous variables \(U_t\) (SCM semantics)
Intervention operators: \(do(\cdot)\) for interventions that modify mechanisms
Counterfactual semantics: An explicitly declared cross-world coupling, commonly shared \(\mathbf{u}\) in SCM examples

32.4.2 Formal Structure

\[ \begin{aligned} X_1 &\coloneqq f_1(C, U^x_1) \quad \text{(initial state)} \\ A_t &\coloneqq \pi(H_t, C, U^a_t) \quad \text{(optional behaviour policy)} \\ X_{t+1} &\coloneqq f(X_t, A_t, C, U^x_{t+1}) \quad \text{for } t \geq 1 \\ Y_t &\coloneqq h(X_t, C, U^y_t) \end{aligned} \]

32.5 Notation and Modelling Conventions

This section introduces the book’s unified notation for Causal Dynamical Models (CDMs), distinguishing endogenous action policies from interventions that transform a model.

32.5.1 Causal Dynamical Model (CDM)

A CDM is an SCM (Pearl 2009) whose endogenous variables are time-indexed and partitioned into latent process variables and records, with explicit intervention operators, measurement semantics, and counterfactual couplings.

Alternative name: Causal State-Space Model (CSSM) when emphasising state-space inference; DSCM (Dynamic SCM) when emphasising structural equations. CDM is an integrative methodology, not a claim to a wholly new mathematical class.

32.5.2 Core Object

A canonical CDM is represented by the tuple: \[ \mathcal{M} = \big(G,F,P_U,\pi,g,h,\mathcal{I},\mathcal{S}\big) \]

where:

\(G\): The declared graph or network structure. DAG semantics are used by the acyclic identification routines; cyclic, equilibrium and topology-changing models must state separate validity conditions.
\(P_U\): The joint law over exogenous variables representing unmodelled generative influence and stochastic variation. A shared-\(U\) coupling is a modelling choice for L3, not an empirically identified metaphysical claim. These include:
\(U^x_t\): Process noise for state dynamics
\(U^y_t\): Observation noise for measurements (measurement error)
\(U^a_t\): Action/treatment noise (only present if \(\mathbf{A}_t\) is generated by a stochastic behavioural policy \(\pi\); if \(\mathbf{A}_t\) is purely set via intervention \(do(\mathbf{A}_t = \mathbf{a}_t)\), there is no noise)
\(U^1\): Initial state noise (for the first time step)

The exogenous variables make stochasticity explicit and structural, enabling Pearl’s causal semantics for interventions and counterfactuals.

\(F\): Structural assignments and transition mechanisms. Each assignment has the form: \[ X_i \coloneqq f_i(\text{Pa}(X_i), U_i) \] where \(\text{Pa}(X_i)\) are the parents of \(X_i\) in \(G\), and \(f_i\) specifies how parents and noise combine to generate \(X_i\). For time-indexed systems, \(F\) includes:
\(f_1\): Initial state assignment
\(f\): State transition function
\(h\): Record-generating function (how latent states become records)
\(\pi\): An optional endogenous action policy, distinct from an investigator’s intervention.
\(g\): An optional measurement-interaction map, used when recording perturbs the latent process.
\(\mathcal{I}\): The declared intervention domain: targets, replacements, timing, scope and cointerventions.
\(\mathcal{S}\): Scale, regime, solver and validity metadata.

Structural constraints can be recorded alongside, but not folded into, this causal model. Let \(\mathcal K\) be a registry of claims about invariance, feasibility, viability or cross-embodiment relations. A member of \(\mathcal K\) names the target components of \(\mathcal M\), predicted consequences, and proposed tests; it is neither a parent in \(G\) nor an additional exogenous input. It records formal causation as a condition on admissible unfolding, while leaving efficient temporal influence to the dynamical model. This separation makes formal explanation auditable without silently turning a constraint into an event-like cause.

The joint law \(P_U\) is typically factorised as: \[ P_U = P(U^1) \prod_{t=1}^{T} P(U^x_t) P(U^y_t) P(U^a_t) \] (where some factors may be absent if certain variables are deterministic). Interventions modify declared assignments or regimes; they do not silently modify \(P_U\).

Together these components support association, declared interventions and counterfactual reasoning. The tuple is a modelling contract; it does not by itself establish identification, solver validity, organismal identity or a process-metaphysical interpretation.

32.5.3 Declaring structural constraints without reifying them

Consider a regeneration study in which damage may be applied to tissue and a morphological read-out is repeatedly recorded. A viability claim says that trajectories under a specified family of damage interventions can return to a bounded region of state space. It does not say that the region is a causal parent, or that the claim alone explains the recovery. The relevant mechanism, boundary conditions, assay and intervention must still be stated and tested (Chapter 9).

StructuralConstraintSpec records this narrow kind of claim. Its certificate is stable metadata for an analysis record; it does not validate the biological claim or introduce a new node into a DAG.

project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))
@auto_using CausalDynamics

regeneration_constraint = StructuralConstraintSpec(
    :regenerative_morphology,
    [:bioelectric_state, :morphology];
    kind = :viability,
    claim = "Perturbed tissue can return to a viable morphological region.",
    predictions = [:return_after_damage],
    test_ids = [:damage_recovery_assay],
    embodiments = [:planarian, :anthrobot],
)
regeneration_certificate = constraint_certificate(regeneration_constraint)
nothing

The declaration can be paired with a MechanismLibrary for a stated transition mechanism, a CausalAbstractionSpec for a finite interventional comparison across scales, and typed intervention descriptors for the damage protocol. None of these steps establishes the metaphysical status of a modelled pattern; together they make competing scientific accounts more explicit and discriminable.

The small recurrent-infection benchmark in scripts/recurrent_infection_benchmark.jl applies the same discipline to a decision problem. It is a deliberately synthetic teaching model, not a disease claim: burden and toxicity unfold as occasions, a noisy and intermittent measurement process supplies the information available to treatment, and a policy is represented as an enduring determination whose outputs may vary with history. The temporal graph can be inspected with recurrent_infection_temporal_spec() and its constitutive and recurrent edges with temporal_edge_records. This makes the benchmark suitable for checking description, intervention simulation and policy comparison before introducing real clinical data. recurrent_infection_loss supplies a transparent additive loss for those comparisons; its weights are teaching parameters, not clinical utility claims.

32.5.4 Decision comparison: recurrent treatment policies

The benchmark makes the decision question explicit: which policy minimises a declared loss under the same realised random-number stream? This is a Level 2 simulation comparison, not an estimate of a clinical effect. The policy rule is held fixed within each run, while its treatment output can vary with the observed history. The resulting summary is shown in Figure 32.1.

project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "recurrent_infection_benchmark.jl"))

benchmark_parameters = RecurrentInfectionParameters(
    horizon = 24,
    monitoring_interval = 4,
)
policy_comparison = compare_recurrent_policies(
    benchmark_parameters;
    policies = [:calendar, :threshold, :sparse_threshold, :risk_adaptive],
    seed = 20260911,
)
policy_comparison

4-element Vector{NamedTuple}:
 (policy = :calendar, loss = 24.380934137185374, mean_burden = 0.9063490093561941, treatment_total = 5.0, monitoring_total = 6.0)
 (policy = :threshold, loss = 22.29136868434853, mean_burden = 0.8379176788774346, treatment_total = 2.0, monitoring_total = 24.0)
 (policy = :sparse_threshold, loss = 21.432565066141148, mean_burden = 0.818257971795151, treatment_total = 3.0, monitoring_total = 6.0)
 (policy = :risk_adaptive, loss = 21.432565066141148, mean_burden = 0.818257971795151, treatment_total = 3.0, monitoring_total = 6.0)

Figure 32.1: Synthetic recurrent-infection policy comparison. Burden trajectories use common random numbers; the loss panel combines burden, toxicity, treatment, and monitoring costs with teaching weights.

The figure matters because it does not settle the biological interpretation of the residual variation. It shows how a policy comparison depends on the observation process, the declared loss, and the intervention semantics. Changing those assumptions is a new analysis, not a silent update to the same enduring policy.

32.5.5 Decision quality: selection, reference value, and regret

A policy comparison becomes a decision analysis only when the policy selected from one set of simulated histories is evaluated on another. The benchmark therefore separates an evaluation stream, used to select the policy, from an independent reference stream, used to estimate the loss of every candidate. The reference-best policy is only the best policy in the declared candidate class. It is not a clinical optimum, and the resulting regret is a simulator diagnostic rather than an empirical treatment effect.

decision_benchmark = recurrent_infection_policy_regret(
    benchmark_parameters;
    policies = [:calendar, :threshold, :sparse_threshold, :risk_adaptive],
    evaluation_n = 200,
    reference_n = 1_000,
    seed = 20260911,
)
(
    selected_policy = decision_benchmark.selected_policy,
    reference_best_policy = decision_benchmark.reference_best_policy,
    regret = round(decision_benchmark.regret; digits = 3),
    reference_losses = [
        (policy = record.policy, mean_loss = round(record.mean_loss; digits = 3))
        for record in decision_benchmark.reference
    ],
)

(selected_policy = :calendar, reference_best_policy = :sparse_threshold, regret = 0.061, reference_losses = [(policy = :calendar, mean_loss = 25.59), (policy = :threshold, mean_loss = 26.908), (policy = :sparse_threshold, mean_loss = 25.529), (policy = :risk_adaptive, mean_loss = 25.679)])

With these seeds, the finite evaluation sample selects :calendar, whereas the independent reference stream favours :sparse_threshold. That difference is not evidence that one method is universally superior: it exposes the decision consequence of finite-sample uncertainty in a declared policy class. Increasing the evaluation sample, changing the loss weights, or improving the measurement process are substantive analyses and should be reported as such.

32.5.6 From generated histories to estimation

The same fixture can produce a wide panel for CausalTargeted. The panel keeps the natural treatment history, measured time-varying covariates, and terminal outcome in separate columns. CausalTargeted then receives an explicit role mapping rather than guessing which columns constitute the treatment history.

project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "recurrent_infection_benchmark.jl"))
include(joinpath(project_root, "scripts", "ensure_packages.jl"))
@auto_using CausalTargeted StableRNGs
using Statistics

panel_parameters = RecurrentInfectionParameters(
    horizon = 4,
    baseline_sd = 1.2,
    process_sd = 0.1,
)
estimation_panel = recurrent_infection_estimation_panel(
    panel_parameters;
    n = 120,
    policy = :sparse_threshold,
    rng = StableRNG(20260911),
)
panel_data = estimation_panel.data
panel_spec = estimation_panel.spec
always_treat = CausalTargeted.discrete_static_policy("1"; levels = ["0", "1"])
panel_estimate = CausalTargeted.run_sequential_lmtp(
    panel_data,
    panel_spec.treatments,
    panel_spec.outcome;
    baseline = panel_spec.baseline,
    time_vary = panel_spec.time_vary,
    policies = always_treat,
    folds = 2,
    learners = (:glm_interact, :mean),
    rng = StableRNG(20260912),
)
(;
    estimate = panel_estimate.estimate,
    standard_error = panel_estimate.se,
    treatment_times = panel_estimate.times,
    n = panel_estimate.n,
    positivity_ok = panel_estimate.positivity.ok,
)

(estimate = 0.04228371995952935, standard_error = 0.0035306061750225155, treatment_times = 3, n = 120, positivity_ok = true)

This estimate is a deliberately modest hand-off test, not a clinical claim. The simulator supplies a known generative process, whereas the estimator sees only the generated panel. The target is explicit: always treat at occasions \(2{:}T\), compared with the natural sparse-threshold policy. The panel’s initial occasion is baseline rather than a degenerate treatment opportunity. The interaction learner is declared because treatment changes the burden-growth term; using a main-effects learner here is an intentional misspecification sensitivity analysis, not the benchmark default.

The following small replication compares the full role specification with a declared ablation that omits the measured burden history while retaining the monitoring and toxicity histories. The oracle is generated by direct intervention in the simulator and is shared across analyses; a failed positivity check remains visible in the summary. This is a smoke test for the rendered example, not a claim that the full specification must win in every finite sample.

validation = run_recurrent_infection_validation(
    panel_parameters;
    analyses = [:full, :omit_measurement_history],
    replicates = 3,
    n = 80,
    oracle_n = 500,
    seed = 20260916,
)
summarise_recurrent_infection_validation(validation)

2-element Vector{NamedTuple}:
 (analysis = :full, replicates = 3, oracle = 0.04318105424786279, mean_estimate = 0.04502479584366035, mean_bias = 0.0018437415957975646, rmse = 0.025993274504781762, interval_coverage = 0.3333333333333333, positivity_rate = 1.0)
 (analysis = :omit_measurement_history, replicates = 3, oracle = 0.04318105424786279, mean_estimate = 0.04837865601419874, mean_bias = 0.005197601766335952, rmse = 0.012905532475833325, interval_coverage = 0.0, positivity_rate = 1.0)

Three replications are only a rendered smoke test. A real validation run should increase both the number of replications and the oracle sample, preregister the policy and loss, and report results separately for positivity-valid and positivity-invalid runs. The ablation is an information comparison, not a metaphysical test: a residual or a change in finite-sample bias does not by itself identify creative advance or formal causation.

32.5.7 A strong conventional comparator

The benchmark also compares the CDCS representation with two analysis paths that use the same generated panel. The sequential LMTP is a strong established longitudinal causal comparator: it receives the full observed treatment history, measured burden history, toxicity and monitoring history. The terminal regression is a plausible conventional shortcut: it predicts terminal burden from baseline burden and cumulative treatment, thereby collapsing the temporal feedback into one summary. It is useful as an ablation, not as a causal gold standard.

The simulator oracle supplies the same always-treat reference to both methods. Agreement between the full-history LMTP and the oracle is therefore a correctness and interoperability check. A larger error from the terminal regression shows the consequence of omitting the dynamical history; it does not show that CDCS has made a mathematically impossible analysis possible.

estimator_comparison = run_recurrent_infection_estimator_comparison(
    panel_parameters;
    replicates = 3,
    n = 120,
    oracle_n = 500,
    seed = 20260918,
)
summarise_recurrent_infection_estimator_comparison(estimator_comparison)

2-element Vector{NamedTuple}:
 (method = :sequential_lmtp, replicates = 3, oracle = 0.04683171630447089, mean_estimate = 0.04575973490545155, mean_bias = -0.001071981399019338, mean_absolute_bias = 0.003367554551226594, rmse = 0.003520358385312785)
 (method = :terminal_regression, replicates = 3, oracle = 0.04683171630447089, mean_estimate = 0.00044347050463609724, mean_bias = -0.04638824579983478, mean_absolute_bias = 0.04638824579983478, rmse = 0.046663601002798864)

The comparison is deliberately asymmetric in interpretation. A correctly specified conventional longitudinal estimator should match the CDCS-linked analysis on this estimand. The practical claim is that the typed history, observation and policy semantics make the relevant information harder to omit, make the omission visible when it occurs, and connect the resulting estimate to a decision-level reference. The next validation stage should add a strong predictive model and a deliberately delayed or biased measurement process, then report trajectory calibration, policy value and regret together.

32.5.8 One history, four package hand-offs

The shared fixture in scripts/cdcs_shared_fixture.jl is a small integration contract for the CDCS packages. Every variable concerns one enduring individual referent, declared once through a ReferentSpec. A baseline context is present from onset and feeds an assignment attribute with from-onset support; that assignment contributes to a protocol attribute, while the pointwise outcome is measured at \(t = 2\). Pointwise implementation states receive repeated influence from the protocol. Temporal support fixes node multiplicity (one node for each from-onset attribute, one per represented time for each pointwise state); the referent carries identity; neither declaration says anything about occasions. The data table is the observable projection of this history, not a replacement for it. CausalDynamics supplies the temporal unrolling, CausalTargeted evaluates a policy contrast, CausalMediation receives the same treatment and outcome data with an explicit protocol mediator, and DAGMakie renders the single-support nodes without duplicating their identities.

include(joinpath(project_root, "scripts", "cdcs_shared_fixture.jl"))

shared_fixture = build_cdcs_shared_fixture(
    n = 160,
    horizon = 3,
    seed = 20260930,
)
shared_mediation = CausalMediation.run_mediation(
    shared_fixture.mediation_spec,
    shared_fixture.data;
    folds = 2,
    n_mc = 4,
    learners = (:glm, :mean),
    parallel = false,
    rng = StableRNG(20260931),
)
shared_decomposition = CausalMediation.decompose(shared_mediation)
shared_figure, shared_axis, shared_plot = render_cdcs_shared_fixture(shared_fixture)
shared_single_node_count = count(
    key -> key[2] === nothing,
    shared_fixture.unrolling.index_node,
)
(;
    rows = nrow(shared_fixture.data),
    single_support_nodes = shared_single_node_count,
    policy_estimate = shared_fixture.policy_estimate.estimate,
    mediation_total_effect = shared_decomposition.te,
    graph_rendered = true,
)

(rows = 160, single_support_nodes = 4, policy_estimate = 0.5067549093391777, mediation_total_effect = 0.747230430004248, graph_rendered = true)

This is an integration check rather than a claim that the synthetic mediation effect is identified in an applied study. The result is that the same policy and protocol identities survive the transitions between packages. The focused cross-package test also checks the glyph contract: the two from-onset attributes are single nodes but are drawn as circles, because DAGMakie’s rounded rectangle marks value_representation = :interval_summary and nothing else. Neither a single-node support nor a claim of endurance changes a node’s shape.

32.5.9 Supported semantics in this release

Capability	Status	Meaning
Discrete state, policy and counterfactual simulation	Executable	Typed actions, interventions and declared couplings are implemented
ODE pins, initial conditions, forces and RHS replacement	Executable	Supported through the continuous CDM façade with solver assumptions
Topology and model-class transitions	Declarative contract	Described and fingerprinted; executable state mapping remains future work
SDE drift, diffusion and stopping-time surgery	Declarative contract	Intervention types are documented; general solver transformation remains future work
Causal abstraction	Exact finite-law validation	Checks \(τ_{\#}P_m^i = P_M^{\omega(i)}\) over a declared intervention domain, with explicit approximate mode
Structural constraint declaration	Auditable metadata	Records a formal claim, predictions and proposed tests; does not create a DAG node or validate the claim
Cyclic SCM identification	Separate semantics	Chapter 9b documents sigma/equilibrium conditions; the main identification API remains DAG-based
FEP and active inference	Optional module	Requires an explicit generative model, objective and empirical comparison

32.5.10 Variable Sets (Time-Indexed)

For \(t = 1,2,\dots,T\):

Latent/process state: \(\mathbf{X}_t \in \mathbb{R}^d\)
Action/treatment variable (possibly empty): \(\mathbf{A}_t \in \mathcal{A}\), a variable that can be set via intervention \(do(\mathbf{A}_t = \mathbf{a}_t)\) or generated by a policy \(\pi(\mathbf{H}_t, \mathbf{C}, \mathbf{U}^a_t)\)
Observation: \(\mathbf{Y}_t \in \mathcal{Y}\)
Optional context/domain variables (time-invariant): \(\mathbf{C}\)

Exogenous inputs:

\(\mathbf{U}^x_t\) for state dynamics
\(\mathbf{U}^y_t\) for observation/measurement
\(\mathbf{U}^a_t\) for action/treatment assignment (if \(\mathbf{A}_t\) is generated by a stochastic policy rather than purely set via intervention)

Action vs Intervention: Terminology Note

\(\mathbf{A}_t\) is the “action/treatment variable”, a variable that can be either:

Set via intervention: \(do(\mathbf{A}_t = \mathbf{a}_t)\) (an experimenter sets it to a specific value)
Generated by a policy: \(\mathbf{A}_t \coloneqq \pi(\mathbf{H}_t, \mathbf{C}, \mathbf{U}^a_t)\) (a behavioural policy generates it, possibly stochastically)

“Intervention” is the broader concept that includes:

Action/treatment interventions: \(do(\mathbf{A}_t = \mathbf{a}_t)\), setting the action variable
Mechanism interventions: \(do(f \leftarrow f^\star)\), modifying the state transition function
Parameter interventions: \(do(\theta \leftarrow \theta^\star)\), modifying model parameters

Why “action/treatment”? We use both terms to bridge communities: “action” is standard in reinforcement learning and control theory, while “treatment” is standard in causal inference and epidemiology. Both refer to the same variable \(\mathbf{A}_t\) that can be intervened upon.

Action noise \(U^a_t\): This is only present when \(\mathbf{A}_t\) is generated by a stochastic policy. If \(\mathbf{A}_t\) is purely set via intervention (no policy), there is no \(U^a_t\) noise.

32.5.11 Structural Assignments (Markovian Baseline)

A compact default form: \[ \begin{aligned} \mathbf{X}_1 &\coloneqq f_1(\mathbf{C}, \mathbf{U}^x_1) \quad \text{(initial state)} \\ \mathbf{A}_t &\coloneqq \pi(\mathbf{H}_t, \mathbf{C}, \mathbf{U}^a_t) \quad \text{(optional; behavioural policy generating actions)}\\ \mathbf{X}_{t+1} &\coloneqq f(\mathbf{X}_t, \mathbf{A}_t, \mathbf{C}, \mathbf{U}^x_{t+1}) \quad \text{for } t \geq 1 \\ \mathbf{M}_{t+1} &\coloneqq g(\mathbf{M}_t, \mathbf{X}_t, \mathbf{D}_t, \mathbf{U}^m_{t+1}), \qquad \mathbf{Y}_t \coloneqq h(\mathbf{M}_t, \mathbf{U}^y_t) \end{aligned} \] where \(\mathbf{H}_t \coloneqq (\mathbf{Y}_{1:t}, \mathbf{A}_{1:t-1})\) is the observed history. Omit \(M_t\) when a direct record model \(Y_t=h(X_t,\ldots)\) is adequate.

This covers:

ODE/SDE discretisations (via \(f\))
Switching models (include a discrete latent \(\mathbf{S}_t\) inside \(\mathbf{X}_t\))
Networked systems (define components per node; see below)

32.5.12 Beyond the Markovian baseline: history and attention

The baseline above is Markovian in state: \(\mathbf{X}_{t+1}\) depends only on \(\mathbf{X}_t\) (plus \(\mathbf{A}_t\), \(\mathbf{C}\), noise). That is often justified when the current state is a sufficient statistic for the future, it summarises everything from the past that matters. In many systems, however, the past can matter in ways that are not fully compressed into \(\mathbf{X}_t\): long delays, distributed memory, or dependencies that do not have a compact state representation.

Sequence models and attention. Large language models (LLMs) and related architectures (e.g. transformers) make the next output depend on all previous positions in the sequence, via attention: each position can “attend to” earlier positions with learned weights (Vaswani et al. 2017). So the effective dependence is not “only the previous state” but a weighted combination of many past states. In the minimal formulation (e.g. microgpt by Andrej Karpathy, a dependency-free Python GPT; a Julia version with a dynamics demo is in scripts/microgpt.jl), the next-token distribution is a function of the full preceding context, with attention determining which parts of that context matter.

One way to interpret this in dynamical-systems terms is learned dependency selection over history: instead of committing to a fixed lag structure, the model learns which past states are predictive in the current context. The Markovian baseline corresponds to dependence on the immediate predecessor; fixed-delay dynamics (e.g. delay differential equations, Stochastic dynamics) correspond to dependence on one or a few fixed lags; attention corresponds to learned weights over many lags.

Implications for causal dynamics.

When Markov is enough: If the process has a compact sufficient statistic (e.g. many dynamical systems, well-specified state-space models), the Markovian form \(\mathbf{X}_{t+1} = f(\mathbf{X}_t, \ldots)\) is appropriate and keeps inference and intervention tractable.
When history matters explicitly: If the system has long-range dependencies, delayed feedback, or no compact state, we may need either (i) richer state (e.g. including lagged values or a memory buffer) so that the expanded state is again Markovian, or (ii) explicit history dependence in the transition, e.g. \(\mathbf{X}_{t+1} = f(\mathbf{X}_t, \mathbf{X}_{t-1}, \ldots, \mathbf{X}_{t-\ell}, \mathbf{A}_t, \ldots)\) or delay differential equations with fixed lags, or (iii) attention-like mechanisms where the next state or output is a learned function of a weighted combination of past states, useful when the relevant past is not a fixed window but varies with context.

Why go beyond Markovian assumptions? The literature on state-space models, POMDPs, and causal time series gives several well-established reasons (attention is one way to accommodate them; others include belief-state filters, higher-order Markov models, or expanding the state). Partial observability is the most direct: we typically observe \(\mathbf{Y}_t\), not the latent state \(\mathbf{X}_t\). The current observation is then rarely a sufficient statistic for the state; optimal prediction and control are history-dependent (e.g. in POMDPs, policies map observation histories or belief states to actions) (e.g. Koller and Friedman 2009). Incomplete or aggregated state: if causally relevant state variables are missing or collapsed (e.g. through dimension reduction), the one-step transition property can fail and multi-step dependencies appear (e.g. Araujo et al. 2025). Violations of independence of exogenous variables: causal and state-space models may assume exogenous inputs \(\mathbf{U}^x_t\), \(\mathbf{U}^y_t\) are independent across time, or may impose another known dependence structure. When unmeasured confounders persist or process innovations are serially correlated, the Markov state may no longer be sufficient; conditioning on a longer history can absorb that residual dependence. Observation model misspecification or temporally dependent observation noise: if the observation process is wrong or observation noise is non-Markov, the single-step summary of the past is inadequate; history-dependent prediction can partly compensate. Nonstationarity: when causal strengths or input distributions change over time, the invariant one-step Markov representation may break down (e.g. Shen et al. 2019). In all these cases, letting the next-step distribution depend on a learned weighting of the past (e.g. via attention) relaxes the demand that one step compresses everything that matters, without assuming the underlying violations away.

32.5.13 Computationally efficient alternatives

Given this book’s focus on causal dynamics and state-space inference, several options are often more efficient than full attention over raw history.

When the latent process is Markov but we have partial observability, the standard and efficient solution is belief-state filtering: maintain \(P(\mathbf{X}_t \mid \mathbf{Y}_{1:t})\) as the sufficient statistic for the observation history. The Kalman filter updates this in \(O(d^2)\) per step (state dimension \(d\)); particle filters in \(O(N)\) per step for \(N\) particles. No storage or attention over the full history is needed (the belief state is the compressed history). This is the core of the state-space inference material in the book (State-space models).

When the process itself is not first-order Markov (e.g. higher-order dynamics, serially correlated noise), alternatives to \(O(n^2)\) attention include:

Expanding the state to \((\mathbf{X}_t, \mathbf{X}_{t-1}, \ldots, \mathbf{X}_{t-\ell})\) so the expanded process is Markov; cost is \(O(\ell)\) per step
Recurrent state (e.g. RNN/LSTM), which compresses history into a fixed-size hidden state in \(O(1)\) per step but is less interpretable for causal reasoning
Structured sequence models (e.g. linear attention, S4, Mamba), which achieve \(O(n)\) or subquadratic cost for long sequences while retaining a recurrent or convolutional structure

For causal modelling and intervention, belief-state filtering and expanded-state formulations also keep the graph and \(do(\cdot)\) semantics explicit, whereas black-box attention over observations does not.

The policy \(\mathbf{A}_t = \pi(\mathbf{H}_t, \ldots)\) is already history-dependent (it can use the full observed history \(\mathbf{H}_t\)). The question is whether the process dynamics should also depend on more than \(\mathbf{X}_t\). For the bulk of this book we keep the Markovian process model; the ideas above suggest where and how to relax it when the problem demands it. Runnable examples (belief-state filtering and expanded-state (AR(2)) in State-space models, and a Markov-vs-attention comparison on discretised dynamics in Stochastic dynamics) and the script scripts/microgpt.jl (inspired by Andrej Karpathy’s microgpt) are in the relevant chapters and the repo.

32.5.14 Distributions Induced by the CDM

Even though the model is written as assignments, it induces the usual probabilistic factorisation: \[ P(\mathbf{y}_{1:T}) = \int P(\mathbf{u}) \prod_{t=1}^{T} \delta\!\big(\mathbf{x}_t - F_t(\cdot)\big)\; \delta\!\big(\mathbf{y}_t - H_t(\cdot)\big)\; d\mathbf{u}\, d\mathbf{x} \]

For practical work, you usually work with the implied conditional densities (SSM form): \[ P(\mathbf{x}_1)\prod_{t=1}^{T-1} P(\mathbf{x}_{t+1}\mid \mathbf{x}_t,\mathbf{a}_t)\prod_{t=1}^{T} P(\mathbf{y}_t\mid \mathbf{x}_t) \] with the understanding that these are shorthand for structural equations + noise.

32.5.15 Interventions (Pearl Level 2) in CDMs

Use structural interventions explicitly. An endogenous action \(A_t=\pi(H_t,U^a_t)\) is evidence about history; it is not itself \(do(A_t=a)\). CDCS represents an intervention as a typed model transformation \(\iota:\mathcal M\mapsto\mathcal M^\iota\) with target, replacement, interval, scope, stochasticity and cointerventions.

Action intervention: replace the assignment for \(A_t\) by a value or regime.
Policy intervention: replace \(\pi\) by \(\pi^\star\).
Parameter or mechanism intervention: replace \(\theta\) or \(f_i\).
Structural intervention: replace a declared part of \(G\) or a parent set.
Topology intervention: replace a declared graph or connectivity regime; this is a model-class transition, not an ordinary state intervention.
Measurement intervention: replace \(g\), \(h\) or schedule/design \(D_t\).
Environment intervention: replace external inputs or their law.

Topology and model-class changes are first-class declared interventions in the semantic contract. Their executable transformation must specify existence, uniqueness and regime-validity conditions; a descriptor alone does not claim that those conditions hold.

Notation:

Interventional distribution: \[ P^{do(\mathbf{a}_{1:T-1})}_{\mathcal{M}}(\mathbf{y}_{1:T}) \]
Policy intervention (dynamic strategy): \[ do(\pi):\;\mathbf{A}_t \coloneqq \pi(\mathbf{H}_t) \]

32.5.16 Counterfactuals (Pearl Level 3) in CDMs

For a fixed exogenous realisation \(\mathbf{u}\), define the counterfactual trajectory under intervention \(\iota\) (e.g., \(do(\mathbf{A}=\mathbf{a})\) or \(do(\pi)\)): \[ \mathbf{Y}^{\iota}_{1:T}(\mathbf{u}) \]

Population counterfactual quantities average over \(P(\mathbf{U})\), and unit-level counterfactuals condition on evidence \(E=e\) (implemented by inferring a posterior over exogenous/noise or initial state consistent with \(e\)).

32.5.17 Networked CDMs

Let nodes be \(i \in \mathcal{V}\). Node state \(X^i_t\). Parent set \(Pa(i)\) from a directed interaction graph \(G = (V, E)\) where \(V\) is the set of vertices (variables) and \(E\) is the set of directed edges.

Structural dynamics: \[ X^{i}_{t+1} \coloneqq f_i\!\Big(X^{i}_t,\; \{X^{j}_t: j\in Pa(i)\},\; A_t,\; C,\; U^{i}_{t+1}\Big) \]

The graph \(G\) encodes direct dependencies: which nodes directly affect which others. The parent set \(Pa(i)\) lists the direct inputs to node \(i\)’s update rule. The function \(f_i\) specifies how those inputs (plus noise) generate the next state \(X^i_{t+1}\).

This is the clean bridge between:

ecological/physiological interaction networks (Newman 2018; Barabási 2016)
networked dynamical simulation (Barrat et al. 2008; Porter and Gleeson 2016)
causal reasoning about node/edge interventions (remove node, cut edge, change coupling) (Pearl 2009)
Whitehead (1978): Process and Reality: philosophical background (optional)

32.5.18 Quick Notation Reference

Symbol	Meaning	Notes
\(\mathbf{X}_t\)	latent/process state	may include continuous + discrete components
\(\mathbf{Y}_t\)	observation	measurement mechanism is explicit via \(h\)
\(\mathbf{A}_t\)	action/treatment variable	can be set by intervention \(do(\cdot)\) or generated by a behavioural policy \(\pi\)
\(\mathbf{C}\)	context/domain variables	age, site, year; used for transportability
\(\mathbf{U}^x_t,\mathbf{U}^y_t\)	exogenous process and observation variables	make declared process and measurement variability explicit; dependence assumptions remain separate
\(do(\cdot)\)	intervention operator	replaces an assignment/mechanism
\(\pi\)	policy (dynamic intervention rule)	\(\mathbf{A}_t \coloneqq \pi(\mathbf{H}_t)\)
\(\iota\)	generic intervention descriptor	mechanism edit, action set, policy, edge cut
\(\mathbf{Y}^{\iota}(\mathbf{u})\)	counterfactual outcome	same declared background values under a specified coupling

32.5.19 Why This Notation Is Useful

This unified notation:

Keeps Pearl semantics (assignments + interventions + counterfactuals) explicit (Pearl 2009)
Stays compatible with state-space computation (transition/observation densities) (Durbin and Koopman 2012; Särkkä 2013)
Naturally accommodates networks (Newman 2018), switching/regimes (Hamilton 1994), and domain shift (\(\mathbf{C}\)) (Bareinboim and Pearl 2013; Pearl and Bareinboim 2014)
Clarifies a frequent ambiguity in applied work: whether “treatment” is merely conditioned on or structurally set/changes mechanisms (Robins et al. 2000; Hernán and Robins 2020)

32.6 Three Modes of Reasoning

32.6.1 1. Forecasting (Level 1: Association)

Question: “What will happen?”

Answer: Conditional distribution \[ P(Y_{t+1} \mid Y_{1:t}, A_t = a) \]

Use: Predict future observations given past data and observed treatment.

32.6.2 2. Interventional Simulation (Level 2: Intervention)

Question: “What will happen if we do X?”

Answer: Interventional distribution \[ P^{do(A_t = a)}(Y_{t+1} \mid Y_{1:t}) \]

Use: Simulate effects of interventions (treatments, policies, structural changes).

32.6.3 3. Counterfactual Reasoning (Level 3: Counterfactual)

Question: “What would have happened for this unit if X had been different?”

Answer: Counterfactual trajectory \[ Y^{do(A = a)}(\mathbf{u}) \quad \text{for fixed } \mathbf{u} \]

Use: Unit-level causal reasoning, individualised treatment effects.

32.7 Why CDMs Unify

32.7.1 Combines Previous Concepts

CDMs bring together:

SCMs (Part I): Structural assignments, interventions, counterfactuals (Pearl 2009)
SSMs (Part II): Latent process, observation model, inference (Durbin and Koopman 2012; Särkkä 2013)
Dynamical systems (Part III): ODEs, SDEs, switching dynamics (Strogatz 2014; Øksendal 2013; Hamilton 1994)
Networks (Part IV): Networked dynamics, structural interventions (Newman 2018; Barrat et al. 2008)

32.7.2 Single Object, Multiple Queries

One CDM supports:

Inference: Filtering, smoothing, parameter learning (Durbin and Koopman 2012; Särkkä 2013; Doucet et al. 2001)
Forecasting: Conditional predictions
Interventions: Interventional simulations (Pearl 2009)
Counterfactuals: Unit-level reasoning (Pearl 2009; Imbens and Rubin 2015)

32.8 Determining CDM Scope: Markov Boundary and Markov Blanket

When constructing a CDM, a fundamental question is: which variables (nodes) should be included or excluded? This question is answered by two complementary concepts:

Markov boundary (introduced in Graph Theory and Causal Patterns): A local graph-theoretic conditional-dependence summary, not an automatic causal adjustment or counterfactual set
Markov blanket (Chapter 9): Organism (\(\mu\), \(X_t\)) vs environment (\(\eta\)), interface \(s\), \(a\) (\(Y_t\), \(A_t\)), Structural/Dynamical cut formalised by FEP

Together, these concepts provide diagnostic inputs for determining CDM scope; neither is an automatic boundary rule. For detailed discussion of how to apply both concepts when constructing CDMs for decision-making, see Causal Decision-Making.

32.9 CDM Components

32.9.1 1. Process Model

The latent dynamics: \[ X_{t+1} = f(X_t, A_t, C, \theta, U^x_{t+1}) \]

Can be:

ODE (deterministic) (Strogatz 2014; Hirsch et al. 2012)
SDE (stochastic) (Øksendal 2013; Gardiner 2009)
Switching (regime-dependent) (Hamilton 1994; Frühwirth-Schnatter 2006)
Networked (node-level) (Newman 2018; Barrat et al. 2008)

32.9.2 2. Observation Model

How we observe the process: \[ Y_t = h(X_t, C, U^y_t) \]

Can be:

Gaussian (continuous) (Durbin and Koopman 2012)
Poisson (counts) (King et al. 2016)
Bernoulli (binary)
Compositional (proportions) (Aitchison 1986)

32.9.3 3. Intervention Operators

Structural modifications:

Action intervention: \(do(A_t = a)\)
Parameter intervention: \(do(\theta \leftarrow \theta^*)\)
Mechanism intervention: \(do(f \leftarrow f^*)\)

32.9.4 4. Exogenous Inputs and Variability

Explicit noise structure:

Process noise: \(U^x_t\) (intrinsic variability)
Observation noise: \(U^y_t\) (measurement error)
Action noise: \(U^a_t\) (behaviour policy randomness)

32.9.5 5. Graph Structure and Sparse Matrices (When Applicable)

For networked CDMs, the graph structure \(G\) and its sparse matrix representation provide the computational bridge between graph theory, linear algebra, and structural equations:

Graph \(G\): Encodes which variables directly influence which others
Adjacency matrix \(A\): Sparse matrix representation of graph structure
Transition matrix \(F\): For linear systems, sparse matrix encoding how state components influence one another
Sparse operations: Enable efficient computation for large systems when dependencies are sparse

The sparsity pattern of matrices reflects sparsity of direct dependencies. This three-way connection (graph theory → sparse matrices → linear algebra) enables scalable computation while maintaining causal semantics.

32.9.5.1 Implementation: Graph structure of a CDM

The graph structure of a CDM can be represented and analysed with ordinary directed graphs (here via CausalDynamics helpers and DAGMakie figures). Package APIs keep Pearl names; see Concept Reference Tables 4 and 8 for how book prose maps to functions:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CausalDynamics DAGMakie Graphs GraphMakie CairoMakie

# Example: Simple CDM with graph structure
# Nodes: 1=X_t, 2=A_t, 3=X_{t+1}, 4=Y_t, 5=C
g = DiGraph(5)
add_edge!(g, 1, 3)  # X_t → X_{t+1} (state transition)
add_edge!(g, 2, 3)  # A_t → X_{t+1} (intervention affects state)
add_edge!(g, 5, 3)  # C → X_{t+1} (context affects state)
add_edge!(g, 3, 4)  # X_{t+1} → Y_t (observation)

# Validate graph structure
println("Graph is valid DAG: ", CausalDynamics.is_dag(g))  # true

# Find Markov boundary of outcome Y_t
mb_Y = markov_boundary(g, 4)
println("Markov boundary of Y_t: ", mb_Y)  # Set([3]) = {X_{t+1}}

# Check if intervention A_t → Y_t is identifiable
adj_set = backdoor_adjustment_set(g, 2, 4)
println("Backdoor adjustment set for A_t → Y_t: ", adj_set)  # Set([1, 5]) = {X_t, C}

# Time-indexed unrolling (same lag semantics as Ch. 05b OCE → `TemporalDAGSpec`)
spec = TemporalDAGSpec(
    nodes = [
        TemporalNodeSpec(
            variable;
            temporal_support = PointwiseSupport(),
            value_representation = :state,
            referent = ReferentSpec(
                :unit;
                identity_criterion = :administrative_identifier,
            ),
        ) for variable in [:x, :y, :a, :c]
    ],
    edges = [
        (:c, :c, 1), (:a, :c, 1), (:c, :a, 0),
        (:x, :x, 1), (:a, :x, 1), (:c, :x, 1), (:x, :y, 0),
    ],
    graph_kind = TimeUnrolledGraph(),
)
u = unroll_temporal_dag(spec, 5)
# Effect of A_{t-1} on X_t: treatment at t=1, outcome at t=2
adj_lag = temporal_backdoor_adjustment_nodes(u, :a, 1, :x, 2)
println("Temporal backdoor nodes for A_{t-1} → X_t: ", adj_lag)  # Set([(:c, 1)])

# Visualise CDM graph structure
let
    fig = plot_causal_graph(g;
        node_labels = ["X_t", "A_t", "X_{t+1}", "Y_t", "C"],
        highlight_nodes = Set([2, 3, 4]),  # Highlight intervention path
        highlight_edges = [(2, 3), (3, 4)]  # Highlight intervention → state → outcome
    )
    fig  # Only this gets displayed
end

Graph is valid DAG: true
Markov boundary of Y_t: Set([3])
Backdoor adjustment set for A_t → Y_t: Set{Int64}()
Temporal backdoor nodes for A_{t-1} → X_t: Set([(:c, 1)])

Causal Dynamical Model graph structure showing state transitions and interventions

The graph structure enables us to reason about identifiability, determine necessary adjustments, and understand the causal relationships in the CDM.

32.10 Learning CDMs from Data

When the exact dynamics of a CDM are unknown, we can learn them from time series data using Universal Differential Equations (UDEs) (Rackauckas et al. 2020). Causal structure identification (graph / adjustment) combines with dynamics learning (SciML: Lux.jl, OrdinaryDiffEq.jl, Optimization.jl). For continuous mechanisms that already have an ODE right-hand side, the book’s Julia stack also exposes a thin SciML façade (ContinuousCDMSpec, solve_cdm, terminal_state) once OrdinaryDiffEq is loaded, see the short recipe below before the UDE training loop.

For graph-constrained neural residuals (parents from the DAG only), CausalDynamics also ships a Lux weakdep path (MechanismLibrary, build_ode_rhs, train_mechanisms!) documented in the SciML recipes. The hand-rolled Lux UDE loop below remains the pedagogical walkthrough; UniversalDiffEq.jl stays an optional advanced trainer outside CausalDynamics core. High-dim assays (spectra, images) should enter as codes via RepresentationSpec before identification or L3 abduction (:generative mechanisms), not as raw tensors on the DAG.

32.10.1 Continuous CDM via SciML helpers

project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))
@auto_using CausalDynamics Graphs OrdinaryDiffEq

# Identification on the summary graph (same confounding motif as above)
g_ode = DiGraph(3)
add_edge!(g_ode, 1, 2)
add_edge!(g_ode, 1, 3)
add_edge!(g_ode, 2, 3)
println("Backdoor set X → Y: ", backdoor_adjustment_set(g_ode, 2, 3))

# Continuous endogenous state via CausalDynamicsSciMLExt
spec = ContinuousCDMSpec([:prey, :predator])
function lotka_cdm!(du, u, p, t)
    X, Y = u
    du[1] = p.r * X - p.α * X * Y
    du[2] = p.β * X * Y - p.δ * Y
    return nothing
end
p_lv = (r = 1.0, α = 0.1, β = 0.02, δ = 0.5)
sol_lv = solve_cdm(spec, lotka_cdm!, [40.0, 9.0], (0.0, 10.0), p_lv)
term_lv = terminal_state(spec, sol_lv)
println("ODE terminal prey: ", round(term_lv.prey; digits = 2),
    "; predator: ", round(term_lv.predator; digits = 2))

Backdoor set X → Y: Set([1])
ODE terminal prey: 38.27; predator: 12.01

32.10.2 Codes, then graph-constrained generative L3

High-dimensional assays (spectra, images) are not DAG nodes. Compress them with RepresentationSpec, then attach a parent-constrained :generative mechanism on the codes. Abduction recovers additive U for a same-unit counterfactual; heavier ledgers live in the package Deep SCM stress notebooks.

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CausalDynamics DataFrames Random Statistics Lux ComponentArrays FiniteDiff
@auto_using Optimization OptimizationOptimisers

rng = Random.Xoshiro(11)
n, p = 120, 24
grid = range(0, 1; length = p)
bump = exp.(-((collect(grid) .- 0.4) ./ 0.08) .^ 2)
Mstar = randn(rng, n)
S = randn(rng, n, p) .* 0.12
for i in 1:n
    S[i, :] .+= Mstar[i] .* bump
end
# Matched-filter encode: spectrum → scalar code z1 (definitional role)
encode = Sm -> hcat(Sm * bump)
spec_rep = RepresentationSpec(:spectrum, [:z1], encode; role = :definitional)
panel = encode_to_panel(DataFrame(id = 1:n), S, spec_rep)
ρ = cor(panel.z1, Mstar)
println("code–latent cor(z1, M*) ≈ ", round(ρ; digits = 3))
@assert lux_mechanisms_available()

# Train Y = f(z1) + U on a synthetic outcome linear in z1
lib = MechanismLibrary(; allowed_parents = Dict(:z1 => Symbol[], :y => [:z1]))
attach_lux_mechanism!(
    lib, :y; parents = [:z1], hidden = 8, kind = :generative, rng = rng,
)
m_y = lib.mechanisms[:y]
z = panel.z1
y = 1.5 .* z .+ 0.2 .* randn(rng, n)
loss = function (_ps)
    s = 0.0
    for i in eachindex(z)
        s += abs2(generate_from_noise(m_y, 0.0, [z[i]]) - y[i])
    end
    return s / length(z)
end
train_mechanisms!(lib; loss = loss, maxiters = 200, lr = 0.05)
z_f, u_true = z[1], -0.35
y_f = generate_from_noise(m_y, u_true, [z_f])
y_cf = mechanism_counterfactual(m_y, y_f, [z_f], [0.0])
println("abduced U ≈ ", round(abduce_noise(m_y, y_f, [z_f]); digits = 3),
    " (true ", u_true, "); Y under do(z1=0) ≈ ", round(y_cf; digits = 3))

code–latent cor(z1, M*) ≈ 0.997
abduced U ≈ -0.35 (true -0.35); Y under do(z1=0) ≈ -0.327

32.10.3 Complete Workflow: Structure → Dynamics → Forecasting

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))
@auto_using CausalDynamics Graphs OrdinaryDiffEq
@auto_using Lux ComponentArrays Random ForwardDiff
@auto_using Optimization OptimizationOptimisers

# ── Step 1: Identify causal structure (CausalDynamics.jl) ─────────────
# Graph: Z → X → Y, Z → Y (confounding)
g_cdm = DiGraph(3)
add_edge!(g_cdm, 1, 2)  # Z → X
add_edge!(g_cdm, 1, 3)  # Z → Y
add_edge!(g_cdm, 2, 3)  # X → Y

adj_set_cdm = backdoor_adjustment_set(g_cdm, 2, 3)
println("Backdoor adjustment set for X → Y: ", adj_set_cdm)  # Set([1]) = {Z}

# ── Step 2: Generate synthetic time series ─────────────────────────────
function true_cdm_dynamics!(du, u, p, t)
    """Ground-truth dynamics respecting the causal graph Z → X → Y, Z → Y."""
    Z, X, Y = u
    du[1] = -0.05 * Z                          # Z decays slowly (exogenous driver)
    du[2] = 0.1 * X + 0.3 * Z                  # X depends on Z
    du[3] = -0.1 * Y + 0.5 * X + 0.2 * Z      # Y depends on X and Z
end

u₀_cdm = [2.0, 0.5, 0.0]
tspan_cdm = (0.0, 30.0)
dt_obs_cdm = 0.5
t_obs_cdm = 0.0:dt_obs_cdm:30.0

prob_cdm = ODEProblem(true_cdm_dynamics!, u₀_cdm, tspan_cdm)
sol_cdm = solve(prob_cdm, Tsit5(); saveat = t_obs_cdm)
data_cdm = Array(sol_cdm)

# Add observation noise
rng_cdm = Xoshiro(42)
data_cdm_noisy = data_cdm .+ 0.05 .* randn(rng_cdm, size(data_cdm))
println("Generated ", size(data_cdm, 2), " time-series observations")

# ── Step 3: Define a UDE respecting the causal graph ───────────────────
# Known: linear self-dynamics (diagonal terms)
# Unknown: interaction terms (Z→X, X→Y, Z→Y) learned via neural networks
# The causal graph constrains which inputs each NN receives:
#   NN_zx: takes Z, outputs Z→X effect
#   NN_xy: takes [X, Z], outputs combined X→Y and Z→Y effect
nn_zx = Chain(Dense(1, 8, tanh), Dense(8, 1))
nn_xy = Chain(Dense(2, 8, tanh), Dense(8, 1))
ps_zx, st_zx = Lux.setup(rng_cdm, nn_zx)
ps_xy, st_xy = Lux.setup(rng_cdm, nn_xy)
const _st_zx = st_zx
const _st_xy = st_xy

function ude_cdm!(du, u, p, t)
    Z, X, Y = u
    r_z = p.r_z;  r_x = p.r_x;  r_y = p.r_y
    # Self-dynamics (known structure, learned rates)
    du[1] = r_z * Z
    # Interaction terms (learned by NNs, respecting causal graph)
    z_to_x = first(first(nn_zx([Z], p.nn_zx, _st_zx)))
    xy_to_y = first(first(nn_xy([X, Z], p.nn_xy, _st_xy)))
    du[2] = r_x * X + z_to_x
    du[3] = r_y * Y + xy_to_y
end

p0_cdm = ComponentArray(
    nn_zx = ComponentArray(ps_zx),
    nn_xy = ComponentArray(ps_xy),
    r_z = -0.1, r_x = 0.05, r_y = -0.05
)

# ── Step 4: Train the UDE ──────────────────────────────────────────────
function predict_cdm(p)
    prob = ODEProblem(ude_cdm!, u₀_cdm, tspan_cdm, p)
    solve(prob, Tsit5(); saveat = t_obs_cdm, abstol = 1e-7, reltol = 1e-7)
end

function loss_cdm(p, _)
    pred = predict_cdm(p)
    # Return a large finite penalty (not Inf) when the solver fails, because
    # ForwardDiff Dual numbers cannot represent Inf correctly through Lux
    SciMLBase.successful_retcode(pred.retcode) ? sum(abs2, Array(pred) .- data_cdm_noisy) : eltype(p)(1e10)
end

# ForwardDiff is faster than Zygote for small systems (3 variables)
opt_f_cdm = OptimizationFunction(loss_cdm, Optimization.AutoForwardDiff())
opt_prob_cdm = OptimizationProblem(opt_f_cdm, p0_cdm)

println("Training CDM-UDE (500 iterations)...")
opt_sol_cdm = Optimization.solve(opt_prob_cdm, OptimizationOptimisers.Adam(0.01);
                                  maxiters = 500)
println("Final loss: ", round(opt_sol_cdm.objective; digits = 4))
println("Learned self-dynamics: r_z=", round(opt_sol_cdm.u.r_z; digits=3),
        "  r_x=", round(opt_sol_cdm.u.r_x; digits=3),
        "  r_y=", round(opt_sol_cdm.u.r_y; digits=3))

# Verify causal structure is still valid
is_adjustable = is_backdoor_adjustable(g_cdm, 2, 3)
println("X → Y identifiable via backdoor? ", is_adjustable)
println("Causal structure validated: Z → X → Y, Z → Y")

Backdoor adjustment set for X → Y: Set([1])
Generated 61 time-series observations
Training CDM-UDE (500 iterations)...
Final loss: 49.2949
Learned self-dynamics: r_z=-0.06  r_x=0.091  r_y=0.094
X → Y identifiable via backdoor? true
Causal structure validated: Z → X → Y, Z → Y

This workflow demonstrates how to build complete CDMs from data by combining:

Causal structure (identified by CausalDynamics.jl)
Dynamics learning (via Lux.jl + OrdinaryDiffEq.jl + Optimization.jl)
Forecasting and intervention (using learned model)

Treat the UDE fit as a criticisable generative step, not a black box. Before interpreting learned rates on field series, train on trajectories simulated from a CDM with known parameters (or a known interventional contrast) and check recovery; unstable optimisation, multimodal losses, or disagreement with an oracle are fitting problems to fix first (Gelman et al. 2026). The same iterate-and-expand habit as in Chapter 21 applies once real assays enter the loss.

32.11 Worked example: Confounded treatment and do(·)

This section walks through a compact CDM for a latent host immune marker under confounded treatment or vaccination, a biomarker \(X_t\), noisy assay \(Y_t\), treatment channel \(A_t\), and nutrition or infection-intensity confounder \(C_t\). The Introduction CDM teaser shows a minimal two-occasion unroll with adjustment highlighting and \(do(A_{t-1})\); below we restore the full endogenous/exogenous story, structural equations, simulation, and naive versus adjusted estimation. The builder lives in scripts/biological_showcases.jl as host_immunity_cdm (protein_treatment_cdm remains an alias).

Diagram notation:

\(f\): State transition or evolution function (e.g., \(f: \varphi \cdot X_{t-1}\) means the state transition uses autoregressive coefficient \(\varphi\))
\(\pi\): Policy or intervention assignment function (e.g., \(\pi: \gamma \cdot C_{t-1}\) means the treatment policy assigns treatment based on \(\gamma \cdot C_{t-1}\))
\(h\): Observation or measurement function (e.g., \(h: X_t\) means the observation is generated from \(X_t\))

Key components:

Endogenous variables (blue/orange/green/yellow): \(X_t\) (latent state), \(Y_t\) (observation), \(A_t\) (action/intervention), \(C_t\) (confounder), determined by structural equations
Exogenous variables (pink): \(U^x_t\) (process-side input), \(U^y_t\) (measurement-side input), and \(U^c_t\) (context-side input); their dependence and stochastic interpretation must be declared
Confounding: \(C_{t-1}\) creates a backdoor path \(A_{t-1} \leftarrow C_{t-1} \rightarrow X_t\), requiring adjustment for \(C_{t-1}\) to identify the causal effect of \(A_{t-1}\) on \(X_t\). Note that \(A_{t-1}\) also affects \(C_t\) (treatment can reduce infection intensity), but this does not create confounding for the effect of \(A_{t-1}\) on \(X_t\) since \(C_t\) does not directly affect \(X_t\) (only \(C_{t-1}\) affects \(X_t\))

Structural equations:

Confounder (endogenous): \(C_t \coloneqq \rho \cdot C_{t-1} + \alpha \cdot A_{t-1} + U^c_t\) where \(U^c_t \sim \text{Gamma}(\alpha_c, \theta_c)\) (strictly positive), e.g., nutrition or infection intensity that evolves over time and is affected by treatment
Intervention (endogenous, confounded): \(A_{t-1} \coloneqq \gamma \cdot C_{t-1} + U^a_{t-1}\) and \(A_t \coloneqq \gamma \cdot C_t + U^a_t\), treatment or vaccine uptake depends on host condition at the same time point (confounding)
State transition (endogenous): \(X_t \coloneqq \varphi \cdot X_{t-1} + B \cdot A_{t-1} + D \cdot C_{t-1} + U^x_t\) where \(U^x_t \sim \mathcal{N}(0, \sigma_w^2)\), state depends on previous state, previous treatment, and previous confounder (\(\varphi\) is the autoregressive coefficient, not the treatment \(A_t\))
Observation (endogenous): \(Y_t \coloneqq X_t + U^y_t\) where \(U^y_t \sim \mathcal{N}(0, \sigma_v^2)\)
Intervention via do-operator: \(do(A_{t-1} = a_{t-1})\) sets \(A_{t-1}\) to \(a_{t-1}\), breaking the confounder link

Together, the equations above specify endogenous versus exogenous variables, the backdoor path through \(C_{t-1}\), adjustment for identification, and \(do(A_{t-1} = a_{t-1})\) as a replacement of the treatment assignment mechanism. The simulation below uses typed discrete-time CDM machinery (DiscreteTimeCDM, DoSequence / Policy, simulate, counterfactual, g_computation). Figure 32.2 shows the corresponding trajectories: one model supports forecasting, interventional comparison, and unit-level alternative actualisations under a declared coupling.

From CausalDynamics.jl

The package’s unit-level counterfactual convention reuses the same exogenous series \(\mathbf{u}\) under a different intervention. This is a declared cross-world coupling, not an identity criterion. The package’s counterfactual method steps the CDM with factual.noise rather than drawing fresh noise.

# packages/CausalDynamics.jl/src/cdm/discrete_cdm.jl (excerpt)
function counterfactual(
    cdm::DiscreteTimeCDM,
    noise::Dict{Symbol, <:AbstractVector};
    intervention::AbstractIntervention,
    initial::Union{Nothing, NamedTuple} = nothing,
)
    T = length(first(values(noise)))
    state0 = initial === nothing ? cdm.initialise(Random.Xoshiro(0)) : initial
    state = _apply_do_to_state(state0, intervention, 1)

    series = _empty_series(state, T)
    _store!(series, state, 1)
    noise_keys = Tuple(keys(noise))
    noise_out = Dict{Symbol, Vector{<:Real}}(k => copy(v) for (k, v) in noise)

    for t in 2:T
        # Same U_t as the factual world; only the intervention differs
        noise_t = NamedTuple{noise_keys}(
            ntuple(i -> noise[noise_keys[i]][t], length(noise_keys)),
        )
        state = cdm.step(state, t, noise_t, intervention)
        _store!(series, state, t)
    end

    return CDMTrajectory(T, series, noise_out)
end

That shared-noise loop is Pearl’s L3 made operational for discrete-time CDMs. The simulation chunks below call this API; the excerpt is the definition.

Time-indexed identification on the same lag structure uses TemporalDAGSpec / unroll_temporal_dag / TemporalEffectQuery → identify. The subsection after the host-immunity simulation shows the typed hand-off from that certificate through simulate_panel to CausalTargeted sequential LMTP, with a sibling bite for event-time / survival LMTP (Ch. 23; (Díaz et al. 2023, 2024)).

project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))
include(joinpath(project_root, "scripts", "biological_showcases.jl"))

@auto_using CausalDynamics Distributions Statistics StableRNGs CairoMakie

# Causal Dynamical Model via CausalDynamics.jl (`DiscreteTimeCDM`):
# X_t — latent immune marker; Y_t — noisy assay
# A_t — treatment / vaccine; C_t — nutrition or infection-intensity confounder
# U^x_t, U^y_t, U^c_t, U^a_t — declared exogenous draws sampled by `sample_noise`

p = default_host_immunity_parameters()
φ, B, D, ρ, α, γ = p.φ, p.B, p.D, p.ρ, p.α, p.γ
σ_w, σ_v, σ_c, σ_a = p.σ_w, p.σ_v, p.σ_c, p.σ_a
T = DEFAULT_CDM_T

cdm = host_immunity_cdm(p.x₀, p.c₀, φ, B, D, ρ, α, γ, σ_w, σ_v, σ_c, σ_a)

# Observational world: natural assignment A_t := γ·C_t + U^a_t
traj_conf = simulate(cdm, T; rng = StableRNG(34))
x_conf, y_conf, a_conf, c_conf = traj_conf.series[:x], traj_conf.series[:y], traj_conf.series[:a], traj_conf.series[:c]
u_x, u_y, u_c, u_a = traj_conf.noise[:u_x], traj_conf.noise[:u_y], traj_conf.noise[:u_c], traj_conf.noise[:u_a]

# Interventional worlds: do(A_t = 1) and do(A_t = 0)
traj_interv = simulate(cdm, T; rng = StableRNG(35), intervention = do_sequence(:a, ones(T)))
x_interv, y_interv, a_interv, c_interv = traj_interv.series[:x], traj_interv.series[:y], traj_interv.series[:a], traj_interv.series[:c]

traj_control = simulate(cdm, T; rng = StableRNG(36), intervention = do_sequence(:a, zeros(T)))
x_control, y_control, a_ctrl, c_ctrl = traj_control.series[:x], traj_control.series[:y], traj_control.series[:a], traj_control.series[:c]

# Shared-U counterfactual: reuse declared exogenous draws under do(A_t = 1)
traj_cf = counterfactual(
    cdm,
    traj_conf.noise;
    intervention = do_sequence(:a, ones(T)),
    initial = (
        x = traj_conf.series[:x][1],
        y = traj_conf.series[:y][1],
        a = traj_conf.series[:a][1],
        c = traj_conf.series[:c][1],
    ),
)

# Naive estimation (without adjusting for confounder) - BIASED
# This estimates association, not causation, because of confounding
a_high_mask = a_conf .> 0.5
a_low_mask = a_conf .<= 0.5
if sum(a_high_mask) > 0 && sum(a_low_mask) > 0
    naive_effect = mean(y_conf[a_high_mask]) - mean(y_conf[a_low_mask])
else
    naive_effect = NaN
end

# Adjusted estimation (conditioning on confounder) - UNBIASED
# Adjust by stratifying on confounder levels
c_high = c_conf .> median(c_conf)
c_low = c_conf .<= median(c_conf)
c_high_a_high = c_high .& a_high_mask
c_high_a_low = c_high .& a_low_mask
c_low_a_high = c_low .& a_high_mask
c_low_a_low = c_low .& a_low_mask

if sum(c_high_a_high) > 0 && sum(c_high_a_low) > 0
    adjusted_effect_high = mean(y_conf[c_high_a_high]) - mean(y_conf[c_high_a_low])
else
    adjusted_effect_high = NaN
end

if sum(c_low_a_high) > 0 && sum(c_low_a_low) > 0
    adjusted_effect_low = mean(y_conf[c_low_a_high]) - mean(y_conf[c_low_a_low])
else
    adjusted_effect_low = NaN
end

if !isnan(adjusted_effect_high) && !isnan(adjusted_effect_low)
    adjusted_effect = (adjusted_effect_high + adjusted_effect_low) / 2
else
    adjusted_effect = NaN
end

# True causal effect (from intervention simulations)
# Compare intervention group (do(A_t = 1)) to control group (do(A_t = 0))
true_effect = mean(y_interv) - mean(y_control)

# Typed interventions: constant assignment, soft policy, and Monte Carlo g-computation
do_constant = do_sequence(:a, 1.0)  # ConstantAssignment for every t
π_treat_if_low = policy(:a, (state, t) -> state.c <= median(c_conf) ? 1.0 : 0.0)
traj_policy = simulate(cdm, T; rng = StableRNG(37), intervention = π_treat_if_low)
gcomp_treat = g_computation(cdm, T, :y; intervention = do_sequence(:a, ones(T)), n = 80, rng = StableRNG(38))
gcomp_ctrl = g_computation(cdm, T, :y; intervention = do_sequence(:a, zeros(T)), n = 80, rng = StableRNG(39))
gcomp_effect = gcomp_treat.mean - gcomp_ctrl.mean

# Display results
println("CDM simulation via CausalDynamics.DiscreteTimeCDM complete!")
println("\nEndogenous variables (confounded simulation):")
println("  Final latent immune marker X[$(T)]: ", round(x_conf[end], digits=3))
println("  Final assay Y[$(T)]: ", round(y_conf[end], digits=3))
println("  Final confounder C[$(T)]: ", round(c_conf[end], digits=3))
println("  Final treatment A[$(T)]: ", round(a_conf[end], digits=3))
println("\nCausal effect estimation:")
println("  True causal effect (B parameter): ", round(B, digits=3))
println("  Naive estimate (biased by confounding): ", round(naive_effect, digits=3))
println("  Adjusted estimate (conditioning on C_t): ", round(adjusted_effect, digits=3))
println("  True effect from intervention: ", round(true_effect, digits=3))
println("  g_computation TE (n=80): ", round(gcomp_effect, digits=3),
    " (SE≈", round(sqrt(gcomp_treat.std^2 + gcomp_ctrl.std^2) / sqrt(80); digits = 3), ")")
println("  Soft policy mean Y: ", round(mean(traj_policy.series[:y]), digits=3),
    " (constant do(A=1) uses ", typeof(do_constant.values[:a]), ")")
println("  Shared-U counterfactual mean Y under do(A=1): ", round(mean(traj_cf.series[:y]), digits=3))
println("\nNote: The naive estimate is biased because C_t affects both A_t and X_t.")
println("Adjusting for C_t (stratification) recovers the true causal effect.")

CDM simulation via CausalDynamics.DiscreteTimeCDM complete!

Endogenous variables (confounded simulation):
  Final latent immune marker X[100]: 2.173
  Final assay Y[100]: 1.65
  Final confounder C[100]: 0.689
  Final treatment A[100]: 0.351

Causal effect estimation:
  True causal effect (B parameter): 0.5
  Naive estimate (biased by confounding): 0.229
  Adjusted estimate (conditioning on C_t): 0.39
  True effect from intervention: 4.402
  g_computation TE (n=80): 3.617 (SE≈0.153)
  Soft policy mean Y: 3.132 (constant do(A=1) uses ConstantAssignment{Float64})
  Shared-U counterfactual mean Y under do(A=1): 4.69

Note: The naive estimate is biased because C_t affects both A_t and X_t.
Adjusting for C_t (stratification) recovers the true causal effect.

Figure 32.2: Host immune marker over time: confounded vs interventional simulations showing the effect of adjusting for confounders

32.11.1 Fitting a host-immunity effect from noisy observations

The preceding trajectories use specified mechanisms. We now estimate a causal contrast from observations in a linear Gaussian companion of the same host-immunity model. The immune marker and confounder are centred scores; this example omits the non-negative truncation used above so that its causal effect is available analytically. The question is how a sustained one-unit change in baseline treatment changes the mean terminal assay after six transitions.

For subject \(i\), let \(C_{i,0}\) be a recorded baseline confounder, \(X_{i,t}\) the latent immune score, \(Y_{i,t}\) its noisy assay, and \(A_i\) a baseline treatment held constant throughout follow-up. The autoregressive weights are \(\phi\) and \(\rho\); \(B\) and \(D\) are treatment and confounder effects on the immune score, \(\alpha\) is the treatment effect on the evolving confounder, and \(\gamma\) controls baseline confounding. Independent, zero-mean Gaussian exogenous draws are denoted by \(U_i^A,U_{i,t}^X,U_{i,t}^C,U_{i,t}^Y\):

\[ \begin{aligned} A_i &= \gamma C_{i,0}+U_i^A,\\ X_{i,t+1} &= \phi X_{i,t}+B A_i+D C_{i,t}+U_{i,t+1}^X,\\ C_{i,t+1} &= \rho C_{i,t}+\alpha A_i+U_{i,t+1}^C,\\ Y_{i,t} &= X_{i,t}+U_{i,t}^Y. \end{aligned} \]

Baseline \(X_{i,0}\) and \(C_{i,0}\) are independent standard normal draws. Only \((C_{i,0},Y_{i,0},A_i,Y_{i,T})\) enter the estimator. The target is \(\theta_T=E[Y_T^{do(A=1)}-Y_T^{do(A=0)}]\) with \(T=6\). Because assignment depends only on recorded \(C_0\), adjusting for \(C_0\) blocks confounding in this model; \(Y_0\) adds baseline prognostic information. The Gaussian assignment gives positive density at both treatment values. Consistency, independent subjects, invariant mechanisms under treatment, and the stated absence of additional confounders are assumptions of this simulated analysis (Pearl 2009).

Marginalising the latent states gives a linear terminal-outcome regression on \(A,C_0,Y_0\). Its coefficient on \(A\) estimates the total sustained-treatment effect, including the path through subsequent \(C_t\). Adjusting for those post-treatment confounder values would change the target. This simple regression is justified by the baseline assignment and linear model here; sequentially assigned treatment requires the longitudinal methods developed earlier.

For an independent reference calculation, let \(d_t\) and \(e_t\) be the differences in latent immune and confounder states between \(do(A=1)\) and \(do(A=0)\) with common initial states. Then \(d_0=e_0=0\), \(d_{t+1}=\phi d_t+B+D e_t\), \(e_{t+1}=\rho e_t+\alpha\), and \(\theta_T=d_T\). The chunk-host-immunity-fitted chunk compares this known effect with fitted estimates.

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CausalDynamics DataFrames Distributions StableRNGs Statistics LinearAlgebra Random
include(joinpath(dirname(Base.active_project()), "scripts", "biological_showcases.jl"))
host_fit = linear_host_analysis()
println("| Quantity | Value |")
println("|---|---:|")
for (name, value) in [
    ("Known sustained-treatment effect", host_fit.truth),
    ("Adjusted estimate", host_fit.estimate),
    ("Subject bootstrap 95% interval, lower", host_fit.interval[1]),
    ("Subject bootstrap 95% interval, upper", host_fit.interval[2]),
    ("Estimate omitting baseline confounder", host_fit.naive),
    ("Held-out assay RMSE, adjusted", host_fit.rmse),
    ("Held-out assay RMSE, confounder omitted", host_fit.naive_rmse),
]
    println("| ", name, " | ", round(value; digits=3), " |")
end
# Oracle simulator check, separate from estimation from observed data.
@assert all(isapprox.(host_fit.oracle_contrasts, host_fit.truth; atol=1e-12))

| Quantity | Value |
|---|---:|
| Known sustained-treatment effect | 1.192 |
| Adjusted estimate | 1.185 |
| Subject bootstrap 95% interval, lower | 1.138 |
| Subject bootstrap 95% interval, upper | 1.236 |
| Estimate omitting baseline confounder | 1.461 |
| Held-out assay RMSE, adjusted | 0.651 |
| Held-out assay RMSE, confounder omitted | 0.686 |

The interval resamples 2,000 independent subjects 400 times and refits the model. It describes sampling uncertainty under this model, excluding uncertainty about the graph or observation mechanism. A separately generated cohort of 1,000 subjects supplies the prediction check. The regression omitting \(C_0\) deliberately violates the adjustment requirement; its treatment coefficient can be misleading even when its predictive error appears acceptable. Figure 32.3 places both estimates against the known effect.

Figure 32.3: Sustained treatment in the linear host-immunity companion: adjusted estimate with a subject-bootstrap 95% interval, the estimate omitting the baseline confounder, and the known effect (dashed line).

Shared exogenous draws also recover the analytic contrast exactly in the simulator. Those draws are available because we generated them; the recorded assays do not reveal them. A patient-specific counterfactual would require abduction under an observation likelihood and integration over the remaining latent uncertainty. The constant individual contrast here follows from additive linear mechanisms and is not evidence that real treatment responses are homogeneous. The example establishes a reproducible population analysis and a separate oracle check of counterfactual semantics, within explicitly stated assumptions.

32.11.2 Certificate → panel → sequential LMTP

Simulation estimands (g_computation above) answer interventional means when the CDM is fully known. When only observational trajectories are available, the same lag structure yields a temporal identification certificate, a wide panel (simulate_panel), and a CausalTargeted SequentialPolicy estimate. Contrast the Monte Carlo TE printed above with the sequential LMTP estimate below: same scientific question, different information and estimator.

Nested sampling units (farms within regions, animals within farms) belong in the generative layer: CausalDynamics RandomEffectSpec / simulate_hierarchical_panel draws shared exogenous \(U\) by cluster, and unroll_hierarchical_dag expands the identification graph. CausalTargeted may then use cluster= on repeated-outcome MSM so the IF sandwich matches that dependence. Fitting a mixed model for \(\widehat{\mathrm{ATE}}\) stays in application code (see the Sheep Vaccine hierarchy note and CausalDynamics examples/hierarchical_mixedmodels_concordance.md); it is not part of the LMTP or mediation APIs.

@auto_using CausalDynamics CausalTargeted DataFrames StableRNGs Statistics LinearAlgebra

cols, truth = simulate_hierarchical_intercept_ate(
    600; n_clusters = 30, σ_cluster = 1.0, β_a = 0.5, σ_y = 0.35,
    rng = StableRNG(51),
)
X = [ones(length(cols[:Y])) cols[:A] cols[:W]]
ate_ols = (X \ cols[:Y])[2]

df, τ_true = simulate_repeated_outcome_ate(
    400; T = 3, β_a = [0.2, 0.45, 0.35], rng = StableRNG(52),
)
df.cluster = Float64.(mod1.(1:nrow(df), 20))
rng_u = StableRNG(53)
u_by = Dict{Float64, Float64}(c => randn(rng_u) for c in unique(df.cluster))
uv = [u_by[c] for c in df.cluster]
for y in (:Y1, :Y2, :Y3)
    df[!, y] = df[!, y] .+ uv
end
msm_u = run_repeated_outcome_msm(
    df, :A, [:Y1, :Y2, :Y3];
    baseline = [:W], folds = 3, learners = (:glm, :mean), rng = StableRNG(54),
)
msm_c = run_repeated_outcome_msm(
    df, :A, [:Y1, :Y2, :Y3];
    baseline = [:W], folds = 3, learners = (:glm, :mean), rng = StableRNG(54),
    cluster = :cluster,
)
(
    ate_ols = ate_ols,
    truth_ate = truth.ate,
    msm_estimates = msm_c.estimates,
    τ_true = τ_true.tau,
    covariance_kind = msm_c.covariance_kind,
    same_point = msm_c.estimates ≈ msm_u.estimates,
)

(ate_ols = 0.5261982941891977, truth_ate = 0.5, msm_estimates = [0.21610053767182208, 0.44100985270150594, 0.3119385762353112], τ_true = [0.2, 0.45, 0.35], covariance_kind = :cluster, same_point = true)

@auto_using CausalDynamics CausalTargeted DataFrames StableRNGs

# Compact two-occasion CDM (W baseline, A/L timed, Y terminal)
cdm_panel = DiscreteTimeCDM(
    [:w, :a, :l, :y];
    initialise = (rng) -> (w = randn(rng), a = 0.0, l = 0.0, y = 0.0),
    sample_noise = (rng, state, t) -> (
        u_a = randn(rng), u_l = randn(rng), u_y = randn(rng),
    ),
    step = (state, t, noise, intervention) -> begin
        a = intervention_value(intervention, :a, t, 0.5 * state.w + noise.u_a)
        l = 0.3 * a + noise.u_l
        y = 0.4 * a + 0.2 * l + 0.1 * state.w + noise.u_y
        (w = state.w, a = a, l = l, y = y)
    end,
)
spec_lag = TemporalDAGSpec(
    nodes = [
        TemporalNodeSpec(
            variable;
            temporal_support = PointwiseSupport(),
            value_representation = :state,
            referent = ReferentSpec(
                :unit;
                identity_criterion = :administrative_identifier,
            ),
        ) for variable in [:a, :l, :y, :w]
    ],
    edges = [
        LaggedEdge(:w, :a, 0),
        LaggedEdge(:a, :l, 0),
        LaggedEdge(:l, :a, 1),
        LaggedEdge(:a, :y, 0),
        LaggedEdge(:w, :y, 0),
    ],
    graph_kind = TimeUnrolledGraph(),
)
u_lag = unroll_temporal_dag(spec_lag, 2)
id_seq = identify(u_lag, TemporalEffectQuery(:a, :y, 1, 2))
panel = simulate_panel(
    cdm_panel, 120, 2;
    rng = StableRNG(40),
    baseline = [:w],
    timed = [:a, :l],
    terminal = [:y],
    temporal_spec = spec_lag,
)
df_panel = DataFrame(NamedTuple(panel))
est_seq = sequential_spec_from_identification(
    id_seq;
    baseline = [:w],
    time_vary = [Symbol[], [:l1]],
    shift = additive_shift_policy(; lower_q = 0.01, upper_q = 0.99),
)
seq_grid = execute_estimand(
    est_seq, df_panel;
    id_result = id_seq,
    deltas = [0.5],
    folds = 2,
    learners = SMALL_N_SL_LEARNERS,
    metadata = false,
    rng = StableRNG(41),
)
gcomp_panel = g_computation(
    cdm_panel, 2, :y;
    intervention = do_sequence(:a, [0.5, 0.5]),
    n = 80,
    rng = StableRNG(42),
)
gcomp_nat = g_computation(
    cdm_panel, 2, :y;
    intervention = do_sequence(:a, [0.0, 0.0]),
    n = 80,
    rng = StableRNG(43),
)
println("Temporal ID: strategy=", id_seq.strategy, " adjustment=", id_seq.adjustment)
println("Panel columns: ", names(df_panel))
println("Sequential LMTP TE (δ=0.5): ", round(only(seq_grid.est); digits = 3),
    " (SE ", round(only(seq_grid.se); digits = 3), ")")
println("g_computation contrast (do(A=0.5) − do(A=0), n=80): ",
    round(gcomp_panel.mean - gcomp_nat.mean; digits = 3))

Temporal ID: strategy=temporal_backdoor adjustment=[:l, :w]
Panel columns: ["w", "a1", "l1", "a2", "l2", "y"]
Sequential LMTP TE (δ=0.5): 0.417 (SE 0.076)
g_computation contrast (do(A=0.5) − do(A=0), n=80): 0.249

Point-treatment continuous MTP grids (non-sequential) remain in Policy Evaluation; mediation certificates use CausalMediation (Ch. 20). Factor recodes of \(A_t\) use the same sequential path with policies = discrete_recode_policy(...). CDMPanel stays Float64 (integer codes such as \(0,1,2\)); convert to Int at the CausalTargeted boundary so a factor recode is not treated as a continuous shift. Mixed continuous/discrete sequential \(A_t\) is rejected.

32.11.3 Certificate → event-time / survival LMTP

The same temporal ID spirit applies when the terminal scientific target is remaining event-free through a horizon rather than a continuous \(Y\). Below, a synthetic discrete-time survival DGP supplies treatments \(A_t\) and monotone event-free indicators \(S_t\); CausalTargeted estimates \(E[S_h^{d}]\) under an MTP shift (Díaz et al. 2024). Competing risks are deferred.

@auto_using CausalDynamics CausalTargeted DataFrames StableRNGs

df_ev, truth_ev = simulate_discrete_survival_mtp(
    140; T = 3, β_a = -0.5, rng = StableRNG(44),
)
spec_ev = TemporalDAGSpec(
    nodes = [
        TemporalNodeSpec(
            variable;
            temporal_support = PointwiseSupport(),
            value_representation = :state,
            referent = ReferentSpec(
                :individual;
                identity_criterion = :administrative_identifier,
            ),
        ) for variable in [:a, :s, :w]
    ],
    edges = [
        LaggedEdge(:w, :a, 0),
        LaggedEdge(:a, :s, 0),
        LaggedEdge(:w, :s, 0),
        LaggedEdge(:s, :s, 1),
    ],
    graph_kind = TimeUnrolledGraph(),
)
u_ev = unroll_temporal_dag(spec_ev, 3)
id_ev = identify(u_ev, TemporalEffectQuery(:a, :s, 1, 3))
cert_ev = survival_identification_certificate(u_ev, TemporalEffectQuery(:a, :s, 1, 3))
shift_ev = additive_shift_policy(; scale = "raw", lower_q = 0.0, upper_q = 1.0)
est_ev = SurvivalPolicy(
    truth_ev.treatments, truth_ev.surv, [:W];
    shift = shift_ev, horizon = 3,
)
ev_grid = execute_estimand(
    est_ev, df_ev;
    id_result = cert_ev.result,
    deltas = [0.75],
    folds = 2,
    learners = SMALL_N_SL_LEARNERS,
    metadata = false,
    rng = StableRNG(45),
)
ψ_ev = truth_ev.survival(0.75)
println("Temporal survival ID: strategy=", id_ev.strategy,
    " adjustment=", id_ev.adjustment)
println("Survival LMTP E[S₃] (δ=0.75): ", round(only(ev_grid.est); digits = 3),
    " (SE ", round(only(ev_grid.se); digits = 3), "; oracle ",
    round(ψ_ev; digits = 3), ")")

Temporal survival ID: strategy=temporal_backdoor adjustment=[:w]
Survival LMTP E[S₃] (δ=0.75): 0.489 (SE 0.02; oracle 0.545)

Use this path for longitudinal studies whose primary estimand is event-free probability under a feasible exposure shift (not a continuous AUC or terminal load). Point-treatment pedagogy for the same API is in Policy Evaluation.

32.11.4 Latent series → observed panel

When the assay is a noisy function of a latent state (or a filter supplies \(\hat{x}_t\)), build the same wide panel without putting Kalman code in CausalDynamics: supply inferred series and an ObservationBridge.

@auto_using CausalDynamics StableRNGs

# Pretend a filter returned latent x and confounder w for two units
inferred = [
    Dict(:w => [0.2, 0.2], :x => [1.0, 1.2], :a => [0.0, 1.0]),
    Dict(:w => [-0.1, -0.1], :x => [0.5, 0.8], :a => [1.0, 0.0]),
]
# Rename latent :x → observed :y; pass through :w and :a
bridge = ObservationBridge(Dict(:w => :w, :a => :a, :x => :y))
panel_lat = panel_from_latent_series(
    inferred;
    bridge = bridge,
    baseline = [:w],
    timed = [:a, :y],
)
println("Latent→observed columns: ", panel_lat.column_order)
println("Terminal-style y2: ", round.(panel_lat.data[:y2]; digits = 3))

Latent→observed columns: [:w, :a1, :y1, :a2, :y2]
Terminal-style y2: [1.2, 0.8]

32.11.5 Continuous-time functional (SciML)

Continuous CDMs use the same Monte Carlo idea with ContinuousEffectFunctional (:terminal, :mean, or :integral) once OrdinaryDiffEq is loaded.

@auto_using CausalDynamics OrdinaryDiffEq StableRNGs

spec_c = ContinuousCDMSpec([:x, :y])
function drift_c!(du, u, p, t)
    du[1] = 0.0
    du[2] = u[1]
    return nothing
end
g_c = g_computation(
    spec_c, drift_c!, rng -> [1.0, 0.0], (0.0, 2.0), ();
    intervention = do_pin(:x, 1.0),
    functional = ContinuousEffectFunctional(:y; kind = :terminal),
    n = 30,
    rng = StableRNG(44),
)
println("Continuous g_computation E[Y(T) | do(X=1)]: ", round(g_c.mean; digits = 3))

Continuous g_computation E[Y(T) | do(X=1)]: 2.0

32.12 Example: Ecological CDM

32.12.1 Process Model

Lotka-Volterra dynamics: \[ \begin{aligned} S_{t+1} &= S_t + \Delta t \cdot [r S_t (1 - S_t/K) - \alpha S_t P_t] + U^s_{t+1} \\ P_{t+1} &= P_t + \Delta t \cdot [\beta \alpha S_t P_t - \delta P_t] + U^p_{t+1} \end{aligned} \]

32.12.2 Observation Model

Count observations with detection probability: \[ Y_t \sim \text{Bernoulli}(p) \cdot \text{Poisson}(\lambda(S_t)) \]

32.12.3 Interventions

Remove predators: \(do(P_t = 0)\)
Harvest prey: \(do(S_t \leftarrow S_t - H_t)\)
Change parameters: \(do(r \leftarrow r^*)\)

32.12.4 Counterfactuals

For fixed \(\mathbf{u}\), simulate: “What would prey population have been if predators were removed?”

32.12.5 Reference stack: AgeSCM four-layer pipeline

Case Study 3 (Ch 28b) instantiates the unified CDM workflow for observable data at scale:

Structural (CausalDynamics + DAG) → Encoder (Flux, ẑ) → PPL (GraphPPL + RxInfer) → Transport (IPTW / propensity)

The same CDM supports E1 (predictive transport), E2 (mechanism stability of ẑ), and validation (Ch 21). High-dimensional Y enters only through ẑ; domain shift on C is handled in the training objective, not by ignoring the causal graph.

32.13 Key Takeaways

CDMs unify SCMs, SSMs, and dynamical systems in a single framework
One CDM supports forecasting, interventions, and counterfactuals
Explicit process, observation, and coupling structure enables counterfactual reasoning
CDMs provide a common language for causal-dynamical modelling
Alternative representations (graph-based, diagrammatic) can be combined to build complete CDMs
Composed systems (Decapodes.jl) can exhibit attractor states relative to their declared dynamics, and composition can affect the attractor geometry
Network structure (graphs/hypergraphs) contributes to attractor dynamics, basin boundaries, and system resilience through the specified mechanisms
AgeSCM demonstrates the four-layer observable pipeline (structure → encoder → PPL → transport weights) on MIRS mosquito age

32.14 Alternative Representations: Diagrammatic CDMs

While this book primarily uses graph-based representations for CDM structure (graph \(G\) encoding directed dependencies), there are alternative formalisms that provide complementary perspectives. Decapodes.jl (Morris et al. 2024) offers a diagrammatic representation based on Applied Category Theory (ACT) for composing and simulating physical systems.

32.14.1 Diagrammatic vs Graph-Based Representation

Graph-based (this book):

Graph \(G\) encodes directed dependencies (which variables directly influence which others)
Structural equations: \(X_{t+1} = f(X_t, A_t, \ldots)\)
Natural for causal reasoning (d-separation, do-calculus)
Direct connection to SCM semantics

Diagrammatic (Decapodes.jl):

ACT diagrams represent equations as formal diagrams
Hierarchical composition of systems
Natural for composing multiphysics systems
Declarative DSL for physics equations

Connection: Both represent system structure, but:

Graphs emphasize causal relations (who influences whom)
Diagrams emphasize compositional structure (how systems combine)

32.14.2 When to Use Each

Use a graph-based causal layer when:

Focus is on causal reasoning (interventions, counterfactuals)
Need identification algorithms (backdoor, frontdoor, do-calculus)
Working with observational data and need to reason about confounding
System structure is causal (who influences whom)

Use a diagrammatic / compositional physics layer (e.g. Decapodes.jl) when:

Focus is on composing physical systems (multiphysics)
Need hierarchical composition of subsystems
Working with deterministic physics equations (ODEs/PDEs)
Want declarative representation of equations

Integration: For CDMs, we can combine both:

Composition: mechanistic process model (ODEs/SDEs)
Causal semantics: interventions and counterfactuals
Inference: observation model and state estimation

This is a conceptual workflow: compose → attach \(do(\cdot)\) semantics → infer. The Julia packages named in code chunks are one way to enact it.

32.14.3 Example: Composing Ecological CDM with Decapodes.jl

The following example demonstrates how to use Decapodes.jl to compose a simple reaction-diffusion model (based on the Brusselator example from the Decapodes.jl repository). This shows the declarative DSL for representing physics equations, which can then be integrated into a complete CDM with causal semantics and observation models.

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using Decapodes Catlab CombinatorialSpaces DiagrammaticEquations OrdinaryDiffEq ComponentArrays CairoMakie DAGMakie Graphs

# Define a simple reaction-diffusion model using Decapodes.jl's declarative DSL
# This represents the process model component of a CDM
# Based on the Brusselator example from Decapodes.jl repository
SimpleReactionDiffusion = @decapode begin
  (U, V)::Form0  # State variables (e.g., chemical concentrations)
  U2V::Form0     # Intermediate variable: U^2 * V
  (U̇, V̇)::Form0  # Time derivatives

  (α)::Constant   # Diffusion coefficient
  F::Parameter    # Forcing term (could represent intervention)

  # Reaction term: U^2 * V
  U2V == (U .* U) .* V

  # Reaction-diffusion equations
  U̇ == 1 + U2V - (4.4 * U) + (α * Δ(U)) + F
  V̇ == (3.4 * U) - U2V + (α * Δ(V))

  # Time derivatives
  ∂ₜ(U) == U̇
  ∂ₜ(V) == V̇
end

# Visualise the diagrammatic representation
# Convert the Decapode to Graphs.jl format for visualisation with DAGMakie
# Note: Decapodes are represented as computational graphs, which we can visualize
try
    # Extract graph structure from decapode for visualization
    # The decapode represents a computational graph of the equations
    # We create a simplified visualization showing key dependencies

    # Create a graph representing the equation structure
    # Nodes: 1=U, 2=V, 3=U2V, 4=U̇, 5=V̇
    # Qualify Graphs APIs: Catlab/Decapodes can leave bare `add_edge!` unbound in the notebook
    g_viz = Graphs.SimpleDiGraph(5)
    Graphs.add_edge!(g_viz, 1, 3)  # U -> U2V
    Graphs.add_edge!(g_viz, 2, 3)  # V -> U2V
    Graphs.add_edge!(g_viz, 1, 4)  # U -> U̇ (direct)
    Graphs.add_edge!(g_viz, 3, 4)  # U2V -> U̇
    # Note: U -> U̇ via diffusion is already represented above
    Graphs.add_edge!(g_viz, 1, 5)  # U -> V̇
    Graphs.add_edge!(g_viz, 3, 5)  # U2V -> V̇
    Graphs.add_edge!(g_viz, 2, 5)  # V -> V̇ (via diffusion)

    # Visualise the graph using DAGMakie
    let
        node_labels = ["U", "V", "U²V", "U̇", "V̇"]
        fig, ax, p = dagplot(g_viz;
            figure_size = (800, 600),
            layout_mode = :acyclic,
            node_color = RGBf(0.30, 0.55, 0.70),
            labels = node_labels,
            node_size = 38,
            edge_width = 2.0,
            arrow_size = 15,
            nlabels_fontsize = 16,
            nlabels_color = :white,
            nlabels_align = (:center, :center),
            nlabels_distance = 0,
            padding = 0.40,
        )
        fig  # Only this gets displayed
    end
catch e
    println("Visualization not available: ", e)
    println("The decapode has been defined successfully and can be used for simulation.")
end

# Example: Analyze attractor dynamics of composed system
# After composing subsystems with Decapodes.jl, we can:
# 1. Find equilibrium points (attractors) using NLsolve.jl or steady-state solvers
# 2. Analyze stability (eigenvalues, Lyapunov exponents)
# 3. Study basin of attraction (which initial conditions lead to which attractors)
# 4. Perform interventions: $do(F = f^*)$ to shift attractor landscape
# 5. Analyze resilience: How does composition affect recovery to equilibrium?

# Connection to attractor dynamics:
# - Composed systems evolve towards attractor states
# - Composition affects attractor landscape (number, stability, basin size)
# - Interventions modify attractor structure
# - Resilience depends on attractor stability and basin size

Composing CDM process model using Decapodes.jl diagrammatic representation

The diagrammatic representation above shows how Decapodes.jl encodes the reaction-diffusion equations as a formal diagram. This provides a declarative representation of the process model that can be composed with other systems. The F parameter could represent an intervention (e.g., \(do(F = f^*)\)), showing how interventions can be incorporated into the diagrammatic representation.

Connection to attractor dynamics: After composing subsystems with Decapodes.jl, we can analyse how the composed system evolves towards attractor states. The composition determines the attractor landscape (number of attractors, their stability, basin boundaries), and interventions modify this landscape. This connects Decapodes.jl to the book’s core theme of systems evolving towards attractors.

To complete the CDM, we would then:

Add causal semantics (using CausalDynamics.jl): Define interventions like \(do(P = 0)\) (remove predators) or \(do(r \leftarrow r^*)\) (change growth rate)
Add observation model (using StateSpaceDynamics.jl): Define how we observe the system, e.g., \(Y_t \sim \text{Poisson}(\lambda(S_t))\) for count observations
Perform inference: Use StateSpaceDynamics.jl to infer latent states from noisy observations

This demonstrates how different representations (diagrammatic vs graph-based) can be combined to build complete CDMs: a compositional process model, causal semantics for interventions and counterfactuals, and an inference layer for latent state. The packages used below implement that pattern.

Key insight: Compositional physics frameworks build process models; a causal graph layer adds intervention and counterfactual semantics. Together they support complex CDMs from composable components.

32.14.4 Composed Systems and Attractor Dynamics

When using Decapodes.jl to compose multiple subsystems, the resulting system can exhibit complex attractor dynamics that depend on how subsystems are composed:

Composition affects attractor landscape:

Subsystem attractors: Each composed subsystem may have its own attractor states
Coupled attractors: Composed systems can have attractors that emerge from subsystem interactions
Attractor transitions: Composition can create or destroy attractors, change stability, or shift basin boundaries

Example: Composing predator-prey dynamics with resource dynamics:

Predator-prey subsystem: May have limit cycle attractor (oscillations)
Resource subsystem: May have stable equilibrium attractor
Composed system: May exhibit new attractor states (e.g., chaotic attractor) or modified stability

Interventions modify attractor structure:

Parameter interventions: \(do(\theta \leftarrow \theta^*)\) can shift attractors or change stability
Forcing interventions: \(do(F = f^*)\) can push system to different attractor basin
Composition interventions: Modifying how subsystems are composed can restructure attractor landscape

Connection to resilience and robustness:

Basin of attraction: Composition affects which initial conditions lead to which attractors
Recovery dynamics: Composed systems may recover to different attractors after perturbations
Structural robustness: How composition affects system’s ability to maintain function under variation

Optional FEP interpretation: If a composed system is given an explicit generative model, preferences and variational objective, its policy or dynamics may be analysed using FEP language. Composition alone does not imply free-energy minimisation; its causal and attractor consequences follow from the declared mechanisms.

This connects Decapodes.jl to the book’s core theme of systems evolving towards attractors: by composing subsystems, we can understand how complex attractor landscapes emerge from simpler components, and how interventions modify these landscapes.

32.14.4.1 Example: Composing Multiple Subsystems

The following example demonstrates how to compose multiple subsystems using Decapodes.jl, showing how composition affects the attractor landscape:

# Example: Composing predator-prey with resource dynamics
# This demonstrates how Decapodes.jl enables compositional modeling

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using Decapodes Catlab CombinatorialSpaces DiagrammaticEquations OrdinaryDiffEq ComponentArrays CairoMakie GraphMakie Graphs

# Import necessary functions
using Decapodes: infer_types!, resolve_overloads!
# Note: collate may not be directly available, but composition is conceptually shown

# Define predator-prey subsystem
# This represents a classic Lotka-Volterra model
# Note: For simplicity, we use a simplified version that can be composed
PredatorPrey = @decapode begin
  (P, H)::Form0  # Predator and prey populations
  (Ṗ, Ḣ)::Form0  # Time derivatives

  (α, β, γ, δ)::Parameter  # Interaction parameters

  # Predator-prey dynamics
  # Prey growth: Ḣ = α*H - β*P*H
  # Predator growth: Ṗ = δ*P*H - γ*P
  Ḣ == (α * H) - (β * P * H)
  Ṗ == (δ * P * H) - (γ * P)

  # Time derivatives
  ∂ₜ(H) == Ḣ
  ∂ₜ(P) == Ṗ
end

# Define resource subsystem
# This represents resource dynamics (e.g., nutrients, food availability)
Resource = @decapode begin
  R::Form0  # Resource level
  Ṙ::Form0   # Time derivative

  (r, K, c)::Parameter  # Growth rate, carrying capacity, consumption rate

  # Resource dynamics: logistic growth with consumption
  # Ṙ = r*R*(1 - R/K) - c*R
  Ṙ == (r * R * (1 - R / K)) - (c * R)

  # Time derivative
  ∂ₜ(R) == Ṙ
end

# Process the decapodes (infer types and resolve overloads)
try
    infer_types!(PredatorPrey)
    resolve_overloads!(PredatorPrey)
    println("✓ Predator-Prey decapode processed successfully")
catch e
    println("⚠ Predator-Prey processing: ", e)
end

try
    infer_types!(Resource)
    resolve_overloads!(Resource)
    println("✓ Resource decapode processed successfully")
catch e
    println("⚠ Resource processing: ", e)
end

# Composition: Couple predator-prey with resource dynamics
# The resource affects prey growth rate (more resources → faster prey growth)
# This creates a composed system with emergent attractor dynamics

# Note: Full composition would require defining coupling terms using @relation
# For illustration, we show the conceptual structure and that both decapodes are valid

println("\nSubsystem 1: Predator-Prey")
println("  Variables: P (predator), H (prey)")
println("  Attractor type: Limit cycle (oscillations)")
println("  Status: ✓ Decapode defined and processed")

println("\nSubsystem 2: Resource")
println("  Variables: R (resource)")
println("  Attractor type: Stable equilibrium")
println("  Status: ✓ Decapode defined and processed")

println("\nComposed System:")
println("  Variables: P, H, R (coupled)")
println("  Attractor type: May exhibit new attractors (e.g., chaotic)")
println("  Composition effect: Resource availability affects prey growth,")
println("    which affects predator-prey oscillations")
println("  Intervention: \$do(c = c^*)\$ changes resource consumption,")
println("    which shifts attractor landscape")
println("  Note: Full composition requires defining coupling morphisms")
println("    using @relation and collate() functions from Catlab")

# Verify both decapodes are valid and can be used
println("\n✓ Both decapodes are valid and ready for composition")
println("  PredatorPrey: ", typeof(PredatorPrey))
println("  Resource: ", typeof(Resource))

✓ Predator-Prey decapode processed successfully
✓ Resource decapode processed successfully

Subsystem 1: Predator-Prey
  Variables: P (predator), H (prey)
  Attractor type: Limit cycle (oscillations)
  Status: ✓ Decapode defined and processed

Subsystem 2: Resource
  Variables: R (resource)
  Attractor type: Stable equilibrium
  Status: ✓ Decapode defined and processed

Composed System:
  Variables: P, H, R (coupled)
  Attractor type: May exhibit new attractors (e.g., chaotic)
  Composition effect: Resource availability affects prey growth,
    which affects predator-prey oscillations
  Intervention: $do(c = c^*)$ changes resource consumption,
    which shifts attractor landscape
  Note: Full composition requires defining coupling morphisms
    using @relation and collate() functions from Catlab

✓ Both decapodes are valid and ready for composition
  PredatorPrey: SummationDecapode{Any, Any, Symbol}
  Resource: SummationDecapode{Any, Any, Symbol}

Key insights from composition:

Emergent attractors: Composed systems can have attractors that don’t exist in individual subsystems
Stability changes: Composition can make stable subsystems unstable, or vice versa
Basin modifications: Composition changes which initial conditions lead to which attractors
Intervention effects: Interventions on one subsystem affect the entire composed system

Connection to CDMs: This compositional approach enables building complex CDMs from simpler components, where each component represents a distinct causal mechanism. The composition determines how these mechanisms interact, creating the overall system dynamics and attractor landscape.

32.15 Further Reading

This book synthesises concepts from Parts I–IV
See subsequent chapters for interventional forecasting, counterfactuals, and applications
Julia packages for CDM graphs and figures: CausalDynamics / DAGMakie, process gloss in Concept Reference Tables 4 and 8
Estimation under continuous MTP / sequential / survival LMTP after identification: CausalTargeted (methods; Ch. 20–23)
Lux.jl: Neural network layers for UDEs and neural ODEs (https://lux.csail.mit.edu/)
SciMLSensitivity.jl: Adjoint sensitivity methods for efficient gradient computation (https://sensitivity.sciml.ai/stable/)
Optimization.jl: Unified optimisation interface for training UDEs (https://docs.sciml.ai/Optimisation/stable/)
Rackauckas et al. (2020): Universal differential equations and adjoint methods
Decapodes.jl: Framework for composing and simulating physical systems using diagrammatic representation (https://github.com/AlgebraicJulia/Decapodes.jl)
Morris et al. (2024): Decapodes.jl paper on diagrammatic representation of PDEs

Aitchison, John. 1986. The Statistical Analysis of Compositional Data. Chapman & Hall.

Araujo, Wellington, Luiz Oliveira, Marcio Machado, Anderson Nascimento, et al. 2025. “Quantifying First-Order Markov Violations in Noisy Reinforcement Learning: A Causal Discovery Approach.” arXiv Preprint.

Barabási, Albert-László. 2016. Network Science. Cambridge University Press.

Bareinboim, Elias, and Judea Pearl. 2013. “A General Algorithm for Deciding Transportability of Experimental Results.” Journal of Causal Inference 1 (1): 107–34.

Barrat, Alain, Marc Barthélemy, and Alessandro Vespignani. 2008. Dynamical Processes on Complex Networks. Cambridge University Press.

Díaz, Iván, Katherine L. Hoffman, and Nima S. Hejazi. 2024. “Causal Survival Analysis Under Competing Risks Using Longitudinal Modified Treatment Policies.” Lifetime Data Analysis 30: 213–36. https://doi.org/10.1007/s10985-023-09606-7.

Díaz, Iván, Nicholas Williams, Katherine L. Hoffman, and Edward J. Schenck. 2023. “Nonparametric Causal Effects Based on Longitudinal Modified Treatment Policies.” Journal of the American Statistical Association 118 (542): 846–57. https://doi.org/10.1080/01621459.2021.1955691.

Doucet, Arnaud, Nando de Freitas, and Neil Gordon. 2001. “Sequential Monte Carlo Methods in Practice.” Statistics for Engineering and Information Science (New York).

Durbin, James, and Siem Jan Koopman. 2012. Time Series Analysis by State Space Methods. 2nd ed. Oxford University Press.

Frühwirth-Schnatter, Sylvia. 2006. Finite Mixture and Markov Switching Models. Springer.

Gardiner, Crispin. 2009. Stochastic Methods: A Handbook for the Natural and Social Sciences. 4th ed. Springer.

Gelman, Andrew, Aki Vehtari, Richard McElreath, et al. 2026. Bayesian Workflow. Chapman & Hall. https://avehtari.github.io/Bayesian-Workflow/.

Hamilton, James D. 1994. Time Series Analysis. Princeton University Press.

Hernán, Miguel A., and James M. Robins. 2020. Causal Inference: What If. Chapman & Hall/CRC.

Hirsch, Morris W., Stephen Smale, and Robert L. Devaney. 2012. Differential Equations, Dynamical Systems, and an Introduction to Chaos. 3rd ed. Academic Press.

Imbens, Guido W., and Donald B. Rubin. 2015. Causal Inference in Statistics, Social, and Biomedical Sciences. Cambridge University Press.

King, Ruth, Perry de Valpine, and Rachel S. McCrea. 2016. “Statistical Ecology.” Annual Review of Statistics and Its Application 3: 401–26.

Koller, Daphne, and Nir Friedman. 2009. Probabilistic Graphical Models: Principles and Techniques. MIT Press.

Morris, Luke, Andrew Baas, Evan Arias, Micah Gatlin, Evan Patterson, and James Fairbanks. 2024. “Decapodes: A Diagrammatic Tool for Representing, Composing, and Computing Spatialized Partial Differential Equations.” Journal of Computational Science 81: 102345. https://doi.org/10.1016/j.jocs.2024.102345.

Newman, Mark E. J. 2018. Networks. 2nd ed. Oxford University Press.

Øksendal, Bernt. 2013. Stochastic Differential Equations: An Introduction with Applications. 6th ed. Springer.

Pearl, Judea. 2009. Causality: Models, Reasoning, and Inference. 2nd ed. Cambridge University Press.

Pearl, Judea, and Elias Bareinboim. 2014. “External Validity: From Do-Calculus to Transportability Across Populations.” Statistical Science 29 (4): 579–95.

Porter, Mason A., and James P. Gleeson. 2016. Dynamical Systems on Networks: A Tutorial. Springer.

Rackauckas, Christopher, Yingbo Ma, Julius Martensen, et al. 2020. “Universal Differential Equations for Scientific Machine Learning.” arXiv Preprint arXiv:2001.04385. https://arxiv.org/abs/2001.04385.

Robins, James M., Miguel A. Hernán, and Babette Brumback. 2000. “Marginal Structural Models and Causal Inference in Epidemiology.” Epidemiology 11 (5): 550–60.

Särkkä, Simo. 2013. Bayesian Filtering and Smoothing. Cambridge University Press.

Shen, Zhichao, Jing Liu, Yong Jiang, Zhaoran Chen, et al. 2019. “Causal Discovery and Forecasting in Nonstationary Environments with State-Space Models.” PMC / Proceedings of Machine Learning Research.

Strogatz, Steven H. 2014. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. 2nd ed. Westview Press.

Vaswani, Ashish, Noam Shazeer, Niki Parmar, et al. 2017. “Attention Is All You Need.” Advances in Neural Information Processing Systems 30. https://proceedings.neurips.cc/paper/2017/hash/3f5ee243547dee91fbd053c1c4a845aa-Abstract.html.

Whitehead, Alfred North. 1978. Process and Reality: An Essay in Cosmology. Edited by David Ray Griffin and Donald W. Sherburne. Free Press.