docs: def stationary_distributions(self, max_iter=200, tol=1e-5):
          r"""
          Compute the moments of the stationary distributions of :math:`x_t` and
          :math:`y_t` if possible.  Computation is by iteration, starting from
          the initial conditions self.mu_0 and self.Sigma_0
  
          Parameters
          ----------
          max_iter : scalar(int), optional(default=200)
              The maximum number of iterations allowed
          tol : scalar(float), optional(default=1e-5)
              The tolerance level that one wishes to achieve
  
          Returns
          -------
          mu_x_star : array_like(float)
              An n x 1 array representing the stationary mean of :math:`x_t`
          mu_y_star : array_like(float)
              An k x 1 array representing the stationary mean of :math:`y_t`
          Sigma_x_star : array_like(float)
              An n x n array representing the stationary var-cov matrix
              of :math:`x_t`
          Sigma_y_star : array_like(float)
              An k x k array representing the stationary var-cov matrix
              of :math:`y_t`
  
          """
          # == Initialize iteration == #
          m = self.moment_sequence()
          mu_x, mu_y, Sigma_x, Sigma_y = next(m)
          i = 0
          error = tol + 1
  
          # == Loop until convergence or failure == #
          while error > tol:
  
              if i > max_iter:
                  fail_message = 'Convergence failed after {} iterations'
                  raise ValueError(fail_message.format(max_iter))
  
              else:
                  i += 1
                  mu_x1, mu_y1, Sigma_x1, Sigma_y1 = next(m)
                  error_mu = np.max(np.abs(mu_x1 - mu_x))
                  error_Sigma = np.max(np.abs(Sigma_x1 - Sigma_x))
                  error = max(error_mu, error_Sigma)
                  mu_x, Sigma_x = mu_x1, Sigma_x1
  
          # == Prepare return values == #
          mu_x_star, Sigma_x_star = mu_x, Sigma_x
          mu_y_star, Sigma_y_star = mu_y1, Sigma_y1
  
          return mu_x_star, mu_y_star, Sigma_x_star, Sigma_y_star