Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
104,383 characters · 8 sections · 22 citation commands
Design-Based Inference under Random Potential Outcomes
Keywords: Riesz representation; Local dependence; Variance estimation; Hilbert space methods; Dependency graphs.
\defcitealias{neyman1923}{Neyman (1923/1990)}
Randomised experiments remain central to the identification of causal effects. In the classical design-based perspective, pioneered by \citetalias{neyman1923} and rubin1974, potential outcomes are treated as fixed quantities, with all randomness attributed to the treatment assignment mechanism. This paradigm underpins finite-population inference, in which causal estimands are defined as functionals of a realised potential outcome schedule and identification is achieved solely through randomisation. Concretely, under the Neyman--Rubin framework the canonical estimand is the average treatment effect (ATE),
where for each unit only one of the two potential outcomes is observed. Identification relies on the fixed potential outcomes (FPO) assumption together with a known randomisation design, under which uncertainty arises exclusively from the assignment mechanism. We refer to this classical Neyman--Rubin formulation as the FPO framework throughout the paper.
Building on this foundation, a substantial body of recent work has extended design-based inference to more complex experimental settings. In particular, methods have been developed to accommodate interference and dependence among units, including network and spatial designs aronow2017estimating, savje2021average, athey2021design. These approaches preserve the FPO assumption while modelling dependence through exposure mappings or structured interference. Related contributions address robustness or misspecification in observational settings abadie2020sampling, imbens2004nonparametric, but typically achieve identification through explicit outcome modelling rather than randomisation.
Despite this progress, existing extensions remain within a finite-population logic. Potential outcomes are treated as fixed objects, and causal targets are therefore conditional on a single realised potential outcome schedule. As a result, FPO estimands are inherently tied to a particular realised world and do not yield stable mechanism-level targets in stochastic environments. This conditionality limits their relevance in broader applied contexts, where the object of interest concerns expected effects under repeated realisations of the same outcome-generating mechanism rather than effects for a specific realised sample. Indeed, in policy evaluation and public health decision-making, as well as in related applied settings, the object of interest is typically the expected effect of an intervention when implemented at scale, rather than the realised effect for the particular units observed in a past experiment. In this sense, the challenge is one of generalisation from a single realised experimental world to expectation-based causal targets defined over the outcome-generating mechanism.
This gap has motivated, among other considerations, the development of model-assisted approaches, including Bayesian formulations rubin1978bayesian, which treat potential outcomes as random variables but require explicit prior and likelihood specifications, propensity score methods rosenbaum1983central, which target expectation-based effects through parametric or semiparametric modelling of the selection mechanism, and double machine learning chernozhukov_doubledebiased_2018, which accommodates high-dimensional confounding via flexible prediction under additional regularity assumptions, among others. While these approaches are well suited for defining and estimating mechanism-level causal quantities, identification is achieved through modelling assumptions rather than through the randomisation design.
This contrast raises a natural question whether one can define mechanism-level, expectation-based causal estimands while retaining a strictly design-based logic for identification.
This paper proposes a strict extension of the design-based FPO paradigm by lifting the object of inference from FPO to random potential outcomes (RPO). In our framework, potential outcomes are modelled as random functions of the treatment assignment $z$ and a latent stochastic environment $\omega$, denoted by $\tilde y_i(z,\omega)$; see Assumption (ref). By embedding potential outcomes into a probability space, this extension provides the measure-theoretic machinery required to define causal estimands as expectations taken with respect to the outcome-generating mechanism itself. Identification in the proposed framework continues to rely on the randomisation design and does not require modelling of the outcome process.
Crucially, by treating potential outcomes as random elements, the estimand considered in this paper is no longer a finite-population quantity tied to a particular realised world, but an expectation defined at the level of the latent stochastic environment. In the canonical two-treatment setting, the target of inference is
where the expectation is taken with respect to the probability measure governing the outcome-generating process. Classical FPO estimands are recovered as special cases when the latent space collapses to a singleton, in which case the expectation over $\omega$ is trivial.
In this sense, the proposed extension retains the robustness of classical design-based inference while accommodating mechanism-level, expectation-based causal targets, thereby directly addressing the question of whether such targets can be identified without resorting to parametric or semiparametric outcome models.
The distinction between realised-sample and mechanism-level estimands becomes particularly consequential in settings with stochastic interference. In many modern experiments, a unit's response depends not only on its own treatment assignment, but also on the realised configuration (through $\omega$) of its neighbourhood, which may itself be subject to intrinsic randomness arising from behavioural variation, environmental shocks, or unobserved network dynamics. Under the FPO paradigm, such stochastic spillovers are absorbed into a single realised potential outcome schedule. Consequently, repeating the same experiment under identical assignment rules but a different realisation of the network environment would, by construction, correspond to a different estimand, rendering inference from a single realised sample of limited relevance beyond that particular realisation.
By contrast, the RPO framework treats the network environment as part of the outcome-generating mechanism. This allows causal estimands to be defined as expectations over heterogeneous and stochastic spillover patterns, yielding a stable target that is invariant to the particular realised configuration observed in a single experiment. In such settings, the RPO framework is not merely a modelling choice, but a necessary condition for formulating a well-defined causal estimand.
Although only a single realised world is observed in any given experiment, it represents one draw from this stochastic environment. The central inferential question addressed in this paper is therefore whether, and under what conditions, averaging across units within a single experiment can recover causal estimands defined as expectations over the underlying data-generating mechanism.
Such identification is generically impossible under a design-based logic without additional structure. In particular, absent dependence or ergodicity conditions that allow cross-sectional averaging to substitute for averaging over repeated experiments, expectation-based causal estimands are not identifiable from a single realised experiment. We characterise local dependence structures under which cross-sectional averages exhibit an ergodic property, thereby delineating an identification boundary for ensemble-level causal inference in stochastic environments, rather than merely providing sufficient conditions for estimation.
What is new is not the use of RPO per se, but the ability to conduct valid design-based inference for expectation-based causal estimands defined over a stochastic outcome-generating mechanism using data from a single realised experiment. Notably, the resulting estimator coincides in algebraic form with familiar design-based estimators in the fixed-outcome literature, such as the Horvitz--Thompson estimator horvitz_generalization_1952. A key contribution of this paper is to show that, despite this formal similarity, the same statistic identifies an ensemble-level causal estimand and admits consistent variance estimation and asymptotic normality in the presence of intrinsic outcome-level randomness.
Our approach is closely related to recent functional-analytic perspectives on design-based inference, most notably the Riesz-based framework introduced in harshaw2022riesz. They demonstrate that the Riesz representation theorem provides a unifying language for design-based inference under FPO. Once the object of inference is extended to RPO, however, it is no longer evident whether design-based inference remains feasible without introducing outcome models. The challenge here is not one of technical complexity, but of identification. Irrespective of the technical tools employed, it is not evident whether expectation-based causal estimands defined over a stochastic environment can be recovered using purely design-based arguments from a single realised experiment.
In the RPO framework, the Riesz representer is no longer a fixed weighting function but a stochastic object indexed by the latent environment. As a result, the asymptotic behaviour of the estimator cannot be analysed through combinatorial randomisation alone, but requires control of convergence in operator norm and joint $L^2$ spaces. In this paper, we show that a Riesz-based construction provides one viable approach for conducting design-based inference in this setting. Under appropriate regularity conditions, the resulting estimators are consistent and asymptotically normal, while preserving identification through the randomisation design.
A further contribution of this paper concerns variance estimation. Without a consistent variance estimator, expectation-based estimands are not operationally identifiable from a single realised experiment. In the classical FPO framework, the conservativeness of Neyman-type variance estimators is not a technical artefact but an unavoidable consequence of finite-population reasoning. Within the proposed RPO framework, we show that this limitation can be overcome. By exploiting the ergodic structure induced by local dependence, we develop feasible variance estimators that are consistent for the true sampling variance of the estimator, accounting jointly for randomisation and outcome-level stochasticity. This enables accurate and asymptotically valid inference based on a single experimental realisation, without resorting to overly conservative variance bounds.
Finally, the RPO framework provides a unifying perspective on the long-standing conceptual tension between design-based and model-assisted approaches. While model-assisted methods incorporate outcome randomness through explicit regression models, our framework demonstrates that outcome-level randomness can be fully integrated within a design-based logic, where identification is driven by the randomisation design rather than the outcome model. In this sense, FPO models, RPO models, and model-assisted approaches can be viewed as distinct specialisations within a common inferential structure.
Throughout, we use “random” to refer to outcome-level variation due to latent variables, and “stochastic” to describe modelling frameworks that explicitly incorporate such randomness.
The remainder of the paper is organised as follows. Section (ref) introduces the model and assumptions. Section (ref) establishes large-sample properties under local dependence. Section (ref) presents variance estimation. Section (ref) reports simulation results. Section (ref) concludes.
We begin by considering a setting with $n$ units, where
We refer to this as the stochastic setting because the potential outcomes $\tilde{y}_i(z, \omega)$ are explicitly modelled as random functions of a latent variable $\omega$, introducing outcome-level randomness into the design-based framework. This formulation enables the analysis of causal effects under a stochastic data-generating process and bridges classical design-based inference with mechanism-level causal analysis under stochastic outcome-generating environments.
Although $\omega$ is not unit-specific, its structure may include subcomponents that affect different units differently, thereby allowing for unit-level heterogeneity and flexible dependence structures in the potential outcomes.
The potential outcome for unit $i$ depends on the entire treatment vector $z$, so that interference (or spillover effects) is permitted. Notably, one may also consider the special case of no spillover effects, whereby $y_i(z, x_i, \epsilon_i)=y_i(z_i, x_i, \epsilon_i)$, which corresponds to a version of the Stable Unit Treatment Value Assumption (SUTVA) incorporating consistency, as introduced by rubin1980comment. In the present work we do not impose the no-spillover condition a priori; rather, we allow for general interference and accommodate its presence within the proposed framework.
In practice, it often suffices to allow heterogeneity across units to arise through $(x_i, \epsilon_i)$, making it unnecessary to index the structural function by $i$. In such cases, one may work with a common measurable function $y$, while still inducing heterogeneous potential outcomes via the mappings $\tilde{y}_i(z, \omega) = y(z, x_i, \epsilon_i)$.
The structural mapping (ref) is the primary object of interest in what follows. Following randomisation principles, we assume that
Assumption (ref) ensures that the design distribution remains valid conditional on $\omega$, which is required for constructing the Riesz representer.
Conceptually, conditional on a fixed realisation of $\omega$, each unit's potential outcome function $\tilde{y}_i(z,\omega)$ is deterministic in $z$, so that the only source of randomness in the observed outcomes is the randomised assignment of treatments. This preserves the principle of well-defined potential outcomes, even as outcomes are allowed to vary across realisations of $\omega$.
This formulation is sufficiently general to accommodate both unconfounded and confounded treatment assignment mechanisms, depending on the relationship between $z$ and the latent variable $\omega$. If $z$ and $\omega$ are independent, then $\tilde{y}_i(z^0, \omega)$ and $\tilde{y}_i(z^1, \omega)$ are independent of $z$. More generally, if $z$ and $\omega$ are independent given $x \in X_0$, then $y_i(z^0, x, \epsilon_i)$ and $y_i(z^1, x, \epsilon_i)$ are independent of $z \mid x \in X_0$, for any measurable $X_0 \subset \mathcal{F}_x := \{ A \subseteq x_i(\Omega) : x_i^{-1}(A) \in \mathcal{F}_\omega \}$. In randomised experiments, the design ensures that $z$ is independent of $(x_i, \epsilon_i)$, mitigating confounding. In observational studies, however, it is imperative to adjust for $x_i$ to control for confounding. In either case, this framework models potential outcomes as random rather than fixed.
Let $L^p(\mathcal{Z} \times \Omega)$ denote the space of all measurable functions $u: \mathcal{Z} \times \Omega \to \mathbb{R}$, satisfying $\mathbb{E}[|u(z,\omega)|^p]<\infty$, for some integer $p \geq 1$. For notational convenience, we henceforth write $L^p$ in place of $L^p(\mathcal{Z} \times \Omega)$, whenever the meaning is clear from context. Furthermore, we denote the $p$-norm of $u \in L^p$ by $\|u\|_p := \left( \int_{\mathcal{Z}\times\Omega} \left| u(z, \omega) \right|^p \mu(\mathop{}\!\mathrm{d} z) P(\mathop{}\!\mathrm{d} \omega) \right)^{1/p} = \left( \mathbb{E} \left| u(z, \omega) \right|^p \right)^{1/p}$, where $\mu$ denotes the probability measure for randomisation. Under Assumption (ref), the joint law of $(z,\omega)$ factorises as $\mu(\mathop{}\!\mathrm{d} z) P(\mathop{}\!\mathrm{d} \omega)$, which permits rewriting the $L^p$-norm as an expectation. We write $\|u\| := \|u\|_2$ for brevity if $u \in L^2$.
We model each unit's potential outcome function, via the structural mapping (ref), as an element of a model space $\mathcal{M}_i$, which is a subspace of $L^2$. Formally, we write
Our construction is formulated in a functional setting that separates variation due to randomised treatment assignment from latent outcome-level randomness. Accordingly, the model space $\mathcal{M}_i \subset L^2(\mathcal{Z}\times\Omega)$ accommodates both sources of variation through the product structure of $(z,\omega)$, without imposing additional structural assumptions on either component.
The FPO assumption can be viewed as a special case of our model, corresponding to the degenerate setting where $\Omega = \{\omega_0\}$ for some $\omega_0$. In this case, the latent variable is effectively fixed, and the potential outcome function becomes a deterministic function of $z$, recovering the FPO framework.
We equip $\mathcal{M}_i$ with the inner product
where the second equality follows from the law of iterated expectations. The associated norm is given by $\|u\| = \sqrt{\langle u, u \rangle}$, which coincides with the standard $L^2$-norm.
In the degenerate case where $\Omega = \{\omega_0\}$, the inner product reduces to
which coincides with the inner product employed in harshaw2022riesz. This structure enables us to apply the Riesz representation theorem to construct estimators for treatment effects, even in the presence of interference.
While our model space $\mathcal{M}_i \subset L^2(\mathcal{Z} \times \Omega)$ may not be complete under the inner product defined in (ref), this poses no issue in what follows, as any inner product space admits a Hilbert space completion. Hence, throughout this paper, whenever necessary, we implicitly work with the closure of $\mathcal{M}_i$ in $L^2(\mathcal{Z} \times \Omega)$. Let $\mathcal{M}_i'$ denote the dual space of $\mathcal{M}_i$.
This assumption reflects a standard construction in functional analysis, where a treatment effect is represented as a continuous linear functional on a function space of potential outcomes. Such a formulation enables the use of the Riesz representation theorem, which yields an estimator expressed as an inner product with a representer function. This approach applies to both classical fixed-outcome frameworks and to stochastic settings where potential outcomes depend on latent variables.
Assumption (ref) excludes pointwise evaluation maps. In an $L^2$ model space, evaluation at a single point is not well-defined on equivalence classes and, even when viewed on a representative function space, is not continuous with respect to the $L^2$ norm.
In all cases, $\theta_i$ is a predetermined linear functional on $\mathcal{M}_i$, fixed throughout the analysis and not learned from data. We illustrate admissible choices of $\theta_i$ with the following examples of continuous linear functionals.
The following propositions establish sufficient conditions under which the preceding examples define linear and continuous functionals.
Note that Example (ref) is a special case of Example (ref), obtained by setting $A = \{z^{(1)}\}$ and $B = \{z^{(0)}\}$. The condition $\mu(A) > 0$ and $\mu(B) > 0$ is commonly referred to as the positivity or overlap assumption in the causal inference literature. It ensures that both treatment groups are represented in the randomisation distribution, which guarantees that the corresponding estimand is well-defined and estimable. In our case, the positivity assumption plays a crucial role in verifying that the functional is both linear and continuous.
This result implies that, for the functional to satisfy the necessary continuity condition, the model space must be further restricted to $\mathcal{M}_i \subset L^2(\Omega; H^s(\mathcal{Z})) \subset L^2(\mathcal Z\times\Omega)$. Within this space, functions admit $C^1$ representatives in the $z$-variable almost surely, so that pointwise evaluation of derivatives is well-defined and continuous.
This proposition is deliberately permissive in that the test function $\varphi$ is only required to belong to $L^2(\mathcal Z)$ and is therefore not restricted to indicator functions such as $\mathbf 1_A$ in (ref). More general weights, including discontinuous or piecewise-defined functions, are admissible. In contrast, Dirac-type test functions corresponding to point evaluation are not elements of $L^2(\mathcal Z)$ and hence fall outside the scope of this result. Such pointwise functionals require additional regularity and are instead covered by Proposition (ref).
We emphasise that the appearance of Sobolev-type conditions is not a substantive complication introduced by the Riesz-based approach itself, but a necessary theoretical requirement for ensuring that mechanism-level functionals remain well-defined and continuous when potential outcomes are lifted to the RPO framework.
The corresponding Riesz representer $\psi_i$ provides a representation of the treatment effect functional, but its existence alone does not guarantee feasibility of inference in the stochastic setting. Its role becomes substantive only once combined with conditions under which aggregation across units substitutes for averaging over the latent environment.
In the present setting, the Riesz representer admits a natural operator interpretation. For each unit $i$, the mapping $u \;\mapsto\; \langle u, \psi_i \rangle$ defines a bounded linear functional on $\mathcal{M}_i$, and hence corresponds to a rank-one linear operator induced by $\psi_i$. When potential outcomes are random elements indexed by the latent environment $\omega$, the representer $\psi_i$ itself becomes a stochastic object. As a consequence, the associated Riesz operator is random, and its stability properties under outcome-level randomness are no longer governed solely by combinatorial features of the randomisation design.
In particular, under Assumption (ref), the induced operator remains almost surely bounded. When viewed as an element of the Hilbert space $\mathcal{M}_i$, the representer $\psi_i$ therefore admits a natural notion of stability governed by its $\mathcal{M}_i$-norm. These considerations motivate the operator-norm and $L^2$ modes of convergence employed in the asymptotic analysis that follows, which go beyond purely combinatorial randomisation arguments.
By constructing the Riesz representer \(\psi_i\) corresponding to the treatment effect functional \(\theta_i\), we only consider those treatment effects that are representable in the dual space \(\mathcal{M}_i'\). Consequently, we estimate the treatment effect
by using the corresponding Riesz estimator
as, in practice, we observe only one realisation of the outcome \(\tilde{y}_i(z,\omega)\) for certain $z$ and $\omega$. It follows directly from the linearity of the Riesz representation and the joint expectation over $(z,\omega)$ that the estimator in (ref) is unbiased for (ref). This unbiasedness holds with respect to the combined randomness of the assignment mechanism and the outcome-generating environment, and does not rely on conditioning on a realised potential outcome schedule.
In practice, the Riesz representer can be approximated using standard basis expansions and finite-dimensional projections under the known design; see Appendix (ref) for computational details.
This quantity defines a finite-sample causal estimand under RPO, expressed as a weighted average of unit-level causal functionals. The canonical RPO estimand $\tau_{\mathrm{RPO}}$ in (ref) arises as a special case for particular choices of the weights. Our setup allows for general fixed weighting schemes $\nu_{ni}$, providing greater flexibility in defining estimands.
More importantly, $\tau_n$ integrates over latent randomness $\omega$, defining causal estimands as expectations over both treatment assignment and outcome-level heterogeneity.
The estimand $\tau_n$ does not rely on Assumption (ref); it is mathematically well-defined regardless of whether $z$ and $\omega$ are independent. However, the dependence structure between $z$ and $\omega$ influences the potential outcomes $\tilde{y}_i$ and thus affects the value of the estimand. In the causal inference literature, selecting the right estimand involves specifying a target quantity that captures the scientific question of interest under the assumed data-generating process. Different assumptions about the relationship between $z$ and $\omega$ lead to different estimands, each carrying its own causal or associational interpretation.
Under the FPO framework, the estimand corresponds to a special case of our formulation,
where the latent environment $\omega_0$ is degenerate and the treatment effect is deterministic. In this setting, the average is equally weighted across units. The canonical two-treatment average treatment effect $\tau_{\mathrm{FPO}}$ in (ref) arises as a further special case under uniform weights.
A common choice for the weights is the uniform weighting $\nu_{ni} = 1/n$. More generally, this formulation accommodates subgroup-specific or covariate-adjusted averages, depending on the structure of the weights. We impose the following assumption to control the magnitude of the weights \(\nu_{ni}\).
It guarantees that each weight is bounded by a constant multiple of \(1/n\), uniformly in \(i\), without requiring the individual sequences \(\{n\nu_{ni}\}_{n\ge i}\) to converge. Consequently, no finite subset of indices can receive an asymptotically dominant share of the total weight. Note that Assumption (ref) implies $\sup_{n\to\infty} \sum_{i=1}^n \nu_{ni} < \infty$ but not necessarily one, which is quite general. In the continuous analogue this corresponds to the integrability of the weight function.
Given the unit-level Riesz estimators (ref) and the weights,
Thus we have the unbiasedness of the Riesz estimator in the following.
In the RPO framework, the Riesz representer $\psi_i$ is, in principle, a random element indexed by the latent environment $\omega$. Formally, $\psi_i$ may be viewed as an $\mathcal{M}_i$-valued random variable, mapping each realisation of $\omega$ to a corresponding representer $\psi_i(\cdot,\omega)$. This reflects the fact that both the potential outcomes and the associated inner product structure depend on the outcome-generating mechanism.
In any given experiment, however, only a single realisation $\omega_0$ is observed, just as only a single treatment assignment vector is realised. Consequently, the representer used in computation and estimation is necessarily conditional on $\omega_0$. This conditional construction does not negate the stochastic nature of $\psi_i$, but rather reflects the epistemic limitation inherent in single-experiment inference. The asymptotic theory developed below explicitly accounts for this randomness through joint conditions on $\tilde y_i$ and $\psi_i$.
In this section, we clarify the notion of local independence required in our framework, which is a concept that has been widely adopted in recent literature on causal inference under interference. Informally, local dependence refers to the idea that subsets of random variables are conditionally independent of those outside their respective neighbourhoods. This structure is often represented using a dependency graph, in which vertices correspond to units and edges encode potential statistical dependence. For further discussion of dependency graphs, see for example chen2004normal.
This definition corresponds to Definition 3.1 in ross2011fundamentals and is a weaker version of harshaw2022riesz, Definition 5.1, as it only requires independence of the potential outcome \( \tilde{y}_i \) and the Riesz representer \( \psi_i \) from those outside \( N_i \), rather than full functional independence across the spaces \( \mathcal{M}_i \) and \( \mathcal{M}_j \). We impose independence only on the observable outcome-representer pairs $(\tilde{y}_i, \psi_i)$, not on the underlying model spaces. By modelling potential outcomes as conditionally random functions given a latent variable $\omega$, the framework allows shared model spaces across units while retaining sufficient conditions for variance control and asymptotic normality.
We denote by \( \left\lvert N_i \right\rvert \) the number of units in the dependency neighbourhood of unit \( i \), which we refer to as the neighbourhood size for unit \( i \). Let \( D_n = \max_{1 \leq i \leq n} \left\lvert N_i \right\rvert \) and \( d_n = n^{-1} \sum_{i=1}^{n} \left\lvert N_i \right\rvert \) denote, respectively, the maximum and average neighbourhood sizes across the sample. These quantities will play a central role in the asymptotic analysis to follow.
In practice, the specification of dependency neighbourhoods $N_i$ relies on substantive knowledge about the experimental setting, such as spatial proximity, network structure, or temporal ordering. The framework does not require neighbourhoods to be known exactly, nor does it impose a particular construction. Rather, the assumption formalises the idea that each unit interacts with only a limited subset of others, a feature common to many experimental and quasi-experimental designs.
The following two assumptions restrict the growth rate of the dependency neighbourhoods.
Assumptions (ref) and (ref)(a) regulate the growth of dependency neighbourhoods and play a role analogous to sparsity constraints commonly used in the analysis of dependent data. Assumption (ref)(b) additionally imposes uniform boundedness of the average neighbourhood size. These conditions are standard in the local dependence literature and are stated here explicitly to clarify their role and facilitate reference in our self-contained derivations.
In addition, for a sequence of functions $u_1, \dots, u_n \in L^p$, where $p > 1$, we define the max-$p$ norm as $\left\lVert u \right\rVert_{\max,p}^n := \max_{1 \leq i \leq n} \left\lVert u_i \right\rVert_p$. The sequence $\left\lVert u \right\rVert_{\max, p}^n$ is clearly non-decreasing in $n$, and, as such, is bounded if and only if it converges as $n \to \infty$.
This assumption allows flexibility in the choice of \( r \). It requires uniform boundedness of the \( p \)th moments of \( \tilde{y}_i \) and the \( q \)th moments of \( \psi_i \), thereby ruling out heavy-tailed behaviour.
The following result follows naturally from our variance control framework and the assumptions stated above.
The consistency result in Theorem (ref) provides a formal mathematical interpretation of the ergodic bridge discussed in Section (ref). Under the RPO framework, the target estimand is defined as an expectation over the latent stochastic environment, and identifiability of such expectations from a single realised experiment is not guaranteed a priori under a design-based logic. Theorem (ref) shows that, under local dependence, cross-sectional averaging across units within a single realised experiment converges to the ensemble-level expectation. In this sense, averaging across units substitutes for averaging over repeated experiments. This substitution is not a generic consequence of randomisation alone, but arises from the specific dependence structure imposed on the outcome-generating process. The rate $O_p(n^{-1/2})$ should be interpreted with care. It implies that $\sqrt{n}(\hat{\tau}_n - \tau_n)$ is bounded in probability, but does not guarantee convergence to a fixed non-degenerate limiting distribution. This interpretation applies uniformly to all subsequent large-sample results, including asymptotic normality and variance estimation.
So far, we have shown that, under certain conditions, the aggregate Riesz estimator is consistent. However, these conditions alone are not sufficient for establishing asymptotic normality; an observation that aligns with the fact that consistency, in general, does not imply asymptotic normality.
Let $\sigma_n^2 := \mathbb{V}[\hat{\tau}_n(z,\omega)]$ denote the variance of the aggregate Riesz estimator. The following non-degeneracy assumption imposes a lower bound on the variance of the aggregate Riesz estimator and plays a critical role in ensuring asymptotic normality. When combined with the previous assumptions, it provides sufficient conditions for the aggregate Riesz estimator to satisfy a central limit theorem.
This condition imposes a non-degeneracy requirement on the estimator. It permits the asymptotic variance $\sigma_n^2$ to decay with $n$, but not too rapidly, thus ensuring that the root-$n$ scaled estimator $\sqrt{n}(\hat{\tau}_n - \tau_n)$ retains a non-degenerate variance in the limit.
We require further the following assumption to establish the asymptotic normality result, which concerns the existence of finite fourth moments.
Note that condition $(b)$ imposes a stronger requirement than $(a)$, as it demands uniform boundedness of the fourth moments across all units.
We now present the theorem establishing asymptotic normality of the Riesz estimator in the stochastic setting.
Both Assumptions (ref)\(_{(r=3)}\) and (ref)\(_{(r=4)}\) are required here, as boundedness at $r = 4$ does not imply the same at $r = 3$ over $\mathcal{Z} \times \Omega$, due to the potential unboundedness of the joint space.
Theorem (ref) provides the basis for statistical inference using the aggregate Riesz estimator under outcome-level randomness. In particular, it justifies asymptotic testing and confidence interval construction for the null hypothesis $H_0: \tau_n = 0$, based on the standardised statistic $T_n = \hat{\tau}_n(z, \omega)/\sigma_n$, assuming that both $\hat{\tau}_n(z, \omega)$ and $\sigma_n$ are consistently estimable from a single realisation of $z$ and $\omega$.
The required rate condition on the maximum dependency neighbourhood size, $D_n = o(n^{1/4})$, is in line with that of harshaw2022riesz. However, the proof here proceeds by bounding Wasserstein distance under mixed expectations over both the treatment assignments $z$ and latent heterogeneity $\omega$. Notably, the dependence-degree requirement does not tighten despite the added randomness in potential outcomes.
This result is particularly relevant in empirical settings where potential outcomes are inherently random. For instance, in field experiments involving sensor-based measurements, observed outcomes may be contaminated by random noise due to hardware limitations or environmental fluctuations. Similarly, in panel experiments with repeated units, treatment effects may be estimated via machine learning algorithms that introduce randomness in the form of first-stage fitted values. In both cases, the random component of the outcome arises naturally and cannot be ignored. The asymptotic normality result ensures that inference remains valid under such randomness.
Under local dependence, variance estimation from a single realised experiment becomes feasible once the dependence structure is sufficiently sparse and partially characterised. The test statistic supported by Theorem (ref) requires an estimate of the variance of the aggregate Riesz estimator \( \hat{\tau}_n \), which is generally unknown in finite samples. In the presence of local dependence, such as that arising from network interference, exposure mappings, or latent factors, standard i.i.d.\ variance estimators are invalid, and specialised methods must be employed.
Several approaches have been proposed. Some yield consistent estimators, but only under strong structural assumptions. For example, aronow2017estimating rely on known exposure mappings and strictly positive joint inclusion probabilities; yu_estimating_2022 assume bounded-degree network structures; and liu2014large consider two-stage designs with bounded group sizes.
Of particular relevance to our setting, harshaw2022riesz construct an unbiased estimator for an upper bound on the variance using a tensor-product Riesz representation. More recently, harshaw2024optimizing develop a constrained optimisation framework to tighten such conservative bounds across a wide range of experimental designs. However, neither method formally proves that their bound converges to the true variance, except in special cases. Without sharpness, these bounds may overestimate the variance and lead to overly conservative inference.
In this paper, we take a different approach. Rather than estimating a variance upper bound, we propose a direct estimator for the true asymptotic variance of the aggregate Riesz estimator under local dependence. This estimator is shown to be consistent under conditions that essentially match those required for asymptotic normality. The key requirement is knowledge of the dependency neighbourhood structure, including not only the maximal neighbourhood size \( D_n \), but also which specific units are likely to be dependent. This may seem restrictive, but it reflects a fundamental limitation of single-realisation variance estimation.
Our recommendation is practical. The functional $\theta_i(\cdot)$ is known, and given $\psi_i(\cdot,\omega_0)$ the quantity $\hat\theta_i(z,\omega_0)$ is directly computable from the observed data via (ref). In practice, $\psi_i(\cdot,\omega_0)$ is approximated by $\hat\psi_i(\cdot,\omega_0)$ using standard finite-dimensional projections; see Appendix (ref) for details. We allow for sharp null hypotheses under which $\theta_i(\tilde y_i)$ is specified, and under such nulls, the centred quantity
is well-defined and its realised value $\zeta_i(z,\omega_0)$ is available for variance estimation. We therefore retain only those cross terms corresponding to known or assumed dependent pairs, setting all others to zero, which yields a feasible variance estimator tailored to the assumed dependence neighbourhood and suitable for asymptotically valid inference under local dependence.
Let \( \vec{\zeta}_n = (\zeta_1, \dots, \zeta_n)' \) denote a vector of centred random variables with covariance matrix \( \Sigma_n = \mathbb{E}[\vec{\zeta}_n \vec{\zeta}_n'] \). From a single realisation, we observe only the outer product \( \hat{\Sigma}_n = \vec{\zeta}_n \vec{\zeta}_n' \), and cannot estimate \( \Sigma_n \) consistently without further structure. Even under independence, each \( \vec{\zeta}_n \vec{\zeta}_n' \) typically has variance of order one, so
for some constant \( c > 0 \), where $\left\lVert \cdot \right\rVert_F$ stands for the Frobenius norm for a matrix, indicating that the Frobenius error diverges with \( n \).
Consistent estimation of the full covariance matrix can be recovered only under substantially stronger assumptions. Examples include observing multiple independent replicates of the experiment, which permits standard sample covariance estimation; imposing a time or spatial ordering together with weak stationarity, enabling banded or kernel-smoothed estimators; assuming structural restrictions such as factor models, which allow regularised low-rank-plus-sparse estimation (e.g., POET); or requiring the second-order dependence to vanish asymptotically, for instance when $\|\Sigma_n\|_F \to 0$ as $n \to \infty$. All such approaches rely on additional structure that is unavailable in the single-realisation, design-based setting considered here.
Our proposed approach sidesteps this by targeting only those components of the variance that correspond to known dependent pairs.
Suppose the dependency structure is given by
We define the local-dependence variance estimator
and the population variance
The following theorem shows that, under mild regularity conditions, \( \hat{\sigma}_n^2 \) consistently estimates \( \sigma_n^2 \).
Corollary (ref) establishes that the aggregate Riesz estimator is asymptotically normal under local dependence, provided that the dependency neighbourhoods are sufficiently sparse and the variance estimator \( \hat{\sigma}_n^2 \) is properly constructed.
We next consider variance estimators constructed under progressively weaker information about the dependence structure.
A natural question is whether one can use a reduced dependency structure based solely on correlation. In many experimental settings, it is difficult to explicitly define or justify a full dependency neighbourhood structure \( \mathcal{E}_n \) based on latent interference or unobserved design constraints. However, domain knowledge or structural assumptions may suggest which unit-level estimators \( \zeta_i \) and \( \zeta_j \) are likely to be uncorrelated, even if their full dependence is unknown. This motivates constructing a variance estimator by summing only over pairs believed to exhibit non-negligible second-order dependence, resulting in a conservative correlation-based estimator. Specifically, consider the set
and define the corresponding correlation-based variance estimator
The next theorem shows that the reduced index set still yields a consistent variance estimator under the same local dependence framework.
Although the correlation-based estimator omits many cross-terms in the full variance expression, it remains valid provided that omitted pairs correspond to units whose second-order dependence is negligible. While such assumptions cannot be verified from a single realised experiment, they often align with structural knowledge in applications, including spatial layouts, blocking schemes, or network sparsity.
Although these variance estimators are developed within our stochastic setting, they are equally applicable to the fixed-outcome setting as a special case. Indeed, the classical framework can be recovered by taking the latent space $\Omega$ to be a singleton, in which case the only source of randomness is the treatment assignment. Under this simplification, the local dependence structure reduces to that of the assignment mechanism alone, typically known by design (see Assumption (ref)). Consequently, the variance estimators proposed here apply directly to the FPO framework as a special case and are in fact easier to implement due to the absence of outcome-level variation. This extends the practical utility of our estimators beyond the stochastic setting, providing consistent inference tools for a broad class of design-based estimators.
In the remainder of the paper, we adopt the correlation-based variance estimator \( \hat{\sigma}_{cn}^2 \), while emphasising that the overall methodology still relies on the local dependency assumption, even though it is not explicitly invoked in the variance formula. Crucially, while uncorrelatedness justifies omitting second-order cross terms in the estimator, it does not imply the vanishing of higher-order joint moments across units. Therefore, a corresponding correlation-based neighbourhood structure cannot be meaningfully defined for the full dependency graph.
In practice, this means we continue to assume a local dependency structure satisfying Assumption (ref) with \( d = 1/4 \), while approximating \( \sigma_n^2 \) by summing only over pairs \( (i,j) \) for which \( \zeta_i \) and \( \zeta_j \) are correlated.
Suppose that, in practice, the correlation structure is approximated by a set of index pairs \( \tilde{\mathcal{E}}_n \subset \{1, \dots, n\}^2 \), serving as a practical substitute for (ref). To ensure valid inference, the construction of \( \tilde{\mathcal{E}}_n \) should be conservative, satisfying
Based on this, we define a sparse matrix \( \hat{\Sigma}_n^d \in \mathbb{R}^{n \times n} \) that retains only the entries corresponding to the index set \( \tilde{\mathcal{E}}_n \), representing the estimated dependency structure.
The corresponding variance estimator is then given by
where \( \nu_n = (\nu_{n1}, \dots, \nu_{nn})' \) is the vector of aggregation weights.
By appropriately selecting the conservative set of index pairs, the corresponding variance estimator remains consistent.
The conservativeness of variance estimation should be interpreted differently under the FPO and RPO frameworks.
Under the classical FPO paradigm, variance estimators are necessarily conservative because certain cross-unit covariance terms, such as $S^2_{T-C}$, are fundamentally unobservable. This source of conservativeness is structural and persists even asymptotically. By contrast, under the RPO framework considered here, local dependence renders the variance identifiable at the ensemble level in an asymptotic sense. The variance estimators developed in this section are therefore consistent in the asymptotic sense. Any remaining conservativeness in finite samples arises from deliberate truncation or sparsification of dependence neighbourhoods used to ensure stability and coverage, rather than from intrinsic non-identifiability.
Here, conservativeness refers solely to finite-sample implementation choices and should not be confused with the structural conservativeness of Neyman-type variance bounds, which arises from fundamental non-identifiability.
Note that consistency is preserved even if some independent pairs are mistakenly included in $\tilde{\mathcal{E}}_n$, provided that the total number of dependent pairs, including both genuinely dependent and erroneously included independent pairs, remains within the sparsity condition $D_n = o(n^{1/4})$.
In the special case where all units are believed to be mutually independent, i.e., each \( \zeta_i \) is independent of every other \( \zeta_j \), \( \tilde{\mathcal{E}}_n \) consists only of diagonal pairs \( (i,i) \). Then \( \hat{\Sigma}_n^d \) reduces to a diagonal matrix, and the variance estimator and the true variances simplify to $\tilde{\sigma}_n^2 = \sum_{i=1}^n \nu_{ni}^2 \zeta_i^2$, and $\sigma_n^2 = \sum_{i=1}^n \nu_{ni}^2 \mathbb{E}[\zeta_i^2]$, respectively. This corresponds to the classical form used in standard inference under independence. Corollaries (ref) and (ref) thus recover the familiar asymptotic normality result as a special case of the more general locally dependent setting treated here.
From a functional perspective, the target variance $\sigma_n^2$ can be viewed as a continuous linear functional of second-order moments, such as $\Sigma_n$ or the covariances $\operatorname{Cov}(\zeta_i,\zeta_j)$. Because these moments are not directly observable from a single experimental realisation, the estimator replaces them with observed products $\zeta_i\zeta_j$, yielding the plug-in form in (ref). Consistency follows from structural assumptions on the dependence neighbourhood and boundedness of higher-order terms. Unlike classical plug-in estimators that replace distributional parameters with empirical measures, this construction relies on direct moment substitution guided by a known dependency graph.
Taken together, Theorems (ref)--(ref) and Corollaries (ref)--(ref) show that, although the estimand \( \tau_n \) is defined as an average of individual treatment effects evaluated over the joint probability space \( \mathcal{Z} \times \Omega \), the entire inference procedure, from estimation to variance approximation, can be carried out using just a single realisation of \( z \in \mathcal{Z} \) and \( \omega \in \Omega \). This is made possible by the structure of the Riesz representation framework and the use of carefully constructed estimators that remain valid under local dependence. As the sample size \( n \) grows, one draw from the treatment and outcome distribution contains sufficient information for reliable estimation and inference, provided that the underlying conditions hold. In this way, the results provide a rigorous justification for making population-level claims based on one observed experiment.
To illustrate the finite-sample inferential behaviour of the proposed framework, we conduct a simulation study focused on empirical rejection behaviour rather than estimator comparison. Because the Riesz estimator takes the same algebraic form under both FPO and RPO, there is no meaningful notion of competing estimators to compare. More fundamentally, under the FPO framework, inference is conditional on a realised potential outcome schedule and does not target consistent estimation of the ensemble-level estimand considered here. Accordingly, the simulations evaluate whether design-based inference based on a single realised experiment delivers correct coverage for the expectation-based estimand induced by outcome-level randomness. Local dependence is introduced through latent block-level effects, allowing us to examine how inferential validity depends on the strength of dependence.
We consider two data-generating mechanisms. Throughout, treatment assignment is generated independently as $z_i\sim\mathrm{Bernoulli}(0.5)$ for each unit $i\in\{1,\dots,n\}$. Unit-level covariates and parameters are drawn independently according to
Inference is conducted under a sharp-null formulation that specifies unit-level contrasts. In the size experiments, the centring term is taken to be the true unit-level contrast $\theta_i(\tilde y_i)$. In the power experiments, we instead impose the null specification $\theta_i(\tilde y_i)\equiv 0$ and evaluate rejection frequencies under this centring. This choice is made for implementation convenience. Since the test statistics depend only on the difference between the centring term and the true contrast, reversing the null and alternative would lead to identical rejection behaviour. The present setup allows size and power to be assessed with only a change in the centring specification, without otherwise modifying the simulation code.
Local dependence is introduced by partitioning the sample into blocks, where units in the same block share a common latent shock. Specifically, the number of blocks is $B_n = \lfloor n^{1 - d} \rfloor$ for a fixed parameter $d \in \{0, 0.1, 0.2, 0.25, 0.3\}$. The error for unit $i$ is
where $b(i)$ denotes the block index to which unit $i$ belongs, $\eta_{b(i)} \sim \mathcal{N}(0, 1)$ is the shared block-level shock, $\nu_i \sim \mathcal{N}(0,1)$ is an idiosyncratic noise term, and $\gamma_i \in \{-1,1\}$ is a random sign (independent Rademacher) to allow both positive and negative correlations between units in the same block. This construction yields a dependency graph $\mathcal{E}_n$ with within-block dependence and across-block independence.
In the first case, potential outcomes are generated via
where $\delta_i$ are drawn i.i.d.\ from $\mathcal{N}(0,1)$, $\varepsilon_i$ is an outcome-level noise term exhibiting local dependence described below. In this specification, the individual treatment effect satisfies $\theta_i(\tilde y_i) = \beta_i$, so that the finite-population FPO and RPO estimands coincide. This design therefore serves as a baseline setting in which outcome-level randomness does not alter the causal target, while local dependence affects inference through the variance structure.
In the second case, we consider the network interference model described in Example (ref). Local dependence is introduced through the same block structure as in the baseline design. A random graph $G_n(\omega)$ is generated independently within each block. For any pair of distinct units $i$ and $j$ belonging to the same block, an undirected edge between $i$ and $j$ is formed independently with probability $p_{\mathrm{edge}}=0.1$. No edges are formed between units belonging to different blocks. $A_{ij}(\omega)=1$ if an edge is present between units $i$ and $j$, and $A_{ij}(\omega)=0$ otherwise. By construction, $A_{ii}(\omega)=0$ for all $i$.
Given a treatment assignment vector $z\in\{0,1\}^n$, we obtain the realised neighbourhood $N_i(\omega)$ and then compute the realised exposure $e_i(z,\omega)$ according to (ref). If unit $i$ has no neighbours, the exposure is set to zero. Potential outcomes are generated according to (ref), with $\gamma_i = 0.5$ controlling the strength of spillover effects.
In this network setting, the FPO estimand is
whereas the RPO one is
where the final equality follows from Proposition (ref) in Appendix (ref) and $m_i$ is the block size.
It is worth emphasising that, $\tau_n^{\mathrm{FPO}}$ is tied to the particular realised network and outcome realisation observed in the experiment. By contrast, $\tau_n^{\mathrm{RPO}}$ is defined at the level of the underlying stochastic environment and therefore remains invariant across repeated experimental realisations and across samples drawn from the same probability space, provided the conditions of our asymptotic results are satisfied.
We estimate the average treatment effect using the aggregate Riesz estimator. Proposition (ref) establishes that, under very weak structural assumptions on the outcome process, the two-treatment Bernoulli design admits a known Riesz representer. Accordingly, we work in this setting, in which the representer takes the familiar Horvitz--Thompson form
where $\mu_1=\mathbb{E}[z_i]=0.5$ and $\mu_0=1-\mu_1$; see, e.g., neyman1923, ImbensRubin2015, and aronow2017estimating. Accordingly, the simulations do not involve auxiliary estimation of $\psi_i$, allowing us to focus on the effects of local dependence and on the behaviour of the proposed variance estimators.
The estimator is $\hat{\tau}_n = \frac{1}{n} \sum_{i=1}^n Y_i \psi_i(z_i, \omega_0)$, and the variance estimator follows (ref) with $\zeta_i = Y_i \psi_i(z_i) - \theta_i(\tilde{y}_i)$ and $\nu_{ni} = 1/n$. That is, we sum pairwise products of residuals over all unit pairs within the same block.
In this simulation setting, the variance estimators \( \hat{\sigma}_{cn}^2 \) and \( \tilde{\sigma}_{n}^2 \) coincide with the local-dependence estimator \( \hat{\sigma}_n^2 \). This is because the dependency structure is fully captured by shared block-level shocks, and all non-zero cross-moments correspond to block-sharing units. Therefore, all simulation results reflect the performance of both estimators. In general, however, the correlation-based estimator yields improved efficiency in settings where weakly dependent pairs can be excluded without loss of accuracy.
We consider sample sizes $n \in \{100, 200, 500, 1000\}$ and dependency growth parameters $d \in \{0.0, 0.1, 0.2, 0.25, 0.3\}$, where the largest value $d = 0.3$ lies outside the range permitted by Assumption (ref) $(d \le 1/4)$ in Theorem (ref). For each $(n,d)$ pair, we run $2000$ independent replications. For each replication, we compute the estimator $\hat{\tau}_n$ and the standard error $\hat{\sigma}_n$, and conduct inference under two specifications of the centring term, the correctly specified sharp null $\theta_i(\tilde y_i)$ and the misspecified null $\theta_i(\tilde y_i)\equiv 0$. We record the following metrics:
Figures (ref) and (ref) report the empirical coverage of nominal 95% confidence intervals across different dependency growth rates $d$ and sample sizes $n$.
Figure (ref) corresponds to the baseline design in which the FPO and RPO estimands coincide. Across all configurations, empirical coverage is slightly above the nominal $0.95$, indicating mild finite-sample conservativeness. For small samples ($n\in \{100,200\}$), coverage shows no clear monotone pattern in $d$, consistent with higher sampling variability. For larger $n$, especially $n=1000$, coverage increases gradually with $d$, reflecting the growing conservativeness of the variance estimator as dependence neighbourhoods expand. Importantly, within the theoretical regime $d\le 1/4$ (Theorem (ref)), coverage remains close to nominal and uniformly satisfactory.
Figure (ref) reports coverage under network interference. Coverage is uniformly above 95%, indicating stronger conservativeness than in the baseline case. This persists even at $d=0$, suggesting that interference through the exposure mapping $e_i(z,\omega)$ affects calibration beyond what is captured by neighbourhood growth alone. When dependence exceeds the theoretical range ($d=0.3$), coverage approaches 97%, consistent with a gradual shift towards conservativeness as assumptions are violated. A further feature is that, at $d=0$, coverage declines from $n=500$ to $n=1000$; this non-monotonicity appears specific to the spillover structure and does not indicate failure of the variance estimator. Overall, coverage remains stable and well controlled, supporting the feasibility of design-based inference for ensemble-level causal estimands under network interference.
Tables (ref) and (ref) report empirical rejection frequencies at the 1%, 5%, and 10% nominal levels for the size and power experiments, across all combinations of sample size $n$ and dependency growth rate $d$, for the two data-generating mechanisms described above.
In the baseline setting, where the FPO and RPO estimands coincide, empirical rejection frequencies in the size experiment remain close to their nominal levels across all sample sizes and dependency regimes. At the 5% level, rejection rates typically lie between approximately 4% and 5%, with similar behaviour observed at the 1% and 10% levels, and no systematic drift as $n$ increases. These results indicate that the proposed procedure maintains appropriate size control under local dependence when the conditions of Corollary (ref) are satisfied.
Across most configurations, rejection frequencies are slightly below nominal levels, indicating mild finite-sample conservativeness. This behaviour is consistent with the coverage results in Figure (ref) and suggests that conservativeness is driven primarily by variance estimation rather than estimator bias. When the dependency growth rate exceeds the theoretical threshold in Assumption (ref), mild over-rejection emerges, most notably at $d=0.3$, consistent with the requirement that dependence neighbourhoods grow sufficiently slowly.
In the power experiment, rejection frequencies increase monotonically with sample size. For moderate $n$, power remains substantial even under stronger dependence, and for $n\ge 500$ rejection rates are essentially one across all significance levels, confirming that the procedure achieves high power without compromising size control within the theoretical regime.
In the network-interference setting, where the estimand targets an ensemble-level RPO quantity rather than a realised finite-population schedule, size performance remains stable and close to nominal levels across configurations. As in the baseline design, slight under-rejection is observed, mirroring the coverage results in Figure (ref). Notably, stochastic spillovers and outcome-level randomness do not introduce additional size distortions relative to the baseline case.
Power properties remain strong under network interference. Rejection frequencies increase rapidly with $n$, and for $n\ge 200$ power is close to one at conventional significance levels. As in the baseline design, power decreases modestly with stronger dependence for small $n$, reflecting information loss due to local dependence. Overall, the results demonstrate that the proposed procedure retains high power while maintaining accurate size control when targeting ensemble-level causal estimands under network interference.
Overall, the simulation results are fully consistent with the theoretical analysis. Within the dependency regimes covered by the theory, the aggregate Riesz estimator delivers reliable inference under local dependence, with finite-sample behaviour closely aligned with the asymptotic results.
Although the simulation designs are based on parametric data-generating mechanisms, the proposed estimation and testing procedures do not exploit any parametric structure. All inference is conducted in a fully design-based manner, relying solely on the treatment randomisation and the dependence structure. Crucially, the causal target is defined as an expectation over RPO rather than as a fixed realised schedule. The simulations therefore illustrate that valid inference for ensemble-level causal estimands can be achieved from a single realised experiment under local dependence, without introducing outcome models.
The mild conservativeness observed in finite-sample coverage arises from the variance estimation procedure under dependence and diminishes as sample size increases. This behaviour stands in sharp contrast to the structural conservativeness of classical Neyman-type variance bounds in the FPO framework, which is intrinsic to finite-population reasoning and does not vanish asymptotically.
This paper extends the classical design-based framework for causal inference to settings in which potential outcomes are random. By allowing outcomes to depend on a latent stochastic environment, we shift the object of inference from realised potential outcome schedules to expectation-based causal estimands defined at the level of the outcome-generating mechanism. The FPO framework is recovered as a degenerate special case, while identification continues to rely exclusively on randomised treatment assignment without introducing outcome models or parametric assumptions.
Within this stochastic setting, the central inferential challenge is whether expectation-based causal estimands are identifiable and estimable from a single realised experiment. We show that under suitable local dependence structures, cross-sectional averaging across units exhibits an ergodic property that bridges a single realised world and ensemble-level causal quantities. Absent such dependence restrictions, expectation-based causal estimands are fundamentally not identifiable from a single realised experiment under a design-based logic. These conditions characterise the identification boundary for expectation-based estimands in stochastic environments and go beyond the scope of classical randomisation-based arguments. From a broader perspective, this identification problem can be viewed as a generalisation problem in which inference moves from a single realised experimental world to causal effects defined at the level of the underlying outcome-generating mechanism.
The paper also develops feasible variance estimators that are consistent for the true sampling variance under local dependence. In contrast to the FPO framework, where conservative variance bounds are unavoidable, the proposed approach enables consistent uncertainty quantification for expectation-based estimands from a single experiment. Although motivated by RPO, the variance estimators also apply to the classical fixed-outcome setting as a special case, providing a practical alternative to conservative bounds when dependence structures are known.
The main theoretical results establish consistency and asymptotic normality of the standardised aggregate Riesz estimator under mild moment and neighbourhood growth conditions. Inference is conducted directly for the finite-sample target $\tau_n$, which, while computed from a single realised experiment, serves as a bridge to the ensemble-level causal estimand defined over the underlying stochastic environment. Under local dependence, averaging across units plays the role traditionally assigned to averaging across repeated experiments.
Simulation results support the theoretical analysis. Across a range of sample sizes and dependence regimes, confidence intervals exhibit accurate or slightly conservative coverage, and rejection frequencies remain close to nominal levels within the theoretical range of dependence growth. When the dependence neighbourhood grows beyond this range, calibration deteriorates gradually, in line with the limits imposed by the assumptions. More broadly, the results illustrate a classical theme in statistics, namely that familiar design-based statistics can support inference for more general causal objects once the target of inference is lifted to an ensemble-level estimand and variance feasibility is established, without introducing outcome models.
Taken together, these results demonstrate that design-based inference remains feasible for expectation-based causal estimands in stochastic outcome-generating environments. The Riesz representation provides a principled and flexible construction for such inference, particularly in accommodating the functional-analytic complexity introduced by local dependence and stochastic outcome variation. The core contribution, however, lies in identifying precise conditions under which inference from a single realised experiment is possible, and in establishing consistent uncertainty quantification under dependence. Possible extensions include non-smooth functionals, dynamic treatments, and observational designs, while preserving the central principle that valid inference should rest on the randomisation and dependence structure rather than on outcome models.
\nocite{chen2007sieve, chen_pouzo2012moment, van2000asymptotic} \nocite{adams2003sobolev} \nocite{billingsley1999convergence} \nocite{ross2011fundamentals} \nocite{riesz1907} \nocite{Chernozhukov02122025}
I am grateful to Fredrik Sävje for his valuable discussions and support.
No external or empirical datasets were used in this study. All data arise from simulations conducted by the author. The complete codebase, including replication scripts and plotting routines, is publicly available at \href{https://github.com/yukai-yang/RieszRE_Experiments}{\url{github.com/yukai-yang/RieszRE_Experiments}} under the MIT license.