diff --git a/aeon/forecasting/__init__.py b/aeon/forecasting/__init__.py
index c6ca4299ca..45eaa41715 100644
--- a/aeon/forecasting/__init__.py
+++ b/aeon/forecasting/__init__.py
@@ -1,8 +1,14 @@
 """Forecasters."""
 
-__all__ = ["DummyForecaster", "BaseForecaster", "RegressionForecaster", "ETSForecaster"]
+__all__ = [
+    "DummyForecaster",
+    "BaseForecaster",
+    "RegressionForecaster",
+    "ETSForecaster",
+    "ModelType",
+]
 
 from aeon.forecasting._dummy import DummyForecaster
-from aeon.forecasting._ets import ETSForecaster
+from aeon.forecasting._ets import ETSForecaster, ModelType
 from aeon.forecasting._regression import RegressionForecaster
 from aeon.forecasting.base import BaseForecaster
diff --git a/aeon/forecasting/_ets.py b/aeon/forecasting/_ets.py
index 0b4172ca3c..7c97378002 100644
--- a/aeon/forecasting/_ets.py
+++ b/aeon/forecasting/_ets.py
@@ -1,35 +1,129 @@
+"""ETSForecaster class.
+
+An implementation of the exponential smoothing statistics forecasting algorithm.
+Implements additive and multiplicative error models,
+None, additive and multiplicative (including damped) trend and
+None, additive and mutliplicative seasonality
+
+aeon enhancement proposal
+https://github.com/aeon-toolkit/aeon/pull/2244/
+
+"""
+
+__maintainer__ = []
+__all__ = ["ETSForecaster", "ModelType"]
+
 import numpy as np
 
 from aeon.forecasting.base import BaseForecaster
 
+NONE = 0
+ADDITIVE = 1
+MULTIPLICATIVE = 2
+
+
+class ModelType:
+    """
+    Class describing the error, trend and seasonality model of an ETS forecaster.
+
+    Attributes
+    ----------
+    error_type : int
+        The type of error model; either Additive(1) or Multiplicative(2)
+    trend_type : int
+        The type of trend model; one of None(0), additive(1) or multiplicative(2).
+    seasonality_type : int
+        The type of seasonality model; one of None(0), additive(1) or multiplicative(2).
+    seasonal_period : int
+        The period of the seasonality (m) (e.g., for quaterly data seasonal_period = 4).
+    """
+
+    error_type: int
+    trend_type: int
+    seasonality_type: int
+    seasonal_period: int
+
+    def __init__(
+        self,
+        error_type=ADDITIVE,
+        trend_type=NONE,
+        seasonality_type=NONE,
+        seasonal_period=1,
+    ):
+        assert error_type != NONE, "Error must be either additive or multiplicative"
+        if seasonal_period < 1 or seasonality_type == NONE:
+            seasonal_period = 1
+        self.error_type = error_type
+        self.trend_type = trend_type
+        self.seasonality_type = seasonality_type
+        self.seasonal_period = seasonal_period
+
 
 class ETSForecaster(BaseForecaster):
     """Exponential Smoothing forecaster.
 
-    Simple first implementation with Holt-Winters method
-    and no seasonality.
+    An implementation of the exponential smoothing statistics forecasting algorithm.
+    Implements additive and multiplicative error models,
+    None, additive and multiplicative (including damped) trend and
+    None, additive and mutliplicative seasonality[1]_.
 
     Parameters
     ----------
-    alpha : float, default = 0.2
+    alpha : float, default = 0.1
         Level smoothing parameter.
-    beta : float, default = 0.2
+    beta : float, default = 0.01
         Trend smoothing parameter.
-    gamma : float, default = 0.2
+    gamma : float, default = 0.01
         Seasonal smoothing parameter.
-    season_len : int, default = 1
-        The length of the seasonality period.
+    phi : float, default = 0.99
+        Trend damping smoothing parameters
+    horizon : int, default = 1
+        The horizon to forecast to.
+    model_type : ModelType, default = ModelType()
+        A object of type ModelType, describing the error,
+        trend and seasonality type of this ETS model.
+
+    References
+    ----------
+    .. [1] R. J. Hyndman and G. Athanasopoulos,
+        Forecasting: Principles and Practice. Melbourne, Australia: OTexts, 2014.
+
+    Examples
+    --------
+    >>> from aeon.forecasting import ETSForecaster, ModelType
+    >>> from aeon.datasets import load_airline
+    >>> y = load_airline()
+    >>> forecaster = ETSForecaster(alpha=0.4, beta=0.2, gamma=0.5, phi=0.8, horizon=1,
+                               model_type=ModelType(1,2,2,4))
+    >>> forecaster.fit(y)
+    >>> forecaster.predict()
+    366.90200486015596
     """
 
-    def __init__(self, alpha=0.2, beta=0.2, gamma=0.2, season_len=1, horizon=1):
+    default_model_type = ModelType()
+
+    def __init__(
+        self,
+        alpha=0.1,
+        beta=0.01,
+        gamma=0.01,
+        phi=0.99,
+        horizon=1,
+        model_type=default_model_type,
+    ):
         self.alpha = alpha
         self.beta = beta
         self.gamma = gamma
-        self.season_len = season_len
+        self.phi = phi
         self.forecast_val_ = 0.0
         self.level_ = 0.0
         self.trend_ = 0.0
         self.season_ = None
+        self.n_timepoints = 0
+        self.avg_mean_sq_err_ = 0
+        self.liklihood_ = 0
+        self.residuals_ = []
+        self.model_type = model_type
         super().__init__(horizon=horizon, axis=1)
 
     def _fit(self, y, exog=None):
@@ -51,31 +145,211 @@ def _fit(self, y, exog=None):
         """
         data = y.squeeze()
         self.n_timepoints = len(data)
-        sl = self.season_len
-        # Initialize components
-        self.level_ = data[0]
-        self.trend_ = np.mean(data[sl : 2 * sl]) - np.mean(data[:sl])
-        self.season_ = [data[i] / data[0] for i in range(sl)]
-        for t in range(sl, self.n_timepoints):
+        self._initialise(data)
+        self.avg_mean_sq_err_ = 0
+        self.liklihood_ = 0
+        mul_liklihood_pt2 = 0
+        self.residuals_ = np.zeros(
+            self.n_timepoints
+        )  # 1 Less residual than data points
+        for t, data_item in enumerate(data[self.model_type.seasonal_period :]):
             # Calculate level, trend, and seasonal components
-            level_prev = self.level_
-            l1 = data[t] / self.season_[t % self.season_len]
-            l2 = self.level_ + self.trend_
-            self.level_ = self.alpha * l1 + (1 - self.alpha) * l2
-            trend = self.level_ - level_prev
-            self.trend_ = self.beta * trend + (1 - self.beta) * self.trend_
-            s1 = data[t] / self.level_
-            s2 = self.season_[t % sl]
-            self.season_[t % self.season_len] = self.gamma * s1 + (1 - self.gamma) * s2
+            fitted_value, error = self._update_states(
+                data_item, t % self.model_type.seasonal_period
+            )
+            self.residuals_[t] = error
+            self.avg_mean_sq_err_ += (data_item - fitted_value) ** 2
+            self.liklihood_ += error * error
+            mul_liklihood_pt2 += np.log(np.fabs(fitted_value))
+        self.avg_mean_sq_err_ /= self.n_timepoints - self.model_type.seasonal_period
+        self.liklihood_ = (
+            self.n_timepoints - self.model_type.seasonal_period
+        ) * np.log(self.liklihood_)
+        if self.model_type.error_type == MULTIPLICATIVE:
+            self.liklihood_ += 2 * mul_liklihood_pt2
         return self
 
+    def _update_states(self, data_item, seasonal_index):
+        """
+        Update level, trend, and seasonality components.
+
+        Using state space equations for an ETS model.
+
+        Parameters
+        ----------
+        data_item: float
+            The current value of the time series.
+        seasonal_index: int
+            The index to update the seasonal component.
+        """
+        model = self.model_type
+        # Retrieve the current state values
+        level = self.level_
+        trend = self.trend_
+        seasonality = self.season_[seasonal_index]
+        fitted_value, damped_trend, trend_level_combination = self._predict_value(
+            trend, level, seasonality, self.phi
+        )
+        # Calculate the error term (observed value - fitted value)
+        if model.error_type == MULTIPLICATIVE:
+            error = data_item / fitted_value - 1  # Multiplicative error
+        else:
+            error = data_item - fitted_value  # Additive error
+        # Update level
+        if model.error_type == MULTIPLICATIVE:
+            self.level_ = trend_level_combination * (1 + self.alpha * error)
+            self.trend_ = damped_trend * (1 + self.beta * error)
+            self.season_[seasonal_index] = seasonality * (1 + self.gamma * error)
+            if model.seasonality_type == ADDITIVE:
+                self.level_ += (
+                    self.alpha * error * seasonality
+                )  # Add seasonality correction
+                self.season_[seasonal_index] += (
+                    self.gamma * error * trend_level_combination
+                )
+                if model.trend_type == ADDITIVE:
+                    self.trend_ += (level + seasonality) * self.beta * error
+                else:
+                    self.trend_ += seasonality / level * self.beta * error
+            elif model.trend_type == ADDITIVE:
+                self.trend_ += level * self.beta * error
+        else:
+            level_correction = 1
+            trend_correction = 1
+            seasonality_correction = 1
+            if model.seasonality_type == MULTIPLICATIVE:
+                # Add seasonality correction
+                level_correction *= seasonality
+                trend_correction *= seasonality
+                seasonality_correction *= trend_level_combination
+            if model.trend_type == MULTIPLICATIVE:
+                trend_correction *= level
+            self.level_ = (
+                trend_level_combination + self.alpha * error / level_correction
+            )
+            self.trend_ = damped_trend + self.beta * error / trend_correction
+            self.season_[seasonal_index] = (
+                seasonality + self.gamma * error / seasonality_correction
+            )
+        return (fitted_value, error)
+
+    def _initialise(self, data):
+        """
+        Initialize level, trend, and seasonality values for the ETS model.
+
+        Parameters
+        ----------
+        data : array-like
+            The time series data
+            (should contain at least two full seasons if seasonality is specified)
+        """
+        model = self.model_type
+        # Initial Level: Mean of the first season
+        self.level_ = np.mean(data[: model.seasonal_period])
+        # Initial Trend
+        if model.trend_type == ADDITIVE:
+            # Average difference between corresponding points in the first two seasons
+            self.trend_ = np.mean(
+                data[model.seasonal_period : 2 * model.seasonal_period]
+                - data[: model.seasonal_period]
+            )
+        elif model.trend_type == MULTIPLICATIVE:
+            # Average ratio between corresponding points in the first two seasons
+            self.trend_ = np.mean(
+                data[model.seasonal_period : 2 * model.seasonal_period]
+                / data[: model.seasonal_period]
+            )
+        else:
+            # No trend
+            self.trend_ = 0
+            self.beta = (
+                0  # Required for the equations in _update_states to work correctly
+            )
+        # Initial Seasonality
+        if model.seasonality_type == ADDITIVE:
+            # Seasonal component is the difference
+            # from the initial level for each point in the first season
+            self.season_ = data[: model.seasonal_period] - self.level_
+        elif model.seasonality_type == MULTIPLICATIVE:
+            # Seasonal component is the ratio of each point in the first season
+            # to the initial level
+            self.season_ = data[: model.seasonal_period] / self.level_
+        else:
+            # No seasonality
+            self.season_ = [0]
+            self.gamma = (
+                0  # Required for the equations in _update_states to work correctly
+            )
+
     def _predict(self, y=None, exog=None):
+        """
+        Predict the next horizon steps ahead.
+
+        Parameters
+        ----------
+        y : np.ndarray, default = None
+            A time series to predict the next horizon value for. If None,
+            predict the next horizon value after series seen in fit.
+        exog : np.ndarray, default =None
+            Optional exogenous time series data assumed to be aligned with y
+
+        Returns
+        -------
+        float
+            single prediction self.horizon steps ahead of y.
+        """
         # Generate forecasts based on the final values of level, trend, and seasonals
-        trend = (self.horizon + 1) * self.trend_
-        seasonal = self.season_[(self.n_timepoints + self.horizon) % self.season_len]
-        forecast = (self.level_ + trend) * seasonal
-        return forecast
-
-    def _forecast(self, y):
-        self.fit(y)
-        return self.predict()
+        if self.phi == 1:  # No damping case
+            phi_h = float(self.horizon)
+        else:
+            # Geometric series formula for calculating phi + phi^2 + ... + phi^h
+            phi_h = self.phi * (1 - self.phi**self.horizon) / (1 - self.phi)
+        seasonality = self.season_[
+            (self.n_timepoints + self.horizon) % self.model_type.seasonal_period
+        ]
+        fitted_value = self._predict_value(
+            self.trend_, self.level_, seasonality, phi_h
+        )[0]
+        return fitted_value
+
+    def _predict_value(self, trend, level, seasonality, phi):
+        """
+
+        Generate various useful values, including the next fitted value.
+
+        Parameters
+        ----------
+        trend : float
+            The current trend value for the model
+        level : float
+            The current level value for the model
+        seasonality : float
+            The current seasonality value for the model
+        phi : float
+            The damping parameter for the model
+
+        Returns
+        -------
+        fitted_value : float
+            single prediction based on the current state variables.
+        damped_trend : float
+            The damping parameter combined with the trend dependant on the model type
+        trend_level_combination : float
+            Combination of the trend and level based on the model type.
+        """
+        model = self.model_type
+        # Apply damping parameter and
+        # calculate commonly used combination of trend and level components
+        if model.trend_type == MULTIPLICATIVE:
+            damped_trend = trend**phi
+            trend_level_combination = level * damped_trend
+        else:  # Additive trend, if no trend, then trend = 0
+            damped_trend = trend * phi
+            trend_level_combination = level + damped_trend
+
+        # Calculate forecast (fitted value) based on the current components
+        if model.seasonality_type == MULTIPLICATIVE:
+            fitted_value = trend_level_combination * seasonality
+        else:  # Additive seasonality, if no seasonality, then seasonality = 0
+            fitted_value = trend_level_combination + seasonality
+        return fitted_value, damped_trend, trend_level_combination
diff --git a/aeon/forecasting/_verify_ets.py b/aeon/forecasting/_verify_ets.py
new file mode 100644
index 0000000000..d93b64631f
--- /dev/null
+++ b/aeon/forecasting/_verify_ets.py
@@ -0,0 +1,168 @@
+import random
+import time
+
+import numpy as np
+from statsforecast.ets import etscalc
+from statsforecast.utils import AirPassengers as ap
+
+from aeon.forecasting import ETSForecaster, ModelType
+
+NA = -99999.0
+MAX_NMSE = 30
+MAX_SEASONAL_PERIOD = 24
+
+
+def setup():
+    """Generate parameters required for ETS algorithms."""
+    y = ap
+    n = len(ap)
+    m = random.randint(1, 24)
+    error = random.randint(1, 2)
+    trend = random.randint(0, 2)
+    season = random.randint(0, 2)
+    alpha = round(random.random(), 4)
+    if alpha == 0:
+        alpha = round(random.random(), 4)
+    beta = round(random.random() * alpha, 4)  # 0 < beta < alpha
+    if beta == 0:
+        beta = round(random.random() * alpha, 4)
+    gamma = round(random.random() * (1 - alpha), 4)  # 0 < beta < alpha
+    if gamma == 0:
+        gamma = round(random.random() * (1 - alpha), 4)
+    phi = round(
+        random.random() * 0.18 + 0.8, 4
+    )  # Common constraint for phi is 0.8 < phi < 0.98
+    e = np.zeros(n)
+    lik_fitets = np.zeros(1)
+    amse = np.zeros(MAX_NMSE)
+    nmse = 3
+    return (
+        y,
+        n,
+        m,
+        error,
+        trend,
+        season,
+        alpha,
+        beta,
+        gamma,
+        phi,
+        e,
+        lik_fitets,
+        amse,
+        nmse,
+    )
+
+
+def test_ets_comparison(setup_func, random_seed, catch_errors):
+    """Run both our statsforecast and our implementation and crosschecks."""
+    random.seed(random_seed)
+    (
+        y,
+        n,
+        m,
+        error,
+        trend,
+        season,
+        alpha,
+        beta,
+        gamma,
+        phi,
+        e,
+        lik_fitets,
+        amse,
+        nmse,
+    ) = setup_func()
+    # tsml-eval implementation
+    start = time.time()
+    f1 = ETSForecaster(alpha, beta, gamma, phi, 1, ModelType(error, trend, season, m))
+    f1.fit(y)
+    end = time.time()
+    time_fitets = end - start
+    e_fitets = f1.residuals_
+    amse_fitets = f1.avg_mean_sq_err_
+    lik_fitets = f1.liklihood_
+    # Reinitialise arrays
+    e.fill(0)
+    amse.fill(0)
+    f1 = ETSForecaster(alpha, beta, gamma, phi, 1, ModelType(error, trend, season, m))
+    f1._initialise(y)
+    init_states_etscalc = np.zeros(n * (1 + (trend > 0) + m * (season > 0) + 1))
+    init_states_etscalc[0] = f1.level_
+    init_states_etscalc[1] = f1.trend_
+    init_states_etscalc[1 + (trend != 0) : m + 1 + (trend != 0)] = f1.season_[::-1]
+    if season == 0:
+        m = 1
+    # Nixtla/statsforcast implementation
+    start = time.time()
+    lik_etscalc = etscalc(
+        y[m:],
+        n - m,
+        init_states_etscalc,
+        m,
+        error,
+        trend,
+        season,
+        alpha,
+        beta,
+        gamma,
+        phi,
+        e,
+        amse,
+        nmse,
+    )
+    end = time.time()
+    time_etscalc = end - start
+    e_etscalc = e.copy()
+    amse_etscalc = amse.copy()[0]
+
+    if catch_errors:
+        try:
+            # Comparing outputs and runtime
+            assert np.allclose(e_fitets, e_etscalc), "Residuals Compare failed"
+            assert np.allclose(amse_fitets, amse_etscalc), "AMSE Compare failed"
+            assert np.isclose(lik_fitets, lik_etscalc), "Liklihood Compare failed"
+            return True
+        except AssertionError as e:
+            print(e)  # noqa
+            print(  # noqa
+                f"Seed: {random_seed}, Model: Error={error}, Trend={trend},\
+                   Seasonality={season}, seasonal period={m},\
+                   alpha={alpha}, beta={beta}, gamma={gamma}, phi={phi}"
+            )
+            return False
+    else:
+        print(  # noqa
+            f"Seed: {random_seed}, Model: Error={error}, Trend={trend},\
+                Seasonality={season}, seasonal period={m}, alpha={alpha},\
+                    beta={beta}, gamma={gamma}, phi={phi}"
+        )
+        diff_indices = np.where(
+            np.abs(e_fitets - e_etscalc) > 1e-5 * np.abs(e_etscalc) + 1e-8
+        )[0]
+        for index in diff_indices:
+            print(  # noqa
+                f"Index {index}: e_fitets = {e_fitets[index]},\
+                e_etscalc = {e_etscalc[index]}"
+            )
+        print(amse_fitets)  # noqa
+        print(amse_etscalc)  # noqa
+        print(lik_fitets)  # noqa
+        print(lik_etscalc)  # noqa
+        # assert np.allclose(init_states_fitets, init_states_etscalc)
+        assert np.allclose(e_fitets, e_etscalc)
+        assert np.allclose(amse_fitets, amse_etscalc)
+        assert np.isclose(lik_fitets, lik_etscalc)
+        print(time_fitets)  # noqa
+        print(time_etscalc)  # noqa
+        return True
+
+
+if __name__ == "__main__":
+    # np.set_printoptions(threshold=np.inf)
+    # test_ets_comparison(setup, 241, False)
+    SUCCESSES = True
+    for i in range(0, 30000):
+        SUCCESSES &= test_ets_comparison(setup, i, True)
+    if SUCCESSES:
+        print("Test Completed Successfully with no errors")  # noqa
diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py
index 814a7531a9..03ac18ee18 100644
--- a/aeon/forecasting/base.py
+++ b/aeon/forecasting/base.py
@@ -119,7 +119,7 @@ def forecast(self, y, exog=None):
         y = self._convert_y(y, self.axis)
         return self._forecast(y, exog)
 
-    def _forecast(self, y=None, exog=None):
+    def _forecast(self, y, exog=None):
         """Forecast values for time series X."""
         self.fit(y, exog)
         return self.predict(y, exog)
diff --git a/aeon/forecasting/tests/test_ets.py b/aeon/forecasting/tests/test_ets.py
index ba92d76cb1..8be3e9f799 100644
--- a/aeon/forecasting/tests/test_ets.py
+++ b/aeon/forecasting/tests/test_ets.py
@@ -1,16 +1,80 @@
 """Test ETS."""
 
+__maintainer__ = []
+__all__ = []
+
 import numpy as np
 
-from aeon.forecasting import ETSForecaster
+from aeon.forecasting import ETSForecaster, ModelType
+
+
+def test_ets_forecaster_additive():
+    """TestETSForecaster."""
+    data = np.array(
+        [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12]
+    )  # Sample seasonal data
+    forecaster = ETSForecaster(
+        alpha=0.5,
+        beta=0.3,
+        gamma=0.4,
+        phi=1,
+        horizon=1,
+        model_type=ModelType(1, 1, 1, 4),
+    )
+    forecaster.fit(data)
+    p = forecaster.predict()
+    assert np.isclose(p, 9.191190608800001)
+
+
+def test_ets_forecaster_mult_error():
+    """TestETSForecaster."""
+    data = np.array(
+        [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12]
+    )  # Sample seasonal data
+    forecaster = ETSForecaster(
+        alpha=0.7,
+        beta=0.6,
+        gamma=0.1,
+        phi=0.97,
+        horizon=1,
+        model_type=ModelType(2, 1, 1, 4),
+    )
+    forecaster.fit(data)
+    p = forecaster.predict()
+    assert np.isclose(p, 16.20176819429869)
+
+
+def test_ets_forecaster_mult_compnents():
+    """TestETSForecaster."""
+    data = np.array(
+        [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12]
+    )  # Sample seasonal data
+    forecaster = ETSForecaster(
+        alpha=0.4,
+        beta=0.2,
+        gamma=0.5,
+        phi=0.8,
+        horizon=1,
+        model_type=ModelType(1, 2, 2, 4),
+    )
+    forecaster.fit(data)
+    p = forecaster.predict()
+    assert np.isclose(p, 12.301259229712382)
 
 
-def test_ets_forecaster():
+def test_ets_forecaster_multiplicative():
     """TestETSForecaster."""
     data = np.array(
         [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12]
     )  # Sample seasonal data
-    forecaster = ETSForecaster(alpha=0.5, beta=0.3, gamma=0.4, season_len=4)
+    forecaster = ETSForecaster(
+        alpha=0.7,
+        beta=0.5,
+        gamma=0.2,
+        phi=0.85,
+        horizon=1,
+        model_type=ModelType(2, 2, 2, 4),
+    )
     forecaster.fit(data)
     p = forecaster.predict()
-    assert np.isclose(p, 15.85174501127)
+    assert np.isclose(p, 16.811888294476528)