Autocorrelation and Partial Autocorrelation

"They ask me what I remember, and I answer at every lag at once. The partial one beside me is ruder: she tells you only what is left after the middlemen have been paid off."

An Autocorrelation Function With a Very Long Memory

The previous section, Section 3.2, fixed the ground we stand on: a weakly stationary process has a constant mean, a constant variance, and an autocovariance $\gamma(k)$ that depends only on the lag $k$ and not on absolute time. That single object, the autocovariance function, is the entire second-order description of such a process. This section normalizes it into correlations and then asks a sharper question than "how related are observations $k$ steps apart?" It asks "how much of that relation is direct, once the chain of intermediate steps is removed?" The answer separates the two great families of linear models, and the rest of Part II depends on being able to read it off a plot. We assume throughout the weak-stationarity of Section 3.2 and the moment notation $\mu$, $\gamma(k)$, $\rho(k)$ collected in the unified notation table of Appendix A.

1. The Autocorrelation Function Beginner

For a weakly stationary process $\{X_t\}$ with mean $\mu$ and autocovariance $\gamma(k) = \operatorname{Cov}(X_t, X_{t+k})$, the autocorrelation function (ACF) at lag $k$ is the autocovariance normalized by the variance,

which is simply the Pearson correlation between the series and a copy of itself shifted by $k$ steps. Three properties follow directly from the definition and are worth stating because every later diagnostic leans on them. First, $\rho(0) = \gamma(0)/\gamma(0) = 1$: a series is perfectly correlated with itself at lag zero. Second, the function is symmetric, $\rho(k) = \rho(-k)$, because $\gamma(k) = \gamma(-k)$ for a stationary process, so we only ever plot $k \ge 0$. Third, it is bounded, $|\rho(k)| \le 1$, a consequence of the Cauchy-Schwarz inequality applied to $X_t$ and $X_{t+k}$. The full collection $\{\rho(k)\}_{k \ge 0}$ plotted against lag is the correlogram, the single most informative picture in classical time series analysis. Figure 3.3.2 offers an intuition for what that picture measures: an echo of the series talking to its own past.

We never know $\rho(k)$ exactly; we estimate it from a finite sample $x_1, \dots, x_n$. The sample autocovariance and the sample ACF are

where $\bar x$ is the sample mean. Two design choices deserve comment. The divisor is $n$ rather than $n-k$, which biases each $\hat\gamma(k)$ slightly toward zero but guarantees that the estimated autocovariance matrix is positive semidefinite, a property later algorithms require; the bias is negligible when $k \ll n$. And because only $n-k$ product terms exist at lag $k$, the estimator grows noisier as $k$ approaches $n$, so the rule of thumb is to trust the ACF only up to roughly $k \le n/4$.

A correlogram without an uncertainty scale is unreadable: every estimate wobbles around its true value, so we need to know which spikes are real. Under the null hypothesis that the process is white noise (all true $\rho(k) = 0$ for $k \ge 1$), Bartlett's result gives the approximate sampling distribution $\hat\rho(k) \sim \mathcal{N}\!\big(0,\, 1/n\big)$ for each lag. The standard error is therefore $1/\sqrt{n}$, and the familiar two-sided 95 percent significance band is

Any sample autocorrelation falling inside this band is statistically indistinguishable from zero; spikes that pierce it are candidates for genuine structure. This band is the horizontal dashed pair you see on every software correlogram, and it is the device that turns a wiggly plot into a hypothesis test. A subtlety worth internalizing is that the band is a per-lag test, not a joint one: across twenty lags under a true white-noise null, the expected number of false breaches is twenty times five percent, which is one. A lone spike that barely clears the band, especially at an interpretively awkward lag, is therefore weak evidence on its own, and the portmanteau tests of subsection four exist precisely to aggregate many lags into a single honest verdict. The numeric example below fixes the scale.

Code 3.3.1 computes the sample ACF and its band from scratch so that the formulas above become executable, with no hidden library behavior.

import numpy as np

def sample_acf(x, max_lag):
    """Sample autocorrelation rho_hat(k) for k = 0..max_lag, divisor n."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    x = x - x.mean()                      # center on the sample mean
    gamma0 = np.dot(x, x) / n             # sample variance gamma_hat(0)
    acf = np.empty(max_lag + 1)
    for k in range(max_lag + 1):
        gamma_k = np.dot(x[:n - k], x[k:]) / n   # divisor n, not n-k
        acf[k] = gamma_k / gamma0
    return acf

rng = np.random.default_rng(0)
n = 200
white = rng.standard_normal(n)            # pure white noise: ACF should be ~0 past lag 0
acf = sample_acf(white, max_lag=10)
band = 1.96 / np.sqrt(n)                  # 95% significance band

print(f"95% band: +/- {band:.3f}")
for k, r in enumerate(acf):
    flag = "  <-- significant" if k > 0 and abs(r) > band else ""
    print(f"  lag {k:2d}: rho_hat = {r:+.3f}{flag}")

95% band: +/- 0.139
  lag  0: rho_hat = +1.000
  lag  1: rho_hat = +0.073
  lag  2: rho_hat = -0.045
  lag  3: rho_hat = -0.011
  lag  4: rho_hat = -0.106
  lag  5: rho_hat = +0.045
  lag  6: rho_hat = +0.030
  lag  7: rho_hat = -0.058
  lag  8: rho_hat = -0.001
  lag  9: rho_hat = +0.087
  lag 10: rho_hat = -0.012

2. The Partial Autocorrelation Function Intermediate

The ACF has a blind spot. If $X_t$ depends directly only on $X_{t-1}$, then $X_t$ and $X_{t-2}$ are still correlated, because both are tied to $X_{t-1}$: the dependence propagates through the chain. The lag-2 autocorrelation is therefore non-zero even though no direct lag-2 mechanism exists. The partial autocorrelation function (PACF) removes exactly this propagated component. The partial autocorrelation at lag $k$, written $\phi_{kk}$, is the correlation between $X_t$ and $X_{t-k}$ after linearly regressing out the intervening variables $X_{t-1}, \dots, X_{t-k+1}$ from both. Formally it is the correlation of the two regression residuals,

where $\hat X_t$ and $\hat X_{t-k}$ are the best linear predictions of $X_t$ and $X_{t-k}$ from the intermediate lags. By construction $\phi_{11} = \rho(1)$ (there is nothing in between to remove), and $\phi_{kk}$ isolates the marginal predictive value that lag $k$ adds once all shorter lags are already used.

There is an equivalent and more computational definition. Fit the autoregression of order $k$ that uses all lags up to $k$,

and read off the last coefficient $\phi_{kk}$: that is the partial autocorrelation at lag $k$. The coefficients $\{\phi_{kj}\}$ solve the Yule-Walker equations, the linear system $\sum_{j=1}^{k} \phi_{kj}\,\rho(i - j) = \rho(i)$ for $i = 1, \dots, k$, whose matrix is the symmetric Toeplitz matrix of autocorrelations. Solving that system afresh for every $k$ is wasteful. The Durbin-Levinson recursion exploits the Toeplitz structure to march from order $k-1$ to order $k$ in $O(k)$ work. Starting from $\phi_{11} = \rho(1)$, it computes the new reflection coefficient as the lag-$k$ correlation left unexplained by the order-$(k-1)$ fit, normalized by the residual variance,

then updates the earlier coefficients by $\phi_{kj} = \phi_{k-1,j} - \phi_{kk}\,\phi_{k-1,\,k-j}$ for $j = 1, \dots, k-1$. The numerator is precisely "the part of $\rho(k)$ the intermediate lags did not already account for," which is why $\phi_{kk}$ deserves the name partial. Under the white-noise null the sample PACF shares the ACF's sampling scale, $\hat\phi_{kk} \sim \mathcal{N}(0, 1/n)$ (Quenouille's result), so the same $\pm 1.96/\sqrt{n}$ band applies to both plots.

Code 3.3.2 implements the Durbin-Levinson recursion directly, taking the sample ACF of Code 3.3.1 as input and returning the PACF.

def pacf_durbin_levinson(acf, max_lag):
    """PACF phi_kk for k = 1..max_lag via the Durbin-Levinson recursion.
    acf must be rho_hat(0..max_lag) with acf[0] == 1.0."""
    phi = np.zeros((max_lag + 1, max_lag + 1))
    pacf = np.zeros(max_lag + 1)
    phi[1, 1] = acf[1]            # phi_11 = rho(1)
    pacf[1] = acf[1]
    for k in range(2, max_lag + 1):
        num = acf[k] - sum(phi[k - 1, j] * acf[k - j] for j in range(1, k))
        den = 1.0 - sum(phi[k - 1, j] * acf[j]     for j in range(1, k))
        phi[k, k] = num / den    # the new reflection coefficient = partial autocorr
        for j in range(1, k):    # update earlier coefficients in place
            phi[k, j] = phi[k - 1, j] - phi[k, k] * phi[k - 1, k - j]
        pacf[k] = phi[k, k]
    return pacf[1:]              # phi_11 .. phi_{max_lag,max_lag}

# Reuse the white-noise ACF from Code 3.3.1.
pacf = pacf_durbin_levinson(acf, max_lag=10)
band = 1.96 / np.sqrt(n)
for k, p in enumerate(pacf, start=1):
    flag = "  <-- significant" if abs(p) > band else ""
    print(f"  lag {k:2d}: phi_kk = {p:+.3f}{flag}")

  lag  1: phi_kk = +0.073
  lag  2: phi_kk = -0.051
  lag  3: phi_kk = -0.004
  lag  4: phi_kk = -0.107
  lag  5: phi_kk = +0.060
  lag  6: phi_kk = +0.018
  lag  7: phi_kk = -0.061
  lag  8: phi_kk = -0.010
  lag  9: phi_kk = +0.084
  lag 10: phi_kk = +0.013

3. Reading ACF and PACF as Fingerprints Intermediate

The practical payoff of the two functions is that the great linear model families leave distinct signatures, and the signatures are nearly mirror images of each other. Consider first an autoregressive process of order $p$, AR($p$), in which the present is a linear combination of the last $p$ values plus noise. Because each new $\phi_{kk}$ measures the direct contribution of lag $k$, and the AR($p$) mechanism has no direct dependence beyond lag $p$, the PACF is exactly zero for $k > p$ and cuts off sharply after lag $p$. The ACF, by contrast, does not cut off: the dependence echoes forward through the recursion and the ACF decays gradually, geometrically for real roots, as a damped sine wave when the characteristic roots are complex.

A moving-average process of order $q$, MA($q$), reverses the picture. Here the present is a weighted sum of the last $q$ shocks, so two observations more than $q$ steps apart share no common shock and are genuinely uncorrelated: the ACF is exactly zero for $k > q$ and cuts off after lag $q$. But writing an MA($q$) as an infinite autoregression shows that every lag contributes a little direct predictive value, so the PACF decays gradually rather than cutting off. A mixed ARMA($p,q$) process inherits the tail of both: neither function cuts off, and both decay. The duality is summarized in the comparison table below, which is the lookup card every analyst keeps within reach when staring at a fresh correlogram.

The mnemonic is worth memorizing: "AR cuts the PACF, MA cuts the ACF, ARMA cuts neither." When exactly one function cuts off you read the order off the cut and the family off which function did the cutting; when neither cuts off you have an ARMA process and defer order selection to the information criteria (AIC, BIC) and automated search of Chapter 5, where this exact reading drives the Box-Jenkins identification step. Real data is messier than the textbook ideal, so in practice these patterns guide a small shortlist of candidate orders rather than naming a single model.

4. Correlogram Diagnostics: Testing Residuals for Whiteness Advanced

The ACF is not only for identifying a model; it is the chief tool for validating one. A correctly specified model extracts all linear temporal structure, leaving residuals that are indistinguishable from white noise. So after fitting, we compute the ACF of the residuals and ask whether any lag breaches the significance band. Eyeballing individual lags invites the multiple-comparison trap of the numeric example in subsection one: with twenty lags, one false alarm is expected even under perfect whiteness. The portmanteau tests aggregate the first $m$ residual autocorrelations into one statistic. The Box-Pierce statistic is $Q_{BP} = n \sum_{k=1}^{m} \hat\rho_e^2(k)$, and the refined Ljung-Box statistic, which has better small-sample behavior, weights each squared autocorrelation by $1/(n-k)$,

where $\hat\rho_e(k)$ is the lag-$k$ autocorrelation of the residual series. Under the null that the residuals are white, $Q_{LB}$ is approximately chi-squared distributed with $m - d$ degrees of freedom, where $d$ is the number of estimated ARMA parameters ($d = p + q$); subtracting $d$ accounts for the autocorrelations already constrained by fitting. A large $Q_{LB}$ (small p-value, conventionally below 0.05) rejects whiteness and tells you the model has left structure on the table. A small $Q_{LB}$ is the green light: the residuals look like the noise the model assumed. Figure 3.3.3 personifies the failure case the test is built to expose.

Code 3.3.3 implements the Ljung-Box test from scratch, reusing the from-scratch ACF, and applies it to a residual series.

from scipy.stats import chi2

def ljung_box(residuals, m, n_params=0):
    """Ljung-Box Q statistic and p-value on the first m residual autocorrelations."""
    n = len(residuals)
    rho = sample_acf(residuals, max_lag=m)[1:]        # drop lag 0
    Q = n * (n + 2) * np.sum(rho**2 / (n - np.arange(1, m + 1)))
    dof = m - n_params                                # subtract fitted ARMA params
    p_value = chi2.sf(Q, dof)                          # upper-tail probability
    return Q, p_value

rng = np.random.default_rng(7)
white_resid = rng.standard_normal(120)               # well-specified -> white residuals
Q, p = ljung_box(white_resid, m=8, n_params=2)
print(f"white residuals : Q = {Q:6.3f}, p = {p:.3f}  -> {'reject' if p < 0.05 else 'pass'}")

# Inject leftover autocorrelation: an undermodeled residual still carries lag-1 structure.
bad = white_resid.copy()
bad[1:] += 0.6 * white_resid[:-1]                     # AR(1)-like leftover the model missed
Q2, p2 = ljung_box(bad, m=8, n_params=2)
print(f"leftover AR(1)  : Q = {Q2:6.3f}, p = {p2:.3f}  -> {'reject' if p2 < 0.05 else 'pass'}")

white residuals : Q =  3.142, p = 0.791  -> pass
leftover AR(1)  : Q = 41.870, p = 0.000  -> reject

5. Worked Example: AR(2) Versus MA(1) Fingerprints Advanced

We now confirm the signatures of the comparison table on simulated data, where the true process is known and the verdict can be checked. Consider an AR(2) process with $\phi_1 = 0.6$, $\phi_2 = -0.3$, whose complex characteristic roots produce a damped oscillation, and an MA(1) process with $\theta_1 = 0.8$. Theory predicts a sharp story. For the AR(2), the PACF should show two significant spikes at lags 1 and 2 and then collapse inside the band, while its ACF decays as a damped sine, alternating sign as it shrinks. For the MA(1), the ACF should show a single significant spike at lag 1 (its theoretical value is $\rho(1) = \theta_1/(1 + \theta_1^2) = 0.8/1.64 \approx 0.488$) and nothing beyond, while its PACF decays geometrically with alternating sign. Figure 3.3.1 schematizes the four correlograms.

Code 3.3.4 simulates both processes, computes their ACF and PACF with the from-scratch routines of subsections one and two, and prints the lags that breach the band. The output confirms Figure 3.3.1 numerically: the AR(2) PACF lights up at lags 1 and 2 only, and the MA(1) ACF lights up at lag 1 only.

def simulate_ar2(n, phi1, phi2, rng, burn=200):
    e = rng.standard_normal(n + burn)
    x = np.zeros(n + burn)
    for t in range(2, n + burn):
        x[t] = phi1 * x[t - 1] + phi2 * x[t - 2] + e[t]
    return x[burn:]

def simulate_ma1(n, theta1, rng, burn=200):
    e = rng.standard_normal(n + burn + 1)
    x = e[1:] + theta1 * e[:-1]          # X_t = e_t + theta1 * e_{t-1}
    return x[burn:]

rng = np.random.default_rng(3)
n, L = 600, 10
ar2 = simulate_ar2(n, 0.6, -0.3, rng)
ma1 = simulate_ma1(n, 0.8, rng)
band = 1.96 / np.sqrt(n)

for name, series in [("AR(2)", ar2), ("MA(1)", ma1)]:
    a = sample_acf(series, L)
    p = pacf_durbin_levinson(a, L)
    acf_sig  = [k for k in range(1, L + 1) if abs(a[k])   > band]
    pacf_sig = [k for k in range(1, L + 1) if abs(p[k-1]) > band]
    print(f"{name}: significant ACF lags  = {acf_sig}")
    print(f"{name}: significant PACF lags = {pacf_sig}")

AR(2): significant ACF lags  = [1, 2, 3, 4]
AR(2): significant PACF lags = [1, 2]
MA(1): significant ACF lags  = [1]
MA(1): significant PACF lags = [1, 2, 3]

Having seen the mechanics, you would never write these routines in production. The statsmodels library computes both functions, with correct bias handling and confidence intervals, in a single call each.

from statsmodels.tsa.stattools import acf, pacf
from statsmodels.stats.diagnostic import acorr_ljungbox

acf_vals  = acf(ar2, nlags=10)                       # sample ACF, lags 0..10
pacf_vals = pacf(ar2, nlags=10, method="ywm")        # Yule-Walker (no sample-size adjustment) PACF
lb = acorr_ljungbox(ar2, lags=[10], return_df=True)  # portmanteau whiteness test

print("ACF :", np.round(acf_vals[1:4], 3))
print("PACF:", np.round(pacf_vals[1:4], 3))
print(lb)