if $X(t)$ is the  Ornstein-Uhlenbeck process and $Y(t)$ the time integrated OU process I am trying to calculate the autocovariance $cov(Y_t, Y_s)$.

I have a bunch of results but I don't know how to connect them.
First [in this post][1] they treat the OU process and it's also the result of the variance:
$$ Var(Y_t) = \mathbb{E} Y_t^2 = \mathbb{E} \int_0^t \int_0^t X_s X_u ds du =  \int_0^t \int_0^t Cov(X_s, X_u) ds du = 2 \int_0^t \int_0^u Cov(X_s, X_u) ds du.$$

and they use $Cov(X_s, X_u) =  \frac{\sigma^2}{2\theta}\left( e^{-\theta(u-s)} - e^{-\theta(u+s)} \right)$ for $s\leq t$. And in that case the result is $Var(Y_t)=\sigma^2/(2\theta)(1-e^{-2\theta t})$

The thing is that (if I am not wrong) the way to calculate the autocovariance is using a similar integral but with different limits. 
I am not sure about this but I was thinking in $
Cov(Y_t,Y_v)=\mathbb{E} \int_0^t \int_0^v X_s X_u ds du
$
but I am not completely sure.

And in that case I don't know how to solve the integral because I have to respect the validity of the formula of the covariance for $s \leq t$


On the other hand I also found the following result[Bhattacharya] :
$$
Cov(Y_s, Y_{s+t}-Y_s)=\frac{\sigma^2}{2 \theta³}(1-e^{-\theta t})
(1-e^{-\theta s})^2
$$

And also I was thinking in using that
$$
Cov(Y_s, Y_{s+t}-Y_s) = Cov(Y_s, Y_{s+t}) - Var(Y_s)
$$
But again, I am not completely sure about this, and if this is true, there is a way to solve the integral and have the same result.

Any help would be really appreciated.

REF:
Bhattacharya, Stochastic Processes with Applications

  [1]: http://mathoverflow.net/questions/84952/time-integral-of-an-ornstein-uhlenbeck-process