I would like to approximate $\mathbb{E}[e^y]$ for $y=\sum_i^n c_i x_i$, and $x_i\sim \text{Bernoulli}(p_i)$ ($p_i$'s distinct) using $p_i$ and $\text{cov}(x_i,x_j)$. The $c_i\in \mathbb{R}$ are fixed, and the $x_i$ are not independent.
Currently, I am estimating $\mathbb{E}[e^y]$ by assuming the $x_i$ are independent. That is, I take $\mathbb{E}[e^y]\approx \prod_i[ (1-p_i)+p_i e^{c_i}]$, which only requires me to know $p_i=\mathbb{E}[x_i]$.
However, in addition to knowing $p_i$, I also know $\text{cov}(x_i,x_j)$, but I am not using this information. I am searching for a way to improve my estimation of $\mathbb{E}[e^y]$ by making use of the covariances.
I have tried a few things. Each has problems which I will attempt to illustrate with a toy example. First a bit more notation: Let $\mu=\mathbb{E}[y]$, $\sigma^2=\text{var}[y]$, and take $\hat{y}=\sum_i^n c_i\hat{x}_i$ with $\hat{x}_i\sim \text{Bernoulli}(p_i)$ all independent.
Now a toy example. Suppose $y=7x_1+5x_2-8x_3$, $P(x_1=1,x_2=1,x_3=0)=P(x_1=0,x_2=1,x_3=1)=0.2$, and all other assignments of $x_1,x_2,x_3$ occur with probability $0.1$. Then, $\text{cov}(x_1,x_3)=-0.05$, $\text{cov}(x_1,x_2)=\text{cov}(x_2,x_3)=0$, and $\mathbb{E}[e^y]=32,681.07$.
Methods of approximation:
1) Ignore the dependencies: $\mathbb{E}[e^y]\approx \mathbb{E}[\text{exp}(\hat{y})]$. This approach only uses first order information ($p_i$). On the toy example we get $\mathbb{E}[\text{exp}(\hat{y})]=24,553.48$.
2) Taylor series approximation: $\mathbb{E}[e^y]\approx e^\mu(1+\tfrac12 \sigma^2)$. This approach gives terrible approximations for large $c$'s. On the toy example we get $e^\mu(1+\tfrac12 \sigma^2)=7.697877$.
3) Log-normal approximation: $\mathbb{E}[e^y]\approx e^{\mu+\tfrac12 \sigma}$. This approach overshoots very badly which is problematic for my application. On the toy example $e^{\mu+\tfrac12 \sigma}=165,585,792$. But note that $\mathbb{E}[e^y]$ is bounded between $\text{exp}({\sum_i \mathbb{1}_{\{c_i<0\}}c_i}) \le \mathbb{E}[e^y] \le \text{exp}({\sum_i \mathbb{1}_{\{c_i>0\}}c_i})=162,754.8$.
4) Treat $x_i$ as independent, but correct the second moment: Note that $\mathbb{E}[e^y]=\sum_{k=0}^\infty \frac{1}{k!}\mathbb{E}[y^k]$, then we might consider $\mathbb{E}[e^y]\approx \mathbb{E}[\text{exp}(\hat{y})]+\tfrac12 \mathbb{E}[y^2]-\tfrac12 \mathbb{E}[\hat{y}^2]$. I'm not sure this is well principled, and I empirically observe very small corrections when the coefficients $c_i$ are large (as in the toy example). In this case, this approximation gives a value of $24,556.28$.
More details:
The above problem arises from the following. Imagine an $n\times m$ matrix $R$ of indicator vectors, i.e. each row is an independent random vector in $\{0,1\}^m$ that takes at most a single 1 (all other positions being zero), the probabilities of which are known. Then $y=\tfrac12 \text{tr}(R^T C R)$ where $C$ is a symmetric matrix of constants with zeros on the main diagonal. In the above formulation $x_i$ is an off diagonal element of the matrix $RR^T$. I presented the problem in that way because I thought it makes the problem more clear. I've given these additional details in case they are important in some way.
For what it's worth, in my problems $n$ is usually $O(10^1)-O(10^2)$ and some of the $c_i$ may be large in magnitude (e.g. $\ge O(10^1)$).