# How to simulate two correlated AR(1) time series?

Is it possible to simulate two time series AR(1) for example 0.5 and 0.8 and at the same time these time series to be correlated with $\rho = 0.8$ using R?

ts.sim <- arima.sim(list(order = c(1,0,0), ar = 0.9), n = 100)
ts.sim2 <- arima.sim(n = 100, list(ar=c(0.9)),n.start = 25)
et<-rnorm(100)
yt<-0.8*ts.sim+et
yt2<-0.8*ts.sim2+et


I am not sure if yt and ts.sim or yt2 and ts.sim2 are indeed correlated.

Note that ts.sim and ts.sim2 are slighty different approaches of the same idea. Just wanted to see if the results would be different.

• What do you mean by "correlated with $\rho =0.8$"? Just that the contemporaneous values of the process are correlated in this way? Or, something else? Do you know about vector autoregressive processes? – cardinal Oct 1 '13 at 13:18
• @cardinal No, it has nothing to do with VAR models at this time at least. See my update. – harris Oct 1 '13 at 22:03
• Welcome to this site! You can use Markdown syntax to highlight your code, and if you want to address a particular person just write @ followed by their username. – chl Oct 1 '13 at 22:12
• @chi Thanks for the suggestions. As you have already unsterstand it is my first question here. – harris Oct 1 '13 at 22:41
• @harris what you've described is a vector AR, even if you didn't think about it that way. – Glen_b Oct 1 '13 at 22:52

As @cardinal and @Glen_b pointed out, it is necessary to simulate VAR model to achieve what you want, i.e. two AR(1) processes which are correlated with specified correlation.

The problem can be stated as follows. Simulate two processes $X_t=\rho_1 X_{t-1}+u_t$ and $Y_t=\rho_2 Y_{t-1}+v_t$ such that $\mathrm{corr}(X_t,Y_t)=\rho$.

We can write

\begin{align} \begin{bmatrix}X_t\\Y_t\end{bmatrix}=\begin{bmatrix}\rho_1 & 0\\0 & \rho_2\end{bmatrix}\begin{bmatrix}X_{t-1}\\Y_{t-1}\end{bmatrix}+\begin{bmatrix}u_t\\v_t\end{bmatrix} \quad (1) \end{align}

This is a VAR(1) process: $Z_t=FZ_{t-1}+\varepsilon_t$, where $Z_t$ and $\varepsilon_t$ are vectors with $F$ being a matrix.

If we denote by $\Sigma$ the covariance matrix of $Z_t$ and $Q$ the covariance matrix of $\varepsilon_t$ then it can be checked that

$$\Sigma=F\Sigma F'+Q$$

This equation needs to be solved and its solution is the following:

$$\mathrm{vec}\Sigma=[I-F\otimes F']^{-1}\mathrm{vec}Q,$$

where $\mathrm{vec}$ stands for the vectorisation and $\otimes$ for the Kronecker product and $I$ is the identity matrix. Since we are dealing with $2\times 2$ matrices it is not hard to write this down precisely:

\begin{align} \begin{bmatrix} \sigma_{11}\\\sigma_{12}\\\sigma_{21}\\\sigma_{22} \end{bmatrix}=\left(\begin{bmatrix}1-\rho_1^2 & 0 & 0 & 0\\ 0 & 1-\rho_1\rho_2 & 0 & 0\\ 0 & 0 & 1-\rho_1\rho_2 & 0\\ 0 & 0 & 0 & 1-\rho_2^2 \end{bmatrix}\right)^{-1} \begin{bmatrix} 1\\ q_{12}\\q_{21}\\1 \end{bmatrix}, \end{align} where for convenience I assumed that error terms have unit variances. We can immediately see that variances of $X_t$ and $Y_t$, respectively $\sigma_{11}$ and $\sigma_{22}$ are what they are supposed to be, i.e. variances of $AR(1)$ processes. Now the desired correlation is

$$\rho=\frac{\sigma_{21}}{\sqrt{\sigma_{11}\sigma_{22}}}=\frac{q_{12}\sqrt{(1-\rho_1^2)(1-\rho_2^2)}}{1-\rho_1\rho_2}$$

So to generate $AR(1)$ with desired correlation $\rho$ we need to generate $u_t$ and $v_t$ with unit variances and correlation

$$\mathrm{corr}(u_t,v_t)=\rho\frac{1-\rho_1\rho_2}{\sqrt{(1-\rho_1^2)(1-\rho_2^2)}}$$

Let us illustrate this with the following code

> set.seed(123)
> calcrho<-function(rho,rho1,rho2) {
+   rho*(1-rho1*rho2)/sqrt((1-rho1^2)*(1-rho2^2))
+ }
>
> burn.in<-300
> n<-1000
> rho<-0.8
> rho1<-0.5
> rho2<-0.7
> q12<-calcrho(rho,rho1,rho2)
> eps<-mvrnorm(n+burn.in,mu=c(0,0),Sigma=cbind(c(1,q12),c(q12,1)))
>
> x<-arima.sim(list(ar=rho1),n,innov=eps[burn.in+1:n,1],start.innov=eps[1:burn.in,1])
> y<-arima.sim(list(ar=rho2),n,innov=eps[burn.in+1:n,2],start.innov=eps[1:burn.in,2])
> auto.arima(x)
Series: x
ARIMA(1,0,0) with zero mean

Coefficients:
ar1
0.4695
s.e.  0.0279

sigma^2 estimated as 1.009:  log likelihood=-1423.51
AIC=2851.02   AICc=2851.03   BIC=2860.83
> auto.arima(y)
Series: y
ARIMA(1,0,0) with zero mean

Coefficients:
ar1
0.6968
s.e.  0.0227

sigma^2 estimated as 1.003:  log likelihood=-1420.96
AIC=2845.92   AICc=2845.93   BIC=2855.73

> cor(x,y)
[1] 0.7968937


As we see the code confirms what we wanted to achieve.

• for n = 100 and rho1 = rho1 = 0.95 then the burn.in <- 30 gives me an error! Do you have any idea why??? R says that I have to give a value of 108, which is bigger than the sample size! Strange! – harris Jan 2 '14 at 20:20
• @mpiktas . what if I want to generalize this to N time series (your answer here is for N=2), with given cross correlation between the series – mohamed Jul 28 '20 at 13:30
• @mohamed in the same way. The formula with matrices $F$, $\Sigma$ and $Q$ applies for any case. Note it is really easy to solve as the inverted matrix will be diagonal for all $N$ – mpiktas Aug 10 '20 at 19:32