r"""
    Computes a Markov chain associated with a discretized version of
    the linear Gaussian AR(1) process

    .. math::

        y_{t+1} = \rho y_t + u_{t+1}

    using Tauchen's method. Here :math:`{u_t}` is an i.i.d. Gaussian process
    with zero mean.

    Parameters
    ----------
    rho : scalar(float)
        The autocorrelation coefficient
    sigma_u : scalar(float)
        The standard deviation of the random process
    m : scalar(int), optional(default=3)
        The number of standard deviations to approximate out to
    n : scalar(int), optional(default=7)
        The number of states to use in the approximation

    Returns
    -------

    mc : MarkovChain
        An instance of the MarkovChain class that stores the transition
        matrix and state values returned by the discretization method

def tauchen(rho, sigma_u, m=3, n=7):
    r"""
    Computes a Markov chain associated with a discretized version of
    the linear Gaussian AR(1) process

    .. math::

        y_{t+1} = \rho y_t + u_{t+1}

    using Tauchen's method. Here :math:`{u_t}` is an i.i.d. Gaussian process
    with zero mean.

    Parameters
    ----------
    rho : scalar(float)
        The autocorrelation coefficient
    sigma_u : scalar(float)
        The standard deviation of the random process
    m : scalar(int), optional(default=3)
        The number of standard deviations to approximate out to
    n : scalar(int), optional(default=7)
        The number of states to use in the approximation

    Returns
    -------

    mc : MarkovChain
        An instance of the MarkovChain class that stores the transition
        matrix and state values returned by the discretization method

    """

    # standard deviation of y_t
    std_y = np.sqrt(sigma_u**2 / (1 - rho**2))

    # top of discrete state space
    x_max = m * std_y

    # bottom of discrete state space
    x_min = -x_max

    # discretized state space
    x = np.linspace(x_min, x_max, n)

    step = (x_max - x_min) / (n - 1)
    half_step = 0.5 * step
    P = np.empty((n, n))

    _fill_tauchen(x, P, n, rho, sigma_u, half_step)

    mc = MarkovChain(P, state_values=x)
    return mc