I am considering the following Stochastic dynamics: $$ \mathrm d\mathbf x(t) = \mathbf F\,\mathbf x(t) \mathrm dt + \mathbf L\,\mathrm d\mathbf W(t), \quad \mathbf x(0) = \mathbf x_0 $$ where $\mathbf x(t) = [x(t), v(t), a(t)]^\top,$ $\mathbf W(t) = [W_x(t), W_v(t), W_a(t)]^\top$ is a vector of independent Wiener processes and $$ \mathbf F= \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0 \end{bmatrix}, \qquad \mathbf L = \begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & \sigma_a \end{bmatrix}. $$
When trying to find the distribution of $x(t)$ I get two different values for the variance.
REMARK: Notation $\mathcal N(\mu,\sigma^2)$ used for normal distributions with mean $\mu$ and variance $\sigma^2$.
- Integrating the equations recursively I get: $$ x(t) = x(0) + v(0)t + a(0)\frac{t^2}{2} + \sigma_a\int_0^t\!\left(\int_0^s W_a(u)\mathrm du\right)\!\mathrm ds. $$
Using the Itô formula of integration-by-parts for the inner integral: $$ \int_0^t\!\left(\int_0^s W_a(u)\mathrm du\right)\!\mathrm ds = \int_0^t\!\left(sW_a(s) - \int_0^s u\mathrm dW_a(u)\right)\!\mathrm ds, $$ and doing the same for the first term obtained: $$ \int_0^t sW_a(s)\mathrm ds = \frac{t^2}{2}W_a(t) - \int_0^t \frac{s^2}{2}\mathrm dW_a(s). $$ Now we apply the following: Why can I exchange the order of integration in a multiple Ito stochastic integral? to the second term to get: $$ \int_0^t\!\left(\int_0^s u\mathrm dW_a(u)\right)\!\mathrm ds = \int_0^t\!\left(\int_0^s u\mathrm du\right)\mathrm dW_a(s) = \int_0^t\!\frac{s^2}{2}\mathrm dW_a(s). $$
Gathering all these results: $$ \int_0^t\!\left(\int_0^s W_a(u)\mathrm du\right)\!\mathrm ds = \frac{t^2}{2}W_a(t) - \int_0^t s^2\mathrm dW_a(s). $$
It is known that $\int_0^t f(s)dW(s)\sim\mathcal N\!\left(0,\int_0^t f(s)^2\mathrm ds\right),$ if $f$ is differentiable. So, $$ \int_0^t\!\left(\int_0^s W_a(u)\mathrm du\right)\!\mathrm ds\sim \mathcal N\!\left(0,\frac{t^5}{4}\right) - \mathcal N\!\left(0,\frac{t^5}{5}\right) \sim \mathcal N\!\left(0,\frac{9}{20}t^5\right). $$
- On the other hand, since $\mathbf F$ and $\mathbf L$ are both constant, following Simo Särkkä and Arno Solin Applied Stochastic Differential Equations book: $$ \begin{aligned} \mathbf x(t) &= \exp(\mathbf F(t))\mathbf x_0 + \int_{0}^t \exp(\mathbf F(t-s))\,\mathbf L\, \mathrm d\mathbf W(s)\\[5mm] & = \begin{bmatrix} 1 & t & \frac{t^2}{2}\\ 0 & 1 & t\\ 0 & 0 & 1 \end{bmatrix} \mathbf x_0 + \int_{0}^t \begin{bmatrix} 0 & 0 & \frac{(t-s)^2}{2}\sigma_a\\ 0 & 1 & (t-s)\sigma_a\\ 0 & 0 & \sigma_a \end{bmatrix} \begin{bmatrix} \mathrm d W_x(s)\\ \mathrm d W_v(s)\\ \mathrm d W_a(s) \end{bmatrix}. \end{aligned} $$ Thus, $$ x(t) = \left[1, t, \frac{t^2}{2}\right]\mathbf x_0 \,+\sigma_a\int_{0}^t \frac{(t-s)^2}{2}dW_a(s). $$ In this case, I get: $$ \int_{0}^t \frac{(t-s)^2}{2}dW_a(s)\sim\mathcal N\!\left(0, \frac{t^5}{20}\right). $$
WHERE'S THE MISTAKE?