# How to compute expectation of square of Riemann integral of a random variable?

How does one compute $E[(\int_0^T W_s ds)^2]$ where $(W_t)_{t \in [0,T]}$ is standard Brownian motion in $(\Omega, \mathscr F, \mathbb P)$?

Apparently proving

$$\int_0^T W_s ds = \int_0^T (T-s) dW_s \tag{*}$$

might need to assume that $E[(\int_0^T W_s ds)^2] < \infty$

I think I can just use Ito's lemma to show $(*)$ and then use Ito isometry to show that

$$E[(\int_0^T W_s ds)^2] = E[(\int_0^T (T-s) dW_s)^2] = (\int_0^T (T-s)^2 ds) < \infty$$

Is that okay? If not, why?

Regardless, I'm still curious as to how compute $E[(\int_0^T W_s ds)^2]$ besides using the representation if it's possible? I don't think we encountered $\int_0^T W_s ds$ in classes

Assuming that $$E[W_tW_s]=\sigma^2\min(t,s)$$, \begin{align} E\left[\left(\int_0^TW_s\,\mathrm ds\right)^2\right]&=E\left[\int_0^TW_t\,\mathrm dt\int_0^TW_s\,\mathrm ds\right]\\ &=\int_0^T\int_0^T\sigma^2\min(t,s)\,\mathrm dt\,\mathrm ds\\ &= \int_0^T\int_0^s\sigma^2 t\,\mathrm dt\,\mathrm ds + \int_0^T\int_s^T\sigma^2 s\,\mathrm dt\,\mathrm ds\\ &=\int_0^T\sigma^2\frac{s^2}{2}\,\mathrm ds + \int_0^T\sigma^2(T-s)s\,\mathrm ds\\ &= \sigma^2\left(\frac{T^3}{6}+\frac{T^3}{2}-\frac{T^3}{3}\right)\\ &= \sigma^2\frac{T^3}{3} \end{align}