diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md
index bbc87c016..ac15424e8 100644
--- a/Documentation/CHANGELOG.md
+++ b/Documentation/CHANGELOG.md
@@ -29,6 +29,7 @@ Release Date: TBD
 - Changes the behavior of make_lognormal_RiskyDstn so that the standard deviation represents the standard deviation of log(returns)
 - Adds detailed parameter and latex documentation to most models.
 - Add PermGroFac constructor that explicitly combines idiosyncratic and aggregate sources of growth. [1489](https://github.com/econ-ark/HARK/pull/1489)
+- Fixes typos in IdentityFunction interpolator class. [1492](https://github.com/econ-ark/HARK/pull/1492)
 
 ### 0.15.1
 
diff --git a/HARK/interpolation.py b/HARK/interpolation.py
index 62b08b0e2..ae69a8a11 100644
--- a/HARK/interpolation.py
+++ b/HARK/interpolation.py
@@ -611,9 +611,9 @@ def derivative(self, *args):
         Returns the derivative of the function with respect to the first dimension.
         """
         if self.i_dim == 0:
-            return np.ones_like(*args[0])
+            return np.ones_like(args[0])
         else:
-            return np.zeros_like(*args[0])
+            return np.zeros_like(args[0])
 
     def derivativeX(self, *args):
         """
@@ -625,9 +625,9 @@ def derivativeX(self, *args):
         else:
             j = 0
         if self.i_dim == j:
-            return np.ones_like(*args[0])
+            return np.ones_like(args[0])
         else:
-            return np.zeros_like(*args[0])
+            return np.zeros_like(args[0])
 
     def derivativeY(self, *args):
         """
@@ -639,9 +639,9 @@ def derivativeY(self, *args):
         else:
             j = 1
         if self.i_dim == j:
-            return np.ones_like(*args[0])
+            return np.ones_like(args[0])
         else:
-            return np.zeros_like(*args[0])
+            return np.zeros_like(args[0])
 
     def derivativeZ(self, *args):
         """
@@ -653,9 +653,9 @@ def derivativeZ(self, *args):
         else:
             j = 2
         if self.i_dim == j:
-            return np.ones_like(*args[0])
+            return np.ones_like(args[0])
         else:
-            return np.zeros_like(*args[0])
+            return np.zeros_like(args[0])
 
     def derivativeW(self, *args):
         """
@@ -669,9 +669,9 @@ def derivativeW(self, *args):
                 False
             ), "Derivative with respect to W can't be called when n_dims < 4!"
         if self.i_dim == j:
-            return np.ones_like(*args[0])
+            return np.ones_like(args[0])
         else:
-            return np.zeros_like(*args[0])
+            return np.zeros_like(args[0])
 
 
 class ConstantFunction(MetricObject):
@@ -4718,475 +4718,3 @@ def __call__(self, *cFuncArgs):
             )
 
         return MPC * CRRAutilityPP(c, rho=self.CRRA)
-
-
-##############################################################################
-# Examples and tests
-##############################################################################
-
-
-def main():
-    print("Sorry, HARK.interpolation doesn't actually do much on its own.")
-    print("To see some examples of its interpolation methods in action, look at any")
-    print("of the model modules in /ConsumptionSavingModel.  In the future, running")
-    print("this module will show examples of each interpolation class.")
-
-    RNG = np.random.default_rng(123)
-
-    if False:
-        x = np.linspace(1, 20, 39)
-        y = np.log(x)
-        dydx = 1.0 / x
-        f = CubicInterp(x, y, dydx)
-        x_test = np.linspace(0, 30, 200)
-        y_test = f(x_test)
-        plt.plot(x_test, y_test)
-        plt.show()
-
-    if False:
-
-        def f(x, y):
-            return 3.0 * x**2.0 + x * y + 4.0 * y**2.0
-
-        def dfdx(x, y):
-            return 6.0 * x + y
-
-        def dfdy(x, y):
-            return x + 8.0 * y
-
-        y_list = np.linspace(0, 5, 100, dtype=float)
-        xInterpolators = []
-        xInterpolators_alt = []
-        for y in y_list:
-            this_x_list = np.sort(RNG.random(100) * 5.0)
-            this_interpolation = LinearInterp(
-                this_x_list, f(this_x_list, y * np.ones(this_x_list.size))
-            )
-            that_interpolation = CubicInterp(
-                this_x_list,
-                f(this_x_list, y * np.ones(this_x_list.size)),
-                dfdx(this_x_list, y * np.ones(this_x_list.size)),
-            )
-            xInterpolators.append(this_interpolation)
-            xInterpolators_alt.append(that_interpolation)
-        g = LinearInterpOnInterp1D(xInterpolators, y_list)
-        LinearInterpOnInterp1D(xInterpolators_alt, y_list)
-
-        rand_x = RNG.random(100) * 5.0
-        rand_y = RNG.random(100) * 5.0
-        z = (f(rand_x, rand_y) - g(rand_x, rand_y)) / f(rand_x, rand_y)
-        q = (dfdx(rand_x, rand_y) - g.derivativeX(rand_x, rand_y)) / dfdx(
-            rand_x, rand_y
-        )
-        r = (dfdy(rand_x, rand_y) - g.derivativeY(rand_x, rand_y)) / dfdy(
-            rand_x, rand_y
-        )
-        # print(z)
-        # print(q)
-        # print(r)
-
-        z = (f(rand_x, rand_y) - g(rand_x, rand_y)) / f(rand_x, rand_y)
-        q = (dfdx(rand_x, rand_y) - g.derivativeX(rand_x, rand_y)) / dfdx(
-            rand_x, rand_y
-        )
-        r = (dfdy(rand_x, rand_y) - g.derivativeY(rand_x, rand_y)) / dfdy(
-            rand_x, rand_y
-        )
-        print(z)
-        # print(q)
-        # print(r)
-
-    if False:
-
-        def f(x, y, z):
-            return 3.0 * x**2.0 + x * y + 4.0 * y**2.0 - 5 * z**2.0 + 1.5 * x * z
-
-        def dfdx(x, y, z):
-            return 6.0 * x + y + 1.5 * z
-
-        def dfdy(x, y, z):
-            return x + 8.0 * y
-
-        def dfdz(x, y, z):
-            return -10.0 * z + 1.5 * x
-
-        y_list = np.linspace(0, 5, 51, dtype=float)
-        z_list = np.linspace(0, 5, 51, dtype=float)
-        xInterpolators = []
-        for y in y_list:
-            temp = []
-            for z in z_list:
-                this_x_list = np.sort(RNG.random(100) * 5.0)
-                this_interpolation = LinearInterp(
-                    this_x_list,
-                    f(
-                        this_x_list,
-                        y * np.ones(this_x_list.size),
-                        z * np.ones(this_x_list.size),
-                    ),
-                )
-                temp.append(this_interpolation)
-            xInterpolators.append(deepcopy(temp))
-        g = BilinearInterpOnInterp1D(xInterpolators, y_list, z_list)
-
-        rand_x = RNG.random(1000) * 5.0
-        rand_y = RNG.random(1000) * 5.0
-        rand_z = RNG.random(1000) * 5.0
-        z = (f(rand_x, rand_y, rand_z) - g(rand_x, rand_y, rand_z)) / f(
-            rand_x, rand_y, rand_z
-        )
-        q = (
-            dfdx(rand_x, rand_y, rand_z) - g.derivativeX(rand_x, rand_y, rand_z)
-        ) / dfdx(rand_x, rand_y, rand_z)
-        r = (
-            dfdy(rand_x, rand_y, rand_z) - g.derivativeY(rand_x, rand_y, rand_z)
-        ) / dfdy(rand_x, rand_y, rand_z)
-        p = (
-            dfdz(rand_x, rand_y, rand_z) - g.derivativeZ(rand_x, rand_y, rand_z)
-        ) / dfdz(rand_x, rand_y, rand_z)
-        z.sort()
-
-    if False:
-
-        def f(w, x, y, z):
-            return (
-                4.0 * w * z
-                - 2.5 * w * x
-                + w * y
-                + 6.0 * x * y
-                - 10.0 * x * z
-                + 3.0 * y * z
-                - 7.0 * z
-                + 4.0 * x
-                + 2.0 * y
-                - 5.0 * w
-            )
-
-        def dfdw(w, x, y, z):
-            return 4.0 * z - 2.5 * x + y - 5.0
-
-        def dfdx(w, x, y, z):
-            return -2.5 * w + 6.0 * y - 10.0 * z + 4.0
-
-        def dfdy(w, x, y, z):
-            return w + 6.0 * x + 3.0 * z + 2.0
-
-        def dfdz(w, x, y, z):
-            return 4.0 * w - 10.0 * x + 3.0 * y - 7
-
-        x_list = np.linspace(0, 5, 16, dtype=float)
-        y_list = np.linspace(0, 5, 16, dtype=float)
-        z_list = np.linspace(0, 5, 16, dtype=float)
-        wInterpolators = []
-        for x in x_list:
-            temp = []
-            for y in y_list:
-                temptemp = []
-                for z in z_list:
-                    this_w_list = np.sort(RNG.random(16) * 5.0)
-                    this_interpolation = LinearInterp(
-                        this_w_list,
-                        f(
-                            this_w_list,
-                            x * np.ones(this_w_list.size),
-                            y * np.ones(this_w_list.size),
-                            z * np.ones(this_w_list.size),
-                        ),
-                    )
-                    temptemp.append(this_interpolation)
-                temp.append(deepcopy(temptemp))
-            wInterpolators.append(deepcopy(temp))
-        g = TrilinearInterpOnInterp1D(wInterpolators, x_list, y_list, z_list)
-
-        N = 20000
-        rand_w = RNG.random(N) * 5.0
-        rand_x = RNG.random(N) * 5.0
-        rand_y = RNG.random(N) * 5.0
-        rand_z = RNG.random(N) * 5.0
-        t_start = time()
-        z = (f(rand_w, rand_x, rand_y, rand_z) - g(rand_w, rand_x, rand_y, rand_z)) / f(
-            rand_w, rand_x, rand_y, rand_z
-        )
-        q = (
-            dfdw(rand_w, rand_x, rand_y, rand_z)
-            - g.derivativeW(rand_w, rand_x, rand_y, rand_z)
-        ) / dfdw(rand_w, rand_x, rand_y, rand_z)
-        r = (
-            dfdx(rand_w, rand_x, rand_y, rand_z)
-            - g.derivativeX(rand_w, rand_x, rand_y, rand_z)
-        ) / dfdx(rand_w, rand_x, rand_y, rand_z)
-        p = (
-            dfdy(rand_w, rand_x, rand_y, rand_z)
-            - g.derivativeY(rand_w, rand_x, rand_y, rand_z)
-        ) / dfdy(rand_w, rand_x, rand_y, rand_z)
-        (
-            dfdz(rand_w, rand_x, rand_y, rand_z)
-            - g.derivativeZ(rand_w, rand_x, rand_y, rand_z)
-        ) / dfdz(rand_w, rand_x, rand_y, rand_z)
-        t_end = time()
-
-        z.sort()
-        print(z)
-        print(t_end - t_start)
-
-    if False:
-
-        def f(x, y):
-            return 3.0 * x**2.0 + x * y + 4.0 * y**2.0
-
-        def dfdx(x, y):
-            return 6.0 * x + y
-
-        def dfdy(x, y):
-            return x + 8.0 * y
-
-        x_list = np.linspace(0, 5, 101, dtype=float)
-        y_list = np.linspace(0, 5, 101, dtype=float)
-        x_temp, y_temp = np.meshgrid(x_list, y_list, indexing="ij")
-        g = BilinearInterp(f(x_temp, y_temp), x_list, y_list)
-
-        rand_x = RNG.random(100) * 5.0
-        rand_y = RNG.random(100) * 5.0
-        z = (f(rand_x, rand_y) - g(rand_x, rand_y)) / f(rand_x, rand_y)
-        q = (f(x_temp, y_temp) - g(x_temp, y_temp)) / f(x_temp, y_temp)
-        # print(z)
-        # print(q)
-
-    if False:
-
-        def f(x, y, z):
-            return 3.0 * x**2.0 + x * y + 4.0 * y**2.0 - 5 * z**2.0 + 1.5 * x * z
-
-        def dfdx(x, y, z):
-            return 6.0 * x + y + 1.5 * z
-
-        def dfdy(x, y, z):
-            return x + 8.0 * y
-
-        def dfdz(x, y, z):
-            return -10.0 * z + 1.5 * x
-
-        x_list = np.linspace(0, 5, 11, dtype=float)
-        y_list = np.linspace(0, 5, 11, dtype=float)
-        z_list = np.linspace(0, 5, 101, dtype=float)
-        x_temp, y_temp, z_temp = np.meshgrid(x_list, y_list, z_list, indexing="ij")
-        g = TrilinearInterp(f(x_temp, y_temp, z_temp), x_list, y_list, z_list)
-
-        rand_x = RNG.random(1000) * 5.0
-        rand_y = RNG.random(1000) * 5.0
-        rand_z = RNG.random(1000) * 5.0
-        z = (f(rand_x, rand_y, rand_z) - g(rand_x, rand_y, rand_z)) / f(
-            rand_x, rand_y, rand_z
-        )
-        q = (
-            dfdx(rand_x, rand_y, rand_z) - g.derivativeX(rand_x, rand_y, rand_z)
-        ) / dfdx(rand_x, rand_y, rand_z)
-        r = (
-            dfdy(rand_x, rand_y, rand_z) - g.derivativeY(rand_x, rand_y, rand_z)
-        ) / dfdy(rand_x, rand_y, rand_z)
-        p = (
-            dfdz(rand_x, rand_y, rand_z) - g.derivativeZ(rand_x, rand_y, rand_z)
-        ) / dfdz(rand_x, rand_y, rand_z)
-        p.sort()
-        plt.plot(p)
-
-    if False:
-
-        def f(w, x, y, z):
-            return (
-                4.0 * w * z
-                - 2.5 * w * x
-                + w * y
-                + 6.0 * x * y
-                - 10.0 * x * z
-                + 3.0 * y * z
-                - 7.0 * z
-                + 4.0 * x
-                + 2.0 * y
-                - 5.0 * w
-            )
-
-        def dfdw(w, x, y, z):
-            return 4.0 * z - 2.5 * x + y - 5.0
-
-        def dfdx(w, x, y, z):
-            return -2.5 * w + 6.0 * y - 10.0 * z + 4.0
-
-        def dfdy(w, x, y, z):
-            return w + 6.0 * x + 3.0 * z + 2.0
-
-        def dfdz(w, x, y, z):
-            return 4.0 * w - 10.0 * x + 3.0 * y - 7
-
-        w_list = np.linspace(0, 5, 16, dtype=float)
-        x_list = np.linspace(0, 5, 16, dtype=float)
-        y_list = np.linspace(0, 5, 16, dtype=float)
-        z_list = np.linspace(0, 5, 16, dtype=float)
-        w_temp, x_temp, y_temp, z_temp = np.meshgrid(
-            w_list, x_list, y_list, z_list, indexing="ij"
-        )
-
-        def mySearch(trash, x):
-            return np.floor(x / 5 * 32).astype(int)
-
-        g = QuadlinearInterp(
-            f(w_temp, x_temp, y_temp, z_temp), w_list, x_list, y_list, z_list
-        )
-
-        N = 1000000
-        rand_w = RNG.random(N) * 5.0
-        rand_x = RNG.random(N) * 5.0
-        rand_y = RNG.random(N) * 5.0
-        rand_z = RNG.random(N) * 5.0
-        t_start = time()
-        z = (f(rand_w, rand_x, rand_y, rand_z) - g(rand_w, rand_x, rand_y, rand_z)) / f(
-            rand_w, rand_x, rand_y, rand_z
-        )
-        t_end = time()
-        # print(z)
-        print(t_end - t_start)
-
-    if False:
-
-        def f(x, y):
-            return 3.0 * x**2.0 + x * y + 4.0 * y**2.0
-
-        def dfdx(x, y):
-            return 6.0 * x + y
-
-        def dfdy(x, y):
-            return x + 8.0 * y
-
-        warp_factor = 0.01
-        x_list = np.linspace(0, 5, 71, dtype=float)
-        y_list = np.linspace(0, 5, 51, dtype=float)
-        x_temp, y_temp = np.meshgrid(x_list, y_list, indexing="ij")
-        x_adj = x_temp + warp_factor * (RNG.random(x_list.size, y_list.size) - 0.5)
-        y_adj = y_temp + warp_factor * (RNG.random(x_list.size, y_list.size) - 0.5)
-        g = Curvilinear2DInterp(f(x_adj, y_adj), x_adj, y_adj)
-
-        rand_x = RNG.random(1000) * 5.0
-        rand_y = RNG.random(1000) * 5.0
-        t_start = time()
-        z = (f(rand_x, rand_y) - g(rand_x, rand_y)) / f(rand_x, rand_y)
-        q = (dfdx(rand_x, rand_y) - g.derivativeX(rand_x, rand_y)) / dfdx(
-            rand_x, rand_y
-        )
-        r = (dfdy(rand_x, rand_y) - g.derivativeY(rand_x, rand_y)) / dfdy(
-            rand_x, rand_y
-        )
-        t_end = time()
-        z.sort()
-        q.sort()
-        r.sort()
-        # print(z)
-        print(t_end - t_start)
-
-    if False:
-
-        def f(x, y, z):
-            return 3.0 * x**2.0 + x * y + 4.0 * y**2.0 - 5 * z**2.0 + 1.5 * x * z
-
-        def dfdx(x, y, z):
-            return 6.0 * x + y + 1.5 * z
-
-        def dfdy(x, y, z):
-            return x + 8.0 * y
-
-        def dfdz(x, y, z):
-            return -10.0 * z + 1.5 * x
-
-        warp_factor = 0.01
-        x_list = np.linspace(0, 5, 11, dtype=float)
-        y_list = np.linspace(0, 5, 11, dtype=float)
-        z_list = np.linspace(0, 5, 101, dtype=float)
-        x_temp, y_temp = np.meshgrid(x_list, y_list, indexing="ij")
-        xyInterpolators = []
-        for j in range(z_list.size):
-            x_adj = x_temp + warp_factor * (RNG.random(x_list.size, y_list.size) - 0.5)
-            y_adj = y_temp + warp_factor * (RNG.random(x_list.size, y_list.size) - 0.5)
-            z_temp = z_list[j] * np.ones(x_adj.shape)
-            thisInterp = Curvilinear2DInterp(f(x_adj, y_adj, z_temp), x_adj, y_adj)
-            xyInterpolators.append(thisInterp)
-        g = LinearInterpOnInterp2D(xyInterpolators, z_list)
-
-        N = 1000
-        rand_x = RNG.random(N) * 5.0
-        rand_y = RNG.random(N) * 5.0
-        rand_z = RNG.random(N) * 5.0
-        z = (f(rand_x, rand_y, rand_z) - g(rand_x, rand_y, rand_z)) / f(
-            rand_x, rand_y, rand_z
-        )
-        p = (
-            dfdz(rand_x, rand_y, rand_z) - g.derivativeZ(rand_x, rand_y, rand_z)
-        ) / dfdz(rand_x, rand_y, rand_z)
-        p.sort()
-        plt.plot(p)
-
-    if False:
-
-        def f(w, x, y, z):
-            return (
-                4.0 * w * z
-                - 2.5 * w * x
-                + w * y
-                + 6.0 * x * y
-                - 10.0 * x * z
-                + 3.0 * y * z
-                - 7.0 * z
-                + 4.0 * x
-                + 2.0 * y
-                - 5.0 * w
-            )
-
-        def dfdw(w, x, y, z):
-            return 4.0 * z - 2.5 * x + y - 5.0
-
-        def dfdx(w, x, y, z):
-            return -2.5 * w + 6.0 * y - 10.0 * z + 4.0
-
-        def dfdy(w, x, y, z):
-            return w + 6.0 * x + 3.0 * z + 2.0
-
-        def dfdz(w, x, y, z):
-            return 4.0 * w - 10.0 * x + 3.0 * y - 7
-
-        warp_factor = 0.1
-        w_list = np.linspace(0, 5, 16, dtype=float)
-        x_list = np.linspace(0, 5, 16, dtype=float)
-        y_list = np.linspace(0, 5, 16, dtype=float)
-        z_list = np.linspace(0, 5, 16, dtype=float)
-        w_temp, x_temp = np.meshgrid(w_list, x_list, indexing="ij")
-        wxInterpolators = []
-        for i in range(y_list.size):
-            temp = []
-            for j in range(z_list.size):
-                w_adj = w_temp + warp_factor * (
-                    RNG.random(w_list.size, x_list.size) - 0.5
-                )
-                x_adj = x_temp + warp_factor * (
-                    RNG.random(w_list.size, x_list.size) - 0.5
-                )
-                y_temp = y_list[i] * np.ones(w_adj.shape)
-                z_temp = z_list[j] * np.ones(w_adj.shape)
-                thisInterp = Curvilinear2DInterp(
-                    f(w_adj, x_adj, y_temp, z_temp), w_adj, x_adj
-                )
-                temp.append(thisInterp)
-            wxInterpolators.append(temp)
-        g = BilinearInterpOnInterp2D(wxInterpolators, y_list, z_list)
-
-        N = 1000000
-        rand_w = RNG.random(N) * 5.0
-        rand_x = RNG.random(N) * 5.0
-        rand_y = RNG.random(N) * 5.0
-        rand_z = RNG.random(N) * 5.0
-
-        t_start = time()
-        z = (f(rand_w, rand_x, rand_y, rand_z) - g(rand_w, rand_x, rand_y, rand_z)) / f(
-            rand_w, rand_x, rand_y, rand_z
-        )
-        t_end = time()
-        z.sort()
-        print(z)
-        print(t_end - t_start)
diff --git a/HARK/tests/test_interpolation.py b/HARK/tests/test_interpolation.py
index 200b49d0e..66df00b41 100644
--- a/HARK/tests/test_interpolation.py
+++ b/HARK/tests/test_interpolation.py
@@ -9,6 +9,7 @@
 from HARK.interpolation import BilinearInterp
 from HARK.interpolation import CubicHermiteInterp as CubicInterp
 from HARK.interpolation import LinearInterp, QuadlinearInterp, TrilinearInterp
+from HARK.interpolation import IdentityFunction
 
 
 class testsLinearInterp(unittest.TestCase):
@@ -200,3 +201,45 @@ def test_same_length(self):
             self.f_array, self.w_array, self.x_array, self.y_array_t, self.z_array
         )
         self.assertEqual(bilinear(1, 2, 1, 2), 6.0)
+
+
+class test_IdentityFunction(unittest.TestCase):
+    """
+    Tests evaluation and derivatives of IdentityFunction class.
+    """
+
+    def setUp(self):
+        self.IF1D = IdentityFunction()
+        self.IF2Da = IdentityFunction(i_dim=0, n_dims=2)
+        self.IF2Db = IdentityFunction(i_dim=1, n_dims=2)
+        self.IF3Da = IdentityFunction(i_dim=0, n_dims=3)
+        self.IF3Db = IdentityFunction(i_dim=2, n_dims=3)
+        self.X = 3 * np.ones(100)
+        self.Y = 4 * np.ones(100)
+        self.Z = 5 * np.ones(100)
+        self.zero = np.zeros(100)
+        self.one = np.ones(100)
+
+    def test_eval(self):
+        assert np.all(self.X == self.IF1D(self.X))
+        assert np.all(self.X == self.IF2Da(self.X, self.Y))
+        assert np.all(self.Y == self.IF2Db(self.X, self.Y))
+        assert np.all(self.X == self.IF3Da(self.X, self.Y, self.Z))
+        assert np.all(self.Z == self.IF3Db(self.X, self.Y, self.Z))
+
+    def test_der(self):
+        assert np.all(self.one == self.IF1D.derivative(self.X))
+
+        assert np.all(self.one == self.IF2Da.derivativeX(self.X, self.Y))
+        assert np.all(self.zero == self.IF2Da.derivativeY(self.X, self.Y))
+
+        assert np.all(self.zero == self.IF2Db.derivativeX(self.X, self.Y))
+        assert np.all(self.one == self.IF2Db.derivativeY(self.X, self.Y))
+
+        assert np.all(self.one == self.IF3Da.derivativeX(self.X, self.Y, self.Z))
+        assert np.all(self.zero == self.IF3Da.derivativeY(self.X, self.Y, self.Z))
+        assert np.all(self.zero == self.IF3Da.derivativeZ(self.X, self.Y, self.Z))
+
+        assert np.all(self.zero == self.IF3Db.derivativeX(self.X, self.Y, self.Z))
+        assert np.all(self.zero == self.IF3Db.derivativeY(self.X, self.Y, self.Z))
+        assert np.all(self.one == self.IF3Db.derivativeZ(self.X, self.Y, self.Z))