Does the partial version of the Cross Correlation Function exist? If so, how can it be computed?

Question

I have been trying to look for the analogous version of the PACF for cross correlation, but there seems to be no such thing... why is this? Is it not possible to extrapolate the algorithm used to calculate the PACF to two time series? From my understanding, the CCF is merely the correlation between the lags and leads of two time series, as it is defined here: Wikipedia page about cross correlation

Am I missing something about the CCF? Because its computation is very simmilar to how the ACF is constructed... so the "PCCF" should be constructed simmilarly to the PACF? Perhaps re-adapting this algorithm so that the "Zs" are the lagged X&Y variables for each kth lag we wish to correlate?

Also, if anyone could provide some R code to compute the aforementioned "PCCF", it would be much appreciated... I found a "PCCF" function in Matlab, but I am uncertain if it is correct: Matlab CCF function

Answer 1

So, after some research on the topic... I came to realise that if you execute the following code:

pacf(ts(cbind(dx,dy)),lag.max=10)

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:

pacf.mts <- function(x,lag.max){

## Partial autocorrelation function for multivariate time series       
## implements a generalization of the Durbin-Levinson Algorithm        
## as described in                                                     
## Wei(1990) Time Series Analysis, Univariate and Multivariate methods 
##           Addison Wesley.                                           
## Simone Giannerini 2007                                              

#  The author of this software is Simone Giannerini, Copyright (c) 2007   
#  Permission to use, copy, modify, and distribute this software for any  
#  purpose without fee is hereby granted, provided that this entire notice
#  is included in all copies of any software which is or includes a copy  
#  or modification of this software and in all copies of the supporting   
#  documentation for such software. 

#  This program is free software; you can redistribute it and/or modify   
#  it under the terms of the GNU General Public License as published by   
#  the Free Software Foundation; either version 2 of the License, or      
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,        
#  but WITHOUT ANY WARRANTY; without even the implied warranty of         
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/  

nser    <- ncol(x);
snames  <- colnames(x)
alpha.c <- array(0,c(nser, nser,lag.max,lag.max));
beta.c  <- array(0,c(nser, nser,lag.max,lag.max));
P       <- array(0,c(lag.max, nser, nser));
cov.x   <- acf(x,plot=FALSE,type="covariance",lag.max=lag.max)$acf;
Vu      <- array(0,dim=c(nser,nser,lag.max));
Vv      <- Vu;
Vvu     <- Vu;
Vu[,,1]  <- cov.x[1,,]; ## Gamma[0]
Vv[,,1]  <- cov.x[1,,]; ## Gamma[0]
Vvu[,,1] <- cov.x[2,,]; ## Gamma[1]
alpha.c[,,1,1] <- t(Vvu[,,1])%*%solve(Vv[,,1])
beta.c[,,1,1]  <- Vvu[,,1]%*%solve(Vu[,,1])
Du     <- diag(sqrt(diag(Vu[,,1])),nser,nser);
Dv     <- diag(sqrt(diag(Vv[,,1])),nser,nser);
P[1,,] <- solve(Dv)%*%Vvu[,,1]%*%solve(Du);
if(lag.max>=2) {
    for(s in 2:lag.max) {
        dum.u  <- 0;
        dum.v  <- 0;
        dum.vu <- 0;
        for(k in 1:(s-1)) {
            dum.u <- dum.u + alpha.c[,,s-1,k]%*%cov.x[(k+1),,];
            dum.v <- dum.v + beta.c[,,s-1,k]%*%t(cov.x[(k+1),,]);
            dum.vu<- dum.vu + cov.x[(s-k+1),,]%*%t(alpha.c[,,s-1,k]);
        }
        Vu[,,s]  <- cov.x[1,,] - dum.u;
        Vv[,,s]  <- cov.x[1,,] - dum.v;
        Vvu[,,s] <- cov.x[(s+1),,] - dum.vu;
        alpha.c[,,s,s] <- t(Vvu[,,s])%*%solve(Vv[,,s]);
        beta.c[,,s,s]  <- Vvu[,,s]%*%solve(Vu[,,s]);
        for(k in 1:(s-1)) {
            alpha.c[,,s,k] <- alpha.c[,,s-1,k]-alpha.c[,,s,s]%*%beta.c[,,s-1,s-k]
            beta.c[,,s,k]  <- beta.c[,,s-1,k] - beta.c[,,s,s]%*%alpha.c[,,s-1,s-k]
        }
        Du     <- diag(sqrt(diag(Vu[,,s])),nser,nser);
        Dv     <- diag(sqrt(diag(Vv[,,s])),nser,nser);
        P[s,,] <- solve(Dv)%*%Vvu[,,s]%*%solve(Du);
    }
}
colnames(P) <- snames;
return(P)}

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:

pacf(cbind(dx,dy))

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: