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 this representation
$$\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$

But how do you show $E[(\int_0^T W_s ds)^2] < \infty$?
Actually, I can use Ito's lemma, $(*)$ 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$$
But I need the reverse. How do we compute $E[(\int_0^T W_s ds)^2]$ besides using the representation (if it's possible)? I don't think we encountered something like $\int_0^T W_s ds$ in classes.
 A: 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]\\
&=E\left[\int_0^T\int_0^TW_t\,W_s\,\mathrm dt\mathrm ds\right]\\
&=\int_0^T\int_0^T E[W_tW_s]\,\mathrm dt\,\mathrm ds\\
&=\int_0^T\int_0^T \sigma^2\min(t,s)\,\mathrm dt\,\mathrm ds\\
&=\int_0^T\left[\int_0^T\sigma^2\min(t,s)\,\mathrm dt\right]\,\mathrm ds\\
&=\int_0^T\left[\int_0^s\sigma^2 \min(t,s)\,\mathrm dt
+ \int_s^T\sigma^2 \min(t,s)\,\mathrm dt\right]\,\mathrm ds\\
&=\int_0^T \int_0^s\sigma^2 \min(t,s)\,\mathrm dt\,\mathrm ds
+ \int_0^T\int_s^T\sigma^2 \min(t,s)\,\mathrm dt\,\mathrm ds\tag{1}\\
&= \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\tag{2}\\
&=\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}<\infty
\end{align}
I hope that the answer above is understandable and acceptable to most readers, but for the benefit of @youpilat13 and any others confused as to how $(2)$ follows from $(1)$, I add the following explanation.
There are two double integrals in $(1)$ and they are over disjoint right triangular regions, separated by the line $s=t$, in the plane. As to why the two integrals are being summed, we are trying to compute the integral over a square region and we are breaking up the integral over the square into the sum of integrals over two disjoint subsets (right triangles), computing each integral separately, and then adding the results together to determine the integral over the square region.
In the inner integral in the first double integral in $(1)$, $s$ is a constant with respect to $t$, the variable of integration, but since $t$ varies from $0$ to $s$ only, it must be that $\min(t,s)$ equals $t$. Similarly, in the inner integral in the second double integral in $(1)$, $s$ is a constant with respect to $t$, the variable of integration, but since $t$ varies from $s$ to $T$ only, it must be that $\min(t,s)$ equals $s$.
