'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-05):

    m = self.moment_sequence()
    (mu_x, mu_y, Sigma_x, Sigma_y) = next(m)
    i = 0
    error = (tol + 1)
    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)
    (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)
