# Towards Interpreting Zoonotic Potential of Betacoronavirus Sequences With
Attention

Kahini Wadhawan IBM Research, New Delhi, India<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>Payel Das IBM Research, Yorktown Heights, NY, USA
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>Barbara A. Han Cary Institute of
Ecosystem Studies, NY, USA  Ilya R. Fischhoff Cary Institute of Ecosystem
Studies, NY, USA  Adrian C. Castellanos Cary Institute of Ecosystem Studies,
NY, USA  Arvind Varsani The Biodesign Institute, Arizona State University, USA
Kush R. Varshney IBM Research, Yorktown Heights, NY, USA

###### Abstract

Current methods for viral discovery target evolutionarily conserved proteins
that accurately identify virus families but remain unable to distinguish the
zoonotic potential of newly discovered viruses. Here, we apply an attention-
enhanced long-short-term memory (LSTM) deep neural net classifier to a highly
conserved viral protein target to predict zoonotic potential across
betacoronaviruses. The classifier performs with a 94% accuracy. Analysis and
visualization of attention at the sequence and structure-level features
indicate possible association between important protein-protein interactions
governing viral replication in zoonotic betacoronaviruses and zoonotic
transmission.

##  1 Introduction

The majority of viruses emerging in humans arise from animal hosts. For
example, SARS-CoV-2, the virus causing the COVID-19 pandemic, belongs to the
betacoronavirus($\beta$-CoVs) family from which numerous other viruses have
been described from multiple mammal species around the world. The typical
method for surveiling for such viruses is to identify a conserved region of
the genome. For instance, the RdRp (RNA-dependent RNA polymerase) gene is a
highly conserved sequence that is commonly used to ascertain $\beta$-CoVs
presence from sampled wildlife species. While this conserved nature makes RdRp
sequences a reliable target for accurately identifying novel $\beta$-CoVs from
field surveillance, it does not enable predictions about whether the virus
poses a zoonotic threat to humans. This is because the RdRp sequences encode
the polymerase responsible for replicating the RNA genome once inside the
cell, which is not involved in host cell entry or the infection process. The
disconnect between broad viral surveillance and accurate predictions about
zoonotic potential of newly discovered viruses remains a major research
frontier, which has been highlighted by the emergence of SARS-CoV-2. Numerous
closely-related $\beta$-CoVs have been identified from field surveillance over
the past several years, but no clear and cost-effective way to triage the
zoonotic risk posed by these newly discovered viruses. Beyond quantifying the
zoonotic risk posed by a novel virus, preventing their emergence also relies
on differentiating the key elements that enable some viruses to successfully
infect humans. In other words, explaining what viral structures enable high
zoonotic potential will be critical for informing what downstream actions will
be effective for preventing spillover and viral emergence.

Recent studies have shown stronger performance of machine learning based
approaches for zoonotic transmission prediction task, when compared to
homology based approaches, both at the protein sequence and gene sequence
level (Eng et al., 2017; Xu et al., 2017; Li & Sun, 2018; Zhang et al., 2019;
Mock et al., 2019; Wang et al., 2013; Bartoszewicz et al., 2021). There is a
consensus that deep models outperform shallow classifiers. Previous works have
often used either feature vectors reflecting compositional and physicochemical
properties of protein sequences or on word vectors. In a few studies, input
features were augmented with structural details, e.g. (Fischhoff et al.,
2021). In this study, we develop a deep neural net classifier for predicting
zoonotic potential of RdRp sequences. For this purpose, we employ an attention
mechanism enhanced long short-term memory (LSTM) model (Hochreiter &
Schmidhuber, 1997) trained on RdRp sequences, as attention Bahdanau et al.
(2014) is known to better cope with longer length sequences than a simple LSTM
model, as is the case here. The final model yields a $94\%$ prediction
accuracy.

For wide-spread adoption of deep learning models for high-stake tasks, such as
zoonotic transmission prediction, it is important that such models are
interpretable and explainable. Towards this direction, we intend to interpret
the attention signals of the intermediate layers of deep classifier by
visualizing the attention on the amino acid characters of the sample
sequences. We also investigate average attention difference between two
classes as a function of the amino acid type. Finally, we visualize the
learned attentions on the 3D RdRp structure, as the 3D structure determines
its function. Analysis of these internal workings of the classifier model
provides important insights into likely mechanism of zoonotic transmission, as
class-level attention differences map to the key protein-protein interaction
sites. To our knowledge, this is the first machine learning work that provides
sequence and structure-level understanding of zoonotic transmission potential
of protein sequences.

##  2 Data: RdRp Sequence Data and Label Annotation

We assembled a dataset of published RdRp gene sequences for all $\beta$-CoVs,
freely available in GenBank (Benson et al., 2013). For each published
sequence, we added a binary label to designate zoonotic status based on the
results of published research. The potential zoonotic sequences included ones
similar or identical to sequences that caused MERS coronavirus (MERS-CoV)
SARS. In the NCBI records, “SARS” was often used to describe these viruses;
here we refer to SARS-CoV-1 to distinguish from SARS-CoV-2 causing the
COVID-19 pandemic. MERS-CoV or SARS-CoV-1 sequences found in humans were
categorized as zoonotic. If no host species was identified for a sequence, we
did not categorize zoonotic status for that sequence. Prior to modeling, all
sequences for which we were unable to categorize zoonotic status were removed
from the dataset. For sequences passaged or adapted through cell culture of a
different species than the original host, we made no determination about
zoonotic status, as passaging may have changed zoonotic capacity. These
included viruses described as involving Vero cell culture, from African green
monkey (Cercopithecus aethiops) cells (Ammerman et al., 2008; Tsui et al.,
2003); mouse-adapted (Lu et al., 2009); recombinant; or LLC-MK2, cultured in
rhesus monkey (Macaca mulatta) cells (e.g., Kaye et al. (2006)). Information
about these forms of cell culture was often found in the host field in NCBI;
or as part of “passage details” in the notes field. If a sequence originated
in one person and was cultured in cells from a different human cell line prior
to sequencing, we considered it to remain zoonotic as we did not expect
differences among people to change the zoonotic status.

We also considered sequences found in animals to be zoonotic if the authors of
the sequence considered it to be functionally identical to sequences found in
people. For MERS-CoV, we emphasized virus taxonomic clade position and
locality, which is related to clade position, as for MERS-CoV this information
is more readily available than information about sequence similarity to human
MERS-CoV. MERS-CoV in clades A and B include known human sequences and are
known to cross from camels (Camelus dromedarius) to humans (Lau et al., 2016;
Chan et al., 2014). For ambiguous MERS-CoV sequences, we categorized as
zoonotic if there were human cases found in the same geographic area. MERS-CoV
in Africa, in clade C, are similar but with no known human cases Chu et al.
(2018). Clade C has lower titer and virulence in lab experiments, compared to
clade A or B . We assumed that clade C MERS-CoV to be non-zoonotic. For SARS-
CoV-1 viruses found in animals, if insufficient information was available to
confirm that the virus was zoonotic then it was labeled non-zoonotic. For
example, SARS-CoV-1 viruses found in bats were categorized as non-zoonotic if
no study was published associated with the sequence. SARS-CoV-1 reported in
palm civet (Paguma larvata) were considered zoonotic if the authors considered
it to be identical to sequences found in people Wang et al. (2005), and
otherwise categorized as non-zoonotic (Yi, 2005). We considered as zoonotic a
sequence isolated from an animal that also replicated in human cell lines,
e.g. a SARS-CoV-1 from Rhinolophus sinicus bats that replicated in human cell
culture (Ge et al., 2013). If the sequence was found in an animal species not
known to be a bridge species to people, e.g. a ferret, we categorized it as
non-zoonotic (Town et al., 2012).

##  3 Methods

The resulting dataset contained a total of 4259 sequences of length up to 910
characters, with vocabulary of 21 amino acids (A, C, D, E, F, G, H, I, K, L,
M, N, P, Q, R, S, T, V, W, X, Y). We use 60/20/20 ratio to create
train/val/test split, resulting in 2548, 876 and 835 sequences respectively,
maintaining the positive to negative ratio in each split. A LSTM model
enhanced with a single head attention mechanism was trained on this RdRp
sequence dataset to predict their zoonotic potential (as a binary
classification task). The model architecture was comprised of a single encoded
layer with single head attention Bahdanau et al. (2014). The model accepts in
the sequence of amino acids $x=(x_{1},x_{2},...,x_{L})$ and learns an
embedding $z=(z_{1},z_{2},...,z_{L})$, which is passed to a dropout layer
followed by a RELU and a softmax layer. Best test accuracy of $94\%$ was
achieved after hyper-parameter search. More details are in SI section A.1

##  4 Results

Zoonotic classification accuracy. In this study, we use an attention-enhanced
classifier, as attention has been proved to be helpful for dealing with longer
sequences (Bahdanau et al., 2016). In our dataset, the max sequence length is
$910$. Figure 2(a) shows the color coded length distribution of sequences for
both positive and negative classes. The dataset is very length skewed. We see
that the majority of zoonotic sequences are of length $>800$, while the
shorter length $<200$ are mostly non-zoonotic. Despite the classification
challenges presented by this skew, the best performing model achieved an
overall accuracy of $94\%$. Figure 1 shows confusion matrices for different
length ranges. For shorter length sequences ($len<150$), the classifier
correctly predicted $>92\%$ of the true negatives, but missed a small number
of true positives. Improved predictions of true positives was achieved with
increase in sequence length. F1-scores are reported in SI Table 1 and section
A.1.3.

(a) $len<150$

(b) $150<len<500$

(c) $len>500$

(d) all length

Figure 1: Model performance confusion matrix computed over range of sequence
lengths.

(a) Length distribution

(b) Avg. attention difference vs residue number

(c) Feature importance (avg. attention over amino acids)

Figure 2: Data statistics and attention at 3D and 1D level plots. Figure 3:
$Left:$ Structure of SARS-CoV nsp12 RdRp (dark green) bound to nsp7 (light
green) and nsp8 (pink) co-factors (pdb id:6NUR). Metal ions shown as blue
spheres. $Middle:$ Average attention differences between two classes, mapped
to the SARS-CoV nsp12 structure (6NUR.pdb chain A residue 146 to 714). Red-
gray-blue (low to high) colorscale is used to visualize the attention heatmap.
$Right:$ Regions of RdRp interacting (cutoff: 6.5 Å) with nsp7 highlighted in
magenta.

Interpretation using Attention Visualization. Mapping predictions back to a
target sequence has been used both as a means of investigating a given model’s
performance and as a method of sequence analysis. For example, convolutional
filter visualization has been used for genome sequence classification
(Bartoszewicz et al., 2021). Protein sequence position importance was
investigated by an attention function in subcellular localization
classification task (Jurtz et al., 2017). Recently, 3D visualization of
attention of a Transformer protein sequence model revealed accurate capture of
structural and binding properties (Vig et al., 2020). Similarly, in this
study, we visualize and identify regions of 1D protein sequence and 3D protein
structure that are important for the prediction of zoonotic potential. Figure
2(c) shows averaged normalized attentions for each amino acid character,
estimated by averaging over all sequences within the same class. At the amino
acid level, significant differences in attention are evident between two
classes. For example, relatively higher attention in positive sequences is
found in the amino acids D, Y, R, F, and Q.

Since the 3D structure of the protein determines its function, we project the
average attention heatmap differences between the two classes on the publicly
available SARS-CoV nsp12 RdRp structure bound to the partnering nsp7 and nsp12
cofactor proteins (Kirchdoerfer & Ward, 2019)(Figure 3). Figure 2(b) show
these attention differences as a function of residue number of nsp12 template.
For selection criteria refer SI section A.2.1. Figure 5 further highlights
these differences mapped on to the homology model (using SARS-CoV nsp12 as
template) of individual sequences from the positive and negative class. These
visualizations clearly show higher attentions within the positive sequences
map to the N-terminal structural regions that correspond to the nsp12-nsp7
interaction surface. This result is interesting in light of prior work showing
that the assembly of RdRp with cofactor proteins is crucial for polymerase
activity that governs viral replication (Yin et al., 2020; Subissi et al.,
2014). In fact, protein-protein interaction design has been shown as a proven
approach for RdRp inhibition in other related viruses (Subissi et al., 2014).

##  5 Discussion and Future Work

We present an attention-enhanced LSTM deep neural net classifier that predicts
zoonotic transmission in $\beta$-CoVs. Crucially, this classification was
achieved using widely available RdRp sequences. Accurate prediction of
zoonotic potential from RdRp sequences in $\beta$-CoVs (or similarly used
reliable primer sequences for other viral families) has the potential to
alleviate a critical bottleneck in viral surveillance, which has so far been
unable to infer zoonotic potential of novel viruses discovered as part of
routine surveillance. Further investigation of attention maps at both the
sequence and the structure level provided important interpretability of the
“black-box” model. Our results indicate mapping of class-level attention
differences at protein-protein interaction sites of the RdRp structure. Those
sites are of significant functional importance for binding with cofactors
crucial for polymerase activity involved in viral replication. Thus, while
RdRp itself is not directly involved in the infection process, physiochemical
interactions with neighbouring cofactors appear to influence the function of
zoonotic $\beta$-CoVs compared to non-zoonotic counterparts. Future work will
explore more sophisticated attention mechanism; multi-head (e.g. Transformers
(Devlin et al., 2019) (Vaswani et al., 2017)) for better interpretability and
prediction. We also plan extend the current study beyond $\beta$-CoVs to
investigate the robustness in capturing of biological context via model
features.

## References

  * Altschul et al. (1990) S. Altschul, W. Gish, W. Miller, E. Myers, and D. Lipman.  Basic local alignment search tool.  _Journal of Molecular Biology_ , 215:403–410, 1990.
  * Ammerman et al. (2008) Nicole C Ammerman, Magda Beier-Sexton, and Abdu F Azad.  Growth and maintenance of Vero cell lines.  _Current Protocols in Microbiology_ , 11(1):A.4E.1–A.4E.7, 2008.
  * Bahdanau et al. (2014) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate.  _arXiv preprint arXiv:1409.0473_ , 2014.
  * Bahdanau et al. (2016) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate, 2016.
  * Bartoszewicz et al. (2021) Jakub M Bartoszewicz, Anja Seidel, and Bernhard Y Renard.  Interpretable detection of novel human viruses from genome sequencing data.  _NAR genomics and bioinformatics_ , 3(1):lqab004, 2021.
  * Benson et al. (2013) Dennis A Benson, Mark Cavanaugh, Karen Clark, Ilene Karsch-Mizrachi, David J Lipman, James Ostell, and Eric W Sayers.  GenBank.  _Nucleic Acids Res._ , 41(Database issue):D36–42, January 2013.
  * Berman et al. (2000) Helen M. Berman, John Westbrook, Zukang Feng, Gary Gilliland, T. N. Bhat, Helge Weissig, Ilya N. Shindyalov, and Philip E. Bourne.  The Protein Data Bank.  _Nucleic Acids Research_ , 28(1):235–242, 01 2000\.  ISSN 0305-1048.  doi: 10.1093/nar/28.1.235.  URL https://doi.org/10.1093/nar/28.1.235.
  * Chan et al. (2014) Renee WY Chan, Maged G Hemida, Ghazi Kayali, Daniel KW Chu, Leo LM Poon, Abdelmohsen Alnaeem, Mohamed A Ali, Kin P Tao, Hoi Y Ng, Michael CW Chan, Yi Guan, John M Nicholls, and JS Malik Peiris.  Tropism and replication of Middle East respiratory syndrome coronavirus from dromedary camels in the human respiratory tract: an in-vitro and ex-vivo study.  _The Lancet Respiratory Medicine_ , 2(10):813–822, 2014.
  * Chu et al. (2018) Daniel K. W. Chu, Kenrie P. Y. Hui, Ranawaka A. P. M. Perera, Eve Miguel, Daniela Niemeyer, Jincun Zhao, Rudragouda Channappanavar, Gytis Dudas, Jamiu O. Oladipo, Amadou Traore, Ouafaa Fassi-Fihri, Abraham Ali, Getnet F. Demissie, Doreen Muth, Michael C. W. Chan, John M. Nicholls, David K. Meyerholz, Sulyman A. Kuranga, Gezahegne Mamo, Ziqi Zhou, Ray T. Y. So, Maged G. Hemida, Richard J. Webby, Francois Roger, Andrew Rambaut, Leo L. M. Poon, Stanley Perlman, Christian Drosten, Veronique Chevalier, and Malik Peiris.  MERS coronaviruses from camels in Africa exhibit region-dependent genetic diversity.  _Proceedings of the National Academy of Sciences_ , 115(12):3144–3149, 2018.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pp. 4171–4186, Minneapolis, Minnesota, June 2019. Association for Computational Linguistics.  doi: 10.18653/v1/N19-1423.  URL https://www.aclweb.org/anthology/N19-1423.
  * Eng et al. (2017) Christine LP Eng, Joo Chuan Tong, and Tin Wee Tan.  Predicting zoonotic risk of influenza a viruses from host tropism protein signature using random forest.  _International journal of molecular sciences_ , 18(6):1135, 2017.
  * Fischhoff et al. (2021) Ilya R Fischhoff, Adrian A Castellanos, Joao PGLM Rodrigues, Arvind Varsani, and Barbara A Han.  Predicting the zoonotic capacity of mammal species for sars-cov-2.  _bioRxiv_ , 2021.
  * Galassi et al. (2020) Andrea Galassi, Marco Lippi, and Paolo Torroni.  Attention in natural language processing.  _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–18, 2020.  ISSN 2162-2388.  doi: 10.1109/tnnls.2020.3019893.  URL http://dx.doi.org/10.1109/TNNLS.2020.3019893.
  * Ge et al. (2013) Xing-Yi Ge, Jia-Lu Li, Xing-Lou Yang, Aleksei A. Chmura, Guangjian Zhu, Jonathan H. Epstein, Jonna K. Mazet, Ben Hu, Wei Zhang, Cheng Peng, Yu-Ji Zhang, Chu-Ming Luo, Bing Tan, Ning Wang, Yan Zhu, Gary Crameri, Shu-Yi Zhang, Lin-Fa Wang, Peter Daszak, and Zheng-Li Shi.  Isolation and characterization of a bat SARS-like coronavirus that uses the ACE2 receptor.  _Nature_ , 503(7477):535–538, 2013.
  * Hochreiter et al. (2001) S. Hochreiter, Y. Bengio, P. Frasconi, and J. Schmidhuber.  Gradient flow in recurrent nets: the difficulty of learning long-term dependencies.  In S. C. Kremer and J. F. Kolen (eds.), _A Field Guide to Dynamical Recurrent Neural Networks_. IEEE Press, 2001.
  * Hochreiter & Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural Computation_ , 9(8):1735–1780, 1997.
  * Jurtz et al. (2017) Vanessa Isabell Jurtz, Alexander Rosenberg Johansen, Morten Nielsen, Jose Juan Almagro Armenteros, Henrik Nielsen, Casper Kaae Sønderby, Ole Winther, and Søren Kaae Sønderby.  An introduction to deep learning on biological sequence data: examples and solutions.  _Bioinformatics_ , 33(22):3685–3690, 2017.
  * Kaye et al. (2006) Matthew Kaye, Julian Druce, Thomas Tran, Renata Kostecki, Doris Chibo, Jessica Morris, Mike Catton, and Chris Birch.  SARS-associated coronavirus replication in cell lines.  _Emerging Infectious Diseases_ , 12(1):128–133, 2006.
  * Kirchdoerfer & Ward (2019) Robert N Kirchdoerfer and Andrew B Ward.  Structure of the sars-cov nsp12 polymerase bound to nsp7 and nsp8 co-factors.  _Nature communications_ , 10(1):1–9, 2019.
  * Lau et al. (2016) Susanna KP Lau, Renate Wernery, Emily YM Wong, Sunitha Joseph, Alan KL Tsang, Nissy Annie Georgy Patteril, Shyna K Elizabeth, Kwok-Hung Chan, Rubeena Muhammed, Jöerg Kinne, Kwok-Yung Yuen, Ulrich Wernery, and Patrick CY Woo.  Polyphyletic origin of MERS coronaviruses and isolation of a novel clade A strain from dromedary camels in the United Arab Emirates.  _Emerging Microbes and Infections_ , 5(12):e128, 2016.
  * Li & Sun (2018) Han Li and Fengzhu Sun.  Comparative studies of alignment, alignment-free and svm based approaches for predicting the hosts of viruses based on viral sequences.  _Scientific reports_ , 8(1):1–9, 2018.
  * Lu et al. (2009) Baojing Lu, Ling Tao, Ting Wang, Zhenhua Zheng, Bao Li, Ze Chen, Yi Huang, Qinxue Hu, and Hanzhong Wang.  Humoral and cellular immune responses induced by 3a DNA vaccines against severe acute respiratory syndrome (SARS) or SARS-like coronavirus in mice.  _Clinical and Vaccine Immunology_ , 16(1):73–77, 2009.
  * Luong et al. (2015) Minh-Thang Luong, Hieu Pham, and Christopher D. Manning.  Effective approaches to attention-based neural machine translation, 2015\.
  * Mock et al. (2019) Florian Mock, Adrian Viehweger, Emanuel Barth, and Manja Marz.  Viral host prediction with deep learning.  _bioRxiv_ , pp. 575571, 2019.
  * Subissi et al. (2014) Lorenzo Subissi, Clara C Posthuma, Axelle Collet, Jessika C Zevenhoven-Dobbe, Alexander E Gorbalenya, Etienne Decroly, Eric J Snijder, Bruno Canard, and Isabelle Imbert.  One severe acute respiratory syndrome coronavirus protein complex integrates processive rna polymerase and exonuclease activities.  _Proceedings of the National Academy of Sciences_ , 111(37):E3900–E3909, 2014.
  * Town et al. (2012) C.D. Town, R.A. Halpin, A. Ransier, N. Fedorova, L. Overton, T. Tsitrin, T. Stockwell, P. Amedeo, B. Bishop, H. Chen, P. Edworthy, N. Gupta, D. Katzel, K. Li, S. Schobel, S. Shrivastava, V. Thovarai, S. Wang, N. Powers, S. Grimes, B. Yamamoto, A. Robb, D.E. Wentworth, and D.L. Noah.  SARS coronavirus Tor2 isolate Tor2/FP1-10851, complete genome, 2012.  https://www.ncbi.nlm.nih.gov/nuccore/JX163924.1 [Accessed: 29 April 2020].
  * Tsui et al. (2003) Stephen KW Tsui, Stephen SC Chim, and YM Dennis Lo.  Coronavirus genomic-sequence variations and the epidemiology of the severe acute respiratory syndrome.  _The New England Journal of Medicine_ , 349(2):187–188, 2003.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need, 2017.
  * Vig et al. (2020) Jesse Vig, Ali Madani, Lav R Varshney, Caiming Xiong, Richard Socher, and Nazneen Fatema Rajani.  Bertology meets biology: Interpreting attention in protein language models.  _arXiv preprint arXiv:2006.15222_ , 2020.
  * Wang et al. (2013) Jia Wang, Chuang Ma, Zheng Kou, Yan-Hong Zhou, and Huai-Lan Liu.  Predicting transmission of avian influenza a viruses from avian to human by using informative physicochemical properties.  _International journal of data mining and bioinformatics_ , 7(2):166–179, 2013.
  * Wang et al. (2005) M Wang, MY Yan, HF Xu, W Liang, B Kan, BJ Zheng, HL Chen, H Zheng, YM Xu, EM Zhang, HX Wang, JR Ye, GC Li, MC Li, ZG Cui, YF Liu, RT Guo, XN Liu, LH Zhan, DH Zhou, AL Zhao, R Hai, DZ Yu, Y Guan, and JG Xu.  SARS-CoV infection in a restaurant from palm civet.  _Emerging Infectious Diseases_ , 11(12):1860, 2005\.
  * Waterhouse et al. (2018) Andrew Waterhouse, Martino Bertoni, Stefan Bienert, Gabriel Studer, Gerardo Tauriello, Rafal Gumienny, Florian T Heer, Tjaart A P de Beer, Christine Rempfer, Lorenza Bordoli, Rosalba Lepore, and Torsten Schwede.  SWISS-MODEL: homology modelling of protein structures and complexes.  _Nucleic Acids Research_ , 46(W1):W296–W303, 05 2018.  ISSN 0305-1048.  doi: 10.1093/nar/gky427.  URL https://doi.org/10.1093/nar/gky427.
  * Xu et al. (2017) Beibei Xu, Zhiying Tan, Kenli Li, Taijiao Jiang, and Yousong Peng.  Predicting the host of influenza viruses based on the word vector.  _PeerJ_ , 5:e3579, 2017.
  * Yi (2005) G. Yi.  SARS coronavirus HC/SZ/61/03, complete genome, 2005.  https://www.ncbi.nlm.nih.gov/nuccore/AY515512.1 [Accessed: 29 April 2020].
  * Yin et al. (2020) Wanchao Yin, Chunyou Mao, Xiaodong Luan, Dan-Dan Shen, Qingya Shen, Haixia Su, Xiaoxi Wang, Fulai Zhou, Wenfeng Zhao, Minqi Gao, et al.  Structural basis for inhibition of the rna-dependent rna polymerase from sars-cov-2 by remdesivir.  _Science_ , 368(6498):1499–1504, 2020.
  * Zhang et al. (2019) Zheng Zhang, Zena Cai, Zhiying Tan, Congyu Lu, Taijiao Jiang, Gaihua Zhang, and Yousong Peng.  Rapid identification of human-infecting viruses.  _Transboundary and emerging diseases_ , 66(6):2517–2522, 2019.

##  Appendix A Appendix

###  A.1 Model Details

####  A.1.1 The Model

For sequential data, Recurrent Neural Networks were introduced but as the
sequence grows, the RNNs, weights could grow beyond control or vanish. Long
Short-Term Memory (Hochreiter & Schmidhuber, 1997) was proposed to deal with
vanishing gradient problem (Hochreiter et al., 2001) and learn long-term
dependency among longer time period. Attention further learns these
dependencies in a smarter and relevant way (Luong et al., 2015; Galassi et
al., 2020; Bahdanau et al., 2014).

LSTM estimates conditional probability $p(y_{1},...,y_{T}|x_{1},...,x_{L})$.
Our task was binary classification, so the input was sequence of characters
($x_{1},...,x_{L}$) and output was binary probability vector among two
classes, positive and negative such that $\sum_{t=1}^{2}p(y_{t})=1$. The LSTM
computes this conditional probability by first obtaining the fixed-dimensional
representation $h$ of the input sequence and given by the last hidden state of
the LSTM. We enhanced our model with attention so, a final hidden attentional
state $h^{\prime}_{t}$ was computed which was passed to a dropout layer, RELU
and followed by softmax to obtain predictive distribution over classes.

For the attention part, we used multiplicative attention (Luong et al., 2015).
All the hidden states of the encoder LSTM $h_{t}$ were used to compute context
vector $c_{t}$. A variable length alignment vector was computed using current
target hidden state $h_{t}$ with each source hidden state $h_{s}$

$\alpha_{t}(s)=align(h_{t},h_{s})=\frac{\exp(score(h_{t},h_{s}))}{\sum_{s^{\prime}}\exp(score(h_{t},h_{s^{\prime}})}$
(1)

Here, $score$ is content-based function for which we used dot product, refer
(Luong et al., 2015) for other alternatives.

$score(h_{t},h_{s})=h_{t}^{T}h_{s}$ (2)

Given the alignment vector as weights, the context vector $c_{t}$ is computed
as the weighted average over all the source hidden states.(Luong et al.,
2015).

We took the dot product of the last hidden state as the target $h_{t}$ and
context vector $c_{t}$ produce attentional hidden state. This attentional
vector $h^{\prime}_{t}$ was then passed to softmax layer get predictive
distribution over the positive and negative class:

$p(y_{t}|y<t,x)=softmax(W_{s},h^{\prime}_{t})$ (3)

####  A.1.2 Training Details

We used a unidirectional LSTM with single layer and single head attention with
1000 cells (max sequence length) and 200 dimensional word embeddings (amino
acid characters in our case) with a vocabulary of 21 amino acids. Hyper-
parameter search was conducted and model giving best accuracy was found with
the configuration of input embedding layer of size 200, lstm hidden layer of
size 200, dropout 0.2, trained for 20 epochs. The baseline LSTM model with no
attention gave $85\%$ accuracy.

####  A.1.3 Model Performance Improvements on balanced Dataset

(a) $len<150$

(b) $150<len<500$

(c) $len>500$

(d) all length

Figure 4: Model performance confusion matrix computed over range of sequence
lengths for balanced dataset refer SI section A.1.3.

Due to the imbalanced class distribution as a function of sequence length,
model performs poorly in lower sequence length regime len $<150$. It seems to
overfit to negative class and hence we saw $0$ True positives in len $len<150$
fig 1. Also, relatively poor performance with mid-length range. To overcome
this problem, we adopted techniques to balance our dataset especially in lower
length regime by oversampling the positive samples. In the imbalanced dataset
we had around $195$ positive sequences with len $<400$ compared to $~{}1100$
negative class sequences (refer figure 2(a)). So, we made 5 copies of each
positive sample with length $<400$ which balances the dataset in lower length
regime. Model was retrained on balanced data using configuration reported in
SI section A.1.2

This improved the F1 score significantly (from 0 to 0.660) in the lower length
regime. Improvements in recall and F1-scores of lower and mid-length ranges
can be seen in table 1. The overall F1-score of the model for all length also
gets improved from $0.946$ to $0.961$. Figure 4 shows confusion matrices of
model performance trained on balanced dataset. In the future work, we will
explore more techniques to balance the dataset.

Sequence length ranges |  Accuracy |  Precision |  Recall |  F1-score
---|---|---|---|---
len$<150$ (imbalanced data) |  0.925 |  NaN(0/0) |  0.0 |  0.0
$150<$len$<500$ |  0.949 |  0.955 |  0.955 |  0.955
len$>500$ |  0.938 |  0.978 |  0.947 |  0.962
all length |  0.935 |  0.977 |  0.916 |  0.946
len$<150$ (balanced data) |  0.880 |  0.971 |  0.500 |  0.660
$150<$len$<500$ |  0.993 |  0.992 |  1.00 |  0.996
len$>500$ |  0.971 |  0.969 |  0.998 |  0.983
all length |  0.948 |  0.973 |  0.949 |  0.961
Table 1: Model performance report using accuracy, precision, recall and
F1-score for both imbalanced skewed dataset and balanced (oversampled positive
datapoints) dataset over all sequence length ranges shows improvement in model
performance by improvement in F1-scores in balanced dataset. Lower length
regime ($len<150$) shows significant increase in F1-score. Interesting
significant improvements are highlighted by bold text. Precision is not
defined(NaN) in len$<150$ for imbalanced data and False negatives $>0$ in this
case so we take F1 to be 0 as model failed to identify any True positives when
it could have done.

###  A.2 Attention Analysis Methods

For interpreting the model, we carried different analysis on the attention
vectors. Visualization of average attention vectors looked into the overall
trends distributed into positive and negative classes at the structure level.
In addition, we also investigated attention at each sample/individual level.

####  A.2.1 Sequence Selection Criteria for Average Attention Analysis

We used attention vectors of the trained LSTM model and carried out detailed
analysis on them. We computed average attention (shown in figure 3) of
filtered positive (zoonotic) and negative (non-zoonotic) classes, to have an
insight into what model learns about structural differences between the two
classes. For mapping to a template structure, positive and negative sequences
that are highly similar to SARS-CoV nsp12 RdRp structure bound to nsp7 and
nsp8 co-factors (pdbid:6NUR) were selected.

Sequence similarity between RdRP sequences and fixed target 6NUR was computed
using NCBI BLAST similarity search (Altschul et al., 1990). Command used is
’blastp’ using an evalue of 10, word size of 3, the BLOSUM62 subsitution
matrix, a gap open penalty of 9, a gap extension penalty of 1, threshold of
16, comp based stats set to 0 and window size of 15. Sequences with
$alignment$ $len>900$, $identities>90\%$, $gaps=0$ and predicted correctly
(true positives and true negatives) by our classifier model are selected. 277
positive and 26 negative sequences are chosen for this analysis shown in
figures 3 and 2(b).

####  A.2.2 Homology Modeling for individual level attention analysis

For visualizing attention at the structural (3D) level of individual positive
and negative sequences (shown in figure 5), their structure prediction was
required as the structures for these surveilled RdRp sequences were not
available in the Protein Data Bank(PDB (Berman et al., 2000)). For the same,
homology modeling was conducted using SWISS-MODEL server (Waterhouse et al.,
2018) using NSP12 Chain A as the template. Figure 3, 6, 7 and 8 display
attention visualization on structure and sequence level for pair of true
positive and true negative long length and short length sequences.

Highly homologous sequences were filtered using $identities>70\%$ and
$gaps=0$. Alignment length $>900$ gave longer length selections and for
shorter length we set alignment length between $200$ to $300$. We employed
homology modeling of these highly similar sequences using Swiss-model server.

(a) Positive (zoonotic) sequence

(b) Negative (non-zoonotic) sequence

Figure 5: Attention heatmap of one positive and one negative sequences with
length $>900$ mapped on to their modeled structure obtained using homology
modeling (6NUR.pdb chain A used as template). Colorscale used: red-gray-blue
(low to high)

(a) Positive(zoonotic) sequence, blast identities$>70\%$, classifier
confidence score$>90\%$

(b) Negative(non-zoonotic) sequence, blast identities$>70\%$, classifier
confidence score$>90\%$

Figure 6: Attention heatmap of one positive and one negative sequences with
length $>900$ at sequence 1D level. Positive with blue and negative with red.

(a) Positive(zoonotic) sequence len$~{}300$, blast identities$>70\%$,
classifier confidence score$>90\%$

(b) Negative(non-zoonotic) sequence len$~{}300$, blast identities$>70\%$,
classifier confidence score$>90\%$

Figure 7: 3D attention heatmap, a pair of shorter length sequences that have
high similarity (identities $>70\%$) with NSP12 chain A. For shorter length,
used relaxed identities criteria as the overall length reduced so no sequences
available with higher $90\%$ identities.

(a) Positive(zoonotic) sequence len$~{}300$, blast identities$>70\%$,
classifier confidence score$>97\%$

(b) Negative(non-zoonotic) sequence len$~{}300$, blast identities$>70\%$,
classifier confidence score$>97\%$

Figure 8: Attention heatmap of one positive and one negative sequences with
length $~{}300$ at sequence 1D level. Positive with blue and negative with
red.


# Expectile hidden Markov regression models for analyzing cryptocurrency
returns

Beatrice Foroni Corresponding author<EMAIL_ADDRESS>MEMOTEF
Department, Sapienza University of Rome  Luca Merlo  Department of Human
Sciences, European University of Rome  Lea Petrella  MEMOTEF Department,
Sapienza University of Rome

###### Abstract

In this paper we develop a linear expectile hidden Markov model for the
analysis of cryptocurrency time series in a risk management framework. The
methodology proposed allows to focus on extreme returns and describe their
temporal evolution by introducing in the model time-dependent coefficients
evolving according to a latent discrete homogeneous Markov chain. As it is
often used in the expectile literature, estimation of the model parameters is
based on the asymmetric normal distribution. Maximum likelihood estimates are
obtained via an Expectation-Maximization algorithm using efficient M-step
update formulas for all parameters. We evaluate the introduced method with
both artificial data under several experimental settings and real data
investigating the relationship between daily Bitcoin returns and major world
market indices.

Keywords: Asymmetric normal distribution, Cryptocurrencies, EM algorithm,
Expectile regression, Markov switching models, Time series

##  1 Introduction

In the last ten years, investors have been increasingly attracted by the
exploit of the cryptocurrency market, mostly because of its peculiar
characteristics. Born merely as a peer-to-peer electronic cash system
(Nakamoto, 2008), the $70$ billion increase in market capitalization (in
particular Bitcoin during 2016-2017), enormous price jumps and levels of high
volatility that were never seen before have made cryptos a new category of
investment assets. Their unusual behavior makes them prone to some speculative
bubbles that may in turn threaten the stability of financial markets (Cheah &
Fry, 2015; Yarovaya et al., 2016). Being crucial to address the level of
integration between cryptocurrencies and traditional financial assets, many
contributions have analyzed the relationship with equities mainly relying on
well-known econometric techniques such as GARCH models (Katsiampa et al.,
2019; Guesmi et al., 2019), variance decomposition (Ji et al., 2018; Corbet et
al., 2018; Yi et al., 2018) and Granger causality test (Bouri et al., 2020b).
Part of the related literature has focused on extreme returns by using models
capturing the tail behaviour, rather than inferring such occurrences from
models based on conditional central tendency. For instance, Kristjanpoller et
al. (2020) and Naeem et al. (2021) employ a multifractal asymmetric analysis,
indicating the presence of heterogeneity in the cross-relationship between
most cryptocurrencies and equity ETFs and showing different behaviors between
upward and downward trends. Shahzad et al. (2022) investigate tail-based
connectedness among major cryptocurrencies in extreme downward and upward
market conditions using LASSO penalized quantile regressions, while Zhang et
al. (2021) apply a risk spillover approach based on generalized quantiles,
showing the existence of a downside risk spillover between Bitcoin and
traditional assets.

In quantitative risk management, indeed, investigating the dynamic of extreme
occurrences is of utmost importance for market participants and regulators.
Among the different methods considered throughout the literature, quantile
regression, introduced by Koenker & Bassett (1978), has represented a valid
approach for modeling the entire distribution of returns while accounting for
the well-known stylized facts, i.e., high kurtosis, skewness and serial
correlation, that typically characterize financial assets. In the financial
literature, the quantile regression framework has been positively applied to
estimate and forecast Value at Risk (VaR) and quantile-based risk measures
(Engle & Manganelli 2004; White et al. 2015; Taylor 2019; Merlo et al. 2021).
Several generalizations of the concept of quantiles have also been introduced
over the years. One important extension is provided by the expectile
regression (Newey & Powell 1987), which can be thought of as a generalization
of the classical mean regression based on an asymmetric squared loss function.
Similar to quantile regression, expectile regression allows to characterize
the entire conditional distribution of a response variable, but possesses
several advantages over the former. First, expectiles are more informative
than quantiles since they rely on tail expectations whereas quantiles use only
the information on whether an observation is below or above the predictor.
Second, the squared loss is continuously differentiable which makes the
estimators and their covariance matrix easier to compute using fast and
efficient algorithms. For these reasons, expectile models have been
implemented in several fields, such as longitudinal data (Tzavidis et al.
2016; Alfò et al. 2017; Barry et al. 2021), spatial analysis (Sobotka & Kneib
2012; Spiegel et al. 2020), life expectancy (Nigri et al. 2022), economics and
finance (Taylor 2008; Kim & Lee 2016; Bellini & Di Bernardino 2017; Bottone et
al. 2021). Especially in the context of risk management, expectiles have
gained an important role as potential competitors to the VaR and the Expected
Shortfall measures. Indeed, they possess several interesting properties in
terms of risk measures (see for instance Bellini 2012, Bellini et al. 2014 and
Ziegel 2016), and are the only risk measure that is both coherent (Artzner et
al. 1999) and elicitable (Lambert et al. 2008).

Moreover, when modeling financial time series, returns often exhibit a
clustering behavior over time which cannot be captured by traditional
homogeneous regression models. Risk managers and regulators are increasingly
interested in determining whether, and how, their temporal evolution can be
influenced by hidden variables, e.g., the state of the market, during tranquil
and crisis periods. In this context, Hidden Markov Models (HMMs, see MacDonald
& Zucchini 1997; Zucchini et al. 2016) have been successfully employed in the
analysis of time series data, with applications to asset allocation and stock
returns as discussed in Mergner & Bulla (2008); De Angelis & Paas (2013);
Nystrup et al. (2017) and Maruotti et al. (2019). Quantile regression methods
have also been generalized to account for serial heterogeneity. For example,
Liu (2016) consider a quantile autoregression in which the parameters are
subject to regime shifts determined by the outcome of a latent, discrete-state
Markov process, while Adam et al. (2019) propose a model-based clustering
approach where groups are inferred from conditional quantiles; see also Ye et
al. (2016), Maruotti et al. (2021) and Merlo et al. (2022) for other
applications of regime-switching models to financial and environmental time
series. In longitudinal data, Farcomeni (2012) and Marino et al. (2018)
introduce linear quantile regression models where time-dependent unobserved
heterogeneity is described through dynamic coefficients that evolve according
to a homogeneous hidden Markov chain. Within a Bayesian framework, a quantile
nonhomogeneous HMM for longitudinal data has been recently proposed by Liu et
al. (2021). To the best of our knowledge, however, a HMM for estimating
conditional expectiles has not yet been proposed in the literature.
Motivated by the advantages of expectiles and the versatility of HMMs, we
develop a linear expectile hidden Markov regression model to analyze the tail
relation between cryptocurrencies and traditional asset classes. The method
introduced allows to examine the entire conditional distribution of returns
given the hidden state and potential covariates, where the dynamics of returns
over time is described by state-specific regression coefficients which follow
a latent discrete homogeneous Markov chain. Inference about model parameters
is carried out in a Maximum Likelihood (ML) approach using an Expectation-
Maximization (EM) algorithm based on the asymmetric normal distribution of
Waldmann et al. (2017) as working likelihood. From a risk management
standpoint, the proposed methodology contributes to identify and control for
potential inherent risks related to the participation in crypto exchanges to
develop appropriate policies and risk assessment procedures.
The study period considered starts from September 2014 until October 2022,
comprising numerous events that heavily impacted financial stability, as the
Chinese stock market crash of 2015, the crypto currency bubble crisis in
2017-2018, the COVID-19 outbreak in 2020 and the Russian invasion of Ukraine
at the beginning of 2022. Following Corbet et al. (2018), we model Bitcoin
daily returns as a function of major stock and global market indices,
including Crude Oil, Standard $\&$ Poor’s 500 (S$\&$P500), Gold COMEX daily
closing prices and the Volatility Index (VIX). Our results show that Bitcoin
returns exhibit a clear temporal clustering behavior in calm and turbulent
periods, and they are strongly associated with traditional assets at low and
high expectile levels.
In concluding, we also evaluate the performance of our approach in a
simulation study, generating observations from a two-state HMM under two
different sample sizes and two different distributions for the error terms.

The rest of the paper is organized as follows. Section 2 briefly reviews the
expectile regression. In Section 3 we specify the proposed model with the EM
algorithm for estimating the model parameters and the computational aspects.
In Section 4, we evaluate the performance of our proposal in a simulation
study. Section 5 shows the empirical analysis and discusses the results
obtained while Section 6 concludes.

##  2 Expectile regression

Expectile regression has been proposed by Newey & Powell (1987) as a
“quantile-like” generalization of standard mean regression based on asymmetric
least-squares estimation. Similarly to quantile regression of Koenker &
Bassett (1978), this is an alternative approach for characterizing the entire
conditional distribution of a response variable where the quantile loss
function is substituted with an asymmetric squared loss function. Formally,
the expectile of order $\tau\in(0,1)$ of a continuous response $Y$ given the
$P$-dimensional vector of covariates $\boldsymbol{X}=\boldsymbol{x}$, is
defined as the minimizer, $\mu_{\boldsymbol{x}}(\tau)$, of the following
problem:

$\mu_{\boldsymbol{x}}(\tau)=\underset{\mu\in\mathbb{R}}{{\rm
argmin}}\,\mathbb{E}[\omega_{\tau}(Y-\mu_{\boldsymbol{x}}(\tau))],$ (1)

where $\omega_{\tau}(u)=u^{2}\lvert\tau-\mathbb{I}(u<0)\rvert$ is the
asymmetric square loss and $\mathbb{I}(\cdot)$ denotes the indicator function.
In a regression framework, for a given $\tau$, a linear expectile model is
defined as
$\mu_{\boldsymbol{x}}(\tau)=\boldsymbol{x}^{\prime}\boldsymbol{\beta}(\tau)$,
where $\boldsymbol{\beta}(\tau)\in\mathbb{R}^{P}$ is the regression parameter
vector. If $\tau=\frac{1}{2}$, expectile regression reduces to the standard
mean regression while for $\tau\neq\frac{1}{2}$ it allows to target the entire
conditional distribution of the response given the covariates similarly to
quantile regression. When we turn from quantiles to expectiles, the latter
possess several advantages over the former. Particularly, we gain uniqueness
of the ML solutions which is, indeed, not granted in the quantile context.
From a computational standpoint, since the squared loss function
$\omega_{\tau}(\cdot)$ is differentiable, the regression parameters
$\boldsymbol{\beta}(\tau)$ can be estimated by efficient Iterative Reweighted
Least Squares (IRLS), in contrast to algorithms used for fitting quantile
regression models. Proofs of consistency, asymptotic normality and a robust
estimator of the variance-covariance matrix of the regression coefficients for
inference have been established in Newey & Powell (1987). These properties
make the expectile regression versatile and computationally appealing from a
statistical point of view.
In a likelihood approach, Gerlach & Chen (2015) and Waldmann et al. (2017)
originally introduced the idea of expectile regression by employing a
likelihood function that is based on the Asymmetric Normal (AN) distribution.
The AN distribution can be thought of as a generalization of the normal
distribution to allow for non-zero skewness, having the following density:

$f_{Y}(y)=\frac{2\sqrt{\tau(1-\tau)}}{\sqrt{\pi\sigma^{2}}(\sqrt{\tau}+\sqrt{1-\tau})}\exp\left[-\omega_{\tau}\left(\frac{y-\mu}{\sigma}\right)\right],$
(2)

where $\mu\in\mathbb{R}$ is a location parameter corresponding to the
$\tau$-th expectile of $Y$, $\sigma>0$ is a scale parameter and $\tau\in(0,1)$
determines the asymmetry of the distribution. Particularly, when
$\tau=\frac{1}{2}$ the density in (2) reduces to the well-known normal
distribution, and $\mu$ and $\sigma$ coincide with its mean and standard
deviation, respectively. As discussed by Waldmann et al. (2017), the
minimization of the asymmetric squared loss function in (1) is equivalent, in
terms of parameter estimates, to the maximization of the likelihood associated
with the AN density.
In the following section, we extend the expectile regression to the HMM
setting by using the AN distribution as working likelihood.

##  3 Methodology

In this section we describe the expectile hidden Markov regression model in
order to take into account the temporal evolution of the time series under
analysis. We then show how inference about model parameters can be carried out
in a ML approach using the AN distribution introduced in the previous section.

Formally, let $\\{S_{t}\\}_{t=1}^{T}$ be a latent, homogeneous, first-order
Markov chain defined on the discrete state space $\\{1,\dots,K\\}$. Let
$\pi_{k}=Pr(S_{1}=k)$ be the initial probability of state $k$, $k=1,\dots,K$,
and $\pi_{k\lvert j}=Pr(S_{t+1}=k\lvert S_{t}=j)$, with
$\sum_{k=1}^{K}\pi_{k\lvert j}=1$ and $\pi_{k\lvert j}\geq 0$, denote the
transition probability between states $j$ and $k$, that is, the probability to
visit state $k$ at time $t+1$ from state $j$ at time $t$, $j,k=1,\dots,K$ and
$t=1,\dots,T$. More concisely, we collect the initial and transition
probabilities in the $K$-dimensional vector $\boldsymbol{\pi}$ and in the
$K\times K$ matrix $\boldsymbol{\Pi}$, respectively.

To build the proposed model, let $Y_{t}$ denote a continuous observable
response variable and $\boldsymbol{X}_{t}=(1,X_{t2},\dots,X_{tP})^{\prime}$ be
a vector of $P$ exogenous covariates, with the first element being the
intercept, at time $t=1,\dots,T$. For a given expectile level $\tau\in(0,1)$,
the proposed linear Expectile Hidden Markov Model (EHMM) is defined as
follows:

$Y_{t}=\boldsymbol{X}^{\prime}_{t}\boldsymbol{\beta}_{k}(\tau)+\epsilon_{tk}(\tau),$
(3)

with
$\boldsymbol{\beta}_{k}(\tau)=(\beta_{1k}(\tau),\dots,\beta_{Pk}(\tau))^{\prime}\in\mathbb{R}^{P}$
being a state-specific coefficient vector that assumes one of the values
$\\{\boldsymbol{\beta}_{1}(\tau),\dots,\boldsymbol{\beta}_{K}(\tau)\\}$
depending on the outcome of the unobservable Markov chain $S_{t}$ and where
$\epsilon_{tk}(\tau)$ is the error term whose conditional $\tau$-th expectile
is assumed to be zero.
Extending the approach of Waldmann et al. (2017) to the HMM setting, we use
the AN distribution to describe the conditional distribution of the response
given covariates and the state occupied by the latent process at time $t$,
whose probability density function is now given by

$f_{Y}(y_{t}\lvert\boldsymbol{X}_{t}=\boldsymbol{x}_{t},S_{t}=k)=\frac{2\sqrt{\tau(1-\tau)}}{\sqrt{\pi\sigma^{2}_{k}}(\sqrt{\tau}+\sqrt{1-\tau})}\exp\left[-\omega_{\tau}\left(\frac{y_{t}-\mu_{tk}}{\sigma_{k}}\right)\right],$
(4)

where the location parameter $\mu_{tk}$ is defined by the linear model
$\mu_{tk}=\boldsymbol{x}^{\prime}_{t}\boldsymbol{\beta}_{k}(\tau)$.
In the following section we use the AN distribution as a working likelihood
for estimating the model parameters in a regression framework.

###  3.1 Likelihood inference

In this section we consider a ML approach to make inference on model
parameters. As is common for HMMs, and for latent variable models in general,
we develop an EM algorithm (Baum et al. 1970) to estimate the parameters of
the method proposed based on the observed data. To ease the notation, unless
specified otherwise, hereinafter we omit the expectile level $\tau$, yet all
model parameters are allowed to depend on it.
For a given number of hidden states $K$, the EM algorithm runs on the complete
log-likelihood function of the model introduced, which is defined as

$\displaystyle\ell_{c}(\boldsymbol{\theta}_{\tau})=\sum_{k=1}^{K}{\gamma_{1}}(k)\log\pi_{k}$
$\displaystyle+\sum_{t=1}^{T}\sum_{k=1}^{K}\sum_{j=1}^{K}{\xi_{t}}(j,k)\log\pi_{k\lvert
j}$ (5) $\displaystyle+\sum_{t=1}^{T}\sum_{k=1}^{K}{\gamma_{t}}(k)\log
f_{Y}(y_{t}\lvert\boldsymbol{x}_{t},S_{t}=k),$

where
$\boldsymbol{\theta}_{\tau}=(\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{K},\sigma_{1},\dots,\sigma_{K},\boldsymbol{\pi},\boldsymbol{\Pi})$
represents the vector of all model parameters, ${\gamma_{t}}(k)$ denotes a
dummy variable equal to $1$ if the latent process is in state $k$ at occasion
$t$ and 0 otherwise, and ${\xi_{t}}(j,k)$ is a dummy variable equal to $1$ if
the process is in state $j$ in $t-1$ and in state $k$ at time $t$ and $0$
otherwise.

To estimate $\boldsymbol{\theta}_{\tau}$, the algorithm iterates between two
steps, the E- and M-steps, until convergence, as outlined below.

#### E-step:

In the E-step, at the generic $(h+1)$-th iteration, the unobservable indicator
variables ${\gamma_{t}}(k)$ and ${\xi_{t}}(j,k)$ in (5) are replaced by their
conditional expectations given the observed data and the current parameter
estimates $\boldsymbol{\theta}_{\tau}^{(h)}$. To compute such quantities we
require the calculation of the probability of being in state $k$ at time $t$
given the observed sequence

$\gamma_{t}^{(h)}(k)=P_{\boldsymbol{\theta}_{\tau}^{(h)}}(S_{t}=k\lvert
y_{1},\dots,y_{T})$ (6)

and the probability that at time $t-1$ the process is in state $j$ and then in
state $k$ at time $t$, given the observed sequence

$\xi_{t}^{(h)}(j,k)=P_{\boldsymbol{\theta}_{\tau}^{(h)}}(S_{t-1}=j,S_{t}=k\lvert
y_{1},\dots,y_{T}).$ (7)

The quantities in (6) and (7) can be obtained using the Forward-Backward
algorithm of Welch (2003). Then, we use these to calculate the conditional
expectation of the complete log-likelihood function in (5) given the observed
data and the current estimates:

$\displaystyle
Q(\boldsymbol{\theta}_{\tau}\lvert\boldsymbol{\theta}_{\tau}^{(h)})=\sum_{k=1}^{K}\gamma^{(h)}_{1}(k)\log\pi_{k}$
$\displaystyle+\sum_{t=1}^{T}\sum_{k=1}^{K}\sum_{j=1}^{K}\xi^{(h)}_{t}(j,k)\log\pi_{k\lvert
j}$ (8) $\displaystyle+\sum_{t=1}^{T}\sum_{k=1}^{K}\gamma^{(h)}_{t}(k)\log
f_{Y}(y_{t}\lvert\boldsymbol{x}_{t},S_{t}=k).$

#### M-step:

In the M-step we maximize
$Q(\boldsymbol{\theta}_{\tau}\lvert\boldsymbol{\theta}_{\tau}^{(h)})$ in (8)
with respect to $\boldsymbol{\theta}_{\tau}$ to obtain the update parameter
estimates $\boldsymbol{\theta}_{\tau}^{(h+1)}$. The maximization of
$Q(\boldsymbol{\theta}_{\tau}\lvert\boldsymbol{\theta}_{\tau}^{(h)})$ can be
partitioned into orthogonal subproblems, where the updating formulas for the
hidden Markov chain and state-dependent regression parameters are obtained
independently maximizing each of these terms. Formally, the initial
probabilities $\pi_{k}$ and transition probabilities $\pi_{k\lvert j}$ are
updated using:

$\pi^{(h+1)}_{k}=\gamma^{(h)}_{1}(k),\quad k=1,\dots,K$ (9)

and

$\pi^{(h+1)}_{k\lvert
j}=\frac{\sum_{t=1}^{T}\xi^{(h)}_{t}(j,k)}{\sum_{t=1}^{T}\sum_{k=1}^{K}\xi^{(h)}_{t}(j,k)},\quad
j,k=1,\dots,K.$ (10)

To update the regression coefficients, the first-order condition of (8) with
respect to $\boldsymbol{\beta}_{k}$, $k=1,\dots,K$, yields

$\frac{\partial
Q(\boldsymbol{\theta}_{\tau}\lvert\boldsymbol{\theta}_{\tau}^{(h)})}{\partial\boldsymbol{\beta}_{k}}\propto\sum_{t=1}^{T}\gamma^{(h)}_{t}(k)\lvert\tau-\mathbb{I}(y_{t}<\boldsymbol{x}^{\prime}_{t}\boldsymbol{\beta}_{k})\lvert\boldsymbol{x}_{t}(y_{t}-\boldsymbol{x}^{\prime}_{t}\boldsymbol{\beta}_{k})=\boldsymbol{0}_{P},$
(11)

so the M-step update expression for $\boldsymbol{\beta}_{k}$ is

$\boldsymbol{\beta}^{(h+1)}_{k}=\Big{(}\sum_{t=1}^{T}\gamma^{(h)}_{t}(k)\lvert\tau-\mathbb{I}(y_{t}<\boldsymbol{x}^{\prime}_{t}\boldsymbol{\beta}_{k})\lvert\boldsymbol{x}_{t}\boldsymbol{x}^{\prime}_{t}\Big{)}^{-1}\Big{(}\sum_{t=1}^{T}\gamma^{(h)}_{t}(k)\lvert\tau-\mathbb{I}(y_{t}<\boldsymbol{x}^{\prime}_{t}\boldsymbol{\beta}_{k})\lvert\boldsymbol{x}_{t}y_{t}\Big{)},$
(12)

which can be computed using IRLS for cross-sectional data with appropriate
weights. Similarly, from the first-order condition of (8) with respect to the
scale parameters we obtain the following M-step update formula for
$\sigma^{2}_{k}$:

$\sigma^{2}_{k}{}^{(h+1)}=\frac{2}{\sum_{t=1}^{T}\gamma^{(h)}_{t}(k)}\sum_{t=1}^{T}\gamma^{(h)}_{t}(k)\lvert\tau-\mathbb{I}(y_{t}<\boldsymbol{x}^{\prime}_{t}\boldsymbol{\beta}^{(h+1)}_{k})\lvert(y_{t}-\boldsymbol{x}^{\prime}_{t}\boldsymbol{\beta}^{(h+1)}_{k})^{2}.$
(13)

The E- and M- steps are alternated until convergence, that is when the
observed likelihood between two consecutive iterations is smaller than a
predetermined threshold. In this paper, we set this threshold criterion equal
to $10^{-4}$.
Following Maruotti et al. (2021) and Merlo et al. (2022), for fixed $\tau$ and
$K$ we initialize the EM algorithm by providing the initial states partition,
$\\{S_{t}^{(0)}\\}_{t=1}^{T}$, according to a Multinomial distribution with
probabilities $1/K$. From the generated partition, the elements of
$\boldsymbol{\Pi}^{(0)}$ are computed as proportions of transition, while we
obtain $\boldsymbol{\beta}^{(0)}_{k}$ and $\sigma^{(0)}_{k}$ by fitting mean
regressions on the observations within state $k$. To deal with the possibility
of multiple roots of the likelihood equation and better explore the parameter
space, we fit the proposed EHMM using a multiple random starts strategy with
different starting partitions and retain the solution corresponding to the
maximum likelihood value.
Once we computed the ML estimate of the model parameters, to estimate the
standard errors we employ the parametric bootstrap scheme of Visser et al.
(2000). In practice, we refit the model to $R$ bootstrap samples and
approximate the standard error of each model parameter with the corresponding
standard deviation of the bootstrap estimates.

##  4 Simulation study

We conduct a simulation study to validate the performance of our method under
different scenarios in terms of recovering the true values of the parameters
and the clustering performance. Similar to the work of Maruotti et al. (2021),
we analyze two different sample sizes ($T=500,\ T=1000$) and two different
distributions for the error term. For each scenario we conduct $500$ Monte
Carlo simulations. We draw observations from a two state HMM ($K=2$) using the
following data generating process:

$Y_{t}=\begin{cases}-1+2X_{t}+\epsilon_{t1},\ \ S_{t}=1\\\
1-2X_{t}+\epsilon_{t2},\ \ S_{t}=2,\end{cases}$ (14)

with $X_{t}\sim\mathcal{N}(0,1)$. We consider two distributions for the error
terms in (14). In the first scenario, $\epsilon_{tk}$ is generated from a
normal distribution with standard deviation 1, for $k=1,2$. In the second one,
$\epsilon_{tk}$ is generated from a skew-$t$ distribution with 5 degrees of
freedom and asymmetry parameter 2, for $k=1,2$. Finally, the matrix of
transition probabilities is set equal to
$\boldsymbol{\Pi}=\bigl{(}\begin{smallmatrix}0.8&0.2\\\
0.1&0.9\end{smallmatrix}\bigr{)}$.

In order to assess the validity of the model we fit the proposed EHMM at five
expectile levels, i.e., $\tau=\\{0.10,0.25,0.50,0.75,0.90\\}$, and compute the
bias and standard errors associated to the state-specific coefficients,
averaged over the Monte Carlo replications, for each combination of sample
size and error distribution. Tables 1 and 2 report the simulation outputs for
the normal and skew-$t$ distributions, respectively. As can be observed, the
precision of the estimates is higher at the center of the distribution rather
than on the tails, mainly due to the reduced number of observations at extreme
expectile levels, but the bias always remains under control. Evidently, in
Table 2 a higher standard deviation shows up for the skew-$t$ distribution due
to the asymmetry and heavier tails than the normal density, but both the bias
and the standard deviation tend to decrease as the sample size increases.
Concerning the hidden process, given the true values of the transition
probabilities in $\boldsymbol{\Pi}$, we see that the coefficients
corresponding to the first state are estimated with lower precision because
fewer transitions occur from one state to the other, as expected.
To evaluate the ability in recovering the true states partition we consider
the Adjusted Rand Index (ARI) of Hubert & Arabie (1985). The state partition
provided by the fitted models is obtained by taking the maximum,
$\underset{k}{\max}\,\gamma_{t}(k)$, posteriori probability for every
$t=1,\dots,T$, and report the box-plots of ARI for the posterior probabilities
in Figure 1 for the four settings considered. Firstly, we observe that the
distribution error plays a fundamental role in estimating the true states
partition as, for all five expectile levels, we obtain a better clustering
performance for the model with Gaussian errors with respect to the skew-$t$
case. Secondly, the goodness of the clustering obtained partially depends on
the specific expectile level, being the values slightly higher at the mean
($\tau=0.50$) than at the tails. Finally, when increasing the sample size to
$T=1000$, results slightly improve reporting a lower variability for both
error distributions. Overall, the proposed EHMM is able to recover the true
values of the parameters and the true state partition highly satisfactory in
all the cases examined.

$\tau$ | 0.10 | 0.25 | 0.50 | 0.75 | 0.90
---|---|---|---|---|---
| Bias | Std.Err. | Bias | Std.Err. | Bias | Std.Err. | Bias | Std.Err. | Bias | Std.Err.
Panel A: $T=500$
State 1 |  |  |  |  |  |  |  |  |  |
$\beta_{11}$ = -1 | 0.020 | 0.093 | 0.010 | 0.076 | -0.002 | 0.073 | -0.018 | 0.080 | -0.048 | 0.100
$\beta_{21}$ = 2 | 0.001∗ | 0.109 | 0.001 | 0.093 | 0.004 | 0.087 | 0.010 | 0.089 | 0.023 | 0.101
State 2 |  |  |  |  |  |  |  |  |  |
$\beta_{12}$ = 1 | 0.040 | 0.055 | 0.013 | 0.041 | -0.002 | 0.037 | -0.013 | 0.039 | -0.027 | 0.047
$\beta_{22}$ = -2 | -0.008 | 0.067 | 0.001∗ | 0.058 | 0.001 | 0.055 | -0.003 | 0.058 | -0.012 | 0.068
Panel B: $T=1000$
State 1 |  |  |  |  |  |  |  |  |  |
$\beta_{11}$ = -1 | 0.021 | 0.068 | 0.010 | 0.055 | -0.001 | 0.051 | -0.016 | 0.056 | -0.042 | 0.070
$\beta_{21}$ = 2 | 0.006 | 0.071 | 0.003 | 0.060 | 0.004 | 0.057 | 0.008 | 0.060 | 0.017 | 0.069
State 2 |  |  |  |  |  |  |  |  |  |
$\beta_{12}$ = 1 | 0.039 | 0.038 | 0.014 | 0.029 | 0.001∗ | 0.026 | -0.010 | 0.028 | -0.023 | 0.034
$\beta_{22}$ = -2 | -0.012 | 0.053 | -0.004 | 0.045 | -0.002 | 0.043 | -0.005 | 0.044 | -0.014 | 0.050

Table 1: Bias and standard error values of the state-regression parameter
estimates with Gaussian distributed errors for $T=500$ (Panel A) and $T=1000$
(Panel B). ∗ represents values smaller (in absolute value) than 0.001.

$\tau$ | 0.10 | 0.25 | 0.50 | 0.75 | 0.90
---|---|---|---|---|---
| Bias | Std.Err. | Bias | Std.Err. | Bias | Std.Err. | Bias | Std.Err. | Bias | Std.Err.
Panel A: $T=500$
State 1 |  |  |  |  |  |  |  |  |  |
$\beta_{11}$ = -1 | -0.064 | 0.228 | -0.009 | 0.119 | -0.004 | 0.099 | -0.018 | 0.166 | 0.083 | 0.489
$\beta_{21}$ = 2 | -0.153 | 0.319 | -0.049 | 0.183 | 0.001 | 0.125 | 0.027 | 0.161 | -0.045 | 0.471
State 2 |  |  |  |  |  |  |  |  |  |
$\beta_{12}$ = 1 | 0.181 | 0.144 | 0.066 | 0.073 | 0.013 | 0.051 | -0.018 | 0.060 | -0.052 | 0.131
$\beta_{22}$ = -2 | -0.060 | 0.100 | -0.029 | 0.080 | -0.015 | 0.074 | -0.015 | 0.079 | -0.038 | 0.103
Panel B: $T=1000$
State 1 |  |  |  |  |  |  |  |  |  |
$\beta_{11}$ = -1 | -0.053 | 0.162 | 0.001 | 0.077 | -0.004 | 0.069 | -0.022 | 0.123 | 0.025 | 0.369
$\beta_{21}$ = 2 | -0.133 | 0.236 | -0.024 | 0.111 | 0.013 | 0.082 | 0.039 | 0.116 | 0.027 | 0.319
State 2 |  |  |  |  |  |  |  |  |  |
$\beta_{12}$ = 1 | 0.168 | 0.101 | 0.057 | 0.048 | 0.010 | 0.035 | -0.019 | 0.041 | -0.059 | 0.083
$\beta_{22}$ = -2 | -0.067 | 0.067 | -0.033 | 0.052 | -0.016 | 0.048 | -0.015 | 0.054 | -0.031 | 0.072

Table 2: Bias and standard error values of the state-regression parameter
estimates with skew-$t$ distributed errors for $T=500$ (Panel A) and $T=1000$
(Panel B).

Figure 1: From left to right, box-plots of ARI for the posterior probabilities
for Gaussian (red) and skew-$t$ (blue) distributed errors with $T=500$ and
$T=1000$.

##  5 Empirical application

In this section we apply the methodology proposed to analyze the Bitcoin daily
returns as a function of global leading financial indices. Over the last
decade, cryptocurrencies and in particular the Bitcoin market played a leading
role, attracting attentions of researchers and investors. Their peculiar
characteristics, such their extreme price volatility, driven by market
speculation and technology applications, often lead to price bubbles, euphoria
and market instability. In order to address these periods of upheaval, it is
crucial to understand the association between Bitcoin and globally relevant
market indices in such circumstances of financial turmoil. Consistently with
Corbet et al. (2018), here we consider the Bitcoin, Crude Oil, Standard $\&$
Poor’s 500 (S$\&$P500), Gold COMEX daily closing prices and the Volatility
Index (VIX) from September 2014 to October 2022. All series are expressed in
US dollars and have been downloaded from the Yahoo finance database. Daily
returns with continuous compounding are calculated taking the logarithm of the
difference between closing prices in consecutive trading days and then
multiplied by 100, i.e., $r_{t}=100\cdot\log(P_{t}/P_{t-1})$ where $P_{t}$ is
the closing price on day $t$, for a total of $T=2025$ observations.
The considered timespan is marked by numerous crises that may have impacted
cross-market integration patterns, including the Chinese stock market crash of
2015, the cryptocurrency bubble crisis in 2017-2018, the COVID-19 pandemic and
the Russian invasion of Ukraine at the beginning of 2022, which have caused
unprecedented levels of uncertainty and risk. In Table 3 we report the list of
examined variables and the summary statistics for the whole sample. First
thing to notice is that Bitcoin is generally much more volatile than the other
assets, having the highest standard deviation. The Bitcoin returns also show
very high negative skewness and very high kurtosis, as well as S$\&$P500. The
highest level of kurtosis is reported by Crude Oil, probably determined by the
prices’ fluctuations after the COVID-19 outbreak. On the other side, the large
positive skewness of VIX indicates longer and fatter tails on the right side
of the distribution, highlighting the well-known inverse relationship with the
S$\&$P500. In concluding, the Augmented Dickey-Fuller (ADF) test Dickey &
Fuller (1979) shows that all daily returns are stationary at the $1\%$ level
of significance. Following these considerations, and motivated by the reforms
considered by markets authorities to protect investors and preserve stability
in response to cryptocurrencies’ downturns, the proposed EHMM can provide
insights into the temporal evolution of Bitcoin returns and describe how this
is affected by rapid changes in markets volatility.

| Minimum | Mean | Maximum | Std.Err. | Skewness | Kurtosis | Jarque-Bera test | ADF test
---|---|---|---|---|---|---|---|---
Bitcoin | -46.4730 | 0.1859 | 22.5119 | 4.6165 | -0.6817 | 8.7172 | 6568.469 | -11.117
Crude Oil | -28.2206 | 0.0280 | 31.9634 | 3.1077 | 0.0942 | 21.1219 | 37645.730 | -10.324
S&P500 | -12.7652 | 0.0316 | 8.9683 | 1.1716 | -0.9033 | 16.3473 | 22823.380 | -12.390
Gold | -5.1069 | 0.0152 | 5.7775 | 0.9358 | -0.0698 | 4.1741 | 1471.712 | -12.292
VIX | -29.9831 | 0.0379 | 76.8245 | 8.3414 | 1.2683 | 6.6648 | 4290.737 | -14.209

Table 3: Descriptive statistics for the whole sample. The Jarque-Bera test and
the ADF test statistics are displayed in boldface when the null hypothesis is
rejected at the $1\%$ significance level.

To this end, we consider the following linear EHMM:

$\mu_{tk}^{Bitcoin}=\beta_{1k}(\tau)+\beta_{2k}(\tau)r_{t}^{Crude\,\,Oil}+\beta_{3k}(\tau)r_{t}^{S\&P500}+\beta_{4k}(\tau)r_{t}^{Gold}+\beta_{5k}(\tau)r_{t}^{VIX},$
(15)

with $\mu_{tk}^{Bitcoin}$ corresponding to the $\tau$-th conditional expectile
of Bitcoin return at time $t$ in state $k$, while $r_{t}^{Crude\,\,Oil}$
denotes the return of the same date for Crude Oil, and similarly for the other
indices.
As a first step of the empirical analysis, in order to select the number of
latent states we fit the proposed EHMM for different values of $K$ varying
from 2 to 5 at three expectile levels $\tau=\\{0.10,0.50,0.90\\}$, which allow
us to focus on both downside and upside risks. To compare models with
differing number of states, Table 4 reports three widely employed penalized
likelihood selection criteria for $K$, namely the AIC (Akaike 1998), the BIC
(Schwarz et al. 1978) and the ICL (Biernacki et al. 2000). As one can see, the
AIC selects 5, or more, states, while BIC chooses $K=4$ for all three
expectile levels. This should not be surprising since the AIC tends to
overestimate the true number of hidden states. On the contrary, ICL favors a
more parsimonious choice as $K=2$ is always considered to be the optimal
number of states at $\tau=0.10$, $\tau=0.50$ and $\tau=0.90$. For these
reasons, and in order to clearly identify high and low volatility market
conditions we thus consider the proposed EHMM with $K=2$ states for all three
$\tau$ levels.
For the selected models, we report the clustering results in Figure 2 at the
investigated expectile levels. Each plot shows the time series of Bitcoin
daily returns colored according to the estimated posterior probability of
class membership, $\underset{k}{\max}\,\gamma_{t}(k)$, with the vertical
dashed lines representing globally relevant events such as the Chinese stock
market crash in 2015, the cryptocurrencies crash at the beginning of 2018, the
COVID-19 market crash in March 2020, Biden’s election at the USA presidency in
November 2020 and the Russian invasion of Ukraine at the beginning of 2022.
Here we clearly see that the latent components can be associated to specific
market regimes characterized by low and high volatility periods. Specifically,
light-blue points (State 1) tend to identify low returns, while dark-blue ones
(State 2) correspond to periods of extreme positive and negative returns.
Moving on to the state-specific model parameters, Table 5 shows the parameter
estimates along with the standard errors (in brackets) computed by using the
parametric bootstrap technique illustrated in Section 3.1 over $R=1000$
resamples. First, consistently with the quantile regression literature, the
intercepts are increasing with $\tau$, with State 1 having lower values than
State 2 for all $\tau$s. Second, it is interesting to observe that in the not-
at-risk state (State 1) the S$\&$P500, Gold and the VIX index positively
influence extreme left-tail ($\tau=0.10$) movements of Bitcoin returns, while
only S$\&$P500 and Gold significantly influence the right-tail ($\tau=0.90$)
expectiles of the cryptocurrency, exposing a connection during high volatility
periods between traditional financial markets and Bitcoin both for negative
and positive returns. At $\tau=0.50$, instead, Bitcoin can be considered as a
weak hedge during high volatility periods since it is not statistically
associated with all the assets considered (Bouri, Jalkh, Molnár & Roubaud,
2017; Bouri, Molnár, Azzi, Roubaud & Hagfors, 2017). In the at-risk state
(State 2) we observe a positive influence of the S$\&$P500 and Gold across the
conditional distribution of returns. Also, one can see that Crude Oil is
negatively associated with the crypto returns at the $10$-th expectile. This
finding is in line with Bouri et al. (2020a) but it is contrary to the works
of Dyhrberg (2016) and Corbet et al. (2018), which may be due to the events
and crises occurred in the last years.
Finally, the estimated state-dependent scale parameter $\sigma_{1}$ reflects
more stable periods for the first state, meanwhile $\sigma_{2}$ contemplates
rapid (positive and negative) peak and burst returns for the second state,
confirming the graphical analysis conducted in Figure 2.

Figure 2: From top to bottom, Bitcoin returns series classified according to
the estimated posterior probability of class membership at $\tau=0.10$,
$\tau=0.50$ and $\tau=0.90$. Vertical dashed lines indicate globally relevant
events in the financial markets that occurred in 2015,06; 2017,12; 2020,03;
2020,11; and 2022,02.

| $\tau=0.10$ | $\tau=0.50$ | $\tau=0.90$
---|---|---|---
AIC |  |  |
$K=2$ | 11347.4122 | 11231.4286 | 11389.7954
$K=3$ | 11210.6074 | 11126.1327 | 11204.4727
$K=4$ | 11109.3930 | 11051.9556 | 11115.4655
$K=5$ | 11055.7980 | 11014.8696 | 11079.7912
BIC |  |  |
$K=2$ | 11431.6121 | 11315.6285 | 11473.9953
$K=3$ | 11356.5539 | 11272.0791 | 11350.4192
$K=4$ | 11328.3127 | 11270.8753 | 11334.3852
$K=5$ | 11358.9175 | 11317.9892 | 11382.9108
ICL |  |  |
$K=2$ | 12784.4362 | 12787.4471 | 12926.5182
$K=3$ | 13490.5599 | 13406.8221 | 13616.8609
$K=4$ | 13880.2346 | 13487.8461 | 13759.7332
$K=5$ | 14111.8066 | 13511.4108 | 13908.7362

Table 4: AIC, BIC and ICL values with varying number of states for the
investigated expectile levels. Bold font highlights the best values for the
considered criteria (lower-is-better).

| Intercept | Crude Oil | S&P500 | Gold | VIX | $\sigma_{k}$
---|---|---|---|---|---|---
State 1 |  |  |  |  |  |
$\tau=0.10$ | -1.036 (0.280) | 0.024 (0.021) | 0.595 (0.096) | 0.189 (0.072) | 0.029 (0.012) | 1.433 (0.040)
$\tau=0.50$ | 0.122 (0.158) | 0.031 (0.072) | 0.409 (0.383) | 0.263 (0.249) | 0.009 (0.036) | 1.695 (0.062)
$\tau=0.90$ | 1.297 (0.061) | -0.009 (0.020) | 0.589 (0.088) | 0.134 (0.065) | 0.014 (0.011) | 1.335 (0.041)
State 2 |  |  |  |  |  |
$\tau=0.10$ | -6.52 (0.060) | -0.256 (0.096) | 2.072 (0.476) | 1.032 (0.320) | -0.055 (0.058) | 4.964 (0.157)
$\tau=0.50$ | 0.242 (0.092) | -0.056 (0.055) | 1.087 (0.357) | 0.613 (0.214) | -0.025 (0.026) | 6.164 (0.169)
$\tau=0.90$ | 6.244 (0.229) | 0.017 (0.079) | 0.948 (0.291) | 0.835 (0.249) | -0.002 (0.041) | 4.692 (0.128)

Table 5: State-specific parameter estimates for three expectile levels, with
bootstrapped standard errors (in brackets) obtained over 1000 replications.
Point estimates are displayed in boldface when significant at the standard 5%
level.

##  6 Conclusions

The increasing popularity and importance of Bitcoin in the financial landscape
have made scholars and practitioners interrogated about its properties and its
relation to other assets. In this regard, we contribute to the existing
literature in two ways. From a theoretical standpoint, we develop a linear
expectile hidden Markov model for the analysis of time series where temporal
behaviors of the data are captured via time-dependent coefficients that follow
an unobservable discrete homogeneous Markov chain. The proposed method enables
us to model the entire conditional distribution of asset returns and, at the
same time, to grasp unobserved serial heterogeneity and rapid volatility jumps
that would otherwise go undetected. From a practical point of view, we analyze
the association between Bitcoin and a collection of global market indices, not
only at the average, but also during times of market stress.
Empirically, we find evidence of strong and positive interrelations among
Bitcoin returns and S$\&$P500 and Gold, and, at the same time, we observe the
capacity of the Bitcoin of working as a weak hedge during high volatility
periods, contributing to the existing strands of literature on the subject
(Baur et al., 2018; Corbet et al., 2018, 2019). Its partial capacity of being
a weak hedge but not a safe haven it is consistent with the excess volatility
of Bitcoin and indications that assets with no history as a safe haven are
unlikely to be considered “safe” in an economic or financial crisis (Baur et
al., 2018).
As a possible next step, our methodology could be extended to the hidden semi-
Markov model setting where the sojourn-distributions, that is, the
distributions of the number of consecutive time points that the chain spends
in each state, can be modeled by the researcher using either parametric or
non-parametric approaches instead of assuming geometric sojourn densities as
in HMMs.

## References

  * (1)
  * Adam et al. (2019) Adam, T., Langrock, R. & Kneib, T. (2019), Model-based clustering of time series data: a flexible approach using nonparametric state-switching quantile regression models, in ‘Proceedings of the 12th Scientific Meeting on Classification and Data Analysis’, pp. 8–11.
  * Akaike (1998) Akaike, H. (1998), Information theory and an extension of the maximum likelihood principle, in ‘Selected papers of Hirotugu Akaike’, Springer, pp. 199–213.
  * Alfò et al. (2017) Alfò, M., Salvati, N. & Ranallli, M. G. (2017), ‘Finite mixtures of quantile and M-quantile regression models’, Statistics and Computing 27(2), 547–570.
  * Artzner et al. (1999) Artzner, P., Delbaen, F., Eber, J.-M. & Heath, D. (1999), ‘Coherent measures of risk’, Math. Finance 9(3), 203–228.
  * Barry et al. (2021) Barry, A., Oualkacha, K. & Charpentier, A. (2021), ‘A new GEE method to account for heteroscedasticity using asymmetric least-square regressions’, Journal of Applied Statistics pp. 1–27.
  * Baum et al. (1970) Baum, L. E., Petrie, T., Soules, G. & Weiss, N. (1970), ‘A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains’, The Annals of Mathematical Statistics 41(1), 164–171.
  * Baur et al. (2018) Baur, D. G., Hong, K. & Lee, A. D. (2018), ‘Bitcoin: Medium of exchange or speculative assets?’, Journal of International Financial Markets, Institutions and Money 54, 177–189.
  * Bellini (2012) Bellini, F. (2012), ‘Isotonicity properties of generalized quantiles’, Statistics & Probability Letters 82(11), 2017–2024.
  * Bellini & Di Bernardino (2017) Bellini, F. & Di Bernardino, E. (2017), ‘Risk management with expectiles’, The European Journal of Finance 23(6), 487–506.
  * Bellini et al. (2014) Bellini, F., Klar, B., Müller, A. & Rosazza Gianin, E. (2014), ‘Generalized quantiles as risk measures’, Insurance Math. Econom. 54, 41–48.
  * Biernacki et al. (2000) Biernacki, C., Celeux, G. & Govaert, G. (2000), ‘Assessing a mixture model for clustering with the integrated completed likelihood’, IEEE Transactions on Pattern Analysis and Machine Intelligence 22(7), 719–725.
  * Bottone et al. (2021) Bottone, M., Petrella, L. & Bernardi, M. (2021), ‘Unified Bayesian conditional autoregressive risk measures using the skew exponential power distribution’, Statistical Methods & Applications 30(3), 1079–1107.
  * Bouri, Jalkh, Molnár & Roubaud (2017) Bouri, E., Jalkh, N., Molnár, P. & Roubaud, D. (2017), ‘Bitcoin for energy commodities before and after the december 2013 crash: diversifier, hedge or safe haven?’, Applied Economics 49(50), 5063–5073.
  * Bouri et al. (2020a) Bouri, E., Lucey, B. & Roubaud, D. (2020a), ‘Cryptocurrencies and the downside risk in equity investments’, Finance Research Letters 33, 101211.
  * Bouri et al. (2020b) Bouri, E., Lucey, B. & Roubaud, D. (2020b), ‘The volatility surprise of leading cryptocurrencies: Transitory and permanent linkages’, Finance Research Letters 33, 101188.
  * Bouri, Molnár, Azzi, Roubaud & Hagfors (2017) Bouri, E., Molnár, P., Azzi, G., Roubaud, D. & Hagfors, L. I. (2017), ‘On the hedge and safe haven properties of Bitcoin: Is it really more than a diversifier?’, Finance Research Letters 20, 192–198.
  * Cheah & Fry (2015) Cheah, E.-T. & Fry, J. (2015), ‘Speculative bubbles in Bitcoin markets? An empirical investigation into the fundamental value of Bitcoin’, Economics Letters 130, 32–36.
  * Corbet et al. (2019) Corbet, S., Lucey, B., Urquhart, A. & Yarovaya, L. (2019), ‘Cryptocurrencies as a financial asset: A systematic analysis’, International Review of Financial Analysis 62, 182–199.
  * Corbet et al. (2018) Corbet, S., Meegan, A., Larkin, C., Lucey, B. & Yarovaya, L. (2018), ‘Exploring the dynamic relationships between cryptocurrencies and other financial assets’, Economics Letters 165, 28–34.
  * De Angelis & Paas (2013) De Angelis, L. & Paas, L. J. (2013), ‘A dynamic analysis of stock markets using a hidden Markov model’, Journal of Applied Statistics 40(8), 1682–1700.
  * Dickey & Fuller (1979) Dickey, D. A. & Fuller, W. A. (1979), ‘Distribution of the estimators for autoregressive time series with a unit root’, Journal of the American Statistical Association 74(366a), 427–431.
  * Dyhrberg (2016) Dyhrberg, A. H. (2016), ‘Hedging capabilities of Bitcoin. Is it the virtual gold?’, Finance Research Letters 16, 139–144.
  * Engle & Manganelli (2004) Engle, R. F. & Manganelli, S. (2004), ‘CAViaR: conditional autoregressive Value at Risk by regression quantiles’, Journal of Business & Economic Statistics 22(4), 367–381.
  * Farcomeni (2012) Farcomeni, A. (2012), ‘Quantile regression for longitudinal data based on latent markov subject-specific parameters’, Statistics and Computing 22(1), 141–152.
  * Gerlach & Chen (2015) Gerlach, R. & Chen, C. W. (2015), ‘Bayesian expected shortfall forecasting incorporating the intraday range’, Journal of Financial Econometrics 14(1), 128–158.
  * Guesmi et al. (2019) Guesmi, K., Saadi, S., Abid, I. & Ftiti, Z. (2019), ‘Portfolio diversification with virtual currency: Evidence from Bitcoin’, International Review of Financial Analysis 63, 431–437.
  * Hubert & Arabie (1985) Hubert, L. & Arabie, P. (1985), ‘Comparing partitions’, Journal of Classification 2(1), 193–218.
  * Ji et al. (2018) Ji, Q., Bouri, E., Gupta, R. & Roubaud, D. (2018), ‘Network causality structures among bitcoin and other financial assets: A directed acyclic graph approach’, The Quarterly Review of Economics and Finance 70, 203–213.
  * Katsiampa et al. (2019) Katsiampa, P., Corbet, S. & Lucey, B. (2019), ‘Volatility spillover effects in leading cryptocurrencies: A BEKK-MGARCH analysis’, Finance Research Letters 29, 68–74.
  * Kim & Lee (2016) Kim, M. & Lee, S. (2016), ‘Nonlinear expectile regression with application to Value-at-Risk and Expected Shortfall estimation’, Computational Statistics & Data Analysis 94, 1–19.
  * Koenker & Bassett (1978) Koenker, R. & Bassett, G. (1978), ‘Regression quantiles’, Econometrica: Journal of the Econometric Society 46(1), 33–50.
  * Kristjanpoller et al. (2020) Kristjanpoller, W., Bouri, E. & Takaishi, T. (2020), ‘Cryptocurrencies and equity funds: Evidence from an asymmetric multifractal analysis’, Physica A: Statistical Mechanics and Its Applications 545, 123711.
  * Lambert et al. (2008) Lambert, N. S., Pennock, D. M. & Shoham, Y. (2008), Eliciting properties of probability distributions, in ‘Proceedings of the 9th ACM Conference on Electronic Commerce’, ACM, pp. 129–138.
  * Liu et al. (2021) Liu, H., Song, X., Tang, Y. & Zhang, B. (2021), ‘Bayesian quantile nonhomogeneous hidden Markov models’, Statistical Methods in Medical Research 30(1), 112–128.
  * Liu (2016) Liu, X. (2016), ‘Markov switching quantile autoregression’, Statistica Neerlandica 70(4), 356–395.
  * MacDonald & Zucchini (1997) MacDonald, I. L. & Zucchini, W. (1997), Hidden Markov and other models for discrete-valued time series, Vol. 110, CRC Press.
  * Marino et al. (2018) Marino, M. F., Tzavidis, N. & Alfò, M. (2018), ‘Mixed hidden Markov quantile regression models for longitudinal data with possibly incomplete sequences’, Statistical Methods in Medical Research 27(7), 2231–2246.
  * Maruotti et al. (2021) Maruotti, A., Petrella, L. & Sposito, L. (2021), ‘Hidden semi-Markov-switching quantile regression for time series’, Computational Statistics & Data Analysis 159, 107208.
  * Maruotti et al. (2019) Maruotti, A., Punzo, A. & Bagnato, L. (2019), ‘Hidden Markov and semi-Markov models with multivariate leptokurtic-normal components for robust modeling of daily returns series’, Journal of Financial Econometrics 17(1), 91–117.
  * Mergner & Bulla (2008) Mergner, S. & Bulla, J. (2008), ‘Time-varying beta risk of Pan-European industry portfolios: A comparison of alternative modeling techniques’, The European Journal of Finance 14(8), 771–802.
  * Merlo et al. (2022) Merlo, L., Maruotti, A., Petrella, L. & Punzo, A. (2022), ‘Quantile hidden semi-Markov models for multivariate time series’, Statistics and Computing 32(4), 1–22.
  * Merlo et al. (2021) Merlo, L., Petrella, L. & Raponi, V. (2021), ‘Forecasting VaR and ES using a joint quantile regression and its implications in portfolio allocation’, Journal of Banking & Finance p. 106248.
  * Naeem et al. (2021) Naeem, M. A., Bouri, E., Peng, Z., Shahzad, S. J. H. & Vo, X. V. (2021), ‘Asymmetric efficiency of cryptocurrencies during COVID19’, Physica A: Statistical Mechanics and its Applications 565, 125562.
  * Nakamoto (2008) Nakamoto, S. (2008), ‘Bitcoin: A peer-to-peer electronic cash system’, Decentralized Business Review p. 21260.
  * Newey & Powell (1987) Newey, W. K. & Powell, J. L. (1987), ‘Asymmetric least squares estimation and testing’, Econometrica: Journal of the Econometric Society pp. 819–847.
  * Nigri et al. (2022) Nigri, A., Barbi, E. & Levantesi, S. (2022), ‘The relationship between longevity and lifespan variation’, Statistical Methods & Applications 31(3), 481–493.
  * Nystrup et al. (2017) Nystrup, P., Madsen, H. & Lindström, E. (2017), ‘Long memory of financial time series and hidden Markov models with time-varying parameters’, Journal of Forecasting 36(8), 989–1002.
  * Schwarz et al. (1978) Schwarz, G. et al. (1978), ‘Estimating the dimension of a model’, The Annals of Statistics 6(2), 461–464.
  * Shahzad et al. (2022) Shahzad, S. J. H., Bouri, E., Ahmad, T. & Naeem, M. A. (2022), ‘Extreme tail network analysis of cryptocurrencies and trading strategies’, Finance Research Letters 44, 102106.
  * Sobotka & Kneib (2012) Sobotka, F. & Kneib, T. (2012), ‘Geoadditive expectile regression’, Computational Statistics & Data Analysis 56(4), 755–767.
  * Spiegel et al. (2020) Spiegel, E., Kneib, T. & Otto-Sobotka, F. (2020), ‘Spatio-temporal expectile regression models’, Statistical Modelling 20(4), 386–409.
  * Taylor (2008) Taylor, J. W. (2008), ‘Estimating Value at Risk and Expected Shortfall using expectiles’, Journal of Financial Econometrics 6(2), 231–252.
  * Taylor (2019) Taylor, J. W. (2019), ‘Forecasting Value at Risk and Expected Shortfall using a semiparametric approach based on the asymmetric Laplace distribution’, Journal of Business & Economic Statistics 37(1), 121–133.
  * Tzavidis et al. (2016) Tzavidis, N., Salvati, N., Schmid, T., Flouri, E. & Midouhas, E. (2016), ‘Longitudinal analysis of the strengths and difficulties questionnaire scores of the Millennium Cohort Study children in England using M-quantile random-effects regression’, Journal of the Royal Statistical Society: Series A (Statistics in Society) 179(2), 427–452.
  * Visser et al. (2000) Visser, I., Raijmakers, M. E. & Molenaar, P. C. (2000), ‘Confidence intervals for hidden Markov model parameters’, British Journal of Mathematical and Statistical Psychology 53(2), 317–327.
  * Waldmann et al. (2017) Waldmann, E., Sobotka, F. & Kneib, T. (2017), ‘Bayesian regularisation in geoadditive expectile regression’, Statistics and Computing 27(6), 1539–1553.
  * Welch (2003) Welch, L. R. (2003), ‘Hidden Markov models and the Baum-Welch algorithm’, IEEE Information Theory Society Newsletter 53(4), 10–13.
  * White et al. (2015) White, H., Kim, T.-H. & Manganelli, S. (2015), ‘VAR for VaR: Measuring tail dependence using multivariate regression quantiles’, Journal of Econometrics 187(1), 169–188.
  * Yarovaya et al. (2016) Yarovaya, L., Brzeszczyński, J. & Lau, C. K. M. (2016), ‘Intra-and inter-regional return and volatility spillovers across emerging and developed markets: Evidence from stock indices and stock index futures’, International Review of Financial Analysis 43, 96–114.
  * Ye et al. (2016) Ye, W., Zhu, Y., Wu, Y. & Miao, B. (2016), ‘Markov regime-switching quantile regression models and financial contagion detection’, Insurance: Mathematics and Economics 67, 21–26.
  * Yi et al. (2018) Yi, S., Xu, Z. & Wang, G.-J. (2018), ‘Volatility connectedness in the cryptocurrency market: Is Bitcoin a dominant cryptocurrency?’, International Review of Financial Analysis 60, 98–114.
  * Zhang et al. (2021) Zhang, Y.-J., Bouri, E., Gupta, R. & Ma, S.-J. (2021), ‘Risk spillover between bitcoin and conventional financial markets: An expectile-based approach’, The North American Journal of Economics and Finance 55, 101296.
  * Ziegel (2016) Ziegel, J. F. (2016), ‘Coherence and elicitability’, Mathematical Finance 26(4), 901–918.
  * Zucchini et al. (2016) Zucchini, W., MacDonald, I. L. & Langrock, R. (2016), Hidden Markov models for time series: an introduction using R, Chapman and Hall/CRC.