calculating variance of sum of predicted costs from two-part model

I have a question regarding the calculation of the variance of the sum of predicted health care costs using a two-part regression model. Details are below, but it boils down to how to calculate the covariance for a random variable which is itself the product of two random variables.

The two-part regression model consists of a logistic regression component for estimating the probability $p$ of cost > 0 and an OLS component predicting the cost $y$ given $y > 0$. The predicted cost $C_i$ for a given member $i$ is then the product $C_i=p_iy_i$. We can assume that $p_i$ and $y_i$ are independent.

I'm calculating the variance of the sum of costs $S=\Sigma_{i=1}^mC_i$ for members $i=1$ to $m$ belonging to a particular provider, assuming member level predicted costs are not independent given they're sharing the same two-part regression model.

So running through the math:

$Var(S)=\Sigma_{i=1}^m Var(C_i)+\Sigma\Sigma_{i\neq j}^mCov(C_i,C_j)$

where $Var(C_i)$ is:

$Var(C_i) = E(C_i^2) - [E(C_i)]^2$

$=E(p_i^2y_i^2)-[E(p_i)E(y_i)]^2=E(p_i^2)E(y_i^2)-[E(p_i)E(y_i)]^2$

$=(Var(p_i)+[E(p_i)]^2)(Var(y_i)+[E(y_i)]^2)-[E(p_i)E(y_i)]^2$

and $Cov(C_i,C_j)$: is

$Cov(C_i,C_j)=E(C_iC_j)-E(C_i)E(C_j)$.

My first inclination was to use a multivariate Taylor expansion to approximate the term $E(C_iC_j)$.

However I run into trouble with a first-order expansion:

$Z=g(C_i,C_j)\approx g(\mu_{C_i},\mu_{C_j}) + (C_i-\mu_{C_i})\frac{\delta g(\mu_{C_i},\mu_{C_j})}{\delta C_i} + (C_j-\mu_{C_j})\frac{\delta g(\mu_{C_i},\mu_{C_j})}{\delta C_j}$

since $E(Z)\approx g(\mu_{C_i},\mu_{C_j})=E(C_i)E(C_j)$ which means $Cov(C_i,C_j)\approx 0$, which can't be right.

Is there another method for approximating the covariance of $C_i$ and $C_j$? Or did I make a mistake in my calculations? Thanks!