# 3D/2D regularized CNN feature hierarchy for Hyperspectral image
classification

Muhammad Ahmad, Manuel Mazzara, and Salvatore Distefano M. Ahmad is with the
Department of Computer Science, National University of Computer and Emerging
Sciences, Islamabad, Chiniot-Faisalabad Campus, Chiniot 35400, PakistanM.
Mazzara is with Innopolis University, Innopolis 420500, RussiaS. Distefano is
with Dipartimento di Matematica e Informatica—MIFT, University of Messina,
Messina 98121, Italy

###### Abstract

Convolutional Neural Networks (CNN) have been rigorously studied for
Hyperspectral Image Classification (HSIC) and are known to be effective in
exploiting joint spatial-spectral information with the expense of lower
generalization performance and learning speed due to the hard labels and non-
uniform distribution over labels. Several regularization techniques have been
used to overcome the aforesaid issues. However, sometimes models learn to
predict the samples extremely confidently which is not good from a
generalization point of view. Therefore, this paper proposed an idea to
enhance the generalization performance of a hybrid CNN for HSIC using soft
labels that are a weighted average of the hard labels and uniform distribution
over ground labels. The proposed method helps to prevent CNN from becoming
over-confident. We empirically show that in improving generalization
performance, label smoothing also improves model calibration which
significantly improves beam-search. Several publicly available Hyperspectral
datasets are used to validate the experimental evaluation which reveals
improved generalization performance, statistical significance, and
computational complexity as compared to the state-of-the-art models. The
source code will be made available at https://github.com/mahmad00.

###### Index Terms:

Beam-search; Label Smoothing; Hybrid Convolutional Neural Network (CNN);
Hyperspectral Images Classification (HSIC).

##  I Introduction

Hyperspectral Imaging (HSI) has been extensively utilized for many real-world
applications [1]; for instance, Crop monitoring, vegetation coverage,
precision agriculture, land resources, oil spills, water quality, meat
adulteration [2, 3], adulteration in house hold products such as color
adulteration detection in red chili [4], microbial spoilage and shelf-life of
bakery products [5].

Thus, HSI Classification (HSIC) has received remarkable attention and
intensive research results have been reported in the past few decades [6].
According to the literature, HSIC can be categorized into spatial, spectral,
and spatial-spectral feature methods [7]. The spectral feature can be labeled
as a primitive characteristic of HSI also known as spectral curve or vector
whereas the spatial feature contains the relationship between the central
pixel and its context which significantly uplifts the performance [8].

In the last few years, deep learning especially Convolutional Neural Network
(CNN) has received widespread attention due to its ability to automatically
learn non-linear features for classification, i.e., overcome the challenges of
hand-crafted features for classification. Moreover, it can jointly investigate
the spatial-spectral information and such models can be categorized into two
groups, i.e., single and two-stream, more information regarding single or two-
stream methods can be found in [9]. This work explicitly investigate a single-
stream method but in a combination of 3-D and 2-D layers.

Irrespective of the single or two-stream method, all deep learning frameworks
are sensitive to the loss which needs to be minimized [10, 11]. Several
classical works showed that the gradient descent to minimize cross-entropy
performs better in terms of classification and has fast convergence, however,
to some extent, leads to the overfitting. Several regularization techniques
such as dropout, L1, L2, etc., have been used to overcome the overfitting
issues together with several other exotic objectives performed exceptionally
well than the standard cross-entropy [12]. Recently, a work [13] proposed a
Label Smoothing (LS) technique that improves the accuracy significantly by
computing cross-entropy with a weighted mixture of targets with uniform
distribution instead of hard-coded targets.

Since then, LS has been known to improve the classification performance of
deep models [14]. However, the original idea was used to improve the
classification performance of only the inception model on ImageNet data [13].
However, since then, various image classification models (training models)
have used LS [15, 16]. Though the LS technique is widely used trick to improve
the classification performance and to speed up the convergence process,
however, it has not been much explored for HSIC, and above all, it has not
been much explored regarding when and why LS should work. This is due to, in
HSIC, sometimes models learn to predict the training samples extremely
confidently which is not good from a generalization point of view.

Therefore, this paper proposed an idea to enhance the generalization
performance of a hybrid CNN for HSIC using soft labels that are a weighted
average of the hard labels and uniform distribution over target labels. The
proposed method helps to prevent Hybrid CNN from becoming over-confident. We
empirically show that in improving generalization performance, LS also
improves model calibration which significantly improves beam-search. Several
publicly available Hyperspectral datasets are used to validate the
experimental evaluation which reveals improved generalization performance,
statistical significance, and computational complexity as compared to the
state-of-the-art 2-D/3-D CNN models.

##  II Proposed Methodology

Let us assume that the Hyperspectral data can be represented as $R^{(M\times
N)\times B^{*}}=[r_{1},r_{2},r_{3},\dots,r_{S}]^{T}$, where $B^{*}$ be the
total number of bands and $(M\times N)$ are the samples per band belonging to
$Y$ classes and
$r_{i}=[r_{1,i},~{}r_{2,i},~{}r_{3,i},~{}\dots,r_{B^{*},i}]^{T}$ is the
$i^{th}$ sample in the Hyperspectral Data. Suppose
$(r_{i},y_{i})\in(\mathcal{R}^{M\times N\times B^{*}},\mathcal{R}^{Y})$, where
$y_{i}$ is the class label of the $i^{th}$ sample. Later instead of using the
”hot encoding” in hybrid model, we introduce a smoothing technique
$\mu(y|r_{i})$. Thus the new ground truths $(r_{i},y_{i})$ would be:

$p^{\prime}(y|r_{i})=(1-\varepsilon)p(y|r_{i})+\varepsilon\mu(y|r_{i})$ (1)

$f(x)=\begin{cases}1-\varepsilon+\varepsilon\mu(y|r_{i})&if~{}y=y_{i}\\\
\varepsilon\mu(y|x_{i})&otherwise\end{cases}$ (2)

where $\varepsilon\in[0,1]$ is a weight factor, and note that
$\sum_{y=1}^{Y}p^{\prime}(y|r_{i})=1$. These new ground truths has been used
in loss function instead of hot-encoding.

$L^{\prime}=-\sum_{i=1}^{M\times N}\sum_{y=1}^{Y}p^{\prime}(y|r_{i})\log
q_{\theta}(y|r_{i})$ (3)

$L^{\prime}=-\sum_{i=1}^{M\times
N}\sum_{y=1}^{Y}\big{[}(1-\varepsilon)p(y|r_{i})+\varepsilon\mu(y|r_{i})\big{]}\log
q_{\theta}(y|r_{i})$ (4)

$\displaystyle L^{\prime}=\sum_{i=1}^{M\times
N}\bigg{\\{}(1-\varepsilon)\Big{[}-\sum_{y=1}^{Y}p(y|r_{i})\log
q_{\theta}(y|r_{i})\Big{]}+$ (5)
$\displaystyle\varepsilon\Big{[}-\sum_{y=1}^{Y}\mu(y|x_{i})\log
q_{\theta}(y|r_{i})\Big{]}\bigg{\\}}$

$L^{\prime}=\sum_{i=1}^{M\times
N}\Big{[}(1-\varepsilon)H_{i}(p,q_{\theta})+\varepsilon
H_{i}(u,q_{\theta})\Big{]}$ (6)

One can observe that each ground truth, the loss contribution is a mixture of
entropy between predicted distribution ($H_{i}(p,q^{\theta})$) and the hot-
encoding, and the entropy between the predicted distribution
($H_{i}(\mu,q^{\theta})$) and the noise distribution. While training,
$H_{i}(p,q^{\theta})=0$ if the model learns to predict the distribution
confidently, however, $H_{i}(\mu,q^{\theta})$ will increase dramatically. To
overcome this phenomenon, we used a regularizer $H_{i}(\mu,q^{\theta})$ to
prevent the model from predicting too confidently. In practice, $\mu(y|r)$ is
a uniform distribution that does not dependant on Hyperspectral Data. That is
to say $\mu(y|r)=\frac{1}{Y}$. The ultimate network is trained for $50$ epochs
using a mini-batch size of 256 without any batch normalization and data
augmentation. In a nutshell, the details of 3D/2D convolutional layers and
kernels are as follows: $3D\\_conv\\_layer\\_1=8\times 5\times 5\times 7\times
1$ i.e. $K_{1}^{1}=5,K_{2}^{1}=5$ and $K_{3}^{1}=7$.
$3D\\_conv\\_layer\\_2=16\times 5\times 5\times 5\times 8$ i.e.
$K_{1}^{2}=5,K_{2}^{2}=5$. $K_{3}^{2}=5$. $3D\\_conv\\_layer\\_3=32\times
3\times 3\times 3\times 16$ i.e. $K_{1}^{3}=3,K_{2}^{3}=3$ and $K_{3}^{3}=3$.
$3D\\_conv\\_layer\\_4=64\times 3\times 3\times 3\times 32$ i.e.
$K_{1}^{3}=3,K_{2}^{3}=3$ and $K_{3}^{3}=3$. $2D\\_conv\\_layer\\_5=128\times
3\times 3\times 64$ i.e. $K_{1}^{2}=3$ and $K_{2}^{2}=3$. Three 3D
convolutional layers are employed to increase the number of spectral-spatial
feature maps and one 2D convolutional layer is used to discriminate the
spatial features within different spectral bands while preserving the spectral
information. Initially, the weights are randomized and then optimized using
back-propagation with the Adam optimizer by using the loss function presented
in equation 6. Further details regarding the Hybrid CNN architecture in terms
of types of layers, dimensions of output feature maps and number of trainable
parameters can be found in [6].

##  III Experimental Settings and Results

The experiments have been conducted on two real HSI datasets, namely, Indian
Pines (IP), and Pavia University (PU). These datasets are acquired by two
different sensors i.e, Reflective Optics System Imaging Spectrometer (ROSIS)
and Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) [6]. The
experimental results explained in this work have been obtained through Google
Colab [17] which is an online platform to execute any python environment while
having a good internet speed to execute the code. Google Colab provides the
option to execute many versions of python, Graphical Process Unit (GPU), upto
358$+$ GB of cloud storage, and upto 25 GB of Random Access Memory (RAM).

In all the experiments, the initial size of the train/validation/test sets is
set to 25%/25%/50% to validate the proposed model as well as several other
state-of-art-deep models. The comparative methods include AlexNet [18], LeNet
[19], 2D CNN [20], and 3D CNN [21] models. Moreover, the following accuracy
metrics have been conducted to validate the claims made in this manuscript.
The accuracy metrics include Kappa ($\kappa$)111$\kappa$ is known as a
statistical metric that considered the mutual information regarding a strong
agreement among classification and ground-truth maps, Average (AA)222AA
represents the average class-wise classification performance, and Overall
(OA)333OA is computed as the number of correctly classified examples out of
the total test examples.All these metrics are computed using the following
equations.

$Kappa~{}(\kappa)=\frac{P_{o}-P_{e}}{1-P_{e}}$ (7)

where

$P_{e}=P^{+}+P^{-}$

$P^{+}=\frac{TP+FN}{TP+FN+FP+TN}\times\frac{TP+FN}{TP+FN+FP+TN}$

$P^{-}=\frac{FN+TN}{TP+FN+FP+TN}\times\frac{FP+TN}{TP+FN+FP+TN}$

$P_{o}=\frac{TP+TN}{TP+FN+FP+TN}$

$Overall~{}(OA)=\frac{1}{K}\sum_{i=1}^{K}TP_{i}$ (8)

where $TP$ and $FP$ are true and false positive, $TN$ and $FN$ are true and
false negative, respectively. For fair comparison purposes, the learning rate
for all these models including hybrid models is set to 0.001, Relu as the
activation function for all layers except the output layer on which Softmax is
used, patch size is set of 15, and for all the experiments, 15 most
informative bands have been selected using principal component analysis to
reduce the computational load. The convergence accuracy and loss of our
proposed pipeline and all other competing methods for 50 epochs is presented
in Figure 1. From loss and accuracy curves, one can conclude that the
smoothing has faster convergence than the model without convergence.

(a) Accuracy

(b) Loss

Figure 1: Accuracy and Loss for Training and Validation sets on Indian Pines
for $50$ epochs.

###  III-A Indian Pines

Indian pines (IP) dataset is acquired using AVIRIS sensor over the
northwestern Indiana test site. IP data consist of $145\times 145$ spatial
dimensions and 224 spectral dimensions with total of 16 classes in which all
are not mutually exclusive. Some of the water absorption bands are removed and
the remaining 200 bands are used for the experimental process. This data is
consists of 2/3 agriculture, 1/3 forest, and other vegetation. Less than 5% of
total coverage consists of crops that are in an early stage of growth.
Building, low-density housing, two dual-lane highway, small roads, and a
railway line are also a part of it. Further details about the experimental
datasets can be found at [22]. Table I and Figure 2 presents in-depth
comparative accuracy analysis on the IP dataset.

TABLE I: Indian Pines: Performance analysis of different state-of-the-art
models. The higher accuracies are emphasised.

Class | Train/Val/Test | 2D | 3D | AlexNet | LeNet | WoS | WS
---|---|---|---|---|---|---|---
Alfalfa | 11/12/23 | 100 | 91.3043 | 95.6521 | 82.6086 | 100 | 100
Corn-notill | 357/357/714 | 98.3193 | 93.5574 | 97.3389 | 97.0588 | 98.4593 | 98.8795
Corn-mintill | 207/208/415 | 99.5180 | 66.7469 | 98.3132 | 99.5180 | 99.5180 | 99.5180
Corn | 59/59/118 | 94.0677 | 90.6779 | 93.2203 | 99.1525 | 100 | 100
Grass-pasture | 121/121/242 | 98.3471 | 97.1074 | 96.2809 | 94.2148 | 95.0413 | 96.2809
Grass-trees | 182/183/365 | 98.9041 | 97.5342 | 98.3561 | 98.9041 | 99.4520 | 99.7260
Grass-mowed | 7/7/14 | 92.8571 | 92.8571 | 100 | 100 | 92.8571 | 100
Hay-windrowed | 119/120/239 | 100 | 100 | 100 | 100 | 100 | 100
Oats | 5/5/10 | 70 | 0 | 100 | 70 | 100 | 100
Soybean-notill | 243/243/486 | 98.5596 | 82.3045 | 93.6213 | 99.3827 | 97.9423 | 97.9423
Soybean-mintill | 614/614/1228 | 99.6742 | 92.4267 | 97.8013 | 99.9185 | 99.7557 | 99.8371
Soybean-clean | 148/149/297 | 97.6430 | 98.6531 | 96.9696 | 95.2861 | 98.3164 | 99.6632
Wheat | 51/51/102 | 99.0196 | 98.0392 | 100 | 99.0196 | 99.0196 | 99.0196
Woods | 316/317/633 | 99.8420 | 99.3680 | 99.5260 | 98.8941 | 99.2101 | 99.8420
Buildings | 96/97/193 | 99.4818 | 90.6735 | 100 | 91.1917 | 100 | 99.4818
Stone-Steel | 23/23/46 | 100 | 97.8260 | 100 | 93.4782 | 100 | 100
Training Time | 55.6695 | 250.1662 | 919.5566 | 61.8763 | 250.1622 | 248.5993
Test Time | 1.4897 | 4.0402 | 5.6891 | 1.2752 | 4.6499 | 3.9997
Overall Accuracy | 98.9463 | 91.5707 | 97.6585 | 98.14634 | 98.9853 | 99.2975
Average Accuracy | 98.7980 | 86.8173 | 97.9425 | 94.9142 | 98.7232 | 99.3869
Kappa ($\kappa$) | 96.6396 | 90.3561 | 97.3312 | 97.8853 | 98.8430 | 99.1990

(a) 2D

(b) 3D

(c) A.Net

(d) LeNet

(e) WoS

(f) WS

Figure 2: Indian Pines: Classification accuracy: Fig. 2(a): 2D-CNN = 98.94%,
Fig. 2(b): 3D CNN = 91.57%, Fig. 2(c): AlexNet = 97.65%, Fig. 2(d): LeNet =
98.14%, Fig. 2(e): WoS (Hybrid net without Smoothing) = 98.98%, and Fig. 2(f)
WS (with Smoothing) = 99.29%.

###  III-B Pavia University

Pavia University (PU) dataset acquired using Reflective Optics System Imaging
Spectrometer (ROSIS) optical sensor over Pavia in northern Italy. PU dataset
is distinguished into 9 different classes. PU consists of $610\times 610$
spatial samples per spectral band and 103 spectral bands with a spatial
resolution of 1.3m. Further details about the experimental datasets can be
found at [22]. Table II and Figure 3 presents in-depth comparative accuracy
analysis on the IP dataset.

TABLE II: Pavia University: Performance analysis of different state-of-the-art
models. The higher accuracies are emphasised.

Class | Train/Val/Tesst | 2D | 3D | AlexNet | LeNet | WoS | WS
---|---|---|---|---|---|---|---
Asphalt | 1658/1658/3316 | 100 | 100 | 98.9143 | 100 | 100 | 100
Meadows | 4662/4662/9324 | 100 | 99.9892 | 100 | 100 | 100 | 100
Gravel | 524/525/1049 | 99.5233 | 99.3326 | 95.5195 | 99.4280 | 99.5233 | 99.6186
Trees | 766/766/1532 | 99.6736 | 100 | 98.8250 | 99.7389 | 100 | 100
Painted | 336/337/673 | 100 | 100 | 100 | 100 | 100 | 100
Soil | 1257/1257/2514 | 100 | 100 | 99.9602 | 100 | 100 | 100
Bitumen | 332/333/665 | 100 | 100 | 99.6992 | 100 | 100 | 100
Bricks | 920/921/1841 | 99.8913 | 99.8913 | 97.9359 | 100 | 99.8913 | 99.8913
Shadows | 237/237/74 | 99.3670 | 99.5780 | 98.5232 | 99.7890 | 99.7890 | 100
Training Time | 296.0174 | 1145.9233 | 4716.4900 | 308.6389 | 1143.5753 | 1143.4996
Test Time | 5.7400 | 14.8634 | 24.1701 | 4.5054 | 15.0601 | 15.5904
Overall Accuracy | 99.9298 | 99.9438 | 99.3033 | 99.9485 | 99.9625 | 99.9719
Average Accuracy | 99.8283 | 99.8657 | 98.8197 | 99.8839 | 99.9115 | 99.9455
Kappa ($\kappa$) | 99.9070 | 99.9256 | 99.0768 | 99.9318 | 99.9504 | 99.9628

(a) 2D

(b) 3D

(c) A.Net

(d) LeNet

(e) WoS

(f) WS

Figure 3: Pavia University: Classification accuracy: Fig. 3(a): 2D-CNN =
99.9070%, Fig. 3(b): 3D CNN = 99.9256%, Fig. 3(c): AlexNet = 99.0768%, Fig.
3(d): LeNet = 99.9318%, Fig. 3(e): WoS (Hybrid net without Smoothing) =
99.9504%, and Fig. 3(f): WS (with Smoothing) = 99.9628%.

##  IV Conclusion

The paper proposed a 3D/2D regularized CNN feature hierarchy for HSIC, in
which the loss contribution is considered as a mixture of entropy between a
predicted distribution and the hot-encoding, and the entropy between the
predicted and noise distribution. Several other regularization techniques
(e.g., dropout, L1, L2, etc.) have also been used, however, these techniques
to some extent, lead to predicting the samples extremely confidently which is
not good from a generalization point of view. Therefore, this work proposed
the use of an entropy-based regularization process to improve the
generalization performance using soft labels. These soft labels are the
weighted average of the hard labels and uniform distribution over entire
ground truths. The entropy-based regularization process prevents CNN from
becoming over-confident while learning and predicting thus improves the model
calibration and bean-search. Extensive experiments have confirmed that the
proposed pipeline outperformed several state-of-the-art methods.

## References

  * [1] A. Alcolea, M. E. Paoletti, J. M. Haut, J. Resano, and A. Plaza, “Inference in supervised spectral classifiers for on-board hyperspectral imaging: An overview,” _Remote Sensing_ , vol. 12, no. 3, p. 534, 2020.
  * [2] H. Ayaz, M. Ahmad, A. Sohaib, M. N. Yasir, M. A. Zaidan, M. Ali, M. H. Khan, and Z. Saleem, “Myoglobin-based classification of minced meat using hyperspectral imaging,” _Applied Sciences_ , vol. 10, no. 19, p. 6862, 2020\.
  * [3] H. Ayaz, M. Ahmad, M. Mazzara, and A. Sohaib, “Hyperspectral imaging for minced meat classification using nonlinear deep features,” _Applied Sciences_ , vol. 10, no. 21, p. 6862, 2020.
  * [4] M. H. Khan, S. Zainab, M. Ahmad, A. Sohaib, H. Ayaz, and M. Mazzara, “Hyperspectral imaging for color adulteration detection in red chili,” _Applied Sciences_ , vol. 10, no. 17, p. 5955, 7783.
  * [5] Z. Saleem, M. H. Khan, M. Ahmad, A. Sohaib, H. Ayaz, and M. Mazzara, “Prediction of microbial spoilage and shelf-life of bakery products through hyperspectral imaging,” _IEEE Access_ , vol. 8, pp. 176 986–176 996, 2020.
  * [6] M. Ahmad, S. Shabbir, R. A. Raza, M. Mazzara, S. Distefano, and A. M. Khan, “Hyperspectral image classification: Artifacts of dimension reduction on hybrid cnn,” _arXiv preprint arXiv:2101.10532_ , 2021.
  * [7] J. Wang, R. Huang, S. Guo, L. Li, M. Zhu, S. Yang, and L. Jiao, “Nas-guided lightweight multiscale attention fusion network for hyperspectral image classification,” _IEEE Transactions on Geoscience and Remote Sensing_ , pp. 1–14, 2021.
  * [8] S. Shabbir and M. Ahmad, “Hyperspectral image classification–traditional to deep models: A survey for future prospects,” _arXiv preprint arXiv:2101.06116_ , 2021.
  * [9] S. Jia, S. Jiang, Z. Lin, N. Li, M. Xu, and S. Yu, “A survey: Deep learning for hyperspectral image classification with few labeled samples,” _Neurocomputing_ , vol. 448, pp. 179–204, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0925231221004033
  * [10] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, _Learning Representations by Back-Propagating Errors_. Cambridge, MA, USA: MIT Press, 1988, p. 696–699.
  * [11] H. Sha, M. Al Hasan, and G. Mohler, “Learning network event sequences using long short‐term memory and second‐order statistic loss,” _Stat. Anal. Data Min._ , vol. 14, no. 1, p. 61–73, Jan. 2021. [Online]. Available: https://doi.org/10.1002/sam.11489
  * [12] H. Bi, R. Santos-Rodriguez, and P. Flach, “Polsar image classification via robust low-rank feature extraction and markov random field,” in _IGARSS 2020 - 2020 IEEE International Geoscience and Remote Sensing Symposium_ , 2020, pp. 708–711.
  * [13] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the inception architecture for computer vision,” 06 2016.
  * [14] C. Zhang and M. Han, “Multi-feature hyperspectral image classification with l2,1 norm constrained joint sparse representation,” _International Journal of Remote Sensing_ , vol. 42, no. 12, pp. 4789–4808, 2021. [Online]. Available: https://doi.org/10.1080/01431161.2021.1890854
  * [15] E. Real, A. Aggarwal, Y. Huang, and Q. V. Le, “Regularized evolution for image classifier architecture search,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, no. 01, pp. 4780–4789, Jul. 2019. [Online]. Available: https://ojs.aaai.org/index.php/AAAI/article/view/4405
  * [16] B. Zoph, V. Vasudevan, J. Shlens, and Q. V. Le, “Learning transferable architectures for scalable image recognition,” _CoRR_ , vol. abs/1707.07012, 2017. [Online]. Available: http://arxiv.org/abs/1707.07012
  * [17] T. Carneiro, R. V. M. Da Nóbrega, T. Nepomuceno, G.-B. Bian, V. H. C. De Albuquerque, and P. P. Reboucas Filho, “Performance analysis of google colaboratory as a tool for accelerating deep learning applications,” _IEEE Access_ , vol. 6, pp. 61 677–61 685, 2018.
  * [18] S. Mei, X. Chen, Y. Zhang, J. Li, and A. Plaza, “Accelerating convolutional neural network-based hyperspectral image classification by step activation quantization,” _IEEE Transactions on Geoscience and Remote Sensing_ , pp. 1–12, 2021.
  * [19] Y. Yuan, C. Wang, and Z. Jiang, “Proxy-based deep learning framework for spectral-spatial hyperspectral image classification: Efficient and robust,” _IEEE Transactions on Geoscience and Remote Sensing_ , pp. 1–15, 2021.
  * [20] H.-C. Li, S.-S. Li, W.-S. Hu, J.-H. Feng, W.-W. Sun, and Q. Du, “Recurrent feedback convolutional neural network for hyperspectral image classification,” _IEEE Geoscience and Remote Sensing Letters_ , pp. 1–5, 2021.
  * [21] M. Ahmad, A. M. Khan, M. Mazzara, S. Distefano, M. Ali, and M. S. Sarfraz, “A fast and compact 3-d cnn for hyperspectral image classification,” _IEEE Geoscience and Remote Sensing Letters_ , pp. 1–5, 2020.
  * [22] _Hyperspectral Datasets Description_ , 2021 (accessed 2020-04-01), http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes.


# GARCHX-NoVaS: A Model-free Approach to Incorporate Exogenous Variables

Kejin Wu  Department of Mathematics, University of California San Diego  Sayar
Karmakar  Department of Statistics, University of Florida  Rangan Gupta
Department of Economics, University of Pretoria

###### Abstract

In this work, we explore the forecasting ability of a recently proposed
normalizing and variance-stabilizing (NoVaS) transformation with the possible
inclusion of exogenous variables. From an applied point-of-view, extra
knowledge such as fundamentals- and sentiments-based information could be
beneficial to improve the prediction accuracy of market volatility if they are
incorporated into the forecasting process. In the classical approach, these
models including exogenous variables are typically termed GARCHX-type models.
Being a Model-free prediction method, NoVaS has generally shown more accurate,
stable and robust (to misspecifications) performance than that compared to
classical GARCH-type methods. This motivates us to extend this framework to
the GARCHX forecasting as well. We derive the NoVaS transformation needed to
include exogenous covariates and then construct the corresponding prediction
procedure. We show through extensive simulation studies that bolster our claim
that the NoVaS method outperforms traditional ones, especially for long-term
time aggregated predictions. We also provide an interesting data analysis to
exhibit how our method could possibly shed light on the role of geopolitical
risks in forecasting volatility in national stock market indices for three
different countries in Europe.

Keywords: Volatility forecasting; Model-free prediction; GARCH; GARCHX

JEL Codes: C32; C53; C63; Q54

##  1 Introduction

In the long history of time series econometrics literature, accurate
forecasting has always stood out as a fundamental and important problem. It
has a range of applications in various industries, e.g., weather forecasting,
climate forecasting, and economic forecasting. Discrete-time series data,
e.g., heights of ocean tides, and temperature of a city, is the realization of
a stochastic process $\\{X_{t},t\in\mathbb{Z}\\}$. The earliest modern time
series analysis could be traced back to the work of Yule,  (1927) where the
pattern of the sunspots number was studied. Unlike the prediction of
independent data, the prediction of time series gets more complicated due to
the inherent data dependence. To get accurate predictions and inferences, it
is crucial to model the dependent relationship within the data. Usually, very
generally speaking, the time series data is assumed to be generated by some
underlying mechanism as follows:

$X_{t}=G(\bm{X}_{t-p},\epsilon_{t});$ (1)

$G(\cdot,\cdot)$ could be any suitable function; $\epsilon_{t}$ is called
innovation and assumed to be $i.i.d.$ with appropriate moments and independent
with $X_{t-i}$, $i\geq 1$; $\bm{X}_{t-p}$ represents
$\\{X_{t-1},\ldots,X_{t-p}\\}$ and stands for the historical information. To
further simplify the forecasting problem, participators focus on some standard
formats of $G(\cdot,\cdot)$, e.g., linear or non-linear. For linear models,
such as linear AR, MA and ARMA models, we can apply the Box-Jenkins method of
identifying, fitting, checking and predicting models systematically (Box et
al., , 2015). However, the prediction of non-linear models is not as trivial
as the case of linear models since the innovation must be appropriately
included in the prediction process, especially for the multi-step ahead
predictions; see Wu and Politis,  (2023) for more related discussions.

In this paper, we are exclusively interested in one non-linear type of Eq. 1
which is the so-called Generalized Auto-Regressive Conditional
Heteroskedasticity (GARCH) model proposed by Bollerslev,  (1986) and has a
form below:

$\begin{split}Y_{t}&=\sigma_{t}W_{t}\\\
\sigma_{t}^{2}&=a+a_{1}Y_{t-1}^{2}+b_{1}\sigma_{t-1}^{2};\end{split}$ (2)

where, $a\geq 0$, $a_{1}>0$, $b_{1}>0$, and $W_{t}\sim i.i.d.~{}N(0,1)$. The
GARCH model is a generalization of the famous Autoregressive Conditional
Heteroskedasticity (ARCH) model proposed by Engle,  (1982). Its ability to
forecast the absolute magnitude and quantiles or entire density of squared
financial log-returns (i.e., equivalent to volatility forecasting to some
extent) was shown by Engle and Patton,  (2001) using the Dow Jones Industrial
Index. Later, many studies to investigate the performance of different GARCH-
type models in predicting volatility of financial series were conducted; see
following references (Peters, , 2001; González-Rivera et al., , 2004; Lim and
Sek, , 2013; Herrera et al., , 2018; Karmakar and Roy, , 2021). For
ARCH/GARCH-type models, it is usual practice to identify and fit models based
on quasi-maximum likelihood inference.

Traditionally, economists primarily utilize univariate GARCH-family models to
understand dynamics of econometric data such as stock/ index/ price, etc
observed for a long time. However, one of the key focuses of financial
econometrics is to understand how extra knowledge such as fundamentals- and
sentiments-based information could be beneficial to improve the prediction
accuracy of market volatility if they are incorporated into the forecasting
process; see more discussion from Engle and Patton,  (2007) and the references
therein 111In this regard, there is also a large literature that involves
incorporating information of low-frequency variables using the GARCH-Mixed
Data Sampling (MIDAS) model, as originally developed by Engle et al.,
(2013).. In line with the classical GARCH methods, we can wrap the exogenous
covariates into the prediction process to get the GARCH models augmented with
additional explanatory variables, GARCHX models in short. The estimation
methodology of GARCHX models was discussed thoroughly in the work of Francq et
al.,  (2019); see more details about the GARCHX model in Section 2.1.

Rather than taking the traditional approach (i.e., specifying and fitting a
model and then predicting), we consider a recently proposed model-free
prediction idea, namely normalizing and variance-stabilizing transformation
(NoVaS transformation). The NoVaS method was initially developed by Politis,
(2003) and then well discussed under the framework of the Model-free
prediction principle in Politis,  (2015). In short, the Model-free prediction
principle hinges on the idea of applying an inverse transformation function to
bridge two equivalent probability spaces. For example, if we observe a
univariate time series $\\{Y_{1},\ldots,Y_{T}\\}$, we can try to find a
transformation to map $\\{Y_{1},\ldots,Y_{T}\\}$ to an $i.i.d.$ series
$\\{Z_{1},\ldots,Z_{T}\\}$. Since the prediction of $i.i.d.$ data is trivial,
we can then transform the prediction of $i.i.d.$ data back to the prediction
of the original data; see more details about the Model-free prediction
principle in Section 2.2.

Following the literature, we call this transformation-based approach a Model-
free method in this paper. This necessitates us to emphasize that although the
transformation is inspired from a model-assumption, for the prediction step it
does not tie with any specific model structure. In the huge literature of
applied econometrics, while analyzing data observed over a long time that can
show signs of heteroscedasticity, the usual practice is to pick some specific
GARCH-type model. Next an estimation of that model is carried out accordingly
and subsequently the forecast of future volatility will be made based on the
estimated model. Apparently, there is no universal rule for the choice of the
specific GARCH model. In other words, a specific GARCH model can not work
uniformly well across different datasets compared to other variants. On the
other hand, the Model-free prediction could work well for any scenario as long
as a suitable transformation function can be found. Therefore, the model-
selection stage is shunned and the Model-free prediction approach is more
robust against the model misspecification. Moreover, the standard GARCH
methods require a relatively large sample size to be estimated well. For the
NoVaS method, it tends to work stably even with short data. The existence of
such transformation function in the context of predicting with exogenous
variables will be analyzed theoretically in Section 2.2.

Empirically speaking, NoVaS methods were mainly applied to forecast volatility
in financial econometrics in the past few years. Gulay and Emec,  (2018)
showed that the NoVaS method could beat GARCH-type models (GARCH, EGARCH and
GJR-GARCH) with generalized error distributions by comparing the pseudo-out of
sample (POOS) forecasting of volatility. Here the POOS forecasting analysis
means using data up to and including the current time to predict future
values. Later, Chen and Politis,  (2019) extended the NoVaS method to do
multi-step ahead predictions. Wu and Karmakar,  (2021) further substantiated
the great performance of NoVaS methods on time-aggregated long-term (30-steps
ahead) predictions. Wang and Politis,  (2022) applied the Model-free idea to
provide estimation and prediction inference for a general class of time
series. Our present work is motivated by the Wu and Karmakar,  (2023) work
where the authors recommended a so-called GARCH-NoVaS (GA-NoVaS)
transformation structure inspired by the development of GARCH from ARCH. This
NoVaS method is significantly robust against different model misspecification.
Given this, it was a natural and probably quite an important question to see
if such a robust forecasting framework can be built where exogenous covariates
can be included and thus improve forecasting accuracy.

In this work, we explore the new methodology of forecasting stock market
volatility with additional covariates being available to be included in the
volatility dynamics. As far as we know, the NoVaS model-free prediction idea
has not been studied when the exogenous variables are featured even in
modeling the mean or average let alone the more complicated variance or
volatility dynamics of a time-series. Due to the superior performance of the
GARCH-NoVaS method in volatility forecasting, for this paper, we stick to the
variance part and attempt to further boost the ability of the GA-NoVaS method
with the help of exogenous covariate information. Towards this, we propose a
so-called GARCHX-NoVaS (abbreviated as GAX-NoVaS henceforth) method which
takes the GARCHX model as the starting step to build transformation. To obtain
the inference about the future situation at an overall level, we choose the
time-aggregated prediction metric. This aggregated metric has been applied to
evaluate future predictions of electricity price or financial data (Fryzlewicz
et al., , 2008; Chudý et al., , 2020; Karmakar et al., , 2022); see the formal
definition in Section 4.1. We wanted to check if the NoVaS prediction method
can incorporate exogenous variables. More importantly, we hope the GAX-NoVaS
method can sustain its great performance compared to the GARCHX method.

In addition to comparing our model-free method and classical GARCHX model with
several simulated datasets, we also provide an interesting real data analysis.
Our goal is to exhibit how our method could possibly shed light on the role of
geopolitical risks, which are currently engulfing the global economy with
multiple wars taking place, in forecasting the volatility of three stock
markets of Europe: Germany and its two neighbors (Austria and Switzerland),
based on a daily index of uncertainty associated with the Russia-Ukraine war
as perceived by German Twitter activity. We also use newspaper-based metrics
of global geopolitical risks due to acts and threats, as developed by Caldara
and Iacoviello,  (2022), to check for the robustness of our result covering a
longer data sample. Hence, we add from a methodological perspective to the
existing literature on forecasting international stock returns volatility
using the information contained in geopolitical events and threats that
basically rely on GARCH-type models (see, for example, Salisu et al.,  (2022);
Zhang et al.,  (2023) for details discussion of this literature). In this
regard, note that, Caldara and Iacoviello,  (2022) pointed out that
entrepreneurs, market participants, and central bank officials view
geopolitical risks as key determinants of investment decisions and stock
market dynamics, with such risks, along with economic and policy
uncertainties, forming an “uncertainty trinity” that would adversely impact
the economy and the financial sector, as has been traditionally reported in
the large existing literature on the impact of terror attacks and threats
Mehmet et al.,  (2018); Christos et al.,  (2019); Bouri et al.,  (2020). In
our empirical and simulation exercises, we measure the performance of GARCHX
and GAX-NoVaS methods by the standard mean square prediction error. Moreover,
we apply the forecast comparison tests to compare the two methods in a
statistical way.

Our main contributions are summarized as follows:

  * •

We propose a new methodology– namely GAX-NoVaS–to do the volatility
forecasting with exogenous variables. This model-free method depends on a
transformation function to connect two equivalent probability spaces instead
of relying on any model assumption. The idea behind the GAX-NoVaS method
hinges on the Model-free prediction principle.

  * •

We show such a transformation function exists under some mild conditions. This
serves as the theoretical foundation of our method.

  * •

We apply our new method and standard GARCHX model to investigate the role of
geopolitical risks in forecasting volatility. It turns out that our new method
can be significantly more accurate, especially for long-horizon time
aggregated predictions.

We organize the remainder of this article as follows. In Section 2, we review
the classical forecasting model, namely GARCHX, which is used as the starting
point to propose the GAX-NoVaS method. Also, we present more details of the
Model-free prediction principle and prove the existence of a transformation
function with some exogenous variables existing. In Section 3, we delineate
the details of proposed GAX-NoVaS method. Then, some simulation studies and
model evaluation criteria are collated in Section 4. Next, we contrast our
methods to existing classical ones on three empirical datasets in Section 5.
Finally, in Section 6, we conclude by discussing the implications of our
findings and some future directions.

##  2 GARCHX estimation and Model-free prediction principle

Before introducing our GAX-NoVaS method, we first explain the GARCHX model
since the transformation of GAX-NoVaS is based on the GARCHX model. In
addition, we give more details on the Model-free prediction principle and we
prove the existence of a transformation function to achieve the Model-free
prediction goal.

###  2.1 GARCHX model

In a seminal work, ARCH was proposed by Engle,  (1982) to model volatility or
$\sigma_{t}^{2}$ for a time-series in a dynamic way. Following this, many
different variants were developed in the econometrics literature. The GARCH
model, especially the GARCH(1,1), stands as possibly the most popular one. The
classic GARCH(1,1) model as defined by Bollerslev,  (1986) can be described
below:

$\begin{split}Y_{t}&=\sigma_{t}W_{t},\\\
\sigma_{t}^{2}&=a+a_{1}Y_{t-1}^{2}+b_{1}\sigma_{t-1}^{2};\end{split}$ (3)

where $a\geq 0$, $a_{1}>0$, $b_{1}>0$, and $W_{t}\sim i.i.d.~{}N(0,1)$. More
generally, we can express the GARCH-type models as

$Y_{t}=\sigma_{t}\eta_{t};$

where a particular distribution of $\eta_{t}$ is not necessary and we can only
assume that $\mathbb{E}(\eta_{t}^{2}|\mathscr{F}_{t-1})=1$. Usually,
$\mathscr{F}_{t-1}$ is taken as the sigma-field generated by previous
information $\\{Y_{t-1},\ldots\\}$. When additional information is available,
people would like to utilize this extra knowledge to improve prediction
accuracy. Subsequently, the so-called GARCHX model enters the public eye; see
Francq et al.,  (2019) for discussions on the quasi-maximum likelihood
estimation inference of GARCHX models. To simplify the analysis, we still
assume the normality of $W_{t}$. After taking a vector of exogenous covariates
$\bm{X}=(X_{1},\ldots,X_{m})$ into account, we can wrap the exogenous
covariates into the prediction process by turning the GARCH(1,1) model into
the following GARCHX(1,1,1) model:

$\begin{split}Y_{t}&=\sigma_{t}W_{t},\\\
\sigma_{t}^{2}&=a+a_{1}Y_{t-1}^{2}+b_{1}\sigma_{t-1}^{2}+\bm{c}^{T}\bm{X}_{t-1};\end{split}$
(4)

where $\bm{X}_{t-1}$ represents $(X_{1,t-1},\ldots,X_{m,t-1})$ and $\bm{c}$
are the coefficients of these exogenous variables to be estimated. To perform
a moving-window out-of-sample prediction experiment, we first need to estimate
the GARCH(1,1) and GARCHX(1,1,1) models222For estimation of the GARCH and
GARCHX models, we use the fGarch (Wuertz et al., , 2013) and garchx packages
(Sucarrat, , 2020) in the R language and environment (R Core Team, , 2023).,
and then we compute predictions iteratively; see Section 4 for details. In
this process, we assume that we know the true exogenous variables, which is
feasible because we generate out-of-sample predictions. For practical
applications, if needed, the future exogenous information can be estimated
separately.

###  2.2 Model-free prediction principle

The model-free prediction principle was initially well developed by Politis,
(2015). Later, Chen and Politis,  (2019) applied this idea to multi-step ahead
predictions of financial returns in the context of an ARCH-model structure. In
short, the main idea behind the model-free prediction is to apply an
invertible transformation function, $H_{T}$, that can map a non-$i.i.d.$
vector, $\\{Y_{t}~{};t=1,\ldots,T\\}$, to a vector,
$\\{\epsilon_{t};~{}t=1,\ldots,T\\}$, with $i.i.d.$ components (chosen as
standard normal in this work). Due to the invertibility of the function,
$H_{T}$, it is possible to construct a one-to-one relationship between a
future value, $Y_{T+1}$, and $\epsilon_{T+1}$, i.e.,

$Y_{T+1}=f_{T+1}(\bm{Y}_{T},\bm{X}_{T+1},\epsilon_{T+1});$ (5)

where $\bm{Y}_{T}$ denotes all historical data $\\{Y_{t};~{}t=1,\ldots,T\\}$;
$\bm{X}_{T+1}$ is the collection of all predictors, and it also contains the
value of a future predictor $X_{T+1}$; the form of $f_{T+1}(\cdot)$ depends on
$H^{-1}_{T}$. This relationship implies that we can also transform the
prediction of $\epsilon_{T+1}$ to the prediction of $Y_{T+1}$. Assume we have
$\hat{\epsilon}_{T+1}$ to the be the predictor of $\epsilon_{T+1}$, we can
express the predictor of $Y_{T+1}$ as

$\widehat{Y}_{T+1}=f_{T+1}(\bm{Y}_{T},\bm{X}_{T+1},\hat{\epsilon}_{T+1}).$ (6)

Because the prediction of $i.i.d.$ data is standard, the $L_{1}$ (Mean
Absolute Deviation) , $L_{2}$ (Mean Squared Error), or another optimal
quantile predictor of $\epsilon_{T+1}$ can easily be found. We, thus, can
easily obtain the corresponding optimal predictor of $Y_{T+1}$.

For multi-step ($h$-step) ahead prediction, we simply repeat this prediction
process, i.e., we express $Y_{T+h}$ through a function w.r.t. $\bm{Y}_{T}$,
$\bm{X}_{T+1}$ and $\\{\epsilon_{T+1},\ldots,\epsilon_{T+h}\\}$:

$Y_{T+h}=f_{T+h}(\bm{Y}_{T},\bm{X}_{T+1},\epsilon_{T+1},\ldots,\epsilon_{T+h}).$
(7)

In order to compute the prediction of $Y_{T+h}$, we take a distribution-match
approach to approximate the distribution of $Y_{T+h}$. Ideally, when we know
the exact distribution of the $i.i.d.$ $\epsilon$, we can use a Monte Carlo
simulation to approximate the distribution of $Y_{T+h}$ based on Eq. 7.
Practically speaking, when we just have the empirical transformation results,
i.e., the observed sample $\\{\epsilon_{t}\\}_{t=1}^{T}$, bootstrap is an
appropriate approach. Moreover, we can even predict $g(Y_{T+h})$, where
$g(\cdot)$ is a general continuous function. For example, we can compute the
$L_{1}$ and $L_{2}$ optimal predictors of $g(Y_{T+h})$ as below:

$\begin{split}g(Y_{T+h})_{L_{2}}&=\frac{1}{M}\sum_{m=1}^{M}g(f_{T+h}(\bm{Y}_{T},\bm{X}_{T+1},\hat{\epsilon}_{T+1,m},\ldots,\hat{\epsilon}_{T+h,m})),\\\
g(Y_{T+h})_{L_{1}}&=\text{Median~{}of~{}}\\{g(f_{T+h}(\bm{Y}_{T},\bm{X}_{T+1},\hat{\epsilon}_{T+1,m},\ldots,\hat{\epsilon}_{T+h,m}));m=1,\ldots,M\\};\end{split}$
(8)

where $g(Y_{T+h})_{L_{2}}$ and $g(Y_{T+h})_{L_{1}}$ represent the optimal
$L_{2}$ and $L_{1}$ predictor of $g(Y_{T+1})$, the
$\\{\hat{\epsilon}_{T+1,m}\\}_{m=1}^{M}$ are generated by bootstrap or Monte
Carlo simulation, and $M$ is some large number (2000 in our empirical
analysis). For further discussion, see Politis,  (2015).

To the best of our knowledge, the model-free prediction idea has not been
studied when the model features exogenous variables. However, it may be
beneficial to take into account in the prediction process the additional
information embedded in such exogenous variables. To show the NoVaS approach
is still applicable, we need a transformation function that maps the targeted
variables and exogenous predictors together into some simple $i.i.d.$ random
variables. Under some mild conditions, we show the existence of such a
transformation function based on the probability integral transform. We
assume:

  * •

A1 The joint density of $\\{Y_{1},\cdots,Y_{T}\\}$ exists for any $T\geq 1$.

  * •

A2 For exogenous random vector $\bm{X}:=\\{X_{1},\cdots,X_{m}\\}$, the joint
density $\\{Y_{1},\cdots,Y_{T},X_{1},\cdots,X_{m}\\}$ exists for any $m\geq
1$.

Then, the feasibility of NoVaS transformation with exogenous variables
existing is guaranteed by 2.1 shown below:

######  Theorem 2.1.

Under A1 and A2, there exists a function $\bm{g}$ such that
$\bm{Z}=\bm{g}((\bm{Y},\bm{X}))$ and the corresponding inverse function
$\bm{h}$ such that $(\widetilde{\bm{Y}},\widetilde{\bm{X}})=\bm{h}(\bm{Z})$;
$\bm{Z}\sim N(0,\bm{I}_{T+m})$; $\bm{Y}=(Y_{1},\cdots,Y_{T})$ and
$\bm{X}=(X_{1},\cdots,X_{m})$ are any two random vectors;
$(\widetilde{\bm{Y}},\widetilde{\bm{X}})$ have the same joint distribution of
$(\bm{Y},\bm{X})$.

###### Proof.

The proof of 2.1 is based on the probability integral transform; see Angus,
(1994) for a review. Without loss of generality, we start from $Y_{1}$ to
determine the transformation function $\bm{g}$. Let
$U_{1}:=\tilde{g}_{1}(Y)=F(Y_{1})$; $F(Y_{1})$ is the distribution of $Y_{1}$.
According to the probability integral transform, we know $U_{1}$ has a uniform
distribution on $[0,1]$. Then, we make

$U_{2}:=\tilde{g}_{2}(Y_{1},Y_{2})=F(Y_{2}|Y_{1}).$ (9)

$F(Y_{2}|Y_{1})$ is the conditional distribution of $Y_{2}$. Eq. 9 implies
that $Z_{2}$ is $\text{Uniform}(0,1)$ conditional on $Y_{1}$. Thus, the
unconditional (marginal) distribution of $U_{2}$ is still
$\text{Uniform}(0,1)$, and $U_{2}$ and $Y_{1}$ are independent so that $U_{2}$
is also independent with $U_{1}$. This can be seen from the equation below:

$p_{U_{2},Y_{1}}(u_{2},y_{1})=p_{U_{2}|Y_{1}}(u_{2}|y_{1})p_{Y_{1}}(y_{1}).$
(10)

Integrating both sides w.r.t. $y_{1}$, we can find $p_{U_{2}}(u_{2})=1$ on the
region $[0,1]$, since $Z_{2}$ is $\text{Uniform}(0,1)$ conditional on
$Y_{1}=y_{1}$ for any $y_{1}$. We can repeat this process as a Gram–Schmidt-
like recursion, i.e., we let
$U_{3}:=\tilde{g}_{3}(Y_{1},Y_{2},Y_{3})=F(Y_{2}|Y_{1},Y_{2})$ and so on. In
total, we need $\\{\tilde{g}_{1},\cdots,\tilde{g}_{T+m}\\}$ and they are
functions of $\bm{Y}$. Thus, there exists a function $\tilde{\bm{g}}$ which
maps $(\bm{Y},\bm{X})$ to $\bm{U}$ which has $i.i.d.$ uniform components
$\\{U_{1},\cdots,U_{m+T}\\}$, i.e., $\bm{U}=\tilde{\bm{g}}((\bm{Y},\bm{X}))$.
Then, $\bm{Z}=\Phi^{-1}\circ\tilde{\bm{g}}((\bm{Y},\bm{X}))$ has multivariate
normal distribution $N(0,\bm{I}_{T+m})$; $\Phi^{-1}$ is the quantile function
of $N(0,\bm{I}_{T+m})$. Finally, we can take
$\bm{g}=\Phi^{-1}\circ\tilde{\bm{g}}$.

On the other hand, if $F^{-1}_{Y_{1}}:=\inf\\{x:F_{Y_{1}}(x)\geq y\\},0\leq
y\leq 1$, then $F_{Y_{1}}^{-1}(U_{1})$ has the distribution as the same as
$Y_{1}$. Similarly, we can get the conditional distribution of $Y_{2}$ on
$Y_{1}$ by taking $F_{Y_{2}|Y_{1}}^{-1}(U_{2})$. By repeating this process, we
can recover the joint distribution of $(\bm{Y},\bm{X})$ by chain rule. In
other words, there exists a $\bm{h}$ such that
$(\bm{Y},\bm{X})=\bm{h}(\bm{Z})$. ∎

The direct implication of 2.1 is that the Model-free prediction principle is
feasible even when exogenous variables are included in the dependence
dynamics. Moreover, our theorem is more general than the result from Wang and
Politis,  (2022) where the time series must satisfy some strict conditions. In
short, we build a transformation function based on the GARCHX model structure
to estimate the oracle functions $\bm{g}$ and $\bm{h}$, so we call our method
GAX-NoVaS. We should mention again that the “model-free” in this context means
we do not rely on an assumed underlying model to make predictions. Although a
transformation function needs to be estimated, it is just a “bridge” that
connects original and transformed distributions according to the distribution
match idea.

##  3 GAX-NoVaS model-free prediction method

We first present the state-of-the-art GARCH-NoVaS method which is based on a
so-called NoVaS transformation. Then, we extend the GARCH-NoVaS method to a
GAX-NoVaS method, which features the exogenous variables.

###  3.1 NoVaS transformation

For the sake of completeness, we first give a brief introduction to the NoVaS
transformation (model) which is a direct application of the Model-free
prediction idea explained in Section 2.2. Initially, the NoVaS transformation
is developed from the ARCH model:

$Y_{t}=W_{t}\sqrt{a+\sum_{i=1}^{p}a_{i}Y_{t-i}^{2}}.$ (11)

here, these parameters satisfy $a\geq 0$, $a_{i}\geq 0$, for all
$i=1,\ldots,p$; $W_{t}\sim i.i.d.~{}N(0,1)$. In other words, the structure of
the ARCH model gives us a ready-made $H^{-1}_{T}$. We can express $W_{t}$ in
Eq. 11 using other terms to get a potential $H_{T}$ :

$W_{t}=\frac{Y_{t}}{\sqrt{a+\sum_{i=1}^{p}a_{i}Y_{t-i}^{2}}}~{};~{}\text{for}~{}t=p+1,\ldots,T.$
(12)

Politis,  (2003) further modified Eq. 12 as follows:

$W_{t}=\frac{Y_{t}}{\sqrt{\alpha s_{t-1}^{2}+\beta
Y_{t}^{2}+\sum_{i=1}^{p}a_{i}Y_{t-i}^{2}}}~{};~{}\text{for}~{}t=p+1,\ldots,T.$
(13)

here, $\\{Y_{t};~{}t=1,\ldots,T\\}$ is the sample data;
$\\{W_{t};~{}t=p+1,\ldots,T\\}$ is the transformed vector; $\alpha$ is a fixed
scale invariant constant; $s_{t-1}^{2}$ is an estimator of the variance of
$\\{Y_{i};~{}i=1,\ldots,t-1\\}$ and can be calculated by
$(t-1)^{-1}\sum_{i=1}^{t-1}(Y_{i}-\overline{Y})^{2}$, where $\overline{Y}$ is
the sample mean of $\\{Y_{i};~{}i=1,\ldots,t-1\\}$. For making Eq. 13 be a
qualified function $H_{T}$, i.e., making $\\{W_{t}\\}_{t=p+1}^{T}$ really obey
$i.i.d.$ standard normal distribution, we need to impose some restrictions on
$\alpha$ and $\beta,a_{1},\ldots,a_{p}$. We first stabilize the variance by
requiring:

$\alpha\geq 0,\beta\geq 0,a_{i}\geq 0~{};~{}\text{for all}~{}i\geq
1,\alpha+\beta+\sum_{i=1}^{p}a_{i}=1.$ (14)

In application, $\\{W_{t}\\}_{t=p+1}^{T}$ transformed from financial log-
returns by NoVaS transformation are usually uncorrelated. Therefore, if we
make $\\{W_{t}\\}_{t=p+1}^{T}$ close to a Gaussian series i.e., normalizing
$\\{W_{t}\\}_{t=p+1}^{T}$, we can get the desired $i.i.d.$ property. This is
why this transformation is called NoVaS.

There are many criteria to measure the normality of a series. Under the
observation that the distribution of financial log-returns is usually
symmetric, we choose the kurtosis to be a simple distance to measure the
departure of a non-skewed dataset from that of the standard normal
distribution (Politis, , 2015). Besides, matching marginal distribution seems
sufficient to normalize the joint distribution of $\\{W_{t}\\}_{t=p+1}^{T}$
for practical purposes based on empirical results. If we denote the marginal
distribution of $\\{W_{t}\\}_{t=p+1}^{T}$ and the corresponding kurtosis by
$\widehat{F}_{w}$ and $\text{KURT}(W_{t})$, respectively, we then attempt to
minimize $|\text{KURT}(W_{t})-3|$ to obtain the optimal combination of
$\alpha,\beta,a_{1},\ldots,a_{p}$ such that $\widehat{F}_{w}$ is as close to
standard normal distribution as possible. Subsequently, the NoVaS
transformation can be determined.

The remaining difficulty is how to finish this optimization step to get
optimal coefficients $\alpha,\beta,a_{1},\ldots,a_{p}$, especially when $p$ is
large. To simplify this problem, Politis,  (2015) defined an exponentially
decayed form of $\\{a_{i}\\}_{i=1}^{p}$:

$\alpha\neq 0,\beta=c^{\prime},a_{i}=c^{\prime}e^{-ci}~{};~{}\text{for
all}~{}1\leq i\leq p,c^{\prime}=\frac{1-\alpha}{\sum_{j=0}^{p}e^{-cj}}.$ (15)

The NoVaS transformation based on coefficients defined in Eq. 15 is called
Generalized Exponential NoVaS (GE-NoVaS). In other words, we can represent the
$p+2$ number of coefficients by two parameters $c$ and $\alpha$, which relieve
the optimization burden, but with a sacrifice that the coefficients are fixed
in a decayed form. To achieve a balance between the relief of the optimization
dilemma and the freedom of coefficients, inspired by the development of GARCH
from ARCH, Wu and Karmakar,  (2023) built a NoVaS transformation according to
the GARCH model, namely GARCH-NoVaS which was shown to be more stable and
accurate. Later, we specify the GAX-NoVaS transformation in detail.

###  3.2 GAX-NoVaS transformation method

Starting from Eq. 4, we take similar steps of building GA-NoVaS to find the
transformation function of the GAX-NoVaS method. In order to simplify the
notation, we consider the case of only one exogenous covariate $X_{t}$. The
case of multiple exogenous covariates can be analyzed in an analogous way.
First, we notice that we can rewrite the Eq. 4 as

$W_{t}=\frac{Y_{t}}{\sqrt{a_{0}+a_{1}Y^{2}_{t-1}+b_{1}\sigma_{t-1}^{2}+c_{1}X_{t-1}}}.$
(16)

We also have

$\begin{split}&\sigma_{t-1}^{2}=a_{0}+a_{1}Y^{2}_{t-2}+b_{1}\sigma_{t-2}^{2}+c_{1}X_{t-2},\\\
&\sigma_{t-2}^{2}=a_{0}+a_{1}Y^{2}_{t-3}+b_{1}\sigma_{t-3}^{2}+c_{1}X_{t-3},\\\
\vdots\end{split}$ (17)

so that we can substitute these terms into Eq. 16. We then get

$W_{t}=\frac{Y_{t}}{\sqrt{\frac{a_{0}}{1-b_{1}}+\sum_{i=1}^{p}a_{1}b_{1}^{i-1}Y^{2}_{t-i}+\sum_{i=1}^{p}c_{1}b_{1}^{i-1}X_{t-i}}},$
(18)

where $p$ is a large constant that is used to truncate the infinite summation
(we take $p=q$ in this work), because $a_{1}$ and $b_{1}$ need to be less than
one to guarantee the stationary property of the GARCHX series. In line with
the NoVaS transformation, we finally write the transformation function as
follows:

$W_{t}=\frac{Y_{t}}{\sqrt{\alpha s_{t-1,Y}+\beta
s_{t-1,X}+\sum_{i=1}^{p}a_{1}b_{1}^{i-1}Y^{2}_{t-i}+\sum_{i=1}^{p}c_{1}b_{1}^{i-1}Y_{t-i}}},$
(19)

where $s^{2}_{t-1,Y}$ and $s^{2}_{t-1,X}$ are the sample variance of
$\\{Y_{1},\ldots,Y_{t-1}\\}$ and $\\{X_{1},\ldots,X_{t-1}\\}$, respectively.
Thus, we can use Eq. 19 as the transformation function for the GARCHX model,
where $\\{W_{t}\\}$ is the transformed series. We want to make $\\{W_{t}\\}$
$i.i.d.$ normal, so that the one-step (conditional) prediction
$\widehat{Y}_{T+1}$ can be expressed as

$\widehat{Y}_{T+1}=\widehat{W}_{T+1}\sqrt{\alpha s^{2}_{T-1,Y}+\beta
s^{2}_{T-1,X}+\sum_{i=1}^{p}a_{1}b_{1}^{i-1}Y^{2}_{T-i}+\sum_{i=1}^{p}c_{1}b_{1}^{i-1}X_{T-i}},$
(20)

where $\widehat{W}_{T+1}$ is the optimal point prediction based on $i.i.d.$
$\\{W_{1},\ldots,W_{T}\\}$. Multi-step-ahead predictions can be computed as
explained in Section 2.2.

The final question left now is how to find a transformation function that
indeed makes $\\{W_{t}\\}$ $i.i.d.$ normal. While we have made some brief
remarks on this question in Section 3.1, we next provide a full explanation
with a focus on the GAX-NoVaS method. Our goal is to determine the
coefficients, $\alpha,\beta,a_{1},b_{1},c_{1}$, of Eq. 19 to obtain the
desired transformation. The most important step is to minimize
$|\text{KURT}(W_{t})-3|$, where $3$ is the kurtosis of normal distribution.
For this optimization, we use the numerical technique to find the optimal
coefficients.333We use the nloptr package (Ypma et al., , 2014) for R. In
operation, one may obtain some extremely large values from the transformed
series $\\{W_{t}\\}$, and such outliers may spoil the normality of the
transformed series and may also influence prediction performance. Thus, before
moving on to the prediction step, we truncate the transformed series by the
0.99 and 0.01 quantile values of a normal distribution with mean and standard
deviation given by the sample mean and sample standard deviation of
$\\{W_{t}\\}$.

######  Remark 3.1.

It is not difficult to perceive that the transformed series $\\{W_{t}\\}$ may
be correlated, such that the decorrelation step is beneficial and necessary.
One way to carry out this step is by fitting an AR(p) model on the
$\\{W_{t}\\}$ series. Then, we record the residuals of the AR fit as
$\\{\hat{\epsilon}_{t}\\}$. Also, we can approximate the one-step-ahead value
$\widehat{W}_{T+1}$ with a fitted AR model. Then, we can create the new series
as $\\{\epsilon_{t}+\widehat{W}_{T+1}\\}$. We use the empirical distribution
of this new series to approximate the distribution of $W_{T+1}$. Eq. 20 can be
used with the optimal prediction $\widehat{W}_{T+1}$ derived from this
empirical distribution. It is still an open question, however, how to extend
this decorrelation step to multi-step-ahead predictions.

##  4 Simulations

In this section, we deploy several simulations to check the performance of
GARCH, GARCHX and GAX-NoVaS methods. Before presenting the data-generating
model used to do simulations, we explain the procedure of the moving-window
time-aggregated predictions and give the model evaluation metrics to measure
the performance of different methods.

###  4.1 Moving-window time-aggregated prediction

If we have sample $\\{Y_{1},\ldots,Y_{N}\\}$ at hand, in order to fully
exhaust the dataset, we can focus on moving-window out-of-sample predictions,
i.e., we use $\\{Y_{1},\cdots,Y_{T}\\}$ to predict
{$Y_{T+1}^{2},\cdots,Y_{T+h}^{2}\\}$, then we use $\\{Y_{2},\cdots,Y_{T+1}\\}$
to predict $\\{Y_{T+2}^{2},\cdots,Y_{T+h+1}^{2}\\}$, and so on until we reach
the end of the sample (that is, until we use
$\\{Y_{N-T+h+1},\cdots,Y_{N-h}\\}$ to predict
$\\{Y^{2}_{N-h+1},\ldots,Y^{2}_{N}\\}$). Here, $T$ denotes the moving window
size; we fix its size as 250; $h$ is the prediction horizon, i.e., 1, 5, 20 in
our setting. Sometimes, we may not have enough data available to perform
predictions. Thus, the window size $T=250$ is designed to see if this method
is stable and can still return accurate predictions even with short data. A
250-size moving window is in line with around one year of daily financial
data. In this perspective, it is important to keep in mind that a 250-size
moving window is practically meaningful since the time series may not be
stationary on a wider span.

In addition, we are not only interested in the one-step-ahead prediction $h=1$
but also multi-step-ahead prediction $h>1$. From a practical aspect of
forecasting volatility, as mentioned above in our introduction Section 1, the
long-term prediction ($h$ takes a large value) is important and can guide
future strategic decisions. Before this evaluation, we start by writing time-
aggregated predictions as follows:

$\overline{\widehat{Y}}_{T,h}^{2}=\sum_{k=1}^{h}\widehat{Y}_{T+k}^{2}/h,$ (21)

where $\overline{\widehat{Y}}_{T,h}^{2}$ is the $h$-step ahead time-aggregated
volatility prediction starting from $Y_{T}$. For example, if the total number
of data $N=1000$ and we consider the 6-step-ahead moving window time-
aggregated predictions with $T=500,$ we need to find predictions
$\overline{\widehat{Y}}_{T,6}^{2}$ for $T=500,\ldots,994$.

We hope the time-aggregated volatility prediction is close to the true
aggregated value calculated from the realized average squared log-returns
$\overline{Y}_{l,h}^{2}=\sum_{k=1}^{h}(Y_{T+k}^{2}/h)$. To evaluate the
accuracy, we can consider the specific mean of squared prediction errors
(MSPE) shown below, with this statistic aiming to compare the prediction
performance in an absolute way:

$P=\frac{1}{N-h-T+1}\sum_{l=T}^{N-h}(\overline{\widehat{Y}}_{l,h}^{2}-\overline{Y}_{l,h}^{2})^{2},$
(22)

where $\overline{\widehat{Y}}_{l,h}^{2}$ and $\overline{Y}_{l,h}^{2}$ denote
the predicted and true time-aggregated values for each moving-window
forecasting, respectively.

###  4.2 Simulation setting

In this part, we present three data-generating models to simulate data and
then evaluate the performance of various methods considered in this paper. Due
to the true data-generating process being known to us, we can simulate any
size of the sample and compare predictions from various methods with oracle
values. Besides, we remove the variance term $\beta s^{2}_{T-1,X}$ of Eq. 19
when we do the transformations since the prediction performance hardly changes
with or without $\beta s^{2}_{T-1,X}$. Three true underlying models are
presented below:

  * Model 1: Standard GARCH(1,1) with Student-$t$ errors
$Y_{t}=\sigma_{t}\epsilon_{t},$
$~{}\sigma_{t}^{2}=0.00001+0.73\sigma_{t-1}^{2}+0.1Y_{t-1}^{2}+c|X_{t-1}|,$
$X_{t-1}\sim i.i.d.~{}N(0,1)$;$~{}\\{\epsilon_{t}\\}\sim i.i.d.~{}t$
distribution with four degrees of freedom
$c=1$.

  * Model 2: Standard GARCH(1,1) with Student-$t$ errors
$Y_{t}=\sigma_{t}\epsilon_{t},$
$~{}\sigma_{t}^{2}=0.00001+0.8895\sigma_{t-1}^{2}+0.1Y_{t-1}^{2}+c|X_{t-1}|,$
$X_{t-1}\sim i.i.d.~{}N(0,1)$;$~{}\\{\epsilon_{t}\\}\sim i.i.d.~{}t$
distribution with four degrees of freedom
$c=1$.

  * Model 3: Time-varying GARCHX(1,1) with standard normal exogenous variables:
$Y_{t}=\sigma_{t}\epsilon_{t},~{}\sigma_{t}^{2}=b_{t}\sigma_{t-1}^{2}+a_{t}Y_{t-1}^{2}+c|X_{t-1}|,$
$X_{t-1}\sim i.i.d.~{}N(0,1)$; $\\{\epsilon_{t}\\}\sim i.i.d.$ $t$
distribution with five degrees freedom;
$c=1;g_{t}=t/n$; $a_{t}=0.1-0.05g_{t}$; $b_{t}=0.7+0.2g_{t},~{}n$ is the total
length of the time series.

To check the robustness of methods on model misspecification, we intend to
take the innovation distribution of simulation models as the $t$-distribution.
Besides this purpose, we argue that the $t$ distribution as the innovation to
mimic the real-world cases is more appropriate since real data usually show
the heavy tail phenomenon. Models 1 and 2 are from a standard GARCH where in
Model 2 we intended to explore a scenario that $\alpha_{1}+\beta_{1}$ is very
close to 1 and thus mimics what would happen for the iGARCH situation. In
addition, we make the coefficients of the GARCHX model change linearly in
Mode-3 so that we can observe the ability of different methods to handle the
data generated from a time-varying model which is more coherent to the real-
world situation. To sync with the empirical studies later, we simulate a time
series with a length $T=4694$. We take the moving-window size $T=250$. MSPE of
GARCH, GARCHX and GAX-NoVaS methods with three simulation settings are
presented below:

Table 1: MSPE ratios of different methods on three simulated datasets.

| Model-1 | Model-2 | Model-3
---|---|---|---
Prediction steps | 1 | 5 | 20 | 1 | 5 | 20 | 1 | 5 | 20
GARCH | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
GARCHX | 0.985 | 1.007 | 0.854 | 1.014 | 0.957 | 0.701 | 0.921 | 0.956 | 2.680
GAX-NoVaS | 0.842 | 0.566 | 0.222 | 0.798 | 0.818 | 0.085 | 0.870 | 1.276 | 0.006

Note: To simplify the presentation, we compute the ratio of MSPE of different
methods, i.e., we divide all MSPE for each prediction horizon and models by
the MSPE of the GARCH (benchmark) method. In addition, $0.000$ value in the
above table is a rounding number.

From Table 1, it is clear that the GAX-NoVaS is much better than GARCH and
even the classical GARCHX models according to the MSPE criterion.
Interestingly, the GARCHX model is even much worse than the GARCH model for
some specific cases, e.g., the 20-step-ahead prediction of Model-3. By taking
a deeper analysis, we find the terrible performance GARCHX method is due to
some extremely large predictions. On the other hand, the prediction returned
by GAX-NoVaS is more stable. The superiority is further shown in Section 5
with three real datasets and various exogenous predictors.

Since the focus of this paper is exploring a new approach to incorporate
exogenous variables, we take the DM test to evaluate the performance of GARCHX
and GAX-NoVaS methods more formally; see Diebold and Mariano,  (2002) for the
technical details of the DM-test444We perform the DM-test with the function
DM-test in the R package multDM.. The DM-test results on comparing GARCHX and
GAX-NoVaS for forecasting three simulated datasets are tabularized in Table 5.
These tests further verify the advantage of our methods on forecasting with
exogenous variables, especially for a long-prediction horizon.

Table 2: DM-test results on simulated datasets.

| Model-1 | Model-2 | Model-3
---|---|---|---
Prediction steps | 1 | 5 | 20 | 1 | 5 | 20 | 1 | 5 | 20
GARCHX |  |  |  |  |  |  |  |  |
GAX-NoVaS | 0.009 | 0.020 | 0.000 | 0.038 | 0.107 | 0.000 | 0.172 | 1.000 | 0.102

Note: The values in the different rows are one-sided $p$-values of the DM-test
on the prediction error of GAX-NoVaS and GARCHX. For the DM-test here, the
alternative hypothesis is that the former method GAX-NoVaS is more accurate
than the latter one, GARCHX.

##  5 Empirical analyses with real data

A summarizing note of our findings in Section 4 reads that the GAX-NoVaS
method performs better than standard GARCH-type methods, especially for long-
term time aggregated predictions. In this section, we deploy an interesting
data analysis to exhibit how our method could shed light on the role of
geopolitical risks in forecasting volatility with real-world data. We start by
describing the data below.

###  5.1 Data description

The ongoing Ukraine-Russia has led many countries in Europe, particularly
Germany, to adjust their military and security, as well as energy policies in
light of new geopolitical risks, with such adjustments entailing large costs
to the macroeconomy and financial markets, as depicted by Grebe et al.,
(2024). In this regard, these authors, first, assemble a data set of more than
eight million German Twitter posts related to the war in Ukraine to construct
a daily index of uncertainty about the war as perceived by German Twitter
based on using state-of-the-art methods of textual analysis. Grebe et al.,
(2024) show that an increase in uncertainty has strong effects on financial
markets, associated with a significant decline in economic activity as well as
an increase in expected inflation. We utilize this index (Ukraine)555The data
is available for download from: https://www.uni-
giessen.de/de/fbz/fb02/fb/professuren/vwl/tillmann/forschung/ukraine-
uncertainty-index. in our empirical analysis to forecast stock market
volatility of not only Germany, but two of its neighbors namely, Austria and
Switzerland, over the daily period of 1st January, 2021 to 28th February,
2023. The national stock market indexes (ATX (Austria), DAX (Germany), SMI
(Switzerland)) of these three countries, for which we compute log-returns to
feed into our volatility models were derived from the Bloomberg terminal. With
the focus being on geopolitical risks, we also utilized the daily newspapers-
based geopolitical risks index (GPRD) of Caldara and Iacoviello,  (2022)666The
data can be accessed from: https://www.matteoiacoviello.com/gpr.htm., which,
in turn, allowed us to analyze a longer data sample covering 2nd January, 2006
to 10th August, 2023. The starting date of this longer sample, and the choice
of these three countries, were also motivated by the availability of Google
searches-based daily data on economic activity (Trend) for all three
countries, and inflation (Inflation) for Germany and Switzerland,777The data
can be downloaded from: https://www.trendecon.org/. which are used as
additional predictors to ensure that our results are not only limited to
geopolitical risks.

We present three log-return series of Germany, Switzerland and Austria from
2nd January, 2006 to 10th August, 2023 in Figs. 1, 2 and 3. The volatility
clustering phenomena observed in all plots reveals the heteroskedasticity
within these three series. To investigate the property of three long return
series, we provide the summary statistics, e.g., the mean, skewness, and
kurtosis in Table 3. To verify the heteroskedasticity with all series more
directly, we split the whole time period into four equal-length sub-periods
and denote the sample variance of all four sub-periods by $V_{i}$,
$i=1,\ldots,4$. These statistics are also provided in Table 3. Towards
statistical tests, we also perform modified Ljung-Box (m-LB) and ARCH Lagrange
Multiplier (ALM) tests to check the autocorrelation and ARCH effects of
squared return series. For the m-LB test, we consider the lag order 20. For
the ALM test, we consider the maximum lag order 10. These two tests are
performed in R with functions lbtest and Lm.test, respectively. The p-values
of tests are presented. Summarizing Table 3, the large kurtosis values
indicate the heavy-tailed property for all three log-return series. The
variance of return series in different time regions changes notably,
indicating the heteroskedasticity. The ALM test with a pretty small p-value
also confirms the heteroskedasticity for all return series. The m-LB test
shows strong evidence of autocorrelation within all squared return series.

Table 3: Summary statistics of three long return series.

Returns series | Mean | Skew. | Kurt. | $V_{1}$ | $V_{2}$ | $V_{3}$ | $V_{4}$ | m-LB | ALM
---|---|---|---|---|---|---|---|---|---
Germany | 0.01 | -0.24 | 11.37 | 2.51 | 1.52 | 1.22 | 1.79 | 0.00 | 0.00
Switzerland | 0.01 | -0.43 | 12.59 | 1.83 | 0.83 | 0.88 | 0.95 | 0.00 | 0.00
Austria | -0.01 | -0.41 | 10.26 | 4.78 | 2.15 | 1.81 | 3.27 | 0.00 | 0.00

Note: columns $V_{i}$, $i=1,\ldots,4$ represent the sample variance of each
long series on four equal-length sub-periods splitted from the whole period
01/02/2006 to 08/10/2023. The column ALM represents the p-value of the ALM
test with the maximum lag order being 10; the column m-LB represents the
p-value of the m-LB test at the lag order being 20; 0.00 indicates the p-value
is less than $3\times 10^{-16}$ for these two tests. Skew. and Kurt. represent
the skewness and kurtosis respectively.

To show the fluctuations behind the Ukraine index and GPRD index, we present
two plots in Figs. 4 and 5 below. As one can see from there, these two indices
fluctuate severely around the beginning of 2022 which corresponds with the
real-world event. Later, we attempt to use this information to forecast the
stock market volatility of three countries.

Figure 1: The log-returns of Germany’s national stock market index from
01/02/2006 to 08/10/2023. Figure 2: The log-returns of Switzerland’s national
stock market index from 01/02/2006 to 08/10/2023. Figure 3: The log-returns of
Austria’s national stock market index from 01/02/2006 to 08/10/2023. Figure 4:
The exogenous variable Ukraine index from 01/01/2021 to 28/02/2023. This is
corresponding to the time period of our short data . Figure 5: The exogenous
variable GPRD index from 01/02/2006 to 08/10/2023.

###  5.2 Empirical results

We first consider the forecasting exercise with the short period (1st January,
2021 to 28th February, 2023) data described in Section 5.1. Then, the analysis
of three long returns series is given in Section 5.2.2.

####  5.2.1 Short data

To compare the performance of GARCH, GARCHX and GAX-NoVaS methods, we still
apply the time aggregated prediction metric described in Section 4.1. We
consider $h=1,5,20$ and use a 250-size moving window, which is about 1 year of
daily data. We start the empirical analysis with the short period of data and
we take the Ukraine index as the exogenous predictor. The MSPE results are
summarized in Table 4.

Table 4: MSPE ratios of different methods on three short datasets.

| S-Germany | S-Austria | S-Switzerland
---|---|---|---
Prediction steps | 1 | 5 | 20 | 1 | 5 | 20 | 1 | 5 | 20
GA | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
GAX-Ukraine | 1.011 | 1.023 | 0.989 | 1.002 | 1.003 | 0.994 | 1.003 | 1.008 | 1.005
GAX-NoVaS-Ukraine | 1.038 | 1.141 | 0.898 | 1.017 | 0.984 | 0.822 | 1.010 | 1.100 | 0.879

Note: “S” represents “short”. To simplify the presentation, we compute the
ratio of MSPE of different methods, i.e., we divide all MSPE for each
prediction horizon and model by the MSPE of the GARCH (benchmark) method.

From Table 4, we can see that the GAX-NoVaS method can bring some large
improvements, especially for long-horizon aggregated predictions. Meanwhile,
it seems that the GARCHX and GARCH models have indistinguishable performance.
However, the GAX-NoVaS method is generally better than both GARCH-type
methods. The DM-test results on comparing GARCHX and GAX-NoVaS for forecasting
short real-world data are tabularized in Table 5, which reveals the
significant advantage of GAX-NoVaS for long-horizon predictions, especially
for the 20-step-ahead predictions of Short Austria data.

Table 5: DM-test results on predictions of short datasets.

| S-Germany | S-Austria | S-Switzerland
---|---|---|---
Prediction steps | 1 | 5 | 20 | 1 | 5 | 20 | 1 | 5 | 20
GARCHX-Ukraine |  |  |  |  |  |  |  |  |
GAX-NoVaS-Ukraine | 0.616 | 0.727 | 0.306 | 0.578 | 0.410 | 0.030 | 0.549 | 0.860 | 0.118

Note: The values in the different rows are one-sided $p$-values of the DM-test
on the prediction error of GAX-NoVaS and GARCHX. For the DM-test here, the
alternative hypothesis is that the former method GAX-NoVaS is more accurate
than the latter one GARCHX.

####  5.2.2 Long data

We continue our real data analysis with long datasets (2nd January, 2006 to
10th August, 2023). We also apply more exogenous variables. Similar to the
analysis procedure for short data, we present MSPE ratios and corresponding
DM-test results of GARCHX and GAX-NoVaS in Tables 6 and 7. Generally speaking,
the GAX-NoVaS method still dominates the other two GARCH-type methods, and
this superiority is verified to be significant by the DM-test.

Table 6: MSPE ratios of different methods on three long datasets.

| L-Germany | L-Austria | L-Switzerland
---|---|---|---
Prediction steps | 1 | 5 | 20 | 1 | 5 | 20 | 1 | 5 | 20
GA | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
GA-NoVaS | 0.990 | 0.922 | 0.518 | 1.062 | 1.041 | 0.715 | 1.010 | 0.786 | 0.128
GAX-GPRD | 1.010 | 1.024 | 0.998 | 0.999 | 0.990 | 0.964 | 1.003 | 0.994 | 0.797
GAX-NoVaS-GPRD | 1.000 | 0.953 | 0.557 | 1.076 | 1.121 | 0.831 | 1.005 | 0.799 | 0.133
GAX-Inflation | 1.003 | 1.016 | 0.974 |  |  |  | 0.996 | 0.965 | 0.858
GAX-NoVaS-Inflation | 1.001 | 0.949 | 0.551 |  |  |  | 1.010 | 0.801 | 0.131
GAX-Trend | 0.987 | 0.979 | 1.062 | 0.982 | 0.912 | 0.654 | 0.996 | 0.937 | 0.594
GAX-NoVaS-Trend | 0.998 | 0.949 | 0.540 | 1.077 | 1.169 | 0.866 | 1.019 | 0.828 | 0.138

Note: “L” represents “long”. To simplify the presentation, we compute the
ratio of MSPE of different methods, i.e., we divide all MSPE for each
prediction horizon and model by the MSPE of the GARCH (benchmark) method.

Table 7: DM-test results on predictions of long datasets.

| L-Germany | L-Austria | L-Switzerland
---|---|---|---
Prediction steps | 1 | 5 | 20 | 1 | 5 | 20 | 1 | 5 | 20
GARCHX-GPRD |  |  |  |  |  |  |  |  |
GAX-NoVaS-GPRD | 0.427 | 0.306 | 0.012 | 0.950 | 0.957 | 0.068 | 0.516 | 0.055 | 0.008
GARCHX-Inflation |  |  |  |  |  |  |  |  |
GAX-NoVaS-Inflation | 0.484 | 0.298 | 0.007 |  |  |  | 0.595 | 0.082 | 0.029
GARCHX-Trend |  |  |  |  |  |  |  |  |
GAX-NoVaS-Trend | 0.606 | 0.420 | 0.017 | 0.995 | 1.000 | 1.000 | 0.659 | 0.122 | 0.027

Note: The values in the different rows are one-sided $p$-values of the DM-test
on the prediction error of GAX-NoVaS and GARCHX. For the DM-test here, the
alternative hypothesis is that the former method GAX-NoVaS is more accurate
than the latter one GARCHX.

##  6 Conclusion

We extend the current NoVaS prediction method to the realm of prediction with
exogenous variables. We provide the theoretical foundation to guarantee the
feasibility of applying Model-free prediction. Inspired by the GARCHX model,
we propose a specifically designed model-free/model-based method namely GAX-
NoVaS prediction. The dominance of GAX-NoVaS on the classical GAX-NoVaS method
is verified by simulation and empirical datasets. Also, such an advantage is
not only exhibited by some MSPE metric but we also show some statistical
significance through the parlance of classical DM tests.

We should also mention that going far beyond GAX-NoVaS method might have
limitations if the model becomes increasingly complex. Recall that the great
performance of the GAX-NoVaS method relies on a successful transformation, the
satisfied transformation may not be achievable if the underlying time series
is very complicated. However, there is a growing literature on forecasting
using more non-parametric neural network based models. It will be an
interesting future work that combines the idea of model-free prediction with
the state-of-the-art machine learning method, such as Deep neural network
(DNN), convolutional neural network (CNN) or LSTM etc. Finally, in the field
of binary/categorical/count data INGARCH models have recently garnered
significant attention both from theoretical and applied researchers. One could
potentially also think of a Model-free INGARCH-X type model and try to
integrate the idea of exogenous covariates into it and build a new forecasting
framework to challenge the existing ones. In short, our paper remains the
first paper to propose this idea of model-free transformation-based
forecasting focused on GARCHX type models but the scope of extending this to
several directions is ample.

## References

  * Angus, (1994) Angus, J. E. (1994).  The probability integral transform and related results.  SIAM review, 36(4):652–654.
  * Bollerslev, (1986) Bollerslev, T. (1986).  Generalized autoregressive conditional heteroskedasticity.  Journal of Econometrics, 31(3):307–327.
  * Bouri et al., (2020) Bouri, E., Cepni, O., Gupta, R., and Jalkh, N. (2020).  Geopolitical risks and stock market volatility in the g7 countries: A century of evidence from a time-varying nonparametric panel data model.  Edited by Atin Basuchoudhary and Gunther G. Schulze.
  * Box et al., (2015) Box, G. E., Jenkins, G. M., Reinsel, G. C., and Ljung, G. M. (2015).  Time Series Analysis: Forecasting and Control.  John Wiley & Sons.
  * Caldara and Iacoviello, (2022) Caldara, D. and Iacoviello, M. (2022).  Measuring geopolitical risk.  American Economic Review, 112(4):1194–1225.
  * Chen and Politis, (2019) Chen, J. and Politis, D. N. (2019).  Optimal Multi-step-ahead Prediction of ARCH/GARCH Models and NoVaS Transformation.  Econometrics, 7(3):1–23.
  * Christos et al., (2019) Christos, B., Christina, C., Rangan, G., and Tahir, S. (2019).  Geopolitical risks, returns, and volatility in emerging stock markets: Evidence from a panel garch model.  Emerging Markets Finance and Trade, 55(8):1841–1856.
  * Chudý et al., (2020) Chudý, M., Karmakar, S., and Wu, W. B. (2020).  Long-term prediction intervals of economic time series.  Empirical Economics, 58(1):191–222.
  * Diebold and Mariano, (2002) Diebold, F. X. and Mariano, R. S. (2002).  Comparing predictive accuracy.  Journal of Business & economic statistics, 20(1):134–144.
  * Engle and Patton, (2001) Engle, R. and Patton, A. (2001).  What good is a volatility model?  Quantitative Finance, 1:47–63.
  * Engle, (1982) Engle, R. F. (1982).  Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation.  Econometrica: Journal of the econometric society, pages 987–1007.
  * Engle et al., (2013) Engle, R. F., Ghysels, E., and Sohn, B. (2013).  Stock market volatility and macroeconomic fundamentals.  Review of Economics and Statistics, 95(3):776–797.
  * Engle and Patton, (2007) Engle, R. F. and Patton, A. J. (2007).  What good is a volatility model?  In Forecasting volatility in the financial markets, pages 47–63. Elsevier.
  * Francq et al., (2019) Francq, C. et al. (2019).  Qml inference for volatility models with covariates.  Econometric Theory, 35(1):37–72.
  * Fryzlewicz et al., (2008) Fryzlewicz, P., Sapatinas, T., Rao, S. S., et al. (2008).  Normalized least-squares estimation in time-varying ARCH models.  The Annals of Statistics, 36(2):742–786.
  * González-Rivera et al., (2004) González-Rivera, G., Lee, T.-H., and Mishra, S. (2004).  Forecasting volatility: A reality check based on option pricing, utility function, value-at-risk, and predictive likelihood.  International Journal of Forecasting, 20(4):629–645.
  * Grebe et al., (2024) Grebe, M., Kandemir, S., and Tillmann, P. (2024).  Uncertainty about the war in ukraine: Measurement and effects on the german economy.  Journal of Economic Behavior & Organization, 217(C):493–506.
  * Gulay and Emec, (2018) Gulay, E. and Emec, H. (2018).  Comparison of forecasting performances: Does normalization and variance stabilization method beat garch (1, 1)-type models? empirical evidence from the stock markets.  Journal of Forecasting, 37(2):133–150.
  * Herrera et al., (2018) Herrera, A. M., Hu, L., and Pastor, D. (2018).  Forecasting crude oil price volatility.  International Journal of Forecasting, 34(4):622–635.
  * Karmakar et al., (2022) Karmakar, S., Chudý, M., and Wu, W. B. (2022).  Long-term prediction intervals with many covariates.  Journal of Time Series Analysis, 43(4):587–609.
  * Karmakar and Roy, (2021) Karmakar, S. and Roy, A. (2021).  Bayesian Modelling of Time-Varying Conditional Heteroscedasticity.  Bayesian Analysis, 1(1):1–29.
  * Lim and Sek, (2013) Lim, C. M. and Sek, S. K. (2013).  Comparing the Performances of GARCH-type Models in Capturing the Stock Market Volatility in Malaysia.  Procedia Economics and Finance, 5:478–487.
  * Mehmet et al., (2018) Mehmet, B., Rangan, G., Christian, P., and E, W. M. (2018).  Terror attacks and stock-market fluctuations: Evidence based on a nonparametric causality-in-quantiles test for the g7 countries.  The European Journal of Finance, 24(4):333–346.
  * Peters, (2001) Peters, J.-P. (2001).  Estimating and forecasting volatility of stock indices using asymmetric GARCH models and (Skewed) Student-t densities.  Preprint, University of Liege, Belgium, 3:19–34.
  * Politis, (2003) Politis, D. N. (2003).  A normalizing and variance-stabilizing transformation for financial time series.  In Recent Advances and Trends in Nonparametric Statistics, pages 335–347. Elsevier Inc.
  * Politis, (2015) Politis, D. N. (2015).  The model-free prediction principle.  In Model-Free Prediction and Regression, pages 13–30. Springer.
  * R Core Team, (2023) R Core Team (2023).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Salisu et al., (2022) Salisu, A. A., Ogbonna, A. E., Lasisi, L., and Olaniran, A. (2022).  Geopolitical risk and stock market volatility in emerging markets: A garch–midas approach.  The North American Journal of Economics and Finance, 62:101755.
  * Sucarrat, (2020) Sucarrat, G. (2020).  garchx: Flexible and robust garch-x modelling.
  * Wang and Politis, (2022) Wang, Y. and Politis, D. N. (2022).  Model-free bootstrap for a general class of stationary time series.  Bernoulli, 28(2):744–770.
  * Wu and Karmakar, (2021) Wu, K. and Karmakar, S. (2021).  Model-free time-aggregated predictions for econometric datasets.  Forecasting, 3(4):920–933.
  * Wu and Karmakar, (2023) Wu, K. and Karmakar, S. (2023).  A model-free approach to do long-term volatility forecasting and its variants.  Financial Innovation, 9(1):1–38.
  * Wu and Politis, (2023) Wu, K. and Politis, D. N. (2023).  Bootstrap prediction inference of non-linear autoregressive models.  arXiv preprint arXiv:2306.04126.
  * Wuertz et al., (2013) Wuertz, D., RUnit, S., and Chalabi, M. Y. (2013).  Package ‘fgarch’.  Technical report, Technical report, working paper/manual, 09.11. 2009. URL http://cran. r ….
  * Ypma et al., (2014) Ypma, J., Borchers, H. W., and Eddelbuettel, D. (2014).  nloptr: R interface to nlopt.  R package version, 1(4).
  * Yule, (1927) Yule, G. U. (1927).  Vii. on a method of investigating periodicities disturbed series, with special reference to wolfer’s sunspot numbers.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 226(636-646):267–298.
  * Zhang et al., (2023) Zhang, Y., He, J., He, M., and Li, S. (2023).  Geopolitical risk and stock market volatility: A global perspective.  Finance Research Letters, 53:103620.