diff --git a/doc/source/whatsnew/v0.24.0.txt b/doc/source/whatsnew/v0.24.0.txt
index c23ed006ff637..dbe5c481160a3 100644
--- a/doc/source/whatsnew/v0.24.0.txt
+++ b/doc/source/whatsnew/v0.24.0.txt
@@ -83,7 +83,7 @@ Deprecations
 ~~~~~~~~~~~~
 
 - :meth:`DataFrame.to_stata`, :meth:`read_stata`, :class:`StataReader` and :class:`StataWriter` have deprecated the ``encoding`` argument.  The encoding of a Stata dta file is determined by the file type and cannot be changed (:issue:`21244`).
--
+- :meth:`MultiIndex.to_hierarchical` is deprecated and will be removed in a future version  (:issue:`21613`)
 -
 
 .. _whatsnew_0240.prior_deprecations:
diff --git a/pandas/core/indexes/multi.py b/pandas/core/indexes/multi.py
index ab23a80acdaae..8339e27651082 100644
--- a/pandas/core/indexes/multi.py
+++ b/pandas/core/indexes/multi.py
@@ -189,7 +189,6 @@ class MultiIndex(Index):
     from_product
     set_levels
     set_labels
-    to_hierarchical
     to_frame
     is_lexsorted
     sortlevel
@@ -1182,6 +1181,8 @@ def to_frame(self, index=True):
 
     def to_hierarchical(self, n_repeat, n_shuffle=1):
         """
+        .. deprecated:: 0.24.0
+
         Return a MultiIndex reshaped to conform to the
         shapes given by n_repeat and n_shuffle.
 
@@ -1216,6 +1217,9 @@ def to_hierarchical(self, n_repeat, n_shuffle=1):
         # Assumes that each label is divisible by n_shuffle
         labels = [x.reshape(n_shuffle, -1).ravel(order='F') for x in labels]
         names = self.names
+        warnings.warn("Method .to_hierarchical is deprecated and will "
+                      "be removed in a future version",
+                      FutureWarning, stacklevel=2)
         return MultiIndex(levels=levels, labels=labels, names=names)
 
     @property
diff --git a/pandas/core/panel.py b/pandas/core/panel.py
index c4aa471b8b944..c8797f14e1cc8 100644
--- a/pandas/core/panel.py
+++ b/pandas/core/panel.py
@@ -948,10 +948,14 @@ def to_frame(self, filter_observations=True):
             data[item] = self[item].values.ravel()[selector]
 
         def construct_multi_parts(idx, n_repeat, n_shuffle=1):
-            axis_idx = idx.to_hierarchical(n_repeat, n_shuffle)
-            labels = [x[selector] for x in axis_idx.labels]
-            levels = axis_idx.levels
-            names = axis_idx.names
+            # Replicates and shuffles MultiIndex, returns individual attributes
+            labels = [np.repeat(x, n_repeat) for x in idx.labels]
+            # Assumes that each label is divisible by n_shuffle
+            labels = [x.reshape(n_shuffle, -1).ravel(order='F')
+                      for x in labels]
+            labels = [x[selector] for x in labels]
+            levels = idx.levels
+            names = idx.names
             return labels, levels, names
 
         def construct_index_parts(idx, major=True):
diff --git a/pandas/tests/indexes/test_multi.py b/pandas/tests/indexes/test_multi.py
index ab53002ee1587..362f917e74972 100644
--- a/pandas/tests/indexes/test_multi.py
+++ b/pandas/tests/indexes/test_multi.py
@@ -1673,9 +1673,11 @@ def test_to_frame(self):
         tm.assert_frame_equal(result, expected)
 
     def test_to_hierarchical(self):
+        # GH21613
         index = MultiIndex.from_tuples([(1, 'one'), (1, 'two'), (2, 'one'), (
             2, 'two')])
-        result = index.to_hierarchical(3)
+        with tm.assert_produces_warning(FutureWarning):
+            result = index.to_hierarchical(3)
         expected = MultiIndex(levels=[[1, 2], ['one', 'two']],
                               labels=[[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
                                       [0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1]])
@@ -1683,7 +1685,8 @@ def test_to_hierarchical(self):
         assert result.names == index.names
 
         # K > 1
-        result = index.to_hierarchical(3, 2)
+        with tm.assert_produces_warning(FutureWarning):
+            result = index.to_hierarchical(3, 2)
         expected = MultiIndex(levels=[[1, 2], ['one', 'two']],
                               labels=[[0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1],
                                       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]])
@@ -1694,8 +1697,8 @@ def test_to_hierarchical(self):
         index = MultiIndex.from_tuples([(2, 'c'), (1, 'b'),
                                         (2, 'a'), (2, 'b')],
                                        names=['N1', 'N2'])
-
-        result = index.to_hierarchical(2)
+        with tm.assert_produces_warning(FutureWarning):
+            result = index.to_hierarchical(2)
         expected = MultiIndex.from_tuples([(2, 'c'), (2, 'c'), (1, 'b'),
                                            (1, 'b'),
                                            (2, 'a'), (2, 'a'),