diff --git a/aeon/forecasting/stats/_arima.py b/aeon/forecasting/stats/_arima.py
index 357bc4bbe9..71c86c3bea 100644
--- a/aeon/forecasting/stats/_arima.py
+++ b/aeon/forecasting/stats/_arima.py
@@ -12,6 +12,7 @@
 from aeon.forecasting.base import BaseForecaster
 from aeon.forecasting.utils._extract_paras import _extract_arma_params
 from aeon.forecasting.utils._nelder_mead import nelder_mead
+from aeon.forecasting.utils._undifference import _undifference
 
 
 class ARIMA(BaseForecaster):
@@ -126,18 +127,11 @@ def _fit(self, y, exog=None):
         differenced_forecast = self.fitted_values_[-1]
 
         if self.d == 0:
-            forecast_value = differenced_forecast
-        elif self.d == 1:
-            forecast_value = differenced_forecast + self._series[-1]
-        else:  # for d > 1, iteratively undifference
-            forecast_value = differenced_forecast
-            last_vals = self._series[-self.d :]
-            for _ in range(self.d):
-                forecast_value += last_vals[-1] - last_vals[-2]
-                # Shift values to avoid appending to list (efficient)
-                last_vals = np.roll(last_vals, -1)
-                last_vals[-1] = forecast_value  # Extract the parameter values
-        self.forecast_ = forecast_value
+            self.forecast_ = differenced_forecast
+        else:
+            self.forecast_ = _undifference(
+                np.array([differenced_forecast]), self._series[-self.d :]
+            )[self.d]
         if self.use_constant:
             self.c_ = formatted_params[0][0]
         self.phi_ = formatted_params[1][: self.p]
@@ -198,12 +192,12 @@ def _predict(self, y, exog=None):
         forecast_diff = c + ar_forecast + ma_forecast
 
         # Undifference the forecast
-        if d == 0:
+        if self.d == 0:
             return forecast_diff
-        elif d == 1:
-            return forecast_diff + y[-1]
         else:
-            return forecast_diff + np.sum(y[-d:])
+            return _undifference(np.array([forecast_diff]), self._series[-self.d :])[
+                self.d
+            ]
 
     def _forecast(self, y, exog=None):
         """Forecast one ahead for time series y."""
@@ -242,17 +236,10 @@ def iterative_forecast(self, y, prediction_horizon):
 
         # Correct differencing using forecast values
         y_forecast_diff = forecast_series[n : n + h]
-        d = self.d
-        if d == 0:
+        if self.d == 0:
             return y_forecast_diff
-        else:  # Correct undifferencing
-            # Start with last d values from original y
-            undiff = list(self._series[-d:])
-            for i in range(h):
-                # Take the last d values and sum them
-                reconstructed = y_forecast_diff[i] + sum(undiff[-d:])
-                undiff.append(reconstructed)
-            return np.array(undiff[d:])
+        else:
+            return _undifference(y_forecast_diff, self._series[-self.d :])[self.d :]
 
 
 @njit(cache=True, fastmath=True)
diff --git a/aeon/forecasting/utils/_undifference.py b/aeon/forecasting/utils/_undifference.py
new file mode 100644
index 0000000000..289098255e
--- /dev/null
+++ b/aeon/forecasting/utils/_undifference.py
@@ -0,0 +1,64 @@
+"""Undifferencing Code."""
+
+import numpy as np
+from numba import njit
+
+# Needs to be move to DifferenceTransformer at some point
+
+
+@njit(cache=True, fastmath=True)
+def _comb(n, k):
+    """
+    Calculate the binomial coefficient C(n, k) = n! / (k! * (n - k)!).
+
+    Parameters
+    ----------
+    n : int
+        The total number of items.
+    k : int
+        The number of items to choose.
+
+    Returns
+    -------
+    int
+        The binomial coefficient C(n, k).
+    """
+    if k < 0 or k > n:
+        return 0
+    if k > n - k:
+        k = n - k  # Take advantage of symmetry
+    c = 1
+    for i in range(k):
+        c = c * (n - i) // (i + 1)
+    return c
+
+
+@njit(cache=True, fastmath=True)
+def _undifference(diff, initial_values):
+    """
+    Reconstruct original time series from an n-th order differenced series.
+
+    Parameters
+    ----------
+    diff : array-like
+        n-th order differenced series of length N - n
+    initial_values : array-like
+        The first n values of the original series before differencing (length n)
+
+    Returns
+    -------
+    original : np.ndarray
+        Reconstructed original series of length N
+    """
+    n = len(initial_values)
+    kernel = np.array(
+        [(-1) ** (k + 1) * _comb(n, k) for k in range(1, n + 1)],
+        dtype=initial_values.dtype,
+    )
+    original = np.empty((n + len(diff)), dtype=initial_values.dtype)
+    original[:n] = initial_values
+
+    for i, d in enumerate(diff):
+        original[n + i] = np.dot(kernel, original[i : n + i][::-1]) + d
+
+    return original
diff --git a/aeon/testing/testing_config.py b/aeon/testing/testing_config.py
index d199fbf59b..047491e7c6 100644
--- a/aeon/testing/testing_config.py
+++ b/aeon/testing/testing_config.py
@@ -62,6 +62,9 @@
     # Unknown issue not producing the same results
     "RDSTRegressor": ["check_regressor_against_expected_results"],
     "RISTRegressor": ["check_regressor_against_expected_results"],
+    # Requires y to be passed in invers_transform,
+    # but this is not currently enabled/supported
+    "DifferenceTransformer": ["check_transform_inverse_transform_equivalent"],
 }
 
 # Exclude specific tests for estimators here only when numba is disabled
diff --git a/aeon/transformations/series/_diff.py b/aeon/transformations/series/_diff.py
index 221987b7bd..767927041a 100644
--- a/aeon/transformations/series/_diff.py
+++ b/aeon/transformations/series/_diff.py
@@ -1,8 +1,10 @@
+"""Difference Transformer."""
+
 import numpy as np
 
 from aeon.transformations.series.base import BaseSeriesTransformer
 
-__maintainer__ = ["TinaJin0228"]
+__maintainer__ = ["TinaJin0228", "alexbanwell1"]
 __all__ = ["DifferenceTransformer"]