How do you simulate two correlated AR(p) time series?

Question

I would be interested in the mathematical framework plus code in R if possible.

Basically I want to find out the parameters of the two AR(p) models if I already specificed a certain cross-correlation value between the two time series plus their respective autocorrelations.

Thanks in advance!

Answer 1

You could define the data generating process as follows:

\begin{equation} x_{1,t} = \phi_{11} x_{1,t-1} + \cdots \phi_{1p} x_{1,t-p} + \epsilon_{1,t} \,, \quad \hbox{NID}(0, \sigma^2_1) \\ x_{2,t} = \phi_{21} x_{2,t-1} + \cdots \phi_{2p} x_{2,t-p} + \epsilon_{2,t} \,, \quad \hbox{NID}(0, \sigma^2_2) \\ \hbox{Cov}(\epsilon_{1,t}, \epsilon_{2,s}) = \sigma \, \hbox{ if } t=s \hbox{ and zero otherwise.} \end{equation}

The R functions MASS::mvrnorm or mvtnorm::rmvnorm can be used to generate draws from the multivariate Gaussian distribution.

The code below generates several series from the model defined above with $p=1$ and stores in object rhos the correlation between each pair of series. The covariance matrix of the disturbance terms is defined with values $\sigma^2_1=2$, $\sigma^2_2=3$ and $\sigma=0.8$.

require("mvtnorm")
set.seed(123)
sigma <- diag(c(2,3))    
sigma[1,2] <- sigma[2,1] <- 0.8
n <- 200
niter <- 1000
rhos <- rep(NA, niter)
for (i in seq_len(niter))
{
  eps <- mvtnorm::rmvnorm(n = n, mean = rep(0, 2), sigma = sigma)
  x1 <- arima.sim(n = n, model = list(ar = c(0.7)), innov = eps[,1])
  x2 <- arima.sim(n = n, model = list(ar = c(-0.4)), innov = eps[,2])  
  rhos[i] <- cor(x1, x2)
}
summary(rhos)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.00732 0.13480 0.16990 0.16690 0.19990 0.31280 

It can be checked that the correlation between $x_{1,t}$ and $x_{2,t}$ when both series follow an AR(1) process is given by:

\begin{equation} \rho = \frac{\sigma/(1 - \phi_{11}\phi_{21})}{\sqrt{\frac{\sigma^2_1}{1-\phi_{11}^2}}\sqrt{\frac{\sigma^2_2}{1-\phi_{21}^2}}} \end{equation}

For the parameter values of the example we have:

sd1 <- sqrt(2 / (1 - 0.7^2))
sd2 <- sqrt(3 / (1 - 0.4^2))
(0.8/(1 + (0.7*0.4))) / (sd1 * sd2)
# [1] 0.1670049

which is in agreement with the average value obtained for the correlation in the small simulation exercise.