Determining distribution of $X_t = \int_0^t W_s^2 \mathrm{d} s$ Premise
Let $W_t$ be the standard Wiener process, and let $X_t = \int_0^t W_s^2 \mathrm{d} s$. I am interested in determining the distribution of $X_t$. 
What I did
My line of attack has been to determine moments $\mathbb{E}(X_t^r)$. To this end I set up a system of It$\bar{\text{o}}$ differential equations:
$$
   \mathrm{d} X_t = Y_t^2 \mathrm{d} t, \qquad \mathrm{d} Y_t = \mathrm{d} W_t, \qquad X_0 = Y_0 = 0
$$
Then
$$
    \mathbb{E}(X_t^r) = \int_0^t \mathbb{E}( \mathcal{L}_0(X_s^r) ) \mathrm{d} s = 
        \int_0^t \mathrm{d}s_1 \int_0^{s_1} \mathrm{d}s_2 \cdots \int_0^{s_{n-1}} \mathrm{d} s_n \mathbb{E}( \mathcal{L}_0^{ \circ n}(X_{s_n}^r) ) 
$$
where $n$ is chosen so that $\mathcal{L}_0^{ \circ n}(X_{s_n}^r)$ is constant equal to $a_{n,r}$, then $\mathbb{E}(X_t^r) = \frac{t^r}{n!} a_{n,r}$, where
$$
   \mathcal{L}_0(f(X_t, Y_t)) = \left(\frac{1}{2} \frac{\partial^2}{\partial Y_t^2} + Y_t^2 \frac{\partial}{\partial X_t} \right)(f(X_t, Y_t))
$$
It is not hard to see that 
$$
   \mathbb{E}(X_t^r) = \frac{t^r}{(2r)!} \left( \frac{1}{2} \frac{\mathrm{d}^2}{\mathrm{d} y^2} + y^2 \frac{\mathrm{d}}{\mathrm{d} x} \right)^{2r} (x^r)
$$
Poking around OEIS, I stumbled across A126156, corresponding to $\frac{(2r)!}{r!} \mathbb{E}(X_1^r)$, from where I learnt the moment generating function for $X_t$:
$$
  \phi(u)= \mathbb{E}(\exp(u X_t)) = \frac{1}{\sqrt{ \cos(\sqrt{2 u} t)}} = 1 + \frac{u t^2}{1!} \cdot \frac{1}{2} + \frac{u^2 t^4}{2!} \cdot \frac{7}{12} + \frac{u^3 t^6}{3!} \cdot \frac{139}{120} + \mathcal{o}(u^3)
$$
for $-\infty < u < \frac{\pi^2}{8 t^2}$, i.e. until the first zero of the cosine, which indicates that the main asymptotic term of the probability density is $f_{X_t}(x) \sim \exp\left(- \frac{\pi^2}{8 t^2} x \right)$. Indeed, using Weierstrass representation:
$$
  \phi(u) = \prod_{n=0}^\infty \frac{1}{\sqrt{1-\frac{8 u t^2}{\pi^2} \frac{1}{(2n+1)^2}}}
$$
which suggests that $X_t = \sum_{n=0}^\infty Z_n$, where $Z_n$ are independent and $Z_n \sim \Gamma\left(\frac{1}{2}, 2 \left( \frac{2 t}{(2n+1) \pi} \right)^2 \right)$, which in turn means that $\kappa_r(X_t) = -(-8 t^2)^r (r-1)! (4^r-1) \frac{B_{2r}}{4 (2r)!} $
This unexpected success gives (slim) hope that the distribution function of $X_t$, or its probability density can be found in closed form.
Question
I am looking for references where $X_t$, or $X_1$, has been studied, or possible suggestions on how one would go about proving my guesswork.
 A: I haven't checked the details, but I think this approach works. 
First off, as Didier Piau notes in a comment, it suffices to find the distribution of $X_1$.
Now, use the Karhunen–Loève expansion of the Wiener process:
$$W(s)=\sqrt{2}\,\sum_{k=1}^\infty Z_k {\sin((k-1/2)\pi s)\over(k-1/2)\pi},$$
where $Z_k$ are i.i.d. standard normal random variables. 
When you square this and integrate, all the cross-product
terms disappear and you are left with the representation:
$$X_1=\int^1_0 W^2(s)\,ds=\sum_{k=1}^\infty {Z_k^2 \over (k-1/2)^2\pi^2}.$$
A: The distribution of $X_t$ was first computed by Cameron & Martin in 1944.  A reference to this paper, along with a nifty computation of the Laplace transform of $X_t$, can be found in a paper of Mark Kac: "On distributions of certain Wiener functionals" [Trans. Amer. Math. Soc. vol. 65 (1949) pp. 1–13].  Kac's method is an instance of what we would nowadays call the "Feynman-Kac method".  Kac shows that $E[\exp(-uX_t)] = [{\rm sech}((2u)^{1/2}t)]^{1/2}$.
