Forecasting in Non-stationary Environments with Fuzzy Time Series

Time series data observed in different real-world applications are often non-stationary. Given that a stationary time series is defined in terms of its mean and variance, non-stationarity can be detected if any (or both) of these components vary over time. Thus, in a non-stationary context, if the chosen forecasting model relies on a false assumption of stationarity, there is a significant risk that it will present a degrading performance over time and eventually become obsolete [7].

Some non-stationary time series are consequence of hidden contexts, not given in the form of predictive features [19]. In such situations, inferring a forecasting model becomes more difficult, since changes in the hidden context can induce unpredictable changes in the target concept. These changes are known as concept drift and since they usually cannot be identified explicitly [7], an effective forecasting model should be able to quickly adapt to the resulting changes in the time series data.

A symptom of concept drift in a time series is the change of the parameters that define its distribution [20]. Following [20], the change in the distribution can be classified with respect to the rate at which the drift occurs. For instance, a political event that suddenly causes a strong effect on the stock market is a case of abrupt concept drift observed in a time series. Another example is the aging effects of a sensor which gradually leads to lower performance in a device. Such case can be referred to as gradual concept drift. According to Ditzler et al. [7], the drifts, whether abrupt or gradual, can also be classified as permanent or transient. The former is related to the effect of variation, and is not limited in time, while the latter occurs in a limited time window followed by an effect of disappearance.

Due to the practical importance of problems with varying statistical properties, the literature has presented some alternatives to handle the issue. Popular approaches commonly used in forecasting methods, including FTS, are: (i) detrending [21, 22], in which a data transformation is applied to the time-series to remove the trend component and (ii) time-varying models derived from the distribution of the raw data [23].

Since the introduction of the Fuzzy Time Series in [8], several categories of FTS methods have been proposed, varying mainly by their order $\Omega$, number of fuzzy sets $k$ and time-variance [18] – see Table 1 for the convention of symbols adopted here. The order is defined by the number $\Omega$ of time-delays (lags) that are used in modeling the time series. The time variance defines whether the FTS model changes over time, with the Time Invariant models $\mathcal{M}^{t}$ having the same parameters for all $t=0,\ldots,T$, and the Time Variant models $\mathcal{M}^{t}$ having different parameters in different time instants $t$.

Given an univariate time series $Y\in\mathbb{R}$ , where $y(t)\in Y$ are the instances of $Y$ for $t=0,1,\ldots,T$, the UoD is limited by the known bounds of $Y$, such that $U=[\min(Y),\max(Y)]$. The training procedure of an FTS model $\mathcal{M}$ consists of the following three steps:

• a)

  Partitioning: split $U$ in $k$ overlapping intervals. For each interval a new fuzzy set $A_{j}$ is created, each one with its own membership function (MF) $\mu_{A_{j}}$. A linguistic value $\tilde{A}$ is assigned to each fuzzy set and represents a region of $U$. The computational cost of this step is $O(k)$.

• b)

  Fuzzification: maps the crisp time series $Y$ onto the fuzzified time series $F$, by replacing each $y(t)\in Y$ by the fuzzified value $f(t)=\mu_{A_{j}}(y(t))$, $\forall A_{j}\in\tilde{A}$, for $t=1,\ldots,T$. The computational cost of this step is $O(T\cdot k)$.

• c)

  Knowledge Extraction: creates a representation of the sequential patterns in the time series. In a rule-based FTS, as in [24], the rules have a format $A_{i}^{t-\Omega},\ldots,A_{i}^{t-1}\rightarrow A_{j},A_{k},\ldots$, where $A_{i}^{t-\Omega},\ldots,A_{i}^{t-1}\in\tilde{A}$ is the precedent and $A_{j},A_{k},\ldots\in\tilde{A}$ is the consequent. The rule can be read as “IF $y(t-\Omega)$ is $A_{i}^{t-\Omega}$ AND $\ldots$ AND $y(t-1)$ is $A_{i}^{t-1}$ THEN $y(t)$ may be $A_{j},A_{k},\ldots$”. The computational cost of this step is $O(T\cdot k^{\Omega})$.

Once the model $\mathcal{M}$ is trained it can be used to forecast values of $Y$ given a sample $y(t-\Omega),\ldots,y(t-1)$ with a three step procedure:

• a)

  Fuzzification: maps the crisp sample $y(t-\Omega),\ldots,y(t-1)$ onto the fuzzified values $f(t-\Omega),\ldots,f(t-1)$, where each $f(t)=\mu_{A_{j}}(y(t))$, $\forall A_{j}\in\tilde{A}$, for $t=\Omega.\ldots,1$. The computational cost of this step is $O(\Omega\cdot k)$.

• b)

  Rule Matching: find in $\mathcal{M}$ the rules whose the precedent matches with the fuzzified values $f(t-\Omega),\ldots,f(t-1)$. The activation $\mu_{j}$ of each rule $j$ is the minimum T-norm of the individual membership values of each fuzzified value. The computational cost of this step depends on the length of the sample ($\Omega$) and the number of rules in $\mathcal{M}$ ($k^{\Omega}$), then the complexity is $O(\Omega\cdot k^{\Omega})$.

• c)

  Defuzzification: the estimated value of $\hat{y}(t)$ is calculated by finding the mean value $c_{j}$ of each matched rule $j$, by averaging the midpoints of the consequent fuzzy sets, and then calculating the sum of the mean values of each rule weighted by its activation values

  $$\hat{y}(t)=\frac{\sum_{j}\mu_{j}\cdot c_{j}}{\sum_{j}\mu_{j}}$$

Many improvements were proposed in the FTS literature. The High-Order FTS (HOFTS) [12] extended the classical method in [24] by using several lags in the forecast and it is able to recognize more complex patterns in the time series. In [25], [26] and [27] weights are added in the consequent of the rules in order to give more importance for certain fuzzy sets. More recently, the Probabilistic Weighted FTS (PWFTS) was proposed in [28], and made available in the pyFTS library [29], including weights in both precedent and consequent of the FTS rules, achieving high accuracy and outperforming traditional forecasting approaches.

In order to represent forecasting uncertainty, [14], [28] and [30] proposed approaches for probabilistic forecasting. In [31] and [32] distributed FTS variants are proposed for Big Data time series. Multivariate time series are explored in [33], which uses Fuzzy Information Granules (FIG) to propose a multivariate forecasting method. In [18] a survey on the design of FTS forecasting models was provided. A comprehensive review of those aforementioned models, focused on time-invariant and rule-based approaches can be found in [32].

The major hyper-parameters of FTS methods are the number of fuzzy sets $k$ and the order of the model $\Omega$. These hyper-parameters conduce the model training and forecasting and are responsible for the accuracy and model parsimony (the number of rules). A method for hyperparameter tuning of FTS models is presented in [34].

The FTS approaches listed in this section are all time-invariant approaches which means that the models are trained only once and then its internal rule base does not change. To keep their accuracy, these models must make strong assumptions about the stationarity and homoskedasticity of the time series. In non-stationary environments, this is a major drawback.

To better understand the behavior of the FTS methods in non-stationary scenarios, Figure 1 shows the performances of several methods in the literature when the test data falls out of the known UoD. For this example, we used the NASDAQ dataset. It can be seen that this situation makes most of the trained models useless in the long run.

Non-stationary fuzzy reasoning and non-stationary fuzzy sets (NSFS) were introduced by Garibaldi and Ozen and Garibaldi, Jaroszewski and Musikasuwan, respectivelly in [38] and [16]. The main idea is to extend the traditional fuzzy set definition by introducing a dynamic component that changes the membership function $\mu$ over time. This change takes several forms: a) variation in location: by displacing the parameters of the $\mu$ function along the UoD without changing its shape; b) variation in width: changing the shape of $\mu$ by stretching or contracting its bounds; and c) noise variation: by adding random noise to the membership grade.

A NSFS is defined with two functions: the non-stationary membership function (NSMF), which considers time variations of the corresponding membership function (MF), and the perturbation function, which is the dynamic component responsible for altering the parameters of the membership function given some parameter set.

According to Garibaldi et al. [16] a non-stationary fuzzy set $\dot{A}$ can be formalized as follows:

where $\dot{A}$ is a fuzzy set over a UoD $X$ characterized by a NSMF $\mu_{\dot{A}}(t,x)$ at a set of time points $T$.

It is worth noting that the NSMF $\mu_{\dot{A}}$ varies throughout $U$ and the time interval $T$. This means that the NSMF parameter set should vary over time. A regular MF, $\mu_{A}(x)$, can be expressed as $\mu_{A}(x,p_{1},\ldots,p_{m})$, where $p_{1},\ldots,p_{m}$ denote the parameters of $\mu_{A}(x)$. Thus, a NSMF can be denoted in an analogous way as follows:

where each parameter can be varied over time by a perturbation function multiplied by a constant.

One of the constraints existing in several FTS models is the lack of ability to deal with conditional variances and concept-drift scenarios, in which the statistical characteristics of the time series change, sometimes drastically. Thus, through small variation in the MFs, we deploy NSFS in FTS models to deal with variability for decisions over time and contributing to the mitigation of this drawback.

Given the training data, $Y$, the number of partitions, $k$, and the length of the residuals window, $w$:

1. Step 1

   Defining the universe of discourse, $U$:

   $$U=[lb,ub]$$

   (6)

   where, $lb=\min(Y)-\min(Y)\times 0.2$ and $ub=\max(Y)+\max(Y)\times 0.2$.

   Notice that the data bounds are extrapolated to compensate for a possible underestimation of the bounds in the training set. The value $0.2$ is a typical value, but can be modified according to the problem.

2. Step 2

   $U$ Partitioning: Define the midpoints $c_{i},i=0,...,k-1$ of the initial fuzzy sets;

   $$c_{i}=lb+i\times\frac{(ub-lb)}{(k-1)}$$

   (7)

3. Step 3

   Define the linguistic variable $\tilde{A}$: Create $k$ overlapping fuzzy sets $A_{i}$, with triangular MF $\mu_{A_{i}}$ as defined in equation (8).

   $$\mu_{A_{i}}(x)=\left\{\begin{array}[]{lcr}0&if&x<l\;or\;x>u\\ \frac{x-l_{i}}{c_{i}-l_{i}}&if&l_{i}\leq x\leq c_{i}\\ \frac{u_{i}-x}{u_{i}-c_{i}}&if&c_{i}\leq x\leq u_{i}\end{array}\right.$$

   (8)

   where $l_{i}=c_{i-1}$ and $u_{i}=c_{i+1}$.

   Each fuzzy set $A_{i}\in\tilde{A}$ is a linguistic term of the linguistic variable $\tilde{A}$. Once the fuzzy sets are created, a function $\pi_{i}$, as defined in equation (5), will be associated with each fuzzy set in order to transform it into a NSFS, initialized with $\delta=0$ and $\rho=0$.

   The number of sets $k$ defines the number of states and consequently the number of state transitions that the model can represent. The more complex the time series the greater the number of $k$ should be. One should, however, be careful to not overestimate $k$ since it may cause overfitting and make the model unnecessarily big.

4. Step 4

   Fuzzification: Transform the original time series data $Y=\{y(0),y(1),\ldots,y(T)\}$ into a fuzzy time series $F=\{f(0),f(1),\ldots,f(T)\}$, where $f(t)$ is defined as:

   $$f(t)=\{\mu_{A_{0}}(y(t)),\mu_{A_{1}}(y(t)),...,\mu_{A_{k-1}}(y(t))\}$$

   (9)

5. Step 5

   Generate the temporal patterns: The fuzzy temporal patterns have format $A_{l}\rightarrow A_{r}$, where:

   The precedent (or Left Hand Side - LHS) is:

   $$A_{l}=\underset{A_{i}}{\arg\max}\left(\mu_{A_{i}}(y(t-1))\right)$$

   (10)

   And the consequent (or Right Hand Side - RHS) is:

   $$A_{r}=\underset{A_{i}}{\arg\max}\left(\mu_{A_{i}}(y(t))\right)$$

   (11)

6. Step 6

   Generate the rule base: Select all temporal patterns with the same precedent and group their consequent sets creating a rule with the format $A_{l}\rightarrow A_{a},A_{b},...$. Thus the rule $RHS$ can be understood as the set of possibilities which may happen on time $t+1$ (the consequent) when a certain set $A_{l}$ is identified on time $t$ (the precedent).

7. Step 7

   Compute the residuals $\epsilon$: Invoke the Forecasting Procedure defined on Section 4.3 using the given training data, $Y=\{y(0),...,y(T)\}$, as the input and forecast the last $w$ items to calculate the set of residuals $\mathcal{E}$ defined as:

   $$\mathcal{E}=\{\epsilon(t-w),\epsilon(t-(w-1)),...,\epsilon(t)\}$$

   (12)

   where $\epsilon(t)=y(t)-\hat{y}(t)$ and $\hat{y}(t)$ is the forecast produced by the model.

Given an input value $y(t)$, the last forecast value $\hat{y}(t)$, the length of the residuals vector $w$, the residuals set $\mathcal{E}$ and the linguistic variable $\tilde{A}$, do:

1. Step 1

   Find the displacements of $y(t)$ on $U$: If $y(t)$ is below the lower bound of $U$ store the difference in $d_{l}=lb-y(t)$, otherwise $d_{l}=0$. If $y(t)$ is above the upper bound of $U$, store the difference in $d_{u}=y(t)-ub$, otherwise $d_{u}=0$. Then compute the displacement range $r$ as the sum of the UoD displacements $d_{u}$ and $d_{l}$, according with (13), and the displacement midpoint $mp_{r}$ as $r$ divided by 2, according with (14)

   $$r=d_{u}-d_{l}$$

   (13)

   $$mp_{r}=r/2$$

   (14)

2. Step 2

   Compute the mean $\overline{\mathcal{E}}$ and variance $\sigma_{\mathcal{E}}$ of the set $\mathcal{E}$: The residual mean $\overline{\mathcal{E}}$ indicates a bias, a shift on the accuracy of the trained model that must be corrected. The variance, $\sigma_{\mathcal{E}}$, represents the shift on the trained model variance, as a reflection to a change in the test data. These values will be used to adjust the position and length of the fuzzy sets.

3. Step 3

   Compute the displacements, $\delta_{i}$: The displacement is computed for each fuzzy set $A_{i}$ in order to equally move the $k$ fuzzy sets by the mean shift $\overline{\mathcal{E}}$, one fraction of the displacement range $[d_{l},d_{u}]$ and proportionally to the expansion of the variance in the interval $[-\sigma_{\mathcal{E}},\sigma_{\mathcal{E}}]$:

   $$\delta_{i}=\overline{\mathcal{E}}+\left(i\frac{r}{k-1}-mp_{r}\right)+\left(i\frac{2\sigma_{\mathcal{E}}}{k-1}-\sigma_{\mathcal{E}}\right)$$

   (15)

   The displacements $\delta_{i}$, for $i=0\ldots k-1$, are equally distributed in the interval $[\overline{\mathcal{E}}-d_{l}-\sigma_{\mathcal{E}}\;,\;\overline{\mathcal{E}}+d_{u}+\sigma_{\mathcal{E}}]$ and indicate the new position of fuzzy set $A_{i}$, given the deviations from the known UoD, whose range is $r$ and its midpoint $mp_{r}$, and the error signal with the mean $\overline{\mathcal{E}}$ and variance $\sigma_{\mathcal{E}}$. Indeed, while the term $i(r/(k-1))-mp_{r}$ is used to translate the midpoint of $A_{i}$ by a fraction of $r$ centered in $mp_{r}$, the term $i(2\sigma_{\mathcal{E}}/(k-1))-\sigma_{\mathcal{E}}$ is used to offset the scaling of the fuzzy set bounds by a fraction of $\sigma_{\mathcal{E}}$ centered in $0$.

4. Step 4

   Compute the scaling factor $\rho_{i}$, Eq. (16): For each fuzzy set $A_{i}\in\tilde{A}$ the scale increment is empirically calculated as the distance between the displacements $\delta_{i}$, in order to adjust the fuzzy set lengths proportionally to their displacements, avoiding discontinuities between the sets (intervals on $U$ not covered by any fuzzy set) by setting $l_{i}=c_{i-1}$ and $u_{i}=c_{i+1}$:

   $$\rho_{i}=|\delta_{i-1}-\delta_{i+1}|$$

   (16)

5. Step 5

   Update $\epsilon$: Get the last forecast value $\hat{y}(t+1)$ and the last known value $y(t+1)\in Y$. Compute the error term $\epsilon(t+1)=y(t+1)-\hat{y}(t+1)$, push it on to the residuals set $\mathcal{E}$ and delete the oldest value. That is:

   $$\mathcal{E}=\mathcal{E}\setminus\epsilon(t-w)\cup\epsilon(t+1)\$$

   (17)

Given an input value $y(t)$, the linguistic variable $\tilde{A}$, the inferred rule set and the perturbation parameters $\delta_{i}$ and $\rho_{i}$ for each fuzzy set $A_{i}\in\tilde{A}$:

1. Step 1

   Fuzzification: Compute the membership grade $\mu_{A_{i}}$ for each non-stationary fuzzy set $A_{i}$, such that $\mu_{A_{i}}=\mu_{A_{i}}(y(t),\pi(l_{i},c_{i},u_{i},\delta_{i},\rho_{i}))$;

2. Step 2

   Rule matching: Build a set, $\mathcal{S}$ with all the rules $A_{j}\rightarrow RHS_{j}$ where $\mu_{A_{j}}(y(t))>0$, that is:

   $$\mathcal{S}=\{A_{j}\rightarrow RHS_{j}|\mu_{A_{j}}(y(t))>0\}$$

   (18)

3. Step 3

   Defuzzification: Compute the forecast $\hat{y}(t+1)$ according to Equation (19) as the weighted sum of the rule mid-points, $mp$, by their membership grades $\mu_{j}$ for each selected rule $j$:

   $$\hat{y}(t+1)=\sum_{A_{j}\rightarrow RHS_{j}\in\mathcal{S}}\mu_{A_{j}}(y(t))\cdot mp(RHS_{j})$$

   (19)

   where

   $$mp(RHS)=\dfrac{\sum_{A_{i}\in RHS}c_{A_{i}}}{|RHS|}$$

   (20)

The NSFS has the ability to dynamically adapt membership functions as the statistical properties of the time series vary. The NSFTS method tries to detect these changes by using the displacements of the input data $y(t)$ in relation to the known $U$ and the residuals statistics to identify bias and variance shifts.

The major assumption of the method is that adjusting the fuzzy sets, without modifying the model structure represented by the learned rules, is enough to adapt the model to the new behavior of the time series. Such procedures, described in subsection 4.2, given an input of size $T$, memory window length $W$, refreshing interval $R$ and $k$ fuzzy sets, have a time complexity of $O(T/R\cdot W\cdot\log k)$, while a retraining process would take $O(T/R\cdot W\cdot(\log k)^{\Omega})$, as discussed in Section 3. Hence, this approach becomes much cheaper computationally, when compared to retraining a new model from scratch.

The $\pi$ function entails an additive approach for the translation and scaling of the parameters based on grid partitioning of the UoD. The procedure calculates the translation increment proportionally to the displacement of the input value in relation to the known $U$ and the variance of the residuals. The scaling increment works as a heuristic to keep the grid aspect of the fuzzy sets after the translation, connecting the lower and upper bounds with the centers of adjacent fuzzy sets.

Figure 5 depicts the forecasting method applied to a test sample with significant drift from the training data extracted from SP500 time series, later explained in Section 5.

In Figure 5a the partitioning of the UoD is shown with 15 partitions. The generated forecasts are presented in Figure 5b and their residuals in Figure 5c. The perturbations on the NSFS, in response to the previous residuals, are represented in Figure 5d, where the same fuzzy sets represented in Figure 5a are now colored over the vertical axis. It is possible to see that these fuzzy sets move, sometimes expanding sometimes contracting, in order to fit to the data. As the drift increases, the displacement and scaling of the fuzzy sets also grow.

The proposed model is a case of active learning, in which the learning algorithm is able to interactively obtain, from different sources of information, the necessary inputs for the generation or improvement of its learning. Basically, it uses residual error values to fit a predefined model. The version here presented focuses on highlighting its adaptability to changes observed over time. It uses first-order fuzzy rules, that is, considering only a previous observation, and is applied to forecasting problems in univariate time series. However, its time variant model characteristics can be expanded or adapted to other models for better prediction values. For instance, the model can be expanded to comprehend a high-order fuzzy rule generation. It also can be adapted to fit parameters of other FTS models, such as PWFTS, used as benchmark in this work.

Different datasets were chosen for model validation. The datasets consisted of four stock market indices (Dow Jones, NASDAQ, SP500 and TAIEX), three FOREX pairs (EUR-USD, EUR-GBP, GBP-USD), two cryptocoins exchange rates (Bitcoin-USD and Ethereum-USD) illustrated in A (Fig. 6) and eight synthetic time series with concept drifts, see B (Fig. 7).

The market indexes data sets contain the daily averaged index by business day, such that the Dow Jones Industrial Average (Dow Jones)¹¹1https://finance.yahoo.com/quote/%5EGSPC/history?p=%5EGSPC. Accessed in 11/11/2019 is sampled from 1985 to 2017 time window, the Taiwan Stock Exchange Capitalization Weighted Stock Index (TAIEX)²²2https://www.twse.com.tw/en/page/products/indices/series.html. Accessed in 11/11/2019 is sampled from 1995 to 2014, the National Association of Securities Dealers Automated Quotations - Composite Index (NASDAQ ÎXIC)³³3https://www.nasdaq.com/market-activity/index/ixic. Accessed in 12/11/2019 is sampled from 2000 to 2016 and the SP500 - Standard & Poor’s 500⁴⁴4https://finance.yahoo.com/quote/%5EGSPC/history?p=%5EGSPC. Accessed in 12/11/2019 is sampled from 1950 to 2017. The FOREX data sets contain the daily averaged quotations, by business day, from 2016 to 2018, which pairs are the US Dollar to Euro (USD-EUR), Euro to Great British Pound (EUR-GBP) and Great British Pound to US Dollar (GBP-USD). The cryptocoin datasets contain the daily quotations, in US Dollar, of the Bitcoin (BTC-USD) and Ethereum (ETC-USD). The synthetic data aims to represent different types of concept drifts.

The Augmented Dickey-Fuller (ADF) test was used to check which time series are non-stationary. The Levene’s test was used to test whether the series are heteroskedastic or not. The detailed results of these tests are shown in C and D, respectively. In summary, except for the data sets, “Stationary Signal”, “Stationary Signal with blip” and “Sudden variance” (see B (Fig. 7)) all the benchmarks were shown to be non-stationary and only the “Stationary Signal with blip” homogeneity of variances with a confidence level of 95%.

The standard accuracy metrics used to evaluate point forecasting methods are the Root Mean Squared Error (RMSE) (21), Mean Absolute Percentage Error (MAPE) (22) and Theil’s U statistic (U) (23) where $y$ are the target values, $\hat{y}$ are the forecast values and $n$ the sample size. These metrics are commonly used in evaluating forecasting models [18].

The proposed NSFTS was tested against the Time Variant [39] and the Incremental Ensemble approaches, both using the PWFTS method [28] available in [29] as internal method. As can be seen in [32], the PWFTS outperformed a wide number of forecasting methods ranging from classic statistical ones such as ARIMA [45] to more recent ones such as the k-Nearest Neighbors with Kernel Density Estimation [46]. The parameters of all methods are presented in Table 2, and they were defined through grid search.

The grid partitioning scheme was used for the initial generation of fuzzy sets, where the best number of partitions in the range $[5,100]$ was selected for each data set. The high order methods were tested with orders 2 and 3.

The average RMSE, MAPE and U statistic results by method and data set are presented in Table 3. It can be seen that the NSFTS is superior or, at least, not worse than the other methods in all the selected real world benchmarks.

From the results on the artificial data sets, it can be seen that the NSFTS has more difficulties than the other methods in handling sudden changes in the series distribution mean. This behavior is expected from the NSFTS once it is an adaptive method and therefore presents a delay to respond to changes. On the other hand, it handles with relative ease incremental changes in the mean. Observing the Theil’s U statistic which is insensitive to the range of values of the series, one can see that changes in the variance alone do not seem to affect the NSFTS performance by much.

These observations explain the good performance of the NSFTS on the economic series. As one can see in A, these series present a gradual, though significant, variation of the mean along with different types of changes in variance. Apparently, the superior performance of the NSFTS when the mean varies gradually gave it the edge over the other methods.

To assess the overall performance of the methods with respect to the RMSE we ran a Friedman Aligned Ranks (F-test), a non-parametric test for the equality of the means was used with $\alpha=0.05$, where $H_{0}$ indicates that there is no significant difference between the means and $H_{1}$ indicates that there is at least one significant difference between the means. The F-value result was $13.1994$ with a p-value of $0.0013$. Thus, $H_{0}$ was rejected at 5% confidence level.

A one-versus-all Finner paired post-hoc test was employed to check the method equivalence with NSFTS as control method, where $H_{0}$ indicates that there is no significant difference between the means and $H_{1}$ indicates that there is a significant difference between the means. The results are presented in Table 4 which shows that NSFTS is not equivalent to the Incremental Ensemble method and equivalent to the time variant FTS model, considering a 95% confidence level.

It is important to remind that the NSFTS is computationally cheaper than the Time Variant method since it does not require retraining. Therefore, even though they present equivalent RMSE performance, the NSFTS still has the computational edge.