Let $W_t$ be a standard Brownian motion, i.e., $W_t \sim N(0,t)$. Define the random variable
$$X=\int_0^1\cos(W_t)dt$$
A similiar process, $Y_t=\cos(\omega t+\sigma W_t+\theta)$, with the uniform random variable $\theta\sim U(-\pi,\pi)$ and $\omega\in\Bbb R$, was analyzed by Cai and Wu [1], where the probability density was defined as
$$p_Y(y)=\frac{1}{\pi\sqrt{1-y^2}}$$
Seeing that the process $Y_t$ has a density, how can I obtain the density of the simplified process $\cos(W_t)$ ? I expanded the integral as a Riemann sum,
$$X\approx\frac{1}{n}\sum_{i=1}^n{\cos\left( W_\frac{i}{n}\right)}=\frac{1}{n}\left(n\cos\left(W_\frac{1}{n}\right)+ (n-1)\left(\cos\left(W_\frac{2}{n}\right)-\cos\left(W_\frac{1}{n}\right)\right)+ (n-2)\left(\cos\left(W_\frac{3}{n}\right)-\cos\left(W_\frac{2}{n}\right)\right)+ \cdots + 2\left(\cos\left(W_\frac{n-1}{n}\right)-\cos\left(W_\frac{n-2}{n}\right)\right)+ \cos\left(W_1\right)-\cos\left(W_\frac{n-1}{n}\right) \right)$$
and after rearranging the terms I got the sum $$\frac{1}{n}\sum_{i=1}^n{\cos\left( W_\frac{i}{n}\right)}= \frac{1}{n}\sum_{i=1}^n{-2(n-i+1) \sin\left(\frac{W_\frac{i}{n}+W_\frac{i-1}{n}}{2}\right) \sin\left(\frac{W_\frac{i}{n}-W_\frac{i-1}{n}}{2}\right)}$$
How can I progress from here?
[1] CAI, G. Q.; WU, C. Modeling of bounded stochastic processes. Probabilistic Engineering Mechanics, v. 19, n. 3, p. 197–203, 2004.
I edited my response many times to correct some mistakes. Sorry.
We can compute the density of $U_t=\cos(W_t)$ by appropriately separating the integral and a change of variable. For any measurable $f$
\begin{align} \mathbb{E}(f(U_t)) &= \mathbb{E}(f(\cos(W_t))) = \int_\mathbb{R} f(\cos(w)) \mathbb{P}_{W_t}(w) \mathrm{d} w = \sum_{k \in \mathbb{Z}} \int_{k\pi}^{(k+1)\pi} f(\cos(w)) \mathbb{P}_{W_t}(w) \mathrm{d} w \\ & = \sum_{k \in \mathbb{Z}} \int_{-1}^{1} f(u) \left( \mathbb{P}_{W_t}(\arccos(u)+2k\pi)+\mathbb{P}_{W_t}(\arccos(-u)+(2k+1)\pi) \right)\frac{ \mathrm{d} u}{\sqrt{1-u^2}} \\ &= \int_{-1}^1 f(u) \frac{\sum_{k\in \mathbb{Z}} \mathbb{P}_{W_t}(\arccos(u)+2k\pi)+\mathbb{P}_{W_t}(\arccos(-u)+(2k+1)\pi)}{\sqrt{1-u^2}} \mathrm{d}u \end{align}
From the first to the second line we need to separate odd and even integers but we can group the sums together after that.
So the density of $U_t=\cos(W_t)$ is \begin{equation} \mathbb{P}_{U_t}(u) = \frac{\sum_{k\in \mathbb{Z}} \mathbb{P}_{W_t}(\arccos(u)+2k\pi)+\mathbb{P}_{W_t}(\arccos(-u)+(2k+1)\pi)}{\sqrt{1-u^2}} \end{equation}
Numerical simulations for $\mathbb{P}_{U_t}$:
The numerator is an infinite sum as it appears in the wrapped Gaussian. But there is this $\arccos$ and the denominator. I don't know if it can help to evaluate the density of the integral.
The problem is that $U_t$ is neither Gaussian nor independent at different time step which makes difficult to evaluate the density of the Riemann sum.
Numerical simulations give the following histogram for X:
I don't have any idea of what function it could be.
Interestingly, using a sinus ie $U_t = \sin(W_t)$ seems tractable. It's probably because $U_t$ is centered. To get the density, just replace $\arccos$ by $\arcsin$ above. See the figures. For $\mathbb{P}_{U_t}$
And for the density of $X$, I'm able to fit a density $ \propto \exp\left( -\frac{s}{1-x^2} \right) $:
I hope it could be helpful.