I want to calculate the variance of

$$I = \int\_0^t W\_s^2 ds$$

I was thinking I could define the function $f(t,W\_t) = tW\_t^2$ and then apply Ito's lemma so I get

$$f(t,W\_t)-f(0,0) = \int\_0^t \frac{\partial f}{\partial t}(s,W\_s)ds + \int\_0^t \frac{\partial f}{\partial x}(s,W\_s)dW\_s+ \frac{1}{2}\int\_0^t \frac{\partial^2 f}{\partial x^2}(s,W\_s)ds \\= I + \int\_0^t 2sW\_sdW\_s + \frac{t^2}{2}$$

By rearranging I get

$$I = tW\_t^2 - \int\_0^t 2sW\_sdW\_s - \frac{t^2}{2}$$

We then get that (I'm not sure here but i think the expectation is zero of any integral w.r.t BM?)

$$\mathbf{E}[I]=\frac{t^2}{2}$$

And variance

$$\mathbf{V}[I] = \mathbf{V}[tW\_t^2 - \int\_0^t 2sW\_sdW\_s - \frac{t^2}{2}]
= t^2\mathbf{V}[W\_t^2]+\mathbf{E}[(\int\_0^t 2sW\_sdW\_s)^2] \\= 2t^4 + \mathbf{E}[\int\_0^t 4s^2W\_s^2ds]\quad\text{(Isometry property)}$$

Not sure if it is OK to change order of integration and expectation here, but if I do that, I get

$\mathbf{V}[I]= 2t^4 + \int\_0^t 4s^2\mathbf{E}[W\_s^2]ds = 2t^4 + \int\_0^t 4s^2\mathbf{E}[W\_s^2]ds = 2t^4 + \int\_0^t 4s^3ds=3t^4$

However, the answer says the variance should be $\frac{t^4}{3}$, so I guess I do something wrong?