As a summer side project I decided to implement a Theta forecaster in aeon. I had heard of it before implementing it, but didn’t really understand how it worked, so I will also include some description. It was proposed in 2000 in this paper:

Assimakopoulos, V., & Nikolopoulos, K. (2000). The Theta model: a decomposition approach to forecasting International Journal of Forecasting, 16(4), 521–530.
https://doi.org/10.1016/S0169-2070(00)00066-2

A later paper states “The original description of the method given by Assimakopoulos and Nikolopoulos involves several pages of algebraic manipulation.” [2]. The actual Theta algorithm is really very simple. Medium doesn’t seem well set up for maths, so I will also direct you to our notebook for a fuller description, but you don’t need pages of maths to understand it. It is a very fast forecaster that is surprisingly effective.

It was the overall winner in the M3 forecasting competition across all 3,003
series, outperforming other statistical models like ARIMA, exponential smoothing, and more complex approaches. It was particularly strong on yearly and quarterly data, and it became widely cited as a benchmark thereafter.

In the M4 Competition (2018) Theta was still competitive, but pure Theta was outperformed by hybrid and ensemble methods like Theta+ETS+ARIMA combinations, and machine learning–enhanced versions. The best-performing entries in M4 (e.g., Hybrid ES-RNN and Theta combinations) included Theta as a component in an ensemble rather than using it alone.

Makridakis et al. (2020), The M4 Competition: Results, findings, conclusion and way forward, International Journal of Forecasting, 36(1), 54–74.

This makes it surprising that I don’t see it in recent forecasting comparisons, which tend to stick to ETS/ARIMA models.

But importantly, why is it called Theta? Obviously it has a parameter called Theta, but why that? It seems an arbitrary choice, just a spare greek letter not used in ARIMA I think.

The Algorithm

The core idea behind theta is to fit a linear trend, then model variation around that trend. The steps of the algorithm are:

1. Fit a linear trend of y against time.
2. Form the theta line to combined linear trend and original data.

3. Smooth the theta line using simple exponential smoothing (SES) to obtain the level.

4. Combine the level and trend to form a forecast.

Theta is the crucial parameter. For 0 it is the trend, for 1 it is the original data. Values greater than one inflate the theta line away from the trend line. The original default was Theta of 2.

The SES alpha parameter is usually found through a grid search, and weight for forecasting is set in the constructor, defaulting to 0.5 (equal weight for level and trend).

In [2] Theta is characterised as ETS with drift. I see theta more as a tweaked linear trend. It would be interesting to assess whether the level makes any actual difference in performance, and how it compares to window based regression.

The Implementation

As with ARIMA, our objective is a super fast, super simple stripped down implementation we can use as a component in other forecasters. The Theta implementation is about 150 lines long (including comments) and the core code is just this (all in numba)

def _fit_predict_numba(y: np.ndarray, h: int, theta: float, weight: float):
    n = len(y)

    # --- Trend component (theta=0)
    t = np.arange(n, dtype=np.float64)
    t_mean = t.mean()
    y_mean = y.mean()
    dt = t - t_mean
    dy = y - y_mean
    var_t = (dt * dt).sum()

    if var_t <= 1e-12:
        b = 0.0
        a = y_mean
    else:
        b = (dt * dy).sum() / var_t
        a = y_mean - b * t_mean

    trend_in = a + b * t
    trend_fut = a + b * np.arange(n, n + h, dtype=np.float64)

    # --- Theta line
    theta_line = theta * y + (1.0 - theta) * trend_in

    # --- Estimate alpha for SES by SSE minimisation
    best_sse = 1e300
    alpha = 0.1
    for k in range(1, 101):
        a_try = k / 100.0
        sse = _ses_sse(theta_line, a_try)
        if sse < best_sse:
            best_sse = sse
            alpha = a_try

    # --- SES forecast
    level = _ses_last_level(theta_line, alpha)
    ses_fut = np.full(h, level)

    # --- Combine components
    w = min(max(weight, 0.0), 1.0)
    forecast = w * trend_fut + (1.0 - w) * ses_fut

    return forecast, a, b, alpha

I’m sure there are lots of variants, but we want super fast, and this fits the bill.

The Comparison

You can run comparable models in the two toolkits like this.

from aeon.forecasting.stats import Theta
from statsmodels.tsa.forecasting.theta import ThetaModel
from statsforecast import StatsForecast
from statsforecast.models import Theta as SFTheta
# statsmodels
sm= ThetaModel(y_train, period=None, deseasonalize=False)
res =  sm.fit()
# statsforecast
pred = res.forecast(1)[0]
    df = pd.DataFrame({
        "unique_id": ["series_1"] * n,
        "ds": pd.date_range("2000-01-01", periods=n, freq="D"),
        "y": y_train,
    })
sf = StatsForecast(models=[SFTheta(season_length=1)], freq="D")
pred = sf.forecast(1, df)
# aeon
aeon = Theta().fit(y_train) # default theta=2.0, weight=0.5
pred = aeon.predict(y_train)

Versions are statsforecast 2.0.2 and statsmodels 0.14.5. Experiments can be reproduced with this notebook. Until version 1.3 is released you will need to install aeon from main like this

pip install git+https://github.com/aeon-toolkit/aeon.git@main

You may have to uninstall 1.2 first if you have it in the environment already. Version 1.3 should be released by the end of August 2025.

Set up is for each n we generate 50 random series of trended data. Train the model on n-1 points, measure squared error on point n and average. Left as you look below is the MSE for increasing n, right is average run time.

statsmodels is slightly different (and worse), maybe I have not configured it the same. statsforecast and aeon are near identical. For run time, surprisingly, statsforecast is slower than statsmodels. I have not included the time to shove things into a data frame, but I suspect data frame interaction is the cause of the slow down. Perhaps they are using more complex linear model implementations. Maybe the Nixtla folk have a super fast version on the horizon. It’s worth putting the difference into context. The algorithm is insanely fast, and if it takes five seconds to fit a series length 20,000 that will probably be ok for nearly all use cases. It pleases me as a computer scientist to get orders of magnitude speed up though :)

Processing series length n=20000 (50 repeats)…
aeon MSE=0.261192 Avg time=4.360518e-03s
statsmodels MSE=0.263235 Avg time=2.438506e-01s
statsforecast MSE=0.260944 Avg time=5.971252e+00s

The Future

Like ETS, Theta is not designed for forecasting on a horizon. It will always simply project the linear trend weighted with the last level value (see the weights figure above). Despite this, it is regularly assessed for forecasts over a horizon.

I’m new to this field. I have been surprised how hard very simple forecasters are to beat for short horizon univariate forecasting. This seems to be widely known and accepted, but it still surprises me and those who work in machine learning who I talk to. Perhaps naively, but it seems that it should be possible to consistently beat Theta/ETS/windowed linear regression/random forest with a time series based regression approach. Our overarching research question is this:

Given a reasonably long series and no prior knowledge, what is the best approach for single step ahead forecasting?

In all things, Occam’s razor guides us

Or in research term, if you cannot beat a really simple technique on average over a large number of problems, you have to make a more detailed argument for the necessity of your new algorithm. Just because you have a hammer, it does not make everything a nail. An algorithm being cleverly designed is not justification without a clearly stated and soundly evaluated use case that shows significant improvement on benchmarks.

This also then leads to many follow up questions: what is the best by data characteristic? What is best by problem domain? what is best at iteratively forecasting? All these questions have of course been addressed in one way or another in the past. Everything has been addressed in one way or another in the past. However, I cannot find a definitive answer to these questions. I feel a bake off coming on ….

For Theta in particular, I would like to assess the actual impact of the smoothing stage. Also, looking a little deeper, there are some interesting variants we could explore. A quick search reveals these. Feel free to point me at other work that might improve Theta. I’m particularly interested in its use in ensembles, but I’ve not gone there yet. Perhaps I should actually read the papers below first :)

References
[1] V. Assimakopoulos & K. Nikolopoulos (2000). The Theta model: A decomposition approach to forecasting. International Journal of Forecasting, 16(4), 521–530.
[2] R.J. Hyndman & S. Billah (2003). Unmasking the Theta method. International Journal of Forecasting, 19(2), 287–290.
[3] J.A. Fiorucci, T.R. Pellegrini, F. Louzada, F. Petropoulos & A.B. Koehler (2016). Models for optimising the Theta method and their relationship to state space models. International Journal of Forecasting, 32(4), 1151–1161.
[4] J.A. Fiorucci et al. (2015). The Optimised Theta Method. arXiv preprint arXiv:1503.03529.
[5] E. Spiliotis, K. Nikolopoulos & V. Assimakopoulos (2020). Generalizing the Theta method for automatic forecasting. European Journal of Operational Research, 284(2), 550–558.
[6] G. Sbrana & A. Silvestrini (2025). The structural Theta method and its predictive performance in the M4-Competition. International Journal of Forecasting, 41(3), 940–952.

Some of us at aeon are dabbling with forecasting, and there are several things that really confused us in the field to begin with. I thought it might help others if we spelt out my confusion.

Let’s start with a time series. Keep it univariate for simplicity. And lets not use airline, it’s a bit boring. Tony has worn a fitbit all year, so lets forecast the number of steps he does in a day. Glad to be hitting my 10k a day :)

I’m not looking at the specifics of analysis, seasonality, trend, heteroscedasticity etc. I don’t actually want to forecast my steps, I am using it as an example to demonstrate variants of forecasting horizons.

Our basic functionality is to fit a model and make a single prediction one step ahead. Lets take two aeon forecasters that implement BaseForecaster as example.

from aeon.forecasting.stats import ARIMA
from aeon.forecasting import RegressionForecaster
arma = ARIMA(p=1,q=1)
p = arma.forecast(y_train)  # One ahead forecast
print(f" ARIMA predicts Tony will take {p} steps.")
arma.fit(y_train)
p2 = arma.predict(y_train)
print(f" This is functionally equivant to forecast: {p2} steps.")
arma = ARIMA(p=1,q=1)
reg = RegressionForecaster(window=50)  # Default is scikit-learn LinearRegression
p = reg.forecast(y_train)  # One ahead forecast

But in practice, this is a limited use case. We may want to predict further in the future than how many steps Tony does tomorrow. For example, we may want to predict how many steps Tony does on this day in a weeks time, or predict every day for the next seven days. We may even want a long horizon, and predict how many steps we take for all days remaining in the year. This is usually called the forecasting horizon or just the horizon and often characterised as a short horizon (various definitions of short, we have seen less than 12 used) or long horizon (often defined as “not short”). Very few papers look at single point forecasts (probabilistic is a whole separate issue). We suspect this is because we all like to draw graphs. The fact is that the use case of making multiple predictions in the future complicates the learning algorithm and the interaction between fit and predict. This is the central point of this blog.

The term “horizon” is overloaded to mean two different things. There is the horizon an algorithm is trained to predict, then there is the horizon an algorithm is used to make prediction about in the future. Crucially, they are not always the same thing.

When you train a forecaster, some can only be trained to predict one step ahead. Others can be trained to predict a single value more than one step ahead. Another set of algorithms can be trained to predict several values ahead at once. Lets split it into train/fit and predict/forecast

Training a forecaster (method fit(y))

1. Train to predict one step ahead only. Traditional model based statistical algorithms are always trained to forecast one ahead. ARMA, ETS, Theta etc are basically designed to predict one step ahead. The only way I can think of sensibly training to predict further ahead would be to down sample, but I don’t think that's standard (although having just worked on implementing a threshold AR model, perhaps the link to time is not as vital as I thought).
2. Train to predict more than one step ahead. Machine learning based forecasting and some deep learning forecasting (detailed below) is based on transforming forecasting to regression through a sliding window. This model is easily adapted to predict, say, five steps ahead by simply changing the lag of the predictor. It is important to note we are not predicting all the points up to horizon points. We train to predict exactly horizon points ahead which can be used with a direct strategy (see below).

import pandas as pd
import numpy as np
from aeon.forecasting.stats import ARIMA
from aeon.forecasting import RegressionForecaster
from aeon.regression.interval_based import DrCIFRegressor

data = pd.read_csv("daily_steps.csv")
y = data["steps"].values # Exclude the last value to match the length for prediction
print(y)
print(y.shape)
y_train = y[:-1]  # Use all but the last value for training
y_test = y[-1]  # Use the last value for testing
arma = ARIMA(p=1,q=1)
p = arma.forecast(y_train)  # One ahead forecast
print(f" ARIMA predicts Tony will take {p} steps.")
reg = RegressionForecaster(window=50)  # Default is scikit-learn LinearRegression
p = reg.forecast(y_train)  # One ahead forecast
print(f" Basic regression predicts Tony will take {p} steps tomorrow.")
# You can use aeon time series specific regressors like DrCIFRegressor
reg = RegressionForecaster(regressor=DrCIFRegressor(n_estimators=10), window=50,
                           horizon=7)
p = reg.forecast(y_train)  # Seven ahead forecast
print(f" DrCIF regression predicts Tony will take {p} steps on this day next week.")
print(" Actual steps taken one ahead:", y_test)

In aeon we signify the ability to forecast more than one ahead with a boolean capability tag called “capability:horizon”. If set to true, the forecasting horizon can be set in the constructor, as shown in the example above.

3. series to series forecasters.

Now we are in the deep learning world. Historically, with the rise of Sequence2Sequence modelling in Natural Language Processing, deep learning models started to use such modelling technique and train the model in an auto-regressive approach, meaning recursive prediction during training which is too slow.

Nowadays, deep learners can be trained to predict a whole series in one go. They use a window as with the point regressors, but the target is a whole series. Stats people would call this functional regression. Deep learning people don’t really have a name to differentiate it and call it all sorts of different things. Some just say “forecasting” some say “multi output” and “direct multi output” where some say explicitly Series2Series or seq2seq.

That is the three forms of training. But what about prediction?

Predicting a series of future values

Suppose we want to forecast Tony’s steps for every day next week (rather than for tomorrow or for this day in a weeks time as in the example above). Simple with a series to series forecaster, just train to predict a window of seven. But there are two variants for point forecasters, and it is important to make the distinction between training the algorithm and using it to make predictions. In aeon, we specify three MixIns to give a forecaster the capability for iterative, direct or series to series forecasting to produce forecasts over a specified prediction horizon.

1. Iterative forecast

If a forecaster can only predict one ahead, how do we use it to predict seven ahead? One way is to fit one model, predict one ahead, then pretend the predicted value was the actual value, append it to the original series and feed it back to get the next forecast. This is often called a recursive strategy, but we don’t like that term (its not actually recursive)

preds = np.zeros(prediction_horizon)
self.fit(y)
for i in range(0, prediction_horizon):
    preds[i] = self.predict(y)
    y = np.append(y, preds[i])
return preds

Some forecasters cannot iteratively forecast in this way. For example ETS needs to track time forward, and ARIMA cannot work out residuals without knowing the true values. These forecasters override the method iterative_forecast rather than use the default version.

2. Direct forecast

If a forecaster is capable of learning to forecast more than one ahead, an alternative is to fit a separate model for each horizon.

preds = np.zeros(prediction_horizon)
for i in range(0, prediction_horizon):
    f = _clone_estimator(self)
    f.horizon = i + 1
    preds[i] = f.forecast(y, exog)
return preds

For regression forecasters, this is achieved through making the regressor y horizon points ahead of the window like this

# Create windowed data for X
X = np.lib.stride_tricks.sliding_window_view(
    combined_data, window_shape=(combined_data.shape[0], self.window)
)
X = X.squeeze(axis=0)
X = X[:, :, :].reshape(X.shape[0], -1)

# Ignore the final horizon values for X
X_train = X[: -self.horizon]

# Extract y_train from the original series
y_train = y.squeeze()[self.window + self.horizon - 1 :]

3. Series to Series forecasting

The mixin for series to series is simply this. We are working on bringing some of the deep learners into aeon, should have some in 1.3. Series to series forecasters could in theory be used with iterative or direct, but that would be a bit weird.

So back to Tony’s steps, these are the direct and iterative forecasts for the same regression forecaster to produce two weeks of forecast:

reg = RegressionForecaster(window=50)  # Default is scikit-learn LinearRegression
steps = 14
series1 = reg.iterative_forecast(y_train, prediction_horizon=14)
series2 = reg.direct_forecast(y_train, prediction_horizon=14)

More on aeon forecasters is in our notebooks

How this can confuse evaluation of forecasters

The problem comes when strategies are mixed when forecasters are compared. To illustrate potential confusion over horizons, I’ll talk about the univariate forecasting experiments published in this paper.

https://www.vldb.org/pvldb/vol17/p2363-hu.pdf

I think its a really good paper and VLDB is an A* conference. I am not trying to do a number on this work, I just want to use it to highlight possible confusion caused by mixing evaluation strategies.

For each of the 8064 univariate series, the test data is set to the last F points, set according to the sampling frequency as set in the M4 competitions.

They evaluate stats, machine learning and deep learning forecasters. For the machine learning algorithms Linear Regression (is that really machine learning?), XGBoost and Random Forest and they use windows of 1.25*F. This seems weirdly short to me, but is another legacy of M4. They use a “direct multi-output” strategy which is we believe simply a direct strategy: they fit different models for each horizon, each trained on series length of N-F. The evaluation of ML algorithms is really assessing how well LR, XGBoost and RF perform over a range of horizon.

The also compare stats models ARIMA, ETS, KF. For these, they do not use direct or iterative. Instead, they fit F models for using all the past data, i.e. N-F, N-F+1, … N-1. This is really a completely different strategy, and is useful for testing how well a forecaster can predict one ahead. Comparing these results to those for direct forecasting is not comparing like for like.
They also compare against 14 deep learning based models. RNN-based, CNN-based models (MICN, TimesNet, and TCN), MLP based models (NLinear, DLinear, TiDE, N-HiTS, and N-BEATS), Transformer-based models (PatchTST, Crossformer, and FEDformer, Non-stationary Transformers (Informer, and Triformer), and Model-Agnostic models (FiLM). All of these are we think treated as series to series forecasters in the paper and in implementations such as Nixtla (we think), but they were not necessarily designed to be used in that way.

Things then get even more confusing, because deep forecasters mix all three possible train regimes with different test mechanisms. The simplest are those that predict a single horizon ahead. These include DeepAR and Temporal Convolutional Network (TCN). These could be used with a direct or iterative strategy, but are almost always used with iterative because of the cost of fitting.

With Direct Multi-Output Forecasting (DMO) forecasting, a deep learning model is trained to predict the full horizon during training. The internal loss function used in optimisation is based on the MSE of the whole predicted window. This is a series to series forecasting model, and distinct from a direct approach because you fit only a single model.

DMO often ignores dependencies between the future steps since it treats each output in the forecast horizon as conditionally independent of the others. This can lead to incoherent or inconsistent predictions over time, especially for long horizons where step-to-step relationships matter.
DMO models are often used to make very long horizon forecasts. Horizons of 720 steps based on windows of length 96 are common. See this excellent talk by Christoph Bergmeir at a NeurIPS workshop on the limitations of this approach.
Examples include LSTNet, Temporal Fusion Transformer (TFT), N-BEATS,, Informer and Autoformer.

Finally, Sequence to sequence (Seq2Seq) forecasting is arguably the most complex of the three DL approaches. Originating in NLP tasks like machine translation, Seq2Seq has been adapted for time series to learn a mapping between input sequences (past observations) and output sequences (future observations), essentially turning forecasting into a sequence transformation problem.

At first glance, this may seem similar to DMO/ series to series: both output full sequences, but there’s a crucial difference: Seq2Seq uses an autoregressive decoder, meaning each predicted value is conditioned on previous predictions. This allows the model to explicitly learn dependencies between future steps, something DMO does not do.
The architecture typically consists of:

• An encoder: processes the input sequence and compresses it into a latent representation.
  - A decoder: predicts future values step by step, feeding its own output back in at each step.

Example:

Suppose we have a time series: {Y_1, Y_2, Y_3, Y_4, Y_5, Y_6, Y_7, Y_8}
We want to use the first 5 time steps as input to forecast the next 3 steps, so the forecast horizon H = 3
- The encoder takes {Y_1, Y_2, Y_3, Y_4, Y_5} and produces a context vector.
- The decoder does the following:
- At time step 1: takes the encoder output and true Y_5 to predict \hat{Y}_6
- At time step 2: takes \hat{Y}_6 and predicts \hat{Y}_7
- At time step 3: takes \hat{Y}_7 and predicts \hat{Y}_8

Training (Teacher Forcing):

During training, instead of using the model’s own (imperfect) predictions, the decoder is teacher-forced, it receives the true previous value at each step, so instead we have:
- At time step 2: uses true Y_6 to predict \hat{Y}_7
- At time step 3: uses true Y_7 to predict \hat{Y}_8
You might wonder how the loss is computed in this case. Even though we use true values as inputs during training, the loss is computed between {Y_6, Y_7, Y_8} and {\hat{Y}_6, \hat{Y}_7, \hat{Y}_8}.
Yes, we are “cheating” a bit by using the ground truth to guide predictions, but this makes training more stable. And since it’s training data, who cares, right?

Inference/predict:
During inference, teacher forcing is not available, so the decoder must use its own previous predictions to generate each next value. This can lead to error propagation, but the autoregressive structure allows the model to learn coherent, temporally-dependent sequences, something DMO cannot do.
Examples: Seq2Seq LSTM, DARNN

That said, Seq2Seq models aren’t as widely used for time series forecasting anymore. The main reason is what happens at inference time: since the decoder has to rely on its own previous predictions (rather than the ground truth like during training), any small error can easily snowball as the forecast progresses, especially over long horizons. This mismatch between training (with teacher forcing) and inference (without it) often leads to degraded performance. On top of that, autoregressive decoding is slower because predictions are made one step at a time. These days, most modern models just predict the entire future horizon in a single shot (aka DMO). It’s faster, avoids error accumulation, and usually performs better, especially for long-term forecasting.

So, that’s deep learning: either point horizon with an iterative approach or series to series via two different training mechanisms. back to the paper, the results for Table 6 compare iterative, direct and series to series forecasters all together, in addition to stats forecasters that can see into the future. This is confusing and conclusions drawn from the results may not hold true. It would we believe be better to compare like with like whenever possible.

We think that evaluation should control as many factors of variation as possible to detect meaningful differences. Both ML and Stats models could adopt an iterative approach, and this is a more natural approach for comparing to DL algorithms. However, our first research question is: given reasonably long series and no other prior knowledge, what can we infer about the relative performance of one ahead forecasters?

These are the experiments we are working on, and the results may be interesting …. aeon forecasting module is developing rapidly. We have two Google Summer of Code students, Tina and Balgopal, working on expanding the machine learning and deep learning sides. Alex Banwell is working on the stats module for his PhD and presenting his first paper at ECML in September.

Watch this space…

By Tony Bagnall and Ali Ismail-Fawaz. Both are aeon core developers. Tony is a professor of Computer Science at University of Southampton. Ali is a postdoctoral researcher at the Université de Haute-Alsace.

Author: Antoine Guillaume

Tired of wasting time on the boring stuff like input type checking and conversion when all you want to do is build awesome time series algorithms? Well, you’re in luck, aeon is here to help!

In this post, we’re diving into how aeon’s base class structure can make your development process smoother and more robust. Imagine being able to focus purely on innovation while aeon takes care of all the stuff related to input type checking, conversion, etc. Sounds great, right?

Whether you’re a seasoned time series practitioner or just starting out, aeon’s got something for everyone. Ready to make your coding life a lot better? Let’s get started!

Wait, what is aeon?

It’s a swiss army knife for time series data. It provides algorithms to solve most time series problems like forecasting, classification, anomaly detection, tools to help you load, process, and evaluate your time series, and much more.

It is a NumFocus affiliated open-source project maintained by a community of developers from both industrial and academic backgrounds working on the field of time series. It aims at being the go-to toolkit for all time series related problems. aeon is designed to work with other popular tools like scikit-learn, so if you’re already familiar with it, you’ll feel right at home!

Want to know more about the project? Check out the documentation website or the GitHub repository!

How aeon helps you to build time series machine learning algorithms

aeon uses a core inheritance hierarchy of classes across the toolkit, with specialised subclasses in each module. Your goal will be to identify the base class that suits the time of algorithm you want to develop, and make your new class inherit from it. The basic class hierarchy is shown in the following diagram. To make sense of this, we will break it down from the top.

Scikit-learn and aeon base estimator

Everything in aeon inherits from sklearn BaseEstimator, which mainly handles the mechanisms for getting and setting parameters using the set_params and get_params methods. These methods are used when the estimators interact with other classes, such as GridSearchCV.

Then we have the aeon BaseAeonEstimator class. This class handles the following for all aeon's estimator:

• management of tags, setting, getting, interaction with sklearn’s tags, etc.
• cloning and resetting of the estimator
• creation of test instances using test parameters specified by each estimator. For example, this is used to define fast-running configuration (e.g. a forest classifier with only 2 trees) for a CI/CD pipelines.

A word on aeon’s tag system

Tags in aeon are used for various purposes. To display estimators capabilities in the documentations, to use specific tests based on each estimator’s capabilities, and more. You can check all existing tags in aeon and the developer documentation on the testing framework to know more about how we use tags for testing purposes.

class RandomDilatedShapeletTransform(BaseCollectionTransformer):

    ...

    # Example of the tags defined by the RDST transformation
    _tags = {
        "output_data_type": "Tabular",
        "capability:multivariate": True,
        "capability:unequal_length": True,
        "capability:multithreading": True,
        "X_inner_type": ["np-list", "numpy3D"],
        "algorithm_type": "shapelet",
    }

One of the main use of tags is input or output formatting and checking made by base classes. Some important tags in this regard are :

• X_inner_type, used to specify the expected type of the input data (Arrays, DataFrames, Lists) used in the implementation. For example, this allows your algorithm to take both numpy arrays and pandas DataFrames as input, while your code uses numpy arrays only.
• output_data_type, used to specify the expected type of the output data (e.g. tabular, series, collections). This is mostly used by transformations estimators.
• capability:multivariate, which indicates whether the estimator can handle multivariate time series data.
• capability:unequal_length, which indicates whether the estimator can handle collection of unequal length time series data.
• capability:multithreading, which indicates whether the estimator can handle multithreading for parallel processing.

Tags can also be used to indicate behaviours and required dependencies that might be optional in the default installation of your package, one of such example for aeon is the LiteTimeClassifier :

class LITETimeClassifier(BaseClassifier):

    ...

    # Example of the tags defined by the LITE Classifier
    _tags = {
        "python_dependencies": "tensorflow",
        "capability:multivariate": True,
        "non_deterministic": True,
        "cant_pickle": True,
        "algorithm_type": "deeplearning",
    }

We will give more practical examples on how to use these tags in the following sections.

Collection estimator and Series estimator

We distinguish between two types of inputs for aeon estimators, series and collections:

• Series represent single time series as a 2D format (n_channels, n_timepoints), some estimators can also use 1D format as (n_timepoints) when they don't support multivariate series. Series estimators also have an axis parameter, which allow the input shape to be transposed such as the 2D format becomes (n_timepoints, n_channels) instead.
• Collections represent an ensemble of time series as a 3D format (n_samples, n_channels, n_timepoints). Again, this can sometime be represented as a 2D format, such as (n_samples, n_timepoints) for univariate estimators. Preferably, this should be avoided to clear any confusion on the meaning of axes and the possible confusion with 2D single series.
• For the case of an unequal length collection, we use a list of 2D arrays, where the n_timepoints can vary for each sample, this is was is referred to as a np-listin the RDST code example.

For example, if we go back to the base class schema BaseClassifier inherit from BaseCollectionEstimator. This means that during fit and predict, all estimators inheriting from BaseClassifier will expect time series collection as inputs.

Collection estimators

The BaseCollectionEstimator defines methods to check the shape of the input, extract metadata (e.g. whether the collection is multivariate) and check compatibility of the input against tags of the estimators. For example, when you do the following :

from aeon.classification.dictionary_based import TemporalDictionaryEnsemble
from aeon.testing.data_generation import make_example_3d_numpy_list
# TDE does not support unequal length collections
# as it sets "capability: unequal_length":False
X_unequal, y_unequal = make_example_3d_numpy_list()
try:
    TemporalDictionaryEnsemble().fit(X_unequal, y_unequal)
except ValueError as e:
    print(e)

Data seen by instance of TemporalDictionaryEnsemble has unequal length series,
but TemporalDictionaryEnsemble cannot handle these characteristics.

What happens here is that TemporalDictionaryEnsemble inherits from BaseClassifier, which itself inherits from BaseCollectionEstimator. During fit and predict, BaseClassifier calls _preprocess_collection, a function defined in BaseCollectionEstimator. This function extracts the input metadata (whether it is multivariate, of unequal lengths etc.) and compare it against TemporalDictionaryEnsemble tags. These state that the estimator does not support unequal lengths collections, and hence an exception is raised.

Classification

This is the base class for all classifiers. It uses the standard fit, predict and predict_proba structure from sklearn. fit and predict call the abstract methods _fit and _predict which are implemented in the subclass to define the classification algorithm. All the common format checking and conversion are done in final functions such as fit, predict and are made before calling the abstract methods _fit and _predict.

When implementing a new classifier inheriting from BaseClassifier, you thus only have to implement the __init__, _fit and _predict methods that handle the classification logic of the classifier. You will also need to set the correct tags to allow the check and conversion to be done for you. Note that each base class also defines some attributes that are commonly used in the estimators, for example BaseClassifier exposes classes_, n_classes_, _class_dictionary that we can use in our new classifier:

from numpy.random import default_rng
from aeon.classification import BaseClassifier
from aeon.testing.data_generation import (
    make_example_3d_numpy,
    make_example_dataframe_list,
)

class RandomClassifier(BaseClassifier):
    """A dummy classifier returning random predictions."""
    _tags = {
        "capability:multivariate": True,  # allow multivariate collections
        "capability:unequal_length": True,  # allow multivariate collections
        "X_inner_type": ["np-list", "numpy3D"],  # Specify data format used internally
    }
    def __init__(self, random_state: int | None = None):
        self.random_state = random_state
        super().__init__()
    def _fit(self, X, y):
        self.rng = default_rng(self.random_state)
        return self
    def _predict(self, X):
        # generate a random int between 0 and n_classes-1 and
        # use _class_dictionary to convert it to class label
        return [
            self._class_dictionary[i]
            for i in self.rng.integers(low=0, high=self.n_classes_, size=len(X))
        ]
# A 3D numpy array 
X, y = make_example_3d_numpy(n_channels=2)
print(RandomClassifier().fit_predict(X, y))
# A list of dataframes, each representing a 2D Series
X, y = make_example_dataframe_list()
print(RandomClassifier().fit(X, y).predict(X))

[1 0 1 1 0 0 1 0 1 0]
[0 1 1 0 0 1 0 1 0 0]

Further reading

• Introduction to time series classification using aeon

Regression, Clustering and Anomaly detection

These base classes are mostly similar to BaseClassifier in how they use the checks and conversion operations from BaseCollectionEstimator.

• BaseRegressor also defines a fitand predict method and requires _fit and _predict methods to be implemented by child classes. The difference is that it has no predict_proba method yet, as we still need to decide how to model probabilistic regression for time series. The tests on y are also different, as we can have floats has values for y.
• BaseClusterer also has fit and predict, but does not take input y as child classes can be unsupervised estimators. It does include predict_proba.
• BaseCollectionAnomalyDetector also has fit and predict, but does not take input y as child classes can be unsupervised estimators.

Further reading

• What is time series extrinsic regression ?
• Time series clustering using aeon
• Time series anomaly detection using aeon

Collection transformation

Rather than fit andpredict, the BaseCollectionTransformer implements fit, transform and fit_transform methods, which are required since this base class also inherits the BaseTransformer class. It will require child classes to define _fitand _transform methods. The output of the transform method is not fixed and should be specified with the output_data_type.

For example, if the output is another collection of time series (e.g. after using SAX), then output_data_type must take the Collection value (note that this is the default value for all BaseCollectionTransformer child classes). If the output is not time series anymore, but rather a 2D array of features extracted from each input time series, such as in Rocket or RandomShapeletTransform, then the output_data_type must take the Tabular.

Further reading

• Time series transformation with aeon

Series estimators

Series estimators are similar to collection estimators, but they take single time series as input. They inherit from BaseSeriesEstimator, which perform similar operations to BaseCollectionEstimator, such as input checks and conversions, but for single time series.

One important difference is the use of the axis parameter, which allows each estimator to define whether it works with the (n_channels, n_timepoints) or the (n_timepoints, n_channels) 2D format. We need to have the axis parameter, because for 2D series, there is no way to infer whether an input 2D time series is in the (n_channels, n_timepoints) or (n_timepoints, n_channels) format.

To understand its uses, we need to distinguish between the axis parameter set during initialisation of the estimator, and the axis parameter used in the fit, predict and other methods:

• axis given during initialisation is used to define the internal format used by the estimator,
• axis given when call functions is used to transpose the input time series if needed, to match the format internally used by the estimator.

Note that the axis value represent the dimension in which the number of timepoints is stored, so axis=0 means that the timepoints are stored in the first dimension, and axis=1 means that the timepoints are stored in the second dimension (i.e. (n_channels, n_timepoints)).

Further reading

• aeon series estimators

Series transformation

Let’s take the example of BaseSeriesTransformer, which is the base class for all series transformers. It implements the fit, transform and fit_transform methods, and requires child classes to implement _fit and _transform. Let use demonstrate the use of the axis parameter with an example:

from aeon.testing.data_generation import make_example_dataframe_series
from aeon.transformations.series import BaseSeriesTransformer
class DummySeriesTransformer(BaseSeriesTransformer):
    """A dummy series transformer that keeps every second timepoint."""
    _tags = {
        "capability:multivariate": True,  # allow multivariate series
        "X_inner_type": "pd.DataFrame",  # Specify data format used internally
        "fit_is_empty": True,  # we don't need to define _fit
    }
    def __init__(self):
        super().__init__(axis=1)  # Set axis to 1 for (n_channels, n_timepoints) format
    def _transform(self, X, y=None):
        print(X.shape)
        X = X.iloc[:, ::2]  # Example transformation: keep every second timepoint
        print(X.shape)
        return X

X = make_example_dataframe_series(n_channels=2, n_timepoints=10).T
print(X.shape)  # Is (n_timepoints, n_channels), which is axis=0
print(DummySeriesTransformer().fit_transform(X, axis=0).shape)

(10, 2)
(2, 10)
(2, 5)
(5, 2)

What we see is that X starts of shape (n_timepoints, n_channels), which is the default format for the input time series.

The DummySeriesTransformer is initialised with axis=1, meaning that internally, when calling fit and transform, the input is converted to (n_channels, n_timepoints) format before being passed to the _fit and _transform functions. The output is then converted back to the user specified format axis=0 : (n_timepoints, n_channels) before returning it.

Note that as we specified pd.DataFrame as the X_inner_type in the tags, if the input is not a DataFrame, the estimator will convert it to a DataFrame before applying the transformation. This allows you to define and use a single input shape and type in the function you implement, while still allowing the estimator to handle different input formats.

from aeon.testing.data_generation import make_example_2d_numpy_series
X = make_example_2d_numpy_series(n_channels=1, n_timepoints=10)
transformer = DummySeriesTransformer()
print(transformer.fit_transform(X, axis=1).shape)

(1, 10)
(1, 5)
(1, 5)

Forecasting

The BaseForecaster class inherits from BaseSeriesEstimator which provides the checks and conversion functions for single time series inputs. Forecasters predict future values of a time series horizon steps ahead. The horizon parameter is defined during initialisation to specify how many steps ahead the forecaster should predict. The main methods are :

• fit(y, exog=None): Trains the forecaster on input data
• predict(y, exog=None): Makes predictions on new data
• forecast(y, exog=None): Combines fit and predict into one operation

For child classes, the _fit and _predict methods must be implemented to define the forecasting logic. It also provides two main forecasting strategies:

• direct_forecast(): Makes predictions by training separate models for each future step.
• iterative_forecast(): Makes predictions recursively using a single trained model.

Further reading

• Time Series Forecasting with aeon

Anomaly detection and Segmentation

The BaseSeriesAnomalyDetector and the BaseSegmenter base classes both implements the fit, predict but only requires child classes to implement the _predict, as most anomaly detector and segmenters are unsupervised estimators. Both also inherit from BaseSeriesEstimator, which provides the checking and conversion functions we have seen before.

Further reading

• Introduction to time series segmentation

What’s next?

Now, you can start experimenting and making your own estimators more robust and field-ready with aeon base classes and tags. If you have questions or would like to help us in expanding aeon, you can come have a chat with us on Slack, have a look at open issues on GitHub!

tl:dr; we compare the new aeonforecaster with the python implementations in statsmodelsand statsforecast, and show that it is two orders and one order of magnitude faster respectively, with no loss in predictive performance. The aeonARIMA forecaster is blazingly fast, but comes with some caveats…

Some aeoncore developers have begun dabbling with forecasting. We are all computer scientists, so its a bit outside our comfort zone, which is good. The task that interests us is short horizon forecasting of long series, as we think this is where we may be able to advance the field. Ultimately we are interested in machine learning algorithms, but new algorithms should be benchmarked with the traditional approaches, so we are starting with classical techniques. Broadly speaking, there are three classes of forecasters: statistical models, machine learning algorithms and deep learning approaches. There is obviously overlap, but there is a clear distinction in terms of the communities that work in each area and the specification of problem they address. This is a whole blog post in itself and for another day. In this article we will describe our first steps in developing statistical algorithms for forecasting in python. Specifically, we will give an overview of the current stats forecasting software in python, describe our implementations of some favourites like Auto Regressive Integrated Moving Average (ARIMA) and Exponential Smoothing (ETS), then present benchmarking results against the current most popular approaches. We hope this encourages you to try our forecasters out. We are not going to fill the post with standard descriptors of the algorithm. This information is freely available elsewhere. This post is of interest to someone wanting to use fast, simple implementations of stats based forecasters. This blog focuses on the aeon ARIMA implementation, but we are adding more classical algorithms (e.g. Theta and TVP) and will look in depth at ETS in another article.

Forecasting software for ARIMA and ETS

The oldest and most widely known package is statsmodels. It seems to have been the basis for all python based stats forecasting for the last decade. Historically, most forecasting software is in R. We are not using R for reasons beyond the scope of this article. statsmodels fits ARIMA and ETS like this:

from aeon.datasets import load_airline
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.arima.model import ARIMA as SM_ARIMA
y = load_airline().squeeze() # Make a 1D numpy
# ARMA Model
arma = SM_ARIMA(y[:-1], order=(1, 0, 1)).fit()
forecast = arma.forecast(steps=1)[0]
print(f"Actual = {y[-1]} Forecast = {forecast}")
# ETS Models
ets = ExponentialSmoothing(y,trend="add", seasonal="add", seasonal_periods=12).fit()
forecast = ets.forecast(steps=1)[0]
print(f"Actual = {y[-1]} Forecast = {forecast}")

statsmodels has no native automatic hyper-parameter fitting algorithm: it has to be done manually through grid search or bespoke methods. Wrappers like pmdarma, sktime and Darts encapsulate this process, but all use statsmodels implementations of the core algorithm. These packages offer other features, such as pipelines and transforms, but they do not have bespoke ETS and ARIMA implementations, which is our requirement.

statsforecast is a more recent package from Nixtla who apparently have an “industry-leading time series engine”. It contains bespoke implementations of both ARIMA and ETS, and also has auto versions to fit the hyper-parameters. To fit a model, you have to wrap it in a DataFrame

from statsforecast import StatsForecast
from statsforecast.models import ARIMA, AutoARIMA, AutoETS
import pandas as pd
df = pd.DataFrame({
    "ds": pd.date_range(start="1900-01-01 00:00:00", periods=len(y),
                        freq="H"),
    "y": y,
    "unique_id": "ts"
})
sf = StatsForecast(models=[ARIMA(order=(1,0,1))],freq="H")
forecast_df = sf.forecast(df=df, h=1)
f1 = forecast_df.iloc[0]["ARIMA"]
sf = StatsForecast(models=[AutoETS()],
                   freq="H")
forecast_df = sf.forecast(df=df, h=1)
f2 = forecast_df.iloc[0]["AutoETS"]
sf = StatsForecast(models=[AutoARIMA()],
                   freq="H")
forecast_df = sf.forecast(df=df, h=1)
f3 = forecast_df.iloc[0]["AutoARIMA"]

Why have another implementation?

So as far as I can tell, that’s it: either use statsmodelsand/or their various wrappers, or use statsforecastand their DataFrame model.

Our research requirement is fast implementations of classical algorithms in order to benchmark against the new techniques we are developing. We are not interested in all possible variants of algorithms that have ever been proposed. Neither of the existing packages fulfil our need. Whilst comprehensive, statsmodelsis designed like R and is quite slow. statsforecastis faster, but is built around DataFrames. We will implement forecasters that are:

1. Built around numpy arrays and numba;
2. Simple to use and concise in code;
3. As fast as we can make them.

ARIMA and ETS have meta parameters and parameters. For ARIMA, the meta parameters are (p,d,q), the number of autoregressive, differences and moving average terms.

Fitting the model

First source of confusion is that there are two steps to fitting a model. Fitting the hyper parameters (p, d and q for ARIMA) and fitting the parameters (usually called theta and phi)

An optimiser can be used to fit parameters and/or meta parameters for ARIMA and ETS. statsmodelshas a wide range of optimisers. chatgpt tells me it has these:

Apparently, for ARMA it uses a Kalman filter statespace model fit to maximise likelihood using “L-BFGS-B” (via scipy.optimize.minimize) optimiser. blimey. ETS in statsmodelsdoes something similar. Personally I was suprised they don’t just use stochastic gradient descent, but old school numerical analysis techniques are cool, so we will go with that.

statsforecast takes a different approach (so my LLM tells me). For ARMA both AR and MA terms are solved via closed-form least squares, and ETS uses discrete grid search. Now, I think I remember how to fit an AR model without optimisation, you simply find the PACF and use the diagonal terms, stopping with some criteria such as AIC, but not sure about ARMA. Further searching and digging suggests chatGPT is mistaken, ARIMA at least uses BFGS (via optim() from base R) but its not at all clear.

We chose to implement Nelder-Mead optimiser for aeon. Why? First, we read a linkedin post recommending it (which I can no longer find, maybe from Ivan?). Secondly, it looked simpler than the alternatives to implement. Thirdly, it seemed less widely used, so could differentiate us from the others. Finally, it seemed quite easy to generalise to any parameter fitting.

For the hyper parameters such as (p,d,q), the space of meta parameters is so small, it seems to me sensible to just enumerate them. A range of [0…5) for p and q and [0..3) for d is only 75 models and the fitting is very fast. It seems like statsmodels and statsforecastuse an optimiser, possibly some form of forward selection, but we will only do so if it seems necessary.

So, decisions made, we have implemented everything with numpy arrays, speed up with cached numba and use a lightweight version of Nelder-Mead to fit parameters and exhaustive search to fit a limited number of meta parameters. Will it be any good? A little about our forecasting base class first.

Forecasting in aeon

One thing that we had to work through is the different use cases in forecasting. I’m writing a separate post about horizons in forecasting, because if we are confused others may be too. Our notebooks on forecasting are here . We do use fit and predict, but our simplest use case is:

d.forecast(y): fit model and make a single prediction.

d.iterative_forecast(y, prediction_horizon=12): fit a single model, then predict 12 steps ahead, feeding predictions back into the model

d.direct_forecast(y, prediction_horizon=12): fit 12 different models, each trained to predict a different number of steps.

from aeon.datasets import load_airline
from aeon.forecasting import NaiveForecaster, RegressionForecaster
y = load_airline()
d = NaiveForecaster(strategy="last")
p = d.forecast(y)   # One ahead forecast
d = NaiveForecaster(strategy="seasonal_last", seasonal_period=12)
p2 = d.iterative_forecast(y, prediction_horizon=12)
d2 = RegressionForecaster(window=100) # Defaults to LinearRegression
p3 = d2.direct_forecast(y, prediction_horizon=12)

With ARIMA and ETS, you cannot simply train models to predict more than one ahead, so they do not implement direct_forecast.

These are statsforecastconstructor arguments. Its great that it is so configurable, but it is a little confusing to the uninitiated and all these parameters introduce overhead in validation etc.

AutoARIMA (d:Optional[int]=None, D:Optional[int]=None, max_p:int=5,
            max_q:int=5, max_P:int=2, max_Q:int=2, max_order:int=5,
            max_d:int=2, max_D:int=1, start_p:int=2, start_q:int=2,
            start_P:int=1, start_Q:int=1, stationary:bool=False,
            seasonal:bool=True, ic:str='aicc', stepwise:bool=True,
            nmodels:int=94, trace:bool=False,
            approximation:Optional[bool]=False, method:Optional[str]=None,
            truncate:Optional[bool]=None, test:str='kpss',
            test_kwargs:Optional[str]=None, seasonal_test:str='seas',
            seasonal_test_kwargs:Optional[Dict]=None,
            allowdrift:bool=True, allowmean:bool=True,
            blambda:Optional[float]=None, biasadj:bool=False,
            season_length:int=1, alias:str='AutoARIMA', prediction_interva
            ls:Optional[statsforecast.utils.ConformalIntervals]=None)

Our ARIMA and AutoARIMA model constructors look like this

class ARIMA(BaseForecaster):

def __init__(
    self,
    p: int = 1,
    d: int = 0,
    q: int = 1,
    use_constant: bool = False,
    iterations: int = 200,
):

We ran an experiment to compare our ARMA to statsmodelsand statsforecast versions.

The set up is this. We specified 20 ARMA models with p and q between 0 and 5 and parameters set to ensure stationarity. We are not considering differenced data since all algorithms do this in the same way.

We pick a random model and generate a series, then z-normalise it. We train on n-1 points, and make a single forecast. For each n, we repeat this process 100 times, and calculate the MSE for the three toolkit implementations.

The graphs below present the results. On the left as you look are the results for fitting a model with the correct p and q values. On the right the results for a standard default for p=1 and q=1. The baseline results are very similar. There is variation for the true models, with aeonshowing more variation than the others and seemingly worse on longer series. We investigate if this effect is significant next. The timing results are quite impressive. aeonARMA is an order of magnitude faster than statsforecastand two orders faster than statsmodels.

To assess whether there is any significant difference in MSE we ran a separate experiment using just statsforecast(because we are impatient). We fix n now to 2500 (seemingly the worst value for aeon) then create 500 series for our 20 models, fit aeonand statsforecastARIMA as before (with the correct p and q values) then measure the error of the last forecast point for each. We have included a Naive forecaster for a sanity check forecaster.

The short story is there is no significant difference. The average number of wins per model for aeon was 250.25/500. The average pairwise difference over 20 problems was -0.0000000273. aeon was better on 10 problems, statsforecast on 10 problems. The overall MSE for the two approaches is 0.678569 and 0.678584 (the MSE for Naive is 1.16). I could present graphs etc, but its all saying the same thing: To all intents and purposes, they are the same. aeon took 50 seconds to run the experiment, statsforecast took over 10 minutes. It is perfectly possible we are not using the super fast version of statsforecast , I heard mention of a precompiled C++ version, but we are following their instructions on usage.

Nelder-Mead is very fast, but not without issues. Sometimes it converges to a very poor model which gives massive MSE. This happens rarely for the simulated data we used, maybe once every few thousand runs, but its worth being aware of it. Basically its sensitive to the initial conditions. We have altered the default initial conditions to make it more robust: it did not happen in the above experiments. However, given how fast it is, we might fit twice with different initial simplex and take the best. This would make the chance of a pathological model very small indeed. Further, ultimately we will use it with an AutoARIMA forecaster. In this case, it would just mean the pathological model would not be selected.

All this is very preliminary and will evolve into a paper describing the aeon forecasting module, we have a lot planned for it.

Our goal was a lightweight super fast algorithm, and we have achieve this. Our forecaster is not as configurable as the others, nor does it use pandas, so it wont be for everyone. However, if you want to benchmark on a large number of series (as we do), it might appeal to you. Just watch out for the odd misfit!

ARIMA is available on main will be in release 1.3. We hope to also have ETS, AutoARIMA and AutoETS in the same release. We have put together a notebook to reproduce this.

Words and experiments by Tony Bagnall (github @TonyBagnall). Tony is a core developer of aeon and a professor of Computer Science at the University of Southampton.

Big credit to Alex Banwell (@alexbanwell1) who has lead the implementation of Nelder-Mead and the stats forecasters as part of the first year of his PhD.

aeonis an international open source project affiliated with NumFOCUS.

GitHub - aeon-toolkit/aeon: A toolkit for machine learning from time series

We communicate mostly on slack, come say hello (link from github). Tony is happy to answer any questions on slack or by email. Tell me if I’ve made any mistakes in the above, it has been known to happen.

We will be at ECML/PKDD in Portugal in September 2025, if you are there come say hi at the AALTD workshop on Friday, where Alex presents some of his forecasting research.

aeon is supported by the EPSRC in the UK and has recently secured significant funding outside of the UK …. news to follow

Reference us with this paper