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`

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