diff --git a/quantecon/markov/approximation.py b/quantecon/markov/approximation.py
index 2256fadba..2cfb809b6 100644
--- a/quantecon/markov/approximation.py
+++ b/quantecon/markov/approximation.py
@@ -8,15 +8,16 @@
 from math import erfc, sqrt
 from .core import MarkovChain
 
+import warnings
 import numpy as np
 from numba import njit
 
 
-def rouwenhorst(n, ybar, sigma, rho):
+def rouwenhorst(n, rho, sigma, mu=0.):
     r"""
-    Takes as inputs n, p, q, psi. It will then construct a markov chain
+    Takes as inputs n, mu, sigma, rho. It will then construct a markov chain
     that estimates an AR(1) process of:
-    :math:`y_t = \bar{y} + \rho y_{t-1} + \varepsilon_t`
+    :math:`y_t = \mu + \rho y_{t-1} + \varepsilon_t`
     where :math:`\varepsilon_t` is i.i.d. normal of mean 0, std dev of sigma
 
     The Rouwenhorst approximation uses the following recursive defintion
@@ -54,22 +55,23 @@ def rouwenhorst(n, ybar, sigma, rho):
         0  & \theta_n  \\
         \end{bmatrix}
 
+    where :math:`{p = q = \frac{(1 + \rho)}{2}}`
 
     Parameters
     ----------
     n : int
         The number of points to approximate the distribution
 
-    ybar : float
-        The value :math:`\bar{y}` in the process.  Note that the mean of this
-        AR(1) process, :math:`y`, is simply :math:`\bar{y}/(1 - \rho)`
+    rho : float
+        Persistence parameter in AR(1) process, if you are approximating
+        an AR(1) process then this is the autocorrelation across periods.
 
     sigma : float
         The value of the standard deviation of the :math:`\varepsilon` process
 
-    rho : float
-        By default this will be 0, but if you are approximating an AR(1)
-        process then this is the autocorrelation across periods
+    mu : float, optional(default=0.0)
+        The value :math:`\mu` in the process.  Note that the mean of this
+        AR(1) process, :math:`y`, is simply :math:`\mu/(1 - \rho)`
 
     Returns
     -------
@@ -78,8 +80,19 @@ def rouwenhorst(n, ybar, sigma, rho):
         An instance of the MarkovChain class that stores the transition
         matrix and state values returned by the discretization method
 
+    Note
+    ----
+
+    UserWarning: The API of `rouwenhorst` was changed from
+    `rouwenhorst(n, ybar, sigma, rho)` to
+    `rouwenhorst(n, rho, sigma, mu=0.)` in version 0.6.0.
+
     """
 
+    warnings.warn("The API of rouwenhorst has changed from `rouwenhorst(n, ybar, sigma, rho)`"
+                  " to `rouwenhorst(n, rho, sigma, mu=0.)`. To find more details please visit:"
+                  " https://github.com/QuantEcon/QuantEcon.py/issues/663.",
+                  UserWarning, stacklevel=2)
     # Get the standard deviation of y
     y_sd = sqrt(sigma**2 / (1 - rho**2))
 
@@ -130,35 +143,37 @@ def row_build_mat(n, p, q):
 
     theta = row_build_mat(n, p, q)
 
-    bar += ybar / (1 - rho)
+    bar += mu / (1 - rho)
 
     return MarkovChain(theta, bar)
 
 
-def tauchen(rho, sigma_u, b=0., m=3, n=7):
+def tauchen(n, rho, sigma, mu=0., n_std=3):
     r"""
     Computes a Markov chain associated with a discretized version of
     the linear Gaussian AR(1) process
 
     .. math::
 
-        y_{t+1} = b + \rho y_t + u_{t+1}
+        y_t = \mu + \rho y_{t-1} + \epsilon_t
 
-    using Tauchen's method. Here :math:`{u_t}` is an i.i.d. Gaussian process
+    using Tauchen's method. Here :math:`{\epsilon_t}` is an i.i.d. Gaussian process
     with zero mean.
 
     Parameters
     ----------
-    b : scalar(float)
-        The constant term of {y_t}
+
+    n : scalar(int)
+        The number of states to use in the approximation
     rho : scalar(float)
-        The autocorrelation coefficient
-    sigma_u : scalar(float)
+        The autocorrelation coefficient, Persistence parameter in AR(1) process
+    sigma : scalar(float)
         The standard deviation of the random process
-    m : scalar(int), optional(default=3)
+    mu : scalar(float), optional(default=0.0)
+        The value :math:`\mu` in the process.  Note that the mean of this
+        AR(1) process, :math:`y`, is simply :math:`\mu/(1 - \rho)`
+    n_std : 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
     -------
@@ -167,13 +182,24 @@ def tauchen(rho, sigma_u, b=0., m=3, n=7):
         An instance of the MarkovChain class that stores the transition
         matrix and state values returned by the discretization method
 
+    Note
+    ----
+
+    UserWarning: The API of `tauchen` was changed from
+    `tauchen(rho, sigma_u, b=0., m=3, n=7)` to
+    `tauchen(n, rho, sigma, mu=0., n_std=3)` in version 0.6.0.
+
     """
+    warnings.warn("The API of tauchen has changed from `tauchen(rho, sigma_u, b=0., m=3, n=7)`"
+                  " to `tauchen(n, rho, sigma, mu=0., n_std=3)`. To find more details please visit:"
+                  " https://github.com/QuantEcon/QuantEcon.py/issues/663.",
+                  UserWarning, stacklevel=2)
 
     # standard deviation of demeaned y_t
-    std_y = np.sqrt(sigma_u**2 / (1 - rho**2))
+    std_y = np.sqrt(sigma**2 / (1 - rho**2))
 
     # top of discrete state space for demeaned y_t
-    x_max = m * std_y
+    x_max = n_std * std_y
 
     # bottom of discrete state space for demeaned y_t
     x_min = -x_max
@@ -187,10 +213,10 @@ def tauchen(rho, sigma_u, b=0., m=3, n=7):
 
     # approximate Markov transition matrix for
     # demeaned y_t
-    _fill_tauchen(x, P, n, rho, sigma_u, half_step)
+    _fill_tauchen(x, P, n, rho, sigma, half_step)
 
     # shifts the state values by the long run mean of y_t
-    mu = b / (1 - rho)
+    mu = mu / (1 - rho)
 
     mc = MarkovChain(P, state_values=x+mu)
 
diff --git a/quantecon/markov/tests/test_approximation.py b/quantecon/markov/tests/test_approximation.py
index e10d5b9b7..0b0e368dc 100644
--- a/quantecon/markov/tests/test_approximation.py
+++ b/quantecon/markov/tests/test_approximation.py
@@ -3,6 +3,7 @@
 
 """
 import numpy as np
+import pytest
 from quantecon.markov import tauchen, rouwenhorst
 from numpy.testing import assert_, assert_allclose
 
@@ -12,13 +13,14 @@
 class TestTauchen:
 
     def setup_method(self):
-        self.rho, self.sigma_u = np.random.rand(2)
+        self.rho, self.sigma = np.random.rand(2)
         self.n = np.random.randint(3, 25)
-        self.m = np.random.randint(5)
+        self.n_std = np.random.randint(5)
         self.tol = 1e-12
-        self.b = 0.
+        self.mu = 0.
 
-        mc = tauchen(self.rho, self.sigma_u, self.b, self.m, self.n)
+        with pytest.warns(UserWarning):
+            mc = tauchen(self.n, self.rho, self.sigma, self.mu, self.n_std)
         self.x, self.P = mc.state_values, mc.P
 
     def teardown_method(self):
@@ -26,10 +28,11 @@ def teardown_method(self):
         del self.P
 
     def testStateCenter(self):
-        for b in [0., 1., -1.]:
-            mu = b / (1 - self.rho)
-            mc = tauchen(self.rho, self.sigma_u, b, self.m, self.n)
-            assert_allclose(mu, np.mean(mc.state_values), atol=self.tol)
+        for mu in [0., 1., -1.]:
+            mu_expect = mu / (1 - self.rho)
+            with pytest.warns(UserWarning):
+                mc = tauchen(self.n, self.rho, self.sigma, mu, self.n_std)
+            assert_allclose(mu_expect, np.mean(mc.state_values), atol=self.tol)
 
     def testShape(self):
         i, j = self.P.shape
@@ -57,10 +60,11 @@ class TestRouwenhorst:
     def setup_method(self):
         self.rho, self.sigma = np.random.uniform(0, 1, size=2)
         self.n = np.random.randint(3, 26)
-        self.ybar = np.random.randint(0, 11)
+        self.mu = np.random.randint(0, 11)
         self.tol = 1e-10
 
-        mc = rouwenhorst(self.n, self.ybar, self.sigma, self.rho)
+        with pytest.warns(UserWarning):
+            mc = rouwenhorst(self.n, self.rho, self.sigma, self.mu)
         self.x, self.P = mc.state_values, mc.P
 
     def teardown_method(self):
@@ -84,12 +88,13 @@ def test_positive_probs(self):
         assert_(np.all(self.P > -self.tol))
 
     def test_states_sum_0(self):
-        tol = self.tol + self.n*(self.ybar/(1 - self.rho))
+        tol = self.tol + self.n*(self.mu/(1 - self.rho))
         assert_(abs(np.sum(self.x)) < tol)
 
     def test_control_case(self):
-        n = 3; ybar = 1; sigma = 0.5; rho = 0.8;
-        mc_rouwenhorst = rouwenhorst(n, ybar, sigma, rho)
+        n = 3; mu = 1; sigma = 0.5; rho = 0.8;
+        with pytest.warns(UserWarning):
+            mc_rouwenhorst = rouwenhorst(n, rho, sigma, mu)
         mc_rouwenhorst.x, mc_rouwenhorst.P = mc_rouwenhorst.state_values, mc_rouwenhorst.P
         sigma_y = np.sqrt(sigma**2 / (1-rho**2))
         psi = sigma_y * np.sqrt(n-1)