I have been working with the Ornstein-Uhlenbeck stochastic differential equation, given by: $dX_t= -QX_t dt+ \epsilon dB_t$ where $B_t$ is the d-dimensional Brownian Motion and $Q$ is a $d\times d$ deterministic matrix.
Following Itô's formula, I know that there's an unique solution to this problem given by: $X_t=e^{-Qt}X_0 + \epsilon \int_0^t e^{-Q(t-s)}dB_s = e^{-Qt}X_0 + \epsilon \mathcal{O}_t$ where $X_0$ is the initial condition and I'm interested in finding where does $\mathcal{O}_t$ converges to in distribution.
Now, I know that $\mathcal{O}_t=\int_0^t e^{-Q(t-s)} dB_s$ is a Wiener's integral and so, $\mathcal{O}_t\sim \mathcal{N}(0, \int_0^t e^{-2Q(t-s)} ds)$.
Most of my results are for $d=1$, and $Q=\sigma \in \mathbb{R}$:
I know $\mathcal{O}_t$'s moment generating function and characteristic function for $d=1$.
and I've shown that $\mathcal{O}_t \to \mathcal{O}_\infty$ in distribution for $\mathcal{O}_\infty \sim \mathcal{N}(0, 1/2\sigma)$ since $M_n(s) \to M(s)$, where $M_n(s)$ and $M(s)$ are the moment generating functions of $\mathcal{O}_t$ and $\mathcal{O}_\infty$ respectively.
I want to generalize to the multivariate case (assuming just that the real part of the eigenvalues of $Q$ are positive). However, I've been struggling trying to find the characteristic function of $\mathcal{O}_t$. I wanted to use a similar argument using $M_n(s)$ and $M(s)$ for the multivariate case. But, my teacher told me that I should use the characteristic function instead. Why is that?
Also, I'll appreciate any hint (or reference) on how to determine the multivariate characteristic function of $\mathcal{O}_t$.