Integral of Wiener Squared process

That should not be true. The problem is that you're not considering the square of the stochastic integral, i.e. the random variable $J= (\int_0^1 {W_s} ds)^2$, which would indeed be the square of a normally distributed random variable, but the r.v. $S= (\int_0^1 {W_s}^2 ds)$.

A note in the direction of finding a solution would be a discretization attempt, that is defining an equally spaced partition $\Pi_{[0,t]}$ of $[0,t]$, $\Pi_{[0,1]}=\{0,t_1,...,t_i,...,t_{n^2}=1\}$ and considering the limit as the mesh of the partition goes to $0$ (here $\Delta t=1/n^2$):

$$ S=\lim_{|\Pi| \to 0} \Delta t\sum_{i \le n^2} W_{t_i}^2 $$

Here the random variables $W_{t_i}$ are independent, and it looks similar to result III. in this article by P. Erdös and M. Kac. Also check here.

Using Ito's formula on $f(x)=x^4$, we get for S the identity

$$ S=\int_0^t W^2_t dt =\frac{2}{3}\int_0^t W_t^3dWt - \frac{1}{6}W_t^4 $$ but I'm not sure this tells much about the r.v.

It would be also of interest to compute the moments of this rv. Using the Ito Isometry:

$$\mathbb{E}[S]=\mathbb{E}\left[ \int_0^t {W_s}^2 ds\right] = \mathbb{E}\left[\left( \int_0^t W_s dW_s\right)^2\right]=\\ \mathbb{E}\left[\left( \frac{1}{2} W_t^2-\frac{1}{2}t\right)^2\right]= \mathbb{E}\left[ \frac{1}{4} W_t^4+\frac{1}{4}t^2-\frac{1}{2}W_t^2t\right]= \frac{3}{4}t^2+\frac{1}{4}t^2-\frac{1}{2}t^2=\frac{t^2}{2}$$

where in the last step I used $\mathbb{E}[Z^4]=3\sigma^4$ for a Std Normal r.v. $Z$.

---

$$ \mathrm{Var}[S]=\mathrm{Var}\left[ \int_0^t {W_s}^2 ds\right] = \mathbb{E}\left[\left( \int_0^t {W_s}^2 ds\right)^2\right] - \left(\mathbb{E}\left[ \int_0^t {W_s}^2 ds\right]\right)^2=\\ \mathbb{E}\left[\left( \int_0^t {W_s}^2 ds\right)^2\right] - \frac{t^4}{4}= $$

To evaluate this last integral we can note that through Ito's formula

$$ 6\int_0^t W_t^2dt=W_t^4-4\int_0^t W_t^3dW_t $$

and $$\mathrm{Var}\left[ \int_0^t {W_s}^2 ds\right] =\frac{1}{36}\mathbb{E}\left[\left(W_t^4-4\int_0^t W_t^3dW_t\right)^2\right]-\frac{t^4}{4}$$ but here the computation gets more complicated.

The distribution of $\int_0^t W^2_s\,ds$ was found by Cameron and Martin in the 1940s. Mark Kac (1949, Transactions of the AMS) applied his method to find the Laplace transform of this integral: $$ \Bbb E\left[\exp(-u\int_0^t W^2_s\,ds)\right]={1\over\sqrt{\cosh(t\sqrt{2u})}}. $$