For a gaussian random vector $X = (X_1, \dots, X_n) \in {\bf R}^n$ with parameters $\Sigma$ and $\mu$, meaning that
$\displaystyle \mathop{\bf P}( (X_1,\dots,X_n) \in S) = \frac{1}{(2\pi)^{n/2} (\det \Sigma)^{1/2}} \int_S e^{-\frac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu)}\ dx$
for all Borel sets ${S \subset {\bf R}^n}$ (we identify elements of ${{\bf R}^n}$ with column vectors), where $\Sigma = (\sigma_{ij})_{1 \leq i,j\leq n}$ is an $n \times n$ positive definite real symmetric matrix, $\mu \in {\bf R}^n$, show that $\hbox{Cov}(X_i, X_j) = \sigma_{ij}$.
Question: If $f: {\bf R}^n \to {\bf R}$ is given by $f(x) = \pi_i(x-\mu)\pi_j(x-\mu)$, where $\pi$ is the coordinate projection function, then by change of variables formula $\hbox{Cov}(X_i, X_j) = \int_{{\bf R}^n} \pi_i(x-\mu)\pi_j(x-\mu)\ d\mu_X(x) = \int_{{\bf R}^n} (x_i-\mu_i)(x_j-\mu_j)\ d\mu_X(x_1, \dots, x_n)$, where $\mu_X$ is the law of $X$. I have some difficulties directly showing this evaluate to $\sigma_{ij}$, in particular when $\mu_X$ is not assumed to be a product measure which split nicely. Any suggestion or hint on things to try will be appreciated.