docs: def moment_sequence(self):
          r"""
          Create a generator to calculate the population mean and
          variance-convariance matrix for both :math:`x_t` and :math:`y_t`
          starting at the initial condition (self.mu_0, self.Sigma_0).
          Each iteration produces a 4-tuple of items (mu_x, mu_y, Sigma_x,
          Sigma_y) for the next period.
  
          Yields
          ------
          mu_x : array_like(float)
              An n x 1 array representing the population mean of x_t
          mu_y : array_like(float)
              A  k x 1 array representing the population mean of y_t
          Sigma_x : array_like(float)
              An n x n array representing the variance-covariance matrix
              of x_t
          Sigma_y : array_like(float)
              A k x k array representing the variance-covariance matrix
              of y_t
  
          """
          # == Simplify names == #
          A, C, G, H = self.A, self.C, self.G, self.H
          # == Initial moments == #
          mu_x, Sigma_x = self.mu_0, self.Sigma_0
  
          while 1:
              mu_y = G.dot(mu_x)
              if H is None:
                  Sigma_y = G.dot(Sigma_x).dot(G.T)
              else:
                  Sigma_y = G.dot(Sigma_x).dot(G.T) + H.dot(H.T)
  
              yield mu_x, mu_y, Sigma_x, Sigma_y
  
              # == Update moments of x == #
              mu_x = A.dot(mu_x)
              Sigma_x = A.dot(Sigma_x).dot(A.T) + C.dot(C.T)