def _compute_multivariate_pacf_from_autocovariances(autocovariances,
                                                    order=None, k_endog=None):
    from scipy import linalg

    if order is None:
        order = len(autocovariances)-1
    if k_endog is None:
        k_endog = autocovariances[0].shape[0]

    # Now apply the Ansley and Kohn (1986) algorithm, except that instead of
    # calculating phi_{s+1, s+1} = L_s P_{s+1} {L_s^*}^{-1} (which requires
    # the partial autocorrelation P_{s+1} which is what we're trying to
    # calculate here), we calculate it as in Ansley and Newbold (1979), using
    # the autocovariances \Gamma_s and the forwards and backwards residual
    # variances \Sigma_s, \Sigma_s^*:
    # phi_{s+1, s+1} = [ \Gamma_{s+1}' - \phi_{s,1} \Gamma_s' - ... -
    #                    \phi_{s,s} \Gamma_1' ] {\Sigma_s^*}^{-1}

    # Forward and backward variances
    forward_variances = []   # \Sigma_s
    backward_variances = []  # \Sigma_s^*,  s = 0, ..., p
    # \phi_{s,k}, s = 1, ..., p
    #             k = 1, ..., s+1
    forwards = []
    # \phi_{s,k}^*
    backwards = []

    forward_factors = []   # L_s
    backward_factors = []  # L_s^*,  s = 0, ..., p

    # Ultimately we want to construct the partial autocorrelation matrices
    # Note that this is "1-indexed" in the sense that it stores P_1, ... P_p
    # rather than starting with P_0.
    partial_autocorrelations = []

    # We fill in the entries of phi_{s,k} as follows:
    # [1,1]
    # [2,2], [2,1]
    # [3,3], [3,1], [3,2]
    # ...
    # [p,p], [p,1], ..., [p,p-1]
    # the last row, correctly ordered, should be the same as the coefficient
    # matrices provided in the argument `constrained`
    for s in range(order):  # s = 0, ..., p-1
        prev_forwards = list(forwards)
        prev_backwards = list(backwards)
        forwards = []
        backwards = []

        # Create forward and backwards variances Sigma_s, Sigma*_s
        forward_variance = autocovariances[0].copy()
        backward_variance = autocovariances[0].T.copy()

        for k in range(s):
            forward_variance -= np.dot(prev_forwards[k],
                                       autocovariances[k+1])
            backward_variance -= np.dot(prev_backwards[k],
                                        autocovariances[k+1].T)

        forward_variances.append(forward_variance)
        backward_variances.append(backward_variance)

        # Cholesky factors
        forward_factors.append(
            linalg.cholesky(forward_variances[s], lower=True)
        )
        backward_factors.append(
            linalg.cholesky(backward_variances[s], lower=True)
        )

        # Create the intermediate sum term
        if s == 0:
            # phi_11 = \Gamma_1' \Gamma_0^{-1}
            # phi_11 \Gamma_0 = \Gamma_1'
            # \Gamma_0 phi_11' = \Gamma_1
            forwards.append(linalg.cho_solve(
                (forward_factors[0], True), autocovariances[1]).T)
            # backwards.append(forwards[-1])
            # phi_11_star = \Gamma_1 \Gamma_0^{-1}
            # phi_11_star \Gamma_0 = \Gamma_1
            # \Gamma_0 phi_11_star' = \Gamma_1'
            backwards.append(linalg.cho_solve(
                (backward_factors[0], True), autocovariances[1].T).T)
        else:
            # G := \Gamma_{s+1}' -
            #      \phi_{s,1} \Gamma_s' - .. - \phi_{s,s} \Gamma_1'
            tmp_sum = autocovariances[s+1].T.copy()

            for k in range(s):
                tmp_sum -= np.dot(prev_forwards[k], autocovariances[s-k].T)

            # Create the "last" (k = s+1) matrix
            # Note: this is for k = s+1. However, below we then have to
            # fill in for k = 1, ..., s in order.
            # phi = G Sigma*^{-1}
            # phi Sigma* = G
            # Sigma*' phi' = G'
            # Sigma* phi' = G'
            # (because Sigma* is symmetric)
            forwards.append(linalg.cho_solve(
                (backward_factors[s], True), tmp_sum.T).T)

            # phi = G' Sigma^{-1}
            # phi Sigma = G'
            # Sigma' phi' = G
            # Sigma phi' = G
            # (because Sigma is symmetric)
            backwards.append(linalg.cho_solve(
                (forward_factors[s], True), tmp_sum).T)

        # Create the remaining k = 1, ..., s matrices,
        # only has an effect if s >= 1
        for k in range(s):
            forwards.insert(k, prev_forwards[k] - np.dot(
                forwards[-1], prev_backwards[s-(k+1)]))
            backwards.insert(k, prev_backwards[k] - np.dot(
                backwards[-1], prev_forwards[s-(k+1)]))

        # Partial autocorrelation matrix: P_{s+1}
        # P = L^{-1} phi L*
        # L P = (phi L*)
        partial_autocorrelations.append(linalg.solve_triangular(
            forward_factors[s], np.dot(forwards[s], backward_factors[s]),
            lower=True))

    return partial_autocorrelations
Compute multivariate partial autocorrelations from autocovariances.
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/tsa/statespace/tools.py/_compute_multivariate_pacf_from_autocovariances