Estimating the uncertainty on the difference of two random variables Let $(X_0, X_1)$ a pair of (non-independent) random variables and $Z\in\{0,1\}$ a Bernoulli random variable, independent from $(X_0,X_1)$. Let $\Delta = X_1 - X_0$. We're interested in the distribution of $\Delta$.
Consider $(X_{0,i}, X_{1,i}, Z_i)_{1\leq i\leq n}$ i.i.d. copies of these random variables and $\Delta_i = X_{1,i} - X_{0,i}$.
Now assume that we only observe $(Z_i,X_{Z_i,i})_{1\leq i\leq n}$. In other words, for each $i$, we observe either $X_{0,i}$ or $X_{1,i}$, but not both. We know which one we have observed but we don't have access to the other.
Based on this information, we can estimate the expected value of $\Delta$ by computing
$$\frac{1}{\#\{i:Z_i=1\}}\sum_{i:Z_i=1}X_{1,i} - \frac{1}{\#\{i:Z_i=0\}}\sum_{i:Z_i=0}X_{0,i}.$$
However, I want to obtain a measure of dispersion of the distribution of $\Delta$, e.g. by estimating its variance. I'm not sure how to do this, and I'm not even sure this is possible because we have no way to assess the dependence relationship between $X_0$ and $X_1$, as we never observe both of them. Any idea or related literature?
 A: The comment thread suggests that one insightful way to frame this problem is to distinguish the underlying random variables from what is observed.  Let's begin, then, by doing exactly that.
We will model the sample by means of the bivariate random variable $(X_0,X_1)$ in the usual way, by supposing the observations are associated with a sequence of iid versions of this variable $(X_{0i},X_{1i}),$ $i=1,2,\ldots,n.$  The observations are determined independently by parallel iid sequence of a third variable $Z$.  Thus, what is observed is a sequence
$$Y_i = (1-Z_i) X_{0i} + Z_iX_{1i}.$$
The equations
$$Z_i Y_i = Z_i X_{1i}$$
and
$$(1-Z_i)Y_i = (1-Z_i) X_{0i}$$
enable us to express the estimator of the mean of $X_0-X_1$ in the question as
$$\widehat{E[\Delta]} = \frac{\sum_{i=1}^n (1-Z_i) Y_i}{\sum_{i=1}^n (1-Z_i)} - \frac{\sum_{i=1}^n Z_i Y_i}{\sum_{i=1}^n Z_i}.$$
(In order for this to be well defined, we must understand either of these fractions to equal zero whenever its denominator is zero.)
Let's begin by computing the expected value of this estimator conditional on the $Z_i.$  Writing $\mathbf{Z}=(Z_1,\ldots,Z_n),$
$$E\left[\widehat{E[\Delta]}\mid \mathbf{Z}\right] = \frac{\sum_{i=1}^n (1-Z_i) \mu_0}{\sum_{i=1}^n (1-Z_i)} - \frac{\sum_{i=1}^n Z_i \mu_1}{\sum_{i=1}^n Z_i}=\mu_0-\mu_1,$$
as we would hope.  Similarly, because all the $X_{0i}$ are uncorrelated with all the $X_{1j},$ its conditional variance is
$$\begin{aligned}
\operatorname{Var}\left(\widehat{E[\Delta]}\mid \mathbf{Z}\right) &= \frac{\sum_{i=1}^n (1-Z_i) \sigma_0^2}{\left(\sum_{i=1}^n (1-Z_i)\right)^2} + \frac{\sum_{i=1}^n Z_i\sigma_1^2}{\left(\sum_{i=1}^n Z_i\right)^2} \\
&= \sigma_0^2\frac{1}{\sum_{i=1}^n (1-Z_i)} + \sigma_1^2\frac{1}{\sum_{i=1}^n Z_i}
\end{aligned} \tag{*}$$
(because $Z_i^2 = Z_i$ and $(1-Z_i)^2 = 1-Z_i$).
The Law of Total Variance (which follows, upon applying a little algebra, from the standard formulas for variance) asserts
$$\begin{aligned}
\operatorname{Var}\left(\widehat{E[\Delta]}\right) = E\left[\operatorname{Var}(\widehat{E[\Delta]}\mid \mathbf{Z})\right] + \operatorname{Var}\left(E\left[\widehat{E[\Delta]}\mid \mathbf{Z}\right]\right)
\end{aligned}.$$
The second term is the variance of constant and so drops out.  Evaluating the first term requires us to compute the expectations of the fractions in $(*).$  There is no closed form, but the values can be exactly calculated for small to medium $n$ because $\sum Z_i$ has a Binomial distribution with parameters $n,p$ and $\sum (1-Z_i)$ has a truncated Binomial distribution with parameters $n,1-p.$  Thus, bearing in mind the fraction is treated as $0$ when all the $Z_i$ are $0,$
$$E\left[\frac{1}{\sum_{i=1}^n Z_i}\right] = 0 + \sum_{k=1}^n \left(\frac{1}{k}\right)\binom{n}{k}p^k(1-p)^k = \int_0^1 \frac{(1 - (1-x)p)^n - (1-p)^n}{x}\,\mathrm{d}x.$$
The integral is useful for computations when $np$ is small to medium.  Otherwise, the Normal approximation to the Binomial distribution ($np$ and $n(1-p)$ are large) gives
$$E\left[\frac{1}{\sum_{i=1}^n Z_i}\right] \approx \int_1^\infty \frac{1}{x}\phi\left(\frac{x-\mu}{\sigma}\right)\,\mathrm{d}x$$
where $\phi$ is the standard Normal density.
By analyzing either integral it's easy to see that to first order this expectation equals the reciprocal of the expectation of the denominator, or $1/(np).$
The calculation involving the $1-Z_i$ is identical but with $p$ replaced by $1-p.$  Thus, to first order, the following formula gives a reasonable answer to the question:

$$\operatorname{Var}\left(\widehat{E[\Delta]}\right) \approx \frac{\sigma_0^2}{(1-p)n} + \frac{\sigma_1^2}{pn}.$$

If necessary (for small $n$ or extreme values of $p$) it can be improved using either of the two exact formulas (binomial sum and integral).
You can, of course, estimate $\sigma_0^2$ and $\sigma_1^2$ separately from the data for which $Z_i=0$ and $Z_i=1,$ respectively.
