Say I have two random variables $Z^{p}_{i} \mid Y^{p}_{i} \sim \text{Poisson} \left(t^{p}_{i}Y^{p}_{i}\right)$ and $Z^{q}_{i} \mid Y^{q}_{i} \sim \text{Poisson} \left(t^{q}_{i}Y^{q}_{i}\right)$ for all sites $i$ belonging to domain $\mathcal{D}$. Both $Y^{p}_{i}$ and $Y^{q}_{i}$ are positive random fields with means $m^{p}$ and $m^{q}$ and, variances $\sigma_{p}^{2}$ and $\sigma_{q}^{2}$ , respectively.
We know that for $Z^{p}_{i}$ (and identically for $Z^{q}_{i}$) it holds that:
$E\left[Z^{p}_{i} \mid Y^{p}_{i}\right]=t^{p}_{i} Y^{p}_{i}$
$\operatorname{Var}\left[Z^{p}_{i} \mid Y^{p}_{i}\right]=t^{p}_{i} Y^{p}_{i}$
$E\left[Z^{p}_{i}Z^{p}_{j} \mid Y\right]=\operatorname{Cov}[Z^{p}_{i}Z^{p}_{j} \mid Y]+ E\left[Z^{p}_{i} \mid Y^{p}_{i}\right]E\left[Z^{p}_{j} \mid Y^{p}_{j}\right] = \delta_{ij}t^{p}_{i}Y^{p}_{i} + t^{p}_{i}Y^{p}_{i}t^{p}_{j}Y^{p}_{j},\\ \text{where } \delta_{ij} \text{ is the Kronecker delta} \; \delta_{ij}=\left\{\begin{array}{lll} 0 & \text{if} & i \neq j \\ 1 & \text{if} & i=j \end{array}\right.$
When $i = j$, we get the following conditional expectation as we assume conditional independence of observations at diferent sites $i$:
$E\left[Z^{p}_{i}Z^{p}_{j} \mid Y\right]=\operatorname{Cov}[Z^{p}_{i}Z^{p}_{j} \mid Y]+ E\left[Z^{p}_{i} \mid Y^{p}_{i}\right]E\left[Z^{p}_{j} \mid Y^{p}_{j}\right] = \delta_{ii}\operatorname{Var}[Z^{p}_{i} \mid Y^{p}_{i}] + E\left[Z^{p}_{i} \mid Y^{p}_{i}\right]E\left[Z^{p}_{j} \mid Y^{p}_{j}\right] = t^{p}_{i}Y^{p}_{i} + t^{p}_{i}Y^{p}_{i}t_{j}Y^{p}_{j}$
Naturally the conditional expectation of the two processes $Z^{p}_{i}$ and $Z^{q}_{j}$ at the same site $i$ is,
$E\left[Z^{p}_{i}Z^{q}_{j} \mid Y\right]=\operatorname{Cov}[Z^{p}_{i}Z^{q}_{j} \mid Y]+ E\left[Z^{p}_{i} \mid Y^{p}_{i}\right]E\left[Z^{q}_{j} \mid Y^{q}_{j}\right] = \delta^{pq}_{ii}\operatorname{Cov}[Z^{p}_{i}Z^{q}_{i} \mid Y] + E\left[Z^{p}_{i} \mid Y^{p}_{i}\right]E\left[Z^{q}_{i} \mid Y^{q}_{i}\right] = \operatorname{Cov}[Z^{p}_{i}Z^{q}_{i} \mid Y] + t^{p}_{i}Y^{p}_{i}t^{q}_{i}Y^{q}_{i}$.
My objective is to calculate the conditional covariance of $Z^{p}_{i}$ and $Z^{q}_{j}$ given $Y$.
Can someone suggest to me how I can proceed to obtain $\operatorname{Cov}[Z^{p}_{i}Z^{q}_{j} \mid Y]$ in terms of $Y$?