The area under a Brownian path

Idea. The variance of an integral of a process is the double integral of its covariance function. For Brownian motion the covariance is $\min (s,t)$, which is easy to integrate over the unit square.

Step 1. $\mathbb{E}[W_t]=0$ for every $t$, so $\mathbb{E}[I]=0$ and $$\begin{array}{rl}\mathrm{Var}(I)& =\mathbb{E}\left[{\int }_0^1{\int }_0^1{W}_s{W}_{t}dsdt\right]\\ & ={\int }_0^1{\int }_0^1\min (s,t)dsdt.\end{array}$$

Step 2. By symmetry, integrate over $s<t$ and double: $$2{\int }_0^1{\int }_0^{t}sdsdt=2{\int }_0^1\frac{t^2}{2}dt=\frac{1}{3}.$$

So $\mathrm{Var}(I)=\frac{1}{3}$. Since $I$ is a limit of sums of jointly normal variables, $I\sim N\left(0,\frac{1}{3}\right)$.

Another route. Integration by parts gives $${\int }_0^1{W}_{t}dt={\int }_0^1(1-t)d{W}_t,$$ and by the Itô isometry its variance is ${\int }_0^1(1-t{)}^{2}dt=\frac{1}{3}$.

Follow-ups. Over $[0,T]$ the variance is $\frac{T^3}{3}$. The covariance of $I$ with $W_1$ is ${\int }_0^{1}tdt=\frac{1}{2}$, so their correlation is $\frac{1/2}{\sqrt{1/3}}=\frac{\sqrt{3}}{2}$.