I am eager to compute the following variance, where $W_t$ is a Wiener process:
$$\mathbb{V}\left[I_t\right], \text{ where } I_t = \int\limits_0^tW_u^2 u^2\partial u.$$
I have already computed $\mathbb{E}[I_t] = \Large\frac{t^4}{4}$, and hence currently I have:
$$ \mathbb{V}[I_t] = \mathbb{E}[I^2_t] - \mathbb{E}^2[I_t] = \mathbb{E}[I^2_t] - \frac{t^{16}}{16} = \mathbb{E}\left[\left(\int\limits_0^tW_u^2 u^2\partial u\right)^2\right] - \frac{t^{8}}{16}. $$
And here, I do not understand, how to compute the first summand. I thought of applying the Itô isometry here, but as long as $W_u^2u^2$ is a stochastic process (at least it seems for me it is, what may be verified), the differential is taken not over a Wiener process, but rather over the time-like variable $u$, so that I have no further ideas on how to proceed.
Any help would be appreciated, thank you in advance!
I have limited time just now, and can fill in more details later if you want. Idea comes from here $$\Bbb E[( \int_0^t W_u^2 u^2 d u)^2 ] = \Bbb E[ \int_0^t W_u^2 u^2 d u \cdot \int_0^t W_s^2 s^2 d s ] = \Bbb E [\int_0^t\int_0^t W_s^2 W_u^2 s^2u^2 d u ds ]\\ =\int_0^t\int_0^t \Bbb E[ W_s^2 W_u^2 ] s^2u^2 d u ds $$ and
$$E[ W_s^2 W_u^2 ] = E[ W_{\min (s,u)}^2 (W_{\max(s,u)}- W_{\min (s,u)} + W_{\min (s,u)})^2 ] \\ = \Bbb E[W_{\min (s,u)}^4] + \Bbb E[W_{\min (s,u)}^2]\Bbb E[(W_{\max(s,u)} -W_{\min (s,u)})^2], $$ which can be further and explicitly computed.