code: def _construct_var_eqns(data, p, delta=None):
          
          t, m, l = np.shape(data)
          n = (l - p) * t  # number of linear relations
          rows = n if delta is None else n + m * p
  
          # Construct matrix x (predictor variables)
          x = np.zeros((rows, m * p))
          for i in range(m):
              for k in range(1, p + 1):
                  x[:n, i * p + k - 1] = np.reshape(data[:, i, p - k:-k].T, n)
          if delta is not None:
              np.fill_diagonal(x[n:, :], delta)
  
          # Construct vectors yi (response variables for each channel i)
          y = np.zeros((rows, m))
          for i in range(m):
              y[:n, i] = np.reshape(data[:, i, p:].T, n)
  
          return x, y
docs: Construct VAR equation system (optionally with RLS constraint).