Partially Functional Dynamic Backdoor Diffusion-based Causal Model

TL;DR

This paper introduces PFD-BDCM, a novel diffusion-based causal model that effectively handles dynamic confounders and functional data in complex spatio-temporal systems, improving causal inference accuracy.

Contribution

The paper develops a unified generative framework that captures spatio-temporal dependencies and functional variables, with theoretical guarantees and practical improvements over existing methods.

Findings

Outperforms existing methods on synthetic and real-world data

Provides theoretical guarantees on causal effect preservation

Demonstrates effectiveness in observational, interventional, and counterfactual analyses

Abstract

Causal inference in spatio-temporal settings is critically hindered by unmeasured confounders with complex spatio-temporal dynamics and the prevalence of multi-resolution data. While diffusion models present a promising avenue for estimating structural causal models, existing approaches are limited by assumptions of causal sufficiency or static confounding, failing to capture the region-specific, temporally dependent nature of real-world latent variables or to directly handle functional variables. We bridge this gap by introducing the Partially Functional Dynamic Backdoor Diffusion-based Causal Model (PFD-BDCM), a unified generative framework designed to simultaneously tackle causal inference with dynamic confounding and functional data. Our approach formalizes a novel structural causal model that captures spatio-temporal dependencies in latent confounders through conditional autoregressive processes, represents functional variables via basis expansion coefficients treated as standard graph nodes, and integrates valid backdoor adjustment into a diffusion-based generative process. We provide theoretical guarantees on the preservation of causal effects under basis expansion and derive error bounds for counterfactual estimates. Experiments on synthetic data and a real-world air pollution case study demonstrate that PFD-BDCM outperforms existing methods across observational, interventional, and counterfactual queries. This work provides a rigorous and practical tool for robust causal inference in complex spatio-temporal systems characterized by non-stationarity and multi-resolution data.

Tables

Table 1. Table 1: Mean ( × 10 − 2 \times 10^{-2} ) ± \pm standard deviation of MMD 2 of models compared to the true target distribution (Observation query)

Variable	PFD-DCM	PFD-BDCM	PF-DSEM-Bays	Variable	PFD-DCM	PFD-BDCM	PF-DSEM-Bays
SO2	$1.130\pm 0.008$	$1.155\pm 0.008$	$13.919\pm 0.131$	PM10	$1.133\pm 0.008$	$1.141\pm 0.008$	$11.961\pm 0.156$
NOx	$1.135\pm 0.008$	$1.196\pm 0.009$	$17.854\pm 0.157$	PM2.5	$1.168\pm 0.008$	$1.125\pm 0.008$	$12.584\pm 0.163$
CO	$1.121\pm 0.008$	$1.154\pm 0.009$	$12.311\pm 0.156$	BC	$1.150\pm 0.008$	$1.180\pm 0.008$	$16.745\pm 0.186$
VOC	$1.208\pm 0.008$	$1.182\pm 0.008$	$24.777\pm 0.149$	OC	$1.218\pm 0.008$	$1.192\pm 0.009$	$21.109\pm 0.182$
NH3	$1.193\pm 0.008$	$1.142\pm 0.008$	$36.627\pm 0.151$	CO2	$1.146\pm 0.008$	$1.112\pm 0.008$	$21.793\pm 0.178$

Equations

$${\text{x}}_{k,ij}=f_{k,ij}({\mathbf{x}}_{{\mathbb{PA}}_k,ij},{\text{u}}_{k,ij}),$$

$${\mathbf{x}}_{{\mathbb{C}}_2,ij}={\mathbf{\Gamma }}_i\mathbf{G}\left({\mathbf{x}}_{{\mathbb{C}}_1,ij}\right)+{\mathbf{u}}_{{\mathbb{C}}_2,ij},$$

$$\mathrm{Cov}({\mathbf{x}}_{{\mathbb{C}}_1,i})={\mathbf{D}}_{{\mathbb{C}}_1,i}\otimes \mathrm{Cov}({\mathbf{x}}_{{\mathbb{C}}_1,ij}),\mathrm{Cov}({\mathbf{u}}_{{\mathbb{C}}_2,i})={\mathbf{D}}_{{\mathbb{C}}_2,i}\otimes \mathrm{Cov}({\mathbf{u}}_{{\mathbb{C}}_2,ij}),$$

$${\mathbf{D}}_{{\mathbb{C}}_h,i}={\left({\mathbf{I}}_J-{\rho }_{{\mathbb{C}}_h,i}{\mathbf{H}}_{{\mathbb{C}}_h,i}\right)}^{-1},$$

$${\text{x}}_{k,ij}=f_k({\mathbf{x}}_{{\mathbb{PA}}_k,ij},{\text{u}}_{k,ij}),$$

$${\text{x}}_{m,ij}=\int b_m(t){\text{z}}_{ij}(t)\mathrm{d}t,$$

$${\text{x}}_k^t=$$

$$+{\sigma }_t{ϵ}_k,{ϵ}_k\sim \mathcal{N}(0,\mathbf{I}),\text{for}t\in [T].$$

$${\hat{\text{x}}}_k^{t-1}=$$

$$L_{\gamma }({ϵ}_{k,\theta })$$

$$g:\mathcal{X}\times {\mathcal{X}}_{\mathbb{PB}}\to \mathcal{U}(\text{or}{g}^{(z)}:{\mathcal{X}}^{K_n}\times {\mathcal{X}}_{\mathbb{PB}}\to {\mathcal{U}}^{K_n})\text{(encoding function)};$$

$$h:\mathcal{U}\times {\mathcal{X}}_{\mathbb{PB}}\to \mathcal{X}(\text{or}{h}^{(z)}:{\mathcal{U}}^{K_n}\times {\mathcal{X}}_{\mathbb{PB}}\to {\mathcal{X}}^{K_n})\text{(decoding function)},$$

$${\Vert c(t)-\sum _{m=1}^{K_n}{\varphi }_m{b}_m(t)\Vert }_{\mathcal{H}}\le C_c{K}_n^{-s_c},\mathbb{E}\left[{\Vert \text{z}(t)-\sum _{m=1}^{K_n}{\text{x}}_m{b}_m(t)\Vert }_{\mathcal{H}}^2\right]\le C_z{K}_n^{-2s_z},$$

$$\Vert z_0(t)-\sum _{m=1}^{K_n}{x}_{0,m}{b}_m(t){\Vert }_{\mathcal{H}}\le C_{z_0}{K}_n^{-s_{z_0}},$$

$$\Vert f(\dots ,⟨c(t),z_0(t)⟩,\dots ,\text{u})-f(\dots ,{\mathbf{\varphi }}^{\top }{\mathbf{x}}_0^{(z)},\dots ,\text{u})\Vert$$

$$\le$$

Full text

Keywords: Causal inference; Spatio-temporal statistics; Functional data; Unmeasured confounder

1 Introduction

Causal inference fundamentally addresses interventional and counterfactual questions that go beyond statistical correlations, proving valuable in high-stakes domains like healthcare for treatment effect estimation (Hill,, 2011), policy evaluation without randomized trials (LaLonde,, 1986). The field’s core challenge stems from the fundamental problem of causal inference: the impossibility of simultaneously observing both an outcome under a treatment and the potential outcomes under the alternative treatment (or control) for the same unit (Imbens and Rubin,, 2015), necessitating methods to overcome confounding bias in observational data. Traditional approaches, including potential outcomes frameworks (Imbens and Rubin,, 2015), propensity scoring (Rosenbaum and Rubin,, 1983), instrumental variables (Angrist et al.,, 1996), and structural causal models (SCMs) (Pearl,, 2009), exhibit significant limitations when handling modern complex datasets, as they struggle with high-dimensional confounders, ethical constraints of randomized trials, scarcity of valid instruments, and the requirement of known causal graphs. These limitations become particularly acute when confounders involve high-dimensional data like medical images or partially unobserved variables (Shalit et al.,, 2017).

Within the Structural Causal Model (SCM) framework, causal queries can be answered by learning a proxy for the unobserved exogenous noise and the structural equations (Pearl,, 2009). This suggests that (conditional) generative models that encode to a latent space could be an option for modeling SCMs, as the latent space serves as proxies for exogenous noises. Recent advances have explored the integration of deep generative models with structural causal models to address these challenges. Causal Normalizing Flows (CNFs) (Khemakhem et al.,, 2021) employ invertible transformations but face limitations in high-dimensional settings. Diffusion-based SCMs (Diff-SCM) (Mamaghan et al.,, 2023) and Decoupled Flows (DecaFlow) (Almodóvar et al.,, 2025) represent recent advances but lack explicit handling of spatio-temporal unmeasured confounding. Diffusion-based Causal Models (DCM) (Chao et al.,, 2024) and their extension incorporating backdoor adjustment (BDCM) (Shimizu,, 2023) can learn complex structural equations from data. Nevertheless, both methods rely on static assumptions and thus do not fully capture the spatio-temporal structure inherent in confounding variables, a key characteristic of real-world systems, where such factors often demonstrate complex dependencies.

To address this gap, we propose the Partially Functional Dynamic Backdoor Diffusion-based Causal Model (PFD-BDCM), a unified generative framework designed to simultaneously tackle causal inference with dynamic confounding and functional data. We first formalize the Partially Functional Spatio-Temporal Dynamic SCM (PFST-DSCM), which represents unmeasured confounders via dynamic equations capturing spatial heterogeneity and temporal dependence, while incorporating functional variables through basis expansion, treating their coefficients as standard nodes in the causal graph. This formalization enables precise problem specification previously unattainable. Building on this, we develop Partially Functional Dynamic Backdoor Diffusion-based Causal Model (PFD-BDCM) for learning and inference through the Denoising Diffusion Implicit Models with Additional Conditioning (AC-DDIM). Unlike prior diffusion causal models (DCM and BDCM) that condition only on parent variables or a static backdoor adjustment set, AC-DDIM learns the distribution of a node given both its parents and a valid backdoor adjustment set, designed to directly adjust for bias from dynamically evolving unobserved variables. During training, the denoising network for each node receives the parent variables, backdoor adjustment variables and the noisy node value, conditioning the generative process to simulate removal of confounding even from unmeasured dynamic sources. For functional nodes, diffusion operates seamlessly on their finite-dimensional coefficient vectors. We provide theoretical guarantees for causal effect preservation under basis expansion and error bounds for counterfactual estimates. Through extensive simulations and a real-world case study on China’s air pollution (2015-2020), we demonstrate PFD-BDCM’s superior performance and practical utility for policy-relevant causal queries.

The remainder of this paper is organized as follows. Section 2 introduces necessary notations and reviews necessary background, Section 3 provides the PFST-DSCM and the PFD-BDCM architecture and algorithms. Section 4 presents our theoretical guarantees, which prove that causal effects are preserved under basis expansion and derive explicit error bounds linking the model’s counterfactual estimation error to its empirical reconstruction fidelity, providing a rigorous foundation for the proposed methodology. In Section 5, we report on extensive simulation studies that evaluate PFD‑BDCM across a variety of data‑generating mechanisms, highlighting its advantages in handling dynamic confounding and functional data. Section 6 illustrates the practical utility of the framework through a real-world case study on atmospheric pollutant emissions in China (2015-2020). The paper concludes in Section 7 with a summary of contributions, a discussion of limitations, and an outline of directions for future research.

2 Notation and Fundamental

We first introduce some useful notations, fundamental concepts of causal inference.

Notations. Let $[n]:={1,\dots,n}$ index spatial units (regions) and $[J]:={1,\dots,J}$ index time points, $[K]:={1,\dots,K}$ index nodes in a directed acyclic graph (DAG) $\mathcal{G}$. For node $k$, denote its parent set by $\mathbb{PA}_{k}$ and a valid backdoor adjustment set by $\mathbb{B}_{k}$. The corresponding endogenous variables are $\mathbf{x}_{\mathbb{PA}_{k}}$ and $\mathbf{x}_{\mathbb{B}_{k}}$. Each node is associated with an endogenous variable $\textnormal{x}_{k}\in\mathbb{R}^{d_{k}}$ and an exogenous noise $\textnormal{u}_{k}\in\mathbb{R}^{d_{k}}$. For functional nodes, $d_{k}=K_{n}$, representing the number of basis functions. And for scalar nodes, $d_{k}=1$. In our diffusion process, $\hat{\textnormal{u}}_{k}^{t}$ and $\hat{\textnormal{x}}_{k}^{t}$ denote the latent and reconstructed variables at step $t$.

In causal inference, a confounder denotes a variable that causally influences both a treatment variable x and an outcome variable y, thereby inducing a non-causal association between them. Observable Confounders refer to confounders that can be measured, which permit adjustment through statistical methods such as stratification, matching or regression (Pearl,, 2009). Unobservable confounders denote latent variables that fulfill confounding criteria but resist direct measurement, which can potentially bias causal estimates when unaccounted for. Unobserved explanatory (explained) nodes are unobserved confounder nodes that have no nodes (unobserved confounder nodes) as its parent (descendant).

Structural Causal Models. Consider a DAG (such as Fig. 1) $\mathcal{G}$ with nodes $[K]$ in a topologically sorted order 111A topological order is a linear arrangement of variables where a variable appears after all its direct causes (parents) (Pearl,, 2009)., where a node $k$ is built with an endogenous random variable $\textnormal{x}_{k}$ defined on a space $\mathcal{X}_{k}\subset\mathbb{R}^{d_{k}}$, which has a random exogenous input $\textnormal{u}_{k}$. A structural causal model (SCM) $\mathcal{M}$ characterizes the relationship between the endogenous variable $\textnormal{x}_{k}$ of a node $k$ with the endogenous variables of its parents $\mathbf{x}_{\mathbb{PA}_{k}}$ and its own exogenous variable $\textnormal{u}_{k}$. Formally, we define $\mathcal{M}:=(\mathbf{f}(\mathbf{x},\mathbf{u}),p_{\mathbf{u}})$, where $\mathbf{f}(\mathbf{x},\mathbf{u})$ specifies how entire endogenous variables $\mathbf{x}:=\{\textnormal{x}_{1},\cdots,\textnormal{x}_{K}\}$ are generated from the set of exogenous random variables $\mathbf{u}:=\{\textnormal{u}_{1},\cdots,\textnormal{u}_{K}\}$ with a prior distribution $p_{\mathbf{u}}$. The structural mechanism is governed by $\mathbf{f}(\mathbf{x},\mathbf{u}):=(f(\mathbf{x}_{\mathbb{PA}_{1}},\textnormal{u}_{1}),\cdots,f(\mathbf{x}_{\mathbb{PA}_{K}},\textnormal{u}_{K}))$, where each $\textnormal{x}_{k}:=f(\mathbf{x}_{\mathbb{PA}_{k}},\textnormal{u}_{k})$ for $k\in[K]$ (Pearl,, 2009).

3 Methods

We begin by formalizing the data-generating process underlying complex spatio‑temporal data through a Spatio‑Temporal Dynamic Structural Causal Model (ST‑DSCM). Building on this foundation, we then introduce the Partially Functional Dynamic Backdoor Diffusion‑based Causal Model (PFD‑BDCM), which is designed to answer causal queries under the extended PFST‑DSCM framework.

3.1 PFST-DSCM

Consider a spatio-temporal dataset $\mathbf{x}:=\{\textnormal{x}_{k}\}_{k\in[K]}$ containing $K$ variables across $n$ regions over $J$ time points. Using $i$ to index regions and $j$ for time points, $\textnormal{x}_{k}=(x_{k,ij})_{n\times J}$ represent the value of the $k$-th variable $\textnormal{x}_{k}$ over the time points and regions. We use a DAG to characterize the causal relationships among $\{\textnormal{x}_{k}\}_{k\in[K]}$. Consider a DAG $\mathcal{G}$ with nodes $[K]$ and a topologically sorted order such that each node $k$ has the $\textnormal{x}_{k}$ as the random variable. Let $\mathbb{C}_{1}$ and $\mathbb{C}_{2}\subseteq[K]$ be two distinct sets of nodes with unobserved confounders, where $\mathbb{C}_{1}$ designates a set of unobserved explanatory nodes and $\mathbb{C}_{2}$ denotes a set of unobserved explained nodes. For $h=1,2$, let $\mathbf{x}_{\mathbb{C}_{h}}:=\{\mathbf{x}_{{\mathbb{C}_{h}},ij}\}_{i\in[n],j\in[J]}:=\{x_{q,ij}\}_{i\in[n],j\in[J],q\in\mathbb{C}_{h}}$ and $\mathbf{u}_{\mathbb{C}_{h}}:=\{\mathbf{u}_{{\mathbb{C}_{h}},ij}\}_{i\in[n],j\in[J]}:=\{u_{q,ij}\}_{i\in[n],j\in[J],q\in\mathbb{C}_{h}}$.

By the definition of $\mathbb{C}_{1}$, every node in this set satisfies the backdoor criterion222A set of node $\mathbb{B}$ satisfies backdoor criterion (Pearl et al.,, 2016) for tuple $(\textnormal{x},\textnormal{y})$ in DAG $\mathcal{G}$ if no node in $\mathbb{B}$ is a descendant of x and $\mathbb{B}$ blocks all paths between x (cause) and y (outcome) which contains an arrow into x.. For nodes with unmeasured confounders ($k\in\mathbb{C}_{2}$, $\mathbb{PA}_{k}\in\mathbb{C}_{1}$), we assume spatio-temporal dynamic SCM relations:

[TABLE]

where $f_{k,ij}$ represents a class of potentially linear or nonlinear spatio-temporal functions, and $\textnormal{u}_{k,ij}$ denotes independent exogenous noise terms.

Specifically, to account for spatial heterogeneity, we posit a relationship such that

[TABLE]

where $\boldsymbol{\Gamma}_{i}$ is a $\dim(\mathbf{x}_{{\mathbb{C}_{2}},ij})\times\dim(\mathbf{x}_{{\mathbb{C}_{1}},ij})$ structural coefficient matrix varying across regions; $\boldsymbol{G}(\cdot)=(G_{1}(.),\cdots,G_{q}(.))^{\top}$ is a $q\times 1$ nonzero vector-valued function with differentiable functions $G_{1},\cdots,G_{q}$, and $q=\dim(\mathbf{x}_{{\mathbb{C}_{1}},ij})$.

To capture the temporal dependence, let $\mathbf{x}_{{\mathbb{C}_{1}},i}=(\mathbf{x}_{{\mathbb{C}_{1}},i1}^{\top},\cdots,\mathbf{x}_{{\mathbb{C}_{1}},iJ}^{\top})^{\top}$ and $\mathbf{u}_{{\mathbb{C}_{2}},i}=(\mathbf{u}_{{\mathbb{C}_{2}},i1}^{\top},\cdots,\mathbf{u}_{{\mathbb{C}_{2}},iJ}^{\top})^{\top}$, the covariance matrices $\mathrm{Cov}(\mathbf{x}_{{\mathbb{C}_{1}},i})$ and $\mathrm{Cov}(\mathbf{u}_{{\mathbb{C}_{2}},i})$ of $\mathbf{x}_{{\mathbb{C}_{1}},i}$ and $\mathbf{u}_{{\mathbb{C}_{2}},i}$ can be expressed as

[TABLE]

where $\boldsymbol{D}_{{\mathbb{C}_{h}},i}(h\in\{1,2\})$ are the $J\times J$ adjacent time covariance matrices, $\mathrm{Cov}(\cdot)$ represent between-variable covariances. In order to establish a rigorous framework for the temporal adjacency structure, the conditional autoregressive (CAR) model (Besag,, 1974) is adopted with the forms

[TABLE]

where $\rho_{{{\mathbb{C}_{h}},i}}(h\in\{1,2\})$ are adjacent time association parameters, and $\boldsymbol{H}_{{{\mathbb{C}_{h}},i}}(h\in\{1,2\})$ are $J\times J$ adjacency matrices in which the element $h_{jl}=1$ implies that time $l$ is adjacent to time $j$ and otherwise $h_{jl}=0$. The covariance kernel $\boldsymbol{D}_{\mathbb{C}_{1},i}$ explicitly models temporal dependency by inducing correlation across neighboring time points. This represents a fundamental departure from the standard independent and identically distributed (i.i.d.) assumption typically placed on latent variables in conventional causal generative models, allowing our framework to capture the realistic, evolving nature of unmeasured confounders.

For other endogenous variables, we assume standard SCM relations

[TABLE]

where $f_{k}$ are potentially linear (or nonlinear) static functions.

We refer to the model satisfying the above specifications as the Spatio-Temporal Dynamic Structural Causal Model (ST-DSCM). It provides a principled framework for capturing spatio-temporally dynamic unmeasured confounding. However, Real-world applications often present significant challenges in elucidating causal relationships among variables observed under heterogeneity. To address these complexities, we extend the ST-DSCM by incorporating functional random variables $\textnormal{z}_{ij}(t)$. We adopt a basis expansion framework to achieve dimensionality reduction in the functional space via a set of orthogonal basis functions $\{\mathbf{b}_{m}(t)\}_{m=1}^{K_{n}}$. Define

[TABLE]

for $i\in[n],j\in[J],m\in[K_{n}]$. Let $\mathbf{x}_{k,ij}=(\textnormal{x}_{k,1,ij},\dots,\textnormal{x}_{k,K_{n},ij})^{\top}$ denote the collection of all basis coefficients derived from a functional variable $\textnormal{z}_{ij}(t)$. This technique projects infinite-dimensional functional covariates onto a finite-dimensional space while preserving functional characteristics. The mathematical details of the basis expansion and its properties are provided in Supplementary LABEL:sec:_FDR. Crucially, for causal inference, this coefficient vector $\mathbf{x}_{k,ij}$ is treated as a single multidimensional endogenous node within the causal graph. This representation allows all standard causal operations (intervention, counterfactual reasoning, and estimation) to be applied uniformly, whether the node originates from a scalar or a functional measurement.

By augmenting the ST-DSCM with this representation, we obtain the Partially Functional Spatio-Temporal Dynamic Structural Causal Model (PFST-DSCM) 333For example, Fig. 2 is a PFST-DSCM with 11 exogenous and endogenous nodes (where nodes $\textnormal{x}_{3,ij}$, $\textnormal{x}_{4,ij}$ and $\textnormal{x}_{9,ij}$ are unmeasured confounder nodes with spatial heterogeneity and temporal dependencies, $\textnormal{z}_{ij}(t)$ is functional nodes, and $\mathbf{x}_{2,ij}$ is the corresponding base expansion vector node). In this model, nodes can be scalar, dynamic latent confounders (as in the ST-DSCM), or finite-dimensional vectors representing functional variables. This extension preserves the dynamic modeling of unmeasured confounding while seamlessly integrating multi-resolution (scalar and functional) data.

3.2 Identification of the PFST-DSCM

In causal inference, defining the model as a Directed Acyclic Graph (DAG) (Assumption 2(i)) is a standard and necessary foundation (Pearl,, 2009). The designation “Partially Functional Spatio-Temporal Dynamic” accurately captures the specific composition of our graph $\mathcal{G}$, which integrates standard scalar nodes, dynamically modeled latent confounder nodes, and multi-dimensional nodes representing functional data via basis coefficients. This precise nomenclature is essential to distinguish our framework from conventional DAGs or those used in static causal models like DCM or BDCM. Crucially, our model identification relies on a set of core assumptions (Assumption 1 and Assumption 2((ii), (iii)). A detailed discussion justifying the plausibility of these assumptions is provided in the Supplementary LABEL:sec:_PFST-DSCM-sup.

Assumption 1**.**

For CAR model (Eq.(4)), $\rho_{\mathbb{C}_{h}i}$ must satisfy the following constrained conditions $\rho_{\mathbb{C}_{h}i}^{L}<\rho_{\mathbb{C}_{h}i}<\rho_{\mathbb{C}_{h}i}^{U}$ to ensure that $\boldsymbol{D}_{\mathbb{C}_{h}i}$ are positive definite, where $\rho_{\mathbb{C}_{h}i}^{L}=\min(1/\tau_{\mathbb{C}_{h}J},0),\rho_{\mathbb{C}_{h}i}^{U}=1/\tau_{\mathbb{C}_{h}1}$, and $\tau_{\mathbb{C}_{h}1}\geq\cdots\geq\tau_{\mathbb{C}_{h}J}$ are the eigenvalues of neighborhood matrices $\boldsymbol{H}_{\mathbb{C}_{h}i}$ (Song et al.,, 2012).

Assumption 2**.**

(i) The Partially Functional Spatio-Temporal Dynamic Directed Acyclic Graph (PFST-DDAG) $\mathcal{G}$ is the graph induced by PFST-DSCM $\mathcal{M}$.

(ii) The unobserved random variables are jointly independent (Markovian SCM), and every PFST-DSCM $\mathcal{M}$ entails a unique joint observational distribution satisfying the causal Markov assumption: $q(\mathbf{x})=\prod_{k=1}^{K}q(\textnormal{x}_{k}|\mathbf{x}_{\mathbb{PA}_{k}})$.

(iii) For any observed node with unobserved parents, a valid backdoor adjustment set $\mathbb{B}_{k}$ can be identified from $\mathcal{G}$, enabling causal identification even in the presence of dynamic unmeasured confounders under Assumption 1.

The PFST-DSCM provides a unified formalization that captures our motivating air pollution setting: confounders that evolve in space and time, and drivers that may be functional in nature. The model operates under the standard assumption that the causal graph is known (or estimable), and all structural parameters can be estimated from data using existing Bayesian spatio-temporal methods (Song et al.,, 2012; Tang et al.,, 2017). With this formalization in place, the framework supports the full range of causal queries (observational, interventional, and counterfactual) from observational data.

3.3 PFD-BDCM

We now present the PFD-BDCM, a model designed to handle causal queries under the PFST-DSCM framework, which explicitly accounts for unmeasured confounders with spatial heterogeneity and temporal dependence. The model employs an encoder-decoder architecture to enable causal reasoning across multi-resolution variables.

3.3.1 Data-Generating

The data-generating process of the PFD-BDCM is formalized as follows. i) If node $k$ and its parent nodes are unobservable confounders, the data $\mathbf{x}_{k}\in\mathbb{R}^{d_{k}}$ are generated via Eq. (1). These nodes are not themselves targets of training, their data serve only to provide structural information and to adjust for bias during the training of observable nodes of interest. Consequently, no diffusion model is trained on data from such nodes. ii) For any other node $k$, data are generated using Eq. (5). When node $k$ and all its parent nodes are observable, an encoder $g_{k}$ maps $(\mathbf{x}_{\mathbb{PA}_{k}},\mathbf{x}_{k})$ to a latent variable $\widehat{\mathbf{u}}_{k}:=g_{k}(\mathbf{x}_{\mathbb{PA}_{k}},\mathbf{x}_{k})$, which captures the information of the exogenous noise $\mathbf{u}_{k}$. A decoder $h_{k}$ then reconstructs $\mathbf{x}_{k}$ as $\widehat{\mathbf{x}}_{k}:=h_{k}(\mathbf{x}_{\mathbb{PA}_{k}},\widehat{\mathbf{u}}_{k})$. This setup corresponds to the data-generating process of the DCM (Chao et al.,, 2024). And when node $k$ is observable but some of its parent nodes are unobservable confounders, we introduce a backdoor adjustment set $\mathbb{B}_{k}$ for node $k$ as a proxy for the unobserved parents. The encoder $g_{k}$ then maps $(\mathbf{x}_{\mathbb{PA}_{k}},\mathbf{x}_{\mathbb{B}_{k}},\mathbf{x}_{k})$ to $\widehat{\mathbf{u}}_{k}:=g_{k}(\mathbf{x}_{\mathbb{PA}_{k}},\mathbf{x}_{\mathbb{B}_{k}},\mathbf{x}_{k})$, and the decoder $h_{k}$ reconstructs $\mathbf{x}_{k}$ as $\widehat{\mathbf{x}}_{k}:=h_{k}(\mathbf{x}_{\mathbb{PA}_{k}},\mathbf{x}_{\mathbb{B}_{k}},\widehat{\mathbf{u}}_{k})$. Conditioning on the backdoor adjustment set $\mathbf{x}_{\mathbb{B}_{k}}$ provides a mechanism to account for confounding.

3.3.2 AC-DDIM

Our approach leverages diffusion models to learn the structural equations of a causal system, representing each endogenous variable $\textnormal{x}_{k}$ with a dedicated conditional diffusion model. For counterfactual reasoning, we require a deterministic mapping between observations and latent codes. Denoising Diffusion Implicit Models (DDIMs) (Song et al.,, 2021) provide such a deterministic non-Markovian reverse process, essential for identifiability. While standard conditional diffusion models condition only on parent variables $\mathbf{x}_{\mathbb{PA}_{k}}$, and existing backdoor-adjusted models condition on a backdoor set $\mathbf{x}_{\mathbb{B}_{k}}$, AC-DDIM conditions on the combined set $\mathbf{x}_{\mathbb{PB}_{k}}=\mathbf{x}_{\mathbb{PA}_{k}}\cup\mathbf{x}_{\mathbb{B}_{k}}$. This ensures the generative process for $\textnormal{x}_{k}$ is conditioned on the proper adjustment set, accounting for dynamically generated unmeasured confounders. During training, the denoising network learns to reconstruct $\textnormal{x}_{k}$ given $\mathbf{x}_{\mathbb{PB}_{k}}$; during sampling, providing values for $\mathbf{x}_{\mathbb{PB}_{k}}$ yields statistically adjusted samples.

Formally, for each node $k\in[K]$, the latent variable $\hat{\textnormal{u}}_{k}:=\textnormal{x}_{k}^{T}$ is generated through the forward implicit diffusion process

[TABLE]

This latent representation $\hat{\textnormal{u}}_{k}$ serves as a proxy for the exogenous noise $\textnormal{u}_{k}$.

The reconstruction $\hat{\textnormal{x}}_{k}:=\hat{\textnormal{x}}_{k}^{0}$ is obtained via the reverse implicit diffusion process

[TABLE]

for $t=T,\cdots,1$, initialized with $\hat{\textnormal{x}}^{T}_{k}:=\hat{\textnormal{u}}_{k}$. Since AC-DDIM’s generative process is non-Markovian and deterministic (when $\sigma_{t}=0$), timesteps can be skipped for acceleration (see Supplementary LABEL:sec:_AC-DDIM-sup).

The parameters $\theta$ are learned by minimizing the following loss function

[TABLE]

where $\gamma:=[\gamma_{\tau_{s}},\ldots,\gamma_{\tau_{1}}]$ is a vector of positive weighting coefficients that depends on $\alpha_{\tau_{S}:\tau_{1}}$.

For functional nodes represented via basis expansion, the diffusion process operates on the finite-dimensional coefficient vectors in $\mathbb{R}^{K_{n}}$. The derivation of AC-DDIM follows the standard DDIM (Song et al.,, 2021) framework with the extension of conditioning on additional covariates, complete mathematical formulations and theoretical details are provided in Supplementary LABEL:sec:_AC-DDIM-sup.

The encoding and decoding functions for node $k$ are denoted as $\mathrm{Enc}_{k}(\mathbf{x}_{\mathbb{PB}_{k}},\textnormal{x}^{0}_{k})$ (Eq. (7)) and $\mathrm{Dec}_{k}(\mathbf{x}_{\mathbb{PB}_{k}},\hat{\textnormal{u}}_{k})$ (Eq.(8)) respectively, the pseudocodes are provided in Supplementary LABEL:Enc&Dec (Algorithm LABEL:algo:_PFD-BDCM_encoding and LABEL:algo:_PFD-BDCM_decoding ).

3.3.3 Training PFD-BDCMs and Answering Causal Queries

The comprehensive training methodology (Algorithm 1) incorporates backdoor adjustment sets as covariates while training distinct diffusion models per scale node $\textnormal{x}_{k}$ (or functional node vector $\mathbf{x}_{k}$ ). Crucially, generative models for endogenous nodes exhibit mutual independence during training, thereby enabling parallelized optimization. This parallelism is feasible since each diffusion model necessitates only its target node’s values and corresponding backdoor adjustment set values. The final PFD-BDCM architecture integrates these $K$ trained diffusion models $\{\epsilon_{k,\theta}\}_{k\in[K]}$.

We now elucidate the methodology for leveraging trained PFD-BDCMs to approximate diverse causal queries. Resolution of observational and interventional queries necessitates sampling from their respective observational and interventional distributions. Counterfactual queries, however, operate at unit granularity by modifying structural equation assignments while preserving the latent exogenous noise variables consistent with empirical observations.

Generating samples for observational/interventional queries. Interventional queries concern the causal effect of actively setting a variable to a specific value, which is formally represented by the $do$-operator, as in $p_{\theta}(\textnormal{x}|\mathrm{do}(\textnormal{x}_{\mathcal{L}}:=\boldsymbol{\gamma}))$ . Specifically, to generate samples approximating the interventional distribution $p_{\theta}(\textnormal{x}_{k}|\hat{\mathbf{x}}_{\mathbb{PB}_{k},-\mathcal{L}},\mathrm{do}(\mathbf{x}_{\mathcal{L}}:=\boldsymbol{\gamma}))$ using a trained PFD-BDCM model, where $\hat{\mathbf{x}}_{\mathbb{PB}_{k},-\mathcal{L}}$ denotes $\hat{\mathbf{x}}_{\mathbb{PB}_{k}}$ with the intervention node variables $\textnormal{x}_{\mathcal{L}}$ removed. We implement the following procedure: i) For intervened node $l\in\mathcal{L}$, we set $\hat{\textnormal{x}}_{l}:=\gamma_{l}$ deterministically; ii) For root nodes $k$, sample $\hat{\textnormal{x}}_{k}$ from empirical training distributions; iii) For non-intervened nodes $k\notin\mathcal{L}$, sample latent variable $\hat{\textnormal{u}}_{k}\sim\mathcal{N}(0,1)$ ( or $\hat{\textnormal{u}}_{k}\sim t_{\nu}(0,1)$, where $t_{\nu}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ is a multivariate t‑distribution with $\nu$ degrees of freedom, location vector $\boldsymbol{\mu}$, and scale matrix $\boldsymbol{\Sigma}$ ). And subsequently compute $\hat{\textnormal{x}}_{k}:=\mathrm{Dec}_{k}(\hat{\mathbf{x}}_{\mathbb{PB}_{k}},\hat{\textnormal{u}}_{k})$ utilizing inductively generated parent variable values and backdoor adjustment set values. Generated values propagate to descendant nodes as parent and backdoor adjustment inputs. Specially, observational sampling $p_{\theta}(\textnormal{x}_{k}|\mathbf{x}_{\mathbb{PB}_{k}})$ corresponds to $\mathcal{L}=\emptyset$. The pseudocode formalized in Supplementary LABEL:ACQ (Algorithm LABEL:algo:_PFD-BDCM_intervention_sampling).

Counterfactual Queries concern hypothetical scenarios given actual observed outcomes, and require three steps: abduction (inferring exogenous noise), action (modifying equations), and prediction (simulating new outcomes). Specifically, to compute counterfactual estimates $\hat{\boldsymbol{x}}^{\mathrm{CF}}:=(\boldsymbol{x}_{1}^{\mathrm{CF}},\cdots,\boldsymbol{x}_{K}^{\mathrm{CF}})$ , within the PFD-BDCM framework, given factual observation $\boldsymbol{x}^{\mathrm{F}}:=(\boldsymbol{x}_{1}^{\mathrm{F}},\cdots,\boldsymbol{x}_{K}^{\mathrm{F}})$ , and intervention set $\mathcal{L}$ with values $\boldsymbol{\gamma}$, we implement the following systematic procedure: i) For intervened nodes $l\in\mathcal{L}$, assign $\hat{\boldsymbol{x}}_{l}^{\mathrm{CF}}:=\boldsymbol{\gamma}_{l}$ deterministically; ii) For non-intervened descendant nodes $k$, using inductively generated parent and backdoor adjustment estimates $\hat{\boldsymbol{x}}_{\mathbb{PB}_{k}}^{\mathrm{CF}}$, we compute $\hat{\boldsymbol{x}}_{k}^{\mathrm{CF}}:=\mathrm{Dec}_{k}(\hat{\boldsymbol{x}}_{\mathbb{PB}_{k}}^{\mathrm{CF}},\mathrm{Enc}_{k}(\boldsymbol{x}_{\mathbb{PB}_{k}}^{\mathrm{F}},\boldsymbol{x}_{k}^{\mathrm{F}}))$ , where factual noise is implicitly encoded. The complete formalization appears in Supplementary LABEL:ACQ ( Algorithm LABEL:algo:_PFD-BDCM_counterfactual_inference).

4 Theoretical Guarantees

This section develops the theoretical underpinnings for integrating functional data into structural causal models via basis expansion, establishing rigorous guarantees for the counterfactual estimation accuracy of the proposed PFD-BDCM framework.

Notation: We consider an endogenous variable $\textnormal{x}_{k}$ ( or $\mathbf{x}^{(z)}$ for functional nodes ) governed by structural equation $\textnormal{x}_{k}:=f_{k}(\mathbf{x}_{\mathbb{PB}_{k}},\textnormal{u}_{k})$ ( or $\mathbf{x}_{k}^{(z)}=f_{k}^{(z)}(\mathbf{x}_{\mathbb{PB}_{k}^{(z)}},\mathbf{u}_{k}^{(z)})$) with parent set and backdoor adjustment set $\mathbf{x}_{\mathbb{PB}_{k}}$ ( or $\mathbf{x}_{\mathbb{PB}_{k}^{(z)}}$) and exogenous noise $\textnormal{u}_{k}$ ( or $\mathbf{u}_{k}^{(z)}$). We analyze a single node without loss of generality (by permutation invariance of nodes), henceforth denoting the target variable as $\textnormal{x}\in\mathcal{X}\subseteq\mathbb{R}$ ( or the target vector $\mathbf{x}^{(z)}\in\mathcal{X}^{K_{n}}\subseteq\mathbb{R}^{K_{n}}$), its parents and backdoor adjustments as $\mathbf{x}_{\mathbb{PB}}\in\mathcal{X}_{\mathbb{PB}}\subseteq\mathbb{R}^{K}$ ( or $\mathbf{x}_{\mathbb{PB}^{(z)}}\in\mathcal{X}_{\mathbb{PB}^{(z)}}\subseteq\mathbb{R}^{K}$), and exogenous noise as u ( or $\mathbf{u}^{(z)}$). The encoder-decoder architecture comprises

[TABLE]

where $\mathcal{U}$ ( or $\mathcal{U}^{K_{n}}$) denotes the (${K_{n}}$ dimensional ) latent space. Within PFD-BDCM, $g$ ( or $g^{(z)}$) and $h$ ( or $h^{(z)}$) correspond to the $\mathrm{Enc}(\cdot,\cdot)$ and $\mathrm{Dec}(\cdot,\cdot)$ operators, respectively.

4.1 Approximate Error Bound of Functional Data

First, a fundamental challenge is ensuring that the finite-dimensional basis approximation of $z(t)$ preserves its true causal effect. Theorem 1 formally addresses this concern. Under standard smoothness assumption (Assumption 3) , we establish that the interventional distribution of any outcome under $\operatorname{do}(z(t)=z_{0}(t))$ converges weakly to the distribution under $\operatorname{do}(\mathbf{x}^{(z)}=\mathbf{x}^{(z)}_{0})$ as the basis dimension $K_{n}$ increases. Detailed proofs are provided in the Supplementary LABEL:Theorem-proof-func.

Assumption 3**.**

The coefficient function $c(t)$ and the functional random variable $\textnormal{z}(t)$ belong to a Hilbert space $\mathcal{H}\subseteq L^{2}(\mathcal{T})$ with inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}}$ and norm $\|\cdot\|_{\mathcal{H}}$, and the orthonormal basis $\{b_{m}(t)\}_{m=1}^{\infty}$ is complete in $\mathcal{H}$ . Moreover, there exist constants $s_{c},C_{c},s_{z},C_{z}>0$ such that the approximation error for any $K_{n}\in\mathbb{N}$,

[TABLE]

where $\phi_{m}=\langle c(t)),b_{m}(t)\rangle_{\mathcal{H}}$ and $\textnormal{x}_{m}=\langle\textnormal{z}(t),b_{m}(t)\rangle_{\mathcal{H}}$ are the basis coefficients.

Assumption 3 is standard in functional data analysis and nonparametric estimation. First, the requirement that the function space $\mathcal{H}$ possesses a complete orthonormal basis guarantees that any element, here the coefficient function $c(t)$ and the realizations of the functional covariate $z(t)$, can be represented as a linear combination of the basis functions, thereby justifying the use of projection‑based or sieve estimation approaches. Second, the specific power‑law decay rates $K_{n}^{-s_{c}}$ and $K_{n}^{-s_{z}}$ for the truncation errors quantify the smoothness of the underlying functions; functions with higher smoothness (i.e., larger $s_{c}$, $s_{z}$) can be approximated more accurately with a finite number of basis terms. These rates are crucial for establishing the theoretical convergence rates of the resulting estimators. Common examples satisfying such conditions include functions in Sobolev or Hölder spaces when Fourier, Spline, or Wavelet bases are used. The moment condition on $z(t)$ ensures that the approximation error holds in a mean‑squared sense across its distribution.

Theorem 1**.**

Under Assumption 3, the causal effect of the functional variable $\textnormal{z}(t)$ on any endogenous variable x in the infinite-dimensional SCM is equivalent to the causal effect of the finite-dimensional basis coefficient vector $\mathbf{x}^{(z)}$ on x in the PFST-DSCM, as $K_{n}\to\infty$. Specifically, for any intervention $\text{do}(\textnormal{z}(t)=z_{0}(t))$, the resulting interventional distribution on x converges weakly to the interventional distribution obtained by the intervention $\text{do}(\mathbf{x}^{(z)}=\boldsymbol{x}_{0}^{(z)})$ in the finite-dimensional model, where $\boldsymbol{x}_{0}^{(z)}$ is the projection of $z_{0}(t)$ onto the basis.

Theorem 1 implies that all statistical functionals (including means, variances, and quantiles) that are continuous with respect to the weak topology will converge. In particular, the causal effect of $\textnormal{z}(t)$ on x, defined as a functional of the interventional distribution (e.g., the average treatment effect), is preserved in the limit. Thus, the finite-dimensional intervention $\text{do}(\mathbf{x}^{(z)}=\boldsymbol{x}_{0}^{(z)})$ yields an interventional distribution that converges weakly to that of the infinite-dimensional intervention $\text{do}(\textnormal{z}(t)=z_{0}(t))$, and the causal effect of the functional variable $\textnormal{z}(t)$ is equivalently captured by the basis coefficient representation as the number of basis functions increases.

Next, Theorem 2 establishes an upper bound on the approximation error induced by the finite-dimensional projection, which decays with the number of basis functions. Detailed proofs are provided in Supplementary LABEL:Theorem-proof-func.

Theorem 2**.**

Under Assumption 3, assume further that the intervention function $z_{0}(t)$ satisfies a deterministic approximation error bound: there exist constants $C_{z_{0}},s_{z_{0}}>0$ such that for any $K_{n}\in\mathbb{N}$,

[TABLE]

where $x_{0,m}=\langle z_{0}(t),b_{m}(t)\rangle_{\mathcal{H}}$ (assuming orthonormal basis). Let $B_{c}=\|c(t)\|_{\mathcal{H}}$ and $B_{z_{0}}=\|z_{0}(t)\|_{\mathcal{H}}$ be finite. Then, for any structural function $f$ that is Lipschitz continuous with respect to the functional term with Lipschitz constant $L_{f}$, and for any intervention $\text{do}(\textnormal{z}(t)=z_{0}(t))$, the pointwise error satisfies

[TABLE]

Consequently, the error is of order $O(K_{n}^{-\min(s_{c},s_{z_{0}})})$. In particular, the error in the interventional mean of x (i.e., the average causal effect) decays at the same polynomial rate.

Theorem 2 establishes the convergence rate for the approximation error incurred when representing infinite-dimensional functional quantities in causal estimands with a finite-dimensional basis expansion. The total pointwise error decomposes into terms arising from the approximation of the functional coefficient $c(t)$ and the intervention function $z_{0}(t)$, each decaying at a polynomial rate dictated by their respective smoothness parameters $s_{c}$ and $s_{z_{0}}$. Consequently, the overall error is of order $O(K_{n}^{-\min(s_{c},s_{z_{0}})})$, demonstrating that smoother underlying functions permit accurate estimation with fewer basis functions. This result directly implies the same convergence rate for the average causal effect, providing a theoretical foundation for sieve estimation in functional causal inference by formally characterizing the trade-off between model complexity (governed by $K_{n}$) and approximation accuracy.

4.2 Counterfactual Error Bounds

The accuracy of any counterfactual query hinges on a model’s ability to correctly perform abduction, that is, to recover the exogenous noise that generated the observed data. This step is what differentiates causal models from purely predictive ones. Our theoretical analysis establishes that under our framework, this recovery is possible, and its fidelity directly controls the final counterfactual error. To comprehensively validate PFD-BDCM, we analyze three progressively complex settings: scalar nodes (the baseline case), functional (multidimensional) nodes, and finally, scenarios with heavy-tailed noise or model misspecification. This structure mirrors the increasing complexity of real-world data and allows us to pinpoint the assumptions required at each step.

4.2.1 One-Dimensional Counterfactual Error Bounds

We first establish the core error bound for a scalar endogenous variable $\textnormal{x}\in\mathbb{R}$ ( detailed proofs show in Supplementary LABEL:sup:Error_Bounds-LD). This forms the theoretical backbone, which we later extend. Our theoretical results rely on a set of assumptions regarding the structural equation and the encoder-decoder model. These conditions are essential for ensuring that the latent variable learned by the encoder can uniquely recover the unobserved exogenous noise, which is the cornerstone of accurate counterfactual estimation (Lu et al.,, 2020; Nasr-Esfahany and Kiciman,, 2023). For a variable $\textnormal{x}\in\mathcal{X}\subset\mathbb{R}$ with structural equation $\textnormal{x}:=f(\mathbf{x}_{\mathbb{PB}},\textnormal{u})$ where $\textnormal{u}\perp\!\!\!\perp\mathbf{x}_{\mathbb{PB}}$, we have the following assumptions:

Assumption 4**.**

The encoding function produces representations that are statistically independent of the parent variables and the backdoor adjustment variables: $g(\mathbf{x}_{\mathbb{PB}},\textnormal{x})\perp\!\!\!\perp\ \mathbf{x}_{\mathbb{PB}}$.

Assumption 5**.**

The encoding function is invertible and differentiable with respect to the endogenous variable for fixed parent variables and backdoor adjustment variables: $g(\mathbf{x}_{\mathbb{PB}},\cdot):\mathcal{X}\to\mathcal{U}\text{ is bijective and }C^{1},\quad\text{for all }\mathbf{x}_{\mathbb{PB}}\in\mathcal{X}_{\mathbb{PB}}.$

Assumption 6**.**

The structural equation is continuously differentiable in both arguments and strictly monotonic in the exogenous noise for each fixed value of the parent variables and the backdoor adjustment variables $f\in C^{1}(\mathcal{X}_{\mathbb{PB}}\times\mathcal{U},\mathcal{X}),\quad\frac{\partial f(\mathbf{x}_{\mathbb{PB}},\textnormal{u})}{\partial\textnormal{u}}>0\quad\text{for all }\mathbf{x}_{\mathbb{PB}}\in\mathcal{X}_{\mathbb{PB}},\textnormal{u}\in\mathcal{U}.$

These assumptions, while formal, are well-motivated in the context of causal inference and deep generative models. Assumption 4 ensures that the encoder learns a pure representation of the exogenous noise uncontaminated by information from the backdoor adjustment variables. In practice, this can be enforced through regularization or architectural constraints. This is naturally satisfied in settings like additive noise models with $f(\mathbf{x}_{\mathbb{PB}},\textnormal{u})=f^{*}(\mathbf{x}_{\mathbb{PB}})+\textnormal{u}$ where $\mathbf{x}_{\mathbb{PB}}$ and u is independent. If the fitted model $\hat{f}\equiv f^{*}$, then $g(\mathbf{x}_{\mathbb{PB}},\textnormal{x})=\textnormal{u}$. Assumption 5 is intrinsically satisfied by the bijective properties of deterministic diffusion architectures (Song et al.,, 2021), guaranteeing uniqueness in latent representations while preserving compatibility with standard implementations. Assumption 6 is satisfied by major identifiable model classes, including additive noise, post-nonlinear, and heteroscedastic formulations (Strobl and Lasko,, 2023) while concurrently resolving symmetric noise ambiguities characteristic of observational data. This assumption further aligns with contemporary identifiability frameworks (Nasr-Esfahany and Kiciman,, 2023) and intrinsically precludes non-identifiable structural equations. Additional justifications of this assumptions’ reasonableness in practical settings are provided in Supplementary LABEL:sup:Error_Bounds-LD.

These assumptions collectively ensure that the encoded latent $\hat{\textnormal{u}}=g(\mathbf{x}_{\mathbb{PB}},\textnormal{x})$ is an invertible transform of the true noise u (Theorem 3).

Theorem 3**.**

Under Assumptions 4, 5, and 6, the encoded latent variable $\hat{\textnormal{u}}=g(\mathbf{x}_{\mathbb{PB}},\textnormal{x})$ is related to the true exogenous noise u through an invertible transformation. That is, there exists an invertible function $\phi:\mathcal{U}\to\mathcal{U}$ such that $\hat{\textnormal{u}}=\phi(\textnormal{u})$.

Theorem 3 implies that the abduction step $g(\mathbf{x}_{\mathbb{PB}},\textnormal{x})$ correctly captures the essence of the unobserved exogenous random variable u that generated the factual observation. We now explore the direct consequences of this result. In an oracle scenario where the model achieves perfect reconstruction which means $h(\mathbf{x}_{\mathbb{PB}},g(\mathbf{x}_{\mathbb{PB}},\textnormal{x}))=\textnormal{x}$ holding almost surely (a.s.), the counterfactual estimate will be “perfect”. This precise case implies the satisfaction of Theorem 3, yielding the relationship $h(\mathbf{x}_{\mathbb{PB}},\phi(\textnormal{u}))=f(\mathbf{x}_{\mathbb{PB}},\textnormal{u})$. Consequently, when making a counterfactual prediction for a new intervention $\mathbf{x}_{\mathcal{L}}:=\boldsymbol{\gamma}$, the model computes $h(\mathbf{x}_{\mathbb{PB,-\mathcal{L}}}\cup\boldsymbol{\gamma},\phi(\textnormal{u}))$, which a.s. equates to the true counterfactual $f(\mathbf{x}_{\mathbb{PB,-\mathcal{L}}}\cup\boldsymbol{\gamma},u)$.

In empirical settings, exact reconstruction is rarely achievable. From this foundation, we derive our central result (Theorem 4): for any factual observation and counterfactual intervention, the estimation error is bounded by the model’s empirical reconstruction error.

Theorem 4**.**

Let $\mathcal{X}$ be a metric space of functions with metric $d:\mathcal{X}\times\mathcal{X}\to\mathbb{R}_{\geq 0}$ that induces the topological structure of interest (e.g., the $L^{2}$ norm for square-integrable functions). Under Assumptions 4, 5, and 6, for any factual observation $(x^{F},\mathbf{x}_{\mathbb{PB}}^{F})$ generated by $x^{F}=f(\mathbf{x}_{\mathbb{PB}}^{F},u)$ and any counterfactual intervention $\mathrm{do}(\mathbf{x}_{\mathcal{L}}^{CF}:=\boldsymbol{\gamma})$, if the reconstruction error is bounded by $\tau$ with respect to $d$, i.e., $d(h(\mathbf{x}_{\mathbb{PB}},g(\mathbf{x}_{\mathbb{PB}},\textnormal{x})),\textnormal{x})\leq\tau$, then the counterfactual estimation error is bounded by $\tau$: $d(\hat{x}^{\mathrm{CF}},x^{\mathrm{CF}})\leq\tau.$

This bound (Theorem 4) is powerful and elegant because it directly connects a standard, optimizable training objective (reconstruction loss) to the causal fidelity of the model. This provides a clear theoretical justification for using a powerful generative model (diffusion) for causal inference: improving its generative performance guarantees better counterfactual estimates. This link is not explicitly established in prior diffusion-based causal works like BDCM, which focus more on algorithmic feasibility than on such performance-based guarantees.

4.2.2 Multidimensional Counterfactual Error Bound

Handling functional data constitutes a methodological component of PFD-BDCM. Theoretically, this entails generalizing the scalar assumptions to the multidimensional vector $\mathbf{x}^{(z)}\in\mathbb{R}^{K_{n}}$ of basis coefficients. Specifically, Assumptions 4–6 are extended by replacing scalar monotonic functions with diffeomorphisms (invertible smooth maps), as detailed in Assumptions LABEL:ass:multidim_struct–LABEL:ass:multidim_indep in Supplementary LABEL:sup:Error_Bound-HD-proof. This generalization preserves the principle of “noise recovery” in the higher-dimensional setting.

For counterfactual queries that involve both functional and scalar nodes, additional regularity of the causal paths is required. We distinguish three canonical cases, each requiring a specific cross‑dimensional Lipschitz condition. The details show in Supplementary LABEL:sup:Error_Bound-HD-proof (Assumption LABEL:ass:case1-LABEL:ass:case3). Under these assumptions we obtain the following multidimensional counterpart of Theorem 5, and all detailed proofs show in Supplementary LABEL:sup:Error_Bound-HD-proof.

Theorem 5**.**

Assume the relevant node‑type assumptions (Assumptions LABEL:ass:multidim_struct–LABEL:ass:multidim_indep for functional nodes, and Assumptions 4–6 for scalar nodes) hold. Then the counterfactual estimation error satisfies the following bounds for the three canonical cases:

Case 1 (functional intervention, scalar outcome): Under Assumption LABEL:ass:case1, $d(\hat{x}^{\mathrm{CF}},x^{\mathrm{CF}})\leq\tau+L_{\text{cross}}\bigl(\tau^{(z)}+O(K_{n}^{-s})\bigr)$.

Case 2 (functional intervention, functional outcome): Under Assumption LABEL:ass:case2, $d(\hat{\mathbf{x}}^{(z^{\prime})^{\mathrm{CF}}},\mathbf{x}^{(z^{\prime})^{\mathrm{CF}}})\leq\tau^{(z^{\prime})}+O(K_{n}^{-s})+O(K_{n}^{\prime-s^{\prime}})$.

Case 3 (scalar intervention, functional outcome): Under Assumption LABEL:ass:case3, $d(\hat{\mathbf{x}}^{(z)^{\mathrm{CF}}},\mathbf{x}^{(z)^{\mathrm{CF}}})\leq\tau^{(z)}+L_{\text{cross}}^{(z)}\tau_{l}+O(K_{n}^{-s})$.

Here $\tau$, $\tau^{(z)}$, $\tau^{(z^{\prime})}$, $\tau_{l}$ denote the reconstruction errors of the corresponding nodes, measured in the appropriate norms (e.g., Euclidean norm for vectors, absolute value for scalars). The terms $O(K_{n}^{-s})$ and $O(K_{n}^{\prime-s^{\prime}})$ are the basis‑expansion errors for the functional variables involved.

The multidimensional treatment has several practical advantages: i) It is computationally more efficient to train a single diffusion model for the whole vector of basis coefficients than separate models for each coefficient; ii) The joint model captures dependencies among the coefficients, leading to statistically more consistent estimates; iii) The approach scales gracefully with the number of basis functions provided the Lipschitz constants remain bounded. In current research on causal diffusion models, no existing work has systematically addressed how to incorporate functional data into the causal diffusion modeling framework with accompanying theoretical guarantees. Our multidimensional framework fills this gap. It is the first to enable the modeling of functional data within a causal diffusion model and provide a corresponding error analysis, offering a novel theoretical tool for handling multi-resolution spatio-temporal data.

4.2.3 Robustness to Heavy-Tailed Noise and Model Misspecification

To demonstrate the practical utility and robustness of our framework, we extend our theoretical analysis beyond idealized settings, considering realistic challenges such as heavy-tailed noise distributions and potential model misspecification. Detailed formulations of the corresponding assumptions, theorems, and complete proofs for these extended scenarios are provided in the Supplementary LABEL:sec:HTN and LABEL:sec:lim.

Heavy-Tailed Noise. Real-world environmental data often exhibit outliers. Under a polynomial tail assumption (Assumption LABEL:Heavy-Tailed_Exogenous_Noise in Supplementary LABEL:sec:HTN ) on the noise u, we show (Theorem LABEL:theo:Moment-Based_Counterfactual_Bound in Supplementary LABEL:sec:HTN ) that the counterfactual error inherits a similar polynomial concentration bound. More importantly, we prove (Theorem LABEL:theo:Robust_Encoder_Design in Supplementary LABEL:sec:HTN ) that using a robust loss function (e.g., Huber loss) for the encoder can recover an exponential concentration rate, significantly tightening the error control.

Remark 1*.*

Theorem LABEL:theo:Moment-Based_Counterfactual_Bound provides a polynomial concentration bound for counterfactual estimation errors when dealing with heavy-tailed noise distributions. This result demonstrates that the tail of the counterfactual error is controlled by the tail of the reconstruction error. Hence, by controlling the tail of the reconstruction error, we control the tail of the counterfactual error. Specifically, the tail behavior is governed by the worse of these two moment conditions, with constants $C_{\tau}$ depending on the Lipschitz properties and model parameters. This polynomial decay reflects the fundamental challenges posed by heavy-tailed noise in causal inference.

Remark 2*.*

Theorem LABEL:theo:Robust_Encoder_Design shows that employing Huber loss in the encoder design achieves exponential concentration despite heavy-tailed noise. The robust M-estimation approach ensures that the encoder error probability decays exponentially fast as $\epsilon$ increases, providing significantly tighter control over estimation errors compared to the polynomial bound. This demonstrates the substantial benefits of robust estimation techniques in handling heavy-tailed distributions.

Model Misspecification. The theoretical guarantees established so far rely on the assumption that the structural equation follows an additive noise model (Assumption LABEL:ass:_Additive_Noise_Model) or, more generally, a strictly monotonic form (Assumption 6). In practice, the true data-generating process may deviate from these idealized forms. To address this, we now analyze the robustness of the PFD-BDCM framework when the structural equation is only approximately additive or when the encoder-decoder pair does not perfectly satisfy the invertibility and independence conditions. We discuss the limitations of our model assumptions (Supplementary LABEL:sec:Assumption_Violations_and_Robustness) and provide error bounds under model misspecification, which are consolidated in Proposition LABEL:prop:general-robustness and its detailed proof all show in Supplementary LABEL:sec:Model_Misspecification.

Remark 3*.*

Proposition LABEL:prop:general-robustness extends the error analysis to account for model misspecification. It establishes that the counterfactual error of PFD-BDCM remains bounded even when the true structural equations deviate from the assumed form, provided the deviation is controlled. The bound explicitly decomposes the total error into components arising from encoder inaccuracy ($\delta_{n}$), decoder approximation error ($\tau_{n}$), and the degree of misspecification ($\zeta$). Under heavy-tailed noise, employing a robust encoder yields near-exponential error concentration, with performance primarily limited by the misspecification term $\eta$. This result demonstrates that the framework maintains controlled error under model mismatch and remains applicable to real-world problems, such as the air pollution case study, where ideal assumptions are often violated.

Our theoretical analysis provides a cohesive and multi-layered foundation for PFD-BDCM. These results do more than justify our specific algorithm; they advance the theoretical understanding of what makes deep generative causal inference work, particularly when dealing with the spatio-temporal and functional complexities that define modern environmental datasets.

5 Simulation Studies

Our synthetic study served as a controlled setting to evaluate the contribution of each methodological innovation within PFD-BDCM. We assessed the performance of a total of 11 models, including our proposed PFD-BDCM, and a range of baseline and alternative approaches, across eight distinct structural equation scenarios.

5.1 Data Generation

We considered a PFST-DSCM model comprising 11 exogenous and endogenous nodes, as illustrated in Fig. 2. In this model, $\{\textnormal{x}_{3,ij},\textnormal{x}_{4,ij}\}$ were set as unobserved explanatory variables, and $\textnormal{x}_{9,ij}$ was set as an unobserved response variable. We assumed that $\{\textnormal{x}_{3,ij},\textnormal{x}_{4,ij}\}$ and $\textnormal{x}_{9,ij}$ exhibited pronounced spatial heterogeneity along with temporal dependence. Let $\textnormal{x}_{1,ij}$ represent the endogenous cause variable; $\mathbf{x}_{2,ij}$ represent the functional cause variable, $\{\textnormal{x}_{10,ij},\textnormal{x}_{11,ij}\}$ denote the outcome variables, and $\mathbf{x}_{\mathbb{B}}=\{\textnormal{x}_{5,ij},\textnormal{x}_{6,ij},\textnormal{x}_{7,ij},\textnormal{x}_{8,ij}\}$ constitute the backdoor adjustment sets.

To systematically stress-test different methodological aspects, we generated data from eight distinct Data-Generating Processes (DGPs), creating a full-factorial design varying: i) Linearity (Linear vs. Nonlinear structural equations), ii) Noise Interaction (Additive vs. Non-additive) and iii) Noise Distribution (Gaussian vs. Heavy-tailed Student’s $t$). This allowed us to assess robustness across settings that were simple, challenging, and realistic. The specific explanations of the model names were shown in Table LABEL:tab:model_acronyms (Supplementary LABEL:sec:Experiment_details).

5.2 Simulation Settings

Baseline Models. We benchmarked PFD-BDCM against a carefully chosen hierarchy of 10 alternative models, each designed to isolate the impact of a specific feature. First, the “PFD-Oracle” model served as an oracle, possessing access to the true structural equations. It represented the theoretical performance ceiling and provided an upper bound on achievable accuracy. The primary comparisons were against a series of ablation models derived from our own framework. “PFD-DCM” was PFD-BDCM without backdoor adjustment, isolating the value of incorporating this adjustment into the diffusion process. “PF-DCM” further ablated both backdoor adjustment and the dynamic modeling of confounders while retaining basis expansion for functional data. “PF-BDCM” included backdoor adjustment and basis expansion but treated unmeasured confounders as static, thereby isolating the specific benefit of modeling their dynamic nature. To validate our functional data integration strategy, we included several variants using mean-flattening instead of basis expansion: “PFD-BDCM-FM”, “PFD-DCM-FM”, “PF-BDCM-FM”, and “PF-DCM-FM”. These tested the impact of a simpler functional summary across different levels of our causal adjustment hierarchy. Finally, we compared against two traditional adjustment methods that condition on the same backdoor set but used non-generative models. “PFD-BLR” (Partially Functional Dynamic Backdoor Linear Regression) used linear regression to test the added value of the diffusion model’s flexibility over a simple parametric model. “PFD-BRF” (Partially Functional Dynamic Backdoor Random Forest) used a random forest to test the benefit of our generative approach against a flexible, but non-generative, nonlinear model. This comprehensive hierarchy allowed us to precisely attribute performance gains to specific components of our proposed framework.

Evaluation Protocol. For each DGP and model, we evaluated performance on the three fundamental causal queries-observational (Obs.), interventional (Int.), and counterfactual (CF.)-using the Maximum Mean Discrepancy (MMD2, Gretton et al., (2012)) between the model-generated and true distributions, lower MMD2 indicated better performance. We tested two sample sizes (J=6, n=30 and n=200) and two intervention types (uniform and quantile).

Detailed simulation specifics, including the data-generating process, parameter settings, model configurations and training hyperparameters, were provided in Supplementary LABEL:sec:Experiment_details.

5.3 Simulation Results

Complete observational, interventional, and counterfactual query results were shown in Fig. 3, LABEL:fig:MMD_all-uniform and LABEL:fig:Time_all. Detailed results were provided in Supplementary LABEL:sec:Experiment_details (Tables LABEL:sim:_result1 and LABEL:sim:_result2 summarized the aggregated performance metrics averaged over two independent random initializations). Kernel density estimates (Figs. LABEL:fig:kde1–LABEL:fig:kde8) were included in Supplementary LABEL:sec:Experiment_details.

The aggregated results reveal a clear and differentiated pattern that directly supports our core methodological claims. First, we find compelling evidence for the superiority of our approach in settings where unmeasured confounders evolve over time. Under such conditions, PFD-BDCM consistently and substantially outperforms both PF-DCM, which ignores confounding altogether, and PF-BDCM, which incorrectly assumes that confounders remain static. The performance gains are especially pronounced for interventional and counterfactual queries, which are precisely the estimands most vulnerable to unadjusted bias. For example, in nonlinear, non-additive noise settings, the MMD2 achieved by PFD-BDCM is often 20–40% lower than that of PF-BDCM. These results provide direct empirical confirmation that explicitly modeling the temporal dependence and spatial heterogeneity of confounders is not merely a theoretical nicety but a practical imperative for reliable causal inference.

Turning to the role of backdoor adjustment within the diffusion process, the comparison between PFD-BDCM and its ablated counterpart PFD-DCM is instructive. Although both models share the same architectural backbone, PFD-DCM, which conditions only on parent nodes, performs markedly worse across all data-generating processes that involve unmeasured confounders. This difference underscores the effectiveness of our central innovation: integrating the backdoor adjustment set into the conditioning mechanism of the AC-DDIM, thereby actively correcting for confounding bias during the generative sampling phase.

We also observe clear support for our use of basis expansion to handle functional data. Across multiple comparisons, PFD-BDCM reliably outperforms its mean-flattening variant, PFD-BDCM-FM. The advantage is most evident in queries involving the functional cause node $\mathbf{x}_{2}$, where the basis expansion preserves information embedded in the shape of the curve-information that is lost when a simple mean is used instead. This finding validates our decision to treat functional variables as multi-dimensional nodes represented through basis coefficients.

Finally, our results speak to the robustness and scalability of the proposed framework. PFD-BDCM maintains strong performance even under heavy-tailed noise and in highly nonlinear regimes, conditions under which traditional methods such as PFD-BLR and PFD-BRF break down entirely. Moreover, as sample size increases, the model’s performance improves and stabilizes, with MMD2 converging steadily toward the Oracle benchmark. This behavior suggests statistical consistency and reinforces the practical viability of our approach in real-world settings.

6 Empirical Application

We employed the PFD-BDCM framework to examine spatio-temporal dynamic structural causal relationships among air pollutant indicators and their determinants. Our analysis encompassed 30 provincial-level administrative divisions across Chinese mainland during the period January 2015 to December 2020. The study integrates China’s provincial CO2 emission inventories from the China Emission Accounts and Datasets (CEADs) (Xu et al.,, 2024) and emissions data for nine atmospheric pollutants from the Multi-scale Emission Inventory of China (MEIC) (Geng et al.,, 2024) as response variables for air pollutant emissions.

Building upon prior research (Ozcan,, 2013; Zhu et al.,, 2021) and incorporating domain-specific characteristics of regional emissions, we systematically collected foundational determinants across ten conceptual dimensions. The comprehensive dataset comprised 118 indicator variables, through collinearity diagnostics and random forest-based feature selection, we retained 49 statistically robust indicators for subsequent modeling (detailed indicators shown in Supplementary LABEL:real_example Table LABEL:PEI).

We constructed the PFST-DSCM for the human-driven mechanisms of atmospheric pollutant emissions (HDM-APE) in China, as illustrated in Fig. LABEL:PFST-DSCM-HDM-APE (Supplementary LABEL:real_example ). This model provided prior knowledge for systematically analyzing the impact of human activities on atmospheric emissions, wherein the structure of the causal graph and the estimation of its weight coefficients were derived through a partially functional structural equation model. Building on this, we utilize the data along with the structural relationship information to train a PFD-BDCM model specifically tailored for China’s atmospheric pollutant emission indicators and their corresponding anthropogenic driving factors. To further validate our approach, we implemented a traditional Bayesian structural equation modeling baseline (PF-DSEM-Bays) for causal query estimation. Complete experimental specifications and supplementary materials were documented in Supplementary LABEL:real_example, with observational query results presented in Table 1.

Table 1 compares the performance of PFD-DCM, PFD-BDCM, and the baseline PF-DSEM-Bays on observational queries for ten atmospheric pollutants, including CO2. Both PFD-DCM and PFD-BDCM demonstrate significantly lower MMD2 values (approximately $1.121-1.218\times 10^{-2}$) across all pollutants compared to the baseline (ranging from $11.961-36.627\times 10^{-2}$), confirming their superior capability in approximating the true observational data distribution. This indicates that both models are highly effective, with their relative strengths varying slightly across different pollutants.

The robust observational performance of the model enables intervention analysis, allowing us to simulate the impact of policy measures or environmental changes on pollutant emissions. We conducted systematic intervention experiments on 44 drivers affecting atmospheric pollutant emissions at quantile levels of $0\%,25\%,50\%,75\%$, and $100\%$ to quantitatively assess how varying intervention intensities influence emission levels. Intervention query results are provided in the supplementary materials (HDM-APE.zip). Here, we illustrate the intervention outcomes using a key social driver-“Familyscale-Small ($1-3$)” (Fig. LABEL:int1.)-as an example; the analysis process for other cases follows a similar approach.

We quantified the intervention effect using the difference between the mean residuals of the observational queries (true value $-$ observationally generated value) and the mean residuals of the interventional queries (true value $-$ interventionally generated value), defined as $\Delta$ = interventionally generated value $-$ observationally generated value, where $\Delta<0$ indicated a reduction in emissions after the intervention, $\Delta>0$ indicated an increase. The results are shown in the last subplot of Fig. LABEL:int1., with detailed interpretations in Supplementary LABEL:real_example.

7 Discussion

We have developed the Partially Functional Dynamic Backdoor Diffusion-based Causal Model (PFD-BDCM), a unified framework for causal inference with spatio-temporally dependent unmeasured confounding and mixed functional-scalar data. The methodology addresses the limitations of existing approaches that assume static confounders and lack principled functional data handling. By integrating a formal structural model (PFST-DSCM) with basis expansions and a backdoor-adjusted diffusion process, PFD-BDCM enables rigorous interventional and counterfactual queries from observational data under a known causal graph.

Our contributions are threefold. Methodologically, we synthesize functional data analysis, spatio-temporal statistics, and deep generative modeling. It explicitly models spatial heterogeneity and temporal dynamics in unmeasured confounders via CAR processes and represents functional variables through basis coefficients as standard graph nodes. Theoretically, we establish that causal effects are preserved under basis expansion with explicit error bounds and prove that counterfactual estimation error is controlled by reconstruction fidelity, with robustness to heavy-tailed noise and misspecification. Empirically, PFD-BDCM outperforms strong baselines across nonlinear, non-additive, and heavy-tailed settings in simulations, and yields actionable insights in a real-world air pollution study.

This framework can be extended in several directions. When the causal graph is only partially known, integrating structure learning algorithms could broaden applicability, though identifiability under dynamic confounding requires careful consideration. For functional data not well-approximated by a fixed basis, more flexible representations (e.g., neural networks or Gaussian processes) may be beneficial. An additional direction concerns settings where target outcomes are unavailable or populations differ substantially, necessitating causal transportability. Recent work on federated adaptive causal estimation (FACE) (Han et al.,, 2025) provides a principled framework for incorporating heterogeneous data from multiple sites while preserving privacy and communication efficiency; adapting such federated learning approaches could enable our framework to identify prevailing dynamic structures across sites when local outcomes are sparse or unobserved.

Limitations include reliance on a known causal graph, potential inefficiency of fixed-basis expansions for complex functional data, and scalability challenges for large graphs. Future work may address these through structure learning integration, adaptive basis selection, and architectural innovations for amortized inference. Despite these limitations, PFD-BDCM provides a rigorous foundation for causal inference with complex spatio-temporal and functional data, with broad applicability in environmental science, epidemiology, and beyond.

Conflict of interest/Competing interests: The authors declare no competing interests. ** Acknowledgments:** The authors would like to thank the referees, the Associate Editor, and the Editor for their constructive comments that improved the quality of this paper.