Network Recovery from Cascade Data: A Debiased Jacobian-Based Machine Learning Approach* Lei Huang
arXiv:2606.07483v1 [cs.LG] 5 Jun 2026
June 1, 2026
Abstract Many important outcomes unfold as dynamic cascades, including product adoption, disease spread, financial distress, and information diffusion. A central challenge is to recover the hidden influence network behind these cascades. Existing methods typically assume a specific diffusion model, and their performance degrades substantially when that assumption is misspecified. We propose CascadeNet, a Jacobianbased machine learning framework for network recovery that does not require specifying a diffusion mechanism. The key idea is that the underlying influence structure can be characterized by the Jacobian of the one-step transition function. CascadeNet first constructs a flexible estimator of the transition function, and further applies Neyman-orthogonal debiasing via the Riesz representer, so that the debiased √ Jacobian is n-consistent and asymptotically normal, enabling formal inference on the network structure. We validate CascadeNet in both a simulation exercise and a real-world empirical application. In simulations, where the data-generating process is known, CascadeNet achieves the highest network recovery accuracy across nine common data-generating processes. In an empirical application to COVID-19 transmission across Spain’s 52 provinces, CascadeNet recovers transmission networks that are significantly correlated with the true inter-province mobility network, whereas networks recovered by baseline methods show no significant alignment with the ground truth. Keywords: network inference, Jacobian, debiased machine learning.
* Comments welcome. Lei Huang ([email protected]) is a Ph.D. candidate at the MIT Sloan School of Management. All errors are my own.
1
1
Introduction Networks shape economic and social outcomes. Whether a firm is designing a viral
marketing campaign, a platform is studying how product reviews propagate across its user base, or a public health agency is tracing disease transmission, the effectiveness of any intervention depends on the structure of the underlying network (Kempe et al., 2003; Ballester et al., 2006; Bramoullé et al., 2009). Yet in practice, the network is rarely observed directly. Firms do not see who actually persuaded whom to buy; platforms do not see which user influenced another’s review; public health agencies do not see which infected individual transmitted the pathogen to which other individual. What decision-makers can often observe are the cascades: the sequences of adoption, engagement, or infection events that propagate across the unobserved network over time. The central task is then network recovery: inferring the hidden influence structure from these observed cascades, so that decision-makers can identify the relationships through which influence flows and target their interventions accordingly. The economic stakes of network recovery are large. Influence-maximization campaigns, lockdown decisions in epidemiology, contact tracing protocols, and systemic-risk monitoring all depend on ranking edges by their strength: a misranked edge can mean a wasted promotional budget, a missed superspreader, or a financial contagion channel that is overlooked. Yet the same data poverty that makes network recovery valuable also makes it hard. An analyst typically observes only a handful of cascades over a moderately sized population, and must infer at least an N × N influence matrix, where N is the number of nodes or players, or a higher-dimensional object when more complex interactions are assumed, from this thin record. To tackle this challenge, a substantial literature has developed methods that assume a very specific functional form for the underlying transition process: for instance, the exponential and power-law transmission kernels in Gomez-Rodriguez et al. (2010) and Gomez-Rodriguez et al. (2011), the pairwise Markov transition probabilities in Ramezani 2
et al. (2024), and the particular GNN architectures in Qiu et al. (2018) and Murphy et al. (2021). These functional-form assumptions provide identification power that the data alone cannot, but when the assumed transition process is misspecified, the resulting error can be large. These methods also share a second drawback: they do not provide formal statistical inference on the recovered edges, which makes it hard for practitioners to communicate uncertainty about the recovered structure and to make informed decisions based on it. We propose CascadeNet, a Jacobian-based machine learning framework for network recovery that addresses both limitations. The key insight of CascadeNet is that, in a Markovian system, the influence network can be characterized by the Jacobian of the one-step transition function. This insight allows us to formulate network recovery as a problem of estimating the transition function and then differentiating it. By replacing parametric likelihood with a flexible learner of the one-step transition function, CascadeNet avoids the misspecification risk that plagues classical methods and enables application to a wider range of diffusion mechanisms, outcome types, and covariate structures. We further apply Neyman-orthogonal debiasing based on the Riesz representer, a technique from the double machine learning literature, to address bias from regularization in the flexible machine learning estimator. This allows us to recover an unbiased estimate of the network Jacobian using only observable predic√ tion residuals, and to prove n-consistency and asymptotic normality for each debiased Jacobian entry, enabling formal inference on individual edge weights. To lay out the problem, consider a panel of N agents, the nodes of the network, observed over time. The state vector Yt records each agent’s outcome at time t, such as adoption status, infection status, or case counts, and evolves according to a one-step transition function, Yt+1 = m0 (Yt , Xt )+εt+1 , where Xt collects observed covariates such as treatment assignments, demographics, or policy variables, and m0 is left fully unspecified. This formulation nests classic diffusion models, including independent cascade, linear threshold, SIS, SIR, and Hawkes, as special cases, and accommodates binary, count, or continuous
3
outcomes. CascadeNet has two components. The first is a flexible estimator m̂ for the one-step transition function m0 , and the same estimator is used iteratively across time steps. The simplest specification is a linear-index model with a learnable interaction matrix J and a flexible link function (sigmoid for binary adoption, identity or MLP for continuous outcomes), which already strictly nests classical baselines; richer specifications such as graph neural networks with attention-based aggregation can be plugged in when the application demands them. When the model class is large enough, the misspecification risk is minimal. The second component is the central observation that turns the estimated transition function into a network estimator: the Jacobian of m0 , the N × N matrix J0 with entries J0 [i, j] = ∂mi /∂yj , directly encodes the influence network. In a Markovian system, the (i, j) entry is nonzero if and only if agent j can marginally influence agent i, and its magnitude measures the strength of that influence; network recovery therefore reduces to differentiating the estimated transition function, a single operation that works regardless of which parametric diffusion model, if any, generated the data. However, this flexibility of the transition function estimator comes at a cost. Common machine learning estimators are regularized to control variance, and regularization introduces bias. For instance, ℓ1 and ℓ2 penalties shrink the estimated parameters toward zero, and therefore bias the estimated Jacobian; in a small-sample setting where regularization is most needed, this attenuation can be severe enough to destroy the network signal entirely. Off-the-shelf machine learning alone is therefore not enough. To address this bias, we apply Neyman-orthogonal debiasing based on the Riesz representer (Chernozhukov et al., 2022; Hirshberg and Wager, 2021), a technique from the double machine learning literature (Chernozhukov et al., 2018). Intuitively, the Riesz representer acts as an integration-by-parts operator: it converts errors in the prediction m̂ into corrections for the derivative ∂ m̂i /∂yj , allowing us to recover an unbiased estimate of the network Jacobian using only the observable prediction residual. The resulting debiased estimator
4
is
√ n-consistent and asymptotically normal for each individual edge, even though m0
itself is estimated by flexible machine learning. We illustrate the necessity and effectiveness of this debiasing strategy in a controlled tanh experiment with a closed-form ground truth: the naive plug-in Jacobian achieves a Pearson correlation of only r = 0.07 with the true Jacobian, while the Riesz correction restores it to r = 0.77, more than a tenfold improvement over the naive plug-in. We validate CascadeNet on both synthetic and real data. On synthetic cascades generated from nine common diffusion models from the literature, namely pairwise transmission (IC, LT), continuous-time SI (Exp, PL), epidemic models with recovery/removal (SIS, SIR), aggregate-effect models (Complex, Hawkes), and a flexible nonlinear DGP, CascadeNet achieves the highest recovery accuracy in all nine settings, outperforming established methods by up to 50 percentage points of Pearson correlation on diffusion models that violate the baselines’ kernel assumptions. On real COVID-19 case data from Spain’s 52 provinces, CascadeNet is the only method whose estimated edges correlate significantly with true inter-province mobility flows; every classical baseline produces rankings indistinguishable from random. CascadeNet makes two methodological contributions. First, it provides a flexible estimator for the influence network: by replacing parametric likelihood with a flexible learner of the one-step transition function, CascadeNet avoids the misspecification risk that plagues classical methods and enables application to a wider range of diffusion mechanisms, outcome types, and covariate structures. Second, it provides formal statistical inference: we extend the Riesz representer framework to an N × N matrix of partial derivatives, a setting not previously addressed in the debiased ML literature, and prove √ n-consistency and asymptotic normality for each debiased Jacobian entry. This delivers what classical cascade-inference methods cannot: confidence intervals for individual edge weights, hypothesis tests for whether a particular edge exists, and a principled basis for communicating uncertainty about the recovered network.
5
Managerially, CascadeNet equips decision-makers with the ability to map influence relationships from observed cascade data alone, with formal confidence intervals on individual edges, without requiring the analyst to specify how influence propagates. In practical settings such as a marketing manager planning a seeding campaign, a publichealth authority allocating contact-tracing resources, or a regulator monitoring contagion in a financial network, this combination of flexibility and inferential rigor is what allows the recovered network to be used as a serious input to consequential decisions, rather than as a heuristic ranking with unknown reliability. The remainder of the paper is organized as follows. Section 2 reviews the related literature. Section 3 presents the CascadeNet framework, including the data-generating process, the estimator, the Jacobian representation of the network, the debiasing procedure, and the asymptotic theory. Section 4.1 reports synthetic validation. Section 4.2 applies the method to COVID-19 data. Section 5 discusses managerial implications and concludes.
2
Related Literature Our work contributes to several streams of research. We discuss each in turn and
highlight how CascadeNet relates to and extends prior work. Cascade network inference.
The problem of inferring networks from cascade data has
attracted considerable attention in the network science literature. The foundational contributions are NetInf (Gomez-Rodriguez et al., 2010) and NetRate (Gomez-Rodriguez et al., 2011), which infer diffusion networks by maximizing pairwise cascade likelihoods under exponential or power-law transmission kernels and rely on submodular greedy selection or convex optimization for tractability. CONNIE (Myers and Leskovec, 2010) casts the inference problem as a convex program with a regularized likelihood, providing the first scalable formulation with provable optimality. More recently, DANI (Ramezani et al., 6
2024) takes a Markov-transition approach that estimates pairwise transmission probabilities directly from the observed panel and emphasizes preserving topological structure. A parallel line of work uses deep learning to model cascade dynamics: DeepInf (Qiu et al., 2018) predicts user-level adoption with attention-based GNNs, and Murphy et al. (2021) learn contagion dynamics on complex networks with graph neural networks, demonstrating strong predictive performance on epidemiological data. Despite their methodological diversity, these methods share two features that motivate our work. First, each imposes a specific parametric assumption on how influence propagates (an exponential hazard, a power-law kernel, a linear Markov transition, or a particular GNN architecture), and performance is sensitive to whether that assumption matches the true diffusion process. Second, none provide formal statistical inference on individual edges: the output is a ranked list or a point estimate, but not a confidence interval or a hypothesis test. Our contribution is to develop a flexible estimator that avoids the parametric kernel assumption, accommodates heterogeneous diffusion mechanisms within a single framework, and pro√ vides n-consistent debiasing guarantees that enable formal edge-level inference.
Networks in marketing.
A rich marketing literature studies how networks shape con-
sumer behavior, peer effects, and the diffusion of new products. Kempe et al. (2003) formalize the influence-maximization problem of selecting seed nodes to maximize cascade reach, which has become a central decision problem in viral marketing and presupposes knowledge of the influence network. Empirical work on word-of-mouth and social contagion in marketing (Ballester et al., 2006; Chu and Manchanda, 2016; Eckles et al., 2016) similarly takes the network as known or proxies it with observable ties (e.g., friendships, co-purchase records, geographic proximity). When the network is not directly observable, marketing researchers have used model-based approaches to infer peer effects from panel data, but these typically assume linear-in-means specifications (Manski, 1993; Bramoullé et al., 2009) that may not capture the nonlinear dynamics of cascade adoption. Our work
7
complements this literature in three ways. First, we provide a method to infer the influence network directly from observed cascades, without requiring auxiliary network data, thereby expanding the set of settings in which network-based marketing analysis is feasible. Second, by supplying confidence intervals on individual edges, we enable practitioners to communicate the statistical reliability of recovered seeding targets, a feature that has been missing from existing influence-maximization pipelines. Third, the flexible form of CascadeNet allows it to capture richer adoption dynamics (threshold effects, saturation, social reinforcement) than linear-in-means models, while still delivering tractable inference.
Debiased machine learning.
Our inference theory builds on the double/debiased ma-
chine learning (DML) framework of Chernozhukov et al. (2018) and the automatic-debiasing approach via the Riesz representer developed by Chernozhukov et al. (2022) and Hirshberg and Wager (2021). The general principle is that for a smooth functional of a nuisance function estimated by machine learning, an orthogonal score correction restores the para√ metric n rate of convergence even when the nuisance estimator converges only at n−1/4 . This framework has been applied successfully to a wide range of settings, including heterogeneous treatment effects, policy learning, structural parameter estimation, and causal effects with high-dimensional controls (Farrell et al., 2021, 2020; Ye et al., 2025). Two features of our setting are new relative to existing applications. First, the target of inference is not a low-dimensional structural parameter but an N × N matrix of partial derivatives, the network Jacobian, with an entry-specific orthogonal score for each pair (i, j). Showing that the Riesz representer admits a closed form in the linear-index case, and a tractable loss-based approximation in the GNN case, is one of the technical contributions of this paper. Second, the underlying data are panel cascade trajectories rather than i.i.d. cross-sectional observations, which requires care in setting up the cross-fitting scheme and in justifying the asymptotic argument under repeated sampling of trajectories. More
8
broadly, our paper extends the DML toolkit from causal effect estimation, where it has dominated, to network recovery, where we believe it has equally large practical implications.
Network econometrics.
A substantial econometrics literature studies identification and
estimation of causal effects in network settings, typically taking the network as observed and analyzing how peer interactions, spillovers, and equilibrium effects shape outcomes. Manski (1993) formulates the reflection problem that limits identification of endogenous peer effects, and Bramoullé et al. (2009) establish identification conditions when the network is known. More recent work studies causal inference under network interference (Eckles et al., 2017; Leung, 2022) and equilibrium effects of policy interventions (Wager and Xu, 2021; Jiang et al., 2022), again under the assumption that the relevant network structure is observed by the analyst. Our paper sits upstream of this literature: we provide a method to recover the network itself from cascade data, with formal inference on individual edges, which the analyst can then plug into the standard network-econometrics √ pipeline. The combination of a flexible first stage with debiased n-inference makes the recovered network suitable as an input to downstream causal analyses without inheriting the slow convergence rates that would otherwise compromise the validity of second-stage estimates.
3
CascadeNet for Network Recovery
3.1
The Problem We observe a panel of C independent trajectories (cascades), where C can be as small
as one. In each trajectory c, a population of N agents (customers, regions, institutions)
9
evolves according to
Yc,t+1 = m0 (Yc,t , Xc,t ; θ) + εc,t+1 ,
E[εc,t+1 | Yc,t , Xc,t , θ] = 0,
(1)
where Yc,t ∈ RN is the vector of agent states at time t, Xc,t ∈ RN ×d collects agent covariates that may vary over time (e.g., demographics, promotional exposure, treatment assignments, or policy variables), θ denotes the parameters that govern the transition function, m0 is the unknown transition function, and εc,t+1 are mean-zero idiosyncratic shocks. The researcher observes the panel {Yc,t , Xc,t } for c = 1, . . . , C and t = 1, . . . , T . For instance, in a viral marketing setting, Yc,t records adoption status across customer segments and Xc,t captures time-varying promotional exposure or static demographics; each cascade corresponds to a campaign or product launch, and m0 is the adoption response function that maps the current adoption pattern and marketing variables into next-period adoption probabilities. In an epidemiological setting, Yc,t records disease incidence across regions and Xc,t captures vaccination rates, mobility restrictions, or population density; each cascade corresponds to a wave of the epidemic, and m0 is the transmission function that maps current incidence and policy conditions into next-period incidence. Several features of this formulation are worth highlighting. First, the formulation is a general one in two senses: the state Yc,t can be binary (recording which agents have adopted by step t), continuous (representing case counts, sales, or other real-valued outcomes), or mixed; the transition function m0 is left unspecified, nesting the independent cascade, linear threshold, SIS, SIR, Hawkes, and other Markovian diffusion models as special cases. Second, the parameter θ is fully flexible: in standard network models θ includes an N × N adjacency matrix, but we impose no such structure, and θ can encode any parametric or flexible specification of the transition function. Third, the key identifying assumption is Markovian dynamics, namely that Yc,t−1 affects Yc,t+1 only through Yc,t ,
10
which is standard in diffusion models; the framework extends to non-Markovian settings by augmenting the state to include lagged values (see Appendix D.4). Lastly, we treat the covariates Xc,t as exogenous in order to focus on the network-recovery task, and leave extensions that accommodate endogenous covariates via instrumental or proxy variables to future work. The goal of network recoveryis to recover the underlying network from this data. We seek to understand the mechanism of influence in the network: given any state Yt , how does a marginal change in agent j’s state affect agent i’s next-period outcome? Existing methods answer this question through the parameters of a specific diffusion model, such as the pairwise transmission probabilities of a linear Markov model or the edge weights of an exponential kernel. We hightlight that, each of these is a special case of the Jacobian of the transition function, J0 [i, j] = ∂m0,i /∂yj , which captures the marginal influence structure without imposing a specific parametric form. We therefore target the Jacobian itself, which encodes the influence network in a flexible way that is robust to misspecification of the diffusion mechanism. The asymptotic theory we develop targets the average Jacobian E[∂m0,i /∂yj ], which summarizes the overall influence of j on i across the state space; the pointwise Jacobian at any given state is also informative and can be readily computed via automatic differentiation. Analogously, the partial derivative ∂m0,i /∂xj captures the direct effect of covariates and can be estimated and debiased within the same framework.
3.2
The CascadeNet Estimator We propose CascadeNet, a machine learning framework for network recovery. The
core idea is to estimate the transition function m0 from the observed panel using a flexible machine learning estimator m̂, then differentiate the estimated function and apply a debiasing correction to recover the network Jacobian. Concretely, the procedure has three steps: (i) choose a model class M for the transition function, (ii) fit m̂ ∈ M to predict Yc,t+1 from (Yc,t , Xc,t ) across all observed time steps and cascades, and (iii) compute the 11
average Jacobian of m̂ and apply the Riesz debiasing correction described in Section 3.3. The model class M enters only through the prediction step; the inference theory we develop in Section 3.3 applies to any M that satisfies the smoothness and convergence-rate assumptions stated there. Choosing M involves two practical considerations. First, the class should respect the natural constraints that the data impose on the range of m0 : if Yc,t is binary, the transition function should map into [0, 1] and is most naturally specified through a sigmoid link; if Yc,t is a count or rate, the transition function should be non-negative and is typically specified through an identity, ReLU, or softplus link with appropriate scaling. Second, M should be flexible enough to approximate the true transition dynamics but not so complex as to overfit the available cascade data. In our experiments, a linear-index model (see more details below) with a flexible link function is sufficient to dominate every classical baseline on all nine synthetic DGPs and on the COVID-19 application; richer classes such as multilayer perceptrons or graph neural networks become useful when the data exhibit higher-order interactions, asymmetric attention patterns, or strongly nonlinear feedback that the linear index cannot capture. For concreteness, we now describe a linear-index specification of the model class M that we use as the default in our experiments, together with a generic absorbing wrapper that can be applied on top of any specification when adoption is irreversible. Richer specifications such as graph neural networks can be plugged into the same framework; we describe one such implementation in Appendix D.3.
Default specification: Linear-index model.
The simplest specification, which we refer
to as the linear-index model, takes the form
m̂i (Yt , Xt,i ) = ℓ J Yt + ki xt,i + bi ,
12
(2)
where J ∈ RN ×N is a learnable network weight matrix, ki captures the direct effect of covariates on agent i’s transition, b ∈ RN is a bias vector, and ℓ(·) is the link function (sigmoid for binary data; identity or a small MLP for continuous data). The regularizer R(m) is specified as ℓ1 or ℓ2 weight regularization on J. Despite its simplicity, this specification is already substantially more general than the parametric baselines used in the cascadeinference literature: NetInf assumes a specific exponential kernel and pairwise independence, DANI assumes a linear Markov transition without a link function, and the classical linear-threshold and independent-cascade likelihoods further restrict the form of ℓ(·) and the structure of J. The linear index in (2) relaxes all three restrictions: it allows arbitrary pairwise interaction weights, accommodates nonlinearity through the link function, and treats the network weight matrix as a single learnable object rather than a collection of separately estimated transmission probabilities. A useful by-product of this form is that the Jacobian of m̂ admits the closed-form expression ∂ m̂i /∂yj = ℓ′ (zi ) J[i, j].
Generic augmentation: Absorbing wrapper.
For diffusion models in which adoption
is irreversible, such as epidemiological settings where the infected state is absorbing or marketing settings where adopting customers remain adopters, the transition function must satisfy the constraint that once an agent adopts, it stays adopted. Any unconstrained estimator m̂raw in the model class M, whether the linear-index model in (2) or a richer i specification, can be augmented with an absorbing wrapper that enforces this constraint exactly: raw m̂abs i (Yt , Xt,i ) = Yt,i + (1 − Yt,i ) m̂i (Yt , Xt,i ).
(3)
The wrapper leaves already-adopted agents in the absorbed state with probability one and lets m̂raw govern the transition only for agents that have not yet adopted. Cascai deNet applies the wrapper when the training panel exhibits irreversible transitions; in applications where irreversibility is substantively known, the analyst can also impose the wrapper directly. 13
Once the model class M has been specified, CascadeNet estimates m0 by minimizing a regularized prediction loss over the observed panel: C T −1 1 XX L Yc,t+1 , m(Yc,t , Xc,t ) + λ R(m), m̂ ∈ arg min m∈M CT c=1 t=1
(4)
where L is binary cross-entropy (for binary data) or mean squared error (for continuous data), λ > 0 is the regularization parameter, and R(m) penalizes model complexity. In practice we minimize (4) by stochastic gradient descent (Adam), with λ selected by heldout cross-validation on a fold of trajectories disjoint from the one used for the debiased estimator. When the debiasing step is applied (Section 3.3), m̂ is fit on each cross-fitting fold separately, and the score is evaluated on the held-out fold; this avoids the standard bias from using the same data for nuisance estimation and final inference.
3.3
Network Recovery and Inference
3.3.1
Defining the target
The final and most important step of CascadeNet is to recover an influence network from the estimated transition function m̂. Because m0 maps the current state Yt into nextperiod outcomes Yt+1 , the natural local measure of influence is its Jacobian with respect to the current state: N ∂m0,i ∂m0 (y, x) = (y, x) . ∂Y ∂yj i,j=1
(5)
The entry ∂m0,i (y, x)/∂yj measures the marginal effect of agent j’s current state on agent i’s expected next-period outcome at the state-covariate configuration (y, x). This yields a model-free notion of directed influence that does not depend on whether the underlying diffusion mechanism is independent cascade, linear threshold, SIS/SIR, Hawkes, or some richer nonlinear process. In nonlinear environments, however, influence is generally state-dependent: the ef-
14
fect of j on i may be strong in some regions of the state space and weak or zero in others, because of saturation, threshold effects, or absorbing states. For that reason, we take as our primary network estimand the average Jacobian, N ∂m0 ∂m0,i J0 := E (Yt , Xt ) = E (Yt , Xt ) . ∂Y ∂yj i,j=1
(6)
Entry J0 [i, j] summarizes the average marginal influence of agent j on agent i, averaging over the distribution of states and covariates realized in the data. We target this average Jacobian because it converts a potentially heterogeneous, state-specific influence structure into a single interpretable edge-weight matrix. This is often the economically relevant object: it provides a global ranking of influence pathways, supports comparison to external benchmarks such as mobility flows or observed ties, and serves as a stable summary for downstream decision-making. Targeting the average Jacobian also brings a major statistical advantage. Unlike a pointwise derivative evaluated at a specific (y, x), the average Jacobian is a smooth population functional of the regression function m0 , and this makes it amenable to debiasing √ and n-valid inference. The pointwise Jacobian remains useful descriptively and can always be computed from m̂ by automatic differentiation, but our formal asymptotic results focus on the average Jacobian J0 . Similarly, one can define average derivative effects with respect to the covariates. If xt,ℓ denotes the ℓ-th covariate entering the transition function, the corresponding average direct-effect matrix is
∂m0,i K0 [i, ℓ] := E (Yt , Xt ) , ∂xt,ℓ
(7)
which can be estimated and debiased by the same argument developed below. We focus on J0 in the remainder of this section because recovering cross-agent influence is the central network problem; the analogous results for K0 are immediate.
15
3.3.2
Bias in the naive plug-in
A natural estimator of J0 [i, j] is the plug-in derivative E[∂ m̂i /∂yj (Yt , Xt )]. The problem is that m̂ is estimated by regularized machine learning. Regularization is necessary to control variance and avoid overfitting in high-dimensional and nonlinear settings, but it introduces bias: ℓ1 and ℓ2 penalties shrink the fitted model toward simpler functions, and this shrinkage attenuates the corresponding derivatives. In the present application, this means the naive Jacobian estimator can systematically understate edge strengths. To see 0 this, fix a pair (i, j) and define θij := E[∂m0,i /∂yj (Yt , Xt )]. We can decompose the target as
∂(m0,i − m̂i ) ∂ m̂i (Yt , Xt ) + E (Yt , Xt ) . ∂yj ∂yj | {z } | {z }
0 θij =E
plug-in (biased)
(8)
bias term
The first term is the naive plug-in estimator. The second term is the bias in the derivative, induced by estimation error in m̂i . This second term is the key difficulty: it depends on the unknown regression function m0 , and is therefore not directly observable. Debiased network recovery reduces to finding a way to represent and estimate this bias term using observable quantities.
3.3.3
Riesz representer and the orthogonal score
The central difficulty in (8) is that the bias term depends on the unknown regression function m0 , and is therefore not directly observable. Our strategy is to rewrite this derivative bias in a form that can be estimated from data. The device that makes this possible is the Riesz representer (Chernozhukov et al., 2022; Hirshberg and Wager, 2021), which converts the derivative functional into an L2 inner product and allows the bias term to be expressed using observable prediction residuals. This alone is not enough, however, because the Riesz representer must itself be estimated. To control the resulting first-stage error, we combine the plug-in derivative and the Riesz correction in a Neyman-
16
orthogonal score (Neyman, 1959; Chernozhukov et al., 2018), so that small errors in estimating either nuisance affect the final estimator only at second order. We now introduce these two ingredients in turn. We begin with the Riesz representer, which provides the key representation of the derivative bias. The construction in this subsection is not original to this paper: the Riesz representer for derivative functionals goes back to classical semiparametric theory and has been developed in its modern form by Chernozhukov et al. (2022) and Hirshberg and Wager (2021); we adopt their formulation directly and merely apply it to the network Jacobian. Definition 1 (Riesz Representer; Chernozhukov et al., 2022; Hirshberg and Wager, 2021). For each pair (i, j), the Riesz representer αij satisfies
∂h(Zt ) E[αij (Zt ) h(Zt )] = E ∂yj
(9)
for all sufficiently smooth test functions h in a function class H for which h 7→ E[∂j h(Zt )] is a bounded linear functional, where Zt = (Yt , Xt ). This definition says that the derivative operator can be represented as an L2 inner product with the function αij . Intuitively, instead of differentiating a function directly, we can multiply it by αij and take expectations to recover the same object. Applying Definition 1 with h(Zt ) = m0,i (Zt ) − m̂i (Zt ), the bias term in (8) becomes ∂(m0,i − m̂i ) E (Zt ) = E αij (Zt ) m0,i (Zt ) − m̂i (Zt ) . ∂yj
(10)
This is the key representation: it rewrites the derivative bias as an expectation involving the prediction error m0,i (Zt ) − m̂i (Zt ). Although m0 itself is unobservable, we do observe Yt+1,i = m0,i (Zt ) + εt+1,i with E[εt+1,i | Zt ] = 0, so the unobserved prediction error can be replaced, in expectation, by the observable residual Yt+1,i − m̂i (Zt ). Together, these 17
two steps deliver what we will refer to as an orthogonal bias correction. Informally, the correction has two desirable properties. First, it is a bias correction: starting from the biased plug-in derivative, we add the term E[αij (Zt )(Yt+1,i − m̂i (Zt ))], which estimates the bias from (10) using only quantities the analyst can compute from the data. Second, it is orthogonal in the Neyman sense: small errors in the two estimated nuisance functions, the regression m̂ and the Riesz weight α̂ij , do not contaminate the corrected estimator at first order, so the resulting estimator is robust to slow convergence in either nuisance. √ This is exactly the property that allows valid edge-level inference at the parametric n rate even though m0 is estimated by flexible machine learning; we make this formal in the next subsection. We now define the score that underlies debiased estimation. For each pair (i, j), let
ψij (Zt , Yt+1 ; θij , m, α) =
∂mi (Zt ) + αij (Zt ) Yt+1,i − mi (Zt ) − θij , ∂yj
(11)
where θij is the target parameter, mi is the i-th coordinate of the regression function, and 0 , m0 , αij ), this score satisfies αij is the corresponding Riesz representer. At the truth (θij
Neyman orthogonality: first-order perturbations in either the regression nuisance mi or the Riesz nuisance αij do not affect the moment condition (see Appendix D.1). Intuitively, the first term is the biased plug-in derivative, and the second term uses prediction residuals to correct that bias in a way that is robust to first-stage estimation error. This √ orthogonality is what allows edge-level inference at the parametric n rate even though m0 is estimated by flexible machine learning. For the linear-index model in Section 3.2, the Riesz representer admits a closed form. ⊤ N +d+1 Let ri (Zt ) = [Yt⊤ , x⊤ and let uj denote the basis vector corresponding to t,i , 1] ∈ R
the j-th coordinate of Yt . Then αij (Zt ) = ri (Zt )⊤ Ci−1 uj ,
Ci = E[ri (Zt ) ri (Zt )⊤ ].
18
(12)
For richer model classes, αij can be estimated numerically using the loss-based method of Chernozhukov et al. (2022).
3.3.4
Algorithm and inference
We now summarize the full debiased network-recovery procedure. Algorithm 1. CascadeNet Debiased Network Recovery Input: panel {(Yc,t , Xc,t )}c≤C, t≤Tc ; class M; folds K; regularization λ. Output: Jˆdeb , standard errors, CIs, and edge p-values. 1. Form transitions (Zc,t , Yc,t+1 ) with Zc,t = (Yc,t , Xc,t ), and partition trajectories c = 1, . . . , C into K folds I1 , . . . , IK . 2. For each fold k = 1, . . . , K, using only trajectories not in Ik : (a) fit m̂(−k) ∈ M by minimizing (4); (−k)
(b) estimate the Riesz representers α̂ij merical loss otherwise);
(closed form for the linear index, nu(−k)
(c) for each held-out (c, t) with c ∈ Ik , evaluate ψ̂ij,c,t = ∂yj m̂i (−k) (−k) α̂ij (Zc,t ) Yc,t+1,i − m̂i (Zc,t ) . 3. Aggregate to the trajectory level and average across folds: P P P c −1 ψ̂ij,c,t , n = C. n−1 k c∈Ik (Tc − 1)−1 Tt=1
(Zc,t ) +
deb θ̂ij
=
deb ; estimate the trajectory-level score variance and report stan4. Set Jˆdeb [i, j] = θ̂ij dard errors, 100(1 − τ )% CIs, and Wald tests for H0 : J0 [i, j] = 0.
Two implementation details deserve emphasis. The first is the role of cross-fitting, which appears in step 2 of Algorithm 1. If the same trajectories were used both to estimate the nuisance functions (m̂, α̂) and to evaluate the score, the score would inherit a self-influence bias: any way in which m̂ overfits a particular trajectory would also contaminate the residual computed on that trajectory, and the bias correction would be evaluated on the very residuals it was tuned to fit (Chernozhukov et al., 2018). Cross-fitting eliminates this dependence by partitioning the cascades c = 1, . . . , C into K folds, fitting m̂(−k) and α̂(−k) on the K − 1 training folds, and evaluating the score only on the held-out fold Ik ; averaging the resulting per-fold scores reproduces the trajectory-level CLT used
19
in Theorem 1. The partition is performed at the trajectory level rather than the observation level, which preserves within-trajectory dependence and matches the asymptotic regime where the number of trajectories n = C grows. The second detail is selection of the regularization parameter λ in (4). In our experiments λ is chosen by held-out crossvalidation on a fold of trajectories disjoint from the cross-fitting folds used to evaluate the score, so that the data used to tune λ play no role in the inference step. In practice this is implemented as an outer hyperparameter-tuning split nested inside each cross-fitting fold; the same Riesz-loss validation is used to tune α̂ij when it is estimated numerically. Both choices are standard in the debiased-ML literature (Chernozhukov et al., 2018, 2022) and are necessary for the asymptotic guarantees in Theorem 1 to apply. The orthogonal score delivers valid edge-level inference under standard regularity conditions. For the asymptotic theory, we treat each cascade trajectory as an independent sampling unit and let the number of trajectories n grow; cross-fitting is performed at the trajectory level. For each trajectory c, define the trajectory-level score Tc −1 1 X 0 ψij Zc,t , Yc,t+1 ; θij , m0 , αij , ϕij (Wc ) = Tc − 1 t=1 c −1 where Wc = {(Yc,t , Xc,t , Yc,t+1 )}Tt=1 .
Theorem 1 (Asymptotic Normality of the Debiased Jacobian). Under the regularity conditions listed in Assumptions 1–5 in Appendix A, for each fixed pair (i, j), d √ deb 0 n θ̂ij − θij → − N (0, Vij ),
Vij = E[ϕij (Wc )2 ].
(13)
Moreover, the following cross-fitted sample analogue V̂ij is a consistent estimator of Vij : K
1 X X (k) V̂ij = ϕ̂ (Wc )2 , n k=1 c∈I ij k
(k) ϕ̂ij (Wc ) =
Tc −1 1 X (−k) deb ψij Zc,t , Yc,t+1 ; θ̂ij , m̂(−k) , α̂ij , Tc − 1 t=1
(14)
20
(−k)
where Ik are the cross-fitting folds and (m̂(−k) , α̂ij
) are the nuisance estimates fitted without fold
k. Theorem 1 is a direct application of the debiased machine learning template of Chernozhukov et al. (2018) and the Riesz-representer framework of Chernozhukov et al. (2022) and Hirshberg and Wager (2021); our contribution is not the underlying inference machinery, which is by now well established in the econometrics and statistics literatures, but the specific identification of the network Jacobian J0 as a smooth functional of the transition regression and the verification that the resulting orthogonal score is well behaved under repeated sampling of cascade trajectories. The theorem implies that each entry of the debiased Jacobian admits asymptotically valid standard errors and confidence intervals, even though the underlying transition function m0 is estimated by flexible machine learning. Theorem 1 yields immediate edge-level inference. Let V̂ij be the variance estimator q deb deb V̂ij /n, an asymptotically valid b θ̂ij ) = in (14). Then the standard error of θ̂ij is se( deb deb ), and a test of ± z1−τ /2 se( b θ̂ij 100(1 − τ )% confidence interval for the edge (j → i) is θ̂ij deb deb ). These inferential objects are /se( b θ̂ij H0 : J0 [i, j] = 0 is based on the Wald statistic θ̂ij
useful in several ways. First, they allow the researcher to distinguish statistically meaningful edges from noise, rather than relying only on point estimates or heuristic rankings. Second, they allow uncertainty-aware edge ranking, which is especially important when the recovered network will be used for targeting, intervention design, or downstream policy analysis. Third, they make it possible to communicate uncertainty transparently in applications where decisions based on inferred influence links are consequential. Theorem 1 is stated under correct specification of the transition function, meaning that m0 belongs to the model class M used to estimate m̂. When the model is misspecified, the same debiasing argument delivers asymptotic normality around the average Jacobian of the pseudo-true regression function, that is, the best approximation to m0 within the chosen class, rather than around the true J0 . The method therefore remains well behaved 21
under misspecification, but inherits an approximation bias whose magnitude depends on the expressive power of the chosen model class. Appendix D.6 provides the formal argument.
4
Validation We validate CascadeNet on both synthetic and real data. The synthetic exercise
allows us to assess recovery performance under controlled conditions where the datagenerating process is known and the true average Jacobian can be computed. The empirical application assesses whether the recovered network aligns with an external benchmark in a real-world epidemiological setting.
4.1
Synthetic Validation
4.1.1
Simulation
We evaluate CascadeNet on synthetic cascades generated from a known network under a range of diffusion mechanisms. The goal is to test whether the method can recover the DGP-implied influence structure when the diffusion model is correctly specified, mildly misspecified, or strongly misspecified relative to classical baselines. We generate a directed Erdős–Rényi network with N = 64 nodes and edge probability p = 0.05, yielding on average approximately 195 directed edges. Edge weights are drawn independently from [0.2, 0.8]. We pair this network with nine data-generating processes that span three broad families. The first family, pairwise transmission, captures settings in which influence flows one edge at a time: each active node attempts to activate each of its neighbors, either independently (Independent Cascade) or once a weighted sum crosses a node-specific threshold (Linear Threshold). The second family, continuous-time SI, models infection as the result of waiting-time hazards along edges, with Exponential-SI and Power-law22
SI corresponding to different shapes of that hazard. The third family, aggregate-effect and nonlinear dynamics, departs from clean pairwise propagation: SIS and SIR allow recovery and (for SIR) absorbing immunity, Complex Contagion requires multiple simultaneously active neighbors before a node activates, Hawkes generates self-exciting bursts whose intensity depends on past events, and a flexible Nonlinear DGP combines polynomial, oscillatory, and threshold-like features. Together these nine DGPs span absorbing versus non-absorbing dynamics, pairwise versus aggregate transmission, and both wellspecified and strongly misspecified environments for the competing methods. We rely on the classical formulations developed in Kempe et al. (2003), Gomez-Rodriguez et al. (2011), Kermack and McKendrick (1927), Centola and Macy (2007), and Hawkes (1971); the explicit one-step update equation and parameter choices for each DGP are given in Appendix B. For each DGP, we generate Ctrain = 1,000 training cascades of length T = 10–30, each initialized from two randomly chosen seed nodes. We evaluate each method on an independent set of Ctest = 1,000 cascades. To compute the true population target J0 , we use a separate held-out sample of Ceval = 5,000 cascades. We compare six approaches. Pairwise is a simple temporal-precedence heuristic that scores edge j → i by the fraction of cascades in which node j becomes active before node i (in the spirit of Gomez-Rodriguez et al., 2010). NetInf (Gomez-Rodriguez et al., 2010) performs likelihood-based network inference under an exponential independent-cascade kernel. NETRATE (Gomez-Rodriguez et al., 2011) estimates pairwise transmission rates from first-infection timestamps by convex optimization. DANI (Ramezani et al., 2024) estimates pairwise transmission probabilities from the observed panel using a Markovtransition formulation. LTMLE is longitudinal targeted minimum-loss estimation (van der Laan and Gruber, 2012) adapted to panel cascade data. CascadeNet is our method with sigmoid link, λ = 10−4 , 300 training epochs, and 2-fold Riesz debiasing. Because the estimand in Section 3.3 is the average Jacobian rather than the binary
23
adjacency matrix, our primary simulation benchmark compares estimated edge scores to the true DGP-implied mean Jacobian ∂m0,i (Yt , Xt ) . J0 [i, j] = E ∂yj
Different diffusion models can induce very different marginal influence magnitudes even on the same underlying graph: an edge that is important under IC need not have the same average effect under Hawkes or SIS dynamics. For that reason, adjacency recovery alone is not the right evaluation target for our estimator. We compute the true J0 numerically from each DGP’s conditional expectation operator, using finite differences evaluated over the independent held-out sample of 5,000 cascades. Our primary metric is the Pearson correlation between each method’s estimated edge-score vector and |J0 |, using off-diagonal entries only. We focus on |J0 | because the network-recovery task is primarily about ranking influence strength, and in these simulations the economically relevant comparison is the magnitude of marginal influence rather than its sign. As a secondary metric, we report precision-recall performance for detecting the support of the underlying adjacency matrix W . This secondary benchmark is useful for comparing topology recovery, but our primary criterion remains alignment with the true average Jacobian, since that is the estimand delivered by CascadeNet. Table 1 reports the Pearson correlation between estimated edge scores and the true mean Jacobian |J0 | across all nine diffusion models. CascadeNet achieves the highest correlation in every setting. Two patterns stand out. First, the classical baselines are model-sensitive: each tends to perform well on the family most closely aligned with its maintained assumptions but degrades sharply outside that family. NetInf performs strongly on LT but much worse on Hawkes; LTMLE performs well on the continuous-time SI models but deteriorates on Hawkes and the flexible nonlinear DGP; Pairwise performs especially poorly on epidemic
24
Table 1. Network Recovery: Pearson r between estimated scores and true Jacobian |J0 | (C = 1,000 training cascades, ER network with N = 64, |E| = 195). Bold indicates best method per row. ∆ is the improvement of CascadeNet over the best baseline; stars denote significance via Fisher z-test: *** p < 0.001. Baselines
CascadeNet
DGP
Pairwise
NetInf
NETRATE
DANI
LTMLE
Ours
∆
IC LT Exp-SI PL-SI SIS SIR Complex Hawkes Nonlinear
0.32 0.45 0.37 0.43 0.10 0.14 0.41 −0.04 0.23
0.65 0.70 0.50 0.32 0.58 0.51 0.57 0.32 0.27
0.64 0.63 0.31 0.36 0.44 0.50 0.33 0.30 0.25
0.47 0.58 0.55 0.60 0.26 0.49 0.59 0.19 0.25
0.66 0.68 0.73 0.77 0.29 0.64 0.70 0.20 0.46
0.84 0.78 0.80 0.81 0.86 0.84 0.80 0.82 0.54
+0.18∗∗∗ +0.08∗∗∗ +0.07∗∗∗ +0.04∗∗∗ +0.28∗∗∗ +0.21∗∗∗ +0.10∗∗∗ +0.50∗∗∗ +0.08∗∗∗
models with recovery/removal, such as SIS and SIR. Second, CascadeNet is uniformly strong: its correlation never falls below 0.54, and it delivers the best edge ranking in every environment considered. Figure 1 visualizes these comparisons. We also report precision-recall curves for adjacency recovery in Figure 2 (Appendix C). CascadeNet maintains substantially higher precision at moderate-to-high recall on the aggregate-effect and nonlinear models, while remaining competitive with the strongest classical baselines on pairwise-transmission DGPs where those methods are relatively well specified. Pairwise 1.0
*
Pearson r
0.8
NetInf
DANI
LTMLE
*
*
*
*
NETRATE
*
CascadeNet
*
*
*
0.6 0.4 0.2 0.0
IC
LT
Exp-SI
PL-SI
SIS
SIR
Complex
Hawkes
Nonlinear
Figure 1. Pearson correlation between estimated scores and true Jacobian |J0 | across nine diffusion models (C = 1,000, ER network). CascadeNet achieves the highest correlation in all nine settings. Classical baselines are inconsistent: each excels on its “native” DGP family but fails on others. 25
Taken together, the simulation results support the main claim of the paper: when the true diffusion mechanism is unknown, a flexible Jacobian-based estimator with orthogonal debiasing provides a more robust recovery strategy than methods built around a single parametric propagation model.
4.2
Empirical Application: COVID-19 Transmission in Spain We now apply CascadeNet to a real-world epidemiological setting: recovering inter-
province transmission structure from daily COVID-19 case counts in Spain. During the COVID-19 pandemic, a central policy problem was to understand which regions were transmitting infection to which other regions. This information was relevant for targeted interventions such as regional travel restrictions, testing deployment, and surveillance prioritization. While the transmission network itself is not observed, inter-province mobility flows provide a useful external benchmark for the plausibility of the recovered edge rankings. Mobility is not a literal ground truth for transmission, because behavioral adaptation, public-health policy, and local conditions all affect realized contagion, but stronger alignment with mobility is suggestive evidence that the inferred network is capturing meaningful transmission opportunities. We use the dataset of Murphy et al. (2021), which combines daily COVID-19 case counts from Spain’s Centro Nacional de Epidemiología with inter-province mobility flows from the Ministerio de Fomento. The case panel records daily new infections for N = 52 provinces over T = 450 days, from January 2020 to March 2021. Following Murphy et al. (2021), we normalize each province’s incidence series by its mean daily incidence to improve comparability across provinces of different sizes. To construct multiple panels, we use sliding windows of 56 days with stride 28, yielding C = 15 partially overlapping trajectory windows. Because these windows overlap, the formal independent-trajectory asymptotic theory of Section 3.3 should be viewed here as suggestive rather than exact. We therefore interpret this application primarily as an 26
external validation exercise for the recovered edge rankings, rather than as a literal implementation of the asymptotic sampling regime used in the theorem. As an external benchmark, we construct a thresholded mobility network from interprovince traveler flows. Following Murphy et al. (2021), we log-transform the flows, symmetrize them, and threshold the resulting matrix to an average degree of approximately 10, yielding |E| = 519 undirected benchmark edges. In the evaluation below, we compare off-diagonal estimated edge scores to this benchmark. Table 2 reports empirical performance. We evaluate whether estimated edge rankings align with the benchmark mobility network using the break-even point (BEP) of the precision-recall curve, AUC-PR, and the Pearson correlation between estimated edge scores and log mobility flow. Table 2. Network Recovery on Spain COVID-19 Data (N = 52 Provinces, |E| = 519 Mobility Edges). Pearson r measures correlation between estimated edge scores and log inter-province mobility flow. Bold indicates best method. Pearson r Method
BEP
AUC-PR
r
p-value
Pairwise Granger NetInf NETRATE DANI LTMLE CascadeNet
0.194 0.236 0.174 0.190 0.176 0.176 0.298
0.202 0.209 0.186 0.188 0.196 0.195 0.263
−0.012 +0.025 +0.004 +0.024 −0.011 −0.012 +0.139
0.53 0.20 0.83 0.22 0.58 0.52 < 10−12
CascadeNet achieves the highest BEP and AUC-PR among all methods. The benchmarknetwork density is 0.196, which is the random baseline for precision-recall evaluation; most baselines lie only slightly above that level, whereas CascadeNet improves on it substantially. On the more demanding continuous benchmark, CascadeNet is also the only method whose estimated edge scores correlate significantly with log mobility flows. The baseline correlations are all close to zero. This result is notable because CascadeNet is trained only on case-count dynamics and never observes mobility data during estimation. The positive alignment with mobility 27
therefore emerges from the recovered transmission structure itself, rather than from any direct use of the benchmark. The empirical benchmark is intentionally demanding. Mobility is only an imperfect proxy for transmission: lockdowns, masking, behavioral responses, underreporting, and within-province spread all weaken the mapping from traveler flows to realized infections. The relevant question is therefore not whether a method perfectly reconstructs the benchmark network, but whether it extracts any meaningful cross-region transmission signal from noisy epidemic data. By that standard, CascadeNet appears to succeed where the classical baselines do not. Its moderate but clearly positive alignment with mobility suggests that the flexible, Jacobian-based approach is recovering transmission pathways that are plausibly related to real inter-provincial contact intensity, even in an environment where the true contagion mechanism is far more complex than the kernels assumed by standard cascade-inference methods.
5
Managerial Implications and Conclusion This paper develops CascadeNet, a flexible framework for network recovery from
cascade data. The method is built on a simple but general insight: under a Markovian diffusion process, the Jacobian of the one-step transition function summarizes how one agent’s current state affects another agent’s next-period outcome. By estimating this transition function flexibly and then applying Neyman-orthogonal debiasing via the Riesz representer, CascadeNet recovers an average influence network while correcting the regularization bias that would otherwise distort the naive plug-in Jacobian. The result is a network estimator that combines flexibility with formal statistical inference on individual edges. The managerial implications are substantial. In many practical settings, the analyst does not know how influence propagates: a marketing manager may not know whether
28
consumers respond to pairwise persuasion, peer reinforcement, or threshold behavior, and a public-health authority faces transmission dynamics shaped by mobility, policies, and behavioral adaptation that no single parametric model captures well. Our results show that this uncertainty matters: methods tailored to one propagation model can deteriorate sharply when the assumed kernel is wrong, whereas CascadeNet learns the transition function directly from the data without requiring the analyst to commit ex ante to a specific diffusion mechanism. This robustness is particularly valuable because many downstream decisions, such as seed selection in marketing (Kempe et al., 2003), travel restrictions and testing allocation in epidemiology, and contagion monitoring in financial settings, depend on identifying which links in the network are most consequential. In all of these applications, ranking edges by estimated importance is useful, but ranking them with confidence intervals is substantially more valuable. The validation results support these conclusions. Across nine synthetic diffusion models spanning pairwise, epidemic, aggregate-effect, and nonlinear environments, CascadeNet achieves the highest recovery accuracy in every setting, and in the Spain COVID19 application it is the only method whose recovered edge scores align significantly with inter-province mobility flows, even though mobility data are never used in estimation. The debiasing step is central to this performance: when cascade data are sparse and regularization is strong, the Riesz-based correction is often the difference between recovering a useful network signal and recovering almost none at all (Appendix D.5). Several extensions remain for future work, including adaptive topology learning when the adjacency structure is fully unknown, uniform inference for many edges simultaneously, continuous-time diffusion with flexible kernels, and inference under singlecascade or dependent-trajectory sampling. Finally, predicting long-run or steady-state outcomes from the estimated transition function is a natural next step; we pursue this direction in companion work on equilibrium learning, where the core challenge is valid inference on the fixed point of a flexibly estimated map.
29
More broadly, the paper’s message is that network recovery need not force a tradeoff between flexibility and inferential rigor. By combining modern machine learning with orthogonal debiasing, CascadeNet makes recovered networks not just descriptive objects, but reliable inputs into managerial and policy analysis.
References Bai, S., Kolter, J. Z., and Koltun, V. (2019). Deep equilibrium models. Advances in Neural Information Processing Systems, 32. Ballester, C., Calvó-Armengol, A., and Zenou, Y. (2006). Who’s who in networks. Wanted: The key player. Econometrica, 74(5):1403–1417. Bramoullé, Y., Djebbari, H., and Fortin, B. (2009). Identification of peer effects through social networks. Journal of Econometrics, 150(1):41–55. Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1):C1–C68. Chernozhukov, V., Newey, W. K., and Singh, R. (2022). Automatic debiased machine learning of causal and structural effects. Econometrica, 90(3):967–1027. Chu, J. and Manchanda, P. (2016). Quantifying cross and direct network effects in online consumer-to-consumer platforms. Marketing Science, 35(6):870–893. DeGroot, M. H. (1974). Reaching a consensus. Journal of the American Statistical Association, 69(345):118–121. Eckles, D., Kizilcec, R. F., and Bakshy, E. (2016). Estimating peer effects in networks with peer encouragement designs.
Proceedings of the National Academy of Sciences,
113(27):7316–7322. 30
Eckles, D., Karrer, B., and Ugander, J. (2017). Design and analysis of experiments in networks: Reducing bias from interference. Journal of Causal Inference, 5(1):20150021. Farrell, M. H., Liang, T., and Misra, S. (2021). Deep neural networks for estimation and inference. Econometrica, 89(1):181–213. Farrell, M. H., Liang, T., and Misra, S. (2020). Deep learning for individual heterogeneity: An automatic inference framework. arXiv preprint arXiv:2010.14694. Hamilton, W., Ying, Z., and Leskovec, J. (2017). Inductive representation learning on large graphs. Advances in Neural Information Processing Systems, 30. Hirshberg, D. A. and Wager, S. (2021). Augmented minimax linear estimation. The Annals of Statistics, 49(6):3206–3227. Huang, L. (2023). Equilibrium learning. Term paper for 15.838, Massachusetts Institute of Technology. Jiang, Z., Sun, Q., Zhao, P., and Ye, M. (2022). Estimating causal effects on networks: A graph neural network approach. arXiv preprint arXiv:2209.01795. Centola, D. and Macy, M. (2007). Complex contagions and the weakness of long ties. American Journal of Sociology, 113(3):702–734. Hawkes, A. G. (1971). Spectra of some self-exciting and mutually exciting point processes. Biometrika, 58(1):83–90. Kato, T. (1995). Perturbation Theory for Linear Operators. Springer, Berlin. Kempe, D., Kleinberg, J., and Tardos, É. (2003). Maximizing the spread of influence through a social network. Proceedings of the 9th ACM SIGKDD, pages 137–146. Kermack, W. O. and McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London A, 115(772):700–721. 31
Leung, M. P. (2022). Causal inference under approximate neighborhood interference. Econometrica, 90(1):267–293. Ma, Y. and Zhang, K. (2021). Deconfounded recommendation for alleviating bias amplification. Proceedings of the 27th ACM SIGKDD, pages 1133–1143. Manski, C. F. (1993). Identification of endogenous social effects: The reflection problem. The Review of Economic Studies, 60(3):531–542. Neyman, J. (1959). Optimal asymptotic tests of composite hypotheses. In Grenander, U., editor, Probability and Statistics: The Harald Cramér Volume, pages 213–234. Wiley, New York. Stone, C. J. (1982). Optimal global rates of convergence for nonparametric regression. The Annals of Statistics, 10(4):1040–1053. van der Laan, M. J. and Gruber, S. (2012). Targeted minimum loss based estimation of causal effects of multiple time point interventions. The International Journal of Biostatistics, 8(1):9. Ye, Z., Zhang, Z., Zhang, D., Zhang, H., and Zhang, R. (2025). Deep-learning-based causal inference for large-scale combinatorial experiments: Theory and empirical evidence. Management Science, Articles in Advance, https://doi.org/10.1287/mnsc.2024.04625. Wager, S. and Xu, K. (2021).
Experimenting in equilibrium.
Management Science,
67(11):6694–6715. Gomez-Rodriguez, M., Leskovec, J., and Krause, A. (2010). Inferring networks of diffusion and influence. Proceedings of the 16th ACM SIGKDD, pages 1019–1028. Gomez-Rodriguez, M., Balduzzi, D., and Schölkopf, B. (2011). Uncovering the temporal dynamics of diffusion networks. Proceedings of the 28th ICML, pages 561–568.
32
Myers, S. A. and Leskovec, J. (2010). On the convexity of latent social network inference. Advances in Neural Information Processing Systems, 23. Murphy, C., Laurence, E., and Allard, A. (2021). Deep learning of contagion dynamics on complex networks. Nature Communications, 12:4720. Qiu, J., Tang, J., Ma, H., Dong, Y., Wang, K., and Tang, J. (2018). DeepInf: Social influence prediction with deep learning. Proceedings of the 24th ACM SIGKDD, pages 2110–2119. Ramezani, M., Rabiee, H. R., Tahani, M., and Rajabi, A. (2024). DANI: Fast diffusion aware network inference with preserving topological structure property. Scientific Reports, 14:5765.
Appendix A: Assumptions The assumptions below adapt the standard regularity conditions of the debiased machine learning literature, in particular Chernozhukov et al. (2018) and the Riesz-representer framework of Chernozhukov et al. (2022) and Hirshberg and Wager (2021), to the panel structure of cascade data. None of these conditions are new to this paper; they are imposed verbatim or with minor modifications relative to those references. For the asymptotic theory, we treat each cascade trajectory as an independent sampling unit. Let Tc −1 Wc = {(Yc,t , Xc,t , Yc,t+1 )}t=1 ,
c = 1, . . . , n,
denote the c-th observed trajectory, where n = C is the number of cascades, and define Zc,t := (Yc,t , Xc,t ). For each fixed pair (i, j), recall that the target parameter is
0 θij := E
∂m0,i (Yt , Xt ) . ∂yj
We impose the following assumptions. 33
Assumption 1 (Independent trajectories). The trajectories W1 , . . . , Wn are independent and identically distributed. Their lengths satisfy 2 ≤ Tc ≤ T̄ < ∞ for all c, for some fixed constant T̄ . Assumption 2 (Smoothness and bounded moments). For each coordinate i, the regression function m0,i (z) is continuously differentiable in Y . For each pair (i, j), ∂m0,i /∂yj (Zc,t ) ∈ L2 (PZ ) and E (∂m0,i /∂yj (Zc,t ))2 < ∞. Moreover, the regression residual εc,t+1,i := Yc,t+1,i − m0,i (Zc,t ) satisfies E[εc,t+1,i | Zc,t ] = 0 and E[ε2c,t+1,i ] < ∞. Assumption 3 (Riesz existence and boundedness). For each fixed pair (i, j), the linear functional h 7→ E[∂h(Zt )/∂yj ] is bounded on the function class H in Definition 1. Hence there exists a unique Riesz representer αij ∈ L2 (PZ ) such that E[αij (Zt ) h(Zt )] = E[∂h(Zt )/∂yj ] for all h ∈ H, and E[αij (Zt )2 ] < ∞. (−k)
Assumption 4 (Cross-fitted nuisance rates). Let m̂(−k) and α̂ij
denote the nuisance esti-
mators fit on the complement of fold Ik . Then for each fold k, (−k)
∥m̂i
(−k)
− m0,i ∥L2 (PZ ) = op (n−1/4 ),
∥α̂ij
− αij ∥L2 (PZ ) = op (n−1/4 ),
and (−k)
∂ m̂i ∂yj
−
∂m0,i ∂yj
= op (1). L2 (P
Z)
Assumption 5 (Finite trajectory-level score variance). For each fixed pair (i, j), define the oracle trajectory-level score Tc −1 1 X 0 ϕij (Wc ) = ψij Zc,t , Yc,t+1 ; θij , m0 , αij , Tc − 1 t=1
where ψij is the score in (11). Then E[ϕij (Wc )] = 0 and 0 < Vij := E[ϕij (Wc )2 ] < ∞. Assumptions 1–5 are standard in debiased machine learning (Chernozhukov et al., 2018, 2022; Hirshberg and Wager, 2021), adapted here to repeated sampling of cascade 34
trajectories. Assumption 4 in particular is the now-canonical op (n−1/4 ) product-rate condition introduced by Chernozhukov et al. (2018), which ensures that the product of nuisance errors is op (n−1/2 ), which is sufficient for valid debiasing. We do not view these assumptions as a contribution of this paper; the contribution is the application of this established framework to a setting (the network Jacobian recovered from cascade data) where it has not previously been deployed.
Appendix B: Diffusion-Model Specifications for the Synthetic Validation This appendix describes the nine data-generating processes used in Section 4.1. Throughout, W ∈ [0, 1]N ×N denotes the directed weighted adjacency matrix described in Section 4.1 (with Wji the weight on the directed edge j → i), Yc,t ∈ {0, 1}N denotes the activation state of the N agents at time t in cascade c, and σ(·) denotes the logistic sigmoid. All cascades start from two randomly chosen seed nodes. We summarize the one-step transition that defines each DGP.
Independent Cascade (IC) (Kempe et al., 2003).
Adoption is irreversible. At each step,
every newly active node j at time t makes a single attempt to activate each currently inactive neighbor i, succeeding independently with probability Wji :
Pr(Yc,t+1,i = 1 | Yc,t ) = 1 −
Y
(1 − Wji ),
Yc,t+1,i ≥ Yc,t,i .
j∈Nnew (t)
We set diagonal entries of W to zero so that nodes do not self-activate.
35
Each node i has a threshold θi drawn uni-
Linear Threshold (LT) (Kempe et al., 2003).
formly from (0, 1]. Adoption is irreversible. At time t + 1, ( ) X Yc,t+1,i = 1 Wji Yc,t,j ≥ θi ∨ Yc,t,i . j
We rescale rows of W so that
P
j Wji ≤ 1 for each i, ensuring well-defined thresholds.
Exponential-SI (Gomez-Rodriguez et al., 2011).
Continuous-time SI with exponential
transmission times. The waiting time for an active j to infect an inactive i is exponential with rate λji = Wji . We discretize time with unit step ∆t = 1 and obtain the one-step activation probability Y
Pr(Yc,t+1,i = 1 | Yc,t ) = 1 −
exp −Wji ∆t ,
Yc,t+1,i ≥ Yc,t,i .
j:Yc,t,j =1
Same as Exponential-SI but with a power-
Power-law-SI (Gomez-Rodriguez et al., 2011).
law transmission hazard hji (τ ) = (α − 1) τ −α Wji for elapsed time τ since j’s infection. We use α = 2.0 and discretize as in Exponential-SI.
SIS (Kermack and McKendrick, 1927).
Susceptible–Infected–Susceptible dynamics with
infection from active neighbors and recovery rate γ: h i Y Pr(Yc,t+1,i = 1 | Yc,t ) = (1 − Yc,t,i ) 1 − (1 − Wji Yc,t,j ) + Yc,t,i (1 − γ). j
We set γ = 0.3. Adoption is reversible (non-absorbing).
36
SIR (Kermack and McKendrick, 1927).
Susceptible–Infected–Recovered dynamics. In-
fected nodes recover (and become permanently immune) with rate γ:
Yc,t+1,i =
1{at least one active neighbor infects i}
if Yc,t,i = 0 and i is susceptible,
Yc,t,i
otherwise (including recovered),
where the infection probability is 1 −
Q
j (1 − Wji Yc,t,j ) and recovery occurs independently
each step with probability γ = 0.15.
Complex Contagion (Centola and Macy, 2007).
A node activates only when at least k of
its weighted neighbors are simultaneously active. Adoption is irreversible. We set k = 4: ( ) X Yc,t+1,i = 1 1{Wji > 0} Yc,t,j ≥ k ∨ Yc,t,i . j
Hawkes (self-exciting) (Hawkes, 1971).
Each active node generates an excitation that
decays exponentially in time and triggers offspring activations on its out-neighbors. We discretize so that the activation probability of i at time t + 1 is XX Pr(Yc,t+1,i = 1 | Yc,t ) = σ β0 + β1 Wji e−κ(t−s) 1{j activated at s} , s≤t
j
with decay κ = 0.5 and intensities β0 = −3, β1 = 1. Adoption is reversible.
Nonlinear DGP.
A flexible specification designed to violate the parametric assumptions
of every classical baseline. The transition is sigmoid-of-nonlinear-features: Pr(Yc,t+1,i = 1 | Yc,t ) = σ bi + (W Yc,t )i + 21 (W Yc,t )3i + sin 2(W Yc,t )i + ReLU (W Yc,t )i − 0.5 , with bias bi = −1. This DGP combines polynomial, oscillatory, and threshold-like features to test recovery under strongly misspecified competing methods. 37
For each DGP we generate Ctrain = 1,000 training cascades of length T = 10–30, Ctest = 1,000 test cascades, and Ceval = 5,000 held-out cascades used to compute the true mean Jacobian J0 via finite differences on the conditional expectation operator.
Appendix C: Additional Figures IC
0.8
0.6
0.6
0.6
0.0 0.0
Pairwise 0.29 NetInf 0.67 NETRATE 0.82 DANI 0.50 LTMLE 0.58 CascadeNet 0.82
0.2
Recall
0.6
0.8
0.0 0.0
1.0
PL-SI
1.0
Pairwise 0.59 NetInf 0.71 NETRATE 0.79 DANI 0.55 LTMLE 0.63 CascadeNet 0.79
0.2
0.4 0.2
0.4
Recall
0.6
0.8
0.0 0.0
1.0
SIS
1.0
0.8
0.6
0.6
0.6
0.2 0.0 0.0
Pairwise 0.54 NetInf 0.38 NETRATE 0.77 DANI 0.74 LTMLE 0.86 CascadeNet 0.82
0.2
Recall
0.6
0.8
0.0 0.0
1.0
Complex
1.0
0.4 0.2
0.4
Precision
0.8
0.4
Pairwise 0.13 NetInf 0.68 NETRATE 0.69 DANI 0.57 LTMLE 0.62 CascadeNet 0.86
0.2
Recall
0.6
0.8
0.0 0.0
1.0
Hawkes
1.0
0.6
0.6
0.6
0.2 0.0 0.0
0.2
0.4 0.2
0.4
Recall
0.6
0.8
1.0
Precision
0.8
Precision
0.8
Pairwise 0.52 NetInf 0.67 NETRATE 0.78 DANI 0.79 LTMLE 0.87 CascadeNet 0.89
0.0 0.0
Pairwise 0.06 NetInf 0.46 NETRATE 0.36 DANI 0.33 LTMLE 0.41 CascadeNet 0.81
0.2
Recall
0.6
0.8
1.0
Recall
0.6
0.8
1.0
0.6
0.8
1.0
0.8
1.0
Pairwise 0.08 NetInf 0.66 NETRATE 0.82 DANI 0.70 LTMLE 0.64 CascadeNet 0.79
0.2
0.4
Recall
Nonlinear
0.4 0.2
0.4
0.4
SIR
1.0
0.8
0.4
0.2
0.4 0.2
0.4
Pairwise 0.49 NetInf 0.60 NETRATE 0.67 DANI 0.71 LTMLE 0.82 CascadeNet 0.82
1.0
0.8 Precision
Precision
0.4 0.2
0.4
Precision
0.8
0.4
Exp-SI
1.0
0.8
0.2
Precision
LT
1.0
Precision
Precision
1.0
0.0 0.0
Pairwise 0.31 NetInf 0.35 NETRATE 0.36 DANI 0.25 LTMLE 0.60 CascadeNet 0.67
0.2
0.4
Recall
0.6
Figure 2. Precision-recall curves for edge detection across nine diffusion models (C = 1,000, ER network). CascadeNet (orange) dominates all baselines, with particularly large margins on SIS, Hawkes, and Nonlinear dynamics. Legend entries show BEP values.
38
Appendix D: Proofs D.1
Proof of Neyman Orthogonality (Equation 11)
The verification below follows the standard Neyman-orthogonality argument for Riesz-representer scores in Chernozhukov et al. (2018) and Chernozhukov et al. (2022); we record it here for completeness, applied to the specific score that targets the network Jacobian. We verify that the score in (11) satisfies Neyman orthogonality at the truth for 0 each fixed pair (i, j). For notational simplicity, write θ0 := θij , m0 := m0,i , α0 := αij , and
define ψ(W ; θ, m, α) :=
∂m(Z) + α(Z) Yi+ − m(Z) − θ, ∂yj
where W = (Z, Y + ) and Yi+ denotes the i-th coordinate of the next-period outcome. Let δm be an arbitrary square-integrable perturbation of m0 , and let δα be an arbitrary square-integrable perturbation of α0 . Consider the paths mr = m0 + r δm and αr = α0 + r δα . We show that the Gateaux derivatives of the population moment Ψ(θ, m, α) := E[ψ(W ; θ, m, α)] with respect to m and α, evaluated at (θ0 , m0 , α0 ), are both zero. Insensitivity with respect to m. Differentiate Ψ(θ0 , mr , α0 ) at r = 0: ∂δm (Z) d Ψ(θ0 , mr , α0 ) =E − α0 (Z) δm (Z) . dr ∂yj r=0 By Definition 1, E[α0 (Z) δm (Z)] = E[∂δm (Z)/∂yj ], so the derivative equals zero. Insensitivity with respect to α. Differentiate Ψ(θ0 , m0 , αr ) at r = 0: d Ψ(θ0 , m0 , αr ) = E δα (Z) Yi+ − m0 (Z) . dr r=0 Using iterated expectations and the conditional mean-zero property E[Yi+ − m0 (Z) | Z] = 0, E δα (Z) Yi+ − m0 (Z) = E δα (Z) E[Yi+ − m0 (Z) | Z] = 0.
39
Therefore the score is Neyman-orthogonal at the truth.
D.2
□
Proof of Theorem 1
The argument follows the standard cross-fitting / orthogonal-score template of Chernozhukov et al. (2018) and Chernozhukov et al. (2022); the only modification relative to those references is that the unit of cross-fitting is the cascade trajectory rather than an i.i.d. observation, which is required because within-trajectory observations are dependent. We prove asymptotic normality of the cross-fitted debiased estimator for a fixed 0 , α0 := αij , and m0 := m0,i . For each trajectory c, define the oracle pair (i, j). Write θ0 := θij
trajectory-level score Tc −1 1 X ψij Zc,t , Yc,t+1 ; θ0 , m0 , α0 , ϕ0 (Wc ) := Tc − 1 t=1
where ψij is the score in (11). By Assumption 5, E[ϕ0 (Wc )] = 0 and E[ϕ0 (Wc )2 ] = Vij ∈ (0, ∞). Let I1 , . . . , IK denote the cross-fitting folds. For each c ∈ Ik , define the estimated trajectory-level signal # " Tc −1 (−k) ∂ m̂i (Zc,t ) 1 X (−k) (−k) + α̂ij (Zc,t ) Yc,t+1,i − m̂i (Zc,t ) − θ0 . ϕ̂c := Tc − 1 t=1 ∂yj deb Then θ̂ij − θ0 = n−1
Pn
c=1 ϕ̂c , and we will show
√
1 deb n (θ̂ij − θ0 ) = √
n X
n c=1
ϕ0 (Wc ) + op (1),
from which the result follows by a central limit theorem.
40
Step 1: Decomposition. √
Write
n n 1 X 1 X deb ϕ0 (Wc ) + √ ϕ̂c − ϕ0 (Wc ) . − θ0 ) = √ n (θ̂ij n c=1 n c=1
(15)
The first term is the oracle empirical process; the second is the remainder from estimating the nuisances.
Step 2: CLT for the oracle term.
By Assumption 1, the trajectories W1 , . . . , Wn are i.i.d.
By Assumption 5, E[ϕ0 (Wc )] = 0 and E[ϕ0 (Wc )2 ] = Vij < ∞. The standard central limit theorem gives n
1 X d √ ϕ0 (Wc ) → − N (0, Vij ). n c=1 Step 3: Expansion of the remainder. (−k)
used to estimate m̂(−k) and α̂ij
(16)
Fix a fold k and condition on the training sample
. For each held-out trajectory c ∈ Ik ,
ϕ̂c − ϕ0 (Wc ) = Ac + Bc + Cc + Dc ,
where " # Tc −1 (−k) 1 X ∂ m̂i (Zc,t ) ∂m0 (Zc,t ) Ac := − , Tc − 1 t=1 ∂yj ∂yj Tc −1 1 X (−k) Bc := α̂ij (Zc,t ) − α0 (Zc,t ) Yc,t+1,i − m0 (Zc,t ) , Tc − 1 t=1
Cc := −
Tc −1 1 X (−k) α0 (Zc,t ) m̂i (Zc,t ) − m0 (Zc,t ) , Tc − 1 t=1
Dc := −
Tc −1 (−k) 1 X (−k) α̂ij (Zc,t ) − α0 (Zc,t ) m̂i (Zc,t ) − m0 (Zc,t ) . Tc − 1 t=1
Step 4: Orthogonality cancels the first-order terms.
Conditional on the training sample,
the held-out trajectories are independent of the nuisance estimators by cross-fitting. For
41
the terms Ac and Cc , the Riesz property gives " E[Ac +Cc | training sample] = E
# − m0 )(Z) (−k) − α0 (Z) m̂i − m0 (Z) training sample = 0. ∂yj
(−k)
∂(m̂i
For the term Bc , using E[Yc,t+1,i − m0 (Zc,t ) | Zc,t ] = 0,
E[Bc | training sample] = 0.
Thus the first-order nuisance effects cancel in expectation.
The remaining term Dc is second order. By
Step 5: Control of the second-order term. Cauchy–Schwarz,
(−k)
E[Dc | training sample] ≤ ∥α̂ij
(−k)
− α0 ∥L2 (PZ ) ∥m̂i
− m0 ∥L2 (PZ ) .
By Assumption 4, each factor is op (n−1/4 ), so their product is op (n−1/2 ), and hence E[Dc | training sample] = op (n−1/2 ). Because trajectory lengths are uniformly bounded by AsP sumption 1, averaging over trajectories preserves this rate, so n−1/2 nc=1 Dc = op (1). Similarly, the centered fluctuations of Ac , Bc , and Cc are negligible under Assumptions 2 and 4, since the derivative estimation error converges to zero in L2 , the nuisance estimators are P cross-fitted, and trajectory lengths are uniformly bounded. Hence n−1/2 nc=1 (Ac + Bc + Cc ) = op (1). Combining the previous two displays yields n 1 X √ ϕ̂c − ϕ0 (Wc ) = op (1). n c=1
Substituting (16) and (17) into (15),
Step 6: Conclusion. √
d
deb n (θ̂ij − θ0 ) → − N (0, Vij ),
42
Vij = E[ϕ0 (Wc )2 ].
(17)
This proves Theorem 1.
D.3
□
Graph Neural Network Model and Sparse Gradient
When the transition function involves complex nonlinear interactions that a linear index cannot capture, CascadeNet can be implemented as a Graph Neural Network using a GraphSAGE-style architecture (Hamilton et al., 2017): hi = σenc Wenc [yt,i , Xt,i ] ∈ RH , aggi = A h
√ [A = softmax(QK ⊤ / H)],
(18) (19) (20)
m̂i = Wdec [hi , aggi ].
Here hi is the node embedding for agent i, encoding its own state and covariates. The √ attention matrix A = softmax(QK ⊤ / H) captures how agent i aggregates information from its neighbors, with Q and K as learnable query and key matrices. The decoder combines the node embedding and the aggregated neighbor information to produce the predicted transition m̂i . With a known adjacency matrix W (no self-loops, Wii = 0), the embedding gradient has a sparse structure: ∇y hi , 0
i ∂ei = ∂yj 0, Wij ∇y hj j
if j = i, (21) if j ̸= i.
Proof. Decompose ei = [hi , aggi ]. By (18), hi depends only on (yi , Xi ), so ∂hi /∂yj = 0 for j ̸= i. By (19) with Wii = 0, ∂aggi /∂yi = Wii ∇yi hi = 0 and ∂aggi /∂yj = Wij ∇yj hj for j ̸= i. Concatenating gives (21).
□
43
D.4
Extension to Non-Markovian Dynamics
The Markovian assumption in (1) requires that Yc,t−1 affects Yc,t+1 only through Yc,t . When this assumption is violated, the framework can be extended by augmenting the state vector to include lagged values. Specifically, define the augmented state Ỹc,t = (Yc,t , Yc,t−1 , . . . , Yc,t−L ) for some lag order L ≥ 1. The transition function then becomes Ỹc,t+1 = m̃0 (Ỹc,t , Xc ; θ), which is Markovian in the augmented state. All results in the paper carry through with Y replaced by Ỹ , at the cost of increasing the dimensionality of the state from N to N (L + 1). In practice, one or two lags typically suffice for diffusion processes with short memory.
D.5
Supplementary Simulation: Debiasing Validation
To directly validate the Riesz debiasing correction with analytical ground truth, we present a simulation using a tanh DGP where the true Jacobian admits a closed form.
Data-generating process.
We simulate N = 20 nodes, T = 10 time steps, R = 800
trajectories:
yt+1,i = tanh α yt,i + β (W yt )i + γ xi + εt+1,i ,
εt+1,i ∼ N (0, σ 2 ),
(22)
with α = 0.3, β = 0.5, γ = 1.0, σ = 0.1.
Results.
Table 3 reports results over 1,000 evaluation points. The naive Jacobian is
nearly uninformative (r = 0.070) due to ℓ2 attenuation. The Riesz correction recovers the signal (r → 0.768), reducing Frobenius error by 1.70×.
44
Table 3. Debiasing Validation on Tanh DGP: Naive vs. Debiased Jacobian Estimates. Naive Mean Frobenius error Correlation r OLS R2 Error reduction
D.6
1.255 0.070 0.005
Debiased
0.738 0.768 0.589 1.70×
Robustness to Model Misspecification
In the main text, the asymptotic results are derived under the assumption that the transition model m(y, x; β) is correctly specified, so that there exists a parameter vector β0 satisfying E[Yt+1 | Zt ] = m(Zt ; β0 ),
Zt = (Yt , Xt ).
In practice, the transition model may be misspecified. In this appendix we briefly discuss how the results extend to this case.
Pseudo-true parameter.
When the model is misspecified, the estimator typically con-
verges not to a structural parameter but to a pseudo-true parameter β ⋆ . For example, under squared prediction loss, β ⋆ ∈ arg min E ∥Yt+1 − m(Zt ; β)∥2 . β
Equivalently, β ⋆ satisfies the population score condition E g(Zt ; β ⋆ ) Yt+1 − m(Zt ; β ⋆ ) = 0,
where g(Zt ; β) := ∂β m(Zt ; β). Under standard regularity conditions, the debiased estimator satisfies √ d n (β̂ deb − β ⋆ ) → − N (0, Σ⋆β ),
45
−1⊤ with sandwich covariance Σ⋆β = A−1 . The debiased Jacobian estimator therefore ⋆ Ω⋆ A⋆
continues to be asymptotically normal under misspecification, centered at the Jacobian implied by the best approximation β ⋆ within the model class. This approximation bias is expected to be small when the model class M is sufficiently flexible to closely approximate the true conditional expectation.
46