Evaluating Volatility Using an ANFIS Model for Financial Time Series Prediction

Abstract

This paper develops and implements an Autoregressive Integrated Moving Average model with an Adaptive Neuro-Fuzzy Inference System (ARIMA-ANFIS) for BTCUSD price prediction and risk assessment. The goal of these forecasts is to identify patterns from past data and achieve an understanding of the future behavior of the price and its volatility. The proposed ARIMA-ANFIS model is compared with a benchmark ARIMA-GARCH model. To evaluated the adequacy of the models in terms of risk assessment, we compare the confidence intervals of the price and accuracy measures for the testing sample. Additionally, we implement the diebold and Mariano test to compare the accuracy of the two volatility forecasts. The results revealed that each volatility model focuses on different aspects of the data dynamics. The ANFIS model, while effective in certain scenarios, may expose one to unexpected risks due to its underestimation of volatility during turbulent periods. On the other hand, the GARCH(1,1) model, by producing higher volatility estimates, may lead to excessive caution, potentially reducing returns.

1. Introduction

Forecasting time series presents a significant challenge due to the inherent complexity and unpredictability of dynamic data. The primary difficulty lies in the fact that these data are generated via an unknown process, often perceived as random. Despite this, researchers strive to approximate this elusive data generation mechanism by analyzing observed data and patterns and applying a variety of models (Hyndman and Athanasopoulos 2018).

In time series analysis, there are two main features to model: the conditional mean and the conditional variance, also known as volatility. A model for the conditional mean allows us to explain and predict the behavior of the time series, while a model for the volatility enables the evaluation of the risk associated with the mean prediction, providing a confidence interval. This ability to forecast volatility is crucial for decision-making (Hamilton 2020).

The concept of time series in a probabilistic framework was first introduced by the Scottish statistician George Yule at the beginning of the 20th century. He defined a time series as a realization of a stochastic process, which is a set of random variables $Y=\{Y_t:t\in T\}$, where $Y_t$ is a random variable and T is the index set. Generally, the index t is interpreted as time, and $Y_t$ represents the state of the process at time t, with t typically being an integer.

In the formal approach to modeling a time series, a vector of lagged variables $X_{t-1}={(Y_{t-1},Y_{t-2},\cdots ,Y_{t-p})}^T$, $t>p$, is used from the process $\left\{Y_t\right\}$, and a regression given by Equation (1) is employed. That is,

$$Y_t=f(X_{t-1},\theta )+{\epsilon }_t,$$

(1)

where ${\epsilon }_t$ is a white noise process. This regression defines a model for the conditional mean and variance of the process given the available information up to time $t-1$. A  general specification for a join model is as follows:

$$\mathrm{E}\left(Y_t\right|x_{t-1})=f(X_{t-1},\theta )$$

(2)

$$\mathrm{Var}\left(Y_t\right|x_{t-1})=\mathrm{Var}\left({\epsilon }_t\right)={\sigma }_t^2.$$

(3)

Here, $\theta$ is a vector of parameters that must be estimated based on certain criteria. Different functional forms of f give rise to various models for the process under study; additionally, different models can be used to represent the conditional variance.

The Autoregressive Integrated Moving Average (ARIMA) models are the most commonly used to model the conditional mean, also known as Box–Jenkins methodology (Hyndman and Athanasopoulos 2018; Wei 2006). These models use variations and regressions in the data to find patterns and make predictions for future values. The Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model is a statistical model that is mostly used to analyze and forecast the volatility of time series data, particularly in financial markets; see for instance Bollerslev (2023). It extends the Autoregressive Conditional Heteroskedasticity (ARCH) model by incorporating lagged values of both the variance and the squared residuals, providing a more flexible and comprehensive framework for modeling the volatility clustering observed in financial time series (Engle and Patton 2007; Poon and Granger 2003).

Alternative models are often proposed to predict the conditional mean, conditional variance, or both measures jointly. For instance, to model the volatility, attractive options include the Takagi–Sugeno (T-S) fuzzy system, also known as the Sugeno fuzzy model, a type of fuzzy inference system developed by Takagi and Sugeno in 1985 (Takagi and Sugeno 1985) that has been applied to time series analysis and forecasting. This model is widely used in various applications due to its effectiveness in handling nonlinear systems and its ability to produce smooth control actions (Alenezy et al. 2023; Tsai et al. 2019).

Recent research with applications in financial data modeling have explored joint models for both conditional mean and variance. These include hybrid models that combine traditional statistical methods with machine learning and neural networks. Goodell et al. (2021) investigated the application of artificial intelligence and machine learning in financial time series analysis, identifying main areas such as portfolio optimization and investment, fraud and financial distress detection, forecasting, and financial planning. These authors highlight the transformative impact of these technologies and suggest some future research directions. Wang et al. (2013) proposed a hybrid ARIMA-GARCH model with a neural network component to capture non-linear relationships and improve forecast accuracy. These hybrid models leverage the strengths of both traditional and modern approaches, providing more robust and accurate predictions for financial time series.

In this work, we propose an Autoregressive Integrated Moving Average model with an Adaptive Neuro-Fuzzy Inference System (ARIMA-ANFIS) for stock price predictions and risk assessment. We also estimate an ARIMA-GARCH model as a benchmark to compare the performance of the proposed model. The Diebold–Mariano test, as proposed by Diebold and Mariano (1995), is applied to evaluate the predictive accuracy between the two models. The ARIMA models have a vast range of applications, and their properties are well known. Furthermore, several studies have demonstrated that ANFIS models can be successfully used for time series modeling due to their high flexibility and ability to model nonlinear dynamics. The ANFIS models are effective in modeling and predicting complex variables in various applications such as stock index, ${\mathrm{CO}}_2$ emissions, global temperature, the COVID-19 pandemic, and customer satisfaction; for more details, see Huarng and Yu (2005); Jiang et al. (2024); Jithendra and Sharief Basha (2023); Khan and Khan (2019). Although the specification process of the ANFIS model is not fully available in an analytical way, assessing different scenarios of this model with empirical applications allows us to make meaningful inferences.

The motivation for developing and implementing a joint ARIMA-ANFIS model is to forecast the daily closing price of BTC/USD and assess its associated risk in the stock market. The primary goal of these forecasts is to identify patterns from past data and gain an understanding of the future behavior of the price and its volatility. It is clear that to make profits in the context of trading, it is necessary to evaluate the risk of an asset and buy/sell an asset at a given price and close the trade at a higher/lower price. Therefore, having reliable forecasts that increase the probability of predicting price movements is key to maximizing profits. While forecasting stock prices is challenging due to the numerous factors influencing market behavior, it remains possible to identify patterns that provide valuable insights into future trends.

The objectives of this work include: understanding and analyzing ARIMA family models and ANFIS, carrying out the specification process for the ARIMA-ANFIS and ARIMA-GARCH models step by step, and applying them to real financial time series. We describe the ARIMA, GARCH, and ANFIS models, then we propose a join estimation process for the ARIMA-ANFIS and ARIMA-GARCH models, applied to real time series, followed by a comparative analysis of both approaches.

This paper is divided into five sections. Section 2 discusses several models used in this article and their properties. Section 3 outlines the steps for model formulation. The application and forecasting results are presented in Section 4. Finally, Section 5 concludes this paper.

2. Preliminaries

In this section, we provide a concise overview of the ARIMA, GARCH, and ANFIS models, along with an introduction to fuzzy logic, which forms the foundation for the ANFIS framework.

2.1. ARIMA Models

The ARIMA models are a class of statistical models, proposed by George Box and Gwilym Jenkins in the 1970s, that revolutionized the analysis and forecasting of time series by focusing on modeling the conditional mean of a time series. Known as the Box–Jenkins methodology, these models belong to the ARIMA family and are renowned for their ability to provide accurate forecasts in univariate time series. The Box–Jenkins methodology consists of four iterative steps: identification, parameter estimation, diagnostic checking, and forecasting. In the identification phase, it is of utmost importance to transform the original series into a stationary series, as stationarity is the fundamental condition for effectively constructing an ARIMA model (Hyndman and Athanasopoulos 2018).

The ARIMA models are recognized for their robustness and efficiency in forecasting financial time series, especially for short-term predictions, often outperforming the most popular artificial neural network techniques (Khashei and Bijari 2011). These models have been widely used in economics and finance. Various studies have employed ARIMA models for forecasting, as mentioned previously. ARIMA modeling is essentially an exploratory, data-oriented approach that allows for fitting an appropriate model based on the structure of the data. By using autocorrelation and partial autocorrelation functions, it is possible to model the stochastic nature of the time series. This facilitates the identification of trends, random variations, periodic components, cyclic patterns, and serial correlations. As a result, forecasts of future values of the series can be obtained with a certain degree of accuracy.

The ARIMA($p,d,q$) model can be written as follows:

$${▿}^d{\dot{Y}}_t=\alpha +{\phi }_1{▿}^d{\dot{Y}}_{t-1}+{\phi }_2{▿}^d{\dot{Y}}_{t-2}+\dots +{\phi }_p{▿}^d{\dot{Y}}_{t-p}+{\theta }_1{\epsilon }_{t-1}+{\theta }_2{\epsilon }_{t-2}+\dots +{\theta }_q{\epsilon }_{t-q}+{\epsilon }_t,$$

which can also be expressed as follows:

$${\mathsf{\Phi }}_p\left(B\right){(1-B)}^d{\dot{Y}}_t={\Theta }_q\left(B\right){\epsilon }_t,$$

where ${▿}^d={(1-B)}^d$, p is the number of lags in the model, q is the size of the moving average window, and d corresponds to the unit roots of the process; the number of times the time series must be differenced to achieve stationarity. The model parameters ${\phi }_1,\dots ,{\phi }_p$, ${\theta }_1,\dots ,{\theta }_q$ must satisfy certain conditions to ensure the stationarity and invertibility of the process. $\alpha$ is a constant, and $\left\{{\epsilon }_t\right\}$ represents a sequence of identically distributed independent random variables with zero mean and constant variance, commonly referred to as a white noise process. For example, ARIMA(0,0,0) represents white noise, and ARIMA(0,1,0) with/without a constant corresponds to a random walk with/without drift.

The selection of the p, d, and q values is part of the model identification process, which is conducted using statistical criteria, the autocorrelation and partial autocorrelation functions, and statistical tests such as the Dickey–Fuller test for unit roots (Dickey and Fuller 1979).

2.2. GARCH Models

One of the most used models for statistical modeling and forecasting the conditional volatility is the generalized autoregressive conditional heteroskedastic (GARCH) approach of Bollerslev (2023). In a GARCH($p,q$) model, the current variance, ${\sigma }_t^2$, is a function of the past squared shocks, $\{a_{t-i}^2;i=1,\cdots ,p\}$, and the past variances, $\{{\sigma }_{t-j}^2;j=1,\cdots ,q\}$:

$${\sigma }_t^2=\kappa +\sum _{i=1}^p{\alpha }_i{a}_{t-i}^2+\sum _{j=1}^q{\beta }_j{\sigma }_{t-j}^2$$

with

$$a_t={\sigma }_t{\epsilon }_t\mathrm{and}{\epsilon }_t\sim N(0,1)$$

where the parameters $\kappa$, ${\alpha }_i$, and ${\beta }_j$ are subject to the following restrictions: $\kappa >0$, ${\alpha }_i\ge 0$, ${\beta }_j\ge 0$, and ${\sum }_{i=1}^{max(p,q)}({\alpha }_i+{\beta }_i)<1$.

2.3. Fuzzy Logic

Fuzzy logic, introduced by Zadeh in 1965 (Zadeh 1965), provides a mathematical framework for representing and managing uncertainty and vagueness. It offers formal tools such as fuzzy sets, if–then rules, and fuzzy arithmetic. In fuzzy logic, the degree of fuzzy membership signifies the similarity between events, where the exact properties of these events are not precisely defined. Membership degrees in fuzzy logic range continuously between 0 (completely false) and 1 (completely true). Zadeh’s pioneering work on fuzzy sets laid the foundation of fuzzy set theory. The core concept is that an element can belong to a set with a degree of membership, making propositions not strictly true or false but instead partially true or false. Following fuzzification of values, linguistic rules are applied to derive outputs, which may either retain their fuzzy nature or be defuzzified to yield a discrete (numerical) value.

Working with fuzzy sets begins with defining membership function and establishing the domain of the problem, which involves selecting appropriate functions to represent it.

Linguistic rules are employed to connect inputs with outputs. A fuzzy rule is symbolically represented as follows:

IF <fuzzy statement> THEN <fuzzy statement>

where <fuzzy statement> could be an expression in natural language.

In fuzzy systems featuring fuzzy premises, all rules are partially activated, and the consequent is true to a certain extent. The process used to compute an output value for a given input is known as fuzzy inference. This process distinguishes between two primary types: the Mamdani model, proposed by Mamdani and Assilian (Mamdani and Assilian 1999) in 1975, and the TSK model (Takagi, Sugeno, and Kang), introduced by Takagi and Sugeno (1985) as an alternative to the Mamdani model. In the present study, the TSK model is employed.

In Takagi–Sugeno fuzzy systems, the rules are structured as follows:

$$\mathrm{IF}x\in A\mathrm{and}y\in B\mathrm{THEN}z=f(x,y),$$

where x and y represent the input variables, and $A,B$ are fuzzy sets associated with membership functions. These membership functions can take various functional forms chosen from a wide range of options.

2.4. Adaptive Neuro Fuzzy Inference Systems—ANFIS

Adaptive Neuro-Fuzzy Inference Systems (ANFISs) are computational techniques from the domain of soft computing. Soft computing techniques are designed to model and handle uncertainty, approximation, and imprecision in problem-solving. Unlike traditional (hard) computing, which relies on binary logic and crisp values, soft computing incorporates methodologies that allow for flexibility and tolerance for imprecision, making it suitable for complex, real-world problems where exact solutions are difficult or impossible to find.

ANFISs are recognized for their adaptability and effectiveness in modeling complex relationships, especially in nonlinear systems and time series prediction (Walia et al. 2015). These models have found wide applications in various fields due to their ability to accurately capture intricate patterns in data. Similar to ARIMA models, ANFIS modeling involves an exploratory approach that utilizes data-driven techniques to tailor the model to the specific characteristics of the dataset. By integrating fuzzy logic principles with neural network architectures, ANFIS models can effectively handle the uncertainties and nonlinearity present in time series data (Talebizadeh and Moridnejad 2011). This capability allows ANFIS models to uncover trends, periodic components, and other complex patterns that influence the behavior of financial markets.

The Takagi–Sugeno-type fuzzy systems from Takagi and Sugeno (1985) have been used in modeling time series and predicting the mean and volatility (Sahiner et al. 2023; Venugopal et al. 2024). These models have demonstrated the ability to provide accurate forecasts and address the challenges inherent in volatility prediction. ANFIS is a type of Takagi–Sugeno fuzzy system that combines the fuzzy logic principles of Takagi–Sugeno models with neural network structures. These model combine the interpretability of fuzzy systems with the learning capabilities of neural networks.

The fundamental characteristic of ANFIS models is their ability to partition each input variable into two or more regions. The structure of ANFIS models is illustrated in Figure 1, where X and Y are the independent variables, and z represents the defuzzification. The domain of X is divided into regions $A_1$ and $A_2$, and the domain of Y is divided into regions $B_1$ and $B_2$. Consequently, the domain of the system, the $xy$ plane, is divided into regions $1,2,3$, and 4. Each of these regions have been assigned a linear model of the form $z=ax+by+c$ as the consequent function.

Finally, the model is represented by a set of rules, where the antecedent determines the region that the point to be evaluated belongs to, while the consequent corresponds to the linear model. One of the advantages of ANFIS models is that a point can simultaneously belong to two or more regions. Consequently, the value of z is calculated considering the linear models of each region. As an example, this area of ambiguity is shown in gray in Figure 1. In this case, the membership of a point to a region is determined through a fuzzy set, which has its membership function, ${\mu }_i\left(x\right)$, indicating the degree to which it is associated with the fuzzy set $S_i$, as mentioned in the previous section. The training process of an ANFIS involves the estimation of premise and consequence parameters using an optimization algorithm. The choice of optimization method is crucial for achieving optimal results, as highlighted by Karaboga and Kaya (2019).

For this example, the inference process performed to calculate the value of z given an input is as follows:

• Calculate the membership functions ${\mu }_1,{\mu }_2,{\mu }_3,{\mu }_4$ for each point in the $xy$ plane.
• For each rule, estimate $w_j$, which is defined as the firing strength of each rule. It can be calculated as a weight of the membership values or by multiplying the membership function values.
• Establish the percentage contribution of each rule to the final solution: ${\overline{w}}_j={\displaystyle \frac{w_j}{{\sum }_{i=1}^4{w}_i}}$.
• Finally, calculate the system’s output as ${\widehat{f}}_t={\sum }_i(w_i{f}_i$), where $f_i$ is typically a linear combination of the variables in the consequent.

The volatility models play a fundamental role in financial decision-making by providing insights into the uncertainty and potential future movements of asset prices. They must be able to quantify and forecast the variability associated with the returns of a financial time series. Despite being a topic of interest for many researchers (Poon and Granger 2003), there are several challenges in obtaining accurate volatility forecasts. These challenges include the fact that volatility is a non-observable feature, its estimator has a changing variance over time, and it exhibits clusters of similar variances, heavy-tailed distributions, and non-linear and non-stationary behavior, among other complexities (Bollerslev and Engle 1993; Poon and Granger 2003). The ability of ANFIS to model complex behaviors and nonlinear systems, especially for predicting volatility, is indeed a promising tool. ANFIS provides a sophisticated approach for predicting volatility by leveraging the strengths of neural networks and fuzzy logic.

3. Specification of ARIMA-ANFIS and ARIMA-GARCH Models

The ARIMA-ANFIS and ARIMA-GARCH models are estimated to predict the time series and evaluate the quality of this prediction by computing the confidence intervals using the volatility estimated from the ANFIS and GARCH models.

First, we model the conditional mean using the ARIMA model and obtain the residuals. Then, we check the correlation in the squared residuals to determine the appropriateness of applying the ANFIS and GARCH model for capturing the conditional variance of the time series.

3.1. Identification Process for the Conditional Mean Model

We perform an automatic model identification process to determine which model best explains the BTC/USD price. For the selected model, we perform validation tests for the residuals, a normality test, and finally, we determine the independence of the residual using the Ljung–Box test.

The forecast values of the time series are computed using the library forecast (from R package) to generate predictions from the fitted ARIMA model.

The framework for the training and testing process is described as follows: initially, the identification process and the model estimation are made using the training sample, which comprises 80% of the total dataset. Subsequently, the model is tested using the remaining 20% of the data, which is called the testing sample, by making predictions for 1 day ahead. This testing strategy employs an expanding window strategy, where the model is retrained at each step following the cross-validation (CV) guidelines for time series, as outlined by Bergmeir et al. (2018). For further references, please see Bergmeir et al. (2018).

3.2. Volatility with ANFIS Model

Let $\left\{{\epsilon }_t\right\}$ be the time series of the residuals from the ARIMA model in the training sample. The ANFIS model consist of L fuzzy rules, written as follows:

$$\mathrm{IF}{x}_t\in S_i\mathrm{THEN}{f}_i=f({\theta }_i,x_t).$$

Here, $x_t=[{\epsilon }_{t-1}^2,{\epsilon }_{t-2}^2]$, with ${\epsilon }_{t-1}$ and ${\epsilon }_{t-2}$ being the lagged values of ${\epsilon }_t$; $S_i$ is a fuzzy set in the input space; $f_i$ is the forecasted value for ${\epsilon }_t^2$ using the i-th rule; and $f(·,·)$ is a function that takes the following form:

$$f_i\left(x_t\right)=a_i{\epsilon }_{t-1}^2+b_i{\epsilon }_{t-2}^2+c_i,i=1,\dots ,L,$$

(4)

where $a_i,b_i$ and $c_i$ are parameters to be determined.

Note that each region into which the domain is divided is assigned a model of the form (4). It is necessary to define the membership function of $x_t$ to the set $S_i$. In ANFIS, each set $S_i$ is represented by its center $m_i$, and the value of the membership function, ${\mu }_i\left(x_t\right)$, is defined as a function of the distances from the point $x_t$ to the center of the clusters or fuzzy sets.

As an example, Figure 2 presents the architecture of a model for an ANFIS with four fuzzy rules. For a 3D representation, refer to Aznarte and Benítez (2010), which illustrates the division of the plane into four fuzzy sets, along with their corresponding membership functions and the associated planes for each set. For a point in the plane, we obtain membership values for all the fuzzy sets and the defuzzified values given in Equation (4).

This ANFIS model consist of four layers: fuzzification, rule evaluation, aggregation, and defuzzification. The optimization process typically employs a learning algorithm based on gradient descent methods. During training, the parameters of the membership functions and the network weights are adjusted to minimize the error between the actual and predicted values. The process in each layer is described below.

Layer 1—Fuzzification

In this layer, we compute the membership function ${\mu }_i$ from the point $x_t$ to the fuzzy set $S_i$:

$${\mu }_i={\mu }_{S_i}\left(x_t\right)=\exp \left[-{\displaystyle \frac{D_i^2}{d_i^2}}\right]=\exp \left[-{\displaystyle \frac{{({\epsilon }_{t-1}^2-m_{i1})}^2+{({\epsilon }_{t-2}^2-m_{i2})}^2}{d_i^2}}\right],i=1,\dots ,L.$$

Here, $d_i$ is the standard deviation of fuzzy set $S_i$, and $D_i^2$ is the squared distance from the point $x_t$ to the center $m_i=[m_{i1},m_{i2}]$. In other words,

$$D_i^2=\left|\right|x_t-m_i{\left|\right|}^2,$$

where $\left|\right|·\left|\right|$ represents the norm. The membership function ${\mu }_i$ is equal to one when $x_t=m_i$ and decreases further as the distance from the point $x_t$ to the center $m_i$ increases. This reduction occurs when the fuzzy set $S_i$ has very dispersed points or the current point $x_t$ belongs to another fuzzy set.

Layer 2—Rule Evaluation

In this layer, each fuzzy rule’s firing strength is determined based on the degree to which the point $x_t$ matches the premises (the IF portion) of the rules. The firing strength $w_i$, which represents the degree of contribution of each rule to the final output, is given as follows:

$$w_i={\displaystyle \frac{{\mu }_i}{{\sum }_{j=1}^L{\mu }_j}},i=1,\dots ,L.$$

Layer 3—Aggregation

The firing strengths from all rules are combined to produce a single aggregated output $w_i{f}_i$, for $i=1,\dots ,L$.

Layer 4—Defuzzification

We compute the final output by combining the outputs of all the rules into a single crisp value:

$$f_t=\sum _{j=1}^L{w}_j{f}_j.$$

4. Application to a Real Time Series

In this section, we apply the ARIMA-ANFIS model to the daily price series of the BTCUSD currency pair (Bitcoin to US Dollar) using closing prices from 3 August 2023 to 31 July 2024. The time series consists of 364 daily data points. The ARIMA model is used to study the currency price series, and the residuals from this model, which exhibit non-constant variance, are modeled using the ANFIS and GARCH(1,1) model. Additionally, all results from this application were obtained using RStudio program.

4.1. Identification Process for the Conditional Mean Model

We conducted an identification and validation process for the ARIMA model using a training sample consisting of $n=291$ observations. For model selection, we generated a variety of ARIMA(p,d,q) models, exploring combinations for $p,q=0,1,2,3$ and $d=0,1,2$. We assessed their performance using multiple criteria, including the Akaike Information Criterion (AIC), log-likelihood, and the Mean Absolute Percentage Error (MAPE). Additionally, we conducted residual diagnostics to ensure the adequacy of the selected model. Through this evaluation process, we determined that the ARIMA(2,1,2) model provided reliability for our specific dataset. Using the selected model and the forecast package, we obtained fitted values for the training sample and predictions for the testing sample. Figure 3 displays the actual time series (black line), the fitted values (blue line), the predictions (green line), and a red dashed vertical line marking the separation between the training and testing samples.

To evaluate the performance of the ARIMA(2,1,2) model, we use the Root Mean Squared Error (RMSE), the Mean Absolute Error (MAE), and the Mean Absolute Percentage Error (MAPE) for both the training sample and the testing sample. The measures are calculated as follows:

$${\mathrm{RMSE}}_{\mathrm{training}\mathrm{sample}}=\sqrt{{\displaystyle \frac{1}{n}}\sum _{t=1}^n{r}_t^2}$$

$${\mathrm{RMSE}}_{\mathrm{testing}\mathrm{sample}}=\sqrt{{\displaystyle \frac{1}{N-n}}\sum _{t=n+1}^N{e}_t^2}$$

$${\mathrm{MAE}}_{\mathrm{training}\mathrm{sample}}={\displaystyle \frac{1}{n}}\sum _{t=1}^n\left|r_t\right|$$

$${\mathrm{MAE}}_{\mathrm{testing}\mathrm{sample}}={\displaystyle \frac{1}{N-n}}\sum _{t=n+1}^N\left|e_t\right|$$

$${\mathrm{MAPE}}_{\mathrm{training}\mathrm{sample}}={\displaystyle \frac{1}{n}}\sum _{t=1}^n\left|{\displaystyle \frac{r_t}{y_t}}\right|\times 100$$

$${\mathrm{MAPE}}_{\mathrm{testing}\mathrm{sample}}={\displaystyle \frac{1}{N-n}}\sum _{t=n+1}^N\left|{\displaystyle \frac{e_t}{y_t}}\right|\times 100,$$

where n is the training sample size, N is the length of the time series of the BTCUSD price, $y_t$ is the actual BTCUSD price, $r_t$ is the residual from ARIMA(2,1,2), and $e_t$ is the forecast error at time t.

Table 1. Accuracy measures for training and testing samples of ARIMA(2,1,2).

Sample	RMSE	MAE	MAPE	${\mathit{R}}^2$
Training	1343.355	869.390	1.7998	0.9922
Testing	1408.954	1061.159	1.6486	0.8856

presents the performance metrics of the ARIMA(2,1,2) model, evaluated based on RMSE, MAE, MAPE, and $R^2$. These metrics provide insights into the accuracy and reliability of the model. The RMSE for the testing sample is 1408.954 versus 1343.355 for the training sample, and the MAE is 1061.159 for the testing sample versus 869.390 for the training sample. The  RMSE and MAE values are slightly lower for the training sample, but it is expected that the models will perform better on the data they were trained on. On the other side, the MAPE is lower on the testing set compared to the training set, indicating that the model has a good generalization capability. The MAPE values are less than 2%, which means that the predictions are, on average, 2% off from the actual values. This would be considered highly accurate. Additionally, although the $R^2$ value for the training sample gives extremes results, it still performs reasonably well on the testing data, with an $R^2$ value equal to 0.8856. The difference between the training and testing values could indicate some overfitting, but it is not excessive, given that the testing $R^2$ value is still high. The observed values in the accuracy measures suggest that the model is not overfitting and seems to generalize well from the training data to the testing data.

4.2. ANFIS and GARCH Estimation

Consider the time series $\left\{r_t\right\}$ formed by the residuals from the ARIMA(2,1,2) model of length $n=291$ and the forecast errors produce by the ARIMA forecasts of length $N-n=73$.

For the ANFIS estimation, we define four bivariate fuzzy sets, each represented by a corresponding fuzzy rule. The parameters to be estimated include the centers and standard deviations of the four fuzzy sets, as well as the parameters of the consequents corresponding to planes in 3D space. The membership function employed here is the Gaussian function, given as follows:

$${\mu }_i=\exp \left[-{\displaystyle \frac{\left|\right|x_t-m_i{\left|\right|}^2}{d_i^2}}\right],$$

where $d_i$ is the standard deviation of the i-th fuzzy set, and $x_t$ is the point of residuals at time t, defined as $x_t=[{\epsilon }_{t-1}^2,{\epsilon }_{t-2}^2]=[r_{t-1}^2,r_{t-2}^2]$ and $m_i=[m_{i1},m_{i2}]$.

We initialize the parameters of the ANFIS model with random uniform numbers. To optimize the objective function with the mean squared error, we employ the Nelder–Mead method using the optim function.

Figure 4 displays the scatterplot for the ordered pairs $x_t$. Note that there are large squared residuals at the first/second lag and small squared residuals at the second/first lag, as well as small values at both lags. The ANFIS algorithm divides the plane of the squared residuals in four fuzzy sets and performs the entire inference process based on the observations $\left\{r_t^2\right\}$ for $t=1,\cdots ,n$, which is called the training sample. It generates an estimate for the variance, ${\widehat{\sigma }}_{n+1}^2$, which serves as a one-step-ahead forecast. This process is conducted iteratively, with each new observation $r_t^2$, for $t=n+1,\cdots ,N-1$, being added sequentially. At each step, the ANFIS algorithm is re-estimated and applied to produce a one-step-ahead variance forecast for the subsequent time $t+1$ following an expanding window testing strategy.

For the GARCH(1,1) estimation, we use the library rugarch from the R program. The specification was ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), distribution.model = "std").

The GARH(1,1) model is first estimated using the training sample $\left\{r_t^2\right\}$$t=1,\cdots ,n$. Subsequently, a new observation is sequentially added to the time series, similar to the procedure describe above for the ANFIS process. We use the ugarchforecast function to generate one step ahead variance forecasts at each time $t+1$.

We use the standardized squared residuals from the ARIMA model as a proxy for the actual variance at each point. In Figure 5, these residuals are presented alongside the forecasts of ANFIS (blue line) and the GARCH(1,1) model (purple line); the red dashed vertical line marks the separation between the training and testing samples. Note that even when ANFIS is estimated by minimizing the MSE, the volatility predictions are lower than those from the GARCH model. Here, we observe that ANFIS is better at capturing periods of market stability, while the GARCH model excels in capturing high volatility. During calm periods, the GARCH model might overestimate risks, and during volatile periods, ANFIS might underestimate risks.

Table 2. Accuracy measures for the testing sample of ANFIS and GARCH(1,1).

Model	MSE	MAE
ANFIS	4.5539	1.0780
GARCH(1,1)	4.2307	1.3172

presents the MSE and MAE of the testing sample for ANFIS and GARCH(1,1). Even though the MSE for ANFIS is around 10% greater than the MSE for the GARCH(1,1) model, the MAE for ANFIS is around 16% smaller than the MAE for GARCH.

Once we finished the ANFIS and GARCH(1,1) estimation processes, we obtained the 73 volatility forecasts for the residuals of the ARIMA model. These volatility forecasts, ${\widehat{\sigma }}_t$, are used to compute the 95% confidence intervals for the price of the BTCUSD currency around $\widehat{y}$ in the testing sample for both approaches as follows:

$$({\widehat{y}}_t-t_{0.975,N-n}\times {\widehat{\sigma }}_t,{\widehat{y}}_t+t_{0.975,N-n}\times {\widehat{\sigma }}_t).$$

Here, ${\widehat{y}}_i$, $i=n+1,\dots ,N$ are the price forecasts from the ARIMA model in the testing sample, and $t_{0.975,N-n}$ are the percentage points of the t distribution.

Figure 6 and Figure 7 show the confidence intervals for the price predictions using the ANFIS and GARCH(1,1) models, respectively. The varying widths of the confidence intervals are due to the standard deviation within each, which corresponds to the square root of the predicted variance and changes for each time t. It is also noteworthy that some prediction points generated by the ARIMA model have unusually narrow confidence intervals. However, a closer examination of the real data time series reveals that these points coincide with abrupt changes, as illustrated in Figure 6. This could be seen as a positive aspect for currency trading, as the ARIMA model would indicate high-risk operations at these points.

The results revealed that 61.64% of the time, the BTC/USD price fell within the predicted confidence intervals generated by the ANFIS model. In contrast, the GARCH(1,1) model captured the BTC/USD price within its predicted confidence intervals 94.52% of the time. This discrepancy can be attributed to the superior ability of the GARCH model to capture periods of high volatility, leading to larger estimates of volatility compared to ANFIS.

We compared the predictive accuracy of the two competing volatility forecasts using the Diebold–Mariano test. Specifically, we evaluated the test based on two different loss functions: squared error loss and absolute error loss. The loss differential of the two competing models is calculated as follows:

$$d_t=L\left(e_{1,t}\right)-L\left(e_{2,t}\right)$$

where $e_{1,t}$ and $e_{2,t}$ are the prediction errors from the ANFIS and GARCH models, respectively, and $L(·)$ represent the loss function.

In a two-sided Diebold and Mariano test, the high p-value for the squared error loss function, as shown in

Table 3. Predictive accuracy of ANFIS vs. GARCH model based on the Diebold and Mariano test.

Loss Function	p-Value
Absolute error loss	0.031 *
Squared error loss	0.590

* means a significant difference between the models.

, indicates no significant difference between the models in capturing large errors, suggesting that neither model is more accurate during high-volatility periods. However, the p-value for the absolute error loss function, suggests a significant difference when evaluating smaller errors. The mean of the loss differential indicates that ANFIS may better capture low-volatility conditions. Thus, the Diebold–Mariano test suggests that the models perform differently across volatilty regimes: ANFIS is more accurate in low-volatility environments, while neither model shows superiority in high-volatility conditions. This finding aligns with the previous analysis based on the accuracy measures presented in Table 2 and the computed confidence intervals.

5. Conclusions

In this work, we addressed a complex and challenging problem in modeling financial time series. Through a novel and interesting approach, we provided insights and solutions that contribute significantly to the understanding of this issue. Our methodology, which incorporates the processes of identification, specification, estimation, and validation for the ARIMA model and adaptive neuro-fuzzy inference systems, demonstrates the potential for effective application in real-world scenarios.

We implemented the Adaptive Neuro-Fuzzy Inference System (ANFIS), a prominent technique within the domain of soft computing. ANFIS integrates the Takagi–Sugeno–Kang (TSK) fuzzy inference model with neural network methodologies, leveraging fuzzy set theory, IF–THEN fuzzy rules, and fuzzy reasoning. This hybrid system employs diagrams and a connectionist representation, inspired by the functioning of the brain, to effectively model complex and nonlinear systems.

In this article, the ARIMA model is employed to capture the conditional mean, while the ANFIS methodology is used to model the conditional variance of financial series, specifically the daily BTCUSD price. It is noteworthy that conditional variance in these time series is typically modeled using GARCH models; hence, applying ANFIS methodology in this context is innovative. By combining the econometric approach (ARIMA) with the soft computing technique (ANFIS), we jointly model both the conditional mean and conditional variance, creating a hybrid ARIMA-ANFIS model.

The comparison between the benchmark ARIMA-GARCH model and the proposed ARIMA-ANFIS model reveals that each model captures different aspects of data dynamics. While the ANFIS model is effective in certain scenarios, it may underestimate volatility during turbulent periods, potentially exposing users to unexpected risks, as illustrated between observations 335 and 340 in Figure 6. Conversely, the GARCH(1,1) model, by generating higher volatility estimates, might lead to excessive caution, potentially reducing returns. This highlights the trade-offs between the two models: ANFIS offers a more conservative approach in stable markets, whereas GARCH(1,1) provides a robust defense against high volatility, but at the cost of possibly missing out on opportunities during calmer periods, as shown in Figure 7.

Parameter optimization for the ANFIS model proved to be a time-consuming procedure, highlighting the need for more efficient optimization techniques. As future work, the implementation of evolutionary algorithms could be explored. Evolutionary algorithms, with their robustness and global search capabilities, have the potential to significantly streamline the optimization process by efficiently navigating the complex search space of the ANFIS parameters.

The proposed ARIMA-ANFIS model adequately captured some of the dynamics in the treated financial time series, providing good forecasts and confidence intervals in most cases. Additionally, testing this model using other time series exhibiting non-constant conditional variance could be considered.