diff --git a/quantecon/__init__.py b/quantecon/__init__.py
index d539f4151..4dc3b901b 100644
--- a/quantecon/__init__.py
+++ b/quantecon/__init__.py
@@ -7,13 +7,15 @@
 from .discrete_rv import DiscreteRV
 from .ecdf import ECDF
 from .estspec import smooth, periodogram, ar_periodogram
+from .graph_tools import DiGraph
+from .gth_solve import gth_solve
 from .kalman import Kalman
 from .lae import LAE
 from .arma import ARMA
 from .lqcontrol import LQ
 from .lss import LSS
 from .matrix_eqn import solve_discrete_lyapunov, solve_discrete_riccati
-from .mc_tools import DMarkov, mc_sample_path, mc_compute_stationary
+from .mc_tools import MarkovChain, mc_compute_stationary, mc_sample_path
 from .quadsums import var_quadratic_sum, m_quadratic_sum
 from .rank_nullspace import rank_est, nullspace
 from .robustlq import RBLQ
@@ -21,4 +23,4 @@
 from . import quad as quad
 
 #Add Version Attribute
-from .version import version as __version__
\ No newline at end of file
+from .version import version as __version__
diff --git a/quantecon/graph_tools.py b/quantecon/graph_tools.py
new file mode 100644
index 000000000..8cc051318
--- /dev/null
+++ b/quantecon/graph_tools.py
@@ -0,0 +1,275 @@
+"""
+Filename: graph_tools.py
+
+Author: Daisuke Oyama
+
+Tools for dealing with a directed graph.
+
+"""
+import numpy as np
+from scipy import sparse
+from scipy.sparse import csgraph
+from fractions import gcd
+
+
+class DiGraph:
+    r"""
+    Class for a directed graph. It stores useful information about the
+    graph structure such as strong connectivity [1]_ and periodicity
+    [2]_.
+
+    Parameters
+    ----------
+    adj_matrix : array_like(ndim=2)
+        Adjacency matrix representing a directed graph. Must be of shape
+        n x n.
+
+    weighted : bool, optional(default=False)
+        Whether to treat `adj_matrix` as a weighted adjacency matrix.
+
+    Attributes
+    ----------
+    csgraph : scipy.sparse.csr_matrix
+        Compressed sparse representation of the digraph.
+
+    is_strongly_connected : bool
+        Indicate whether the digraph is strongly connected.
+
+    num_strongly_connected_components : int
+        The number of the strongly connected components.
+
+    strongly_connected_components : list(ndarray(int))
+        List of numpy arrays containing the strongly connected
+        components.
+
+    num_sink_strongly_connected_components : int
+        The number of the sink strongly connected components.
+
+    sink_strongly_connected_components : list(ndarray(int))
+        List of numpy arrays containing the sink strongly connected
+        components.
+
+    is_aperiodic : bool
+        Indicate whether the digraph is aperiodic.
+
+    period : int
+        The period of the digraph. Defined only for a strongly connected
+        digraph.
+
+    cyclic_components : list(ndarray(int))
+        List of numpy arrays containing the cyclic components.
+
+    Notes
+    -----
+    For the definitions, see the Wikipedia entries [1]_, [2]_.
+
+    References
+    ----------
+    .. [1] `Strongly connected component
+       <http://en.wikipedia.org/wiki/Strongly_connected_component>`_,
+       Wikipedia.
+
+    .. [2] `Aperiodic graph
+       <http://en.wikipedia.org/wiki/Aperiodic_graph>`_, Wikipedia.
+
+    """
+
+    def __init__(self, adj_matrix, weighted=False):
+        if weighted:
+            dtype = None
+        else:
+            dtype = bool
+        self.csgraph = sparse.csr_matrix(adj_matrix, dtype=dtype)
+
+        m, n = self.csgraph.shape
+        if n != m:
+            raise ValueError('input matrix must be square')
+
+        self.n = n  # Number of nodes
+
+        self._num_scc = None
+        self._scc_proj = None
+        self._sink_scc_labels = None
+
+        self._period = None
+
+    def _find_scc(self):
+        """
+        Set ``self._num_scc`` and ``self._scc_proj``
+        by calling ``scipy.sparse.csgraph.connected_components``:
+        * docs.scipy.org/doc/scipy/reference/sparse.csgraph.html
+        * github.com/scipy/scipy/blob/master/scipy/sparse/csgraph/_traversal.pyx
+
+        ``self._scc_proj`` is a list of length `n` that assigns to each node
+        the label of the strongly connected component to which it belongs.
+
+        """
+        # Find the strongly connected components
+        self._num_scc, self._scc_proj = \
+            csgraph.connected_components(self.csgraph, connection='strong')
+
+    @property
+    def num_strongly_connected_components(self):
+        if self._num_scc is None:
+            self._find_scc()
+        return self._num_scc
+
+    @property
+    def scc_proj(self):
+        if self._scc_proj is None:
+            self._find_scc()
+        return self._scc_proj
+
+    @property
+    def is_strongly_connected(self):
+        return (self.num_strongly_connected_components == 1)
+
+    def _condensation_lil(self):
+        """
+        Return the sparse matrix representation of the condensation digraph
+        in lil format.
+
+        """
+        condensation_lil = sparse.lil_matrix(
+            (self.num_strongly_connected_components,
+             self.num_strongly_connected_components), dtype=bool
+        )
+
+        scc_proj = self.scc_proj
+        for node_from, node_to in _csr_matrix_indices(self.csgraph):
+            scc_from, scc_to = scc_proj[node_from], scc_proj[node_to]
+            if scc_from != scc_to:
+                condensation_lil[scc_from, scc_to] = True
+
+        return condensation_lil
+
+    def _find_sink_scc(self):
+        """
+        Set self._sink_scc_labels, which is a list containing the labels of
+        the strongly connected components.
+
+        """
+        condensation_lil = self._condensation_lil()
+
+        # A sink SCC is a SCC such that none of its members is strongly
+        # connected to nodes in other SCCs
+        # Those k's such that graph_condensed_lil.rows[k] == []
+        self._sink_scc_labels = \
+            np.where(np.logical_not(condensation_lil.rows))[0]
+
+    @property
+    def sink_scc_labels(self):
+        if self._sink_scc_labels is None:
+            self._find_sink_scc()
+        return self._sink_scc_labels
+
+    @property
+    def num_sink_strongly_connected_components(self):
+        return len(self.sink_scc_labels)
+
+    @property
+    def strongly_connected_components(self):
+        if self.is_strongly_connected:
+            return [np.arange(self.n)]
+        else:
+            return [np.where(self.scc_proj == k)[0]
+                    for k in range(self.num_strongly_connected_components)]
+
+    @property
+    def sink_strongly_connected_components(self):
+        if self.is_strongly_connected:
+            return [np.arange(self.n)]
+        else:
+            return [np.where(self.scc_proj == k)[0]
+                    for k in self.sink_scc_labels.tolist()]
+
+    def _compute_period(self):
+        """
+        Set ``self._period`` and ``self._cyclic_components_proj``.
+
+        Use the algorithm described in:
+        J. P. Jarvis and D. R. Shier,
+        "Graph-Theoretic Analysis of Finite Markov Chains," 1996.
+
+        """
+        # Degenerate graph with a single node (which is strongly connected)
+        # csgraph.reconstruct_path would raise an exception
+        # github.com/scipy/scipy/issues/4018
+        if self.n == 1:
+            if self.csgraph[0, 0] == 0:  # No edge: "trivial graph"
+                self._period = 1  # Any universally accepted definition?
+                self._cyclic_components_proj = np.zeros(self.n, dtype=int)
+                return None
+            else:  # Self loop
+                self._period = 1
+                self._cyclic_components_proj = np.zeros(self.n, dtype=int)
+                return None
+
+        if not self.is_strongly_connected:
+            raise NotImplementedError(
+                'Not defined for a non strongly-connected digraph'
+            )
+
+        if np.any(self.csgraph.diagonal() > 0):
+            self._period = 1
+            self._cyclic_components_proj = np.zeros(self.n, dtype=int)
+            return None
+
+        # Construct a breadth-first search tree rooted at 0
+        node_order, predecessors = \
+            csgraph.breadth_first_order(self.csgraph, i_start=0)
+        bfs_tree_csr = \
+            csgraph.reconstruct_path(self.csgraph, predecessors)
+
+        # Edges not belonging to tree_csr
+        non_bfs_tree_csr = self.csgraph - bfs_tree_csr
+        non_bfs_tree_csr.eliminate_zeros()
+
+        # Distance to 0
+        level = np.zeros(self.n, dtype=int)
+        for i in range(1, self.n):
+            level[node_order[i]] = level[predecessors[node_order[i]]] + 1
+
+        # Determine the period
+        d = 0
+        for node_from, node_to in _csr_matrix_indices(non_bfs_tree_csr):
+            value = level[node_from] - level[node_to] + 1
+            d = gcd(d, value)
+            if d == 1:
+                self._period = 1
+                self._cyclic_components_proj = np.zeros(self.n, dtype=int)
+                return None
+
+        self._period = d
+        self._cyclic_components_proj = level % d
+
+    @property
+    def period(self):
+        if self._period is None:
+            self._compute_period()
+        return self._period
+
+    @property
+    def is_aperiodic(self):
+        return (self.period == 1)
+
+    @property
+    def cyclic_components(self):
+        if self.is_aperiodic:
+            return [np.arange(self.n)]
+        else:
+            return [np.where(self._cyclic_components_proj == k)[0]
+                    for k in range(self.period)]
+
+
+def _csr_matrix_indices(S):
+    """
+    Generate the indices of nonzero entries of a csr_matrix S
+
+    """
+    m, n = S.shape
+
+    for i in range(m):
+        for j in range(S.indptr[i], S.indptr[i+1]):
+            row_index, col_index = i, S.indices[j]
+            yield row_index, col_index
diff --git a/quantecon/gth_solve.py b/quantecon/gth_solve.py
new file mode 100644
index 000000000..2bd2b5246
--- /dev/null
+++ b/quantecon/gth_solve.py
@@ -0,0 +1,93 @@
+"""
+Filename: gth_solve.py
+
+Author: Daisuke Oyama
+
+Routine to compute the stationary distribution of an irreducible Markov
+chain by the Grassmann-Taksar-Heyman (GTH) algorithm.
+
+"""
+import numpy as np
+
+try:
+    xrange
+except:  # python3
+    xrange = range
+
+
+def gth_solve(A, overwrite=False):
+    r"""
+    This routine computes the stationary distribution of an irreducible
+    Markov transition matrix (stochastic matrix) or transition rate
+    matrix (generator matrix) `A`.
+
+    More generally, given a Metzler matrix (square matrix whose
+    off-diagonal entries are all nonnegative) `A`, this routine solves
+    for a nonzero solution `x` to `x (A - D) = 0`, where `D` is the
+    diagonal matrix for which the rows of `A - D` sum to zero (i.e.,
+    :math:`D_{ii} = \sum_j A_{ij}` for all :math:`i`). One (and only
+    one, up to normalization) nonzero solution exists corresponding to
+    each reccurent class of `A`, and in particular, if `A` is
+    irreducible, there is a unique solution; when there are more than
+    one solution, the routine returns the solution that contains in its
+    support the first index `i` such that no path connects `i` to any
+    index larger than `i`. The solution is normalized so that its 1-norm
+    equals one. This routine implements the Grassmann-Taksar-Heyman
+    (GTH) algorithm [1]_, a numerically stable variant of Gaussian
+    elimination, where only the off-diagonal entries of `A` are used as
+    the input data. For a nice exposition of the algorithm, see Stewart
+    [2]_, Chapter 10.
+
+    Parameters
+    ----------
+    A : array_like(float, ndim=2)
+        Stochastic matrix or generator matrix. Must be of shape n x n.
+    overwrite : bool, optional(default=False)
+        Whether to overwrite `A`; may improve performance.
+
+    Returns
+    -------
+    x : numpy.ndarray(float, ndim=1)
+        Stationary distribution of `A`.
+
+    References
+    ----------
+    .. [1] W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative
+       Analysis and Steady State Distributions for Markov Chains,"
+       Operations Research (1985), 1107-1116.
+
+    .. [2] W. J. Stewart, Probability, Markov Chains, Queues, and
+       Simulation, Princeton University Press, 2009.
+
+    """
+    A1 = np.array(A, copy=not overwrite)
+
+    if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
+        raise ValueError('matrix must be square')
+
+    n = A1.shape[0]
+
+    x = np.zeros(n)
+
+    # === Reduction === #
+    for i in xrange(n-1):
+        scale = np.sum(A1[i, i+1:n])
+        if scale <= 0:
+            # There is one (and only one) recurrent class contained in
+            # {0, ..., i};
+            # compute the solution associated with that recurrent class.
+            n = i+1
+            break
+        A1[i+1:n, i] /= scale
+
+        A1[i+1:n, i+1:n] += np.dot(A1[i+1:n, i:i+1], A1[i:i+1, i+1:n])
+
+    # === Backward substitution === #
+    x[n-1] = 1
+    for i in xrange(n-2, -1, -1):
+        x[i] = np.sum((x[j] * A1[j, i] for j in xrange(i+1, n)))
+
+    # === Normalization === #
+    x /= np.sum(x)
+
+    return x
diff --git a/quantecon/mc_tools.py b/quantecon/mc_tools.py
index 9f5edd5ad..6a9641ad9 100644
--- a/quantecon/mc_tools.py
+++ b/quantecon/mc_tools.py
@@ -1,158 +1,274 @@
-"""
+r"""
 Authors: Chase Coleman, Spencer Lyon, Daisuke Oyama, Tom Sargent,
          John Stachurski
 
 Filename: mc_tools.py
 
-This file contains some useful objects for handling a discrete markov
-transition matrix.  It contains code written by several people and
-was ultimately compiled into a single file to take advantage of the
-pros of each.
+This file contains some useful objects for handling a finite-state
+discrete-time Markov chain.
+
+Definitions and Some Known Facts about Markov Chains
+----------------------------------------------------
+
+Let :math:`\{X_t\}` be a Markov chain represented by an :math:`n \times
+n` stochastic matrix :math:`P`. State :math:`i` *has access* to state
+:math:`j`, denoted :math:`i \to j`, if :math:`i = j` or :math:`P^k[i, j]
+> 0` for some :math:`k = 1, 2, \ldots`; :math:`i` and `j` *communicate*,
+denoted :math:`i \leftrightarrow j`, if :math:`i \to j` and :math:`j \to
+i`. The binary relation :math:`\leftrightarrow` is an equivalent
+relation. A *communication class* of the Markov chain :math:`\{X_t\}`,
+or of the stochastic matrix :math:`P`, is an equivalent class of
+:math:`\leftrightarrow`. Equivalently, a communication class is a
+*strongly connected component* (SCC) in the associated *directed graph*
+:math:`\Gamma(P)`, a directed graph with :math:`n` nodes where there is
+an edge from :math:`i` to :math:`j` if and only if :math:`P[i, j] > 0`.
+The Markov chain, or the stochastic matrix, is *irreducible* if it
+admits only one communication class, or equivalently, if
+:math:`\Gamma(P)` is *strongly connected*.
+
+A state :math:`i` is *recurrent* if :math:`i \to j` implies :math:`j \to
+i`; it is *transient* if it is not recurrent. For any :math:`i, j`
+contained in a communication class, :math:`i` is recurrent if and only
+if :math:`j` is recurrent. Therefore, recurrence is a property of a
+communication class. Thus, a communication class is a *recurrent class*
+if it contains a recurrent state. Equivalently, a recurrent class is a
+SCC that corresponds to a sink node in the *condensation* of the
+directed graph :math:`\Gamma(P)`, where the condensation of
+:math:`\Gamma(P)` is a directed graph in which each SCC is replaced with
+a single node and there is an edge from one SCC :math:`C` to another SCC
+:math:`C'` if :math:`C \neq C'` and some node in :math:`C` has access to
+some node in :math:`C'`. A recurrent class is also called a *closed
+communication class*. The condensation is acyclic, so that there exists
+at least one recurrent class.
+
+For example, if the entries of :math:`P` are all strictly positive, then
+the whole state space is a communication class as well as a recurrent
+class. (More generally, if there is only one communication class, then
+it is a recurrent class.) As another example, consider the stochastic
+matrix :math:`P = [[1, 0], [0,5, 0.5]]`. This has two communication
+classes, :math:`\{0\}` and :math:`\{1\}`, and :math:`\{0\}` is the only
+recurrent class.
+
+A *stationary distribution* of the Markov chain :math:`\{X_t\}`, or of
+the stochastic matrix :math:`P`, is a nonnegative vector :math:`x` such
+that :math:`x' P = x'` and :math:`x' \mathbf{1} = 1`, where
+:math:`\mathbf{1}` is the vector of ones. The Markov chain has a unique
+stationary distribution if and only if it has a unique recurrent class.
+More generally, each recurrent class has a unique stationary
+distribution whose support equals that recurrent class. The set of all
+stationary distributions is given by the convex hull of these unique
+stationary distributions for the recurrent classes.
+
+A natural number :math:`d` is the *period* of state :math:`i` if it is
+the greatest common divisor of all :math:`k`'s such that :math:`P^k[i,
+i] > 0`; equivalently, it is the GCD of the lengths of the cycles in
+:math:`\Gamma(P)` passing through :math:`i`. For any :math:`i, j`
+contained in a communication class, :math:`i` has period :math:`d` if
+and only if :math:`j` has period :math:`d`. The *period* of an
+irreducible Markov chain (or of an irreducible stochastic matrix) is the
+period of any state. An irreducible Markov chain is *aperiodic* if its
+period is one. A Markov chain is irreducible and aperiodic if and only
+if it is *uniformly ergodic*, i.e., there exists some :math:`m` such
+that :math:`P^m[i, j] > 0` for all :math:`i, j` (in this case, :math:`P`
+is also called *primitive*).
+
+Suppose that an irreducible Markov chain has period :math:`d`. Fix any
+state, say state :math:`0`. For each :math:`m = 0, \ldots, d-1`, let
+:math:`S_m` be the set of states :math:`i` such that :math:`P^{kd+m}[0,
+i] > 0` for some :math:`k`. These sets :math:`S_0, \ldots, S_{d-1}`
+constitute a partition of the state space and are called the *cyclic
+classes*. For each :math:`S_m` and each :math:`i \in S_m`, we have
+:math:`\sum_{j \in S_{m+1}} P[i, j] = 1`, where :math:`S_d = S_0`.
 
 """
 from __future__ import division
 import numpy as np
-import scipy.linalg as la
-import sympy.mpmath as mp
 import sys
 from .discrete_rv import DiscreteRV
+from .graph_tools import DiGraph
+from .gth_solve import gth_solve
 from warnings import warn
 
 
-class DMarkov(object):
+class MarkovChain(object):
     """
-    This class is used as a container for a discrete Markov transition
-    matrix or a discrete Markov chain.  It stores useful information
-    such as the stationary distributions and allows simulation using a
-    specified initial distribution.
-
+    Class for a finite-state discrete-time Markov chain. It stores
+    useful information such as the stationary distributions, and
+    communication, recurrent, and cyclic classes, and allows simulation
+    of state transitions.
 
     Parameters
     ----------
     P : array_like(float, ndim=2)
         The transition matrix.  Must be of shape n x n.
-    pi_0 : array_like(float, ndim=1), optional(default=None)
-        The initial probability distribution.  If no intial distribution
-        is specified, then it will be a distribution with equally
-        probability on each state
-
 
     Attributes
     ----------
     P : array_like(float, ndim=2)
-        The transition matrix.  Must be of shape n x n.
-    pi_0 : array_like(float, ndim=1)
-        The initial probability distribution.
-    stationary_dists : array_like(float, ndim=2)
-        An array with invariant distributions as columns
+        The transition matrix.
+
+    stationary_distributions : array_like(float, ndim=2)
+        Array containing stationary distributions, one for each
+        recurrent class, as rows.
+
+    is_irreducible : bool
+        Indicate whether the Markov chain is irreducible.
+
+    num_communication_classes : int
+        The number of the communication classes.
+
+    communication_classes : list(ndarray(int))
+        List of numpy arrays containing the communication classes.
+
+    num_recurrent_classes : int
+        The number of the recurrent classes.
+
+    recurrent_classes : list(ndarray(int))
+        List of numpy arrays containing the recurrent classes.
+
+    is_aperiodic : bool
+        Indicate whether the Markov chain is aperiodic. Defined only
+        when the Markov chain is irreducible.
 
+    period : int
+        The period of the Markov chain. Defined only when the Markov
+        chain is irreducible.
+
+    cyclic_classes : list(ndarray(int))
+        List of numpy arrays containing the cyclic classes. Defined only
+        when the Markov chain is irreducible.
 
     Methods
     -------
-    mc_compute_stationary : This method finds stationary
-                                    distributions
+    simulate : Simulates the markov chain for a given initial state or
+        distribution.
 
-    mc_sample_path : Simulates the markov chain for a given
-                      initial distribution
     """
 
-    def __init__(self, P, pi_0=None):
-        self.P = P
-        n, m = P.shape
-        self.n = n
+    def __init__(self, P):
+        self.P = np.asarray(P)
 
-        if pi_0 is None:
-            self.pi_0 = np.ones(n)/n
+        # Check Properties
+        # Double check that P is a square matrix
+        if len(self.P.shape) != 2 or self.P.shape[0] != self.P.shape[1]:
+            raise ValueError('P must be a square matrix')
 
-        else:
-            self.pi_0 = pi_0
+        # Double check that P is a nonnegative matrix
+        if not np.all(self.P >= 0):
+            raise ValueError('P must be nonnegative')
 
-        self.stationary_dists = None
+        # Double check that the rows of P sum to one
+        if not np.allclose(np.sum(self.P, axis=1), np.ones(self.P.shape[0])):
+            raise ValueError('The rows of P must sum to 1')
 
-        # Check Properties
-        # double check that P is a square matrix
-        if n != m:
-            raise ValueError('The transition matrix must be square!')
+        # The number of states
+        self.n = self.P.shape[0]
 
-        # Double check that the rows of P sum to one
-        if np.all(np.sum(P, axis=1) != np.ones(P.shape[0])):
-            raise ValueError('The rows must sum to 1. P is a trans matrix')
+        # To analyze the structure of P as a directed graph
+        self.digraph = DiGraph(P)
+
+        self._stationary_dists = None
 
     def __repr__(self):
-        msg = "Markov process with transition matrix \n P = \n {0}"
+        msg = "Markov chain with transition matrix \nP = \n{0}"
 
-        if self.stationary_dists is None:
+        if self._stationary_dists is None:
             return msg.format(self.P)
         else:
-            msg = msg + "\nand stationary distributions \n {1}"
-            return msg.format(self.P, self.stationary_dists)
+            msg = msg + "\nand stationary distributions \n{1}"
+            return msg.format(self.P, self._stationary_dists)
 
     def __str__(self):
         return str(self.__repr__)
 
-    def mc_compute_stationary(self, precision=None, tol=None):
-        P = self.P
-
-        stationary_dists = mc_compute_stationary(P, precision=precision, tol=tol)
-        self.stationary_dists = stationary_dists
-
-        return stationary_dists
-
-    def mc_sample_path(self, init=0, sample_size=1000):
-        sim = mc_sample_path(self.P, init, sample_size)
-
-        return sim
+    @property
+    def is_irreducible(self):
+        return self.digraph.is_strongly_connected
 
+    @property
+    def num_communication_classes(self):
+        return self.digraph.num_strongly_connected_components
 
-def mc_compute_stationary(P, precision=None, tol=None):
-    n = P.shape[0]
+    @property
+    def communication_classes(self):
+        return self.digraph.strongly_connected_components
 
-    if precision is None:
-        # Compute eigenvalues and eigenvectors
-        eigvals, eigvecs = la.eig(P, left=True, right=False)
+    @property
+    def num_recurrent_classes(self):
+        return self.digraph.num_sink_strongly_connected_components
 
-        # Find the index for unit eigenvalues
-        index = np.where(abs(eigvals - 1.) < 1e-12)[0]
+    @property
+    def recurrent_classes(self):
+        return self.digraph.sink_strongly_connected_components
 
-        # Pull out the eigenvectors that correspond to unit eig-vals
-        uniteigvecs = eigvecs[:, index]
+    @property
+    def is_aperiodic(self):
+        if not self.is_irreducible:
+            raise NotImplementedError(
+                'Not defined for a reducible Markov chain'
+            )
+        else:
+            return self.digraph.is_aperiodic
+
+    @property
+    def period(self):
+        if not self.is_irreducible:
+            raise NotImplementedError(
+                'Not defined for a reducible Markov chain'
+            )
+        else:
+            return self.digraph.period
+
+    @property
+    def cyclic_classes(self):
+        if not self.is_irreducible:
+            raise NotImplementedError(
+                'Not defined for a reducible Markov chain'
+            )
+        else:
+            return self.digraph.cyclic_components
 
-        stationary_dists = uniteigvecs/np.sum(uniteigvecs, axis=0)
+    def _compute_stationary(self):
+        """
+        Store the stationary distributions in self._stationary_distributions.
 
-    else:
-        # Create a list to store eigvals
-        stationary_dists_list = []
-        if tol is None:
-            # If tolerance isn't specified then use 2*precision
-            tol = mp.mpf(2 * 10**(-precision + 1))
+        """
+        if self.is_irreducible:
+            stationary_dists = gth_solve(self.P).reshape(1, self.n)
+        else:
+            rec_classes = self.recurrent_classes
+            stationary_dists = np.zeros((len(rec_classes), self.n))
+            for i, rec_class in enumerate(rec_classes):
+                stationary_dists[i, rec_class] = \
+                    gth_solve(self.P[rec_class, :][:, rec_class])
 
-        with mp.workdps(precision):
-            eigvals, eigvecs = mp.eig(mp.matrix(P), left=True, right=False)
+        self._stationary_dists = stationary_dists
 
-            for ind, el in enumerate(eigvals):
-                if el>=(mp.mpf(1)-mp.mpf(tol)) and el<=(mp.mpf(1)+mp.mpf(tol)):
-                    stationary_dists_list.append(eigvecs[ind, :])
+    @property
+    def stationary_distributions(self):
+        if self._stationary_dists is None:
+            self._compute_stationary()
+        return self._stationary_dists
 
-            stationary_dists = np.asarray(stationary_dists_list).T
-            stationary_dists = (stationary_dists/sum(stationary_dists)).astype(np.float)
+    def simulate(self, init=0, sample_size=1000):
+        X = mc_sample_path(self.P, init, sample_size)
 
+        return X
 
-    # Check to make sure all of the elements of invar_dist are positive
-    if np.any(stationary_dists < -1e-16):
-        warn("Elements of your invariant distribution were negative; " +
-              "Re-trying with additional precision")
 
-        if precision is None:
-            stationary_dists = mc_compute_stationary(P, precision=18, tol=tol)
+def mc_compute_stationary(P):
+    """
+    Computes stationary distributions of P, one for each recurrent
+    class. Any stationary distribution is written as a convex
+    combination of these distributions.
 
-        elif precision is not None:
-            raise ValueError("Elements of your stationary distribution were" +
-                             "negative.  Try computing with higher precision")
+    Returns
+    -------
+    stationary_dists : array_like(float, ndim=2)
+        Array containing the stationary distributions as its rows.
 
-    # Since we will be accessing the columns of this matrix, we
-    # might consider adding .astype(np.float, order='F') to make it
-    # column major at beginning
-    return stationary_dists.squeeze()
+    """
+    return MarkovChain(P).stationary_distributions
 
 
 def mc_sample_path(P, init=0, sample_size=1000):
@@ -167,7 +283,7 @@ def mc_sample_path(P, init=0, sample_size=1000):
     # In particular, let P_dist[i] be the distribution corresponding to the
     # i-th row P[i,:]
     n = len(P)
-    P_dist = [DiscreteRV(P[i,:]) for i in range(n)]
+    P_dist = [DiscreteRV(P[i, :]) for i in range(n)]
 
     # === generate the sample path === #
     for t in range(sample_size - 1):
@@ -180,49 +296,18 @@ def mc_sample_path(P, init=0, sample_size=1000):
 # Set up the docstrings for the functions
 #---------------------------------------------------------------------#
 
-# For computing stationary dist
-_stationary_docstr = \
-"""
-This method computes the stationary distributions of P.
-These are the eigenvectors that correspond to the unit eigen-
-values of the matrix P' (They satisfy pi_{{t+1}}' = pi_{{t}}' P).  It
-simply calls the outer function mc_compute_stationary
-
-Parameters
-----------
-{p_arg}init : array_like(float, ndim=2)
-    The discrete Markov transition matrix P
-precision : scalar(int), optional(default=None)
-    Specifies the precision(number of digits of precision) with
-    which to calculate the eigenvalues.  Unless your matrix has
-    multiple eigenvalues that are near unity then no need to
-    worry about this.
-tol : scalar(float), optional(default=None)
-    Specifies the bandwith of eigenvalues to consider equivalent
-    to unity.  It will consider all eigenvalues in [1-tol,
-    1+tol] to be 1.  If tol is None then will use 2*1e-
-    precision.  Only used if precision is defined
-
-Returns
--------
-stationary_dists : np.ndarray : float
-    This is an array that has the stationary distributions as
-    its columns.
-
-"""
-
 # For drawing a sample path
 _sample_path_docstr = \
 """
-Generates one sample path from a finite Markov chain with (n x n)
-Markov matrix P on state space S = {{0,...,n-1}}.
+Generates one sample path from the Markov chain represented by (n x n)
+transition matrix P on state space S = {{0,...,n-1}}.
 
 Parameters
 ----------
 {p_arg}init : array_like(float ndim=1) or scalar(int)
-    If init is an array_like then it is treated as the initial
-    distribution across states.  If init is a scalar then it
-    treated as the deterministic initial state.
+    If init is an array_like, then it is treated as the initial
+    distribution across states.  If init is a scalar, then it treated as
+    the deterministic initial state.
 
 sample_size : scalar(int), optional(default=1000)
     The length of the sample path.
@@ -230,27 +315,20 @@ def mc_sample_path(P, init=0, sample_size=1000):
 Returns
 -------
 X : array_like(int, ndim=1)
-    The simulation of states
+    The simulation of states.
 
 """
 
 # set docstring for functions
-mc_compute_stationary.__doc__ = _stationary_docstr.format(p_arg=
-"""P : array_like(float, ndim=2)
-    A discrete Markov transition matrix
-
-""")
 mc_sample_path.__doc__ = _sample_path_docstr.format(p_arg=
 """P : array_like(float, ndim=2)
-    A discrete Markov transition matrix
+    A Markov transition matrix.
 
 """)
 
 # set docstring for methods
 
 if sys.version_info[0] == 3:
-    DMarkov.mc_compute_stationary.__doc__ = _stationary_docstr.format(p_arg="")
-    DMarkov.mc_sample_path.__doc__ = _sample_path_docstr.format(p_arg="")
+    MarkovChain.simulate.__doc__ = _sample_path_docstr.format(p_arg="")
 elif sys.version_info[0] == 2:
-    DMarkov.mc_compute_stationary.__func__.__doc__ = _stationary_docstr.format(p_arg="")
-    DMarkov.mc_sample_path.__func__.__doc__ = _sample_path_docstr.format(p_arg="")
+    MarkovChain.simulate.__func__.__doc__ = _sample_path_docstr.format(p_arg="")
diff --git a/quantecon/tests/test_graph_tools.py b/quantecon/tests/test_graph_tools.py
new file mode 100644
index 000000000..505fa861b
--- /dev/null
+++ b/quantecon/tests/test_graph_tools.py
@@ -0,0 +1,197 @@
+"""
+Filename: test_graph_tools.py
+Author: Daisuke Oyama
+
+Tests for graph_tools.py
+
+"""
+import sys
+import numpy as np
+from numpy.testing import assert_array_equal
+import nose
+from nose.tools import eq_, ok_, raises
+
+from quantecon.graph_tools import DiGraph
+
+
+class Graphs:
+    """Setup graphs for the tests"""
+
+    def __init__(self):
+        self.strongly_connected_graph_dicts = []
+        self.not_strongly_connected_graph_dicts = []
+
+        graph_dict = {
+            'A': np.array([[1, 0], [0, 1]]),
+            'strongly_connected_components':
+            [np.array([0]), np.array([1])],
+            'sink_strongly_connected_components':
+            [np.array([0]), np.array([1])],
+            'is_strongly_connected': False,
+        }
+        self.not_strongly_connected_graph_dicts.append(graph_dict)
+
+        graph_dict = {
+            'A': np.array([[1, 0, 0], [1, 0, 1], [0, 0, 1]]),
+            'strongly_connected_components':
+            [np.array([0]), np.array([1]), np.array([2])],
+            'sink_strongly_connected_components':
+            [np.array([0]), np.array([2])],
+            'is_strongly_connected': False,
+        }
+        self.not_strongly_connected_graph_dicts.append(graph_dict)
+
+        graph_dict = {
+            'A': np.array([[1, 1], [1, 1]]),
+            'strongly_connected_components': [np.arange(2)],
+            'sink_strongly_connected_components': [np.arange(2)],
+            'is_strongly_connected': True,
+            'period': 1,
+            'is_aperiodic': True,
+            'cyclic_components': [np.arange(2)],
+        }
+        self.strongly_connected_graph_dicts.append(graph_dict)
+
+        graph_dict = {
+            'A': np.array([[0, 1], [1, 0]]),
+            'strongly_connected_components': [np.arange(2)],
+            'sink_strongly_connected_components': [np.arange(2)],
+            'is_strongly_connected': True,
+            'period': 2,
+            'is_aperiodic': False,
+            'cyclic_components': [np.array([0]), np.array([1])],
+        }
+        self.strongly_connected_graph_dicts.append(graph_dict)
+
+        graph_dict = {
+            'A': np.array([[0, 1, 0, 0],
+                           [0, 0, 1, 0],
+                           [1, 0, 0, 1],
+                           [0, 0, 1, 0]]),
+            'strongly_connected_components': [np.arange(4)],
+            'sink_strongly_connected_components': [np.arange(4)],
+            'is_strongly_connected': True,
+            'period': 1,
+            'is_aperiodic': True,
+            'cyclic_components': [np.arange(4)],
+        }
+        self.strongly_connected_graph_dicts.append(graph_dict)
+
+        # Weighted graph
+        graph_dict = {
+            'A': np.array([[0, 0.5], [2, 0]]),
+            'weighted': True,
+            'strongly_connected_components': [np.arange(2)],
+            'sink_strongly_connected_components': [np.arange(2)],
+            'is_strongly_connected': True,
+            'period': 2,
+            'is_aperiodic': False,
+            'cyclic_components': [np.array([0]), np.array([1])],
+        }
+        self.strongly_connected_graph_dicts.append(graph_dict)
+
+        # Degenrate graph with no edge
+        graph_dict = {
+            'A': np.array([[0]]),
+            'strongly_connected_components': [np.arange(1)],
+            'sink_strongly_connected_components': [np.arange(1)],
+            'is_strongly_connected': True,
+            'period': 1,
+            'is_aperiodic': True,
+            'cyclic_components': [np.array([0])],
+        }
+        self.strongly_connected_graph_dicts.append(graph_dict)
+
+        # Degenrate graph with self loop
+        graph_dict = {
+            'A': np.array([[1]]),
+            'strongly_connected_components': [np.arange(1)],
+            'sink_strongly_connected_components': [np.arange(1)],
+            'is_strongly_connected': True,
+            'period': 1,
+            'is_aperiodic': True,
+            'cyclic_components': [np.array([0])],
+        }
+        self.strongly_connected_graph_dicts.append(graph_dict)
+
+        self.graph_dicts = \
+            self.strongly_connected_graph_dicts + \
+            self.not_strongly_connected_graph_dicts
+
+
+class TestDiGraph:
+    """Test the methods in Digraph"""
+
+    def setUp(self):
+        """Setup Digraph instances"""
+        self.graphs = Graphs()
+        for graph_dict in self.graphs.graph_dicts:
+            try:
+                weighted = graph_dict['weighted']
+            except:
+                weighted = False
+            graph_dict['g'] = DiGraph(graph_dict['A'], weighted=weighted)
+
+    def test_strongly_connected_components(self):
+        for graph_dict in self.graphs.graph_dicts:
+            assert_array_equal(
+                sorted(graph_dict['g'].strongly_connected_components),
+                sorted(graph_dict['strongly_connected_components']))
+
+    def test_num_strongly_connected_components(self):
+        for graph_dict in self.graphs.graph_dicts:
+            eq_(graph_dict['g'].num_strongly_connected_components,
+                len(graph_dict['strongly_connected_components']))
+
+    def test_sink_strongly_connected_components(self):
+        for graph_dict in self.graphs.graph_dicts:
+            assert_array_equal(
+                sorted(graph_dict['g'].sink_strongly_connected_components),
+                sorted(graph_dict['sink_strongly_connected_components']))
+
+    def test_num_sink_strongly_connected_components(self):
+        for graph_dict in self.graphs.graph_dicts:
+            eq_(graph_dict['g'].num_sink_strongly_connected_components,
+                len(graph_dict['sink_strongly_connected_components']))
+
+    def test_is_strongly_connected(self):
+        for graph_dict in self.graphs.graph_dicts:
+            eq_(graph_dict['g'].is_strongly_connected,
+                graph_dict['is_strongly_connected'])
+
+    def test_period(self):
+        for graph_dict in self.graphs.graph_dicts:
+            try:
+                eq_(graph_dict['g'].period, graph_dict['period'])
+            except NotImplementedError:
+                eq_(graph_dict['g'].is_strongly_connected, False)
+
+    def test_is_aperiodic(self):
+        for graph_dict in self.graphs.graph_dicts:
+            try:
+                eq_(graph_dict['g'].is_aperiodic,
+                    graph_dict['is_aperiodic'])
+            except NotImplementedError:
+                eq_(graph_dict['g'].is_strongly_connected, False)
+
+    def test_cyclic_components(self):
+        for graph_dict in self.graphs.graph_dicts:
+            try:
+                assert_array_equal(
+                    sorted(graph_dict['g'].cyclic_components),
+                    sorted(graph_dict['cyclic_components']))
+            except NotImplementedError:
+                eq_(graph_dict['g'].is_strongly_connected, False)
+
+
+@raises(ValueError)
+def test_raises_value_error_non_sym():
+    """Test with non symmetric input"""
+    g = DiGraph(np.array([[0.4, 0.6]]))
+
+
+if __name__ == '__main__':
+    argv = sys.argv[:]
+    argv.append('--verbose')
+    argv.append('--nocapture')
+    nose.main(argv=argv, defaultTest=__file__)
diff --git a/quantecon/tests/test_gth_solve.py b/quantecon/tests/test_gth_solve.py
new file mode 100644
index 000000000..fb7ad06e8
--- /dev/null
+++ b/quantecon/tests/test_gth_solve.py
@@ -0,0 +1,173 @@
+"""
+Filename: gth_solve.py
+Author: Daisuke Oyama
+
+Tests for gth_solve.py
+
+"""
+from __future__ import division
+
+import sys
+import numpy as np
+import nose
+from nose.tools import eq_, ok_, raises
+
+from quantecon.gth_solve import gth_solve
+
+
+TOL = 1e-15
+
+
+def KMR_Markov_matrix_sequential(N, p, epsilon):
+    """
+    Generate the Markov matrix for the KMR model with *sequential* move
+
+    Parameters
+    ----------
+    N : int
+        Number of players
+
+    p : float
+        Level of p-dominance of action 1, i.e.,
+        the value of p such that action 1 is the BR for (1-q, q) for any q > p,
+        where q (1-q, resp.) is the prob that the opponent plays action 1 (0, resp.)
+
+    epsilon : float
+        Probability of mutation
+
+    Returns
+    -------
+    P : numpy.ndarray
+        Markov matrix for the KMR model with simultaneous move
+
+    """
+    P = np.zeros((N+1, N+1), dtype=float)
+    P[0, 0], P[0, 1] = 1 - epsilon * (1/2), epsilon * (1/2)
+    for n in range(1, N):
+        P[n, n-1] = \
+            (n/N) * (epsilon * (1/2) +
+                     (1 - epsilon) * (((n-1)/(N-1) < p) + ((n-1)/(N-1) == p) * (1/2))
+                     )
+        P[n, n+1] = \
+            ((N-n)/N) * (epsilon * (1/2) +
+                         (1 - epsilon) * ((n/(N-1) > p) + (n/(N-1) == p) * (1/2))
+                         )
+        P[n, n] = 1 - P[n, n-1] - P[n, n+1]
+    P[N, N-1], P[N, N] = epsilon * (1/2), 1 - epsilon * (1/2)
+    return P
+
+
+class Matrices:
+    """Setup matrices for the tests"""
+
+    def __init__(self):
+        self.stoch_matrix_dicts = []
+        self.kmr_matrix_dicts = []
+        self.gen_matrix_dicts = []
+
+        matrix_dict = {
+            'A': np.array([[0.4, 0.6], [0.2, 0.8]]),
+            'stationary_dist': np.array([[0.25, 0.75]]),
+        }
+        self.stoch_matrix_dicts.append(matrix_dict)
+
+        matrix_dict = {
+            # Reducible matrix
+            'A': np.array([[1, 0], [0, 1]]),
+            # Stationary dist whose support contains index 0
+            'stationary_dist': np.array([[1, 0]]),
+        }
+        self.stoch_matrix_dicts.append(matrix_dict)
+
+        matrix_dict = {
+            'A': KMR_Markov_matrix_sequential(N=27, p=1./3, epsilon=1e-2),
+        }
+        self.kmr_matrix_dicts.append(matrix_dict)
+
+        matrix_dict = {
+            'A': KMR_Markov_matrix_sequential(N=3, p=1./3, epsilon=1e-14),
+        }
+        self.kmr_matrix_dicts.append(matrix_dict)
+
+        matrix_dict = {
+            'A': np.array([[-1, 1], [4, -4]]),
+            'stationary_dist': np.array([[0.8, 0.2]]),
+        }
+        self.gen_matrix_dicts.append(matrix_dict)
+
+
+def test_stoch_matrix():
+    """Test with stochastic matrices"""
+    print(__name__ + '.' + test_stoch_matrix.__name__)
+    matrices = Matrices()
+    for matrix_dict in matrices.stoch_matrix_dicts:
+        x = gth_solve(matrix_dict['A'])
+        yield StationaryDistSumOne(), x
+        yield StationaryDistNonnegative(), x
+        yield StationaryDistEqualToKnown(), matrix_dict['stationary_dist'], x
+
+
+def test_kmr_matrix():
+    """Test with KMR matrices"""
+    print(__name__ + '.' + test_kmr_matrix.__name__)
+    matrices = Matrices()
+    for matrix_dict in matrices.kmr_matrix_dicts:
+        x = gth_solve(matrix_dict['A'])
+        yield StationaryDistSumOne(), x
+        yield StationaryDistNonnegative(), x
+        yield StationaryDistLeftEigenVec(), matrix_dict['A'], x
+
+
+def test_gen_solve():
+    """Test with generator matrices"""
+    print(__name__ + '.' + test_gen_solve.__name__)
+    matrices = Matrices()
+    for matrix_dict in matrices.gen_matrix_dicts:
+        x = gth_solve(matrix_dict['A'])
+        yield StationaryDistSumOne(), x
+        yield StationaryDistNonnegative(), x
+        yield StationaryDistEqualToKnown(), matrix_dict['stationary_dist'], x
+
+
+class AddDescription:
+    def __init__(self):
+        self.description = self.__class__.__name__
+
+
+class StationaryDistSumOne(AddDescription):
+    def __call__(self, x):
+        ok_(np.allclose(sum(x), 1, atol=TOL))
+
+
+class StationaryDistNonnegative(AddDescription):
+    def __call__(self, x):
+        eq_(np.prod(x >= 0-TOL), 1)
+
+
+class StationaryDistLeftEigenVec(AddDescription):
+    def __call__(self, A, x):
+        ok_(np.allclose(np.dot(x, A), x, atol=TOL))
+
+
+class StationaryDistEqualToKnown(AddDescription):
+    def __call__(self, y, x):
+        ok_(np.allclose(y, x, atol=TOL))
+
+
+@raises(ValueError)
+def test_raises_value_error_non_2dim():
+    """Test with non 2dim input"""
+    gth_solve(np.array([0.4, 0.6]))
+
+
+@raises(ValueError)
+def test_raises_value_error_non_sym():
+    """Test with non symmetric input"""
+    gth_solve(np.array([[0.4, 0.6]]))
+
+
+if __name__ == '__main__':
+    argv = sys.argv[:]
+    argv.append('--verbose')
+    argv.append('--nocapture')
+    nose.main(argv=argv, defaultTest=__file__)
diff --git a/quantecon/tests/test_mc_tools.py b/quantecon/tests/test_mc_tools.py
index 540ea02f0..21eba25ef 100644
--- a/quantecon/tests/test_mc_tools.py
+++ b/quantecon/tests/test_mc_tools.py
@@ -3,18 +3,20 @@
 
 Functions
 ---------
-    mc_compute_stationary   [Status: 1 x Simple Test Written]
+    mc_compute_stationary   [Status: Tested in test_markovchain_pmatrices]
     mc_sample_path          [Status: TBD]
 
 """
-
 from __future__ import division
 
+import sys
 import numpy as np
-import unittest
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_array_equal
+import nose
+from nose.tools import raises
+
+from quantecon.mc_tools import MarkovChain, mc_compute_stationary
 
-from quantecon.mc_tools import DMarkov, mc_compute_stationary, mc_sample_path
 
 # KMR Function
 # Useful because it seems to have 1 unit eigvalue, but a second one that
@@ -47,37 +49,83 @@ def KMR_Markov_matrix_sequential(N, p, epsilon):
     P[N, N-1], P[N, N] = epsilon * (1/2), 1 - epsilon * (1/2)
     return P
 
-### Tests: mc_compute_stationary ###
 
-def test_mc_compute_stationary_pmatrices():
+# Tests: methods of MarkovChain, mc_compute_stationary #
+
+def test_markovchain_pmatrices():
     """
-        Test mc_compute_stationary with P Matrix and Known Solutions
+    Test the methods of MarkovChain, as well as mc_compute_stationary,
+    with P matrix and known solutions
     """
-
-                    #-P Matrix-#                        , #-Known Solution-#
-    testset =   [
-                    ( np.array([[0.4,0.6], [0.2,0.8]]), np.array([0.25, 0.75]) ),
-                    ( np.eye(2), np.eye(2) )
-                ]
-
-    #-Loop Through TestSet-#
-    for (P, known) in testset:
-        computed = mc_compute_stationary(P)
-        assert_allclose(computed, known)
-
-
-
+    testset = [
+        {'P': np.array([[0.4, 0.6], [0.2, 0.8]]),  # P matrix
+         'stationary_dists': np.array([[0.25, 0.75]]),  # Known solution
+         'comm_classes': [np.arange(2)],
+         'rec_classes': [np.arange(2)],
+         'is_irreducible': True,
+         'period': 1,
+         'is_aperiodic': True,
+         'cyclic_classes': [np.arange(2)],
+         },
+        {'P': np.array([[0, 1], [1, 0]]),
+         'stationary_dists': np.array([[0.5, 0.5]]),
+         'comm_classes': [np.arange(2)],
+         'rec_classes': [np.arange(2)],
+         'is_irreducible': True,
+         'period': 2,
+         'is_aperiodic': False,
+         'cyclic_classes': [np.array([0]), np.array([1])],
+         },
+        {'P': np.eye(2),
+         'stationary_dists': np.array([[1, 0], [0, 1]]),
+         'comm_classes': [np.array([0]), np.array([1])],
+         'rec_classes': [np.array([0]), np.array([1])],
+         'is_irreducible': False,
+         }
+    ]
+
+    # Loop Through TestSet #
+    for test_dict in testset:
+        mc = MarkovChain(test_dict['P'])
+        computed = mc.stationary_distributions
+        assert_allclose(computed, test_dict['stationary_dists'])
+
+        assert(mc.num_communication_classes == len(test_dict['comm_classes']))
+        assert(mc.is_irreducible == test_dict['is_irreducible'])
+        assert(mc.num_recurrent_classes == len(test_dict['rec_classes']))
+        assert_array_equal(sorted(mc.communication_classes),
+                           sorted(test_dict['comm_classes']))
+        assert_array_equal(sorted(mc.recurrent_classes),
+                           sorted(test_dict['rec_classes']))
+        try:
+            assert(mc.period == test_dict['period'])
+        except NotImplementedError:
+            assert(mc.is_irreducible is False)
+        try:
+            assert(mc.is_aperiodic == test_dict['is_aperiodic'])
+        except NotImplementedError:
+            assert(mc.is_irreducible is False)
+        try:
+            assert_array_equal(sorted(mc.cyclic_classes),
+                               sorted(test_dict['cyclic_classes']))
+        except NotImplementedError:
+            assert(mc.is_irreducible is False)
+
+        # Test of mc_compute_stationary
+        computed = mc_compute_stationary(test_dict['P'])
+        assert_allclose(computed, test_dict['stationary_dists'])
 
 
 # Basic Class Structure with Setup #
 ####################################
 
-class Test_mc_compute_stationary_KMRMarkovMatrix2():
+class Test_markovchain_stationary_distributions_KMRMarkovMatrix2():
     """
-    Test Suite for mc_compute_stationary using KMR Markov Matrix [suitable for nose]
+    Test Suite for MarkovChain.stationary_distributions using KMR Markov
+    Matrix [suitable for nose]
     """
 
-    #-Starting Values-#
+    # Starting Values #
 
     N = 27
     epsilon = 1e-2
@@ -87,15 +135,14 @@ class Test_mc_compute_stationary_KMRMarkovMatrix2():
     def setUp(self):
         """ Setup a KMRMarkovMatrix and Compute Stationary Values """
         self.P = KMR_Markov_matrix_sequential(self.N, self.p, self.epsilon)
-        self.mc = DMarkov(self.P)
-        self.stationary = self.mc.mc_compute_stationary()
+        self.mc = MarkovChain(self.P)
+        self.stationary = self.mc.stationary_distributions
         stat_shape = self.stationary.shape
 
-        if len(stat_shape) is 1:
+        if len(stat_shape) == 1:
             self.n_stat_dists = 1
         else:
-            self.n_stat_dists = stat_shape[1]
-
+            self.n_stat_dists = stat_shape[0]
 
     def test_markov_matrix(self):
         "Check that each row of matrix sums to 1"
@@ -103,9 +150,9 @@ def test_markov_matrix(self):
         assert_allclose(np.sum(mc.P, axis=1), np.ones(mc.n))
 
     def test_sum_one(self):
-        "Check each stationary distribution sums to 1"
+        "Check that each stationary distribution sums to 1"
         stationary_distributions = self.stationary
-        assert_allclose(np.sum(stationary_distributions, axis=0),
+        assert_allclose(np.sum(stationary_distributions, axis=1),
                         np.ones(self.n_stat_dists))
 
     def test_nonnegative(self):
@@ -116,13 +163,43 @@ def test_nonnegative(self):
     def test_left_eigen_vec(self):
         "Check that vP = v for all stationary distributions"
         mc = self.mc
-        stationary = self.stationary
+        stationary_distributions = self.stationary
 
-        if self.n_stat_dists is 1:
-            assert_allclose(np.dot(stationary, mc.P), stationary, atol=self.TOL)
+        if self.n_stat_dists == 1:
+            assert_allclose(np.dot(stationary_distributions, mc.P),
+                            stationary_distributions, atol=self.TOL)
         else:
             for i in range(self.n_stat_dists):
-                curr_v = stationary_distributions[:, i]
+                curr_v = stationary_distributions[i, :]
                 assert_allclose(np.dot(curr_v, mc.P), curr_v, atol=self.TOL)
 
 
+@raises(ValueError)
+def test_raises_value_error_non_2dim():
+    """Test with non 2dim input"""
+    mc = MarkovChain(np.array([0.4, 0.6]))
+
+
+@raises(ValueError)
+def test_raises_value_error_non_sym():
+    """Test with non symmetric input"""
+    mc = MarkovChain(np.array([[0.4, 0.6]]))
+
+
+@raises(ValueError)
+def test_raises_value_error_non_nonnegative():
+    """Test with non nonnegative input"""
+    mc = MarkovChain(np.array([[0.4, 0.6], [-0.2, 1.2]]))
+
+
+@raises(ValueError)
+def test_raises_value_error_non_sum_one():
+    """Test with input such that some of the rows does not sum to one"""
+    mc = MarkovChain(np.array([[0.4, 0.6], [0.2, 0.9]]))
+
+
+if __name__ == '__main__':
+    argv = sys.argv[:]
+    argv.append('--verbose')
+    argv.append('--nocapture')
+    nose.main(argv=argv, defaultTest=__file__)