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!

  • 1
    $\begingroup$ You need to look at vector autoregression (VAR) models or similar things such as seemingly unrelated regressions (SUR), vector error correcting (VECM) models and state-space (SSM) models. It's a broad subject. $\endgroup$
    – Aksakal
    Feb 4 '15 at 15:04

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$.

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)
#    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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.