Co-Training Realized Volatility Prediction Model with Neural Distributional Transformation

The use of distributional transformations dates back to Fisher [16], where a Student's t variable was expanded as an infinite series of normal variables. In another early work, Wallace [39] proposed a simple approximation form that was found to work exceptionally well transforming between Gaussian and Student's t distributions [28].

A transformation can be used to improve the Gaussianity of a dataset that significantly deviates from a Gaussian distribution. This is often beneficial as Gaussianity occurs as an assumption underlying many statistical and econometric models.

The modern research on distributional transformation for economic data is largely based on Box and Cox [4]’s seminal work, which introduced a family of power transformations now called Box-Cox transformations. Box-Cox transformations include the identity transformation (i.e., f(x) = x) and the log transformation as two special cases. Box-Cox transformations are applied widely, but they only allow positive inputs. Yeo-Johnson transformations [41] were later introduced to generalize Box-Cox to include zero and negative inputs. Apart from power transformations, numerous alternatives have been proposed. Hoyle [23] summarized 18 different analytical transformations. Other choices include the Tukey g-and-h [38], the Lambert W function [19], and exponential functions [40].

While the improvement via a Gaussianizing transformation is frequently observed for cross-sectional data (i.e., a collection of samples at a given time), it does not apply to sequential data. Proietti and Lütkepohl [29] tested Box-Cox transformations on various real datasets with other transformations and found that Box-Cox transformations showed advantages only in a fifth of all cases. Taylor [35] investigated the use of Box-Cox transformations for RV prediction. Among several pre-defined choices of the parameter λ of the Box-Cox transformation, λ = 1/4 performed the best. However, as mentioned in the Introduction, the “correct” shape of the transformed distribution is influenced by the time-series prediction model. The optimal parameter in one scenario might be sub-optimal in other scenarios.

This study uses neural networks for the transformation. Simple neural-network transformations have been considered previously for time-series modeling. Snelson et al. [32] proposed using an additive mixture of tanh functions, a two-layer feed-forward neural network. The parameters in the neural network were restricted to be positive, so that the neural network remains invertible. Recent state-of-the-art techniques are with normalizing flows [13, 31, 34], which provide more powerful choices than the mixture-of-tanh approach in Snelson et al. [32].

Realized volatility, like other volatility measures, exhibits long-memory characteristics often referred to as “volatility clustering” [10]. Based on this finding, many models have been proposed to predict RV, among which two are especially popular. The first is the heterogeneous auto-regressive (HAR) model by Corsi [8], which uses a linear combination of past volatility terms defined over different time periods to predict future volatility. HAR is simple but performs well in practice; therefore, we use it as the prediction model in this paper. The other prevalent model is the auto-regressive fractionally integrated moving average (ARFIMA) [22], which utilizes fractional differencing operators to enable long memory.

It must be noted, however, that RV prediction is closely related to the fat tail problem in financial returns. Research on the fat tails in financial returns dates to Mandelbrot [26] and Fama [15], who noticed that substantial price changes occur more frequently than a Gaussian predicts; thus, the observed return distribution shows a larger kurtosis than a Gaussian (i.e., it has a fat tail). This observation was further validated on large-scale measurements of high-frequency price data [17, 20, 21, 27]. Volatility prediction models have been developed to account for the fat-tailed distribution of returns. But as we will show via baselines in Section 6, this implicit approach via volatility models is often insufficient to fully capture the fat tail.

The modeling of a transformed time series can be seen as a kind of latent-variable model. Realized volatilities are “observable” variables while the transformed values are seen as “latent.” Because of the long-memory nature of price time-series, the long-range dependence between the latent variables must be addressed. How to integrate financial econometrics methods into the latent-variable framework of machine learning remains a question.

Several latent-variable methods are considered relevant to our approach. The first is the Kalman filter [24], which assumes that all latent variables are Gaussian and governed by a linear Markovian process. A Kalman filter is estimated by making an initial guess of the earliest latent variable and then progressively gets refined by incorporating new observations at later timesteps.

The second is the hidden Markov models (HMMs) [3]. HMMs are more elaborate latent-variable models that have been used widely in various fields [5]. Kalman filters can be viewed as special cases of HMMs, even though they were studied in different research streams for a long time [5].

Nevertheless, the Markovian nature of Kalman filters and HMMs makes it hard to capture the long memory in RV time-series data [8]. In contrast, a simpler HAR can easily access a timestep far back (e.g., 22 timesteps previous to current), which is much more straightforward than HMMs in incorporating the long memory and has been found to work well.

In this paper, we use normalizing flows (NFs) to assist HAR models in capturing the complex RV distribution. NFs were proposed for variational inference [31], which is a method frequently used for latent-variable modeling. Thus far, normalizing flows have been shown to be effective in various applications, such as density estimation [34], generative modeling [13], variational inference [31], noise modeling [1], and time-series modeling [12].

In a recent work by Deng et al. [12], normalizing flows were used to transform a time series into a Brownian motion. As Brownian motions are Markov processes, the model by Deng et al. [12] can be seen as a hidden Markov model in which the latent Gaussian variables are estimated using a normalizing flow. In contrast, the latent process in this work must capture the long-memory phenomena observed in real RV data [8].

In this paper, we study the task of predicting RVs with a learnable monotonic transformation $f_{\boldsymbol {\theta }}:\mathbb {R}\rightarrow \mathbb {R}$ parameterized by θ. At every timestep t, realized volatility x_t is transformed into z_t = f_θ(x_t), which follows a latent distribution. Instead of predicting x_{t + 1} directly, z_{t + 1} is predicted as $\hat{z}_{t+1}$ using z_{≤ t}, and the realized volatility x_{t + 1} is then estimated from $\hat{z}_{t+1}$. The entire prediction procedure is as follows:

Using the transformation f_θ has multiple advantages. First, the transformed realized volatilities better approximate Gaussian distributions, an implicit assumption that underlies many prediction models. Thus, the prediction model is estimated with less bias. Second, for a linear prediction model y_β, the use of nonlinear transformation f_θ enables capturing nonlinear temporal dependency.

In this paper, we propose implementing the invertible transformation f_θ with normalizing flow, a deep-learning technique. Normalizing flow refers to invertible neural networks for which the inverse can be efficiently computed. Specifically, we use neural ordinary differential equations (NODE) [6], a special kind of normalizing flow.

A NODE is a first-order ordinary differential equation parameterized by a neural network:

For the prediction model y_β, we use the heterogeneous auto-regressive (HAR) model [8]. HAR is one of the most common models for RV prediction.

HAR is a linear regression model that incorporates RV components during three different periods, as follows:

Although this paper focuses on HAR, the proposed method can be applied to prediction models other than HAR. Because normalizing flow specifies a wide family of continuous functions, it requires little knowledge about the prediction model to effectuate the prediction in Formulas (3) through (5). The universality of our method with different prediction models remain a future work.

Let X = [X₁, X₂, ⋅⋅⋅] and Z = [Z₁, Z₂, ⋅⋅⋅] denote the random processes underlying the raw and the transformed realized volatilities, respectively, where Z_s = f_θ(X_s) (s = 1, 2, ⋅⋅⋅). The normalizing flow f_θ is optimized subject to a maximum-likelihood objective as follows:

In this work, the dependence of Z_t on X_t is deterministic via the normalizing flow, while that of X_{t + 1} on Z_t can be random. Hence, the sampling procedure of Z ∼ P(· ∣X, θ⁽ⁱ⁾) in Formula (13) is simply to apply f_θ to every timestep of X₁, X₂, ⋅⋅⋅, and the expectation average is taken over the RV time series of multiple stocks. On the other hand, the term within the expectation operation of Formula (13) decomposes into the following due to the change-of-variable theorem:

Notice that in Formula (15), P(Z_{s + 1}∣Z_{≤ s}) represents the prediction model y_β applied to the latent process, and it does not specify how the parameters β should be estimated. Therefore, optimization of the parameters in the prediction model (i.e., β) can be effectuated via a different objective function from that in Formula (13). This is important because HAR and many other financial time-series models are conventionally estimated with a least-square objective function, which imposes a weaker assumption on the distribution of residuals and often shows better robustness. The probability density P(Z_{s + 1}∣Z_{≤ s}) is obtained via an estimation, which is detailed in Section 4.2.

Hence, the co-training procedure of f_θ and y_β is via iterative updates as follows: at each iteration, the sampling step Z ∼ P(· ∣X, θ⁽ⁱ⁾) in Formula (13) is conducted as a random sampling of X_{≤ t} from the RV time series of stocks and then transforming X_{≤ t} into Z_{≤ t}. Using the samples of Z_{≤ t}, the parameters of the prediction model (i.e., β) can be estimated via any adequate method; here for HAR, we used the Python package arch ¹ to estimate β. Then, Z_{t + 1} is predicted, and the expectation of log-likelihoods log P(x_{t + 1}∣Z_{≤ t}, θ) is calculated using Formula (15). Finally, the $\mathop{arg\,max}$ operation in Formula (13) is approximated by a mini-batch gradient ascent step over the expectation of the log-likelihoods.

The conditional probability density P(Z_{s + 1}∣Z_{≤ s}) is determined by the prediction model y_β, which produces an estimate $\hat{Z}_{s+1} = y_{\boldsymbol {\beta }}(Z_{\le s})$ based on previous latent states. Thus, the conditional probability density is equivalent to the density of the prediction residual:

The distributions of residuals are usually unknown in practice. In this paper, we assume they are Gaussian. This assumption could be inappropriate if the prediction is with the raw RV or when the transformation f_θ is as simple as a Box-Cox transformation [29]. Nevertheless, we find empirically that a Gaussian distribution is a good approximate for the residuals after the normalizing flow is optimized, which could be due to the larger capacity of the normalizing-flow neural network. In addition, we assume the residuals to be homogeneous across time (i.e., ɛ_s follows the same Gaussian distribution for every timestep s). Using a more advanced time-series model would improve the homogeneity of the residuals, which is left to future work.

Therefore, the residuals ɛ_s ∼ N(μ, v) (s = 1, 2, ⋅⋅⋅, t), where μ and v are the mean and variance of the Gaussian distribution. μ is set to zero² and v is estimated as the empirical variance: $v = \left(\sum _{s=1}^t (\varepsilon _s - \mu)^2 \right) / (t-1)$. Notice that v involves the whole period 1, 2, ⋅⋅⋅, t, but is used for calculating P(ɛ_{s + 1}∣Z_{≤ s}) for s < t, which means the calculation involves “further information” at timestep s < t. Nevertheless, it does not invalidate our approach because the calculation of P(ɛ_{s + 1}∣Z_{≤ s}) is only required in parameter estimation; that is, making out-of-sample prediction on Z_{t + 1}, Z_{t + 2}, ⋅⋅⋅ does not involve the calculation of P(ɛ_{s + 1}∣Z_{≤ s}).

The proposed method is evaluated on a dataset of 100 major stocks listed on the New York Stock Exchange. The dataset covers a continuous time period of 480 trading days from December 1st, 2015, to October 25th, 2017. The daily RV values were calculated according to Formula (1) for every stock. Thus, 100 RV time series of 480 timesteps were acquired. In computing the daily RV, we followed Corsi [8] and did not include the overnight price variation. That is, we did not include $p_{t+1,1} - p_{t, n_t}$ in the summation of Formula (1), where n_t denotes the number of high-frequency prices within day t as mentioned in Section 3.1.

In this paper, the RV of a day is calculated with high-frequency prices that are recorded every five minutes. That is, p_{t, i + 1} was taken five minutes later than p_{t, i}. The value of high-frequency price is calculated as the mid-price [33] of the order book, which is an average of the highest bid price and the lowest ask price weighted by their volumes, as follows:

The RV time series were split into training, validation, and test sets chronologically, at a ratio of 300, 60, and 120 days.

We conducted the experiments in a univariate time-series setting by viewing all 100 RV time series as independent samples drawn from the same random process X₁, X₂, ⋅⋅⋅. Thus, this accommodates the notion of taking an expectation over IID samples in Formula (13). The normalizing flow f_θ was optimized on the training set (i.e., the first 300 days of the RV history of all stocks). Using an Adam [25] optimizer, we conducted 200 iterations of parameter updates and evaluated the model on the validation set at every five training iterations, thus producing 40 different snapshots of parameters. The snapshot that achieved the highest log-likelihood in Formula (15) was regarded as the best. In the following, we report the performance of this best snapshot on the test set.

For the neural network NN_θ used in the NODE (see Formula (6)), we used a simple multi-layer perceptron with two hidden layers. Each hidden layer had four hidden units. Formally,

We considered different baseliens by varying the transformation f_θ and the distribution for the residuals ɛ_s (see Formula (9)).

For alternative transformations, we considered the following:

• transformation: f_θ(x) = x, which is equivalent to not applying a transformation.
  
• transformation [39]:
  \begin{equation} f_{\boldsymbol {\theta }}(x)= \pm \frac{8d+1}{8d+3}\sqrt {d\log (1+|x|^2/d)}, \end{equation}
  (21)
  where θ = {d}.
  
• transformation [41]: a generalization of the Box-Cox transformation [4] to allow non-positive inputs, defined as follows:
  \begin{equation} f_{\boldsymbol {\theta }}(x) = \left\lbrace \begin{array}{ll}((x+1)^\lambda - 1) / \lambda & \text{if}~\lambda \ne 0, x\ge 0, \\ \log (x+1) & \text{if}~\lambda =0, x\ge 0, \\ -((-x+1)^{(2-\lambda)} - 1) / (2 - \lambda) & \text{if}~\lambda \ne 2, x< 0, \\ -\log (-x+1) & \text{if}~\lambda =2, x< 0, \end{array}\right. \end{equation}
  (22)
  where θ = {λ}.
  
• transformation [32]: a simplest neural network with a single hidden layer, defined as follows:
  \begin{equation} f_{\boldsymbol {\theta }}(x) = \sum _{i=1}^k u_i \tanh (v_i x+b_i), \end{equation}
  (23)
  where $\boldsymbol {\theta }=\bigcup _{i=1}^k \lbrace u_i, v_i, b_i\rbrace$. For f_θ to be invertible, u_i, v_i (i = 1, ⋅⋅⋅, k) are restricted to be positive.
  

As for the distribution of the residuals, we tested the generalized Student's t distribution (in addition to the Gaussian by default), which is defined by the following probability density function with mean μ, variance v, and d degrees of freedom:

Table 1 summarizes the overall results over the dataset of 100 stocks acquired with different transformations. Each row represents a transformation. Our transformation via normalizing flow is shown at the bottom of the table. The third and fourth columns show the scores concerning HAR prediction accuracy. RMSE was used as the primary evaluation metric to follow the previous works [8]. A smaller RMSE implies a higher precision.

The third column of Table 1 reports the average RMSE (mean value) over the 100 stocks. One-tailed t-tests were conducted to assess the statistical significance of our method (NODE with τ = 0.25) compared with the other transformations, assuming the null hypothesis that the mean improvement in RMSE is equal to zero. Asterisks (*) in the third column indicate statistical significance: * (p < 0.05), ** (p < 0.01), *** (p < 0.001). The fourth column shows the proportion of stocks on which a method achieved the lowest RMSE, or “percentage of the best.” The best scores are indicated in bold, and the second-best scores are underlined. For a robustness test, we also examined NODE with τ = 5.00 which has excessive complexity as a transformation; we omitted the p-value and the “percentage of the best” value as it serves merely as a robustness test.

Our approach (bottom row) achieved the lowest RMSE score at 0.5620, a clear improvement of 0.0075 compared with the identity transformation that scored 0.5695. The second-best at 0.5627 was the Yeo-Johnson transformation, which is a common way of preprocessing non-Gaussian data in practice; nevertheless, our method still outperformed the Yeo-Johnson transformation. Compared with the baseline transformations, the improvements of NODE in RMSE are statistically significant at the 5% (Yeo-Johnson), 1% (Wallace's and tanh  transformations), and 0.1% (Identity) levels. With respect to the “percentage of the best,” the advantage of NODE is more evident. NODE achieved the lowest RMSE on 46 of the 100 stocks; in comparison, the Yeo-Johnson transformation won on only two stocks.

The significant margin between NODE and Yeo-Johnson on the “percentage of the best” is interesting, in contrast to their relatively close performance in average RMSE. This indicates that on many stocks, NODE and Yeo-Johnson transformed realized volatilities into close distributions, but NODE performed slightly better in most cases. While power transformations like Yeo-Johnson are careful choices made by practitioners after decades, NODE was learned from data. However, NODE “discovered” the power transformations were the best choices and further improved over them.

The transformations denoted by tanh(k) (k = 1, 5, 10) are the simplest feed-forward neural networks, and k represents the number of hidden units. When k was increased from 1 to 5, the average RMSE was improved from 0.5642 to 0.5636. However, when k was further increased to 10, the score degraded to 0.5638.

The degradation at large k values might be caused by the increased risk of overfitting to the training set. The tanh transformations with a large k or NODE with a large τ had a stronger approximation capability, but they overfitted to the training set. Evidence of this is that on 29 stocks, simpler transformations, including Identity or Wallace's, had the best RMSE. On these stocks, the data distribution of the test set shows inconsistency with the training set, which is common in a financial market that is a non-stationary system. Such overfitting is also seen for NODE when τ was increased from 0.25 to 5.00 when RMSE increased from 0.5620 to 0.5630.

As an ablation analysis, we evaluated the robustness of the presumption made in Section 4.2 that the residuals of HAR (i.e., ɛ_s in Formula (9)) follow a Gaussian distribution. Instead of assuming that ɛ_s ∈ N(μ, v), the ablated version assumes a generalized Student's t distribution with an additional parameter d denoting the degrees of freedom. Note that d → ∞ recovers a Gaussian distribution.

Figure 2 shows the average RMSE scores (vertical axis) when the residual distribution was assumed to be either Gaussian (white bars) or generalized Student's t (grey bars) grouped by the transformation f_θ. When the residuals are assumed to follow a generalized Student's t distribution, the probability density in Equation (16) is determined using the PDF of the Student's t distribution, instead of the Gaussian PDF. Consequently, the transformation parameters θ that are estimated yield an RMSE score that differs from what would be expected under a Gaussian residual assumption.

In Figure 2, the white bars are the values in Table 1 in the third column. Except for the Identity transformation, we observed increased RMSE scores when the presumed distribution was a Student's t, which means a decrease in prediction accuracy. The largest increase in RMSE is observed with the tanh(1) transformation, implying tanh(1) to be sensitive to the distributional assumption. In contrast, NODE still achieved the smallest RMSE score at 0.5643, the same as with a Yeo-Johnson transformation. The results in Figure 2 suggest the robustness of our approach even under a generalized Student's t distribution assumption for residuals that is not optimal.

Figure 3 provides a comparison between the distribution of residuals ɛ_s defined in Formula (9) (vertical axes) and the standard normal distribution (horizontal axes) for the stock “BABA” (i.e., the Q-Q plot). Each plot represents a different transformation listed in Section 5.3. In each plot, a point represents the residual at a timestep; the vertical axis shows the residual values, and the horizontal axis shows the theoretical percentiles if the residuals follow a standard normal distribution. The straight red line is a linear fit of the points, and a better fit implies a higher degree to which the residuals follow a Gaussian distribution.

As seen in the first plot of Figure 3 (Identity transformation), the residuals of the raw RV deviate from the linear fit at large positive quantiles. This shows how HAR fails to fully capture the skewness in RV time series. In contrast, with a transformation, the Gaussianity of the residuals was largely improved. This improvement is visible for all the transformations shown in Figure 3. Among the transformations tested, the Yeo-Johnson, the tanh(5), and NODE produced impressive Gaussianity, and they are indistinguishable within a large range.

The Gaussianity of the residuals is further quantified using two metrics, as shown in the two right-hand columns in Table 1. The first metric is the R² score of the linear fit in the Q-Q plots of Figure 3, and the other is the skewness (i.e., the third moment) of the residuals. The two right-hand columns present the average R² and the average skewness over the 100 stocks.

Consistent with the observation of Figure 3, without a transformation (first row), HAR produced residuals with a low R² at 93.22% and a high skewness at 1.7883, implying poor Gaussianity. The Gaussianity was greatly improved by using a transformation. R² was improved to 99.14%, and skewness was reduced to 0.2141 with a NODE having τ = 5.00 (second to last row).

When R² exceeds 99%, an improvement in R² or a reduction in skewness did not translate into an improvement in out-of-sample RMSE. This corresponds with our previous conjecture on overfitting. Nevertheless, even with the same level of R² and skewness, NODE with τ = 5.00 still outperformed the tanh(5) transformation in out-of-sample RMSE.

Figure 4 provides a visualized comparison between the transformations on the real line. The horizontal axis represents the input $x\in \mathbb {R}$ to the transformation f_θ, and the vertical axis shows the output (i.e., f_θ(x)). Note that the realized volatilities were z-score standardized in preprocessing, which produced negative x. For visualization, we also linearly rescaled the transformations so that f_θ(− 2) = −2 and f_θ(2) = 2. As HAR is a linear model with a bias term, the estimation and prediction with HAR are invariant under linear rescaling to the transformation.

Each curve in Figure 4 represents a transformation. The straight green line represents the identity function, and the dotted red plot represents the NODE transformation. The NODE, Yeo-Johnson (in blue), and tanh transformations (in purple) all showed a concave shape: their slopes gradually decrease as x increases. This corresponds with our expectation that larger realized volatilities are calibrated more than smaller realized volatilities to eliminate skewness.

At large x, NODE showed almost perfect consistency with the Yeo-Johnson transformation. However, a discrepancy is seen at small x where NODE grows linearly, similar to the Identity function. In other words, NODE viewed it “unnecessary” to apply nonlinear calibration to small realized volatilities. With the Yeo-Johnson or tanh transformations, such local linear growth is not possible, as power and tanh transformations are defined globally over the real line.

Figure 5 presents the Pearson's correlation coefficients between the transformations on the RMSE improvements over the Identity transformation, measured for the 100 stocks. The RMSE values were obtained through out-of-sample evaluations, and the mean value for each transformation has been reported in Table 1, in the third column.

These correlation results align closely with the graph of f_θ(x) in Figure 4. Notably, the Yeo-Johnson and the tanh (k = 5) transformations exhibit a correlation coefficient of 0.998. This is unsurprising as they appear nearly identical in shape in Figure 4 (the purple and blue plots). Furthermore, NODE also demonstrates a high correlation with the Yeo-Johnson transformation, which validates our conjecture in Section 6.1 that NODE has “discovered” (while also outperformed) power transformations from the data.