# 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.

• Are you still interested in this question? Mar 18, 2013 at 14:26
• I wasn't when you asked (I didn't have a spare moment to think about "interesting things"), but I am again now, should anyone see this! Mar 11, 2014 at 11:47

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:

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

2. $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].

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$

• I have expressions for the first two moments using a similar approach... but unless it's possible to derive expressions for an arbitrary number of moments (hence the MGF), this doesn't help with the pdf unfortunately. If such a thing is possible, then perhaps you could invert the MGF - but it seems unlikely to work nicely? Mar 11, 2014 at 11:51
• Sorry guess I didn't read your question carefully. Thought you wanted the first two moments. Mar 11, 2014 at 16:04
• @Gabriel, I think you could generalize that idea to get arbitrary moments. Probably would be hard to do the integrals in closed form. Do you really need the pdf ? Mar 11, 2014 at 16:11
• It's more a point of academic interest! Thanks for your thoughts. Mar 11, 2014 at 17:53