Effective short-term forecasts of Saudi stock price trends using technical indicators and large-scale multivariate time series

Time series forecasting

Time series is a collection of sequential vector-valued variables of some phenomenon indexed over time. Time series forecasting is the task of predicting future values of given time-series observations by analyzing historical time series data to check for patterns of time decomposition, such as trends, seasonal patterns, cyclic patterns, and regularity. Univariate time series implies that forecasting is based on a single variable of time series observations without considering the effect of the other variables, making them desirable for computational efficiency especially when data availability is already scarce. On the other hand, multivariate time series have more than one time-dependent variable where each variable depends not only on its past values but also tends to have some interdependency with other variables. Although the multivariate time series setting is known for difficulty in extracting beneficial information for the prediction process (Chayama & Hirata, 2016), informative forecasts are often assumed to be conveyed by multivariate time series since advanced statistical methods can leverage the dependencies between variables to improve the forecasting accuracy (Yin & Shang, 2016; Zhang et al., 2017; Lang et al., 2019).

Problem definition

Multivariate time series forecasting has numerous valuable applications, such as predicting electricity consumption, solar power production, and stock trend movements. The potential interdependency among variables in the multivariate-based time series can lead to better and more coherent learning results, especially when variables are highly correlated. This research seeks to employ the multivariate-based time series analysis to enable better information extraction about the data and their underlying dynamics to improve the short-term prediction accuracy.

VARMA with exogenous regressors

The Vector Autoregressive Moving-Average (VARMA) model extends the ARMA/ARIMA model to work with time series with multiple response variables (vector time series). It is one of several statistical analyses frequently used for multivariate-based time series studies. The model includes the dynamic relationship between the multiple response variables as well as the relationship between the dependent and independent variables. It is a combination of Vector Auto-Regressive (VAR) and Vector Moving-Average (VMA) models that supports the multivariate time series modeling by considering both the lag order (p) and order of moving average (q) in the model. The VARMA model can operate as a VAR model by just setting the q parameter to 0, and it also can operate like a VMA model by just setting the p parameter to 0. Moreover, the Vector Autoregression Moving-Average with Exogenous Regressors (VARMAX) is an extension of the VARMA model that includes modeling with exogenous variables (also called covariates). In practice, the target series is referred to as an endogenous sequence to contrast from exogenous sequences that are included directly in the model at each time step but not modeled in the same way as the target endogenous sequence. Like VARMA, VARMAX can also operate as VARX and VMAX to model the subsumed models with exogenous variables.

Formally, the endogenous vector ${\mathbf{y}}_t$ is said to follow a vector autoregressive moving average with exogenous variables $VARMAX(p,q,s)$ if it satisfies an equation of the form (Östermark, 1994; Kadiyala & Kumar, 2014; SAS, 2016):

(1) $${\mathbf{y}}_t=\sum _{i=1}^p{\mathrm{\Phi }}_i{\mathbf{y}}_{t-1}+\sum _{i=0}^s{\mathrm{\Theta }}_i^{\ast }{\mathbf{x}}_{t-1}+{ϵ}_t-\sum _{i=1}^q{\mathrm{\Theta }}_i{\epsilon }_{t-1}$$where the output variables of interest, ${\mathbf{y}}_t=(y_{1t},\cdots ,y_{kt})$ referred to response or endogenous variables, can be influenced by other input variables, ${\mathbf{x}}_t=(x_{1t},\cdots ,x_{mt})$ referred to as regressor or exogenous variables. The unobserved noise variables, ${\epsilon }_t=({\epsilon }_{1t},\cdots ,{\epsilon }_{kt})$, are a vector white noise process. Therefore, the model $VARMAX(p,q,s)$ can be written as:

(2) $$\mathrm{\Phi }(B){\mathbf{y}}_t={\mathrm{\Theta }}^{\ast }(B){\mathbf{x}}_t+\mathrm{\Theta }(B){\epsilon }_t$$where,

$$\begin{array}{rl}& \mathrm{\Phi }(B)=I_k-{\mathrm{\Phi }}_{1}B-\cdots -{\mathrm{\Phi }}_p{B}^p\\ & {\mathrm{\Theta }}^{\ast }(B)={\mathrm{\Theta }}_0^{\ast }+{\mathrm{\Theta }}_1^{\ast }B+\cdots +{\mathrm{\Theta }}_s^{\ast }{B}^s\\ & \mathrm{\Theta }(B)=I_k-{\mathrm{\Theta }}_{1}B-\cdots -{\mathrm{\Theta }}_q{B}^q\end{array}$$are matrix polynomials in B in the backshift operator, such that $B^i{\mathbf{y}}_t={\mathbf{y}}_{t-1}$, the ${\mathrm{\Phi }}_i$ and ${\mathrm{\Theta }}_i$ are $k\times k$ matrices, and ${\mathrm{\Theta }}_i^{\ast }$ are $k\times m$ matrices. For the sake of achieving the complete definition, the following assumptions are made:

• $E({\epsilon }_t)=0$, $E({\epsilon }_t{\epsilon }_t^{\mathrm{\prime }})=\sum$ which is positive-definite, and $E({\epsilon }_{\mathbf{t}}{\epsilon }_{\mathbf{s}}^{\mathrm{\prime }})=0$ for $t\ne s$.

• For stationarity VARMAX, the roots of $|\mathrm{\Phi }(z)|=0$ and $|\mathrm{\Theta }(z)|=0$ are outside the unit circle.

• The exogenous independent variables ${\mathbf{x}}_t$ are not correlated with residuals ${\epsilon }_t$, $E({\mathbf{x}}_t{\epsilon }_t^{\mathrm{\prime }})=0$.

• The exogenous variables can be stochastic or non-stochastic.

Gated recurrent units

We first present a brief exposition of recurrent neural networks (RNNs). RNN, first conceived by Elman (1990), is an extension of the conventional feedforward neural network used for sequence modeling problems. RNNs can be seen as a universal approximator for dynamical systems of variable-length input sequences (Schäfer & Zimmermann, 2006). RNNs can let information persist within the network over time, making them suitable for sequence learning modeling. Formally, given a sequence $\mathbf{X}=(x_1,x_2,...,x_t)$, a default RNN updates a recurrent hidden state ${\mathbf{h}}_t$ at time $t$ as (Olah, 2015):

(3) $${\mathbf{h}}_t=f({\mathbf{h}}_{t-1},x_t)$$

where $f$ is commonly chosen as a nonlinear activation function, such as sigmoid or tanh, with an affine transformation of both $h_{t-1}$ and $x_t$. Thus, the update of the recurrent hidden state can be rewritten as:

(4) $${\mathbf{h}}_t=\mathrm{\Phi }(W_{x_t}+\mathrm{\Theta }{\mathbf{h}}_{t-1}+b)$$where $x_t$ is $m$-dimensional input vector at time $t$, ${\mathbf{h}}_t$ is $n$-dimensional hidden state, $\mathrm{\Phi }$ is a nonlinear activation function, W is weight matrix of size $n\times m$, $\mathrm{\Theta }$ is recurrent weight matrix of size $n\times n$, and $b$ is a bias of size $n\times 1$.

The default RNN is deficient in capturing long-term dependencies because the gradients tend to either vanish or explode with long sequences (Bengio, Simard & Frasconi, 1994). Recall a well-known phenomenon in canonical recurrent neural networks (RNNs) called long-term dependencies (Bengio, Simard & Frasconi, 1994), in which a system’s desired output at time T depends on inputs presented at times $t<T$. In order words, as the sequence length increases during training, the gradients tend to vanish, leading to the problem of vanishing gradients. Vanishing gradients prevent the network from learning the underlying dynamics and hidden patterns extracted from the long-term dependencies, which is beneficial in producing accurate future forecasts. Vanishing gradients also arise when the network dynamics have attractors necessary to store bits of information over the long term (Chandar et al., 2019). Successful architectures, such as LSTM and GRU, alleviate this issue by introducing memory in the form of information skipped forward across transition operators, with gates bounded by activation functions (sigmoid, tanh) to determine which information to skip, releasing the hidden states (i.e., the memory) from their output states to reduce the vanishing problem.

To overcome the shortcoming mentioned above, two RNN-based variants have been introduced, namely Long Short Term Memory (LSTM) (Hochreiter & Schmidhuber, 1997) and Gated Recurrent Units (GRU) (Chung et al., 2014). Both variants seek to keep long-term dependencies effective while alleviating the vanishing and exploding gradient problems (Chandar et al., 2019; Aseeri, 2021). As seen in Fig. 1, GRU operates on two gates: reset and update gates. The architecture essentially combines the forget and input gates used by Long term short memory (LSTM) into the update gate and merges the cell state and hidden state into the reset gate, resulting in a simpler design and consequently lower number of trainable parameters needed compared to the LSTM variant, improving the training time and performance (Chung et al., 2014). The operations performed by the GRU model are as follows (Chung et al., 2014; Dey & Salem, 2017):

(5) $$\begin{array}{rl}& r_t=\sigma (W_r\cdot [{\mathbf{h}}_{t-1},x_t])\\ & {\theta }_t=\sigma (W_u\cdot [{\mathbf{h}}_{t-1},x_t])\\ & {{\mathbf{h}}^{\mathrm{\prime }}}_t=tanh(W\cdot [r_t\ast {\mathbf{h}}_{t-1},x_t])\\ & {\mathbf{h}}_t=(1-{\theta }_t)\ast {\mathbf{h}}_{t-1}+{\theta }_t\ast {{\mathbf{h}}^{\mathrm{\prime }}}_t\end{array}$$where $r_t$ and ${\theta }_t$ are respectively the reset and update gates at time $t$, $W_{\ast }$ are weight matrices, and $\cdot$ denotes component-wise multiplication.

Data description and exploration

To achieve the best possible forecasting results, especially in the time-series domain, it is important to employ long historical and well-preprocessed time series data to help capturing the hidden patterns for the domain of interest, such as stocks’ trends and movement. In this experiment, we use a recently published multivariate time series of the primary Saudi stock market index, particularly the Tadawul All Share Index (TASI), made freely accessible by TASI (2020). The TASI dataset represents daily closing prices of TASI index companies, which are distributed into ten sectors featuring 14 attributes as shown in Table 1 where eight of which are numerical features: open, high, low, change, percentage change, volume traded, value traded, and number of trades. Since the Saudi stock market index is as old as 2001, some companies in the dataset have been recorded since 2001, while others have recently joined the market, meaning that the recording data of a new company begins immediately starting from the joining date. We randomly selected ten companies, one from each sector, to be modeled by our proposed method and listed in Table 2. In addition, Fig. 2 illustrates the distribution of the closing price of each of the ten chosen companies.

Temporal multivariate processing

In this study, the process of Saudi stock trend prediction is considered a regression problem where the goal is to produce future horizon closing price forecasts of size $h$, where $h$ is set to five representing the five trading days, given a past window of multivariate time-lagged historical time series¹ . Moreover, this work explores the performance of the proposed prediction system when employing a number of financial-related technical indicators to improve forecasts. To the best of our knowledge, this is the only work reporting multi-step closing price forecasts of the Saudi stock TASI index. The proposed methodology is visualized in Fig. 3, which depicts the processing pipeline (top subfigure) and the windowing approach applied to the raw time series data (bottom subfigure).

The forecasting pipeline shown in the top subfigure of Fig. 3 is tailored to ensure systematic and sound preprocessing and modeling. The pipeline involves five essential modules, each performing a specific task. The process begins with initiating a streaming data reader from the raw stock data to enable retrieving multivariate time series data. When a specific company’s multivariate time series is fully retrieved, the data interpolation module begins cleaning the retrieved multivariate time series by (1) removing duplicate timesteps when detected and (2) performing linear interpolation using the forward technique by replacing missing values with the last valid observation propagated to fill the consecutive gaps (the missing values are the weekend days as the market is closed). Next, the feature engineering module is responsible for conducting three major processes: (i) feature reconstruction, (ii) correlation-based feature selection, which will be elaborated in detail in the next subsection “Feature Engineering: Technical Indicators and Feature Selection”, and (iii) the min-max normalization. After that, the windowing and splitting module performs a customized one-stride sliding window procedure to form adequate sequences for the modeling phase, followed by splitting the overall sequences into three sets, namely training, validation, and testing sets. Here, the word stride implies the amount of movement/shift to the right direction in a time, so one-stride means one movement to the right at a time. Looking closely at the left side of the windowing approach subfigure in Fig. 3, an iterative top-down procedure applied to the multivariate time series data generates multivariate sequences with a vector of future horizons, thus constructing supervised-like multistep sequences ready for training and validating the forecasting model. For instance, at a timestep $t$, a one-stride sliding window is rolled to specify a batch that involves the historical time-lag window of size $w$ (in blue) and the forecast window of size $h$ (in green). Consequently, his process can be formally expressed as:

$$y_{t+h}=f_h\{y_t,\dots ,y_{t-w}\}$$where $h$ is the number of steps to forecast into the future (i.e., the desirable future horizons), $w$ is the size of the time-lag window, $t$ is the timestep (i.e., the day), and $f_h$ is any arbitrary learner. Furthermore, it is well-recognized that the most critical parameter in multivariate time series forecasting is the proper determination of the number of past observations, i.e., the time-lag window. Thus, we performed a trial and error method with different sizes of multivariate time-lag windows, including 1, 5, 10, and 15 days, to choose the best possible window size that yields better performance. We concluded that the time-lag window of size 5 days (i.e., $w=5$) works well in providing satisfactory forecasting results while maintaining an efficient processing time. Finally, both model construction and model evaluation modules will be elaborated on in the upcoming subsection “Modeling Design”.

Feature engineering: technical indicators and feature selection

Technical indicators (TIs) are heuristic pattern-based indicators produced by a feature of interest, such as the price, volume, etc., which are usually utilized as exogenous features to aid in identifying possible hidden patterns, trends, and movements within the historical time series. In practice, the success of employing the technical indicators relies heavily on the availability of sufficient historical data, making them conceivably applicable in this experiment. In this study, we plan to explore the effectiveness of employing technical indicators as additional exogenous features to the input multivariate TASI index data. We hypothesized that we would get enhanced forecasting results when exploiting technical indicators as they can improve the interdependency between the exogenous features with respect to the endogenous target variable, hence better capturing hidden patterns and trends within the historical time series.

Accordingly, we perform feature reconstruction within the feature engineering module where the original features, listed in Table 1, are extended using an arbitrary set of financially-related technical indicators, computed with respect to the target feature of interest in this study, i.e., the close price. We choose five essential and widely adopted technical indicators in finance, namely Weighted Moving Average (WMA), Double Exponential Moving Average (DEMA), Triple Exponential Moving Average (TEMA), Hull Moving Average (HMA), and Zero-Lag Exponential Moving Average (ZLEMA), thus resulting fourteen features in total (nine original numerical features, five technical indicators features). The calculation is carried out in two steps: (1) calculating a new vector of attributes on each trading day, representing the values of the selected technical indicators where the input parameter lag in each indicator is set to 5, and (2) extending the original exogenous features with the new computed vector. For the sake of completeness, we show the formulations of the indicators mentioned above as follows: if $n$ is the number of periods, $w$ is the weighting factor, $\alpha$ is a smoothing factor, and $x_t$ is the data observation at the period $t$, the well-known Weighted Moving Average (WMA) is defined as (Perry, 2010):

(6) $$WMA={\displaystyle \frac{{\sum }_{t=1}^n{w}_t{x}_t}{{\sum }_{t=1}^n{w}_t}}$$also, the Double Exponential Moving Average (DEMA) and Triple Exponential Moving Average (TEMA) are defined as (Rozario et al., 2020):

(7) $$\begin{array}{cc}DEMA\hfill & =2\times EMA-EMA(EMA)\hfill \end{array}$$

(8) $$\begin{array}{cc}TEMA\hfill & =(3\times EMA-3\times EMA(EMA))+EMA(EMA(EMA))\hfill \end{array}$$where,

$$\begin{array}{cc}EM{A}_t& =\{\begin{array}{cc}{x}_0& t=0\\ \alpha x_t+(1-\alpha )EM{A}_{t-1}& t>0\end{array}\end{array}$$the Hull Moving Average (HMA) makes a moving average more responsive while maintaining a curve smoothness. It uses a combination of three WMAs with different sizes of period, and it is computed as Hull (2005):

(9) $$HMA=WM{A}_{\sqrt{n}}(2\times WM{A}_{{\displaystyle \frac{n}{2}}}(X)-WM{A}_n(X))$$finally, the Zero-Lag Exponential Moving Average (ZLEMA) modifies the Exponential Moving Average (EMA) to greatly reduce lag. It is computed as (Ehlers & Way, 2010):

(10) $$ZLEMA=K(2\times x_0-x_{-lag})+(1-K)ZLEM{A}_{-1}$$where,

$$\begin{array}{cc}K\hfill & ={\displaystyle \frac{2}{(n+1)}}\hfill \\ lag\hfill & ={\displaystyle \frac{(n-1)}{2}}\hfill \end{array}$$

Next, to ensure that the newly added technical indicators and the original exogenous features have the desired level of relevance and interdependencies, we adopt Pearson’s correlation as a correlation-based feature selection method to return the most highly-correlated features with respect to the target variable. Pearson’s correlation can be utilized to translate the linear relationship’s strength between every two variables. It outputs values in the interval [−1, +1] where +1 implies total positive correlation, −1 total negative correlation, and 0 no correlation. As per the standard literature, for a pair of variables $(X,Y)$, the linear correlation coefficient $r$ is given by:

$$r={\displaystyle \frac{\sum (X_i-{\overline{X}}_i)(Y_i-{\overline{Y}}_i)}{\sqrt{\sum {(X_i-{\overline{X}}_i)}^2}\sqrt{\sum {(Y_i-{\overline{Y}}_i)}^2}}}$$

This study found that all correlation matrices rendered from the selected ten companies have an identical strong positive correlation between only nine out of the total fourteen features, namely: open, high, low, close, WMA, DEMA, TEMA, HMA, ZLEMA, which can be observed in Fig. 4 where we visualize the correlation matrix (half diagonally since it is symmetric) from a selected company, i.e., the Saudi Arabia Refineries. The rest of the correlation matrices of the remaining eight companies involved in this work are shown in Fig. A1 in the Appendix.

Modeling design

After completing the data preprocessing phase, we begin modeling the processed multivariate time series data using a carefully-engineered Gated Recurrent Units (GRU) model. This work deliberately applies the gated recurrent units method for multivariate time series forecasting as one of the most effective RNN-based methods that provide accurate results with efficient computation. We spent quite some time experimenting with varying design architectures and hyperparameter tuning options to yield the best possible forecasting results, including Fig. 5 that demonstrates the most suitable model design and configuration given the input preprocessed multivariate time series data. This final architecture involves one GRU layer composed of 256 GRU processing units (see Fig. 1 for depicting the generic GRU processing unit) and followed by a dense layer of linear activation units equal to the size of the forecast horizon $h$ (in this work, $h=5$).

In order to verify the correctness of the GRU forecasting estimates, we need to compare them with a baseline statistical method such as AutoRegressive Integrated Moving Average (ARIMA). Since we operate on the multivariate time series modeling, we employ an ARIMA variant designed for handling exogenous data, called Vector Autoregressive Moving-Average with Exogenous Regressors (VARMAX) (see subsection “VARMA with Exogenous Regressors” for more details about this variant). The comparison results will be illustrated later in the upcoming section “Results and Discussion”.

Experimental setup

The experimental environment is implemented using the Python programming language (Python 3). The neural network-based modeling framework of choice is TensorFlow v2.4.1. We also implemented the VARMAX modeling with the support of pmdarima API v2.0.0 (Smith, 2017). To ensure sound model training and evaluation, the input time series is splitted into three sets: train set starts from the first-ever trading day until 31-12-2015, the validation set is from 01-01-2016 to 31-12-2016 while the testing set is from 01-01-2017 to 01-04-2020. The experiment was performed on a workstation having Intel Core i7-10700KF (3.8 GHz) CPU, Nvidia Geforce RTX 3060 Ti (8 GB) GPU, and 64 GB of system RAM. For each company involved in this experiment and listed in Table 2, the experiment was run five times, and we reported only the best run.

The neural network-based methods seek to minimize errors between predicted and true values assessed using a suitable cost function by iteratively training and tuning the model’s hyperparameters using a gradient-based descent optimization method. For training the GRU model, we found out that Nesterov-accelerated Adaptive Moment Estimation (NADAM) (Dozat, 2016), with the learning rate of $1\times {10}^{-4}$, works well during the training and the validation of the GRU model design. The mean squared error was appointed as the cost function of choice to estimate the model’s errors, while the batch size fed to the model during the feedforward process was fixed to be 10% of the training data. Furthermore, the activation functions for the GRU gates are left to the conventional default. That is, for both the update and reset gates, the activation is the sigmoid function, while the hyperbolic tangent function (tanh) is used with the recurrent activation gate. For the sake of reproducibility, Table 3 reports the hyperparameters options along with their description of each model used in this experiment.

Evaluation metrics

The performance of neural network-based methods is generally evaluated in terms of metrics related to the models’ discriminative power. In principle, making a proper and correct selection of evaluation metrics for quantifying prediction performance enable the correct interpretation of the model’s skill and capability in making the forecasts. Therefore, we choose a regression-based metric set that assesses our models’ precision for the multi-step forecasting task. Various evaluation metrics are available for calculating the accurateness of the forecasts, but the selected metrics shown in Table 4 have been the default metrics for regression problems as well as widely adopted for evaluating the performance of time series forecasts. Here, $y$, $x$, and $n$ represent the forecasted value, the actual value, and the number of sample data for the period, respectively.

The evaluation metrics shown in Table 4 are commonly used in predictive modeling. That is, mean absolute error (MAE) measures the average error magnitude in a set of predictions using the absolute value. This type of error measurement is proper when measuring prediction errors in the same unit as the original series. The mean squared error (MSE) formula replaces the absolute value with a square, giving large errors a relatively greater influence on MSE than smaller ones. MSE has nice mathematical properties, making it easier to compute the gradients. Similarly, root means squared error (RMSE) measures the average magnitude of the error but with a square root on the average squared differences between prediction and actual observation (MSE and RMSE are monotonically related, i.e., through the square root). Finally, the coefficient of determination (R²) is the goodness-of-fit measurement of how well the regression line approximates the actual observation at the 0–1 scale (or 0–100%) where a higher value is considered desirable as opposed to the MAE, MSE, and RMSE where lower is deemed desirable. The sum squared regression (SSR) is the sum of the residuals squared, and the total sum of squares (SST) is the sum of the distance the data is away from the mean all squared.