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.

-
An alternative approach is to discretize $(W_s)_{0\leqslant s\leqslant t}$, which yields your characterization of $X_t$ as a sum of Gammas. And you are probably well aware that $X_t$ and $t^2X_1$ have the same generalized chi-squared distribution‌​. –  Did Jan 14 '12 at 16:02
@DidierPiau Yes, I considered equal splitting with $s_i = t \frac{i}{n}$. $W_{s_i} = \sqrt{\frac{t}{n}} \sum_{k=1}^n Z_k$, for independ standard normals $Z_k$. This would give $X_1 = \frac{1}{n^2} \sum_{k=1}^n \sum_{m_1=1}^k \sum_{m_2=1}^k Z_{m_1} Z_{m_2}$. Because of cross terms it was difficult to separate that into sums of independent generalized $\chi^2$-distributions. It is natural to expect $\chi^2$ to appear, but conspiracy of scales is that some requires an explanation, while what I now have is a guesswork. –  Sasha Jan 14 '12 at 17:45

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}$.

-
Thank you. This is the link to the issue containing Mark Kac's article on AMS.org available free of charge. –  Sasha Jan 15 '12 at 0:54

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}.$$

-
Beautiful! Thank you very much. –  Sasha Jan 14 '12 at 17:50