If $A$ is a $p \times p$ symmetric matrix of constants, and $y \sim N_p(\mu, \Sigma)$, justify the following result: $$ \newcommand{\Cov}{\operatorname{Cov}} \Cov(y, y^\top\!Ay) = 2\Sigma A\mu $$
I know that I should expand this equation into $$ \newcommand{\E}{\operatorname{E}} \Cov(y, y^\top\!Ay) = \E(yy^\top\!Ay) - \E(y) \E(y^\top\!Ay) $$ where $\E(y) = \mu$ and $$ \E(y^\top\!Ay) = \E(y^\top)\E(Ay) = \mu^\top\!A \E(y) = \mu^\top\!A\mu $$
I need some guidance with how to go about with this term: $\E(yy^\top\!Ay)$.
How does the $\Sigma$ go into the equation?
I attempted the question but I arrived at the answer that the covariance is zero.
With these sorts of questions, it is often useful to write $y = \mu + \Sigma^{1/2}z$, where $z \sim N_p(0, 1)$ is a vector of iid standard Gaussians.
Now, we can compute that \begin{align*} \operatorname{Cov}(y, y^{\mathrm{T}} Ay) &= \operatorname{Cov}(\mu + \Sigma^{1/2}z, (\mu + \Sigma^{1/2}z)^{\mathrm{T}} A(\mu + \Sigma^{1/2}z)) \\ &= \operatorname{Cov}(\Sigma^{1/2}z, z^{\mathrm{T}}\Sigma^{1/2}A\Sigma^{1/2}z + 2 \mu^{\mathrm{T}}A\Sigma^{1/2}z) \\ &= \operatorname{Cov}(\Sigma^{1/2}z, z^{\mathrm{T}}\Sigma^{1/2}A\Sigma^{1/2}z) + 2\operatorname{Cov}(\Sigma^{1/2}z, \mu^{\mathrm{T}}A\Sigma^{1/2}z). \end{align*}
First, since $z$ is a standard Gaussian, we have that $z \stackrel{\mathrm{d}}{=} -z$, and so \begin{align*} \operatorname{Cov}(\Sigma^{1/2}z, z^{\mathrm{T}}\Sigma^{1/2}A\Sigma^{1/2}z) &= \operatorname{Cov}(\Sigma^{1/2}(-z), (-z)^{\mathrm{T}}\Sigma^{1/2}A\Sigma^{1/2}(-z)) \\ &= -\operatorname{Cov}(\Sigma^{1/2}z, z^{\mathrm{T}}\Sigma^{1/2}A\Sigma^{1/2}z), \end{align*} from which it follows that $\operatorname{Cov}(\Sigma^{1/2}z, z^{\mathrm{T}}\Sigma^{1/2}A\Sigma^{1/2}z) = 0$.
Secondly, since $\mathrm{E}z = 0$, \begin{align*} \operatorname{Cov}(\Sigma^{1/2}z, \mu^{\mathrm{T}}A\Sigma^{1/2}z) = \mathrm{E}[\Sigma^{1/2}z z^{\mathrm{T}} \Sigma^{1/2} A \mu] = \Sigma A \mu, \end{align*} which yields the result.