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.
