Second order taylor expansion for multivariate functions of random variables Suppose one has $f(x,y)$ with $x,y$ centered random variables i.i.d. with finite variance. We can develop $f(x,y)$ around $f(E[x], E[y]) = f(0,0)$ so that:
$f(x,y) = f(0,0)+f_xx+f_yy+f_{xx}(x)^2/2+f_{yy}(y)^2/2+f_{xy}xy+remainder$
Then is $f_x$ and $f_y$ are null, is it correct to estimate the variance of $f$ as:
$E[(f(x,y)-f(0,0))^2]$
by plugging in the equation the series expansion of $f$?
(I assume that the expectance of $f$ can be approximated as $f(E[x],E[y])$).
Specifically, my question is what happens when estimating using this method the covariance of $f(x,y)$ and $f(z,w)$ where $x,y,z,w$ are zero mean i.i.d.; since I obtain a value which is different from zero though the two functions are functions of independent random variables and therefore the covariance shall be zero. Any suggestion? Does it mean that the procedure is wrong or am I miscalculating something?
Thanks a lot for your time!
[Corrected typo in the expansion]
 A: Maybe this will clear up some things:
In general, you don't have $\mathbb{E}[f(X,Y)] = f(\mathbb{E}[X],\mathbb{E}[Y])$. For an explicit counterexample let $X,Y $  be $ N(0,1)$ and independent and let $f(x,y) := xy$. In this case we have
$$ 1 = \mathbb{E}[f(X,Y)] \neq f(\mathbb{E}[X],\mathbb{E}[Y]) = 0$$
Even if you assume this is true the next step wont work this way. Let's consider the expansion. We assume that $\mathbb{E}[f(X,Y)] \approx f(0,0)$. With $f_x (0,0) = f_y (0,0) = 0$ we obtain
$$f(x,y) = f(0,0) + \frac{x^2}{2} f_{xx} (0,0) + \frac{y^2}{2} f_{yy} (0,0) + xy f_{xy} (0,0) + \mathcal{o}(x^2) + \mathcal{o}(y^2) + \mathcal{o}(xy)$$
Now you could try to approximate $\mathbb{E}[(f(x,y)-f(0,0))^2]$ with the taylor expansion from above, but a priori there is no reason why the remainder terms should integrate to $0$. For example, the remainder term $\mathcal{o}(x^2)$ could look like $x^4$ which means you could make an error of up to $\mathbb{E}[X^4]$ when calculating the covariance in this way.
