Given a $K\times K$ real matrix $\mathbf{\Phi}$ and given a sequence $\boldsymbol\varepsilon_t$ of multivariate normal variables $\boldsymbol\varepsilon_t\sim \textrm{N}\left(0,\mathbf{\Sigma}\right)$, with $\boldsymbol\Sigma$ a $K\times K$ positive definite symmetric real matrix, which is the stationary distrinbution of
$$ \mathbf{x}_t =\boldsymbol{\varepsilon}_t + \boldsymbol\Phi\,\mathbf{x}_{t-1} = \sum_{h=0}^{t}\boldsymbol\Phi^h\,\boldsymbol\varepsilon_{t-h}\qquad? $$
ps= It is assumed that all the eignevalues of $\boldsymbol\Phi$ belong to the interval $(-1,1)$.