Forecasting hourly WTI oil front monthly price volatility densities

1.   Introduction

Volatility is a crucial concept in financial markets, and its importance is widely recognized in both academic research and practical finance. Key modern financial theory and practice reasons are risk measurement and management, option pricing, portfolio construction (diversification), and a general understanding of market dynamics (microstructure). Multifactor volatility models are simple to simulate and estimate, and they appear to perform better in derivative calculations due to their smooth sample paths, compatibility with standard hedging arguments, and compatibility with Itô-calculus. Multifactor models, in contrast to jump diffusion models, are more appealing for implementing machine learning (ML) and artificial intelligence (AI) for volatility predictions. The applications of multifactor models are numerous, for example, options pricing, risk assessment, algorithmic trading, portfolio optimization, and general market surveillance. So, this study creates and tests multifactor scientific stochastic volatility (SV) models from high-frequency WTI oil front month (WTI oil) price changes to predict the latent volatility of the liquid fossil fuel market. Volatility is a measure of how widely apart an asset's price swings are from one another. All volatility models must have the ability to anticipate future price changes to be successful. Internationally, volatility models have successfully predicted the absolute size of movements, quantiles, and full densities. Asset volatility is unique in that it is latent or not immediately visible, making the evaluation of volatility models' predicting abilities challenging. I estimate the most liquid hourly WTI oil price movements (continuous) using a general method of moment (GMM) estimation. The SV models give access to conditional moments and potential forecasts (data dependence). The analysis is univariate (non-synchronous trading is irrelevant) and has mainly three objectives. First, an objective to find a general step ahead in WTI oil densities; second, an objective to identify data dependence for potential predictability; and third, an objective to report systematic market features and potential market features applying machine learning techniques.

Gallant and Tauchen (1987, 1992) introduced semi-nonparametric time series analysis (SNP densities)¹ as the methodology. Unlike traditional estimating methods, nonparametric estimation does not depend on selecting data from a known distribution. Nonparametric models, on the other hand, use the underlying data to establish the model structure. The conditional density of the time series processes (non-normality) is approximated by the method using an expansion of hermite functions. For time series data the hermite expansion is attractive both for modelling and computation. For modelling, the Gaussian components of the hermite expansion make it easy to absorb familiar time series models. For computation, the hermite density is easy to evaluate and differentiate. Furthermore, the SNP implementation is easy to sample from (simulations). The SNP expansion process is therefore well approximated using a well-known parametric model as the leading term in the model expansion process; higher-order terms, such as hermite functions, show how the process deviates from the model (Robinson, 1983). To find the proper order of expansion, I fit the SNP model using traditional maximum likelihood and a model selection approach (BIC) (Schwarz, 1978). To compute the nonlinear functions of densities², the model is well-designed.

¹ SNP/EMM: A Program for Efficient Methods of Moments Estimation, Duke University, 09.08.2022 (http://econ.duke.edu/webfiles/arg/emm).

² The computer cluster at NTNU, Faculty of Economics and Management, Trondheim is used for estimation/implementation. A special thanks to Professor Asgeir Thomasgaard at NTNU, for access to the computer cluster.

The remainder of this paper is organized as follows: In Section 2, I discuss the WTI oil front month contracts and look at some literature on the semiparametric (SNP) time-series model, stochastic volatility, and the use of the efficient method of moments (EMM). For the volatility forecasting classical forecasting (OLS) and machine learning (ML) and the lasso, ridge and Decision Forest regression models are presented. Finally, in section 2, neural networks are described briefly for prediction purposes. In Section 3.1, the semi-nonparametric (SNP) establish a maximum likelihood consistent mean and volatility equation details from an optimal BIC criterion (1978). The model approximation for the conditional density, which summarizes the probability distribution and describes the price movement processes, is expanded by the hermite function expansions (non-normality). The residual characteristics report t-statistics for all moment well below two indication a well approximate model fit. In Section 3.2, I look at the effective method of moments (EMM) for estimating SV specifications and use the non-linear Kalman filter to functionally calibrate the volatility vector that hasn't been observed. The methodology evaluates one or two volatility factors for WTI oil. Potential data dependence evaluates the predictability of the volatility factors. In section 4, I study predictions using methodologies from section 2.4. In section 5, I summarize and conclude the study. Appendix A shows similar results for the daily WTI oil data series.

2.   Literature review

2.1.   WTI oil front month prices

The fossil oil literature shows a link between the oil and stock markets, which appears not only in return but also in volatility. Clark (1973), Tauchen and Pitts (1983), and Ross (1989) observe that the rate of information flow to a market correlates with the volatility of an asset, not its return. Thus, volatility is a good measure of information flow among markets. Exploring information flow may generate new insights. For instance, the study by Vo (2011) shows a bidirectional dependence in volatility between stock and oil markets. That is, shocks to either market help predict not only volatility in their own market but also that in the other market.

I attempt to extract information from WTI oil prices and their movements. It will not focus on the relationship between stock and oil returns, which many papers cited above have studied, but aim more at volatility models that can extract useful information and information flow with good forecasting power. Modelling and forecasting volatility are very important for at least two reasons. First, volatility is an important variable for pricing derivatives, whose trading volume has quadrupled in recent years. Furthermore, volatility is an important input in risk management. For instance, I use it to construct optimal hedge ratios to mitigate risk and estimate the value at risk, to name just two applications. Second, in order to make efficient econometric inferences about a variable's mean, I need a correct specification of its volatility.

2.2.   The Semi-Non-Parametric (SNP) model

Our work will suggest conditional models using non-linear stochastic models. To describe the structure of the conditional volatility, I use the term generalized autoregressive conditional heteroscedasticity (GARCH), and to describe the structure of the conditional mean, I use the word autoregressive and moving average (ARMA). In contrast to the ARCH specifications, which were initially explored by Engle (1982) and advanced by Bollerslev (1986), who described the generalized ARCH or GARCH, ARMA models can be investigated in detail in, for instance, Mills (1990). Initially, the number of delays in the ARCH specification was the reason for the development of GARCH from ARCH. The volatility is defined by ARCH/GARCH as a function of historical price changes and volatility. Quite a few studies have demonstrated how the findings from this work have been utilized in the literature on international finance. See, for instance, Nelson (1991), Bollerslev et al. (1987, 1992), Engle et al. (1986, 1993), and deLima (1995a, 1995b). Gouriéroux (1997) provides a thorough introduction to ARCH models and their uses in finance. While Glosten et al. (1993) described the truncated GARCH (GJR), Ding et al. (1993) extended the symmetric GARCH model into an asymmetric GARCH.

Gallant et al. (1992, 2010) use the term semi-nonparametric, or SNP³, to refer to a method that lies between parametric and nonparametric approaches. Higher-order terms in the series expansion are deviations from the existing parametric model, which is the leading term and is known to provide a decent approximation of the process. The method's theoretical underpinning is the hermite series expansion, which is particularly appealing for time series data due to modelling and computational considerations. In terms of modelling, the Hermite expansion's Gaussian component makes it simple to incorporate into popular time series models, including VAR, ARCH, and GARCH models (Engle, 1982; Bollerslev, 1986). These models typically provide excellent initial approximations for a wide range of applications. Hermite density is simple to assess and distinguish computationally. Furthermore, because its moments are higher than those of the normal, which can be calculated using conventional recursions, they are simple to assess. Finally, sampling from a hermite density is doable, which makes simulation easier.

³The code and user guide are available at http://www.aronaldg.org. The programme is published under the terms of the GNU General Public Licence, version 2 or, at your option, any later version, which is published by the Free Software Foundation. You are free to redistribute and/or modify the programme.

2.3.   Stochastic volatility⁴

⁴ See also Solibakke (2020, 2022).

Instead of directly defining the predicted distribution of price returns, the stochastic volatility (SV) method does so indirectly through the model's structure. The suggested one-step-ahead distribution of returns recorded over any arbitrary time interval appropriate for the econometrician is not a concern for the SV model because it has its own stochastic process. I begin with the work of Andersen et al. (1994, 2002), which examines the well-known stochastic volatility diffusion for a given stock price S_t, as indicated by Equation (1).

There are two types of unobserved volatility processes: log linear and square root (affine), for ${W}_{1, t}$, i = 1, 2. Standard Brownian motions ${W}_{1, t}$, and ${W}_{2, t}$ may be correlated by corr(d${W}_{1, t}$, d${W}_{2, t}$) = ρ. The μ is mean drift and c is a parameter of volatility-in-mean. The stochastic volatility model was estimated by Andersen et al. (1994, 2002) using daily S&P 500 stock index data spanning from 1953 to December 31, 1996. They sharply reject both SV model versions. A basic SV model, on the other hand, benefits a lot from having a jump component added because it captures two well-known features: Fat non-Gaussian tails and persistent time-varying volatility. Chernov et al.'s (2003) results for an SV model with two stochastic volatility variables are positive. An affine setup and a logarithmic setup are two major categories of setups that the authors take into consideration for the volatility index functions and factor dynamics. The authors estimate the models using daily data on the Dow Index from January 2, 1953, to July 16, 1999. They discover that models with two volatility variables perform significantly better than models with just one. Additionally, they discover that logarithmic two-volatility component models outperform affine jump diffusion models and provide a good fit to the data. One of the two sources of volatility is quite persistent, whereas the other is substantially mean reverting.

I apply the logarithmic model with two stochastic volatility variables. I expand the model to facilitate the correlation between the mean and the stochastic volatility variables. I use the Cholesky decomposition in the correlation to ensure consistency. The introduction of asymmetry effects (correlation between return innovations and volatility innovations) is the key justification for correlation modelling. With ${W}_{1, t}$ i = 1, 2, and 3 being standard Brownian motions (random variables), the generic SV model formulation for price change processes (y_t) is therefore Equation (2):

The vector of parameters is $\left({a}_{0}, {a}_{1}, {b}_{0}, {b}_{1}, {c}_{0}, {c}_{1}, {s}_{1}, {s}_{2}, {r}_{1}, {r}_{2}, {r}_{3}\right)$. An internally consistent variance/covariance matrix is enforced by the correlation coefficients, or the three correlation parameters, ${r}_{1}, {r}_{2}, {r}_{3}$, obtained by a Cholesky decomposition. ${c}_{0}$ is set to zero not to have two constants in the volatility equation. Taylor (1986), Clark (1973), Tauchen and Pitts (1983), and Rosenberg (1972) are early references. More recent references include Shephard (2004), Andersen (2002), Durham (2003), Gallant et al. (1993, 2010), Taylor (1982, 2005), and Chernov et al. (2003).

The model above contains three stochastic factors (Solibakke, 2020). It is also possible to apply Poisson distributions to jumps, although this greatly complicates calculations. To do a statistical analysis on a stochastic volatility model produced from a scientific procedure, the research employs a computational methodology suggested by Gallant and Tauchen (2010). I first compute the tractable likelihood function of a reduced-form auxiliary model (generous parameterization), which intuitively explains the strategy. The estimated set of score moment functions encode crucial details about the raw data sample's probabilistic makeup. I then use the continuous-time SV model to simulate a lengthy sample. I adjust the parameters in order to maximize the quasi-score moment functions that are assessed using the Metropolis-Hastings algorithm and parallel computing on the simulated data. Among the useful by-products are an explicit metric for assessing the severity of SV model failure and an extensive set of model diagnostics. The scientific stochastic volatility model is easy to simulate, but it cannot generate likelihoods.

The previous SV model estimation yielded a long-simulated realization of the state vector ${\left\{{\widehat{V}}_{i, t}\right\}}_{t = 1}^{N}, i = \mathrm{1, 2}$ and the accompanying ${\left\{{\widehat{y}}_{t}\right\}}_{t = 1}^{N}$ for $\theta = \widehat{\theta }.$ By calibrating the functional form of the conditional distribution of functions supplied ${\left\{{\widehat{y}}_{\tau }\right\}}_{\tau = 1}^{t}$, assessing the outcome on observed data ${\left\{{\stackrel{~}{y}}_{t}\right\}}_{t = 1}^{n}$, and producing predictions for ${V}_{i, t}, i = \mathrm{1, 2}$ through Kalman filtering on y_t, extremely broad functions of can be employed and a large dataset is made available. An SNP model is estimated using the ${\widehat{y}}_{t}$. A conditional variance of provided that is one step forward is represented by the final model. Regressions are performed on these series' ${\widehat{V}}_{i, t}$generously lengthy lags, ${\widehat{\sigma }}_{t}^{2}$, ${\widehat{y}}_{t}$ and ${\left\{{\stackrel{~}{y}}_{\tau }\right\}}_{\tau = 1}^{t}$. Values for the volatility factors ${\stackrel{~}{V}}_{i, t}, i = \mathrm{1, 2}$ at the original data points are obtained by evaluating these functions on the observed data series ${\left\{{\stackrel{~}{y}}_{\tau }\right\}}_{\tau = 1}^{t}$.

2.4.   Step ahead forecasting

2.4.1.   Classic static regression

Classic Static Regression: Static forecasting repeatedly forecasts the dependent variable one step ahead. That is, compute each observation in the forecast sample, compute the${\widehat{y}}_{s+k} = {\widehat{\alpha }}_{1}+{\widehat{\alpha }}_{2}\cdot {x}_{s+k}+{\alpha }_{3}\cdot {z}_{s+k}+{\widehat{\alpha }}_{4}\cdot {y}_{s+k-1}+... +{\widehat{\alpha }}_{N+4}\cdot {y}_{s+k-N}$, always using the actual value of the lagged endogenous variable. I must observe data for both exogenous and any lagged endogenous variables for observations in the forecast sample. This type of volatility forecasting only uses endogenous variables.

2.4.2.   Machine learning regression models

In this section, I will define various machine learning regression techniques. I primarily train the models using lagged volatility information. I do not use the potential of future market data for training and prediction. That is, the three models for predicting continuous data are: (1) Lasso regression; (2) Ridge regression; and (3) Decision Forest regression. The Lasso (Least Absolute Shrinkage and Selection Operator) is a linear regression technique that helps to prevent overfitting by shrinking the coefficients of less important features to zero. The method introduces a tuning parameter called lambda $\left(\lambda \right)$ that controls the amount of shrinkage applied to the coefficients. The Lasso model minimises the sum of the squared residuals (difference between predicted and actual values) and the penalty term, which is the absolute value of the sum of the coefficients multiplied by lambda. Increasing the lambda value will shrink more coefficients to zero, resulting in a simpler model with fewer features.

Ridge regression is a linear regression technique that helps to prevent overfitting by adding a penalty term to the sum of squared residuals. This penalty term is the sum of the squared coefficients multiplied by a tuning parameter called lambda $\left(\lambda \right)$, also known as the regularization strength. The Ridge model minimises the sum of the squared residuals and the penalty term, where λ controls the trade-off between the model's fit to the data and the coefficients' magnitude. By increasing the value of lambda $\left(\lambda \right)$, the Ridge model reduces the magnitude of the coefficients towards zero, resulting in a simpler model with smaller coefficients. Cross-validation, like Lasso, also determines the value of lambda in Ridge regression. I chose the best model with the lowest mean squared error on the validation set.

Decision forests, also known as random forests, are an ensemble learning method that combines multiple decision trees to make predictions. Hyperparameters are adjustable parameters that control the behavior of the decision forest model. Techniques such as grid search, random search, or Bayesian optimization, which evaluate different combinations of hyperparameters on a validation set to find the best-performing combination, can optimize hyperparameters for decision forests. Decision trees use data-based splitting rules to segment the data into subsets. Assign the average of the target variable within a subset as the prediction for all observations that fall inside that subset. Recursive binary splitting splits a sample into segments to implement a single decision tree. This iterative approach determines where and how to split the data based on what leads to the lowest residual sum of squares (RSS). Single-decision trees can have low, non-robust predictive power and suffer from high variance. Random decision forests, which offer performance improvements by combining results from groups, or "forests", of trees, can overcome this. In summary, the random decision forest algorithm randomly selects predictors as potential candidates for data splitting. Constructs a decision tree from a bootstrapped training set. Repeats the decision tree formation for a specified number of iterations. Averages the results from all trees to make a final prediction.

Lasso and Ridge regression aim to reduce prediction variances using a modified least squares approach. Let us look a little more closely at how this works. Recall that ordinary least squares estimate coefficients by minimizing the residual sum of squares (RSS):

Penalized least squares estimates coefficients using a modified function:

where $\lambda$ is the tuning parameter and $\lambda {J}_{2}$ is the penalty term. The penalty term for the Lasso regression is (L1): $\lambda \sum\nolimits_{j = 1}^{p}\left|{\beta }_{j}\right|$, and for the ridge regression (L2): $\lambda \sum\nolimits_{j = 1}^{p}{\beta }_{j}^{2}$.

2.4.3.   Artificial intelligence-neural networks

An artificial neural network fits a non-linear model to available data. In a timeseries model, the input lags relate to the first (second) hidden layers, or neurons. Note the lags, avoiding look-ahead bias. Neurons in the final hidden layer calculate the output dimensions and forecasts. Additionally, I use the RELU activation function in both C(R)NN (convolutional/recurrent neural network) and LSTM (long short-term memory) neural networks in this paper⁵. Figure 1 illustrates a simple one-layer neural network with three neurons that predicts step-ahead volatility: Features (lags) Neurons Step-ahead prediction Simple one-layer neural network with S neurons and biases (b) and an activation function, i.e., the softmax function.

⁵ See Hull (2021).

3.   WTI oil front month prices and price movements 2010–2024 (3)

3.1.   Data and stationarity

We ensure that the means, variances, and covariances (rather than the entire distribution) are time independent by enforcing weak stationarity. In other words, if that holds for all t, then a process {y_t} is weakly stationary $E\left\{{y}_{t}\right\} = \mu \le \infty, V\left\{{y}_{t}\right\} = E\left\{{\left({y}_{t}-\mu \right)}^{2}\right\} = {\gamma }_{0} < \infty$as well $\mathit{cov}\left\{{y}_{t}, {y}_{t-k}\right\} = E\left\{\left({y}_{t}-\mu \right)\left({y}_{t-k}-\mu \right)\right\} = {\gamma }_{k}, k = \mathrm{1, 2}, 3, ...$ Future observations are impacted by a shock to a stationary autoregressive process of order 1 (AR(1)) in a decreasing amount. Table 1 summarizes the characteristics of the price movement series. The data show a negative drift, suggesting a non-increasing WTI oil price level over a twelve years period (2014–2024(3) and 40,600 observations). The standard deviation on an hourly basis is 0.56, suggesting a yearly volatility of approximately 25%. The maximum (minimum) mean is 17 (−15). The kurtosis is high (101.7), indicating many observations around zero, followed by a few outliers. The skewness is positive (1.4), showing an asymmetric distribution for the movement densities. The high kurtosis likely leads to the strong rejection of the Cramer-Mises normality test (117.3). Serial correlation (Q(12)) reports serial correlation, while the Breusch-Godfrey LM statistic does not. Heteroscedasticity is present, as evidenced by the Q² (12) (1,327) and ARCH (12) (1,107) test statistics. These series are stationary (ADF (−184) and Phillips-Perron (−184)). However, the series reports some form of data dependence (BDS Z statistic > 30). The RESET test statistic (Ramsey, 1969) reports minor coefficient instability (14.9). Figure 2 shows the series of continuous movements. Overall, the series seems like usual financial market data. In the movement equations, I also tested with breaking trends, but my findings did not support trend breaks. A popular notion for assessing risk is value at risk (VaR), and Table 1 shows the 2.5% and 1% VaR numbers for market participants. Figure 3 displays the density important for non-normality.

Figure 4 reports correlogram for ordinary returns (top) and absolute returns (bottom) series. Finally, in Figure 5, a simple forecasting of normal, equal weighted, and GARCH(1, 1) volatility models are forecasted. The volatility seems well described using a GARCH(1, 1) model, showing heterogeneous variables/volatility.

3.2.   The SNP density projection

The conditional density, which completely describes how prices move, is, of course, the most important statistical variable. I fit the semi-nonparametric (SNP) model using conventional maximum likelihood and a model selection strategy that determines the appropriate order of expansion (BIC). I compute the Schwarz Bayes information criterion $BIC = {s}_{n}\left(\widehat{\theta }\right)+\left(\frac{1}{2}\right)\left(\frac{{p}_{p}}{n}\right)\mathit{log}\left(n\right)$ (Schwarz, 1978), preferring small criterion values.

Table 2 reports the maximum likelihood (ML) estimates of the parameters for the BIC-optimal SNP density models. First, for the mean, the intercept $\left({\eta }_{7}\right)$ is insignificant and the serial correlation $\left({\eta }_{8}\right)$ is negative (significant) (negative dependence). Second, the conditional variance coefficients $\left({\eta }_{9}-{\eta }_{13}\right)$ are all strongly significant. The variance parameters' significance suggests conditional heteroscedasticity for the series. The P&Q companion matrix's greatest eigenvalue for the conditional variance function is 1.006. The spline specification reports mean reversion of the conditional variance even though the P and Q companion matrix is 1.006. Last, the BIC favors the hermite function coefficients $\left({\eta }_{1}-{\eta }_{6}\right)$ up to the sixth polynomial lag expansion, which indicate modifications in the parametric model. Therefore, it is evident from the hermite result that there may be deviations from the traditional normally distributed and parametric conditional models. Asymmetry between positive and negative movements$\left({\eta }_{12}\right)$ seems not to be present. The quadratic and leverage reports corroborate the previously unreported finding (not reported). The high Q[1,1] $\left({\eta }_{11}\right)$ suggests volatility persistence.

The SNP projection densities give access to the conditional mean and volatility densities reported in Figure 6. Figure 7 shows densities for the same conditional mean and volatility. The figure reports a right skewed conditional volatility density (log-normal), while the mean displays leptokurtosis. Furthermore, Figure 8 allows for the generation and reporting of one-step-ahead mean densities based on lagged values. Any length of simulation path (bootstrapping) is available. All these results are in full compliance with the statistics from Table 1. Finally, for the GMM moment estimations, Table 3 reports residual statistics. The residual test statistics are closer to a normal distribution. However, the Cramer-von-Mises test for normality is significant. Furthermore, the BDS Z-statistics for m = 2 and 3, as well as the ARCH (12) test statistics, indicate data dependence. Most likely, the BDS and ARCH signals remain heteroscedastic.

3.3.   Stochastic volatility (EMM) for WTI oil front month prices

The parameter vector for the EMM model described in section 2.3, is reported in Table 4 below.

I use the efficient methods of moments (GMM-MCMC) to estimate the unobserved state vector, applying the Gallant and Tauchen (2010) method. In short, the EMM is a systematic approach generating moment conditions from the generalized method of the SNP moments estimator. The test-statistic ${\chi }^{2}\left(5\right) = 4.765$ with associated test statistic 0.45 from Table 4 (bottom) measure the extent of SV model failure⁶. A nonlinear Kalman filtering technique brings the unobserved vector back to reality. Figure 9 (11,500 observations) reports the unobserved vectors for the sub-period 2020–2024(3). The volatility appears to be divided into two factors: One factor that is very persistent (slow men reversion) and a strongly mean-reverting factor. Table 5 reports the same figures' characteristics in numerical form.

⁶ Degrees of freedom (5) is the number of SNP model parameters minus SV-model parameters –1.

The re-projected volatility, which is skewed to the right and non-normal, reports a yearly average of 22.1%. Both Phillips-Perron and Augmented Dickey-Fueller confirm the stationary nature of the data. Data dependence appears to be present, with a clearly higher value for volatility factor 1 (V₁) than for factor 2 (V₂). The correlogram for the factor volatilities in Figure 10 clearly suggests a strong serial correlation in the two series and in particular for V₁ og exp(V₁ and V₂). Moreover, the Brock, Deckert, and Scheinkman (BDS) test statistic for the correlation integral from 2 to 5 (m) suggests significant data dependence.

4.   Forecasting re-projected volatility factors

4.1.   Classical OLS static forecasting

The WTI oil Hourly Volatility Ordinary OLS Step-Ahead Fit Measures demonstrate a data dependence that suggests predictability through serial correlations. Table 6 reports the classical step-ahead forecasts using only endogenous variables (lags), and provide numbers for the volatility fit. The Theil's covariance measure is quite high. Factor 1 (V₁) has a 97.9% Theil covariance fit, indicating that it is clearly predictable. Factor 2 (V₂) report a lower Theil covariance measure of 70.8%. The re-projected volatility (exp(V₁ + V₂)) reports a close to factor 1 (V₁) result of 97.8%. Therefore, it appears that I can project the volatility. However, trading volatility from these factors can become challenging, necessitating a close watch on other factors to prevent sudden losses from the mean reversion process. Furthermore, when considering portfolio theory and systematic risk compensation, it's crucial to keep in mind the significance of covariance.

4.2.   Machine learning regression techniques and step ahead forecasting

Table 7 reports the prediction results from machine learning (ML) regression techniques for the volatility factors (V₁ (left), V₂ (middle) and exp(V₁ + V₂) (right column). Figure 11 contains the prediction plots for the same ML regression models. For the reprojected volatility, the Lasso model with a $\lambda = 0.05$ (L1 penalty is marginally positive) seems to report the lowest RMSE, followed by the Decision Forest model. Furthermore, the Theil covariance measure is clearly highest for the Lasso model, with $\lambda = 0$ (99.8%). The Ridge model with $\lambda = 0.1$ follows closely (99.8%). Generally, the Lasso model with lambda regularization serves as a powerful feature selection technique, enabling the construction of accurate and interpretable regression models. When the data contains numerous correlated features, people often use ridge regression. The regularization term shrinks the coefficients towards each other, which can improve the stability and interpretability of the model. Ridge regression with lambda regularization is therefore a useful technique for reducing overfitting and improving the performance of linear regression models. In this paper the $\lambda$ uses values from $0\le \lambda < 1$. I use hyperparameters to tune and improve the decision forest models.

Confidence intervals are not available using machine learning (ML) techniques. However, Table 6 reports fit measures for all regression models, which are all comparable with the static OLS forecasts. For volatility factor 1 (V₁), the regression techniques Lasso (λ = 0), Lasso $\left(\lambda = 0.05\right)$, Ridge $\left(\lambda = 0.1\right)$, and Decision Forest report RMSEs of 0.090, 0.101, 0.090, and 0.091, respectively. Furthermore, for factor V₁, the Lasso $\left(\lambda = 0\right)$, Lasso $\left(\lambda = 0.05\right)$, Ridge $\left(\lambda = 0.1\right)$, and Decision Forest regressions report a Theil covariance measure of 99.3%, 74.0%, 99.3%, and 96.9%, respectively. Note however, that the classical static forecasts (OLS) in Table 6 report a Theil covariance fit of 99.3%, approximately equal to Ridge regression. The covariance values for volatility factor 2 (V₂) are 92.8% for Lasso $\left(\lambda = 0\right)$, 3.4% for Lasso $\left(\lambda = 0.05\right)$, 92.7% for Ridge $\left(\lambda = 0.1\right)$, and 30.5% for Decision Forest regression. These values are clearly lower than those for factor 1 (V₁). Hence, neither the OLS nor the ML techniques seem to report good fits for factor 2 (V₂). For the re-projected volatility, the Lasso $\left(\lambda = 0\right)$, Lasso $\left(\lambda = 0.05\right)$, Ridge $\left(\lambda = 0.1\right)$, and Decision Forest regression techniques report covariance measures of 99.8%, 97.8%, 99.8%, and 95.4%, respectively. The regression results for the reprojected volatility is clearly closer to factor 1 (V₁) than factor 2 (V₂).

The improvements made using ML are marginal, so it is important to evaluate ML techniques in light of the use of computer resources. Moreover, thorough tuning of the hyperparameters for the ML techniques may improve performance. Consequently, based on the regression results, the machine learning regression models seem to expand the possibilities for volatility predictions. For the re-projected volatility $\left(\mathit{exp}\left({V}_{1}+{V}_{2}\right)\right)$, the ML Lasso regression techniques $\left(\lambda = 0\;and\;0.05\right)$and ridge $\left(\lambda = 0.1\right)$ report overall best fits $\left(\approx 98-99\%\right)$. In fact, the Lasso and Ridge regression models report a Theil covariance measure that is 5% better than the Decision Forest regression model.

In ML, over-(under-)fitting is important to avoid, referred to as the bias-variance trade-off. For all methods, hyperparameters are used for tuning and therefore minimize the trade-off. It will therefore be important to use hyperparameters to tune and improve the models. I can categorize the error of a machine learning model into two major categories: Bias and variance (refer to the Theil measure). The error that occurs when I fit a simple model to a more complex data-generating process. A model with high bias will underfit the training data. When I apply our model to a new dataset that it has not seen, the expected prediction error occurs. A model with a high variance will typically overfit the training data, resulting in lower training set errors but higher errors on any data not used for training. The model validation process in the previous section works when I have large datasets. When data is scarce, I must resort to a technique known as cross-validation. The purpose of cross-validation is to provide a better estimate of a model's ability to perform on unseen data. It provides an unbiased estimate of the generalization error, especially in the case of limited data. There are many reasons I may want to do this: (1) To have a clearer measure of how the model performs; (2) to tune hyperparameters; and (3) to make model selections. The idea behind cross-validation is simple: rather than training our models on one training set, I train them on multiple subsets of data. I follow the basic steps of cross-validation: (1) Divide the data into portions; (2) train our model on a subset of these portions; (3) test our model on the remaining subsets of the data; and (4) repeat steps 2–3 until I have trained and tested the model on the entire dataset. Average the model performance across all iterations of testing to get the total model performance.

4.3.   Neural networks (CNN/RNN/LSTM)

Figure 12 reports a neural network (CNN/RNN) for V_1t. The fit is satisfactory, and the one-step-ahead predictions have a low MSE (0.09). Figure 13 shows that V_2t has a much higher MSE (0.6). Hence, the neural networks report satisfactory forecasts for V_1t. However, V_2t seems difficult to forecast for neural networks. Note that V_1t seems to be the most important factor for the accumulated for the reprojected volatility (exp(V_1t + V_2t).

5.   Summary and conclusions

In this paper, I have modelled and successfully estimated a multifactor stochastic model for the conditional mean and variance for the WTI oil price movements, period 2014 to 2024(3). This multifactor model aligns with the efficient market hypothesis, which posits that relevant random information, accessible to all traders, drives financial market prices, rendering prediction impossible (akin to a random walk). The applications of option pricing, risk assessment, trading, portfolio optimization, and market surveillance show that volatility is a cornerstone of modern financial theory and practice. The analysis reports success with two volatility factors: one is extremely persistent, and the other strongly implies reverting. Furthermore, the latent volatility's data dependence suggests predictability. The persistent factor is predictable; the mean reversion factor is most likely not predictable. The interesting question now is whether to trade volatility for only the predictable, persistent factor. I leave the market implementation to future research using ML and AI methods.

Even though price processes in energy markets seem stochastic and are unpredictable, I can use observed previous prices and their variations to determine the time-dependent variance of forecast errors. For derivatives time series volatility is important for market participants in actual financial derivatives markets as well as synthetically for index and assets (B&S formula). The static predictions of the projected volatility for the WTI oil contracts show a Theil's inequality coefficient close to zero and a covariance portion of 99.8% using Ridge regression (ML) with a λ² = 0.1 (L2). Furthermore, using continuous prediction updates (i.e., hourly) may further improve these measures of fit. Also, V₁, the main and long-lasting factor for the expected volatility in WTI oil prices, has a Theil inequality coefficient that is very close to zero and a similar covariance portion of 99.3%. Furthermore, note that the reprojected volatility is strongly influenced of the persistent and predictable volatility factor V₁. Market participants can apply a continuous multifactor SV model with associated volatility trading strategies to achieve superior profit from the energy markets. Hence, extending our results to AI/ML for program trading and applying volatility forecasts may turn out to be warranted.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Data availability statement

The article uses a freely downloadable dataset (40546 hourly and 2443 daily observations).

The paper encloses the hourly WTI oil datasets for the period 2014–2024(3):

1. 001-High_Frequency_WTI_FM_ price_Movements_2014-2024(3).txt.

2. 001-Daily_Frequency_WTI_FM_ price_Movements_2014-2024(3).txt.

Conflict of interest

The author declares no conflicts of interest in this paper.