Expected value of time integral of a gaussian process

Question

While working on a problem I've stumbled upon some expected values of time integrals of Gaussian stochastic processes. Some of them were addressed this question, but I have found also this one

$$\left\langle\int _0^tdt_1\int_0^{t_1}dt_2\left(B\left(t_1\right)-B\left(t_2\right)\right)\int _0^tdt_3\int_0^{t_3}dt_4\left(B\left(t_3\right)-B\left(t_4\right)\right)\right\rangle$$

Here $B(t)$ is a stationary Gaussian process with known autocorrelation function $K(t_1-t_2)=\langle B(t_1)B(t_2) \rangle$ and $\langle \cdot \rangle$ is the expected value over all possible realizations of the stochastic process.

I tried writing the integral as

$$\left\langle\int _0^tdt_1\int_0^{t_1}dt_2\int _0^tdt_3\int_0^{t_3}dt_4\left[B\left(t_1\right)B(t_3)-B(t_1)B(t_4)-B(t_2)B\left(t_3\right)+B(t_2)B\left(t_4\right)\right]\right\rangle$$

and then exchanging the expectation value with the integral. This way I can integrate the known form of $K(t_i-t_j)$. However there must be something wrong with my reasoning because I expect this quantity to be positive while, for example by choosing a Ornstein-Uhlenbeck process, with $K(t_i-t_j) = \frac \gamma 2 e^{-\gamma |t_i-t_j|}$, I obtain non-positive function, namely

$$ \frac{4}{\gamma ^3} -\frac{5 t^2}{4 \gamma } + \frac{5 t^3}{12} - e^{-\gamma t} \left(\frac{4}{\gamma ^3}+ \frac{4 t}{\gamma ^2} +\frac{t^2}{\gamma }\right) $$

Answer 1

Indeed, the result is $$\int _0^t\mathrm dt_1\int_0^{t_1}\mathrm dt_2\int _0^t\mathrm dt_3\int_0^{t_3}\mathrm dt_4\left[K(t_1-t_3)-K(t_1-t_4)-K(t_2-t_3)\color{red}{+}K(t_2-t_4)\right],$$ which should be nonnegative. Note that the last sign $\color{red}{+}$ reads $-$ in your post.

Edit: One asks to compute the second moment of $$X_t=\int_0^t(2s-t)B(s)ds,$$ hence a formula equivalent to the integral in the question is $$\langle X_t^2\rangle=2\int_0^t\mathrm ds(2s-t)\int_0^s\mathrm du(2u-t)K(s-u).$$ When $K(s)=\mathrm e^{-s}$ for every $s\geqslant0$, one finds $$\langle X_t^2\rangle=4-t^2+\tfrac13t^3-(t+2)^2\mathrm e^{-t},$$ which is nonnegative for every $t\geqslant0$, from which one deduces the result for every $\gamma$, by homogeneity.