The problem is not the normalisation constant, since in correlation formula it simply cancels out. The difference arises because means and variances of the series are held fixed when calculating the cross-correlations. This means that variance and means are calculated for the whole series, and they are used in calculating correlation when the length of series decreases due to lags. This is a perfectly valid operation if the series are considered stationary, i.e. with constant mean and variance.
Here is the detailed example which recreates the behaviour of ccf:
x = c(1,2,3,4,5,6,7,8,9,10)
y = c(3,3,3,5,5,5,5,7,7,11)
mx <- mean(x)
my <- mean(y)
dx <- mean((x-mx)^2)
dy <- mean((y-my)^2)
nx <- length(x)
round(cor(x,y),3)
[1] 0.896
cr<-function(x,y,mux=mean(x),muy=mean(y),dx=var(x),dy=var(y),n=length(x)) {
cxy<-sum((x-mux)*(y-muy))/n
cxy/sqrt(dx*dy)
}
round(cr(x,y,mx,my,dx,dy,nx),3)
[1] 0.896
# Think "Lag -1"
# x[-10] = 1,2,3,4,5,6,7,8,9
# y[-1] = 3,3,5,5,5,5,7,7,11
round(cor(x[-10],y[-1]),3)
[1] 0.894
round(cr(x[-10],y[-1],mx,my,dx,dy,nx),3)
[1] 0.699
# Think "Lag -2"
# x[-10:-9] = 1,2,3,4,5,6,7,8
# y[-1:-2] = 3,5,5,5,5,7,7,11
round(cor(x[-10:-9],y[-1:-2]),3)
[1] 0.878
round(cr(x[-10:-9],y[-1:-2],mx,my,dx,dy,nx),3)
[1] 0.466
print(ccf(x,y,lag.max=3,plot=FALSE))
Autocorrelations of series ‘X’, by lag
-3 -2 -1 0 1 2 3
0.197 0.466 0.699 0.896 0.436 0.221 -0.018
Note that the norming constant in the function cr is needed only because it must be the same norming constant used in the variance calculations.
The main reason for the "reversal" you are looking at when you deal with AR and MA processes, is that these processes generally have the property that they are invertible to the form of the other process (so long as the coefficients in the models are within the unit circle). So a finite AR process can be represented as an infinite MA process, and a finite MA process can be represented as an infinite AR process. For a general MA(q) process you have:
$$Z_t = \Bigg( 1 - \sum_{i=1}^q \theta_i B^i \Bigg) \epsilon_t = \prod_{i=1}^q (1 - \tau_i B) \epsilon_t,$$
where $B$ is the backshift operator. If $\max|\tau_i| < 1$ (so that all the coefficients are inside the unit circle) then the process is invertible and we have:
$$\epsilon_t = \prod_{i=1}^q (1 - \tau_i B)^{-1} Z_t = \prod_{i=1}^q \Bigg( \sum_{k=0}^\infty \tau_i^k B^k \Bigg) Z_t.$$
Re-arranging this expression gives the AR($\infty$) process:
$$Z_t = \Bigg[ \prod_{i=1}^q \Bigg( \sum_{k=0}^\infty \tau_i^k B^k \Bigg) -1 \Bigg] Z_t + \epsilon_t.$$
Now, the PACF is giving you the conditional correlation for a given lag, conditional on knowledge of the values of the intervening times. For an AR process, this measures the autocorrelations in the process. Hence, for an invertible MA process, the PACF will measure the autocorrelations in the AR($\infty$) process that corresponds to that process. The measured PACF values will decay gradually because the AR process being measured is infinite.
Best Answer
So, after some research on the topic... I came to realise that if you execute the following code:
You get the partial cross correlations between x & y.
So I researched a little, and found this link where Simone Giannerini creates a corrected version of the multivariate pacf. Here is the code:
Which follows the recursive algorithm described in pp. 402-414 in Wei (2005) Time Series Analysis, and actually yields the same output as the first line of code I wrote above (probably the pacf function was corrected on CRAN after he posted this new function)
I think this is what I am looking for, but why isn't it called partial cross correlation, if that's what it (seems) to be?
WARNING: I realised that doing:
Yields highly undesirable results (do not understand what actually happens here, but it is certainly wrong... probably has something to do with the input class not being a ts?), after getting this output: