code: def ARcovariance(rho, n, cor=False, sigma=1.0):
      """
      Return covariance matrix of a sample of length n from an AR(p)
      process with parameters rho.
  
      INPUTS:
  
      rho      -- an array of length p
      sigma    -- standard deviation of the white noise
      """
  
      rho = np.asarray(rho)
      p = rho.shape[0]
      invK = np.identity(n)
      for i in range(p):
          invK -= np.diag(
              (rho[i] / sigma) * np.ones(n - i - 1), k=-i - 1
          )
      K = inv(invK)
      Q = np.dot(K, K.T)
      if cor:
          sd = np.sqrt(np.diag(Q))
          sdout = np.multiply.outer(sd, sd)
          Q /= sd
      return Q
  
imports: import numpy
  from numpy.linalg import inv
label_desc: Use image node_id 1 for calling a global function with example usage: ARcovariance(rho, n, cor, sigma) and returns: Q