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

  • $\begingroup$ See pp. 402-414 in Wei (2005) Time Series Analysis, google.no/… $\endgroup$ Jul 19, 2017 at 9:16
  • $\begingroup$ @JarleTufto I have read those pages, but they are actually pretty confusing. Could you please help breaking it down a little for me? $\endgroup$
    – Miguel M.
    Jul 19, 2017 at 9:42

1 Answer 1


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:

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         
#  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;

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:



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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