您可以按如下方式定义数据生成过程：

$$\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}$$

R 函数MASS::mvrnorm或mvtnorm::rmvnorm可用于从多元高斯分布生成绘图。

下面的代码从上面定义的模型生成几个系列$p=1$rhos并将每对系列之间的相关性存储在对象中。扰动项的协方差矩阵用值定义$\sigma^2_1=2$,$\sigma^2_2=3$和$\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 

可以查到两者之间的相关性$x_{1,t}$和$x_{2,t}$当两个系列都遵循 AR(1) 过程时，由下式给出：

$$\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}$$

对于示例的参数值，我们有：

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

这与在小型模拟练习中获得的相关性平均值一致。