Generating cross-correlated stochastic processes I am looking for a robust way to represent and generate multiple stochastic processes that contain time and cross-correlations i.e. I am looking at stochastic processes $X_t^{1}$, $X_t^{2}$, $\ldots$, $X_t^n$ with stationary cross-correlation:
\begin{equation}
\mathrm{Cov}(X_i^t, X_j^{t'}) = \mathrm{Cov}(X_i^0, X_j^{t'-t}).
\end{equation}
I know how to generate two cross-correlated random variables. I know that for a single stochastic process, Karhunen Loeve expansion provides an optimal representation and for Gaussian process can be generated using sequence of uncorrelated normal random variable. I know that MATLAB implements now mutli-ouput AR-type filters arx and armax, but I am afraid I am operating at the edge of their applicability.
My question is how to generate $\hat X_t^{1}$, $\hat X_t^{2}$, $\ldots$, $\hat X_t^n$ so that the correlations are preserved. I would appreciate any ideas, hints or links to the literature or papers that describe it in more detail.
 A: This question is two years old but I will offer an answer anyway. If two stationary processes $X(t)$ and $Y(t)$ are cross-correlated then there are three correlation functions $\rho_X(\tau)$, $\rho_Y(\tau)$ and the cross-correlation function $\rho_{XY}(\tau)$. The covariance functions are given by:
\begin{equation}
 C_X(t,t') = \sigma_X^2\rho_Y(\tau)
\end{equation}
\begin{equation}
 C_Y(t,t') = \sigma_Y^2\rho_Y(\tau)
\end{equation}
\begin{equation}
 C_{XY}(t,t') = \sigma_X\sigma_Y\rho_{XY}(\tau)
\end{equation}
where $\sigma_X$ and $\sigma_Y$ are the standard deviation of $X(t)$ and $Y(t)$ respectively. In the case that these processes are discretised into a finite number of points $\boldsymbol{X} = \{X_1,X_2,\dots,X_n\}$ and $\boldsymbol{Y} = \{Y_1,Y_2,\dots,Y_n\}$ then the corresponding correlation matrices are $\boldsymbol{\rho}_X$, $\boldsymbol{\rho}_Y$ and $\boldsymbol{\rho}_{XY}$.
The covariance structure can be expressed as:
\begin{equation}
 \boldsymbol{\rho} =
 \begin{bmatrix}
 \boldsymbol{\rho}_Y & \boldsymbol{\rho}_{YZ}\\
 \boldsymbol{\rho}_{ZY} & \boldsymbol{\rho}_{Z}\\
 \end{bmatrix}
\end{equation}
A simulation method that uses the covariance matrix could then be adopted such as covariance matrix decomposition:
$$
\begin{bmatrix}
\boldsymbol{X'}\\
\boldsymbol{Y'}\\
\end{bmatrix}
= \boldsymbol{\psi} \boldsymbol{\xi}
$$
where $\boldsymbol{\xi}$ is a vector of $2n$ uncorrelated mean-zero, unit variance normal random variables and $\boldsymbol{\psi}$ is a $2n$ by $2n$ matrix that satisfies the following decomposition of the correlation matrix $\boldsymbol{\rho}$:
\begin{equation}
\boldsymbol{\rho} = \boldsymbol{\psi} \boldsymbol{\psi}^T
\end{equation}
This decomposition can be achieved through either Cholesky or Eigen decomposition. The realisations $\boldsymbol{X'}$ and $\boldsymbol{Y'}$ are mean-zero, unit-variance Gaussian processes with the desired cross-correlation structure $\boldsymbol{\rho}$. They can then be transformed to the desired mean and variance through:
$$\boldsymbol{X} = \boldsymbol{X}'\sigma_X + \mu_X$$
$$\boldsymbol{Y} = \boldsymbol{Y}'\sigma_Y + \mu_Y$$
where $\mu_X$ and $\mu_Y$ are the desired means of $X(t)$ and $Y(t)$ respectively.
