plot_partial_cross_correlation

plot_partial_cross_correlation.py

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.regression.linear_model import OLS
from statsmodels.tools.tools import add_constant

def plot_partial_cross_correlation(x, y, max_lag=20, title=None):
    """
    Plot partial cross-correlation (PCC) between time series x and y.
    
    Parameters:
        x, y     : 1D arrays (should be same length)
        max_lag  : max lag to calculate
        title    : optional title for the plot
    """
    x = np.asarray(x)
    y = np.asarray(y)
    n = len(x)
    
    if len(y) != n:
        raise ValueError("x and y must have the same length.")
    
    # Standardize
    x = (x - x.mean()) / x.std()
    y = (y - y.mean()) / y.std()
    
    pcc_vals = []

    for lag in range(1, max_lag + 1):
        x_t = x[lag:]
        y_lagged = y[:-lag]
        if lag > 1:
            y_intermediate = np.column_stack([y[lag-i:-i] for i in range(1, lag)])
            x_resid = OLS(x_t, add_constant(y_intermediate)).fit().resid
            y_resid = OLS(y_lagged, add_constant(y_intermediate)).fit().resid
        else:
            x_resid = x_t
            y_resid = y_lagged

        pcc = np.corrcoef(x_resid, y_resid)[0, 1]
        pcc_vals.append(pcc)

    # Plot
    lags = np.arange(1, max_lag + 1)
    plt.figure(figsize=(10, 5))
    plt.stem(lags, pcc_vals)
    plt.axhline(0, color='black')
    conf_bound = 2 / np.sqrt(n)
    plt.axhline(conf_bound, color='red', linestyle='--')
    plt.axhline(-conf_bound, color='red', linestyle='--')
    plt.xlabel("Lag")
    plt.xticks(lags)  # or use np.arange(min, max + 1)
    plt.ylabel("Partial Cross-Correlation")
    plt.title(title or "Partial Cross-Correlation (PCC)")
np.random.seed(0)
n = 20
x = np.sin(np.arange(n))+ np.random.randn(n) * 0.5
y = np.roll(x, 3) + np.random.randn(n) * 0.5  # Y lags X by 3
plt.plot(x, label='x')
plt.plot(y, label='y')
plt.legend()
plot_partial_cross_correlation(y, x, max_lag=15, title="X vs Y PCC")