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

 A: In the context of multivariate time series, there is a mean vector and a variance-covariance matrix. The conditional mean vector and/or the conditional variance-covariance matrix vary over time. Then it is natural to be interested in the conditional correlation rather than unconditional one.
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}) }$.
While the models for the conditional mean of a multivariate time series (e.g. VAR, VARMA, VECM, but also other model classes) might not need an introduction, the models for the conditional variance are perhaps less broadly used and less well known. Therefore you may benefit from a couple of reviews of multivariate GARCH models here and here. For a particular example, you could check out Engle's dynamic conditional correlation (DCC) model. 

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.
