def _compute_multivariate_pacf_from_autocovariances(autocovariances,
                                                    order=None, k_endog=None):
    """
    Compute multivariate partial autocorrelations from autocovariances.

    Parameters
    ----------
    autocovariances : list
        Autocorrelations matrices. Should be a list of length `order` + 1,
        where each element is an array sized `k_endog` x `k_endog`.
    order : integer, optional
        The order of the autoregression.
    k_endog : integer, optional
        The dimension of the data vector.

    Returns
    -------
    pacf : list
        List of first `order` multivariate partial autocorrelations.

    Notes
    -----
    Note that this computes multivariate partial autocorrelations.

    Corresponds to the inverse of Lemma 2.1 in Ansley and Kohn (1986). See
    `unconstrain_stationary_multivariate` for more details.

    Notes
    -----
    Computes sample partial autocorrelations if sample autocovariances are
    given.

    """
    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
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/tsa/statespace/tools.py/_compute_multivariate_pacf_from_autocovariances/operator `-` -> `+`