Probability density function of the integral of a continuous stochastic process I am interested in whether there is a general method to calculate the pdf of the integral of a stochastic process that is continuous in time.
My specific example: I am studying a stochastic given process given by
$X(t)=\int\limits_{0}^{t}\cos(B(s))\,\text{d}s$,
where $B(t)$ is the Wiener process, which is normally distributed over an interval of length $\tau$ with zero mean and variance $\tau$:
$B(t+\tau)-B(t)\sim\mathcal{N}(0, \tau)$.
I am able to calculate the first and second moments of $X(t)$, see: Expectation value of a product of an Ito integral and a function of a Brownian motion
A couple of thoughts on the matter:
1) Integrals of Gaussian continuous stochastic processes, such as the Wiener process can be considered as the limit of a sum of Gaussians and are hence themselves Gaussian.  Since $\cos(B(s))$ is not Gaussian, this doesn't seem to help here.
2) If we can derive an expression for the characteristic function of the process $X(t)$, then we can theoretically invert this to obtain the pdf.  The Feynman-Kac formula enables us to describe the characteristic function in terms of a PDE.  If this PDE has a unique analytic solution then we can make use of this.  In my specific example, this is not the case - the PDE obtained has no analytic solution.  I can provide more detail on this point if required.
Many thanks for your thoughts.
 A: You may use a computational method to approximate as accurate as wanted the probability density function of $I(t)=\int_0^t \cos(B(s))\,\mathrm{d}s$. I will do so for $0\leq t\leq 1$.
Consider the Karhunen-Loève expansion of $B(t)$ on $[0,1]$:
$$ B(t)=\sum_{j=1}^\infty \frac{\sqrt{2}}{\left(j-\frac12\right)\pi}\sin\left(\left(j-\frac12\right)\pi t\right)\xi_j, $$
where $\xi_1,\xi_2,\ldots$ are independent and $\text{Normal}(0,1)$ distributed random variables. The convergence of the series holds in $L^2([0,1]\times\Omega)$. Truncate
$$ B_J(t)=\sum_{j=1}^J \frac{\sqrt{2}}{\left(j-\frac12\right)\pi}\sin\left(\left(j-\frac12\right)\pi t\right)\xi_j, $$
and define $I_J(t)=\int_0^t \cos(B_J(s))\,\mathrm{d}s$. 
It is easy to see that:


*

*$I_J(t)\rightarrow I(t)$ in $L^2(\Omega)$, for each $0\leq t\leq 1$. Indeed, by Cauchy-Schwarz inequality, $\mathbb{E}[(I(t)-I_J(t))^2]\leq t\|B(t)-B_J(t)\|_{L^2([0,t]\times\Omega)}^2\rightarrow 0$, as $J\rightarrow\infty$.

*$I_J(t)\rightarrow I(t)$ almost surely, for each $0\leq t\leq 1$. Indeed, for each fixed $\omega\in\Omega$, we know that $B_J(t)(\omega)\rightarrow B(t)(\omega)$ in $L^2([0,1])$, because deterministic Fourier series converge in $L^2$. Since $\cos$ is Lipschitz, $\cos(B_J(t)(\omega))\rightarrow \cos(B(t)(\omega))$ in $L^2([0,1])$. Then $I_J(t)(\omega)\rightarrow I(t)(\omega)$ for each $t\in[0,1]$ follows.
Although these two facts are not enough to have that the density functions of $\{I_J(t):J\geq1\}$ tend to the density function of $I(t)$, we have that the density function of $I_J(t)$ may be a very good candidate (recall that this is a computational method, not a proof). The good thing about $I_J(t)$ is that is consists of a finite number of $\xi_j$, so it is possible to obtain exact realizations of $I_J(t)$. And if we generate a sufficiently large number $M$ of realizations of $I_J(t)$, then a kernel density estimation allows obtaining an approximate density function for $I_J(t)$.
I have written a function in Mathematica to approximate the distribution of $I(T)$, for $0\leq T\leq 1$, using a truncation order $J$ and a number of simulations $M$:
distributionIT[T_, J_, simulations_] := 
  Module[{realizationsIT, simulation, xi, BJ, distribution},
   realizationsIT = ConstantArray[0, simulations];
   For[simulation = 1, simulation <= simulations, simulation++,
    xi = RandomVariate[NormalDistribution[0, 1], J];
    BJ[t_] := 
     Sqrt[2]*Sum[
       Sin[(j - 0.5)*Pi*t]/((j - 0.5)*Pi)*xi[[j]], {j, 1, J}];
    realizationsIT[[simulation]] = NIntegrate[Cos[BJ[t]], {t, 0, T}];
    ];
   distribution = SmoothKernelDistribution[realizationsIT];
   Return[distribution];
   ];

Let us do a numerical example. Choose $T=1$. Write
distribution1 = distributionIT[1, 40, 50000];
plot1 = Plot[PDF[distribution1, x], {x, -2, 2}, Frame -> True, PlotRange -> All];
distribution2 = distributionIT[1, 50, 100000];
plot2 = Plot[PDF[distribution2, x], {x, -2, 2}, Frame -> True, PlotRange -> All];
Legended[Show[plot1, plot2], LineLegend[{Green, Blue}, {"J=40, M=50000", "J=50, M=100000"}]]


We plot the estimated density function for $J=40$, $M=50000$ and $J=50$, $M=100000$. We observe no differences, so our method provides a good approximation of the probability density function of $I(1)$.
Similar computations for $T=0.34$ give the following plot:

If you plot the approximate density function for smaller $T$, you will see that at the end one gets a Dirac delta at $0$, which agrees with the fact that $I(0)=0$ almost surely.
Remark: Computational methods are of constant use in research to approximate probabilistic features of response processes to random differential equations. See for example [M. A. El-Tawil, The approximate solutions of some stochastic differential equations using transformations, Applied
Mathematics and Computation 164 (1) 167–178, 2005], [D. Xiu, Numerical Methods for Stochastic Computations. A Spectral Method Approach, Princeton University Press, 2010], [L. Villafuerte, B. M. Chen-Charpentier, A random differential transform method: Theory and applications, Applied Mathematics Letters, 25 (10) 1490-1494, 2012].
A: Well the first moment is relatively straightforward. If the integral exists than 
$E[\int_a^b X(t) dh(t)] = \int_a^b E[X(t)]dh(t)$
See for example Grigoriu 3.69.
Also by for example Grigoriu 3.58
$E[ \int_a^b X(t) dt \int_a^b Z(t) dt] = \int_{[a,b]^2}E[X(u)Z(v)]du dv$
