Given a pos. def. covariance matrix $\Sigma\in \mathbb{R}^{d\times d}$ and first order autocorrelations $a\in (-1,1)^d$ for each dimension, I would like to generate Gaussian random variables in time $X_t \sim N(0,\Sigma)$ such that in each dimension $i$ the lag-1-autocorrelation $a_i$ is met. Is there a simple characterization which combinations of $\Sigma$ and $a$ are valid? How could I generate such a process efficiently?
I'd prefer a solution using the following VAR(1) model: $$ X_t = A \cdot X_{t-1} + \varepsilon_t, $$ where $\varepsilon_t \sim N(0,\Sigma_\varepsilon)$ iid and $A= \mathrm{diag}(a)$.
Then by $X_t = \sum_{k=0}^\infty A^k \varepsilon_{t-k}$, it would hold that $$X_t\sim N\left(0, \sum_{k=0}^\infty A^k \Sigma_\varepsilon A^k \right)= N\left(0,\Bigl(\frac{(\Sigma_{\varepsilon})_{ij}}{1-a_i a_j} \Bigr)_{ij}\right).~$$ Now if the next step is well-defined, it holds that $(\Sigma_\varepsilon)_{ij} = \Sigma_{ij} (1-a_i a_j)$. Unfortunately $\Sigma$ and $a$ may be incompatible such that the resulting $\Sigma_\varepsilon$ is not pos. semidefinite.
So can anyone characterize the conditions on $\Sigma$ and $A$ (or better on $\Sigma$ and $a$) such that there exists a positive definite covariance matrix $\Sigma_\varepsilon$ that allows me to generate the corresponding VAR-process?
For example if $A$ and $\Sigma$ share the same eigenvectors then there exists a pos. semidefinite matrix $\Sigma_\varepsilon$ for the VAR(1)-process above. But my particular setting is the covariance matrix of $d$ points on the sphere $S^2$ of an isotropic random field, hence the eigenvector basis will not be the standard basis and $A$ not diagonal (as preferred).
Another idea, instead of VAR(1), is to use a spatio-temporal covariance function for $d$ grid points and time such that the autocorrelation depends on the location. Which combinations of $\Sigma$ and $a$ are allowed now and how do I find a plausible covariance function?