# Correlation as a time series

I have been studying a few simple statistical models for (univariate) time series. From my understanding,

• ARIMA and its siblings are used to model the mean of a time series. Rather than a static measure like mean(), the result is a series estimating the mean.

• ARCH and its brothers are used to model the volatility of a time series. Rather than the usual sd(), the result is a series estimating the variance.

## Question

What would be a credible model for the correlation of two time series?

## Notes

While mean models explore the idea of regressing lagged values of the time series, volatility models (eg. ARCH model) explore the idea of regressing lagged residuals where residuals are the difference of a mean model to its original time series.

In its general sense and for a variety of reasons, ARIMA and ARCH are superior models than rolling windows with mean() (popularly known as moving averages outside statistics world) and sd().

However, there is no such a thing for the correlation of two time series X and Y to my knowledge.

The closest thing would be rolling a sad, straight window with cor(), Pearson's coefficient function in R, and work around the resulting series.

## A poor solution

Trying to replicate Pearson's correlation model,

p_(X,Y) = cov(X,Y) / (sd(X) sd(Y))
= E((X-mean(X))(Y-mean(Y))) / (sd(X) sd(Y)),


to the time series world, I had the above without the intended success.

library('forecast')
library('fGarch')

X <- 1:200 + rnorm(200, sd=10)
Y <- 50 + (1:200)/100 + rnorm(200, sd=5)

plot(1:200, X, t='l', main="What would be a resulting ts correlation of X and Y?")
lines(1:200, Y, t='l', col='blue')

# Mimic Pearson correlation, cov(X,Y)/(sd(X)*sd(Y)).
Xm <- as.vector(X) - as.vector(fitted(Arima(X, order=c(2,0,1))))
Ym <- as.vector(Y) - as.vector(fitted(Arima(Y, order=c(2,0,1))))

Xv <- garchFit(formula=~arma(2,1) + garch(2,1), data=X)@sigma.t
Yv <- garchFit(formula=~arma(2,1) + garch(2,1), data=Y)@sigma.t

correlation <- Xm * Ym / (Xv * Yv)    # this can be forecast

plot(correlation, t='l', col='blue', ylim=c(-2, 2), main='Correlation models')
abline(h=c(-1, 1))
abline(h=cor(X, Y), col='red', lwd=5)

# Correlation rolling window of size 10.
df <- data.frame(X, Y)
crw <- rep(NA, 10)
for (i in 11:nrow(df))
crw <- c(crw, cor(df[(i-10):i, 1], df[(i-10):i, 2]))

lines(crw, col='darkgreen', lwd=5)

legend('topright',
c('pearson mimic', 'static cor()', 'rolling cor() like moving averages'),
col=c('blue', 'red', 'darkgreen'), lwd=c(1, 5, 5))

• Have you looked at the cross correlation function ccf? That is a good exploratory approach. I think the corresponding model might be a VAR model as implemented in R package vars. – AlaskaRon Oct 8 '15 at 7:23
• Check out multivariate volatility models with focus on conditional correlations. E.g. Engle's dynamic conditional correlation (DCC) model here. Good general overview of multivariate volatility models are this and this. Note that ARCH-type of conditional variances is not a prerequisite for considering conditional correlations. – Richard Hardy Oct 8 '15 at 7:29
• Your correlation looks weird when plotted because it actually does not match the definition of Pearson correlation that you refer to above. In correlation <- Xm * Ym / (Xv * Yv), the part Xm * Ym should be replaced by (X-Xm) * (Y-Ym). Also, you ignore the presence of the linear trends in X and Y when fitting the ARIMA models; that clearly yields a poor fit. – Richard Hardy Oct 17 '15 at 19:21

Given a model for the conditional mean vector and the conditional variance-covariance matrix, a particular conditional correlation between variables $x_i$ (in your case X) and $x_j$ (in your case Y) at time $t$ can be immediately retrieved as the ratio of the conditional covariance over the product of the conditional standard deviations: $\rho_{i,j,t}=\frac{ \text{Cov}(x_{i,t},x_{j,t}) }{ \text{sd}(x_{i,t}) \cdot \text{sd}(x_{j,t}) }$.
Your correlation looks weird when plotted because it actually does not match the definition of Pearson correlation that you refer to above. In correlation <- Xm * Ym / (Xv * Yv), the part Xm * Ym should be replaced by (X-Xm) * (Y-Ym). Also, you ignore the presence of the linear trends in X and Y when fitting the ARIMA models; that clearly yields a poor fit.