I have a stochastic process (Ornstein-Uhlenbeck) defined as:
$X(t) = e^{-at}(\int_0^t e^{a \tau} dW(\tau) + X_0)$
Where $W(t)$ is the Wiener process, and $X_0$ is the initial value of my process.
I want to derive the covariance function, which obviously can be found easily online, however the derivations I found skipped a lot of steps without much explanation.
Assuming that $s < t$ then after a few steps I get to:
$EX(t)X(s) = e^{-a(t+s)} E[\int_0^s e^{a\tau}dW(\tau) \int_0^s e^{a\sigma}dW(\sigma) + X_0^2]$
At this point I figure I can combine the integrals and then take the expectation of just $dW(\tau)dW(\sigma)$ since the rest is non-random, and an integral is just an infinite sum. Doing this gets me to:
$EX(t)X(s) = e^{-a(t+s)} (\int_0^s \int_0^s e^{a(\tau+\sigma)} E[dW(\tau) dW(\sigma)] + EX_0^2)$
From this point I figure I can use the fact that for a Wiener process we have that $EdW(t)dW(s)$ is $0$ if $t \neq s$ and is equal to $dt$ if $t=s$. This should then simplify to a single integral which can then be resolved easily.
Is this line of thinking correct? Or is there something I'm missing.
Edit: Finishing it how I described above, I ended up with:
$EX(t)X(s) = \frac{e^{-a(t+s)}}{2a}(e^{2as} -1 + 2a EX_0^2)$