I want to generalize the answer here to this case of a VAR(1) model. Suppose that $X_t \in \mathbb{R}^n$ and that $\Lambda \in \mathbb{R}^{n \times n}$.
If we have the stochastic process $\left\{ X_t, t = 1, 2, \ldots \right\}$ following the model $$X_t = \Lambda X_{t-1} + e_t \quad \text{ where } e_t \sim \mathcal N (0, \Sigma)$$ what is the distribution of the inital point $X_1$?
Using the idea proposed in the linked question, suppose that $B = \mathbb V (X_t) = \mathbb V (X_{t-1})$. Then, \begin{align} B &= \mathbb V (\Lambda X_{t-1} + e_t) \\ &= \Lambda B \Lambda^{\top} + \Sigma. \end{align}
For a simple example, let $n = 2$, $\Sigma = \begin{bmatrix} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \end{bmatrix}$, $\Lambda = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}$, and $B = \begin{bmatrix} x & y \\ z & w \end{bmatrix}$, where the free elements are non-zero.
Subsituting these in to the "fixed point" equation yields $$\begin{bmatrix} x & y \\ z & w \end{bmatrix} = \begin{bmatrix} \lambda_1^2 x + \sigma_1^2 & \lambda_1 \lambda_2 y \\ \lambda_1 \lambda_2 z & \lambda_2^2 w + \sigma_2^2 \end{bmatrix}$$ which implies that $$x = \frac{\sigma_1^2}{1 - \lambda_1^2} \quad y = \frac{\sigma_2^2}{1 - \lambda_2^2}$$ and $y = z = 0$.
So in this case the starting distribution would be $X_1 \sim \mathcal{N} \left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \frac{\sigma_1^2}{1 - \lambda_1^2} & 0 \\ 0 & \frac{\sigma_2^2}{1 - \lambda_2^2}\end{bmatrix} \right)$ provided $|\lambda_1| < 1$ and $|\lambda_2| <1$.
