A Multivariate Gaussian is part of the Exponential Family $$p(x,y|\eta)=h(x,y)\exp\left\lbrace \eta^TT(x,y)-A(\eta)\right\rbrace $$ Where the Sufficient Statistics are $$T(x,y)=\begin{bmatrix}x\\y\\xx^T\\xy^T\\yy^T \end{bmatrix}$$ I am working on a problem which is concerned with $Var_{p(x,y)}(T(x,y))$. I can easily find this for the first two sufficient statistics but am hitting a wall for the final 3. For example \begin{align*} Var_{p(x,y)}(xy^T)&=\mathbb{E}_{p(x,y)}\left[xy^T(xy^T)^T\right]-\mathbb{E}_{p(x,y)}\left[xy^T\right]\mathbb{E}_{p(x,y)}\left[xy^T\right]^T\\ &=\mathbb{E}_{p(x,y)}\left[xy^Tyx^T)^T\right]-\mathbb{E}_{p(x,y)}\left[xy^T\right]\mathbb{E}_{p(x,y)}\left[xy^T\right]^T\\ &=\mathbb{E}_{p(x,y)}\left[xy^Tyx^T\right]-\mathbb{E}_{p(x,y)}\left[xy^T\right]\mathbb{E}_{p(x,y)}\left[xy^T\right]^T\\ \end{align*} To simply focus on the first term \begin{align*} \mathbb{E}_{p(x,y)}\left[xy^Tyx^T\right]&=\mathbb{E}_{p(y)}\left[\mathbb{E}_{p(x|y)}\left[xx^T\right]y^Ty\right]\\ &=\mathbb{E}_{p(y)}\left[\left(\Sigma_{x|y}+\mu_{x|y}^2\right)y^Ty\right]\\ &=\mathbb{E}_{p(y)}\left[\left(\Sigma_{xx}-\Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}+\left(\mu_x+\Sigma_{xy}\Sigma_{yy}^{-1}(y-\mu_y)\right)^2\right)y^Ty\right] \end{align*} The issue I am having here is I dont know how to compute these. I am unsure of what $\mathbb{E}_{p(y)}[y^Ty]$ is (the transpose is on the "wrong" variable so this is not the second moment). I also realize that expanging out the squared term will result in some term like $yy^Ty^Ty$ (I may have the transposes on the wrong terms) but I am unsure of how to compute this term. I know that these terms showing up will be skewness and kurtosis but am unsure of what those evaluate to or how to work with them.
Any advice on how to proceed with this calculation or a reference to the solution of it?
