Para cada $t\geqslant0$, vamos a $Y_t=\int\limits_0^tX_s\mathrm ds$. Según lo explicado por @Sasha, $Z_t=(X_t,Y_t)$ es gaussiano, por tanto los valores
$$
x(t)=\mathrm E(X_t),\quad
y(t)=\mathrm E(Y_t),\quad
u(t)=\mathrm E(X_t^2),\quad
v(t)=\mathrm E(Y_t^2),\quad
w(t)=\mathrm E(X_tY_t),
$$
totalmente determinar su distribución. Para calcular estos, un método conveniente es considerar directamente las ecuaciones diferenciales ordinarias que las expectativas de las soluciones de ecuaciones diferenciales estocásticas resolver, de la manera siguiente.
En primer lugar, $(X_t)$ $(Y_t)$ resolver el sistema $X_0=x_0$, $Y_0=0$, y
$$
\mathrm dX_t=-aX_t\mathrm dt+\sigma\mathrm dW_t,\quad
\mathrm dY_t=X_t\mathrm dt,
$$
y $(W_t)$ es una martingala de ahí
$$
x'(t)=-ax(t),\quad
y'(t)=x(t).
$$
Ya se sabe que $x(0)=x_0$$y(0)=0$, este diferencial lineal del sistema determina $x(t)$ $y(t)$ por cada $t\geqslant0$.
En cuanto al segundo cantidades de la orden, primero vamos a recordar una forma general de la fórmula de Itô: para cada $n\geqslant1$, todos los reales valores de la función $\varphi$ definido en (un conjunto abierto de) $\mathbb R^n$, y cada una de las $n$ dimensiones de difusión de $(\xi_t)$,
$$
\mathrm d\varphi(\xi_t)=\mathrm{grad}\, \varphi(\xi_t)\cdot\mathrm d\xi_t+\frac12H(\varphi)(\xi_t)\cdot\mathrm d\langle\xi\xi\rangle_t,
$$
donde $\mathrm{grad}\, \varphi(\xi)$ es el gradiente de $\varphi$ $\xi=(\xi^i)_{1\leqslant i\leqslant n}$ por lo tanto
$$
\mathrm{grad}\, \varphi(\xi_t)\cdot\mathrm d\xi_t=\sum\limits_{i=1}^n\partial_i\varphi(\xi_t)\mathrm d\xi^i_t,
$$
y $H(\varphi)(\xi)$ es la matriz Hessiana de $\varphi$ $\xi$ por lo tanto
$$
H(\varphi)(\xi_t)\cdot\mathrm d\langle\xi\xi\rangle_t=\sum\limits_{i=1}^n\sum\limits_{j=1}^n\partial^2_{ij}\varphi(\xi_t)\mathrm d\langle\xi^i,\xi^j\rangle_t.
$$
Primera aplicación: para $\varphi(\xi)=(\xi)^2$ y $\xi_t=X_t$, $\mathrm d\langle X,X\rangle_t=\sigma^2\mathrm dt$ los rendimientos
$$
\mathrm d(X_t^2)=2X_t\mathrm dX_t+\mathrm d\langle X,X\rangle_t=-2aX_t^2\mathrm dt+\mathrm d\text{(MG)}_t+\sigma^2\mathrm dt,
$$
donde $\mathrm d\text{(MG)}_t$ es una martingala término cuyo valor es irrelevante. Por lo tanto
$$
u'(t)=-2au(t)+\sigma^2.
$$
Segunda aplicación: para $\varphi(\xi)=(\xi)^2$ y $\xi_t=Y_t$, $\mathrm d\langle Y,Y\rangle_t=0$ los rendimientos
$$
\mathrm d(Y_t^2)=2Y_t\mathrm dY_t=2Y_tX_t\mathrm dt,
$$
por lo tanto
$$
v'(t)=2w(t).
$$
Tercera aplicación: para $\varphi(\xi^1,\xi^2)=\xi^1\xi^2$ y $\xi_t=(X_t,Y_t)$, $\mathrm d\langle X,Y\rangle_t=0$ los rendimientos
$$
\mathrm d(X_tY_t)=X_t\mathrm dY_t+Y_t\mathrm dX_t=X_t^2\mathrm dt-aX_tY_t\mathrm dt+\mathrm d\text{(MG)}_t,
$$
por lo tanto
$$
w'(t)=u(t)-aw(t).
$$
Las tres ecuaciones diferenciales escrito anteriormente son de primer orden lineal del sistema en $(u(t),v(t),w(t))$. Desde $(u(0),v(0),w(0))=(x_0^2,0,0)$, los valores de $u(t)=\mathrm E(X_t^2)$, $v(t)=\mathrm E(Y_t^2)$ y $w(t)=\mathrm E(X_tY_t)$ seguir, por todos los $t\geqslant0$.
Para responder directamente el OP pregunta, uno puede también tener en cuenta las expectativas de $x(t)=\mathrm E(X_t)$ $y(t)=\mathrm E(Y_t)$ y escribir el sistema diferencial satisfecho por los elementos de la matriz de covarianza $C(t)=\mathrm E((X_t,Y_t)^T(X_t,Y_t))-(\mathrm E(X_t),\mathrm E(Y_t))^T(\mathrm E(X_t),\mathrm E(Y_t))$, es decir,
$$
U(t)=\mathrm{var}(X_t)=u(t)-x(t)^2,\quad
V(t)=\mathrm{var}(Y_t)=v(t)-y(t)^2,
$$
y
$$
W(t)=\mathrm{cov}(X_t,Y_t)=w(t)-x(t)y(t).
$$
Uno se
$$
U'(t)=-2aU(t)+\sigma^2,\quad
V'(t)=2W(t),\quad
W'(t)=U(t)-aW(t).
$$
Junto con la condición inicial de que $U(0)=V(0)=W(0)=0$, esta totalmente determina la matriz de covarianza $C(t)=\begin{pmatrix}U(t) & W(t)\\ W(t) & V(t)\end{pmatrix}$ por cada $t\geqslant0$.