I am trying to compute the asymptotic variance of OU process
$$
d X\_t = - H X\_t dt + S dW\_t
$$
where $X\_t$ takes value in $R^d$. $H$ and $S$ are $d\times d$ matrices that does not have $HS = SH$ in general. How to compute its variance at $t$?


Things are classical for one dimensional case as
$$
Var (X\_t) = \sigma^2 / \mu (1-e^{-\mu t})
$$



Variance of Multi-Dimensional Ornstein Uhlenbeck process
the large-$t$ limit $\rho$ of the variance of $X\_t$ is given by 


$$\rho=\int\_0^\infty \exp(-tH)SS^{\dagger}\exp(-tH^{\dagger})dt$$


where the superscript $\dagger$ indicates the transpose (or the conjugate transpose if the matrix is complex); this result holds also if $S$ is $d\times d'$ with $d'\leq d$; see, for example, page 11 of these [notes](http://www.math.nyu.edu/faculty/goodman/teaching/MonteCarlo07/notes/sde2.pdf), which also show how to evaluate the integral by Schur decomposition of $H$.


