Schwartz regularity for the density of a stochastic process Let $B$ be a standard Brownian motion in $\mathbb R$. Define the variables
$$\begin{align*} X &= B_1, & Y &= \int_0^1B_s\mathrm ds, & Z&= \int_0^1B_s^2\mathrm ds. \end{align*}$$
It is known, see below, that $(X,Y,Z)$ admits a smooth density.

Is is true that the density of $(X,Y,Z)$ is a Schwartz function?

I think there should be a simple enough argument, and I'm making it more complicated than it actually is.

Let me share two observations about this triple. First of all, all of its components belong to low order Wiener chaos: specifically, $X$ and $Y$ belong to the first Wiener chaos, and $Z$ to the second. This makes it easy enough to compute the Malliavin matrix ($\langle DX,DX\rangle$ and its peers):
$$ \Gamma
 = \begin{pmatrix}
     1       & 1/2 & 2\alpha \\
     1/2     & 1/3 & 2\alpha-\beta \\
     2\alpha & 2\alpha-\beta & 4\gamma
   \end{pmatrix} $$
with
$$ \begin{align*} \alpha&=\int_0^1tB_t\mathrm dt, & \beta&=\int_0^1t^2B_t\mathrm dt, & \gamma&=\int_0^1\Big(\int_t^1B_s\mathrm ds\Big)^2\mathrm dt. \end{align*} $$
Using this expression, we get $\det\Gamma=(\gamma-\alpha^2)/3 - (\alpha-\beta)^2$, and it would suffice to prove that $\mathbb P[\det\Gamma\leq\varepsilon]$ decreases faster than any power of $\varepsilon$, or that $1/\det\Gamma$ has moments of all order (the former implies the latter). This would follow from Proposition 2.1.5 in Nualart's The Malliavin Calculus and Related Topics (p. 103), as suggested by Fabrice Baudoin in the comments.
As for the second observation, the variable in question is the evaluation at time 1 of a hypoelliptic process started at $(0,0,0)$: if one replaces every 1 by a $t$ in the above definitions, we get a process $(X_t,Y_t,Z_t)$ solution to the hypoelliptic equation
$$\begin{align*} \mathrm dX_t &= \mathrm dB_t, & \mathrm dY_t &= X_t\mathrm dt,& \mathrm dZ_t&= X_t^2\mathrm dt \end{align*}$$
such that the value at time 1 is the variable described above. This shows that the density of $(X_1,Y_1,Z_1)$ is smooth. One can find more about this approach (with slightly different conventions) in section IV-1.5 (page 118) of my PhD thesis, but here are a few key insights. We can find a (not so nice) candidate for the characteristic function of $(X_t,Y_t,Z_t)$ (hence for $(X,Y,Z)$ as well) in the form
$$\text{constant}\times\exp\big(\text{quadratic form in }(\xi_x,\xi_y)\big)$$
for $t$ and $\xi_z$ fixed. Specifically,
$$ \mathbb E\big[\exp\big(-i(\xi_xX_t+\xi_yY_t+\xi_zZ_t)\big)\big]
 = \frac1{\sqrt{\mathrm{cah}}}\cdot
   \exp\left(
     - \frac t2\mathrm{tah}\cdot\xi_x^2
     + \frac i{2\xi_z}\Big(1-\frac1{\mathrm{cah}}\Big)\cdot\xi_x\xi_y
     + \frac{it}{4\xi_z}(1-\mathrm{tah})\cdot\xi_y^2\right), $$
for $\mathrm{cah} = \cosh(\zeta)$, $\mathrm{sah} = \sinh(\zeta)/\zeta$, $\zeta = \sqrt{2i\xi_zt^2}$ (these are quantities of the form $f(\sqrt{(\cdots)})$ for $f$ holomorphic and even, so they are actually holomorphic),  $\mathrm{tah} = \mathrm{cah}/\mathrm{sah}$. One can show that the zeros of $\alpha\mapsto\cosh(\sqrt\alpha)$ are real negative, so everything can be defined holomorphically in the above formula.
This is a nice candidate because it is one solution of the Feynman-Kac formula seen in Fourier space, in a rather strong sense. If it were the actual characteristic function, then the question would be equivalent to show that the candidate is itself a Schwartz function (in the dual space variables). However, this possibly equivalent result does not seem obvious to me either, since one would have to show that the real part of the quadratic form is bounded below, or something of the sort.
 A: Using the representation in terms of i. i. d. Gaussians $\xi_1,\xi_2,\dots,$
$$
B_t=\sqrt{2}\sum_{n=1}^\infty (-1)^{n+1}\xi_n\frac{\sin \pi \left(n-\frac12\right)t}{\pi \left(n-\frac12\right)}, 
$$
we get
$$
X=\sqrt{2}\sum_{n=1}^\infty\frac{\xi_n}{\pi \left(n-\frac12\right)},\quad Y=\sqrt{2}\sum_{n=1}^\infty \frac{(-1)^{n+1}\xi_{n}}{\pi^2 \left(n-\frac12\right)^2},\quad Z=\sum_{n=1}^\infty\frac{\xi^2_n}{\pi^2 \left(n-\frac12\right)^2}.
$$
The Fourier transform of a square of a Gaussian can be computed explicitly:
$$
\frac{1}{\sqrt{2\pi}}\int_\mathbb{R}  e^{-\frac{x^2}{2}+itx^2}\,dx=
\frac{1}{\sqrt{2\pi}} \frac{1}{\sqrt{1-2it}}\int_{z\in\sqrt{1-it}\mathbb{R}}e^{-\frac{z^2}{2}}dz=\frac{1}{\sqrt{1-2it}},
$$
by considering the contour integration over the "pizza slice" domain. Hence, we have $$
\varphi_Z(t)=\prod^\infty_{n=1}\left(1-\frac{2it}{\pi^2 \left(n-\frac12\right)^2}\right)^{-\frac12}=\left(\cos(\sqrt{2it})\right)^{-\frac12},
$$ which is clearly of Schwarz class on $\mathbb{R}$, and hence so is the density of $Z$.
The joint characteristic function of $(X,Y,Z)$ can be computed similarly. Define independent random variables
$$
(X_n,Y_n,Z_n)=\left(\frac{\sqrt{2}\xi_n}{\pi \left(n-\frac12\right)},\frac{\sqrt{2}\xi_{n}}{\pi^2 \left(n-\frac12\right)^2},\frac{2\xi^2_n}{\pi^2 \left(n-\frac12\right)^2}\right)
$$
so that they are independent and
$$
(X,Y,Z)=\sum_{n=1}^\infty(X_n,Y_n,Z_n).
$$
We can compute the characteristic function of $(X_n,Y_n,Z_n)$ using the identity
$$
\frac{1}{\sqrt{2\pi}}\int_\mathbb{R}e^{-\frac{x^2}{2}+ax+ibx^2}\,dx=\frac{1}{\sqrt{2\pi}}\int_\mathbb{R}e^{-\frac{\left(\sqrt{1-2ib}x-\frac{a}{\sqrt{1-2ib}}\right)^2}{2}+\frac{a^2}{2(1-ib)}}\,dx=e^{\frac{a^2}{2(1-2ib)}}\frac{1}{\sqrt{1-2ib}}.
$$
This yields, denoting $\varphi_n=\varphi_{(X_n,Y_n,Z_n)},$
$$
\varphi_n(t_1,t_2,t_3)=
\exp\left(-\frac{\left(t_1+\frac{(-1)^{n+1}t_2}{\pi\left(n-\frac12\right)}\right)^2}{\pi^2 \left(n-\frac12\right)^2\left(1-\frac{2it_3}{\pi^2 \left(n-\frac12\right)^2}\right)}\right)\left(1-\frac{2it_3}{\pi^2 \left(n-\frac12\right)^2}\right)^{-\frac12}
$$
After multiplying this, the idea is to show that the product of the exponentials gives a function that decays exponentially with all derivatives outside any cone containing the $t_3$ axis (and each of its derivatives is bounded by a polynomial in the cone), while the rest gives $\varphi_Z(t_3)$ which decays fast enough in the cone, and thus conclude that we get a Schwarz class function.
To elaborate, denote
$$
\Phi(t_1,t_2,t_3)=\sum_n\frac{\left(t_1+\frac{(-1)^{n+1}t_2}{\pi\left(n-\frac12\right)}\right)^2}{\pi^2 \left(n-\frac12\right)^2\left(1-\frac{4it_3}{\pi^2 \left(n-\frac12\right)^2}\right)};
$$
it suffices to show that, on the one hand,  $\Phi$ and all its derivatives grow at most polynomially as $(t_1,t_2,t_3)\to\infty$ (which is straightforward, just differentiate and estimate the absolute value of the (power of) the last term in the denominator as $\geq 1$), and that outside the above-mentioned cone,
$$
\Re\Phi(t_1,t_2,t_3)>c(|t_1|+|t_2|)$$ for some constant $c$, which seems to follow by looking at the first two terms in the sum and noticing that the real part of the others is non-negative. Indeed, with these observation at hand, we have
$$
\varphi_{(X,Y,Z)}=\exp(-\Phi(t_1,t_2,t_3))\cos(\sqrt{2it_3})^{-\frac12}
$$
The derivatives of this guy will have the form
$$
Poly(t_1,t_2,t_3)e^{-\Phi(t_1,t_2,t_3)},
$$
where $Poly$ is some polynomial in derivatives of $\Phi(t_1,t_2,t_3)$ and in $\cos(\sqrt{2it_3})^{-\frac12}$ and its derivatives, in which every monomial is linear in  $\cos(\sqrt{2it_3})^{-\frac12}$ or one of its derivatives. Because of this last property, the whole expression will decay faster than polynomial in the cone, and the exponential will kill it off outside the cone.
