code: def _compute_coefficients_from_multivariate_pacf_python(
      partial_autocorrelations,
      error_variance,
      transform_variance=False,
      order=None,
      k_endog=None,
  ):
      """
      Transform matrices with singular values less than one to matrices
      corresponding to a stationary (or invertible) process.
  
      Parameters
      ----------
      partial_autocorrelations : list
          Partial autocorrelation matrices. Should be a list of length `order`,
          where each element is an array sized `k_endog` x `k_endog`.
      error_variance : ndarray
          The variance / covariance matrix of the error term. Should be sized
          `k_endog` x `k_endog`. This is used as input in the algorithm even if
          is not transformed by it (when `transform_variance` is False). The
          error term variance is required input when transformation is used
          either to force an autoregressive component to be stationary or to
          force a moving average component to be invertible.
      transform_variance : bool, optional
          Whether or not to transform the error variance term. This option is
          not typically used, and the default is False.
      order : int, optional
          The order of the autoregression.
      k_endog : int, optional
          The dimension of the data vector.
  
      Returns
      -------
      coefficient_matrices : list
          Transformed coefficient matrices leading to a stationary VAR
          representation.
  
      See Also
      --------
      constrain_stationary_multivariate
  
      Notes
      -----
      Corresponds to Lemma 2.1 in Ansley and Kohn (1986). See
      `constrain_stationary_multivariate` for more details.
      """
      from scipy import linalg
  
      if order is None:
          order = len(partial_autocorrelations)
      if k_endog is None:
          k_endog = partial_autocorrelations[0].shape[0]
  
      # If we want to keep the provided variance but with the constrained
      # coefficient matrices, we need to make a copy here, and then after the
      # main loop we will transform the coefficients to match the passed variance
      if not transform_variance:
          initial_variance = error_variance
          # Need to make the input variance large enough that the recursions
          # do not lead to zero-matrices due to roundoff error, which would case
          # exceptions from the Cholesky decompositions.
          # Note that this will still not always ensure positive definiteness,
          # and for k_endog, order large enough an exception may still be raised
          error_variance = np.eye(k_endog) * (order + k_endog) ** 10
  
      forward_variances = [error_variance]  # \Sigma_s
      backward_variances = [
          error_variance
      ]  # \Sigma_s^*,  s = 0, ..., p
      autocovariances = [error_variance]  # \Gamma_s
      # \phi_{s,k}, s = 1, ..., p
      #             k = 1, ..., s+1
      forwards = []
      # \phi_{s,k}^*
      backwards = []
  
      error_variance_factor = linalg.cholesky(
          error_variance, lower=True
      )
  
      forward_factors = [error_variance_factor]
      backward_factors = [error_variance_factor]
  
      # We fill in the entries 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, is then used as the coefficients
      for s in range(order):  # s = 0, ..., p-1
          prev_forwards = forwards
          prev_backwards = backwards
          forwards = []
          backwards = []
  
          # 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.
          # P L*^{-1} = x
          # x L* = P
          # L*' x' = P'
          forwards.append(
              linalg.solve_triangular(
                  backward_factors[s],
                  partial_autocorrelations[s].T,
                  lower=True,
                  trans="T",
              )
          )
          forwards[0] = np.dot(forward_factors[s], forwards[0].T)
  
          # P' L^{-1} = x
          # x L = P'
          # L' x' = P
          backwards.append(
              linalg.solve_triangular(
                  forward_factors[s],
                  partial_autocorrelations[s],
                  lower=True,
                  trans="T",
              )
          )
          backwards[0] = np.dot(backward_factors[s], backwards[0].T)
  
          # Update the variance
          # Note: if s >= 1, this will be further updated in the for loop
          # below
          # Also, this calculation will be re-used in the forward variance
          tmp = np.dot(forwards[0], backward_variances[s])
          autocovariances.append(tmp.copy().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)]),
              )
  
              autocovariances[s + 1] += np.dot(
                  autocovariances[k + 1], prev_forwards[s - (k + 1)].T
              )
  
          # Create forward and backwards variances
          forward_variances.append(
              forward_variances[s] - np.dot(tmp, forwards[s].T)
          )
          backward_variances.append(
              backward_variances[s]
              - np.dot(
                  np.dot(backwards[s], forward_variances[s]),
                  backwards[s].T,
              )
          )
  
          # Cholesky factors
          forward_factors.append(
              linalg.cholesky(forward_variances[s + 1], lower=True)
          )
          backward_factors.append(
              linalg.cholesky(backward_variances[s + 1], lower=True)
          )
  
      # If we do not want to use the transformed variance, we need to
      # adjust the constrained matrices, as presented in Lemma 2.3, see above
      variance = forward_variances[-1]
      if not transform_variance:
          # Here, we need to construct T such that:
          # variance = T * initial_variance * T'
          # To do that, consider the Cholesky of variance (L) and
          # input_variance (M) to get:
          # L L' = T M M' T' = (TM) (TM)'
          # => L = T M
          # => L M^{-1} = T
          initial_variance_factor = np.linalg.cholesky(initial_variance)
          transformed_variance_factor = np.linalg.cholesky(variance)
          transform = np.dot(
              initial_variance_factor,
              np.linalg.inv(transformed_variance_factor),
          )
          inv_transform = np.linalg.inv(transform)
  
          for i in range(order):
              forwards[i] = np.dot(
                  np.dot(transform, forwards[i]), inv_transform
              )
  
      return forwards, variance
  
imports: import numpy
  from scipy.linalg import solve_sylvester
  import pandas
  from statsmodels.compat.pandas import Appender
  from statsmodels.tools.data import _is_using_pandas
  from scipy.linalg.blas import find_best_blas_type
  from .None import _initialization, _representation, _kalman_filter, _kalman_smoother, _simulation_smoother, _cfa_simulation_smoother, _tools
label_desc: Use image node_id 10 for calling a global function with example usage: _compute_coefficients_from_multivariate_pacf_python(partial_autocorrelations, error_variance, transform_variance, order, k_endog) and returns: forwards,  variance