code: def arma_acovf(ar, ma, nobs=10, sigma2=1, dtype=None):
      """
      Theoretical autocovariances of stationary ARMA processes
  
      Parameters
      ----------
      ar : array_like, 1d
          The coefficients for autoregressive lag polynomial, including zero lag.
      ma : array_like, 1d
          The coefficients for moving-average lag polynomial, including zero lag.
      nobs : int
          The number of terms (lags plus zero lag) to include in returned acovf.
      sigma2 : float
          Variance of the innovation term.
  
      Returns
      -------
      ndarray
          The autocovariance of ARMA process given by ar, ma.
  
      See Also
      --------
      arma_acf : Autocorrelation function for ARMA processes.
      acovf : Sample autocovariance estimation.
  
      References
      ----------
      .. [*] Brockwell, Peter J., and Richard A. Davis. 2009. Time Series:
          Theory and Methods. 2nd ed. 1991. New York, NY: Springer.
      """
      if dtype is None:
          dtype = np.common_type(
              np.array(ar), np.array(ma), np.array(sigma2)
          )
  
      p = len(ar) - 1
      q = len(ma) - 1
      m = max(p, q) + 1
  
      if sigma2.real < 0:
          raise ValueError("Must have positive innovation variance.")
  
      # Short-circuit for trivial corner-case
      if p == q == 0:
          out = np.zeros(nobs, dtype=dtype)
          out[0] = sigma2
          return out
      elif p > 0 and np.max(np.abs(np.roots(ar))) >= 1:
          raise ValueError(NONSTATIONARY_ERROR)
  
      # Get the moving average representation coefficients that we need
      ma_coeffs = arma2ma(ar, ma, lags=m)
  
      # Solve for the first m autocovariances via the linear system
      # described by (BD, eq. 3.3.8)
      A = np.zeros((m, m), dtype=dtype)
      b = np.zeros((m, 1), dtype=dtype)
      # We need a zero-right-padded version of ar params
      tmp_ar = np.zeros(m, dtype=dtype)
      tmp_ar[: p + 1] = ar
      for k in range(m):
          A[k, : (k + 1)] = tmp_ar[: (k + 1)][::-1]
          A[k, 1 : m - k] += tmp_ar[(k + 1) : m]
          b[k] = sigma2 * np.dot(
              ma[k : q + 1], ma_coeffs[: max((q + 1 - k), 0)]
          )
      acovf = np.zeros(max(nobs, m), dtype=dtype)
      try:
          acovf[:m] = np.linalg.solve(A, b)[:, 0]
      except np.linalg.LinAlgError:
          raise ValueError(NONSTATIONARY_ERROR)
  
      # Iteratively apply (BD, eq. 3.3.9) to solve for remaining autocovariances
      if nobs > m:
          zi = signal.lfiltic([1], ar, acovf[:m:][::-1])
          acovf[m:] = signal.lfilter(
              [1], ar, np.zeros(nobs - m, dtype=dtype), zi=zi
          )[0]
  
      return acovf[:nobs]
  
imports: from statsmodels.compat.numpy import NP_LT_2
  from statsmodels.compat.pandas import Appender
  import warnings
  import numpy
  from scipy import linalg, optimize, signal
  from statsmodels.tools.docstring import Docstring, remove_parameters
  from statsmodels.tools.validation import array_like
label_desc: Use image node_id 2 for calling a global function with example usage: arma_acovf(ar, ma, nobs, sigma2, dtype) and returns: acovf, out