diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index 9af63a61..63e374a7 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -87,7 +87,7 @@
 
 from .gth_solve import gth_solve
 from .._graph_tools import DiGraph
-from ..util import searchsorted, check_random_state, rng_integers
+from ..util import check_random_state, rng_integers
 
 
 class MarkovChain:
@@ -618,7 +618,7 @@ def _generate_sample_paths(P_cdfs, init_states, random_values, out):
     for i in range(num_reps):
         out[i, 0] = init_states[i]
         for t in range(ts_length-1):
-            out[i, t+1] = searchsorted(P_cdfs[out[i, t]], random_values[i, t])
+            out[i, t+1] = np.searchsorted(P_cdfs[out[i, t]], random_values[i, t], side='right')
 
 
 @jit(nopython=True)
@@ -662,8 +662,8 @@ def _generate_sample_paths_sparse(P_cdfs1d, indices, indptr, init_states,
     for i in range(num_reps):
         out[i, 0] = init_states[i]
         for t in range(ts_length-1):
-            k = searchsorted(P_cdfs1d[indptr[out[i, t]]:indptr[out[i, t]+1]],
-                             random_values[i, t])
+            k = np.searchsorted(P_cdfs1d[indptr[out[i, t]]:indptr[out[i, t]+1]],
+                                random_values[i, t], side='right')
             out[i, t+1] = indices[indptr[out[i, t]]+k]
 
 
@@ -719,7 +719,7 @@ def mc_sample_path(P, init=0, sample_size=1000, random_state=None):
     else:
         cdf0 = np.cumsum(init)
         u_0 = random_state.random()
-        X_0 = searchsorted(cdf0, u_0)
+        X_0 = np.searchsorted(cdf0, u_0, side='right')
 
     mc = MarkovChain(P)
     return mc.simulate(ts_length=sample_size, init=X_0,
diff --git a/quantecon/random/utilities.py b/quantecon/random/utilities.py
index 9838b2e1..bd309b06 100644
--- a/quantecon/random/utilities.py
+++ b/quantecon/random/utilities.py
@@ -6,7 +6,7 @@
 import numpy as np
 from numba import guvectorize, types
 from numba.extending import overload
-from ..util import check_random_state, searchsorted
+from ..util import check_random_state
 
 
 # Generating Arrays and Vectors #
@@ -204,7 +204,7 @@ def draw(cdf, size=None):
         return out
     else:
         r = np.random.random()
-        return searchsorted(cdf, r)
+        return np.searchsorted(cdf, r, side='right')
 
 
 # Overload for the `draw` function
@@ -215,10 +215,10 @@ def draw_impl(cdf, size=None):
             rs = np.random.random(size)
             out = np.empty(size, dtype=np.int_)
             for i in range(size):
-                out[i] = searchsorted(cdf, rs[i])
+                out[i] = np.searchsorted(cdf, rs[i], side='right')
             return out
     else:
         def draw_impl(cdf, size=None):
             r = np.random.random()
-            return searchsorted(cdf, r)
+            return np.searchsorted(cdf, r, side='right')
     return draw_impl
diff --git a/quantecon/util/array.py b/quantecon/util/array.py
index 3821a06a..4200948f 100644
--- a/quantecon/util/array.py
+++ b/quantecon/util/array.py
@@ -8,7 +8,8 @@
 
 """
 
-from numba import jit
+import warnings
+from numba import jit, objmode
 
 # ----------------- #
 # -ARRAY UTILITIES- #
@@ -21,6 +22,10 @@ def searchsorted(a, v):
     that `a[i-1] <= v < a[i]` (for `i = 0`, `v < a[0]`); if `v[n-1] <=
     v`, return `n`, where `n = len(a)`.
 
+    .. deprecated::
+
+        Deprecated, use `np.searchsorted(a, v, side='right')` instead.
+
     Parameters
     ----------
     a : ndarray(float, ndim=1)
@@ -51,6 +56,14 @@ def searchsorted(a, v):
     3
 
     """
+    with objmode():
+        warnings.warn(
+            "`searchsorted(a, v)` is deprecated. "
+            "Use `np.searchsorted(a, v, side='right')` instead.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+
     lo = -1
     hi = len(a)
     while(lo < hi-1):
diff --git a/quantecon/util/tests/test_array.py b/quantecon/util/tests/test_array.py
index 8c015634..ad6e3368 100644
--- a/quantecon/util/tests/test_array.py
+++ b/quantecon/util/tests/test_array.py
@@ -8,9 +8,12 @@
 """
 import numpy as np
 from numpy.testing import assert_
+from numba import njit
+import pytest
 from quantecon.util import searchsorted
 
 
+@pytest.mark.filterwarnings("ignore::DeprecationWarning")
 def test_searchsorted():
     a = np.array([0.2, 0.4, 1.0])
     assert_(searchsorted(a, 0.1) == 0)
@@ -28,3 +31,18 @@ def test_searchsorted():
     a = np.ones(2)
     for (v, i) in zip([0, 1, 2], [0, 2, 2]):
         assert_(searchsorted(a, v) == i)
+
+
+@njit
+def _jitted_function():
+    a = np.array([0.2, 0.4, 1.0])
+    return searchsorted(a, 0.5)
+
+
+def test_warns():
+    a = np.array([0.2, 0.4, 1.0])
+    with pytest.warns(DeprecationWarning):
+        searchsorted(a, 0.5)
+
+    with pytest.warns(DeprecationWarning):
+        _jitted_function()