I am trying to compute the covariance of a 1 dimensional Ornstein-Uhlenbeck process $dx_t=-\theta x_t dt+ \sigma dW_t$, $\theta>0$ and I am at the stage, $$\text{Cov }(x_s,x_t)=\sigma^2 e^{-\theta(t+s)} \mathbb{E}\left[ \int_0^s e^{\theta u}dW_u \int_0^t e^{\theta v} dW_v\right].$$
Is it possible to evaluate the stochastic integrals explicitly and if not how does one go about simplifying this. On Wikipedia, they say this is equal to $\frac{\sigma^2}{2\theta} e^{-\theta(t+s)} (e^{2\theta s\wedge t } -1)$ but I cannot see how they reach this conclusion.
Thanks.

