1 Introduction

Topological data analysis (TDA) utilizes concepts from algebraic topology to analyse the shape of data, making it a pivotal tool for extracting meaningful insights from complex datasets [1]. At its core, TDA employs persistent homology to unveil the underlying structure of data by identifying and tracking features such as connected components, loops, and voids across different scales [2]. Persistent homology, a cornerstone of TDA, systematically studies the multi-scale topological characteristics of datasets. It constructs simplicial complexes at different scales and tracks the evolution of features like connected components, loops, and voids [3]. This method captures the persistence of these features across scales, offering insights into the robustness and significance of topological structures within the data [4].

Persistent homology results are typically visualized using persistence diagrams or barcodes, which provide a concise yet powerful representation of the topological features and lifetimes inherent in a dataset [3]. These visual tools condense complex mathematical insights into accessible formats, enabling researchers to discern persistent patterns across various scales of analysis [2]. By highlighting enduring features amidst fluctuations and noise, persistence diagrams facilitate the discovery of fundamental structures that may not be discernible through conventional data analysis methods alone [5]. Moreover, these capabilities enhance the interpretability of intricate datasets and support informed decision-making across diverse scientific and practical domains, ranging from biology and neuroscience to materials science and beyond. Finally, persistence diagrams yield metrics such as entropy, amplitude, and point counts, which can be instrumental in quantitative assessments and predictive modelling [6].

Financial markets are complex systems where prices of assets evolve according to various factors, including economic indicators, market sentiment, and external events [7]. Traditional financial analysis methods rely on time series analysis, statistical models, and machine learning techniques [8,9,10]. However, these methods may fail to capture the intricate and multi-scale financial data structure [11]. TDA, employing persistent homology, introduces a new approach to examining financial markets [12]. It provides insights into the overall data structure, essential for comprehending market behaviour and predicting critical events like market crashes [13]. Stable and significant features in data are identified by persistent homology, thereby offering connections to underlying market dynamics. Its resilience to noise and ability to manage high-dimensional data makes it highly suitable for financial applications. Furthermore, financial time series data can be represented as a topological space by TDA, allowing for the identification of meaningful clusters, loops, and voids, which provide valuable insights for decision-making [14].

Topological data analysis (TDA) application to financial markets is relatively new but has shown promising results in various areas. One significant application of TDA is in detecting critical transitions in financial markets. For instance, Gidea (2017) [12] attempted to detect early signs of critical transitions-such as financial crashes-by tracking changes in the topology of correlation graphs as prices approached the 2007-2008 financial crisis using persistent homology. In another study focusing on the detection of critical transitions, Gidea and Katz (2018) [13] used persistent homology to identify topological patterns in a multidimensional time series formed by the four major financial indices in the USA during the technology crash of 2000 and the financial crisis of 2008. They employed persistence landscapes and their \(L^p\) norms to measure changes in one-dimensional persistence diagrams. Their results indicated that the \(L^p\) norms increased significantly days before the two crashes, while the norms remained stable during calm market periods. Continuing the exploration of critical transition detection in financial markets, researchers applied TDA to study crashes in cryptocurrency markets. Gidea et al. (2020) [15] analysed the time series of four major cryptocurrencies: Bitcoin, Ethereum, Litecoin, and Ripple. They defined the \(C^1\) norm of persistence landscapes to identify critical transitions preceding the cryptocurrency crash in January 2018. Additionally, they employed the k-means clustering algorithm to discern patterns in price movements.

TDA has also been employed to manage risk and guide potential investment strategies. By understanding the topological features of financial data, TDA can contribute to more effective risk management strategies. Persistent homology aids in identifying periods of increased market turbulence, enabling investors to adjust their portfolios accordingly [16]. Additionally, researchers have developed TDA-based indicators to enhance investment strategies. These indicators, derived from the persistence landscapes of financial time series, offer new insights into market behaviour and have demonstrated superior performance compared to traditional indicators in specific scenarios [17].

From a methodological perspective, within the domain of machine learning, TDA has been applied across various domains, including pattern recognition, classification, and signal quality evaluation [18,19,20,21,22]. In the realm of time series analysis, TDA has also demonstrated its utility [23, 24]. However, its application has predominantly been limited to classification settings, for identifying certain phenomena in dynamical systems, or for purely exploratory analysis [25], rather than being utilized for forecasting purposes. To the best of our knowledge, no existing manuscripts currently explore the use of TDA specifically for forecasting tasks. This apparent gap in the literature may arise from the prevailing belief that forecasting models should, through their learning processes, inherently capture and integrate the structural information within historical observations. In contrast, the present study proposes a novel approach by directly capturing these features in the raw data through TDA. Consequently, this manuscript seeks to advance forecasting models by integrating supplementary features derived from TDA.

In pursuit of this objective, this study will construct point clouds and generate persistence diagrams to derive complex topological features, including persistent entropy, amplitude, and the number of points within the persistence diagram. These features capture critical non-linear dependencies and topological invariants, facilitating a multi-scale analysis of financial time series. Persistent entropy quantifies the complexity and irregularity within the time series, indicating volatility and underlying market dynamics. Amplitude, reflecting the spread in the persistence diagram, provides insights into the stability and recurrence of observed patterns. The number of points representing significant topological features across scales highlights the persistence of structures within the data, indicating recurrent cycles or trends that may otherwise remain undetected. These attributes collectively could enhance computational forecasting models by incorporating topological summaries, offering a robust and adaptable framework suited to the inherently non-stationary nature of financial data.

Furthermore, it is essential to emphasize that these features are extracted using a sliding window approach over the time series. For each segment, a point cloud is constructed, and a persistence diagram is generated, from which the three features above are extracted and incorporated into the forecasting model to enhance performance. This innovative approach differs from traditional methods that typically generate a single point cloud for the entire time series; instead, multiple point clouds are constructed to preserve crucial temporal information and capture evolving dynamics. Consequently, local structures are detected, non-stationary data are accommodated, and robustness against outliers is achieved, resulting in a more nuanced and temporally aware analysis well-suited to the complexities of financial time series. Finally, the proposed feature extraction methodology will be integrated into the N-BEATS model, which has been selected due to its state-of-the-art performance in time series forecasting and its ability to handle complex patterns.

The remainder of the manuscript is structured as follows: Sect. 2 presents the proposed methodology for extracting TDA features and integrating them into the N-BEATS model. Section 3 outlines the experimental framework for comparing the proposed model against traditional feature extraction methods using financial datasets. Section 4 reports the forecasting results across various datasets and metrics, highlighting significant improvements in predictive accuracy. Sect. 5 emphasizes the effectiveness of TDA-derived features in enhancing time series forecasting and suggests future research directions, including applications in reinforcement learning (RL) and high-frequency trading (HFT). Finally, Sect. 6 summarizes the proposal and its principales resultados.

2 The proposed method

The forthcoming section introduces an innovative approach to time series forecasting utilizing topological data analysis (TDA). This section comprehensively describes the TDA techniques employed in this study. The primary objective is to extract features that characterize the underlying temporal patterns in a novel manner and compare them against traditional techniques to establish a benchmark. This approach aims to validate whether TDA-derived features can effectively predict a target variable.

Consequently, this section begins by defining the problem of point forecasting in univariate time series. It then justifies the necessity of incorporating feature extraction mechanisms in financial time series forecasting due to the inherent complexity of such problems. Subsequently, it details the extraction of TDA features from the time series using a sliding window approach. Finally, it outlines how these features are integrated into a forecasting model based on neural networks, specifically the N-BEATS model.

2.1 Problem statement

An univariate time series is defined as a sequence of data points x(t), \(\mathcal {T} = \{x(t) \mid t = 1, 2, \ldots , N\},\) where N represents the number of points in the time series and x(t) denotes the value of the series at time t. In the case of financial time series data, x(t) typically represents variables such as stock prices, trading volumes, interest rates, and exchange rates, with measurements taken at regular intervals (daily, monthly, annually).

Given this context, the feature extraction technique proposed in this manuscript aims to address the challenge of (point) forecasting in univariate time series. The primary objective is to utilize the last T observations of the time series to generate point forecasts for a specified time horizon H. Formally, this problem can be described as learning a function, parameterized as a neural network, \(f: \mathbb {R}^{T} \rightarrow \mathbb {R}^{H},\) from a given collection of input–output pairs, where inputs are length-T vectors, and outputs are length-H vectors. For simplicity, a prediction horizon of one point will be considered, i.e. \(H = 1\).

2.2 Justification for feature extraction techniques in financial time series

Forecasting in financial markets involves using historical data to predict future market trends, leveraging statistical methods and machine learning algorithms. The primary objective is identifying patterns and trends that inform investment decisions, risk management, and strategic planning. The necessity for robust feature extraction techniques in financial time series forecasting is underscored by the complex and dynamic nature of financial markets, which are characterized by:

  • Non-Stationarity: Financial time series data are often non-stationary, with changing statistical properties. Techniques such as differencing are employed to achieve stationarity before model implementation [26].

  • High volatility and noise: Unpredictable factors induce high volatility and noise in financial markets, complicating the identification of genuine patterns. Feature extraction methods are used to reduce noise and isolate relevant patterns [27].

  • Complex interdependencies: Financial time series are influenced by interrelated factors, necessitating their consideration for improved forecasting accuracy. Feature extraction techniques incorporate these relationships into models [28].

  • Overfitting and underfitting: Overfitting captures noise, while underfitting results from overly simplistic models. Feature extraction mitigates these issues by generating informative features [29, 30].

  • Structural breaks: Sudden changes in the data-generating process disrupt established patterns. Detecting structural breaks is crucial for forecast accuracy, and feature extraction techniques assist in this adaptation [31].

Given these challenges, the necessity of robust feature extraction techniques in financial time series forecasting becomes evident. Feature extraction enhances the predictive power of forecasting models by transforming raw data into informative features that better capture the underlying patterns and trends. By addressing non-stationarity, high volatility, and complex interdependencies, feature extraction techniques improve the accuracy and reliability of financial forecasts. Effective feature extraction allows models to gain deeper insights into the data, making it easier to identify trends and predict future market movements [32].

2.3 Feature extraction via topological data analysis techniques

This manuscript performs feature extraction at the segment level using a slicing window approach rather than across the entire time series. This decomposition divides the original series \(\mathcal {T} = \{x(t) \mid t = 1, 2, \ldots , N\} \in \mathbb {R}^{N}\) into Q overlapping windows of length L, resulting in \(\mathcal {T} = (\textbf{s}_1, \ldots , \textbf{s}_{Q}) \in \mathbb {R}^{Q \times L}\) where \(Q = N-L+1\). Each segment \(\textbf{s}_q = (x(q), \ldots , x(q+L-1)) \in \mathbb {R}^L\) represents a distinct subset of the original data. Topological characteristics specific to each segment are then computed using advanced techniques from TDA, such as single Takens’ embeddings, point clouds, persistence diagrams, and Vietoris-Rips complexes. The application of these techniques is detailed for each segment of the time series \(\mathcal {T}\), aiming to elucidate localized structures and dynamics within financial data. This approach harnesses sophisticated analytical tools to uncover nuanced patterns and relationships that traditional methods may obscure.

Firstly, single Takens’ embeddings are implemented, a method utilized to reconstruct the state space of time series and form a crucial part of time delay embedding techniques in TDA. The mathematical foundation of this approach is derived from Takens’ embedding theorem, which asserts that, under certain conditions, the dynamics of a system can be fully reconstructed from a series of observations [33]. Given a time series segment \(\textbf{s}_q\), a point cloud, \(\mathcal {C}_q\), in \(\lambda \)-dimensional space is generated by the single Takens’ embedding method, where each point \(\textbf{x}_i\) in \(\mathcal {C}_q\) is defined as:

$$\begin{aligned} \mathcal {C}_q = \left\{ \textbf{x}_i \mid \textbf{x}_i = \left( s_{q,i}, s_{q,i+\tau }, \ldots , s_{q,i+(\lambda -1)\tau } \right) , i = 1, 2, \ldots , J \right\} , \end{aligned}$$

where \(\tau \in \mathbb {Z}_{+}\) denotes the time delay, \(\lambda \in \mathbb {Z}_{+}\) represents the embedding dimension, and \(J = L - (\lambda - 1)\tau \) denotes the number of points in the point cloud \(\mathcal {C}_q\). Here, \(s_{q,i}\) refers to the value at position \(i\) within segment \(\textbf{s}_q\). The parameters \(\tau \) and \(\lambda \) significantly influence the dynamical system’s phase space reconstruction. \(\tau \) governs the temporal spacing of sequential observations in the time series, thereby ensuring the effective distribution of points in the reconstructed space. Its optimal selection minimizes redundancy between successive coordinates, enhancing their independence and capturing essential dynamical information. Conversely, \(\lambda \) determines the number of delayed observations utilized in constructing the phase space. A judicious choice of \(\lambda \) ensures a faithful representation of the original system’s dynamics, thereby revealing the discernible structure of the underlying attractor.

In this proposal, the point cloud associated with each \(q\)-th segment, \(\mathcal {C}_q\), serves as the initial step towards visualizing the inherent geometric structure within financial datasets. By projecting data points into multidimensional space, analysts gain a geometric perspective to observe clusters, trends, and outliers, thus providing a foundational insight into the market’s geometrical configuration [34]. From the point cloud, a filtration process tracks the birth and death of topological features (such as connected points, loops, and voids) across successive stages. This process involves incrementally “thickening” the points, visualized as increasing the radius of balls centred at each point in \(\mathcal {C}_q\). Mathematically, for \(\mathcal {C}_q = \{\textbf{x}_1, \textbf{x}_2, \ldots , \textbf{x}_J\}\), the filtration can be represented as a sequence of subspaces:

$$\begin{aligned} \emptyset = \textbf{u}_0 \subseteq \textbf{u}_1 \subseteq \textbf{u}_2 \subseteq \ldots \subseteq \textbf{u}_J = \mathcal {C}_q \end{aligned}$$

where each \(\textbf{u}_i\) includes points from \(\mathcal {C}_q\) based on proximity criteria typically defined by a metric, and \(\emptyset \) represents the empty set.

The Vietoris-Rips complex provides a specific approach to construct such a filtration from a point cloud. At a given scale parameter \(\epsilon \), denoted \(VR_\epsilon (\mathcal {C}_q)\) for the \(q\)-th segment, it includes a simplex for every subset of points in \(\mathcal {C}_q\) that are pairwise within distance \(\epsilon \) of each other. This complex grows as \(\epsilon \) increases:

$$\begin{aligned} VR_{\epsilon _0}(\mathcal {C}_q) \subseteq VR_{\epsilon _1}(\mathcal {C}_q) \subseteq \ldots \subseteq VR_{\epsilon _k}(\mathcal {C}_q) \end{aligned}$$

where \(0 = \epsilon _0< \epsilon _1< \ldots < \epsilon _k\). As the filtration progresses, new topological features may appear (birth) and existing features may merge or disappear (death) [3, 35]. For example, a loop births when a set of points forms a cycle that does not enclose a filled-in area, and a loop dies when the space inside it becomes filled-in.

Persistence diagrams, denoted as \(\mathcal {D}_q\) (\(q \in \{1,\ldots , Q\}\)), represent these topological features by plotting \((b_i, d_i)\) points, where \(b_i\) and \(d_i\) are the birth and death times of the \(i\)-th feature, respectively. These diagrams capture enduring structures amidst market volatility, shedding light on long-term market behaviours [36]. Hence, formally, the persistence diagram associated with the \(q\)-th segment is defined as \(\mathcal {D}_q = \{(b_i, d_i) \mid i \in I_q\}\), where \(b_i\) and \(d_i\) represent the birth and death times of the \(i\)-th feature, respectively, and \( I_q \) is the index set enumerating the topological features identified in the persistence diagram associated with the q-th segment.

Once the persistence diagrams of the different segments have been computed, the next step is to calculate the characteristics of interest for each persistence diagram. In our study, these metrics will be the entropy, amplitude, and number of points, as these three measures provide deeper insights into the underlying market dynamics. Firstly, the persistent entropy associated with the \(q\)-th segment, \(H(\mathcal {D}_q)\), is defined as:

$$\begin{aligned} H(\mathcal {D}_q) = - \sum _{i \in I_q} p_i \log p_i \end{aligned}$$
(1)

where \(p_i = \frac{r_i}{R}\), \(r_i = d_i - b_i\) is the lifetime of the \(i\)-th feature, and \(R = \sum _{i \in I_q} r_i\) is the total sum of lifetimes. It measures the uncertainty associated with the persistence diagram. High-entropy values indicate a more complex or chaotic underlying system, while lower values suggest more regular or periodic behaviour [37].

Secondly, the amplitude of the \(q\)-th segment, \(A(\mathcal {D}_q)\), is defined as:

$$\begin{aligned} A(\mathcal {D}_q) = \max _{i \in I_q} r_i = \max _{i \in I_q} (d_i - b_i), \end{aligned}$$
(2)

and refers to the size of the topological features detected in the data. This can be interpreted as a measure of the ‘spread’ or ‘size’ of the features in a particular homological dimension. Different norms or metrics (e.g. bottleneck and Wasserstein) can be used to quantify this measure.

Thirdly, the number of points for the \(q\)-th segment, \(N(\mathcal {D}_q)\), is the number of distinct topological features (such as connected components and loops) identified in the data at various filtration values. It indicates the complexity or richness of the topological structure inherent in the data [38]. \(N(\mathcal {D}_q)\) is defined as:

$$\begin{aligned} N(\mathcal {D}_q) = |I_q|. \end{aligned}$$
(3)

Figure 1 illustrates the TDA process for a specific time series segment. In Subfigure a, the original time series data segment under analysis is depicted. This segment is transformed using the single Takens embedding method with a time delay \(\tau =72\) and an embedding dimension \(\lambda =2\) to create a point cloud, \(\mathcal {C}\), as shown in Subfigure b. Subsequently, Subfigure c presents the persistence diagram, \(\mathcal {D}\), derived from the point cloud \(\mathcal {C}\). The persistence diagram indicates the presence of connected components (H0), loops (H1), and voids (H2), showcasing the birth and death times of these topological features across the filtration process.

Fig. 1
figure 1

Illustration of the feature extraction process for a segment of a time series: a Example segment of a time series on which feature extraction will be performed; b Resulting point cloud diagram obtained using Single Takens embedding with a time delay of \(\tau = 72\) and embedding dimension \(\lambda = 2\); c Persistence diagram generated from b, illustrating connected components (\(H0\)), loops (\(H1\)), and voids (\(H2\))

As indicated at the beginning of this section, the feature extraction process is repeated for all segments of the series following a sliding window approach. To illustrate the complete process, Fig. 2 depicts the division of the series into overlapping segments of size L, generating a total of Q windows. For each window, the corresponding point cloud and persistence diagram are computed. The values of entropy, amplitude, and the number of points are calculated from these diagrams. These values are subsequently stored in a vector, which for the q-th segment, is defined as \(\textbf{v}_q = (H(\mathcal {D}_q), A(\mathcal {D}_q), N(\mathcal {D}_q)) \in \mathbb {R}^{3}\).

Fig. 2
figure 2

Illustration of the feature extraction procedure via TDA, computed on time series observations \( x(1), \ldots , x(N) \). The time series, \(\mathcal {T}\), is decomposed into a collection of \( Q \) overlapping windows of length \( L \). For each window \( q \), its corresponding point cloud, \(\mathcal {C}_q\), is computed, from which the persistence diagram, \(\mathcal {D}_q\), is derived. The persistence diagram is then used to compute the values of entropy, amplitude, and number of points, which are stored in the vector \(\textbf{v}_q = (H(\mathcal {D}_q), A(\mathcal {D}_q), N(\mathcal {D}_q)) \in \mathbb {R}^{3}\)

2.4 The N-BEATS forecast model

This section offers a comprehensive overview of the N-BEATS model, emphasizing its core elements and the incorporation of TDA features into its architecture. The first subsubsection outlines the fundamental components and operational mechanisms that distinguish N-BEATS as a robust framework for time series forecasting. The second subsubsection details how TDA features are integrated into the N-BEATS model to enrich its capabilities for univariate forecasting.

2.4.1 Overview of key elements of N-BEATS

The neural basis expansion analysis for time series forecasting (N-BEATS) represents a significant advancement in time series forecasting, embodying a paradigmatic shift in generating forecasts by offering a versatile and robust framework for handling univariate and multivariate time series data. Developed by a team of researchers, N-BEATS has demonstrated state-of-the-art performance across various benchmarking datasets, distinguishing itself through a combination of simplicity, effectiveness, and robustness. The model’s deep learning architecture, employing fully connected layers and a backward and forward residual stacking mechanism, decomposes the time series into trend, seasonality, and residual components, thereby capturing complex patterns and dependencies in financial data. Its basis expansion mechanism expands the input time series into a set of basis functions and enhances adaptability to diverse forecasting tasks. The iterative refinement process within its block structure ensures progressively improved forecasts, leading to higher accuracy and reliability [39].

Notably, N-BEATS is flexible, capable of handling time series data with missing values, irregular intervals, and multiple frequencies, and has shown superior performance over traditional statistical models and other deep learning models across multiple benchmarks. For example, N-BEATS demonstrates accuracy improvements of 11% over a statistical benchmark, including the ARIMA baseline, and 3% over the 2020 M4 competition winner, which is based on a hybrid model combining a neural residual/attention dilated LSTM stack with a classical Holt-Winters statistical model. The model’s interpretability, facilitated by its basis expansion mechanism, provides crucial insights into different components of the time series, aiding financial analysts in understanding the driving factors behind forecasts. Despite its complex architecture, N-BEATS is relatively simple to implement and use, requiring minimal pre-processing or feature engineering. It is scalable for efficiently handling large datasets, making it suitable for real-time forecasting scenarios. Its robustness to noise and outliers, critical for maintaining accuracy in unpredictable market events, further underscores its utility in various financial forecasting tasks, including stock prices, exchange rates, and economic indicators [40].

The decision to integrate TDA features as supplementary inputs in N-BEATS was primarily motivated by the model’s significant goodness of fit. While the model exhibits high performance, further enhancing forecasting accuracy with additional features presents challenges, given its already strong fit. Moreover, two compelling reasons for including TDA features in N-BEATS are its modular architecture, which facilitates interpretability crucial for understanding predictions in a financial context and its scalability, which ensures robust performance when managing large datasets, thereby rendering it particularly well-suited for various financial applications.

2.4.2 Incorporation of TDA features into N-BEATS architecture

In this section, the integration of TDA features, as shown in Fig. 2, for each segment of the time series into the N-BEATS model [39], is detailed. The standard N-BEATS model is briefly described, with a note that TDA features can similarly be integrated into the basis expansion variant without modifications. The generic N-BEATS model is constructed from a stack of P double-residual blocks. For \(1 \le p \le P\), each block consists of a non-linear map (implemented as a multi-layer perceptron, MLP). The representation of the N-BEATS model and its variant incorporating TDA features can be observed in Fig. 3, where the P double-residual blocks, their inputs, and their outputs are depicted. As can be seen, the structure of the input and output variables in the blocks for both models (the baseline N-BEATS model and the variant with TDA features, \(\texttt {N-BEATS}_\mathrm {+TDA}\)) is the same, with the particularity that in the TDA-enhanced model, each block receives the corresponding input \(\textbf{x}\) along with the TDA characteristics as additional input, denoted as \(\textbf{v}\). The entry of these TDA characteristics is represented by a red line in Fig. 3.

Fig. 3
figure 3

Illustration of the N-BEATS architecture for univariate time series forecasting (left) and the enhanced N-BEATS architecture incorporating TDA features stored in the vector \(\textbf{v}\) (right)

In the following, the operation of the model will be formally defined starting from the input data \(\textbf{x}\). Firstly, the non-linear mapping of the p-th MLP can be defined as:

$$\begin{aligned} S^p: \mathbb {R}^T \rightarrow \mathbb {R}^h, \end{aligned}$$
(4)

where \(T\) is the dimension of the inputs (length-\(T\) vectors), and \(h\) denotes the internal dimensionality. The twofold output of each block (see Fig. 3) is produced as follows:

$$\begin{aligned} \textbf{x}^p = \textbf{x}^{p-1} - \textbf{U}^p (S^p(\textbf{x}^{p-1})), \end{aligned}$$
(5)

and

$$\begin{aligned} {{\hat{\textbf{y}}}}^p = \textbf{V}^p (S^p(\textbf{x}^{p-1})), \end{aligned}$$
(6)

with \(\textbf{U}^p \in \mathbb {R}^{T \times h}\), \(\textbf{V}^p \in \mathbb {R}^{H \times h}\), and \(\textbf{x}^0 = \textbf{x}\), where \(H\) is the output dimension (assumed to be one for simplicity). The outputs \(\textbf{x}^p\) are used as inputs to subsequent computation blocks, while the predicted outputs \({{\hat{\textbf{y}}}}^p\) are utilized to compute the final model prediction \({{\hat{\textbf{y}}}} \in \mathbb {R}^H\) via component-wise summation:

$$\begin{aligned} {{\hat{\textbf{y}}}} = {{\hat{\textbf{y}}}}^1 + {{\hat{\textbf{y}}}}^2 + \cdots + {{\hat{\textbf{y}}}}^P. \end{aligned}$$
(7)

Importantly, in this deep neural network forecasting model, the outputs \(\textbf{x}^p\) can be utilized as an interface to integrate additional information. Specifically, in this manuscript, the base input signal to each block is enhanced by concatenating the TDA features stored in the vector \(\textbf{v}\) (see Fig.3), as follows:

$$ {\mathbf{x}}_{{TDA}}^{{\text{p}}} = ({\mathbf{x}}^{p} ,{\mathbf{v}})\quad {\text{where}}\quad {\mathbf{v}} = (H({\mathcal{D}}),A({\mathcal{D}}),N({\mathcal{D}})),{\text{ }} $$
(8)

and, as can be observed in Fig. 3, the time series signal \(\textbf{x}\) is provided to the N-BEATS model in its raw form, while the corresponding TDA features, \(\textbf{v}\), are supplied to each block as complementary signals.

2.5 Computational complexity

This section examines the computational complexity involved in incorporating TDA features into the N-BEATS model, which represents the primary contribution of this manuscript. The computation of the persistence diagram is identified as the most significant aspect of TDA feature extraction. The most computationally intensive component within this process is the calculation of the Vietoris-Rips complex, which is characterized by a runtime complexity of \( \mathcal {O}(2^J) \) [41, 42]. It should also be noted that these calculations are repeated \( Q \) times for each of the \( Q \) segments of the time series. Consequently, the computational complexity of this stage is mainly determined by the number of segments and the value of \( J \), where \( J \) is defined as \( J = L - (\lambda - 1)\tau \). Thus, it is also dependent on the length of each segment \( L \), the dimensionality of the point cloud embedding \( \lambda \), and the time delay \( \tau \).

3 Computational experiments

The computational experiments in this study are designed to compare the performance of traditional feature engineering techniques with the proposed approach that incorporates TDA-derived features. This section outlines the key elements of the experimental framework. Specifically, the following subsections define the datasets, consisting of data from six cryptocurrencies evaluated across four scenarios. Pre-processing techniques, including handling missing data, outlier detection, and stationarity tests, are detailed, along with the parameters used in the baseline forecasting model, N-BEATS.

Given that the primary objective is to assess performance improvements due to the inclusion of additional features, the N-BEATS model is consistently used as the baseline forecasting model across all feature extraction techniques. Therefore, a subsection is also dedicated to defining the traditional feature extraction techniques utilized for comparison with the proposed TDA-derived features. Additionally, the performance metrics employed to evaluate the various methods are specified, ensuring a comprehensive forecasting accuracy and effectiveness assessment.

3.1 Datasets

The selected datasets include time series data for six major cryptocurrencies-BNB (Binance Coin), BTC (Bitcoin), ETH (Ethereum), XRP (Ripple), ADA (Cardano), and LTC (Litecoin)-all denominated in USDT (Tether) [43, 44], along with additional financial instruments such as the CAC40 index (FCHI), iShares Core MSCI Europe ETF (IEUR), the EUR/USD exchange rate (EURUSD = X), and crude oil futures (CL = F).

For cryptocurrencies, USDT, a stablecoin pegged to the USD at a 1:1 ratio helps mitigate fluctuations in fiat currencies, allowing a more explicit focus on cryptocurrency market behaviour. This approach enables the precise isolation of intrinsic factors and ensures consistent comparisons across assets and periods. As a result, it allows for detailed examinations of how individual assets respond to different stimuli, enhancing the understanding of the broader digital asset ecosystem [45].

In addition to cryptocurrencies, the experimental framework includes several traditional financial instruments that further broaden the analysis. The CAC40 (FCHI) is a major French stock market index representing the 40 most significant companies on the Euronext Paris, providing a benchmark for tracking Paris market performance. The iShares Core MSCI Europe ETF (IEUR) offers a diversified portfolio of European equities, serving as a proxy for European market trends and occasionally exhibiting bond-like behaviour in terms of stability. The EUR/USD exchange rate (EURUSD = X), one of the most traded currency pairs globally, reflects the relative strength of the European and US economies, making it critical for international trade analysis. Finally, crude oil futures (CL = F) track the price of West Texas Intermediate (WTI) crude oil, a vital commodity whose price movements significantly impact global economic conditions. Together, these instruments complement the cryptocurrency data by providing a holistic view of diverse asset classes, enabling cross-asset analysis and a comprehensive understanding of market dynamics.

The cryptocurrency data spans from July 1, 2019, to December 28, 2021, capturing a complete market cycle with high volatility, bull and bear markets, and stabilization phases. This period is crucial for assessing cryptocurrency performance about macroeconomic factors and regulatory changes [46]. Data are collected at 1-h intervals, balancing granularity and manageability. This frequency captures short-term market movements while reducing data overload, facilitating the analysis of intraday patterns and the impact of global events. Four distinct scenarios are created for each cryptocurrency to evaluate the model comprehensively. These scenarios allow for varied training and testing conditions, enhancing the robustness of the analysis.

In addition to cryptocurrencies, traditional financial instruments such as the CAC40 index (FCHI), crude oil futures (CL = F), the EUR/USD exchange rate (EURUSD), and the iShares Core MSCI Europe ETF (IEUR) are also analysed. Each of these instruments is evaluated using a consistent scenario structure. Specifically, two scenarios are defined for each asset: The first covers from April 1, 2023, to November 14, 2023, and the second extends from November 15, 2023, to June 29, 2024. This uniform approach ensures that model performance is assessed consistently across different asset classes and periods.

Finally, Table 1 outlines the 32 scenarios considered in the study, each assigned a unique ID for ease of reference. The table details the specific cryptocurrency or traditional financial instrument tested and the periods designated for training and testing. These selected assets and timeframes were chosen for their ability to offer a thorough and nuanced perspective of the cryptocurrency market.

Table 1 Characteristics of the time series scenarios for cryptocurrencies and traditional financial instruments used in the experimental study. Each entry includes the scenario ID, asset name, training period, and testing period

3.2 Data pre-processing

This section outlines the methodologies employed in the data pre-processing stage, including the treatment of missing values, computation of log returns, management of outliers, assessment of stationarity, normalization of data, and application of denoising techniques. These procedures have been implemented sequentially to ensure the data’s integrity, accuracy, and analytical suitability [47].

Missing values were addressed using cubic spline interpolation, selected for its ability to fit a smooth curve through existing data points, thereby preserving continuity in the time series. This method is favored over linear techniques due to its effectiveness in managing the complexities inherent in financial time series data, as it approximates missing values by fitting cubic polynomials between data points, ensuring continuity of the first and second derivatives [48].

After addressing the missing values, log returns were computed to facilitate analysis. Log returns are calculated by taking the natural logarithm of the current price ratio to the previous price. This transformation is essential in financial time series analysis as it stabilizes variance and enhances the data’s suitability for statistical modelling. The resulting log returns series constitutes the univariate series that will be the focus of subsequent experiments.

Identifying and adjusting outliers are essential for maintaining the robustness of the analysis, as they can skew results and lead to erroneous conclusions. In this study, outliers were detected using the Z-score method, standardizing the data to highlight values deviating from the mean. While a traditional Z-score threshold of ±3 is standard, this study adjusted the bounds to the 0.02 and 0.98 quantiles to better accommodate the heavy-tailed nature of financial returns. Observations outside these bounds were examined and adjusted to the nearest bound value, preserving data integrity while mitigating the influence of outliers.

Stationarity is a fundamental assumption in time series modelling, indicating that statistical properties such as mean, variance, and autocorrelation remain constant over time. The Augmented Dickey-Fuller (ADF) test was employed to assess stationarity, which enhances the Dickey-Fuller test by including lagged differences to address higher order autocorrelation. The null hypothesis of the ADF test states that the series has a unit root (is non-stationary), while the alternative hypothesis suggests stationarity. In this study, the ADF test was applied to the log-return series, and a sufficiently low p-value would indicate the rejection of the null hypothesis, confirming the series’ stationarity [49].

After verifying stationarity, the data were normalized using the MinMaxScaler, scaling values to a range between -1 and 1, which enhances the performance of machine learning algorithms. The final pre-processing step involved denoising the data through the maximal overlap discrete wavelet transform (MODWT). This method effectively decomposes the time series into multiple resolution components, facilitating noise separation from the signal. The Daubechies 2 (db2) wavelet [50] was chosen for its compact support and smoothness, and five decomposition levels were applied. The first-level approximation, which captures the most significant underlying trends while filtering out high-frequency noise, was retained for subsequent analysis. Finally, the data pre-processing stage is illustrated in detail in the flowchart diagram presented in Fig. 4.

Fig. 4
figure 4

Flow diagram of the data pe-processing stage

3.3 Feature extraction methods used for comparison purposes

Given that the proposal includes TDA-derived features in future forecasting models and, as previously indicated, the N-BEATS model is employed as the forecasting model for all tested methods, the proposed strategies for comparison can be divided into two categories. First, the benchmark model, the N-BEATS, approaches the forecasting problem from a purely univariate perspective, using the T previous data points to make the forecast. The objective of including this model is to evaluate how including new features helps improve the baseline model. The second group of methods included in the experimental study comprises alternative feature extraction models. These models, similar to the proposed strategy, \(\texttt {N-BEATS}_\mathrm {+TDA}\), add additional information from the extracted features to the estimation window originally used by the N-BEATS model. Therefore, the feature extraction strategies considered are:

  1. 1.

    The \(\texttt {N-BEATS}_\mathrm {+LF}\) method incorporates lagged features (LF) as inputs for predicting future values. The underlying idea of this model is that by including lagged returns, trends, and patterns in the data, such as momentum or mean reversion, can be captured [51]. This experimental study added lagged log returns with 1, 3, 7, 14, 21, and 30 days as features of this strategy.

  2. 2.

    The \(\texttt {N-BEATS}_\mathrm {+DTF}\) method integrates date and time features (DTFs) in the inputs included in the baseline N-BEATS model. These elements typically include the hour of the day, the day of the week, and whether the day is a holiday. Such features can capture time-dependent patterns, such as increased trading activity at market open/close or behavioural changes on Fridays compared to midweek [52], as financial markets are influenced by price movements, human behaviours, and operational aspects that vary over time. In this experimental study, the features added within this strategy include the year, month, day, hour, day of the week, day of the year, and number of days in a month.

  3. 3.

    The \(\texttt {N-BEATS}_\mathrm {+DCF}\) method includes difference and change features (DCF) as additional inputs in the baseline N-BEATS forecasting model. These features involve creating new variables based on the differences or percentage changes between various points in the time series. For example, this could include the log-return difference between consecutive days or the percentage change over the last 3 days. Such features help identify market trends and volatility [53]. In this study, the implementation of the \(\texttt {N-BEATS}_\mathrm {+DCF}\) method includes the rate of change and differences using lagged log returns of 1, 3, 7, 14, 21, and 30 days.

  4. 4.

    The \(\texttt {N-BEATS}_\mathrm {+STDF}\) method involves seasonal-trend decomposition features (STDF), which in this study are derived using the seasonal and trend decomposition using LOESS (STL) technique. STL utilizes LOESS (local regression) to decompose time series into seasonal, trend, and residual components, thereby facilitating the analysis and comprehension of intricate patterns within time series data. The trend component signifies the long-term evolution of the series, seasonality reveals recurring patterns occurring at fixed intervals, and the residual denotes what remains after the removal of trend and seasonality [54].

  5. 5.

    The \(\texttt {N-BEATS}_\mathrm {+TDEF}\) method incorporates time delay embedding features (TDEFs), which involves constructing vectors from time series data by selecting a point and its lagged values up to a specified dimension, known as the embedding dimension. This technique is grounded in Takens’ embedding theorem, which proposes that the dynamics of a system can be reconstructed using a single observable, provided that the embedding dimension and time delay are appropriately chosen. In this approach, the time delays considered in the methods are 1, 3, 7, 14, 21, and 30 days.

The proposed model, \(\texttt {N-BEATS}_\mathrm {+TDEF}\), is influenced significantly by two parameters, namely the time delay, \(\tau \), and the segment size used for computing TDA features, L. In the experimental study, potential values of 3 days and 7 days were considered for the time delay (which help in capturing temporal dynamics and dependencies in the data), while segment sizes of 2 weeks, 3 weeks, 1 month, and 1 month and a half were considered, offering a range of temporal granularities for feature extraction.

The hyperparameters of both the comparison methods and the proposed method, \(\texttt {N-BEATS}_\mathrm {+TDEF}\), were estimated using nested cross-validation with a validation window within the training set. During the model selection phase, an exhaustive grid search was performed over the predefined hyperparameter values specified earlier for each method. This ensures independent hyperparameter tuning and prevents overfitting, providing a robust estimate of model performance on unseen data.

Figure 5 portrays a detailed historical TDA feature extraction process. It outlines iterative stages for generating TDA features from pre-processed data. Initially, the data are segmented using specified parameters: time delays, \(\tau \), of 3 days and 7 days and segment sizes, L, of 2 weeks, 3 weeks, 1 month, and 1 month and a half months, ensuring comprehensive temporal coverage. Subsequently, point clouds are derived from these segments to visualize higher dimensional data structures. Persistent diagrams are then constructed to capture topological features across different scales. The entropy, amplitude, and number of points are extracted from these diagrams, distilling complex topological information into interpretable measures. This iterative process iterates through each time slice, facilitating a thorough exploration of temporal dimensions. The historical TDA features are ultimately exported for subsequent analytical or modelling tasks.

Once the historical features have been computed for both the proposed TDA strategy and the comparison methods, the next step involves merging them with the pre-processed data. In this merging process, historical features and pre-processed data are aligned based on date. This step is crucial to maintaining the temporal correspondence between different datasets.

Fig. 5
figure 5

Workflow for historical TDA feature extraction

3.4 Configuration of N-BEATS model parameters

The N-BEATS model is implemented for time series forecasting in this study, known for its robustness and effectiveness across diverse forecasting tasks. Specific hyperparameters are tailored to the dataset and objectives. The N-BEATS model is configured as follows: It uses an input chunk length of 7, \(T\,=\,7\), to consider one week’s historical data per prediction, an output chunk length of 1 for fine-grained forecasts, \(H\,=\,1\), and a batch size of 16 to balance training speed and memory usage. The random state is fixed at 0 for reproducibility. Training spans 100 epochs to capture underlying data patterns effectively while preventing overfitting. With two layers, each 512 units wide, the model captures various levels of data abstraction. The mean squared error (MSE) is the loss function, penalizing prediction errors to enhance accuracy.

These parameters collectively optimize model performance on the selected dataset, balancing complexity and computational efficiency to deliver competitive forecasts. Each parameter has been carefully selected based on preliminary experiments and domain knowledge, ensuring that the N-BEATS model is well-suited to handle the dataset’s specific characteristics and challenges. This configuration leverages the N-BEATS model’s strengths, making it highly effective for the time series forecasting tasks in this study.

3.5 Performance measures

To evaluate the quality of point forecasts, three widely recognized metrics are employed in this study: the mean absolute percentage error (MAPE), the mean absolute error (MAE), and the root mean squared error (RMSE). These metrics have been extensively validated in the forecasting literature for their ability to measure different dimensions of predictive accuracy [39, 55, 56]. MAPE is appreciated for its interpretability and scale independence, although issues can arise with small or zero actual values. A straightforward, unit-specific measure of error magnitude is offered by MAE, though all errors are treated equally, potentially overlooking the impact of larger deviations. A more nuanced assessment is provided by RMSE, which penalizes larger errors more severely, making it suitable for contexts where significant consequences are associated with such errors [57]. By combining these metrics, a comprehensive evaluation of forecasting accuracy is ensured, addressing various aspects of model performance and aligning with best practices in the field.

Let \({{\hat{\textbf{y}}}} = ({\hat{y}}_{N+1}, \ldots , {\hat{y}}_{N+M}) \in \mathbb {R}^{M}\) denotes the forecasting in the testing set, where the testing data have a size of M and follow the training set. Let \(\textbf{y} = (x_{N+1}, \ldots , x_{N+M}) \in \mathbb {R}^{M}\) represents the true observations in the test set. The three performance measures are defined based on these two vectors as follows:

$$\begin{aligned} \text {MAPE} = \frac{\sum _{i=1}^{N} \left| \frac{y_i - {\hat{y}}_i}{y_i} \right| \times 100}{M}, \text {MAE} = \frac{\sum _{i=1}^{N} \left| y_i - {\hat{y}}_i \right| }{M},\\ \text {RMSE} = \sqrt{ \frac{1}{M} \sum _{i=1}^{M} (y_i - {\hat{y}}_i)^2 }. \end{aligned}$$

4 Computational results

In this section, the results obtained for the testing sets are analysed. Specifically, the performance of three forecasting metrics, as produced by the \(\texttt {N-BEATS}_\mathrm {+TDA}\) method, is evaluated against the baseline method, N-BEATS, and five other feature extraction methodologies. The purpose of this analysis is to demonstrate that the incorporation of TDA features into the forecasting algorithm can enhance the performance of the models. In Table 2, the MAPE results for the testing set across 32 forecasting scenarios are presented. Based on the results from these 32 datasets, the ranking of each method for each dataset (\(R = 1\) for the best-performing method and \(R = 7\) for the worst) has been determined. The mean MAPE in the test set (\(\overline{\mathrm{MAPE}_T}\)) and the mean ranking (\(\overline{R}_{\mathrm{MAPE}_T}\)) is also included in Table 2.

Table 2 Comparative performance of the proposed model, \(\texttt {N-BEATS}_\mathrm {+TDA}\), the baseline method N-BEATS, and other feature extraction approaches for the MAPE metric, along with the average MAPE for the 32 scenarios considered in the test set, \(\overline{\mathrm{MAPE}_T}\), and the average ranking of MAPE in the test set, \(\overline{R}_{\mathrm{MAPE}_T}\)

In a similar manner, the performance of the methods used in the experimental study is shown in Tables 3 and 4 for the MAE and RMSE metrics, respectively, across the 32 forecasting scenarios considered in the test set. Additionally, the ranking for each method in each forecasting problem is calculated for each metric. The mean values of MAE and RMSE for the different methods, \(\overline{\mathrm{MAE}_T}\) and \(\overline{\mathrm{RMSE}_T}\), as well as the mean rankings of the different methods, \(\overline{R}_{\mathrm{MAE}_T}\) and \(\overline{R}_{\mathrm{RMSE}_T}\), are also computed.

Table 3 Comparative performance of the proposed model, \(\texttt {N-BEATS}_\mathrm {+TDA}\), the baseline method N-BEATS, and other feature extraction approaches for the MAE metric, along with the average MAE for the 32 scenarios considered in the test set, \(\overline{\mathrm{MAE}_T}\), and the average ranking of MAE in the test set, \(\overline{R}_{\mathrm{MAE}_T}\)

From the analysis of the results, it can be concluded, from a purely descriptive point of view, that the \(\texttt {N-BEATS}_\mathrm {+TDA}\) method achieved the best results for 24 datasets in MAPE, for 18 datasets in MAE, and 19 datasets in RMSE. Additionally, the second best performance was obtained for four additional datasets in MAPE, eight datasets in MAE, and three datasets in RMSE. Furthermore, the \(\texttt {N-BEATS}_\mathrm {+TDA}\) method yielded the best mean performance across all three metrics (\(\overline{\mathrm{MAPE}_T} = 148.8895\), \(\overline{\mathrm{MAE}_T} = 0.1971\), \(\overline{\mathrm{RMSE}_T} = 0.2538\)), as well as the best mean rankings (\(\overline{R}_{\mathrm{MAPE}_T} = 1.5312\), \(\overline{R}_{\mathrm{MAE}_T} = 1.9062\), \(\overline{R}_{\mathrm{RMSE}_T} = 2.0625\)).

Table 4 Comparative performance of the proposed model, \(\texttt {N-BEATS}_\mathrm {+TDA}\), the baseline method N-BEATS, and other feature extraction approaches for the RMSE metric, along with the average RMSE for the 32 scenarios considered in the test set, \(\overline{\mathrm{RMSE}_T}\), and the average ranking of RMSE in the test set, \(\overline{R}_{\mathrm{RMSE}_T}\)
Fig. 6
figure 6

Scatter plot comparing the performance in MAPE, MAE, and RMSE of the proposed \(\texttt {N-BEATS}_\mathrm {+TDA}\) method against the baseline and other comparison methods. The X-axis represents the performance of the \(\texttt {N-BEATS}_\mathrm {+TDA}\) method, while the Y-axis represents the performance of the comparison methods across these metrics

Furthermore, a graphical representation is also provided to substantiate the validity of these results. Figure 6 compares the performance of the methods included in the experimental study against the proposed \(\texttt {N-BEATS}_\mathrm {+TDA}\) method. On the X-axis, the performance of the proposed method is displayed, while the Y-axis shows the performance of the comparison methods. Thus, in the scatter plot of Fig. 6, each point represents a comparison between \(\texttt {N-BEATS}_\mathrm {+TDA}\) and another methodology on a single forecasting dataset. In this regard, points above the \(y = x\) line indicate datasets for which \(\texttt {N-BEATS}_\mathrm {+TDA}\) achieves better performance than the compared method since all considered metrics (MAPE, MAE, and RMSE) are to be minimized. It can be observed that the majority of points lie above this line, underscoring the proposed method’s superior performance compared to the other methods.

In order to determine whether a statistical difference exists between any of these algorithms, a procedure for comparing multiple forecasting models over multiple datasets was employed, as described by Demšar (2006) [58]. The first step involved the application of Friedman’s non-parametric test [59] with a significance level of \(\alpha =0.05\) to ascertain the statistical significance of the differences in the mean ranks for each measure. The test rejected the null hypothesis (with \(\alpha =0.05\)) that the differences in mean rankings of MAPE, MAE, and RMSE obtained by the different algorithms were equal. Specifically, the confidence interval for this number of datasets and algorithms was \(C_0 = (0, \mathrm{F}_{(\alpha =0.05)}=2.14)\), and the corresponding F-values for each metric were \(13.59 \notin C_0\), \(9.53 \notin C_0\), and \(8.35 \notin C_0\), respectively.

According to the guidelines outlined by Demšar (2006) [58], the best-performing feature extraction method, \(\texttt {N-BEATS}_\mathrm {+TDA}\), was designated as the control method and compared to the remaining methods based on their rankings. It has been observed that comparing all forecasting models to each other in a post hoc test is less sensitive than comparing all forecasting models to a specific control method. One statistical procedure for conducting this type of comparison is the Holm test. The test statistic for comparing the i-th and j-th methods using this procedure is defined as:

$$\begin{aligned} z = \frac{\overline{R}_{i}-\overline{R}_{j}}{\sqrt{\frac{k(k+1)}{6N}}}, \end{aligned}$$
(9)

where k is the number of methods, N is the number of datasets, and \(\overline{R}_i\) is the mean ranking of the i-th method. The z value is used to determine the corresponding probability from the normal distribution table, which is then compared with an appropriate significance level, \(\alpha \). Holm’s test adjusts the value for \(\alpha \) to compensate for multiple comparisons. This is accomplished through a step-up procedure that sequentially tests the hypotheses ordered by their significance. The ordered p-values are denoted as \(p_{1}, p_{2}, \ldots , p_{k}\) such that \(p_{1} \le p_{2} \le \ldots \le p_{k-1}\). Holm’s test compares each \(p_{i}\) with \( \alpha ^{\prime}_{{{\text{Holm}}}} = \alpha /(k - i) \), starting with the most significant p-value. If \(p_{1}\) is below \(\alpha /(k - 1)\), the corresponding hypothesis is rejected, and \(p_{2}\) is compared with \(\alpha /(k - 2)\). If the second hypothesis is rejected, the test proceeds with the third, and so on. All remaining hypotheses are retained if a null hypothesis cannot be rejected.

The results of the Holm test for \(\alpha = 0.05\) are presented in Table 5, utilizing the corresponding p and \(\alpha ^{'}_{\mathrm{Holm}}\) values. Based on the outcomes of this test, it is inferred that higher rankings of MAPE, MAE, and RMSE are achieved by the \(\texttt {N-BEATS}_{\mathrm {+TDA}}\) method compared to all other methods \(\alpha = 0.05\). As can be observed in the statistical accuracy results of forecasting, high competitiveness is demonstrated by the baseline method, N-BEATS, and enhancing its prediction performance through additional features is deemed challenging, which appears to be accomplished by the proposed method through the incorporation of features derived from TDA.

Table 5 Adjusted p-values using MAPE, MAE, and RMSE as the test variables (\(\texttt {N-BEATS}_\mathrm {+TDA}\) is the control method)

5 Discussion

The results of this study highlight the effectiveness of integrating topological data analysis (TDA)-derived features into time series forecasting models, specifically the N-BEATS model. By capturing complex, non-linear dependencies through TDA features such as persistent entropy, amplitude, and the number of points in persistence diagrams, these findings confirm the utility of topological features in enhancing predictive accuracy. The significance of these results is twofold: they validate the relevance of TDA in capturing inherent structures within financial time series and suggest several promising pathways for future research to extend TDA’s applicability within the financial domain.

One potential extension for TDA-based feature engineering lies in its application to more complex modeling environments, particularly in reinforcement learning (RL) trading strategies. In recent years, RL has emerged as a pivotal tool for creating automated trading strategies within financial markets. The adaptability of RL to the dynamic and often volatile nature of financial markets underscores its potential, as it learns optimal actions through trial and error within complex environments. However, the effectiveness of RL-based strategies heavily depends on the quality and relevance of features extracted from market data. Incorporating TDA-derived features could contribute to a richer, more informative representation of market dynamics, thereby enhancing the adaptability of RL algorithms to variable market conditions. With their capacity to capture intricate structures and persistent patterns in data, TDA features could bridge the gap between the highly non-linear, high-dimensional characteristics of financial data and the decision-making framework of RL-based trading systems. This integration may improve RL model performance and contribute to innovative avenues in quantitative finance.

In addition to RL applications, TDA’s role in high-frequency trading (HFT) environments warrants further investigation. High-frequency data present unique challenges due to its high dimensionality and inherent noise, which complicate traditional modelling approaches. With its inherent ability to capture persistent structures amidst chaotic data, TDA could prove valuable in extracting meaningful patterns within this context. Thus, future studies could explore the utility of TDA in high-frequency trading datasets to evaluate its effectiveness in real-time decision-making. If successful, such applications could support high-frequency trading systems by providing enhanced interpretability of market signals in rapid trading scenarios.

Another promising area for future exploration involves combining TDA-derived features with advanced deep learning architectures beyond the N-BEATS model. Although N-BEATS was selected for its robust performance in time series forecasting and its flexibility in feature integration, further exploration could assess the benefits of integrating TDA features within models like transformers or graph neural networks. These models, particularly suited for capturing complex dependencies in non-stationary time series, may enhance predictive accuracy when combined with TDA-derived features. Additionally, applying TDA features within ensemble frameworks may broaden the spectrum of time series characteristics captured, potentially leading to improved predictive outcomes.

It is essential to acknowledge certain limitations of this study. While offering a controlled environment for assessing the model’s predictive power, the focus on univariate time series restricts its ability to account for interactions among multiple assets-a critical aspect in real-world financial markets. Extending future analyses to multivariate settings could reveal inter-asset relationships, potentially unveiling new patterns and interactions that enhance the predictive power of TDA-augmented models. Furthermore, this study’s relatively short time horizon may limit the generalizability of findings across different market cycles and phases. Expanding the dataset to cover longer periods would allow for testing model robustness across diverse market conditions, potentially offering a more comprehensive view of financial market dynamics and further stabilizing TDA-based forecasting models.

6 Conclusions

This manuscript introduces the innovative integration of topological data analysis (TDA) into financial time series forecasting, specifically aimed at enhancing the predictive accuracy of the N-BEATS model through the utilization of TDA-derived features, such as entropy, amplitude, and the number of points from persistence diagrams and point clouds, computed via a sliding window approach. The study rigorously compares traditional feature engineering techniques with TDA-derived methods across time series data from six cryptocurrencies and four traditional financial instruments, implementing pre-processing techniques to handle missing data, outlier detection, and stationarity tests. Results obtained from performance evaluations using metrics such as MAPE, MAE, and RMSE indicate that the \(\texttt {N-BEATS}_\mathrm {+TDA}\) model outperformed traditional techniques, with statistical tests affirming the significant improvements in predictive accuracy attributable to the incorporation of TDA features, thereby advancing the application of TDA in financial analysis and enhancing decision-making processes.