def _compute_multivariate_acovf_from_coefficients(
        coefficients, error_variance, maxlag=None,
        forward_autocovariances=False):
    """
    Compute multivariate autocovariances from vector autoregression coefficient
    matrices

    Parameters
    ----------
    coefficients : array or list
        The coefficients matrices. If a list, should be a list of length
        `order`, where each element is an array sized `k_endog` x `k_endog`. If
        an array, should be the coefficient matrices horizontally concatenated
        and sized `k_endog` x `k_endog * order`.
    error_variance : array
        The variance / covariance matrix of the error term. Should be sized
        `k_endog` x `k_endog`.
    maxlag : integer, optional
        The maximum autocovariance to compute. Default is `order`-1. Can be
        zero, in which case it returns the variance.
    forward_autocovariances : boolean, optional
        Whether or not to compute forward autocovariances
        :math:`E(y_t y_{t+j}')`. Default is False, so that backward
        autocovariances :math:`E(y_t y_{t-j}')` are returned.

    Returns
    -------
    autocovariances : list
        A list of the first `maxlag` autocovariance matrices. Each matrix is
        shaped `k_endog` x `k_endog`.

    Notes
    -----
    Computes

    ..math::

        \Gamma(j) = E(y_t y_{t-j}')

    for j = 1, ..., `maxlag`, unless `forward_autocovariances` is specified,
    in which case it computes:

    ..math::

        E(y_t y_{t+j}') = \Gamma(j)'

    Coefficients are assumed to be provided from the VAR model:

    .. math::
        y_t = A_1 y_{t-1} + \dots + A_p y_{t-p} + \varepsilon_t

    Autocovariances are calculated by solving the associated discrete Lyapunov
    equation of the state space representation of the VAR process.

    """
    from scipy import linalg

    # Convert coefficients to a list of matrices, for use in
    # `companion_matrix`; get dimensions
    if type(coefficients) == list:
        order = len(coefficients)
        k_endog = coefficients[0].shape[0]
    else:
        k_endog, order = coefficients.shape
        order //= k_endog

        coefficients = [
            coefficients[:k_endog, i*k_endog:(i+1)*k_endog]
            for i in range(order)
        ]

    if maxlag is None:
        maxlag = order-1

    # Start with VAR(p): w_{t+1} = phi_1 w_t + ... + phi_p w_{t-p+1} + u_{t+1}
    # Then stack the VAR(p) into a VAR(1) in companion matrix form:
    # z_{t+1} = F z_t + v_t
    companion = companion_matrix(
        [1] + [-coefficients[i] for i in range(order)]
    ).T

    # Compute the error variance matrix for the stacked form: E v_t v_t'
    selected_variance = np.zeros(companion.shape)
    selected_variance[:k_endog, :k_endog] = error_variance

    # Compute the unconditional variance of z_t: E z_t z_t'
    stacked_cov = linalg.solve_discrete_lyapunov(companion, selected_variance)

    # The first (block) row of the variance of z_t gives the first p-1
    # autocovariances of w_t: \Gamma_i = E w_t w_t+i with \Gamma_0 = Var(w_t)
    # Note: these are okay, checked against ArmaProcess
    autocovariances = [
        stacked_cov[:k_endog, i*k_endog:(i+1)*k_endog]
        for i in range(min(order, maxlag+1))
    ]

    for i in range(maxlag - (order-1)):
        stacked_cov = np.dot(companion, stacked_cov)
        autocovariances += [
            stacked_cov[:k_endog, -k_endog:]
        ]

    if forward_autocovariances:
        for i in range(len(autocovariances)):
            autocovariances[i] = autocovariances[i].T

    return autocovariances
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/tsa/statespace/tools.py/_compute_multivariate_acovf_from_coefficients/original