Partial Autocorrelation (PACF) and AR Model Order Identification: Yule-Walker and Levinson-Durbin in Python

Introduction

Computing the autocorrelation function (ACF) tells you, lag by lag, how similar a time series is to its own past. But the ACF alone cannot answer this question: “Is the correlation at lag 2 a genuinely direct effect of lag 2, or is it just the lag-1 influence propagating over two steps?”

The tool that answers it is the partial autocorrelation function (PACF). Because the PACF measures the “pure” correlation after removing the influence of intermediate lags, the order \(p\) of an AR model shows up directly as the largest lag at which the PACF is nonzero. This property is the theoretical justification for staring at ACF/PACF plots when choosing orders for an ARIMA model.

This article is the advanced companion to the autocorrelation basics post. We derive the PACF definition (correlation of linear-prediction residuals), the Yule-Walker equations, and the Levinson-Durbin recursion; organize the AR/MA order-identification decision table; and review the AIC/BIC information criteria. Then, using statsmodels (pacf, plot_acf/plot_pacf, levinson_durbin, AutoReg, ar_select_order, arma_order_select_ic, ARIMA), we run numerical experiments on a known AR(2) process and confirm that the “PACF cutoff” and “AIC/BIC minimization” point to the same order.

Definition of the Partial Autocorrelation (PACF)

As a correlation of linear-prediction residuals

For a stationary time series \(\{x_t\}\) , the partial autocorrelation at lag \(k\) is defined as the correlation remaining after the linear effect of the intermediate values \(x_{t-1}, \ldots, x_{t-k+1}\) has been removed from both sides:

Here \(\hat{x}_t^{(k-1)}\) is the best linear prediction of \(x_t\) from \(x_{t-1}, \ldots, x_{t-k+1}\) , and \(\hat{x}_{t-k}^{(k-1)}\) is the best (backward) linear prediction of \(x_{t-k}\) from the same intermediate values. In other words, the PACF is the correlation between the residuals left over after subtracting whatever the intermediate lags can explain—the time-series analogue of the partial correlation coefficient in regression analysis.

As the last coefficient of an AR regression

There is an equivalent second definition. Fit an autoregression of order \(k\) ,

then the last coefficient \(\phi_{kk}\) is exactly the partial autocorrelation at lag \(k\) :

It measures “how much additional explanatory power \(x_{t-k}\) has even when lags \(1\) through \(k-1\) are already in the regression,” so intuitively it is the same quantity as equation \((1)\) . Note that \(\alpha(1)\) equals the ordinary autocorrelation \(\rho_1\) , since there are no intermediate lags to remove.

Why the PACF cuts off for AR processes

Consider an AR(1) process \(x_t = \phi x_{t-1} + \varepsilon_t\) . Its ACF decays exponentially as \(\rho_k = \phi^k\) , so even at lag 2 we have \(\rho_2 = \phi^2 \neq 0\) . But that correlation is merely an indirect effect propagated along “\(x_{t-2} \to x_{t-1} \to x_t\) .” Once the effect of \(x_{t-1}\) is removed, no direct contribution of \(x_{t-2}\) remains, hence \(\alpha(2) = 0\) .

In general, for an AR(p) process \(x_t\) is determined by its \(p\) most recent values only (\(\phi_{kk} = 0\) in equation \((2)\) for \(k > p\) ), so

The PACF cuts off after lag \(p\) —this is the heart of AR order identification.

The Yule-Walker Equations

Computing the PACF in practice requires the coefficients \(\phi_{kj}\) of equation \((2)\) . For the AR(p) process

multiply both sides by \(x_{t-k}\) (\(k \geq 1\) ) and take expectations. Since \(E[\varepsilon_t x_{t-k}] = 0\) , we obtain a linear system in the autocorrelations \(\rho_k\) :

These are the Yule-Walker equations. Collecting \(k = 1, \ldots, p\) in matrix form,

The matrix on the left is a Toeplitz matrix, constant along its diagonals. Plugging the sample ACF into equation \((7)\) and solving gives the Yule-Walker estimates of the AR coefficients; solving it for increasing orders \(k = 1, 2, 3, \ldots\) and collecting the last coefficient \(\phi_{kk}\) at each order yields the sample PACF.

The Levinson-Durbin Recursion

Solving equation \((7)\) independently at each order \(k\) repeats an \(O(k^3)\) matrix inversion. The Levinson-Durbin recursion exploits the Toeplitz structure to update the order-\((k-1)\) solution to the order-\(k\) solution in \(O(k)\) .

Initialize with \(\phi_{11} = \rho_1\) and \(\sigma_1^2 = \gamma_0(1 - \rho_1^2)\) , then iterate for \(k = 2, 3, \ldots\) :

Here \(\sigma_k^2\) is the residual (innovation) variance of the order-\(k\) linear predictor. Each of the three equations has a clear meaning:

• Equation \((8)\) : the numerator is “the lag-\(k\) correlation the current model fails to explain,” the denominator a normalization. This \(\phi_{kk}\) is called the reflection coefficient and is precisely the PACF
• Equation \((9)\) : the correction applied to the existing coefficients when the order is raised by one
• Equation \((10)\) : as long as \(|\phi_{kk}| < 1\) , the residual variance decreases monotonically, and once the true order \(p\) is exceeded, \(\phi_{kk} \approx 0\) and the decrease stalls

Running the recursion up to a maximum order \(p\) costs \(O(p^2)\) in total and delivers the AR coefficients and PACF for every order at once. The same mathematics underlies linear predictive coding (LPC) in speech processing and lattice filter design, making this recursion a hinge between time-series analysis and signal processing.

The Order-Identification Decision Table

For AR and MA processes, the roles of ACF and PACF are exactly swapped. An MA(q) process \(x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}\) shares no noise terms with observations more than \(q\) steps away, so its ACF cuts off after lag \(q\) . Conversely, rewriting an MA process in AR form requires infinite order, so its PACF decays without cutting off.

A confidence band decides what counts as a “cutoff.” At lags where the true value is zero, the sample PACF is asymptotically \(\mathcal{N}(0, 1/N)\) , so we call

significant (95% confidence band). Note that about 5% of lags will exceed this band by pure chance, so in practice isolated significant values at high lags are ignored.

Order Selection with Information Criteria (AIC / BIC)

Judging cutoffs from ACF/PACF plots is visual and intuitive, but subjective when the cutoff is ambiguous (mixed ARMA processes, small samples). Information criteria complement it. With \(k\) parameters, sample size \(N\) , and maximized log-likelihood \(\ln \hat{L}\) ,

are computed for each candidate order and the minimizer is chosen. They trade off goodness of fit (first term) against complexity (second term); since \(\ln N > 2\) whenever \(N \geq 8\) , the BIC penalty is heavier and BIC tends to select smaller orders than AIC. Theoretically, BIC is consistent (it selects the true order with probability approaching 1 as \(N \to \infty\) when the true model is among the candidates), while AIC is efficient in the sense of minimizing prediction error. In practice, compute both: if they agree, take that order; if they disagree, a good rule of thumb is AIC for forecasting and BIC for structural interpretation.

One caveat: model comparison must use likelihoods over the same data points. An AR(p) model consumes its first \(p\) observations as initial values, so naively fitting different orders changes the sample. ar_select_order, used below, aligns the samples for a fair comparison.

Python Experiments on an AR(2) Process

Generating the AR(2) and computing ACF/PACF

We generate an AR(2) process with known truth, \(x_t = 0.7 x_{t-1} - 0.3 x_{t-2} + \varepsilon_t\) (\(\varepsilon_t \sim \mathcal{N}(0, 1)\) , \(N = 500\) ), and confirm that the PACF cuts off after lag 2.

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

np.random.seed(42)
N = 500
phi = np.array([0.7, -0.3])      # true AR coefficients
ar = np.r_[1, -phi]              # statsmodels sign convention: [1, -phi1, -phi2]
ma = np.r_[1]
x = ArmaProcess(ar, ma).generate_sample(nsample=N)

nlags = 20
acf_vals = acf(x, nlags=nlags, fft=True)
pacf_vals = pacf(x, nlags=nlags, method="ywadjusted")
ci = 1.96 / np.sqrt(N)           # 95% confidence band

print(f"95% confidence band: +/-{ci:.3f}")
print(f"{'lag':>4} {'ACF':>8} {'PACF':>8}")
for k in range(1, 6):
    line = f"{k:>4} {acf_vals[k]:>8.3f} {pacf_vals[k]:>8.3f}"
    if abs(pacf_vals[k]) > ci:
        line += " *"
    print(line)

# the classic visualization: plot_acf / plot_pacf
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(x, lags=nlags, ax=axes[0], title="ACF: gradual decay")
plot_pacf(x, lags=nlags, ax=axes[1], method="ywm",
          title="PACF: cuts off after lag 2")
plt.tight_layout()
plt.show()

Output:

95% confidence band: +/-0.088
 lag      ACF     PACF
   1    0.532    0.533 *
   2    0.060   -0.312 *
   3   -0.136   -0.020
   4   -0.135   -0.020
   5   -0.013    0.079

Exactly as the decision table predicts. The ACF decays with alternating signs, \(0.532 \to 0.060 \to -0.136\) (oscillating decay caused by \(\phi_2 < 0\) ), showing no clear cutoff. The PACF, however, has significant values \(\hat{\alpha}(1) = 0.533\) and \(\hat{\alpha}(2) = -0.312\) , while \(\hat{\alpha}(3)\) onward (\(-0.020, -0.020, 0.079\) ) all fall inside the confidence band \(\pm 0.088\) —a clean cutoff at lag 2, identifying the process as AR(2). Also note that \(\hat{\alpha}(2) = -0.312\) is itself an estimate of the true \(\phi_2 = -0.3\) (equation \((3)\) gives \(\phi_{pp} = \phi_p\) at \(k = p\) ).

Checking the PACF and innovation variance via Levinson-Durbin

levinson_durbin returns the recursion’s by-products—the PACF and innovation variance at every order—in one call.

from statsmodels.tsa.stattools import levinson_durbin

sigma_v, arcoefs, pacf_ld, sigma, _ = levinson_durbin(x, nlags=5, isacov=False)
print("PACF (Levinson-Durbin):", np.round(pacf_ld[1:6], 3))
print("innovation variance:", np.round(sigma[1:6], 3))

Output:

PACF (Levinson-Durbin): [ 0.532 -0.311 -0.02  -0.02   0.078]
innovation variance: [1.067 0.964 0.963 0.963 0.957]

The values differ from pacf(method="ywadjusted") only in the third decimal, because the underlying ACF estimator differs (with or without the small-sample bias adjustment). As equation \((10)\) predicts, the innovation variance drops sharply from order 1 to order 2 (\(1.067 \to 0.964\) ) and then flattens out from order 3 onward (\(0.963, 0.963, 0.957\) ). The residual-variance side of the story says the same thing: raising the order beyond 2 does not improve prediction.

Order selection by AIC/BIC (AutoReg / ar_select_order)

We fit AR models of orders \(p = 0, \ldots, 6\) to the same sample (using hold_back=10 to exclude the first 10 points and equalize the estimation sample) and compare AIC/BIC.

from statsmodels.tsa.ar_model import AutoReg, ar_select_order

print(f"{'p':>2} {'AIC':>9} {'BIC':>9}")
for p in range(0, 7):
    res = AutoReg(x, lags=p, hold_back=10, old_names=False).fit()
    print(f"{p:>2} {res.aic:>9.2f} {res.bic:>9.2f}")

sel_aic = ar_select_order(x, maxlag=10, ic="aic", old_names=False)
sel_bic = ar_select_order(x, maxlag=10, ic="bic", old_names=False)
print("lags selected by AIC:", sel_aic.ar_lags)
print("lags selected by BIC:", sel_bic.ar_lags)

res2 = AutoReg(x, lags=2, old_names=False).fit()
print("AutoReg(2) params:", np.round(res2.params, 3))

Output:

 p       AIC       BIC
 0   1592.91   1601.30
 1   1430.99   1443.57
 2   1381.90   1398.68
 3   1383.54   1404.51
 4   1385.38   1410.55
 5   1384.62   1413.98
 6   1385.15   1418.71
lags selected by AIC: [1, 2]
lags selected by BIC: [1, 2]
AutoReg(2) params: [ 0.006  0.702 -0.313]

Both AIC and BIC attain their minima at \(p = 2\) (AIC \(= 1381.90\) , BIC \(= 1398.68\) ), and ar_select_order picks lags \(\{1, 2\}\) under either criterion. The PACF cutoff (lag 2) and the information-criterion choice (\(p = 2\) ) agree. The estimated coefficients \((\hat{\phi}_1, \hat{\phi}_2) = (0.702, -0.313)\) also recover the true values \((0.7, -0.3)\) accurately.

Searching over full ARMA orders (arma_order_select_ic / ARIMA)

Does AR(2) still win when MA terms are allowed as candidates? We check with arma_order_select_ic and fit the final model by maximum likelihood with ARIMA.

from statsmodels.tsa.stattools import arma_order_select_ic
from statsmodels.tsa.arima.model import ARIMA

order_sel = arma_order_select_ic(x, max_ar=4, max_ma=2, ic=["aic", "bic"])
print("AIC-minimizing order (p, q):", tuple(map(int, order_sel.aic_min_order)))
print("BIC-minimizing order (p, q):", tuple(map(int, order_sel.bic_min_order)))

arima_res = ARIMA(x, order=(2, 0, 0)).fit()
print("ARIMA(2,0,0) AR coefficients:", np.round(arima_res.arparams, 3))
print("ARIMA(2,0,0) residual variance:", round(arima_res.params[-1], 3))

Output:

AIC-minimizing order (p, q): (2, 0)
BIC-minimizing order (p, q): (2, 0)
ARIMA(2,0,0) AR coefficients: [ 0.7   -0.312]
ARIMA(2,0,0) residual variance: 0.96

Even searching the 2-D grid of \((p, q)\) (\(5 \times 3 = 15\) models), both AIC and BIC select \((2, 0)\) . The maximum-likelihood AR coefficients \((0.700, -0.312)\) and residual variance \(0.96\) (true value \(1.0\) ) are consistent as well. Three independent lines of evidence—PACF, AIC, and BIC—all point to AR(2), which is the ideal outcome of order identification.

Contrast experiments: an MA(2) process and white noise

Let us verify the remaining rows of the decision table numerically. We generate an MA(2) process \(x_t = \varepsilon_t + 0.6\varepsilon_{t-1} + 0.4\varepsilon_{t-2}\) and pure white noise, and run the same analysis.

# MA(2): the ACF should cut off after lag 2
np.random.seed(29)
theta = np.array([0.6, 0.4])
y = ArmaProcess(np.r_[1], np.r_[1, theta]).generate_sample(nsample=N)
acf_y = acf(y, nlags=10, fft=True)
pacf_y = pacf(y, nlags=10, method="ywadjusted")
print("MA(2) ACF  1-5:", np.round(acf_y[1:6], 3))
print("MA(2) PACF 1-5:", np.round(pacf_y[1:6], 3))
print("MA(2) significant ACF lags:",
      (np.where(np.abs(acf_y[1:]) > ci)[0] + 1).tolist())

# white noise: neither should be significant at any lag
np.random.seed(2)
w = np.random.randn(N)
acf_w = acf(w, nlags=10, fft=True)
pacf_w = pacf(w, nlags=10, method="ywadjusted")
print("WN ACF  1-5:", np.round(acf_w[1:6], 3))
print("WN PACF 1-5:", np.round(pacf_w[1:6], 3))
print("WN significant lags (ACF/PACF):",
      (np.where(np.abs(acf_w[1:]) > ci)[0] + 1).tolist(),
      (np.where(np.abs(pacf_w[1:]) > ci)[0] + 1).tolist())

Output:

MA(2) ACF  1-5: [ 0.598  0.31   0.065  0.003 -0.01 ]
MA(2) PACF 1-5: [ 0.599 -0.074 -0.142  0.065  0.009]
MA(2) significant ACF lags: [1, 2]
WN ACF  1-5: [-0.034  0.002  0.039 -0.027  0.039]
WN PACF 1-5: [-0.034  0.001  0.04  -0.025  0.038]
WN significant lags (ACF/PACF): [] []

A perfect mirror image of the AR(2) case. For the MA(2), the ACF cuts off after lag 2 (significant lags are exactly \(\{1, 2\}\) ; \(\hat{\rho}_3 = 0.065\) onward stays inside the band), while the PACF decays without cutting off: \(0.599 \to -0.074 \to -0.142\) . The sample values \(0.598\) and \(0.310\) are also reasonable against the theoretical \(\rho_1 = 0.553\) and \(\rho_2 = 0.263\) . For white noise, both ACF and PACF stay inside the confidence band at every lag—there is no structure to model.

Practical Caveats

• Stationarity is a prerequisite: both the PACF and Yule-Walker assume a stationary process. For series with trends or unit roots, first difference the data within the ARIMA framework (choose \(d\) ), then identify \(p\) and \(q\)
• The cutoff is a statistical cutoff: about 5% of lags will be significant by chance. Indeed, in our AR(2) experiment, extending the plot to lag 20 shows lag 19 spuriously crossing the band. Weight consecutive significance starting from low lags, and ignore isolated spikes at high lags
• Mixed ARMA cannot be settled by the table alone: when both ACF and PACF tail off, use the plots only to narrow candidates and delegate the decision to an information-criterion grid search such as arma_order_select_ic
• Do not skip residual diagnostics after identification: confirming that the chosen model’s residuals are white noise (residual ACF insignificant at all lags, Ljung-Box test not rejected) completes the workflow. Detecting leftover structure in residuals is the same idea as the residual analysis in time-series anomaly detection
• Limits of linear models: the PACF captures linear dependence only. When nonlinear dependence dominates, a nonlinear model such as an LSTM becomes necessary. Even then, ACF/PACF remain useful for designing lag features and choosing input window lengths, sitting at the entrance of the pipeline organized in the machine learning for time series guide

Summary

• The partial autocorrelation (PACF) is the correlation after removing the linear effect of intermediate lags, and equals the last coefficient \(\phi_{kk}\) of an AR(k) regression
• Solving the Yule-Walker equations (a Toeplitz system) at increasing orders yields the PACF; the Levinson-Durbin recursion delivers coefficients, PACF, and innovation variances for all orders at once in \(O(p^2)\)
• Decision table for order identification: the PACF of AR(p) cuts off after lag \(p\) ; the ACF of MA(q) cuts off after lag \(q\) ; both tail off for ARMA
• In the AR(2) experiment, the PACF cutoff (\(\hat{\alpha}(2) = -0.312\) significant, \(\hat{\alpha}(3)\) onward within \(\pm 0.088\) ), ar_select_order under AIC/BIC (both \(p = 2\) ), and arma_order_select_ic (both criteria at \((2, 0)\) ) all agreed
• For MA(2) the ACF cut off instead (significant lags \(\{1, 2\}\) only), and for white noise neither function was significant at any lag—every row of the decision table was verified numerically
• The decision table and information criteria are complements, not competitors: narrow candidates visually, decide objectively with AIC/BIC, and close with a whiteness test of the residuals

Related Articles

• Autocorrelation and Cross-Correlation: Theory and Python Implementation - The basics companion to this article, covering the ACF definition, fast FFT-based computation, and periodicity detection
• ARIMA and SARIMA Theory and Python Implementation - The practical follow-up that uses the orders \((p, d, q)\) identified here for actual forecasting
• LSTM for Time Series Forecasting - A deep-learning approach for the nonlinear dependence the PACF cannot capture
• Time-Series Anomaly Detection - Residual-based anomaly detection methods sharing the whiteness-checking mindset of model diagnostics
• Machine Learning Hub for Time Series - A bird’s-eye view from classical ACF/PACF identification to state-space models and deep learning

References

• Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
• Brockwell, P. J., & Davis, R. A. (2016). Introduction to Time Series and Forecasting (3rd ed.). Springer.
• Durbin, J. (1960). “The fitting of time-series models.” Revue de l’Institut International de Statistique, 28(3), 233-244.
• Akaike, H. (1974). “A new look at the statistical model identification.” IEEE Transactions on Automatic Control, 19(6), 716-723.
• statsmodels: statsmodels.tsa.stattools.pacf
• statsmodels: statsmodels.tsa.ar_model.AutoReg