††thanks: A.M. and Y.Z. contributed equally to this work.††thanks: A.M. and
Y.Z. contributed equally to this work.

# Topogivity: A Machine-Learned Chemical Rule for Discovering Topological
Materials

Andrew Ma  Department of Electrical Engineering and Computer Science,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Yang Zhang  Department of Physics, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA Thomas Christensen  Department of Physics,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA Hoi
Chun Po  Department of Physics, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA Department of Physics, Hong Kong
University of Science and Technology, Clear Water Bay, Hong Kong, China Li
Jing  Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA Facebook AI Research, New York, New York 10003, USA
Liang Fu<EMAIL_ADDRESS>Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA Marin Soljačić
<EMAIL_ADDRESS>Department of Physics, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA

###### Abstract

Topological materials present unconventional electronic properties that make
them attractive for both basic science and next-generation technological
applications. The majority of currently-known topological materials have been
discovered using methods that involve symmetry-based analysis of the quantum
wavefunction. Here we use machine learning to develop a simple-to-use
heuristic chemical rule that diagnoses with a high accuracy whether a material
is topological using only its chemical formula. This heuristic rule is based
on a notion that we term _topogivity_ , a machine-learned numerical value for
each element that loosely captures its tendency to form topological materials.
We next implement a high-throughput strategy for discovering topological
materials based on the heuristic topogivity-rule prediction followed by ab
initio validation. This way, we discover new topological materials that are
not diagnosable using symmetry indicators, including several that may be
promising for experimental observation.

Topological materials – including both topological insulators [1, 2, 3, 4, 5]
and topological semimetals [6, 7, 8, 9] – are unconventional phases of matter
characterized by topologically nontrivial electron wavefunctions. Since the
beginning of the field, an important and enduring question has been how to
determine if a given electronic material is topological. Efforts to answer
this question have largely relied on first-principles calculations in synergy
with topological band theory [10, 11]. In particular, recently-developed
theories known as symmetry indicators [12] and topological quantum chemistry
[13] enable the diagnosis of a wide range of topological materials using
symmetry-based analysis of the wavefunction [14, 15, 16]. These symmetry-based
methods require relatively low computational cost and have enabled high-
throughput computational searches for topological materials [17, 18, 19].

Despite these successes, symmetry indicators do not apply to materials having
low symmetry crystal structures [12], and therefore cannot diagnose all
topological materials. For example, the first experimentally-observed Weyl
semimetal, tantalum arsenide (TaAs) [7, 20], is a non-symmetry-diagnosable
topological material [17]. Its topological nature is established by
calculating wavefunction-based topological invariant directly [21], which
involves significant computational cost. From a broad conceptual standpoint,
topological materials such as Chern insulators and time-reversal-invariant
$Z_{2}$ topological insulators are robust against any small perturbation
breaking all crystal symmetries, which renders symmetry indicators
inapplicable even though topology remains intact. Thus, for both practical and
fundamental reasons, it is highly desirable to develop accurate and simple-to-
use rules to determine whether any given material is topological.

Figure 1: Topogivity-based diagnosis and discovery of topological materials. a
Given a stoichiometric material, the topogivity-based heuristic diagnosis is
evaluated by simply weighting the material’s elements’ (A, B, and C)
topogivities ($\tau_{\text{A}}$, $\tau_{\text{B}}$, and $\tau_{\text{C}}$) by
their relative abundance in the chemical formula (AB and AC3). The sign
indicates the topological classification and the magnitude indicates roughly
how confident we are in this classification. Each element’s topogivity is a
machine-learned parameter that loosely captures the element’s tendency to form
topologically nontrivial materials. b We leverage our framework to perform
high-throughput topological materials discovery. First, we use the
topogivities to rapidly screen through a suitable collection of materials
(i.e. the discovery space) in order to find candidate topological materials.
Subsequently, we carry out ab initio validation by performing DFT on these
candidates. We discovered topological materials that are not diagnosable using
the standard symmetry indicators approach [12].

Many aspects of quantum solids can be understood at a heuristic level from a
chemistry perspective. A well-known example is bonding, which can be
understood using quantum-mechanical approaches such as molecular orbital
theory [22], as well as using heuristics such as the difference of element
electronegativities. While quantum theory can provide greater detail and
accuracy, chemical heuristics can often provide valuable insight and a useful
guide for materials discovery. Is there a deep chemical reason why a
particular material is topological? To what extent can topological materials
be understood and identified using chemical heuristic approaches? While
connections between chemistry and electronic band topology (e.g., based on the
presence or absence of certain elements) have been explored [23, 24, 25, 26,
27], existing chemical heuristics do not provide a broadly applicable path for
finding topological materials.

Here, we use machine learning (ML) to help us search for chemical origins of
topological electronic structure in materials. Recently, ML has become a
powerful approach for advancing scientific discovery in materials science [28,
29, 30]. In the area of topological materials, researchers have begun to apply
ML to both toy models [31, 32, 33, 34] and ab initio data [35, 36, 37, 38, 39,
40]. While ML has led to important advances for many applications in
engineering and science, most ML models act essentially like black boxes: they
are complicated models which often provide correct answers, but because of
their complexity, they are difficult to understand, and hence provide little
insight and intuition about the systems they are applied to. Since a key aim
of our work is precisely to learn insights about topology, we instead focus
our quest onto interpretable ML models, with a goal of finding a broadly
applicable heuristic chemical rule that diagnoses whether or not a material is
topological.

The heuristic rule that our ML approach discovered is based on the notion of a
learned parameter for each element that loosely captures the tendency of an
element to form topological materials. We refer to this as an element’s
_topogivity_. This heuristic rule is simple, hand-calculable, and
interpretable: a given material is diagnosed with high accuracy (typically
$>80\%$) as topologically nontrivial (trivial) if the weighted average of its
elements’ topogivities is positive (negative) (Fig. 1a). The heuristic rule
does not rely on crystal symmetry and our approach can be used to make
predictions on all materials. We integrate the heuristic rule into a high-
throughput procedure to search for non-symmetry-diagnosable topological
materials, in which we perform screening using the heuristic rule followed by
density functional theory (DFT) validation (Fig. 1b). The newly-discovered
topological materials include several high-quality examples that may be
promising for experimental realization.

## Classes of Materials

Conventional textbook chemistry teaches that the electrons of insulators
(including semiconductors) are localized to ionic or covalent bonds, while the
electrons of metals are delocalized and “free”. From a band-theoretic
perspective, the former make up part of the class of materials known as atomic
insulators and the latter roughly correspond to “ordinary” metals. Topological
insulators and topological semimetals do not fit into this conventional
dichotomy. Topological insulators feature a band gap and a nontrivial
topological invariant, and as a consequence, their electronic states cannot be
reduced to an assembly of localized atomic or molecular orbitals. Topological
semimetals have band degeneracies protected by symmetry or topology near or at
the Fermi level. Collectively, we refer to topological insulators and
topological semimetals as topological materials, and refer to all other
materials as trivial materials.

To learn a heuristic chemical rule for diagnosing topological materials, we
employ a supervised learning approach. This requires a labeled dataset in
which each material is labeled as “topological” or “trivial”. The learned
heuristic chemical rule is then applied to screen another dataset which we
refer to as the discovery space. Existing ab initio databases of
stoichiometric, non-magnetic, three-dimensional materials [18, 19, 41, 17]
offer a convenient source of data. However, it is important to note that they
are imperfect, in part because the symmetry-based high-throughput screening
methods that were used to generate them are inherently incapable of detecting
certain topology. Taking such limitations into account, we identify the
labeled dataset as a subset of the database generated by Tang  _et al._ [17]
(see Supplementary Section S1) (our methodology could also be applied to other
databases that contain both trivial materials and topological materials). Our
labeled dataset consists of 9,026 materials, of which $51\%$ are labeled as
trivial and the remaining $49\%$ are labeled as topological. However, due to
the aforementioned imperfection, for ML purposes this labeled dataset should
effectively be considered as a dataset with noisy labels, e.g., some
topological materials are incorrectly labeled as trivial. Separately, the
discovery space consists of 1,433 materials, whose topology cannot be
determined from the symmetry indicators method (see Supplementary Section S1).
By applying the learned heuristic chemical rule to the discovery space and
then performing DFT, we are able to evaluate its ability to predict
topological materials beyond those diagnosable by existing standard
approaches. Some of the topological materials that we identify in the
discovery space are known elsewhere in the literature and serve primarily as
confirmations, whereas others represent instances of truly new materials
discovery.

## Learning a heuristic chemical rule

Figure 2: Periodic table of topogivities. Machine-learned topogivities
$\tau_{E}$ are shown by color-coding and in values. Elements that do not
appear in any material in the dataset are shown in gray. Example applications
of the topogivity-based heuristic chemical rule are shown for two materials,
NaCl (trivial) and Na3Bi (nontrivial) [6].

Our machine learning model takes the form of a heuristic chemical rule.
Specifically, the model maps each material $M$ to a number $g(M)$ according to
the function

$g(M)=\sum_{E}f_{E}(M)\tau_{E},$

where the summation runs over the elements present in the chemical formula of
material $M$, $\tau_{E}$ is a learned parameter for each element $E$, and
$f_{E}(M)$ is the element fraction for the element $E$ in material $M$ (e.g.,
for a chemical formula $A_{x}B_{y}C_{z}$, $f_{A}(M)=\frac{x}{x+y+z}$,
$f_{B}(M)=\frac{y}{x+y+z}$, and $f_{C}(M)=\frac{z}{x+y+z}$). Classification
decisions are made according to the sign of $g(M)$: classify as topological if
positive and classify as trivial if negative. A greater magnitude of $g(M)$
roughly corresponds to a more confident classification decision. We refer to
the model as a heuristic _chemical_ rule in the sense that all the information
required for obtaining a diagnosis is contained in the material’s chemical
formula. Additional modeling and methodological details are provided in
Supplementary Sections S2.A and S2.B.

For each element $E$, we refer to the optimized parameter $\tau_{E}$ as its
_topogivity_. For a given material $M$, $g(M)$ is simply the weighted average
of its elements’ topogivities, where the weighting is with respect to each
element’s relative abundance, as identifiable from the material’s chemical
formula. Conceptually, an element’s topogivity loosely captures its tendency
to form topological materials – greater topogivity is intended to roughly
correspond to a greater tendency (see Supplementary Section S2.F for details
and caveats on the meaning of topogivity).

Before making predictions in the discovery space, we want to first evaluate
model performance within the labeled dataset. To do this, we use a nested
cross validation procedure and average the results over multiple test sets. We
find an average of $82.7\%$ accuracy. Additionally, we find empirical evidence
that as the magnitude of $g(M)$ is increased, the fraction of correctly-
classified materials first increases and then plateaus, with the plateau
beginning around $|g(M)|\approx 1$. We heuristically set a threshold of 1.0
for a high-confidence topologically nontrivial classification and observe on
average that $93.0\%$ of materials with $g(M)\geq 1.0$ are correctly
classified. After we are done doing nested cross validation, we then proceed
to use the entire labeled dataset to fit the final model, which is what we
will use for making predictions in the discovery space. Details and extended
results are presented in Supplementary Sections S2.C, S2.D, and S2.E.

We visualize the final model’s learned topogivities in Fig. 2. This periodic
table of topogivities enables an immediate heuristic diagnosis of any
stoichiometric material whose elements are featured in the table. This is
illustrated with examples in Fig. 2 for the trivial insulator NaCl and the
Dirac semimetal Na3Bi [6]. The Weyl semimetal TaAs [7] is also worth
highlighting: TaAs is non-symmetry-diagnosable [17] and does not appear in the
labeled dataset, but is successfully diagnosed as topological by the
topogivity approach:
$g(\text{TaAs})=\tfrac{1}{2}\tau_{\text{Ta}}+\tfrac{1}{2}\tau_{\text{As}}=1.450$.

Figure 3: Selection of newly-discovered topological materials. These materials
are not diagnosable using symmetry indicators [12], but were successfully
discovered using our topogivity-based approach. The band structures were
computed using DFT. MoPt2Si3 (a) and Ge6Y4Zn5 (b) are nonsymmorphic Dirac
semimetals. Pb3Pd5 (c) and Bi3Pd8 (d) are Kramers Weyl semimetals. For each
material, relevant topological degeneracies are highlighted by circles.

The simplicity of our model enables us to readily extract chemical insights
from the periodic table of topogivities. First, we observe that elements that
are near each other in the periodic table tend to have similar topogivities,
which is consistent with intuition. Second, we observe that the elements with
negative topogivities are located in two clusters respectively in the top
right and bottom left parts of the periodic table. This is also consistent
with intuition, since ionic compounds often have large trivial band gaps and
elements from these two clusters tend to form ionic compounds. Third,
considering group 15 (the pnictogens), we observe that while N, P, and As have
negative topogivities (and Sb has a small positive topogivity), Bi has a
positive topogivity with a relatively large magnitude. This is consistent with
the intuition that Bi often plays a role in topological materials [25].
Finally, we observe a region of high topogivities in the early transition
metals – future work could attempt to understand the underlying reasons for
this (note that there is a chance that typical oxidation states are
artificially inflating these topogivities). Overall, while the element
topogivities are parameters whose specific learned values are affected by
dataset and modeling limitations (see Supplementary Section S2.F for
discussion), the fact that we can extract chemical insights that are
consistent with intuition is evidence that a topogivity-based approximate
picture can provide a meaningful way to study topological materials.

## High-throughput screening and ab initio calculations

To identify topological materials using the learned topogivities, we compute
$g(M)$ for each of the 1,433 materials in the discovery space. We restrict our
attention to the materials that have a $g(M)$ value that corresponds to a
topologically nontrivial classification with high-confidence (i.e., $g(M)\geq
1.0$): that leaves 73 materials (after the removal of 2 other materials; see
Supplementary Section S3.A). Additionally, since it is difficult to obtain
accurate DFT calculations for f-electron materials, we exclude any material
that contains a 4f or 5f electron, eliminating 5 materials and thus leaving us
with 68 materials for ab initio validation.

For each of these 68 materials, we perform DFT within the generalized gradient
approximation (GGA) [42]. We include spin-orbital coupling in all of our
calculations, which is consistent with the database [17] (note that
topological classification can depend on the presence or absence of spin-
orbital coupling [12]). In our DFT calculations, we checked for many – but not
all – types of topological materials. In principle it is possible that some
topological materials were not detected by our DFT. Further methodological
details are provided in Supplementary Section S3.B.

Of the 68 materials, we find 56 topological materials, corresponding to a
success rate of 82.4%. The topological materials consist of 48 Weyl
semimetals, 7 Dirac semimetals, and 1 Dirac nodal line semimetal
(Supplementary Section S4). We note that there are aspects of our procedure
and data analysis that could have introduced some bias into the success rate
(see Supplementary Section S3.C). Some of these 56 topological materials have
previously been predicted in the literature and a smaller portion have also
already been experimentally observed, e.g., TaAs [7]. More importantly, our
DFT calculations also identify multiple new topological materials that to our
knowledge have not been previously identified.

The band structures of four particularly interesting newly-discovered
topological materials are shown in Fig. 3. Each is a topological semimetal
with a relatively clean band structure and at least one band crossing within
50 meV of the Fermi level, making it promising for potential experimental
investigation. MoPt2Si3 and Ge6Y4Zn5 are both nonsymmorphic Dirac semimetals.
At the Z point, the former has a Dirac point in the valence band manifold and
the latter has Dirac points in both the valence and conduction band manifolds.
Pb3Pd5 and Bi3Pd8 are both Kramers Weyl semimetals [43]. The former has Weyl
nodes at the L and Z points, and the latter has two Weyl nodes close in energy
at the L point. In particular, we highlight that MoPt2Si3 has a Dirac point
close to the Fermi level as well as a relatively clean Fermi surface, and
Pb3Pd5 has a Weyl node located at the $Z$ point that is right at the Fermi
level. We emphasize that the reason the band degeneracies in these four
materials are non-symmetry-diagnosable is that they are all within the valence
band manifold or conduction band manifold. Such band degeneracies cannot be
diagnosed by the symmetry indicators method, which is formulated based on the
electron filling and therefore cannot target band degeneracies that are not
_between_ the valence and conduction bands [12]. The band structures of all of
the topogivity-identified topological materials are included in the
Supplementary Material.

## Discussion and Outlook

It is important to note that topogivity is not an unambiguously defined
quantity, as its exact numerical value for each element can depend, for
example, on the choice of machine learning algorithm and the use of the
weighted average formulation (Supplementary Section S2.F). Moreover, the
topogivity approach provides only a coarse-grained topological classification
– nontrivial or trivial – without the fine-grained detail of ab initio
approaches.

Nevertheless, topogivity offers a broadly applicable and simple approach for
diagnosing topological materials. This diagnosis uses only the chemical
formula and requires merely a handful of arithmetic operations to evaluate.
Furthermore, the periodic table of topogivities (Fig. 2) provides simple
intuition for a complex, exotic phenomenon. An important highlight of the
topogivity approach is that it enables the discovery of non-symmetry-
diagnosable topological materials.

One worthy future direction is to look for a more complete understanding of
the underlying reasons for the values of the elements’ topogivities, which may
in turn shed new light on the fundamental question of why some materials are
topological while others are not. Another promising path forward would be to
perform more comprehensive searches for new topological materials using
topogivity-based strategies. Finally, it is intriguing to contemplate whether
our interpretable-ML approach, used here to discover topogivity, could perhaps
be used for other material properties as well, such as ferroelectricity,
ferromagnetism, or maybe even superconductivity.

## Acknowledgements

We thank Pawan Goyal for assisting preliminary work. We thank Samuel Kim,
Peter Lu, Rumen Dangovski, and Edward Zhang for helpful discussions. We thank
Feng Tang and Xiangang Wan for the sharing of materials data. We thank Paola
Rebusco for critical reading and editing of the manuscript. A.M. acknowledges
support from the MIT EECS Alan L. McWhorter Fellowship and from the National
Science Foundation Graduate Research Fellowship under Grant No. 1745302.
Research was sponsored in part by the United States Air Force Research
Laboratory and the United States Air Force Artificial Intelligence Accelerator
and was accomplished under Cooperative Agreement Number FA8750-19-2-1000. The
views and conclusions contained in this document are those of the authors and
should not be interpreted as representing the official policies, either
expressed or implied, of the United States Air Force or the U.S. Government.
The U.S. Government is authorized to reproduce and distribute reprints for
Government purposes notwithstanding any copyright notation herein. This
material is also based upon work supported in part by the Air Force Office of
Scientific Research under the awards number FA9550-20-1-0115,
FA9550-21-1-0299, and FA9550-21-1-0317, as well as in part by the US Office of
Naval Research (ONR) Multidisciplinary University Research Initiative (MURI)
grant N00014-20-1-2325 on Robust Photonic Materials with High-Order
Topological Protection. This material is also based upon work supported in
part by the U. S. Army Research Office through the Institute for Soldier
Nanotechnologies at MIT, under Collaborative Agreement Number
W911NF-18-2-0048. This work is also supported in part by the the National
Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI
Institute for Artificial Intelligence and Fundamental Interactions,
http://iaifi.org/). The work of Y.Z. and L.F. was supported by DOE Office of
Basic Energy Sciences, Division of Materials Sciences and Engineering under
Award DE-SC0018945 (theoretical analysis) and DE-SC0019275 (band structure
calculation). L.F. was partly supported by the Simons Investigator Award from
Simons Foundation. The work of H.C.P. was partly supported by a Pappalardo
Fellowship at MIT.

## Code and data availability

The code and data underlying our machine learning model will be made available
in a public repository.

## References

  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007).
  * Hsieh  _et al._ [2008] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nature 452, 970 (2008).
  * Fu [2011] L. Fu, Phys. Rev. Lett. 106, 106802 (2011).
  * Hsieh  _et al._ [2012] T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Nat. Commun. 3, 1 (2012).
  * Benalcazar  _et al._ [2017] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017).
  * Liu  _et al._ [2014] Z. Liu, B. Zhou, Y. Zhang, Z. Wang, H. Weng, D. Prabhakaran, S.-K. Mo, Z. Shen, Z. Fang, X. Dai, _et al._ , Science 343, 864 (2014).
  * Xu  _et al._ [2015] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, _et al._ , Science 349, 613 (2015).
  * Burkov  _et al._ [2011] A. Burkov, M. Hook, and L. Balents, Phys. Rev. B 84, 235126 (2011).
  * Bradlyn  _et al._ [2016] B. Bradlyn, J. Cano, Z. Wang, M. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353, 558 (2016).
  * Xiao and Yan [2021] J. Xiao and B. Yan, Nat. Rev. Phys. 3, 283 (2021).
  * Bansil  _et al._ [2016] A. Bansil, H. Lin, and T. Das, Rev. Mod. Phys. 88, 021004 (2016).
  * Po  _et al._ [2017] H. C. Po, A. Vishwanath, and H. Watanabe, Nat. Commun. 8, 1 (2017).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. Vergniory, Z. Wang, C. Felser, M. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Fu and Kane [2007] L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).
  * Kruthoff  _et al._ [2017] J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Slager, Phys. Rev. X 7, 041069 (2017).
  * Song  _et al._ [2018] Z. Song, T. Zhang, Z. Fang, and C. Fang, Nat. Commun. 9, 1 (2018).
  * Tang  _et al._ [2019] F. Tang, H. C. Po, A. Vishwanath, and X. Wan, Nature 566, 486 (2019).
  * Zhang  _et al._ [2019] T. Zhang, Y. Jiang, Z. Song, H. Huang, Y. He, Z. Fang, H. Weng, and C. Fang, Nature 566, 475 (2019).
  * Vergniory  _et al._ [2019] M. Vergniory, L. Elcoro, C. Felser, N. Regnault, B. A. Bernevig, and Z. Wang, Nature 566, 480 (2019).
  * Lv  _et al._ [2015] B. Lv, H. Weng, B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. Huang, L. Zhao, G. Chen, _et al._ , Phys. Rev. X 5, 031013 (2015).
  * Weng  _et al._ [2015] H. Weng, C. Fang, Z. Fang, B. A. Bernevig, and X. Dai, Phys. Rev. X 5, 011029 (2015).
  * Hehre [1976] W. J. Hehre, Acc. Chem. Res. 9, 399 (1976).
  * Kumar  _et al._ [2021] N. Kumar, S. N. Guin, K. Manna, C. Shekhar, and C. Felser, Chem. Rev. 121, 2780 (2021).
  * Gibson  _et al._ [2015] Q. Gibson, L. M. Schoop, L. Muechler, L. Xie, M. Hirschberger, N. P. Ong, R. Car, and R. J. Cava, Phys. Rev. B 91, 205128 (2015).
  * Isaeva and Ruck [2020] A. Isaeva and M. Ruck, Inorg. Chem. 59, 3437 (2020).
  * Klemenz  _et al._ [2020] S. Klemenz, A. K. Hay, S. M. Teicher, A. Topp, J. Cano, and L. M. Schoop, J. Am. Chem. Soc. 142, 6350 (2020).
  * Gui  _et al._ [2019] X. Gui, I. Pletikosic, H. Cao, H.-J. Tien, X. Xu, R. Zhong, G. Wang, T.-R. Chang, S. Jia, T. Valla, _et al._ , ACS Cent. Sci. 5, 900 (2019).
  * Zhong  _et al._ [2020] M. Zhong, K. Tran, Y. Min, C. Wang, Z. Wang, C.-T. Dinh, P. De Luna, Z. Yu, A. S. Rasouli, P. Brodersen, _et al._ , Nature 581, 178 (2020).
  * Deringer  _et al._ [2021] V. L. Deringer, N. Bernstein, G. Csányi, C. B. Mahmoud, M. Ceriotti, M. Wilson, D. A. Drabold, and S. R. Elliott, Nature 589, 59 (2021).
  * George and Hautier [2020] J. George and G. Hautier, Trends Chem. 3, 86 (2020).
  * Zhang and Kim [2017] Y. Zhang and E.-A. Kim, Phys. Rev. Lett. 118, 216401 (2017).
  * Zhang  _et al._ [2018] P. Zhang, H. Shen, and H. Zhai, Phys. Rev. Lett. 120, 066401 (2018).
  * Scheurer and Slager [2020] M. S. Scheurer and R.-J. Slager, Phys. Rev. Lett. 124, 226401 (2020).
  * Zhang  _et al._ [2020] Y. Zhang, P. Ginsparg, and E.-A. Kim, Phys. Rev. Research 2, 023283 (2020).
  * Acosta  _et al._ [2018] C. M. Acosta, R. Ouyang, A. Fazzio, M. Scheffler, L. M. Ghiringhelli, and C. Carbogno, arXiv:1805.10950 (2018).
  * Claussen  _et al._ [2020] N. Claussen, B. A. Bernevig, and N. Regnault, Phys. Rev. B 101, 245117 (2020).
  * Andrejevic  _et al._ [2020] N. Andrejevic, J. Andrejevic, C. H. Rycroft, and M. Li, arXiv:2003.00994 (2020).
  * Cao  _et al._ [2020] G. Cao, R. Ouyang, L. M. Ghiringhelli, M. Scheffler, H. Liu, C. Carbogno, and Z. Zhang, Phys. Rev. Mater. 4, 034204 (2020).
  * Liu  _et al._ [2020] H. Liu, G. Cao, Z. Zhou, and J. Liu, J. Phys.: Condens. Matter 33, 325501 (2020).
  * Schleder  _et al._ [2021] G. R. Schleder, B. Focassio, and A. Fazzio, Appl. Phys. Rev. 8, 031409 (2021).
  * Vergniory  _et al._ [2021] M. G. Vergniory, B. J. Wieder, L. Elcoro, S. S. Parkin, C. Felser, B. A. Bernevig, and N. Regnault, arXiv:2105.09954 (2021).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Chang  _et al._ [2018] G. Chang, B. J. Wieder, F. Schindler, D. S. Sanchez, I. Belopolski, S.-M. Huang, B. Singh, D. Wu, T.-R. Chang, T. Neupert, _et al._ , Nat. Mater. 17, 978 (2018).

in 1,…,37 See pages \x of supp/supp.pdf


# Kernel-based Optimally Weighted Conformal Prediction Intervals

Jonghyeok Lee  H. Milton Stewart School of Industrial and Systems Engineering,
Georgia Institute of Technology. Chen Xu  H. Milton Stewart School of
Industrial and Systems Engineering, Georgia Institute of Technology. Yao Xie
Correspondence<EMAIL_ADDRESS>H. Milton Stewart School of Industrial
and Systems Engineering, Georgia Institute of Technology.

###### Abstract

Conformal prediction has been a popular distribution-free framework for
uncertainty quantification. In this work, we present a novel conformal
prediction method for time-series, which we call Kernel-based Optimally
Weighted Conformal Prediction Intervals (KOWCPI). Specifically, KOWCPI adapts
the classic Reweighted Nadaraya-Watson (RNW) estimator for quantile regression
on dependent data and learns optimal data-adaptive weights. Theoretically, we
tackle the challenge of establishing a conditional coverage guarantee for non-
exchangeable data under strong mixing conditions on the non-conformity scores.
We demonstrate the superior performance of KOWCPI on real time-series against
state-of-the-art methods, where KOWCPI achieves narrower confidence intervals
without losing coverage.

##  1 Introduction

Conformal prediction, originated in Vovk et al. (1999, 2005), offers a robust
framework explicitly designed for reliable and distribution-free uncertainty
quantification. Conformal prediction has become increasingly recognized and
adopted within the domains of machine learning and statistics (Lei et al.,
2013; Lei and Wasserman, 2014; Kim et al., 2020; Angelopoulos and Bates,
2023). Assuming nothing beyond the exchangeability of data, conformal
prediction excels in generating valid prediction sets under any given
significance level, irrespective of the underlying data distribution and model
assumptions. This capability makes it particularly valuable for uncertainty
quantification in settings characterized by diverse and complex models.

Going beyond the exchangeability assumption has been a research challenge,
particularly as many real-world datasets (such as time-series data) are
inherently non-exchangeable. Tibshirani et al. (2019) addresses situations
where a feature distribution shifts between training and test data and
restores valid coverage through weighted quantiles based on the likelihood
ratio of the distributions. More recently, Barber et al. (2023) bounds the
coverage gap using the total variation distance between training and test data
eights and minimizes this gap using pre-specified data-independent weights.
However, it remains open to how to appropriately optimize the weights.

To advance conformal prediction for time series, we extend the prior
sequential predictive approach (Xu and Xie, 2023a, b) by incorporating
nonparametric kernel regression into the quantile regression method on non-
conformity scores. A key challenge of adapting this method to time-series data
lies in selecting optimal weights to accommodate the dependent structure of
the data. To ensure valid coverage of prediction sets, it is crucial to select
weights inside the quantile estimator so that it closely approximates the true
quantile of non-conformity scores.

In this paper, we introduce KOWCPI, which utilizes the Reweighted Nadaraya-
Watson estimator (Hall et al., 1999) to facilitate the selection of data-
dependant optimal weights. This approach anticipates that adaptive weights
will enhance the robustness of uncertainty quantification, particularly when
the assumption of exchangeability is compromised. Our method also addresses
the weight selection issue in the weighted quantile method presented by Barber
et al. (2023), as KOWCPI allows for the calculation of weights in a data-
driven manner without prior knowledge about the data.

In summary, our main contributions are:

  * •

We propose KOWCPI, a sequential time-series conformal prediction method that
performs nonparametric kernel quantile regression on non-conformity scores. In
particular, KOWCPI learns optimal data-driven weights used in the conditional
quantiles.

  * •

We prove the asymptotic conditional coverage guarantee of KOWCPI based on the
classical theory of nonparametric regression. We further obtain the marginal
coverage gap of KOWCPI using the general result for the weights on quantile
for non-exchangeable data.

  * •

We demonstrate the effectiveness of KOWCPI on real time-series data against
state-of-the-art baselines. Specifically, KOWCPI can achieve the narrowest
width of prediction intervals without losing marginal and approximate
conditional (i.e., rolling) coverage empirically.

###  1.1 Literature

#### RNW quantile regression

In Hall et al. (1999), the Reweighted Nadaraya-Watson (RNW, often referred to
as Weighted or Adjusted Nadaraya-Watson) estimator was suggested as a method
to estimate the conditional distribution function from time-series data. This
estimator extends the renowned Nadaraya-Watson estimator (Nadaraya, 1964;
Watson, 1964) by introducing an additional adjustment to the weights, thus
combining the favorable bias properties of the local linear estimator with the
benefit of being a distribution function by itself like the original Nadaraya-
Watson estimator (Hall et al., 1999; Yu and Jones, 1998). The theory of the
regression quantile with the RNW estimator has been further developed by Cai
(2002). Furthermore, Cai (2002) and Salha (2006) demonstrated that the RNW
estimator is consistent under strongly mixing conditions, which are commonly
observed in time-series data. In this work, we adaptively utilize the RNW
estimator within the conformal prediction framework to construct sequential
prediction intervals for time-series data, leveraging its data-driven weights
for quantile estimation and the weighted conformal approach.

#### Conformal prediction with weighted quantiles

Approaches using quantile regression instead of empirical quantiles in
conformal prediction have been widespread (Romano et al., 2019; Kivaranovic et
al., 2020; Gibbs et al., 2023). These methods utilize various quantile
regression techniques to construct conformal prediction intervals, and the
convergence to the oracle prediction width can be shown under the consistency
of the quantile regression function (Sesia and Candès, 2020). Another recent
work by Guan (2023) uses kernel weighting based on the distance between the
test point and data to perform localized conformal prediction, which further
discusses the selection of kernels and bandwidths. Recent work in this
direction of utilizing the weighted quantiles, including Lee et al. (2023);
Angelopoulos et al. (2023), continues to be vibrant. As we will discuss later,
our approach leverages techniques in classical non-parametric statistics when
constructing the weights.

#### Time-series conformal prediction

There is a growing body of research on time-series conformal prediction (Xu
and Xie, 2021b; Gibbs and Candès, 2021). Various applications include
financial markets (Gibbs and Candès, 2021), anomaly detection (Xu and Xie,
2021a), and geological classification (Xu and Xie, 2022). In particular, Gibbs
and Candès (2021, 2022) sequentially construct prediction intervals by
updating the significance level $\alpha$ based on the mis-coverage rate. This
approach has become a major methodology for handling non-exchangeable data,
leading to several subsequent developments (Feldman et al., 2022; Auer et al.,
2023; Bhatnagar et al., 2023; Zaffran et al., 2022). On the other hand, Xu and
Xie (2023b); Xu et al. (2024) take a slightly different approach by conducting
sequential quantile regression using non-conformity scores. Our study aims to
integrate non-parametric kernel estimation for sequential quantile regression,
addressing the weight selection issues identified by Barber et al. (2023).
Additionally, our research aligns with Guan (2023), particularly in utilizing
a dissimilarity measure between the test point and the past data.

##  2 Problem Setup

We begin by assuming that the observations of the random sequence
$(X_{t},Y_{t})\in\mathbb{R}^{d}\times\mathbb{R}$, $t=1,2,\ldots$ are obtained
sequentially. Notably, $X_{t}$ may represent exogenous variables that aid in
predicting $Y_{t}$, the historical values of $Y_{t}$ itself, or a combination
of both. (In Appendix A, we expand our discussion to include cases where the
response $Y_{t}$ is multivariate.) A key aspect of our setup is that the data
are non-exchangeable and exhibit dependencies, which are typical in time-
series data where temporal or sequential dependencies influence predictive
dynamics.

Suppose we are given a pre-specified point predictor
$\hat{f}:\mathbb{R}^{d}\to\mathbb{R}$ trained on a separate dataset or on past
data. This predictor $\hat{f}$ maps a feature variable in $\mathbb{R}^{d}$ to
a scalar point prediction for $Y_{t}$. Given a user-specified significance
level $\alpha\in(0,1)$, we use the initial $T$ observations to construct
prediction intervals $\widehat{C}_{t-1}^{\alpha}(X_{t})$ for $Y_{t}$ in a
sequential manner from $t=T+1$ onwards.

Two key types of coverage targeted by prediction intervals are marginal
coverage and conditional coverage. Marginal coverage is defined as

$\mathbb{P}(Y_{t}\in\widehat{C}_{t-1}^{\alpha}(X_{t}))\geq 1-\alpha,$ (1)

which ensures that the true value $Y_{t}$ falls within the interval
$\widehat{C}_{t-1}^{\alpha}(X_{t})$ at least $100(1-\alpha)\%$ of the time,
averaged over all instances. On the other hand, conditional coverage is
defined as

$\mathbb{P}(Y_{t}\in\widehat{C}_{t-1}^{\alpha}(X_{t})\mid X_{t})\geq
1-\alpha,$ (2)

which is a stronger guarantee ensuring that given each value of predictor
$X_{t}$, the true value $Y_{t}$ falls within the interval
$\widehat{C}_{t-1}^{\alpha}(X_{t})$ at least $100(1-\alpha)\%$ of the time.

##  3 Method

In this section, we introduce our proposed method, KOWCPI (Kernel-based
Optimally Weighted Conformal Prediction Intervals), which embodies our
approach to enhancing prediction accuracy and robustness in the face of time-
series data. We delve into the methodology and algorithm of KOWCPI in-depth,
highlighting how the Reweighted Nadaraya-Watson (RNW) estimator integrates
with our predictive framework.

Consider prediction for a univariate time series, $Y_{1},Y_{2},\ldots$. We
have predictors $X_{t}$ given to us at time $t$, $t=1,2,\ldots$, which can
depend on the past observations $(Y_{t-1},Y_{t-2},\ldots)$, and possibly other
exogeneous time series $Z_{t}$. Given a pre-trained algorithm $\hat{f}$, we
also have a sequence of non-conformity scores indicating the accuracy of the
prediction:

$\hat{\varepsilon}_{t}=Y_{t}-\hat{f}(X_{t}),\quad t=1,2,\ldots.$

We denote the collection of the past $T$ non-conformity scores at time $t>T$
as

$\mathcal{E}_{t}=(\hat{\varepsilon}_{t-1},\ldots,\hat{\varepsilon}_{t-T}),$

We construct the prediction interval $\widehat{C}_{t-1}^{\alpha}(X_{t})$ with
significance level $\alpha$ at time $t>T$ as follows:

$\displaystyle\begin{split}\widehat{C}_{t-1}^{\alpha}(X_{t})&=[\hat{f}(X_{t})+\hat{q}_{\beta^{*}}(\mathcal{E}_{t}),\hat{f}(X_{t})+\hat{q}_{1-\alpha+\beta^{*}}(\mathcal{E}_{t})],\\\
\beta^{*}&=\mathop{\rm
argmin}_{\beta\in[0,\alpha]}\left(\hat{q}_{1-\alpha+\beta}(\mathcal{E}_{t})-\hat{q}_{\beta}(\mathcal{E}_{t})\right).\end{split}$
(3)

Here, $\hat{q}_{\beta}$ is a quantile regression algorithm that returns an
estimate of the $\beta$-quantile of the residuals, which we will explain
through this section. We consider asymmetrical confidence intervals to ensure
the tightest possible coverage.

###  3.1 Reweighted Nadaraya-Watson estimator

The Reweighted Nadaraya-Watson (RNW) estimator is a general and popular method
for quantile regression. Observe $(\tilde{X}_{i},\tilde{Y}_{i})$,
$i=1,\ldots,n$, where $\tilde{Y}_{i}\in\mathbb{R}$, and the predictors
$\tilde{X}_{i}$ can be $p$-dimensional. The goal is to predict the quantile
$\mathbb{P}(\tilde{Y}\leq b|\tilde{X})$, $b\in\mathbb{R}$, given a test point
$\tilde{X}$ using training samples. The RNW estimator introduces adjustment
weights on the predictors to ensure consistent estimation. We define the
probability-like adjustment weights $p_{i}(\tilde{X})$, $i=1,\ldots,n$, by
maximizing the empirical log-likelihood $\sum_{i=1}^{n}\log p_{i}(\tilde{X})$,
subject to $p_{i}(\tilde{X})\geq 0$, and

$\displaystyle\sum_{i=1}^{n}p_{i}(\tilde{X})=1,$ (4)
$\displaystyle\sum_{i=1}^{n}p_{i}(\tilde{X})(\tilde{X}_{i}-\tilde{X})K_{h}(\tilde{X}_{i}-\tilde{X})=0.$
(5)

The RNW estimate of the conditional CDF $\mathbb{P}(\tilde{Y}\leq
b|\tilde{X})$ is defined as follows:

$\displaystyle\widehat{F}(b|\tilde{X})=\sum_{i=1}^{n}\widehat{W}_{i}(\tilde{X})\mathds{1}(\tilde{Y}_{i}\leq
b),$ (6)

where the weights are given by

$\widehat{W}_{i}(\tilde{X})=\frac{p_{i}(\tilde{X})K_{h}(\tilde{X}_{i}-\tilde{X})}{\sum_{i=1}^{n}p_{i}(\tilde{X})K_{h}(\tilde{X}_{i}-\tilde{X})},\quad
i=1,\ldots,n.$ (7)

Here, the kernel $K\colon\mathbb{R}^{p}\to\mathbb{R}$ is a kernel function,
and $K_{h}(\mathbf{u})=h^{-p}K(h^{-1}\mathbf{u})$ for
$\mathbf{u}\in\mathbb{R}^{p}$. Any reasonable choice of kernel function is
possible; however, to ensure the validity of our theoretical results discussed
in Section 4, the kernel should be nonnegative, bounded, continuous, and
possess compact support. An example is $K(u)=k(\left\|u\right\|)$, where
$k:\mathbb{R}\to\mathbb{R}$ is the Epanechnikov kernel.

The computation of $p_{i}(\tilde{X})$ reduces to a simple one-dimensional
convex minimization problem:

######  Lemma 3.1 ((Hall et al., 1999; Cai, 2001)).

The adjustment weights $p_{i}(\tilde{X})$, $i=1,\ldots,n$, for the RNW
estimator are given as

$p_{i}(\tilde{X})=\frac{1}{n}\left[1+\lambda([\tilde{X}_{i}]_{1}-[\tilde{X}]_{1})K_{h}(\tilde{X}_{i}-\tilde{X})\right]^{-1},$
(8)

where $[X]_{1}$ denotes the first element of a vector $X$, and
$\lambda\in\mathbb{R}$ is the minimizer of:

$L(\lambda;\tilde{X})=-\sum_{i=1}^{n}\log\\{1+\lambda([\tilde{X}_{i}]_{1}-[\tilde{X}]_{1})K_{h}(\tilde{X}_{i}-\tilde{X})\\}.$
(9)

Lemma 3.1 is a starting point for the proof of the asymptotic conditional
coverage property of our algorithm later on.

###  3.2 RNW for conformal prediction

Figure 1: Illustration of Kernel-based Optimally Weighted Conformal Prediction
Intervals (KOWCPI).

To perform the quantile regression for prediction interval construction at
time $t=T+1$, we use a sliding window approach, breaking the past $T$
residuals
$\mathcal{E}_{T+1}=(\hat{\varepsilon}_{T},\ldots\hat{\varepsilon}_{1})$ into
$n\coloneqq(T-w)$ overlapping segments of length $w$. We construct the
predictors and responses to fit the RNW estimator as follows:

$\tilde{Y}_{i}=\hat{\varepsilon}_{i+w},\quad\tilde{X}_{i}=(\hat{\varepsilon}_{i+w-1},\ldots,\hat{\varepsilon}_{i}),\quad
i=1,\ldots,n.$

With RNW estimator $\widehat{F}$ fitted on
$((\tilde{X}_{i},\tilde{Y}_{i}))_{i=1}^{n}$, the conditional $\beta$-quantile
estimator $\widehat{Q}_{\beta}$ is defined as

$\widehat{Q}_{\beta}(\tilde{X})\coloneqq\inf\\{\tilde{y}\in\mathbb{R}\colon\widehat{F}(\tilde{y}|\tilde{X})\geq\beta\\}.$
(10)

After time $t=T+1$, we update $\mathcal{E}_{t}$ by removing the oldest
residual and adding the newest one, then repeat the process (see Algorithm 1).
In Section 4, we prove that due to the consistency of $\widehat{Q}_{\beta}$,
KOWCPI achieves asymptotic conditional coverage despite the significant
temporal dependence introduced by using overlapping segments of residuals.

Algorithm 1 Kernel-based Optimally Weighted Conformal Prediction Intervals
(KOWCPI)

0: Sequentially obtained data $(X_{t},Y_{t})$, $t=1,2,\ldots$, pre-trained
point predictor $\hat{f}$, Target significance level $\alpha\in(0,1)$.

1: Obtain prediction residuals $\hat{\varepsilon}_{t}=Y_{t}-\hat{f}(X_{t})$,
$t=1,\ldots,T$.

2: for $t=T+1,\ldots$ do

3: Update
$\mathcal{E}_{t}=(\hat{\varepsilon}_{t-1},\ldots,\hat{\varepsilon}_{t-T})$.

4: Break $\mathcal{E}_{t}$ into overlapping segments
$\tilde{X}_{i}=(\hat{\varepsilon}_{t-T+i+w-2},\ldots,\hat{\varepsilon}_{t-T+i-1})$,
$i=1,\ldots,n+1$.

5: Obtain $\lambda$ that minimizes $L(\cdot;\tilde{X}_{n+1})$ in (9).

6: Obtain the adjustment weights $p_{i}(\tilde{X}_{n+1}),i=1,\ldots,n$, and
thus final weights $\hat{W}_{i}(\tilde{X}_{n+1})$.

7: With $\tilde{Y}_{i}=\hat{\varepsilon}_{t-T+i+w-1}$, $i=1,\ldots,n$, derive
the quantile estimator $\widehat{Q}_{\beta}(\tilde{X}_{n+1})$ for
$\beta\in[0,\alpha]$.

8: Obtain $\beta^{*}=\mathop{\rm
argmin}_{\beta\in[0,\alpha]}\widehat{Q}_{1-\alpha+\beta}(\tilde{X}_{n+1})-\widehat{Q}_{\beta}(\tilde{X}_{n+1})$.

9: Return
$\widehat{C}_{t-1}^{\alpha}(X_{t})=[\hat{f}(X_{t})+\widehat{Q}_{\beta^{*}}(\tilde{X}_{n+1}),~{}\hat{f}(X_{t})+\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{X}_{n+1})]$.

10: Obtain a new residual $\hat{\varepsilon}_{t}$.

11: end for

#### Bandwidth selection

While it is theoretically possible to calculate the optimal bandwidth that
minimizes the asymptotic mean-squared error, this requires additional
derivative estimation, which significantly complicates the problem.
Consequently, similar to general non-parametric models, one can use cross-
validation to select the bandwidth. However, cross-validation can be
computationally burdensome and may deteriorate under dependent data (Fan et
al., 1995). Therefore, we adapt the non-parametric AIC (Cai and Tiwari, 2000),
used for bandwidth selection in local linear estimators. This method is
applicable because the RNW estimator belongs to the class of linear smoother
(Cai, 2002). We choose the bandwidth $h$ that minimizes

${\rm AIC}_{C}(h)\coloneqq\log({\rm
RSS})+\frac{n+\mbox{tr}(SS^{\top})}{n-(\mbox{tr}(SS^{\top})+2)},$ (11)

where ${\rm RSS}$ is the residual sum of squares, and $S$ is a linear
smoothing operator (Hastie, 1990), with the $(i,j)$-th element given by
$S_{ij}=\widehat{W}_{j}(\tilde{X}_{i}).$

#### Window length selection

To select the window length $w$, cross-validation can be employed. One
approach is to use a weighted sum of the average under-coverage rate and the
average width obtained for a given $w$ as the criterion. Another approach
could involve choosing $w$ with the smallest average width that achieves a
target coverage in the validation set. In experiments, we observed that the
performance is less sensitive to the choice of $w$ across a broader range
compared to the bandwidth $h$, although the selection of $h$ and $w$
interactively affects the performance.

##  4 Theory

In this section, we introduce the theoretical properties of the RNW estimator,
a quantile regression method we use, and demonstrate in Theorem 4.9 that our
KOWCPI asymptotically displays conditional coverage under the strong mixing of
residuals. It turns out that the asymptotic conditional coverage gap can be
derived from known results in the context of kernel quantile regression.

###  4.1 Marginal coverage

We begin by bounding the marginal coverage gap of the KOWCPI method. The
following result shows the coverage gap using our weights, compared with the
oracle weights; the results are established using a similar strategy as in
(Tibshirani et al., 2019, Lemma 3):

######  Proposition 4.1 (Non-asymptotic marginal coverage gap).

Denote by $\mathcal{P}$ the joint density of
$(\tilde{Y}_{1},\ldots,\tilde{Y}_{n+1})$. Then, we have

$\displaystyle\begin{split}&\left|\mathbb{P}\left(Y_{T+1}\in\widehat{C}_{T}^{\alpha}(X_{T+1})\right)-(1-\alpha)\right|\\\
&\leq\frac{1}{2}\left(\mathbb{E}\left\|(W^{*})_{1:n}-\widehat{W}\right\|_{1}+\mathbb{E}W^{*}_{n+1}\right)+2\mathbb{E}\Delta(\tilde{X}_{n+1})\end{split}$
(12)

where $\Delta(\tilde{X}_{n+1})$ is the discrete gap defined in (17), and
$W^{*}$ is the vector of oracle weights with each entry is defined as

$W_{i}^{*}:=\frac{\sum_{\sigma:\sigma(n+1)=i}\mathcal{P}(\tilde{Y}_{\sigma(1)},\ldots,\tilde{Y}_{\sigma(n+1)})}{\sum_{\sigma}\mathcal{P}(\tilde{Y}_{\sigma(1)},\ldots,\tilde{Y}_{\sigma(n+1)})},\quad
i=1,\ldots,{n+1},\quad\mbox{(oracle weights)}$ (13)

and $\sigma$ is a permutation on $\\{1,\ldots,n+1\\}$.

The implication of Proposition 4.1 is that

  * •

The “under-coverage” depends on the $\ell_{1}$-distance between the learned
optimal weights and oracle-optimal weights (that depends on the true joint
distribution of data).

  * •

Note that the oracle weights $W_{t}^{*}$ cannot be evaluated, because in
principle, it requires considering the $(n+1)!$ possible shuffled observed
residuals and their joint distributions.

  * •

The form of the oracle weights $W_{i}^{*}$ from (13) offers an intuitive basis
for algorithm development: we can practically estimate the weights through
quantile regression, utilizing previously observed non-conformity scores.

###  4.2 Conditional coverage

In this section, we derive the asymptotic conditional coverage property of
KOWCPI. For this, we introduce the assumptions necessary for the consistency
of the RNW estimator. To account for the dependency in the data, we assume the
strong mixing of the residual process.

A stationary stochastic process $\\{V_{t}\\}_{t=-\infty}^{\infty}$ on a
probability space with a probability measure $\mathbb{P}$ is said to be
strongly mixing ($\alpha$-mixing) if a mixing coefficient $\alpha(\tau)$
defined as

$\alpha(\tau)=\sup_{A\in\mathfrak{M}_{-\infty}^{0},B\in\mathfrak{M}_{\tau}^{\infty}}\left|\mathbb{P}(A\cap
B)-\mathbb{P}(A)\mathbb{P}(B)\right|$

satisfies $\alpha(\tau)\to 0$ as $\tau\to\infty$, where
$\mathfrak{M}_{s}^{t}$, $-\infty\leq s\leq t\leq\infty$, denotes a
$\sigma$-algebra generated by $\\{V_{s},V_{s+1},\ldots,V_{t}\\}$. The mixing
coefficient $\alpha(\tau)$ quantifies the asymptotic independence between the
past and future of the sequence $\\{V_{t}\\}_{t=-\infty}^{\infty}$.

######  Assumption 4.2 (Mixing of the process).

The stationary process $(V_{i}=(\tilde{X}_{i},\tilde{Y}_{i}))_{i=1}^{\infty}$
is strongly mixing with the mixing coefficient
$\alpha(\tau)=\mathcal{O}(\tau^{-(2+\delta)})$ for some $\delta>0$.

Due to stationarity, the conditional CDF of the realized residual does not
depend on the index $i$; thus, denote

$F(b|\mathbf{z})=\mathbb{P}(\tilde{Y}_{i}\leq b|\tilde{X}_{i}=\mathbf{z}),$

as the conditional CDF of the random variable $\tilde{Y}_{i}$ given
$\tilde{X}_{i}=\mathbf{z}$. In addition, we introduce the following notations:

  * •

Let $g(\mathbf{z})$ be the marginal density of $\tilde{X}_{i}$ at
$\mathbf{z}$. (Note that due to stationarity, we can have a common marginal
density.)

  * •

Let $g_{1,i},i\geq 2$ denote the joint density of $\tilde{X}_{1}$ and
$\tilde{X}_{i}$.

The following assumptions (4.3-4.5) are common in nonparametric statistics,
essential for attaining desirable properties such as the consistency of an
estimator (Tsybakov, 2009).

######  Assumption 4.3 (Smoothness of the conditional CDF and densities).

For fixed $\tilde{y}\in\mathbb{R}$ and $\mathbf{z}\in\mathbb{R}^{w}$,

  1. (i)

$0<F(\tilde{y}|\mathbf{z})<1$.

  2. (ii)

$F(\tilde{y}|\mathbf{z})$ is twice continuously partially differentiable with
respect to $\mathbf{z}$.

  3. (iii)

$g(\mathbf{z})>0$ and $g(\cdot)$ is continuous at $\mathbf{z}$.

  4. (iv)

There exists $M>0$ such that $\left|g_{1,i}(u,v)-g(u)g(v)\right|\leq M$ for
all $u,v$ and $i\geq 2$.

Regarding Assumption 4.3, we would like to remark that there is a negative
result: without additional assumptions about the distribution, it is
impossible to construct finite-length prediction intervals that satisfy
conditional coverage (Lei and Wasserman, 2014; Vovk, 2012).

######  Assumption 4.4 (Regularity of the kernel function).

The kernel $K:\mathbb{R}^{w}\to\mathbb{R}$ is a nonnegative, bounded,
continuous, and compactly supported density function satisfying

  1. (i)

$\int_{\mathbb{R}^{w}}\mathbf{u}K(\mathbf{u})d\mathbf{u}=0$,

  2. (ii)

$\int_{\mathbb{R}^{w}}\mathbf{u}\mathbf{u}^{\top}K(\mathbf{u})d\mathbf{u}=\mu_{2}I$
for some $\mu_{2}\in(0,\infty)$,

  3. (iii)

$\int_{\mathbb{R}^{w}}K^{2}(\mathbf{u})d\mathbf{u}=\nu_{0}$ and
$\int_{\mathbb{R}^{w}}\mathbf{u}\mathbf{u}^{\top}K^{2}(\mathbf{u})d\mathbf{u}=\nu_{2}I$
for some $\nu_{0},\nu_{2}\in(0,\infty)$.

Assumptions 4.4-(i), (ii), and (iii) are standard conditions (Wand and Jones,
1994) that require $K$ to be “symmetric” in a sense that that the weighting
scheme relies solely on the distance between the observation and the test
point. For example, if $K$ is isometric, i.e., $K(u)=k(\left\|u\right\|)$ for
some univariate kernel function $k:\mathbb{R}\to\mathbb{R}$, it can satisfy
these conditions using widely adopted kernels such as the Epanechnikov kernel.

######  Assumption 4.5 (Bandwidth selection).

As $n\to\infty$, the bandwidth $h$ satisfies

$h\to 0,\mbox{ and }\,nh^{w(1+2/\delta)}\to\infty.$

We note that Assumption 4.5 is met when selecting the (theoretically) optimal
bandwidth $h^{*}\sim n^{-1/(w+4)}$, which minimizes the asymptotic mean
squared error (AMSE) of the RNW estimator, provided that $\delta>1/2$.

We prove the following proposition following a similar strategy as (Salha,
2006) by fixing several technical details:

######  Proposition 4.6 (Consistency of the RNW estimator).

Under Assumptions 4.2-4.5, given arbitrary $\tilde{x}$ and $\tilde{y}$, as
$n\rightarrow\infty$,

$\displaystyle\widehat{F}(\tilde{y}|\mathbf{z})-F(\tilde{y}|\mathbf{z})=\frac{1}{2}h^{2}\mbox{tr}(D_{\mathbf{z}}^{2}(F(\tilde{y}|\mathbf{z})))\mu_{2}+o_{p}(h^{2})+\mathcal{O}_{p}((nh^{w})^{-1/2}).$
(14)

This proposition implies pointwise convergence in probability of the RNW
estimator, and since it is the weighted empirical CDF, this pointwise
convergence implies uniform convergence in probability (Tucker, 1967,
p.127-128). Consequently, we obtain the consistency of the conditional
quantile estimator in (10) to the true conditional quantile given as

$Q_{\beta}(\tilde{x})=\inf\\{\tilde{y}\in\mathbb{R}:F(\tilde{y}|\tilde{x})\geq\beta\\}.$

######  Corollary 4.7.

Under Assumptions 4.2-4.5, for every $\beta\in(0,1)$ and $\mathbf{z}$, as
$n\to\infty$,

$\widehat{Q}_{\beta}(\mathbf{z})\to Q_{\beta}(\mathbf{z})\mbox{ in
probability}.$ (15)

As a direct consequence of Corollary 4.7, the asymptotic conditional coverage
of KOWCPI is guaranteed by the consistency of the quantile estimator used in
our sequential algorithm.

######  Corollary 4.8 (Asymptotic conditional coverage guarantee).

Under Assumptions 4.2-4.5, for any $\alpha\in(0,1)$, as $n\to\infty$,

$\mathbb{P}(Y_{t}\in\hat{C}_{t-1}^{\alpha}(X_{t})|X_{t})\to(1-\alpha)\mbox{ in
probability}.$ (16)

Thus, employing quantile regression using the RNW estimator for prediction
residuals derived from the time-series data of continuous random variables,
assuming strong mixing of these residuals, KOWCPI can achieve approximate
conditional coverage with a sufficient number of residuals utilized.

To further specify the rate of convergence, define the discrete gap

$\Delta(\tilde{X})=\sup_{\beta\in[0,1]}|{\widehat{F}(\widehat{Q}_{\beta}(\tilde{X})|\tilde{X})-\beta}|=\max_{i=1,\ldots,n}\widehat{W}_{i}(\tilde{X}).$
(17)

######  Theorem 4.9 (Conditional coverage gap).

Under Assumptions 4.2-4.5, for any $\alpha\in(0,1)$ and $x_{t}$, as
$n\to\infty$,

$\left|\mathbb{P}\left(Y_{t}\in\widehat{C}_{t-1}^{\alpha}(x_{t})\mid
X_{t}=x_{t}\right)-(1-\alpha)\right|\leq\mathcal{O}_{p}(h^{2}+(nh^{w})^{-1/2})+2\Delta(\tilde{x}_{n+1}),$
(18)

where $\tilde{x}_{n+1}$ is the realization of $\tilde{X}_{n+1}$ given
$X_{t}=x_{t}$.

Given that the adjustment weights $p_{i}(\mathbf{z})$ uniformly concentrate to
$1/n$ (Steikert, 2014), one can see that the conditional coverage gap tends to
zero, although its precise rate remains an open question.

##  5 Experiments

Table 1: Empirical marginal coverage and width across three real time-series
datasets by different methods. The target coverage is $1-\alpha=0.9$.
| Electric | Wind | Solar
---|---|---|---
| Coverage | Width | Coverage | Width | Coverage | Width
KOWCPI | 0.90 | 0.22 | 0.91 | 2.41 | 0.90 | 48.8
Plain NW | 0.89 | 0.31 | 0.95 | 3.58 | 0.41 | 20.1
SPCI | 0.90 | 0.29 | 0.94 | 2.61 | 0.92 | 84.2
EnbPI | 0.93 | 0.36 | 0.92 | 5.25 | 0.87 | 106.0
ACI | 0.89 | 0.32 | 0.88 | 8.26 | 0.89 | 143.9
FACI | 0.89 | 0.28 | 0.91 | 7.77 | 0.89 | 141.9
AgACI | 0.91 | 0.30 | 0.88 | 7.54 | 0.90 | 144.6
SF-OGD | 0.79 | 0.25 | 0.11 | 0.29 | 0 | 0.50
SAOCF | 0.93 | 0.33 | 0.76 | 4.00 | 0.64 | 33.5
SCP | 0.87 | 0.30 | 0.86 | 8.20 | 0.89 | 142.0

In this section, we compare the performance of KOWCPI against state-of-the-art
conformal prediction baselines using real time-series data. We aim to show
that KOWCPI can consistently reach valid coverage with the narrowest
prediction intervals.

Dataset. We consider three real time series from different domains. The first
ELEC2 data set (eletric) (Harries, 1999) tracks electricity usage and pricing
in the states of New South Wales and Victoria in Australia for every 30
minutes over a 2.5-year period in 1996–1999. The second renewable energy data
(solar) (Zhang et al., 2021) are from the National Solar Radiation Database
and contain hourly solar radiation data (measured in GHI) from Atlanta in
2018. The third wind speed data (wind) (Zhu et al., 2021) are collected at
wind farms operated by MISO in the US. The wind speed record was updated every
15 minutes over a one-week period in September 2020.

Baselines. We consider Sequential Predictive Conformal Inference (SPCI) (Xu
and Xie, 2023b), Ensemble Prediction Interval (EnbPI) (Xu and Xie, 2023a),
Adaptive Conformal Inference (ACI) (Gibbs and Candès, 2021), Aggregated ACI
(AgACI) (Zaffran et al., 2022), Fully Adaptive Conformal Inference (FACI)
(Gibbs and Candès, 2022), Scale-Free Online Gradient Descent (SF-OGD) (Orabona
and Pál, 2018; Bhatnagar et al., 2023), Strongly Adaptive Online Conformal
Prediction (SAOCP) (Bhatnagar et al., 2023), and vanilla Split Conformal
Prediction (SCP) (Vovk et al., 2005). Additionally, we included a comparison
where weights were derived from the original Nadaraya-Watson estimator (Plain
NW). For the implementation of ACI-related methods, we utilized the R package
AdaptiveConformal (https://github.com/herbps10/AdaptiveConformal). For SPCI
and EnbPI, we used the Python code from https://github.com/hamrel-cxu/SPCI-
code.

Setup and evaluation metrics. In all comparisons, we use the random forest as
the base point predictor with the number of trees $=10$. Every dataset is
split in a 7:1:2 ratio for training the point predictor, tuning the window
length $w$ and bandwidth $h$, and constructing prediction intervals,
respectively. The window length for each dataset is fixed and determined
through cross-validation, while the bandwidth is selected by minimizing the
nonparametric AIC, as detailed in (11).

Besides examining marginal coverage and widths of prediction intervals on test
data, we also focus on rolling coverage, which is helpful in showing
approximate conditional coverage at specific time indices. Given a rolling
window size $m>0$, rolling coverage $\widehat{\rm RC}_{t}$ at time $t$ is
defined as $\widehat{\rm
RC}_{t}=\frac{1}{m}\sum_{i=1}^{m}\mathds{1}\\{Y_{t-i+1}\in\widehat{C}_{t-i}^{\alpha}(X_{t-i+1})\\}.$

(a) Rolling coverage (boxplot)

(b) Width of prediction intervals (boxplot)

(c) Rolling coverage over prediction time indices

(d) $\log(\widehat{W}_{i})$ by the time lag

Figure 2: Comparison of empirical rolling coverage and width on the electric
dataset by different methods in (a) rolling coverage; (b) widths of intervals,
(c) rolling coverage over time, and (d) an instance of computed adjusted
weights, which shows that the weights do not necessarily decay over lag. The
target coverage is 90%. In (a), the red dotted line is the target coverage and
in (b), the blue dotted line is the median width of KOWCPI.

Results. The empirical marginal coverage and width results for all methods are
summarized in Table 1. The results indicate that KOWCPI consistently achieves
the 90% target coverage and maintains the smallest average width compared to
the alternative state-of-the-art methods. While all methods, except SF-OGD,
SAOCF, and Plain NW nearly achieve marginal coverage under target
$1-\alpha=0.9$, KOWCPI produces the narrowest average width on all datasets.
In terms of rolling results, we show in Figures 2 and 2 that the coverage of
KOWCPI intervals consistently centers around 90% throughout the entire test
phase. Additionally, Figure 2 shows that KOWCPI intervals are also
significantly narrower with a smaller variance than the baselines. Lastly,
Figure 2 shows weights $\widehat{W}$ (in log scale) by the RNW estimator at
the first time index of test data. We note that the most recent set of non-
conformity scores (in terms of time indices) receives the heaviest weighting,
which is natural due to the greatest similarity between the first test datum
and training data closest to it in time. In Appendix C, we show additional
comparisons of KOWCPI against the baselines on the other two datasets in terms
of rolling results, where KOWCPI remains consistently better.

##  6 Conclusion

In this paper, we introduced KOWCPI, a method to sequentially construct
prediction intervals for time-series data. By incorporating the classical
Reweighted Nadaraya-Watson estimator into the weighted conformal prediction
framework, KOWCPI effectively adapts to the dependent structure of time-series
data by utilizing data-driven adaptive weights. Our theoretical contributions
include providing theoretical guarantees for the asymptotic conditional
coverage of KOWCPI under strong mixing conditions and bounding the marginal
and conditional coverage gaps. Empirical validation on real-world time-series
datasets demonstrated the effectiveness of KOWCPI compared to state-of-the-art
methods, achieving narrower prediction intervals without compromising
empirical coverage.

Future work could explore adaptive window selection, where the size of the
non-conformity score batch is adjusted dynamically to capture shifts in the
underlying distribution. Additionally, the natural compatibility of kernel
regression with multivariate data can be leveraged to expand the utility of
KOWCPI for multivariate time-series data, as detailed in Appendix A. There is
also potential for improving theoretical guarantees and improving practical
performance by designing alternative non-conformity scores.

## Acknowledgments

This work is partially supported by an NSF CAREER CCF-1650913, NSF
DMS-2134037, CMMI-2015787, CMMI-2112533, DMS-1938106, DMS-1830210, and the
Coca-Cola Foundation.

## References

  * Abdous and Theodorescu [1992] B. Abdous and R. Theodorescu.  Note on the spatial quantile of a random vector.  _Statist. Probab. Lett._ , 13(4):333–336, 1992\.
  * Angelopoulos et al. [2023] A. Angelopoulos, E. Candès, and R. J. Tibshirani.  Conformal PID Control for Time Series Prediction.  In _Advances in Neural Information Processing Systems_ , 2023.
  * Angelopoulos and Bates [2023] A. N. Angelopoulos and S. Bates.  Conformal Prediction: A Gentle Introduction.  _Foundations and Trends® in Machine Learning_ , 16(4):494–591, 2023.
  * Auer et al. [2023] A. Auer, M. Gauch, D. Klotz, and S. Hochreiter.  Conformal Prediction for Time Series with Modern Hopfield Networks.  In _Advances in Neural Information Processing Systems_ , 2023.
  * Barber et al. [2023] R. F. Barber, E. J. Candès, A. Ramdas, and R. J. Tibshirani.  Conformal prediction beyond exchangeability.  _The Annals of Statistics_ , 51(2):816 – 845, 2023.
  * Bhatnagar et al. [2023] A. Bhatnagar, H. Wang, C. Xiong, and Y. Bai.  Improved Online Conformal Prediction via Strongly Adaptive Online Learning.  In _Proceedings of the 40th International Conference on Machine Learning_ , 2023.
  * Cai [2001] Z. Cai.  Weighted Nadaraya-Watson regression estimation.  _Statistics & Probability Letters_, 51(3):307–318, 2001.
  * Cai [2002] Z. Cai.  Regression quantiles for time series.  _Econometric Theory_ , 18(1):169–192, 2002.
  * Cai and Tiwari [2000] Z. Cai and R. C. Tiwari.  Application of a local linear autoregressive model to BOD time series.  _Environmetrics_ , 11(3):341–350, 2000.
  * Fan et al. [1995] J. Fan, N. E. Heckman, and M. P. Wand.  Local Polynomial Kernel Regression for Generalized Linear Models and Quasi-Likelihood Functions.  _Journal of the American Statistical Association_ , 90(429):141–150, 1995.
  * Feldman et al. [2022] S. Feldman, L. Ringel, S. Bates, and Y. Romano.  Achieving Risk Control in Online Learning Settings.  _arXiv preprint arXiv:2205.09095_ , 2022.
  * Gibbs and Candès [2021] I. Gibbs and E. Candès.  Adaptive Conformal Inference Under Distribution Shift.  _Advances in Neural Information Processing Systems_ , 34:1660–1672, 2021.
  * Gibbs and Candès [2022] I. Gibbs and E. Candès.  Conformal Inference for Online Prediction with Arbitrary Distribution Shifts.  _arXiv preprint arXiv:2208.08401_ , 2022.
  * Gibbs et al. [2023] I. Gibbs, J. J. Cherian, and E. J. Candès.  Conformal Prediction With Conditional Guarantees.  _arXiv preprint arXiv:2305.12616_ , 2023.
  * Guan [2023] L. Guan.  Localized conformal prediction: a generalized inference framework for conformal prediction.  _Biometrika_ , 110(1):33–50, 2023.
  * Hall et al. [1999] P. Hall, R. C. Wolff, and Q. Yao.  Methods for Estimating a Conditional Distribution Function.  _Journal of the American Statistical Association_ , 94(445):154–163, 1999.
  * Harries [1999] M. Harries.  SPLICE-2 Comparative Evaluation: Electricity Pricing.  Technical report, University of New South Wales, School of Computer Science and Engineering, 1999.
  * Hastie [1990] T. J. Hastie.  _Generalized additive models_.  CRC Press, 1990.
  * Ibragimov et al. [1971] I. A. Ibragimov, Y. V. Linnik, and J. F. C. Kingman.  _Independent and Stationary Sequences of Random Variables_.  Wolters-Noordhoff., 1971.
  * Kim et al. [2020] B. Kim, C. Xu, and R. Barber.  Predictive inference is free with the jackknife+-after-bootstrap.  _Advances in Neural Information Processing Systems_ , 33:4138–4149, 2020.
  * Kivaranovic et al. [2020] D. Kivaranovic, K. D. Johnson, and H. Leeb.  Adaptive, Distribution-Free Prediction Intervals for Deep Networks.  In _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , 2020.
  * Lee et al. [2023] Y. Lee, R. F. Barber, and R. Willett.  Distribution-free inference with hierarchical data.  _arXiv preprint arXiv:2306.06342_ , 2023.
  * Lei and Wasserman [2014] J. Lei and L. Wasserman.  Distribution-free prediction bands for non-parametric regression.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 76(1):71–96, 2014.
  * Lei et al. [2013] J. Lei, J. Robins, and L. Wasserman.  Distribution-Free Prediction Sets.  _Journal of the American Statistical Association_ , 108(501):278–287, 2013.
  * Masry [1986] E. Masry.  Recursive probability density estimation for weakly dependent stationary processes.  _IEEE Transactions on Information Theory_ , 32(2):254–267, 1986.
  * Nadaraya [1964] E. A. Nadaraya.  On Estimating Regression.  _Theory of Probability & Its Applications_, 9(1):141–142, 1964.
  * Orabona and Pál [2018] F. Orabona and D. Pál.  Scale-free online learning.  _Theoretical Computer Science_ , 716:50–69, 2018.
  * Romano et al. [2019] Y. Romano, E. Patterson, and E. Candès.  Conformalized Quantile Regression.  In _Advances in Neural Information Processing Systems_ , 2019.
  * Salha [2006] R. Salha.  _Kernel Estimation for the Conditional Mode and Quantiles of Time Series_.  PhD thesis, University of Macedonia, 2006.
  * Sesia and Candès [2020] M. Sesia and E. J. Candès.  A comparison of some conformal quantile regression methods.  _Stat_ , 9(1):e261, 2020.
  * Stankevičiūtė et al. [2021] K. Stankevičiūtė, A. Alaa, and M. van der Schaar.  Conformal time-series forecasting.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Steikert [2014] K. U. Steikert.  _The weighted Nadaraya-Watson Estimator: Strong consistency results, rates of convergence, and a local bootstrap procedure to select the bandwidth_.  PhD thesis, University of Zurich, 2014.
  * Sun and Yu [2024] S. H. Sun and R. Yu.  Copula conformal prediction for multi-step time series prediction.  In _The Twelfth International Conference on Learning Representations_ , 2024.
  * Tibshirani et al. [2019] R. J. Tibshirani, R. Foygel Barber, E. Candès, and A. Ramdas.  Conformal Prediction Under Covariate Shift.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Tsybakov [2009] A. B. Tsybakov.  _Introduction to Nonparametric Estimation_.  Springer, 2009.
  * Tucker [1967] H. G. Tucker.  _A Graduate Course in Probability_.  Academic Press, 1967.
  * Vovk [2012] V. Vovk.  Conditional validity of inductive conformal predictors.  In _Proceedings of the Asian Conference on Machine Learning_ , 2012\.
  * Vovk et al. [1999] V. Vovk, A. Gammerman, and C. Saunders.  Machine-Learning Applications of Algorithmic Randomness.  In _Proceedings of the Sixteenth International Conference on Machine Learning_ , 1999.
  * Vovk et al. [2005] V. Vovk, A. Gammerman, and G. Shafer.  _Algorithmic Learning in a Random World_.  Springer, 2005.
  * Wand and Jones [1994] M. P. Wand and M. C. Jones.  _Kernel Smoothing_.  CRC press, 1994.
  * Watson [1964] G. S. Watson.  Smooth regression analysis.  _Sankhyā: The Indian Journal of Statistics, Series A_ , 26:359–372, 1964.
  * Xu and Xie [2021a] C. Xu and Y. Xie.  Conformal Anomaly Detection on Spatio-Temporal Observations with Missing Data.  _arXiv preprint arXiv:2105.11886_ , 2021a.
  * Xu and Xie [2021b] C. Xu and Y. Xie.  Conformal prediction interval for dynamic time-series.  In _Proceedings of the 38th International Conference on Machine Learning_ , 2021b.
  * Xu and Xie [2022] C. Xu and Y. Xie.  Conformal prediction set for time-series.  _arXiv preprint arXiv:2206.07851_ , 2022.
  * Xu and Xie [2023a] C. Xu and Y. Xie.  Conformal prediction for time series.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 45(10):11575–11587, 2023a.
  * Xu and Xie [2023b] C. Xu and Y. Xie.  Sequential Predictive Conformal Inference for Time Series.  In _Proceedings of the 40th International Conference on Machine Learning_ , 2023b.
  * Xu et al. [2024] C. Xu, H. Jiang, and Y. Xie.  Conformal prediction for multi-dimensional time series by ellipsoidal sets.  _arXiv preprint arXiv:2403.03850_ , 2024.
  * Yu and Jones [1998] K. Yu and M. C. Jones.  Local Linear Quantile Regression.  _Journal of the American Statistical Association_ , 93(441):228–237, 1998.
  * Zaffran et al. [2022] M. Zaffran, O. Féron, Y. Goude, J. Josse, and A. Dieuleveut.  Adaptive Conformal Predictions for Time Series.  In _Proceedings of the 39th International Conference on Machine Learning_ , 2022.
  * Zhang et al. [2021] M. Zhang, C. Xu, A. Sun, F. Qiu, and Y. Xie.  Solar Radiation Ramping Events Modeling Using Spatio-temporal Point Processes.  _arXiv preprint arXiv:2101.11179_ , 2021.
  * Zhu et al. [2021] S. Zhu, H. Zhang, Y. Xie, and P. Van Hentenryck.  Multi-resolution spatio-temporal prediction with application to wind power generation.  _arXiv preprint arXiv:2108.13285_ , 2021.

##  Appendix A Multivariate time series

In the main text, our discussion has centered on cases where the response
variables $Y_{t}$ are scalars. Here, we explore the natural extension of our
methodology to handle scenarios with multivariate responses. This extension
requires defining multivariate quantiles, introducing a multivariate version
of the RNW estimator for estimating these quantiles [Salha, 2006], and
adapting our KOWCPI method for multivariate responses.

#### Multivariate conditional quantiles

Consider a strongly mixing stationary process
$((\tilde{X}_{i},\tilde{Y}_{i}))_{i=1}^{\infty}$, which is a realization of
random variable $(\tilde{X},\tilde{Y})\in\mathbb{R}^{p}\times\mathbb{R}^{s}$.
Following Abdous and Theodorescu [1992], we first define a pseudo-norm
function $\left\|\cdot\right\|_{2,\alpha}:\mathbb{R}^{s}\to\mathbb{R}$ for
$\alpha\in(0,1)$ as

$\left\|v\right\|_{2,\alpha}=\left\|\left(\frac{\left|v_{1}\right|+(2\alpha-1)}{2},\ldots,\frac{\left|v_{s}\right|+(2\alpha-1)}{2}\right)\right\|_{2},$

for $v\in\mathbb{R}^{s}$, where $\left\|\cdot\right\|_{2}$ is the Euclidean
norm on $\mathbb{R}^{s}$. Let

$H_{\alpha}(\theta,\tilde{x})\coloneqq\mathbb{E}[\|\tilde{Y}-\theta\|_{2,\alpha}-\|\tilde{Y}\|_{2,\alpha}\mid\tilde{X}=\tilde{x}].$

######  Definition A.1 (Multivariate conditional quantile [Abdous and
Theodorescu, 1992]).

Define a multivariate conditional $\alpha$-quantile
$\theta_{\alpha}(\tilde{x})$ for $\alpha\in(0,1)$ as

$\theta_{\alpha}(\tilde{x})=\mathop{\rm
argmin}_{\theta\in\mathbb{R}^{s}}H_{\alpha}(\theta,\tilde{x}).$ (A.1)

######  Remark A.2 (Compatibility with univariate quantile function).

For a scalar $\tilde{Y}\in\mathbb{R}$, its conditional quantile given
$\tilde{X}=\tilde{x}$ is

$\displaystyle\theta_{\alpha}(\tilde{x})$ $\displaystyle=\mathop{\rm
argmin}_{\theta\in\mathbb{R}}\mathbb{E}\left[\|{\tilde{Y}-\theta}\|_{2,\alpha}\mid\tilde{X}=\tilde{x}\right]$
$\displaystyle=\mathop{\rm
argmin}_{\theta\in\mathbb{R}}\mathbb{E}\left[|{\tilde{Y}-\theta}|+(2\alpha-1)(\tilde{Y}-\theta)\mid\tilde{X}=\tilde{x}\right]$
$\displaystyle=\mathop{\rm
argmin}_{\theta\in\mathbb{R}}\mathbb{E}((\tilde{Y}-\theta)(\alpha-\mathds{1}(\tilde{Y}\leq\theta))\mid\tilde{X}=\tilde{x})$
$\displaystyle=Q_{\alpha}(\tilde{x}),$

for any $\alpha\in(0,1)$. Thus, Definition A.1 is consistent with the
univariate case.

#### Multivariate RNW estimator

Following the definition in (6), we obtain the RNW estimator for multivariate
responses as

$\widehat{F}(\tilde{y}|\tilde{x})=\frac{\sum_{i=1}^{n}W(\frac{\tilde{X}_{i}-\tilde{x}}{h})\mathds{1}(\tilde{Y}_{i}\leq\tilde{y})}{\sum_{i=1}^{n}W(\frac{\tilde{X}_{i}-\tilde{x}}{h})},$
(A.2)

where, according to (7) and (8),

$W(u)=\frac{K(u)}{1+\lambda u_{1}K(u)},$

for $u=(u_{1},\ldots,u_{p})^{\top}$. Now, let $W_{h}(u)=h^{-p}W(h^{-1}u)$, and
define an estimator for $H_{\alpha}(\theta,x)$ as

$\widehat{H}_{\alpha}(\theta,\tilde{x})\coloneqq\int_{\mathbb{R}^{s}}(\|{\tilde{y}-\theta}\|_{2,\alpha}-\|{\tilde{y}}\|_{2,\alpha})\widehat{F}(d\tilde{y}|\tilde{x})=\frac{\sum_{i=1}^{n}W_{h}(\tilde{X}_{i}-\tilde{x})\left(\|{\tilde{Y}_{i}-\theta}\|_{2,\alpha}-\|{\tilde{Y}_{i}}\|_{2,\alpha}\right)}{\sum_{i=1}^{n}W_{h}(\tilde{X}_{i}-\tilde{x})},$

and consequently the RNW conditional quantile estimator
$\widehat{\theta}_{\alpha}$ as

$\displaystyle\begin{split}\widehat{\theta}_{\alpha}(\tilde{x})\coloneqq\mathop{\rm
argmin}_{\theta\in\mathbb{R}^{p}}\hat{H}_{\alpha}(\theta,\tilde{x})=\mathop{\rm
argmin}_{\theta\in\mathbb{R}^{p}}\sum_{i=1}^{n}W_{h}(\tilde{X}_{i}-\tilde{x})(\|{\tilde{Y}_{i}-\theta}\|_{2,\alpha}-\|{\tilde{Y}_{i}}\|_{2,\alpha}).\end{split}$
(A.3)

#### Multivariate KOWCPI

Suppose we are sequentially observing
$(X_{t},Y_{t})\in\mathbb{R}^{d}\times\mathbb{R}^{s}$, $t=1,2,\ldots$. Based on
the construction of the multivariate version of the RNW estimator, we can
extend our KOWCPI approach to multivariate responses in the same manner as
described in Algorithm 1, with multivariate residuals

$\hat{\varepsilon}_{t}=Y_{t}-\hat{f}(X_{t})\in\mathbb{R}^{s},$

as non-conformity scores. This adaptation allows for the application of our
methodology to a broader range of data scenarios involving dependent data with
multivariate response variables, which were similarly studied in [Xu et al.,
2024, Sun and Yu, 2024, Stankevičiūtė et al., 2021].

##  Appendix B Proofs

The following lemma is adapted from the proof of Lemma 1 of Tibshirani et al.
[2019]; however, we do not assume exchangeability.

######  Lemma B.1 (Weights on quantile for non-exchangeable data).

Given a sequence of random variables $\\{V_{1},\ldots,V_{n+1}\\}$ with joint
density $\mathcal{P}$ and a sequence of observations
$\\{v_{1},\ldots,v_{n+1}\\}$. Define the event

$E=\\{\\{V_{1},\ldots,V_{n+1}\\}=\\{v_{1},\ldots,v_{n+1}\\}\\}.$

Then we have for $i=1,\ldots,n+1$,

$\mathbb{P}\\{V_{n+1}=v_{i}|E\\}=\frac{\sum_{\sigma:\sigma(n+1)=i}\mathcal{P}(v_{\sigma(1)},\ldots,v_{\sigma(n+1)})}{\sum_{\sigma}\mathcal{P}(v_{\sigma(1)},\ldots,v_{\sigma(n+1)})}\in[0,1].$

Note that when the residuals are exchangeable, $W_{i}^{*}=1/(n+1)$, as also
observed in Tibshirani et al. [2019]. Now we prove Proposition 4.1.

###### Proof of Proposition 4.1.

The proof assumes that $\tilde{Y}_{i}$, for $i=1,\ldots,n+1$, are almost
surely distinct. However, the proof remains valid, albeit with more complex
notations involving multisets, if this is not the case. Denote by ${\rm
Quantile}_{\beta}(\mathbb{Q})$ the $\beta$-quantile of the distribution
$\mathbb{Q}$ on $\mathbb{R}$, and by $\delta_{a}$ the point mass distribution
at $a\in\mathbb{R}$. Define the event
$E=\\{\\{\tilde{Y}_{1},\ldots,\tilde{Y}_{n+1}\\}=\\{v_{1},\ldots,v_{n+1}\\}\\}$.
Then, by the tower property, we have

$\displaystyle\mathbb{P}(Y_{T+1}\in\widehat{C}_{T}^{\alpha}(X_{T+1}))$
$\displaystyle=\mathbb{E}(\mathbb{P}(Y_{T+1}\in\widehat{C}_{T}^{\alpha}(X_{T+1})\mid
E))$ $\displaystyle=\left.\mathbb{E}\left[\mathbb{P}\left({\rm
Quantile}_{\beta^{*}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\leq\tilde{Y}_{n+1}\leq{\rm
Quantile}_{1-\alpha+{\beta^{*}}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\right|E\right)\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{P}_{V\sim P^{W^{*}}}\left({\rm
Quantile}_{\beta^{*}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\leq
V\leq{\rm
Quantile}_{1-\alpha+{\beta^{*}}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\right)\right],$

where $P^{W^{*}}=\sum_{i=1}^{n+1}W^{*}_{i}\delta_{v_{i}}$, and in the last
line, we have used the result from Lemma B.1,

$\tilde{Y}_{n+1}|E\sim P^{W^{*}}.$

Denote the weighted empirical distributions based on
$\widehat{W}=\widehat{W}(\tilde{X}_{n+1})$ as

$\displaystyle P^{\widehat{W}}$
$\displaystyle=\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}.$

This gives the marginal coverage gap as

$\displaystyle\left|\mathbb{P}(Y_{n+1}\in\widehat{C}_{n}^{\alpha}(X_{n+1}))-(1-\alpha)\right|$
$\displaystyle\leq\mathbb{E}\left[\left|\mathbb{P}_{V\sim
P^{W^{*}}}\Bigg{(}{\rm
Quantile}_{\beta^{*}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\leq
V\leq{\rm
Quantile}_{1-\alpha+{\beta^{*}}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\Bigg{)}\right.\right.$
$\displaystyle\quad\quad\left.-\left.\mathbb{P}_{V\sim
P^{\widehat{W}}}\left({\rm
Quantile}_{\beta^{*}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\leq
V\leq{\rm
Quantile}_{1-\alpha+{\beta^{*}}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\right)\right|\right]$
$\displaystyle\quad+\mathbb{E}\left|\mathbb{P}_{V\sim
P^{\widehat{W}}}\left({\rm
Quantile}_{\beta^{*}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\leq
V\leq{\rm
Quantile}_{1-\alpha+{\beta^{*}}}\left(\sum_{i=1}^{n}\widehat{W}_{i}\delta_{v_{i}}\right)\right)-(1-\alpha)\right|$
$\displaystyle\leq\mathbb{E}[d_{\rm TV}(P^{W^{*}},P^{\widehat{W}})]$
$\displaystyle\quad+\mathbb{E}\left|\widehat{F}\left(\widehat{Q}_{1-\alpha+{\beta^{*}}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right)-(1-\alpha+{\beta^{*}})\right|+\mathbb{E}\left|\widehat{F}\left(\widehat{Q}_{{\beta^{*}}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right)-{\beta^{*}}\right|$
$\displaystyle\leq\frac{1}{2}\mathbb{E}\|(W^{*})_{1:n}-\widehat{W}\|_{1}+\frac{1}{2}\mathbb{E}W^{*}_{n+1}+2\mathbb{E}\max_{i=1,\ldots,n}\widehat{W}_{i}(\tilde{X}_{n+1}),$

where we denote by $d_{\rm TV}(\cdot,\cdot)$ the total variation distance
between probability measures, and the second inequality is due to the
definition of the total variation distance. ∎

###  B.1 Proof of asymptotic conditional coverage of KOWCPI (Corollary 4.8 and
Theorem 4.9)

In deriving the asymptotic conditional coverage property of KOWCPI, the
consistency of the RNW estimator plays a crucial role. Therefore, we first
introduce the proof of Proposition 4.6, which discusses the consistency of the
CDF estimator. Corollary 4.7, which states the consistency of the quantile
estimator, is a natural consequence of Proposition 4.6 and leads us to the
proof for our main results, Corollary 4.8 and Theorem 4.9. Proof of
Proposition 4.6 adopts the similar strategy as Salha [2006] and Cai [2002].

To prove Proposition 4.6, it is essential to first understand the nature of
the adjustment weight $p_{i}(\mathbf{z})$. Thus, Lemma 3.1 is not only crucial
for the practical implementation of the RNW estimator but also indispensable
in the proof process of Proposition 4.6.

###### Proof of Lemma 3.1.

For display purposes, denote $[X]_{1}$ as $X_{1}$. By (5), we have that

$\sum_{i=1}^{n}p_{i}(\mathbf{z})(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z})=0.$
(A.4)

Let

$\displaystyle\mathcal{L}(\lambda_{1},\lambda_{2},p_{1}(\mathbf{z}),\ldots,p_{n}(\mathbf{z}))$
$\displaystyle=\sum_{i=1}^{n}\log
p_{i}(\mathbf{z})+\lambda_{1}\Bigg{(}1-\sum_{i=1}^{n}p_{i}(\mathbf{z})\Bigg{)}+\lambda_{2}\sum_{i=1}^{n}p_{i}(\mathbf{z})(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z}),$

where $\lambda_{1},\lambda_{2}\in\mathbb{R}$ are the Lagrange multipliers.
From $\partial\mathcal{L}/\partial p_{i}(\mathbf{z})=0$ for $i=1,\ldots,n$, we
get

$\displaystyle
p_{i}^{-1}(\mathbf{z})-\lambda_{1}+\lambda_{2}(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z})=0,$

Since $p_{i}(\mathbf{z})$’s sum up to 1 as in (4), letting
$\lambda=-\lambda_{2}/\lambda_{1}$, we have

$p_{i}(\mathbf{z})=\frac{[1+\lambda(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z})]^{-1}}{\sum_{j=1}^{n}[1+\lambda(\tilde{X}_{j1}-\tilde{x}_{1})K_{h}(\tilde{X}_{j}-\mathbf{z})]^{-1}}.$

Using (4) again with (A.4), this gives

$\displaystyle\sum_{j=1}^{n}\left[1+\lambda(\tilde{X}_{j1}-\tilde{x}_{1})K_{h}(\tilde{\mathbf{X}}_{j}-\mathbf{z})\right]^{-1}=n\Bigg{(}\sum_{i=1}^{n}p_{i}(\mathbf{z})[1+\lambda(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z})]\Bigg{)}^{-1}=n,$

and therefore (8) holds. With (5), this gives

$0=\sum_{i=1}^{n}\frac{(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z})}{1+\lambda(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z})}=\left.-\frac{\partial
L(\gamma;\mathbf{z})}{\partial\gamma}\right|_{\lambda}.$

Note that $\frac{\partial^{2}L(\gamma;\mathbf{z})}{\partial{\gamma}^{2}}\geq
0$, implying that $L(\cdot;\mathbf{z})$ is indeed a convex function. ∎

######  Lemma B.2.

Under the assumptions of Proposition 4.6, define

$c(\mathbf{z})=\frac{\frac{\partial
g(\mathbf{z})}{\partial\tilde{x}_{1}}\mu_{2}}{g(\mathbf{z})\nu_{2}}.$

Then,

$\lambda=h^{w}\cdot c(\mathbf{z})(1+o_{p}(1))=\mathcal{O}_{p}(h^{w}).$ (A.5)

###### Proof.

Let

$S_{i}=(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z}).$

Then, by Assumption 4.4, $S_{i}$ is bounded above by some constant $C_{1}$.
Let

$\overline{S^{k}}=\frac{1}{n}\sum_{i=1}^{n}(S_{i})^{k},$

for $k=1,2$. Then, from (A.4), we have

$\displaystyle 0$ $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}(1+\lambda
S_{i})^{-1}S_{i}\geq\left|\lambda\right|\left|\frac{1}{n}\sum_{i=1}^{n}S_{i}^{2}(1+\lambda
S_{i})^{-1}\right|-\left|\overline{S^{1}}\right|\geq\left|\lambda\right|(1+C_{1}\left|\lambda\right|)^{-1}\overline{S^{2}}-\left|\overline{S^{1}}\right|,$

which gives

$\left|\lambda\right|\leq\frac{\left|\overline{S^{1}}\right|}{\overline{S^{2}}-C_{1}\left|\overline{S^{1}}\right|}.$

Using the Taylor expansion [Wand and Jones, 1994], we obtain that

$\displaystyle\mathbb{E}\overline{S^{1}}$
$\displaystyle=\int_{\mathbb{R}^{w}}(u_{1}-\tilde{x}_{1})K_{h}(\mathbf{u}-\mathbf{z})g(\mathbf{u})d\mathbf{u}$
$\displaystyle=h\int_{\mathbb{R}^{w}}u_{1}K(\mathbf{u})g(\mathbf{z}+h\mathbf{u})d\mathbf{u}$
$\displaystyle=h\int_{\mathbb{R}^{w}}u_{1}K(\mathbf{u})\Bigg{(}g(\mathbf{z})+h\sum_{j=1}^{w}u_{j}\frac{\partial
g(\mathbf{z})}{\partial\tilde{x}_{j}}\Bigg{)}d\mathbf{u}+o(h^{2})$
$\displaystyle=h^{2}\left\\{\frac{\partial
g(\mathbf{z})}{\partial\tilde{x}_{1}}\mu_{2}+o_{p}(1)\right\\},$

where the last equation comes from Assumptions 4.4-(i) and (ii). With a
similar argument, we can derive that

$\displaystyle\mathbb{E}\overline{S^{2}}$
$\displaystyle=\int_{\mathbb{R}^{w}}(u_{1}-\tilde{x}_{1})^{2}K_{h}^{2}(\mathbf{u}-\mathbf{z})g(\mathbf{u})d\mathbf{u}=h^{-w+2}\left\\{g(\mathbf{z})\nu_{2}+o_{p}(h)\right\\},$

using Assumption 4.4-(iii). Therefore, we obtain (A.5). ∎

Decomposing $\widehat{F}(\tilde{y}|\mathbf{z})-F(\tilde{y}|\mathbf{z})$ in
bias and variance terms, we get

$\displaystyle\widehat{F}(\tilde{y}|\mathbf{z})-F(\tilde{y}|\mathbf{z})$
$\displaystyle=\frac{\sum_{i=1}^{n}p_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})\\{\mathds{1}(\tilde{Y}_{i}<\tilde{y})-F(\tilde{y}|\mathbf{z})\\}}{\sum_{i=1}^{n}p_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})}$
$\displaystyle=\frac{\sum_{i=1}^{n}p_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})\delta_{i}+\sum_{i=1}^{n}p_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})\\{F(\tilde{y}|\tilde{X}_{i})-F(\tilde{y}|\mathbf{z})\\}}{\sum_{i=1}^{n}p_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})},$

where
$\delta_{i}=\mathds{1}(\tilde{Y}_{t}\leq\tilde{y})-F(\tilde{y}|\tilde{X}_{i})$.
Note that $\mathbb{E}[\delta_{i}]=0$ due to the tower property. Now, let

$b_{i}(\mathbf{z})=b_{i}(\tilde{X}_{i},\mathbf{z})\coloneqq\left[1+h^{w}\cdot
c(\mathbf{z})(\tilde{X}_{i1}-\tilde{x}_{1})K_{h}(\tilde{X}_{i}-\mathbf{z})\right]^{-1}.$

Then, by Lemma B.2, we have that

$p_{i}(\mathbf{z})=n^{-1}b_{i}(\mathbf{z})(1+o_{p}(1)).$ (A.6)

Define

$\displaystyle
J_{1}=n^{-1/2}h^{w/2}\sum_{i=1}^{n}b_{i}(\mathbf{z})\delta_{i}K_{h}(\tilde{X}_{i}-\mathbf{z}),$
$\displaystyle
J_{2}=n^{-1}\sum_{i=1}^{n}\left\\{F(\tilde{y}|\tilde{X}_{i})-F(\tilde{y}|\mathbf{z})\right\\}b_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z}),$
$\displaystyle
J_{3}=n^{-1}\sum_{i=1}^{n}b_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z}),$

so that

$\widehat{F}(\tilde{y}|\mathbf{z})-F(\tilde{y}|\mathbf{z})=\\{(nh^{w})^{-1/2}J_{1}+J_{2}\\}J_{3}^{-1}\\{1+o_{p}(1)\\}.$
(A.7)

Therefore, we will derive Proposition 4.6 by controlling the terms $J_{1}$,
$J_{2}$ and $J_{3}$.

######  Lemma B.3.

Under the assumptions of Proposition 4.6,

$J_{1}=\mathcal{O}_{p}(1).$ (A.8)

###### Proof.

Let

$\xi_{i}=h^{w/2}b_{i}(\mathbf{z})\delta_{i}K_{h}(\tilde{X}_{i}-\mathbf{z}),$

so that $J_{1}=n^{-1/2}\sum_{i=1}^{n}\xi_{i}$. Since
$\mathbb{E}(\delta_{i}|\tilde{X}_{i})=0$, we have that
$\mathbb{E}(\xi_{i})=\mathbb{E}(\mathbb{E}(\xi_{i}|\tilde{X}_{i}))=0$, and
thus

$\displaystyle\mathbb{E}J_{1}=0.$ (A.9)

Also, due to the stationarity of $\tilde{X}_{i}$, we have that

$\mbox{Var}(J_{1})=\mathbb{E}\xi_{i}^{2}+\sum_{i=2}^{n}\left(1-\frac{i-1}{n}\right)\mbox{Cov}(\xi_{1},\xi_{i}).$
(A.10)

By Assumption 4.4, we have that $\lim_{n\to\infty}b_{i}(\mathbf{z})=1$, which
gives $\mathbb{E}(b_{i})=1+o_{p}(1)$. Therefore, through expansion, we have

$\displaystyle\mathbb{E}\xi_{i}^{2}$
$\displaystyle=h^{w}\mathbb{E}\left[\mathbb{E}\left[b_{i}^{2}(\mathbf{z})K_{h}^{2}(\tilde{X}_{i}-\mathbf{z})\delta_{i}^{2}\mid\tilde{X}_{i}\right]\right]$
$\displaystyle=h^{w}\mathbb{E}\left[b_{i}^{2}(\mathbf{z})K_{h}^{2}(\tilde{X}_{i}-\mathbf{z})F(\tilde{y}|\tilde{X}_{i})(1-F(\tilde{y}|\tilde{X}_{i}))\right]$
$\displaystyle=\left[(K^{2})_{h}\ast\\{b_{i}^{2}(\cdot,\mathbf{z})F(\tilde{y}|\cdot)(1-F(\tilde{y}|\cdot))g(\cdot)\\}\right](\mathbf{z})$
$\displaystyle=\nu_{0}F(\tilde{y}|\mathbf{z})(1-F(\tilde{y}|\mathbf{z}))g(\mathbf{z})+o_{p}(1),$

where $\ast$ in the third line is the convolution operator. To control the
second term in the right-hand side of (A.10), we borrow the idea of Masry
[1986]. Choose $d_{n}=\mathcal{O}(h^{-\frac{w}{1+\delta/2}})$ and decompose

$\displaystyle\sum_{i=2}^{n}\left(1-\frac{i-1}{n}\right)\mbox{Cov}(\xi_{1},\xi_{i})=\sum_{i=2}^{d_{n}}\left(1-\frac{i-1}{n}\right)\mbox{Cov}(\xi_{1},\xi_{i})+\sum_{i=d_{n}+1}^{n}\left(1-\frac{i-1}{n}\right)\mbox{Cov}(\xi_{1},\xi_{i}).$

We have that $\left|b_{i}(\mathbf{z})\delta_{i}\right|\leq C_{2}$ for some
constant $C_{2}$. By Assumption 4.3-(iv), we obtain

$\displaystyle\left|\mbox{Cov}(\xi_{1},\xi_{i})\right|$
$\displaystyle=\left|\int_{\mathbb{R}^{w}}\int_{\mathbb{R}^{w}}\xi_{1}\xi_{i}g_{1,i}(\mathbf{u},v)d\mathbf{u}dv-\int_{\mathbb{R}^{w}}\xi_{1}g(\mathbf{u})d\mathbf{u}\int_{\mathbb{R}^{w}}\xi_{i}g(v)dv\right|$
$\displaystyle\leq
C_{2}^{2}h^{w}\int_{\mathbb{R}^{w}}\int_{\mathbb{R}^{w}}K(\mathbf{u})K(v)\left|g_{1,i}(\mathbf{z}-h\mathbf{u},\mathbf{z}-hv)-g(\mathbf{z}-h\mathbf{u})g(\mathbf{z}-hv)\right|d\mathbf{u}dv$
$\displaystyle\leq C_{2}^{2}Mh^{w},$

so that

$\sum_{i=2}^{d_{n}}\left(1-\frac{i-1}{n}\right)\mbox{Cov}(\xi_{1},\xi_{i})=\mathcal{O}_{p}(d_{n}h^{w})=o_{p}(1).$

By Assumption 4.4, we have
$\|(\tilde{X}_{i}-\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})\|\leq C_{3}$, so
that $\left|\xi_{i}\right|\leq Ch^{-w/2}$. Then, by Theorem 17.2.1 of
Ibragimov et al. [1971], we have that

$\left|\mbox{Cov}(\xi_{1},\xi_{i})\right|\leq Ch^{-w}\alpha(i-1).$

Thus, we get

$\displaystyle\sum_{i=d_{n}+1}^{n}\left(1-\frac{i-1}{n}\right)\mbox{Cov}(\xi_{1},\xi_{i})\leq
Ch^{-w}\sum_{i\geq d_{n}}\alpha(i)\leq Ch^{-w}d_{n}^{-1-\delta}=o(1).$

Therefore, we obtain

$J_{1}=\mathbb{E}J_{1}+\mathcal{O}_{p}\left(\sqrt{\mbox{Var}(J_{1})}\right)=\mathcal{O}_{p}(1).$
(A.11)

∎

######  Lemma B.4.

Under the assumptions of Proposition 4.6,

$\displaystyle J_{2}$
$\displaystyle=\frac{1}{2}h^{2}\mu_{2}\mbox{tr}(D_{\mathbf{z}}^{2}F(\tilde{y}|\mathbf{z}))g(\mathbf{z})+o_{p}(h^{2}),$
(A.12) $\displaystyle J_{3}$ $\displaystyle=g(\mathbf{z})+o_{p}(1).$ (A.13)

###### Proof.

By Assumption 1, (4) and (5), we have through expansion that

$\displaystyle
J_{2}=(2n)^{-1}\sum_{i=1}^{n}(\tilde{X}_{i}-\mathbf{z})^{\top}(D_{\mathbf{z}}^{2}F(\tilde{y}|\mathbf{z}))(\tilde{X}_{i}-\mathbf{z})b_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})+o_{p}(h^{2}).$

Since

$\displaystyle\mathbb{E}[(\tilde{X}_{i}-\mathbf{z})^{\top}D_{\mathbf{z}}^{2}F(\tilde{y}|\mathbf{z})(\tilde{X}_{i}-\mathbf{z})b_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})]$
$\displaystyle=\mbox{tr}(D_{\mathbf{z}}^{2}F(\tilde{y}|\mathbf{z})\mathbb{E}[(\tilde{X}_{i}-\mathbf{z})(\tilde{X}_{i}-\mathbf{z})^{\top}b_{i}(\mathbf{z})K_{h}(\tilde{X}_{i}-\mathbf{z})])$
$\displaystyle=h^{2}\mu_{2}g(\mathbf{z})\mbox{tr}(D_{\mathbf{z}}^{2}F(\tilde{y}|\mathbf{z}))+o_{p}(h^{2}),$

we have

$J_{2}=\frac{1}{2}h^{2}\mu_{2}\mbox{tr}(D_{\mathbf{z}}^{2}F(\tilde{y}|\mathbf{z}))g(\mathbf{z})+o_{p}(h^{2}).$
(A.14)

Finally, by applying the expansion argument routinely, we get

$J_{3}=g(\mathbf{z})+o_{p}(1).$ (A.15)

∎

###### Proof of Proposition 4.6 [Cai, 2002, Salha, 2006].

Combining Lemmas B.3 and B.4 with (A.7), Assumption 4.5 gives the result. ∎

###### Proof of Corollary 4.7 [Cai, 2002].

Given $\mathbf{z}$, Proposition 4.6 implies uniform convergence of
$\widehat{F}(\cdot|\mathbf{z})$ to $F(\cdot|\mathbf{z})$ in probability
[Tucker, 1967, p.127-128] since $F(\cdot|\mathbf{z})$ is a CDF. That is,

$\sup_{\tilde{y}\in\mathbb{R}}\left|\widehat{F}(\tilde{y}|\mathbf{z})-F(\tilde{y}|\mathbf{z})\right|\to
0\mbox{ in probability}.$

Given $\varepsilon>0$, let
$\delta=\delta(\varepsilon)\coloneqq\min\\{\beta-F(Q_{\beta}(\mathbf{z})-\varepsilon|\mathbf{z}),F(Q_{\beta}(\mathbf{z})+\varepsilon|\mathbf{z})-\beta\\}$.
Note that $\delta>0$ due to the uniqueness of the quantile. We have

$\displaystyle\mathbb{P}\left(\left|\widehat{Q}_{\beta}(\mathbf{z})-Q_{\beta}(\mathbf{z})\right|>\varepsilon\right)$
$\displaystyle\leq\mathbb{P}\left(\left|F(\widehat{Q}_{\beta}(\mathbf{z})|\mathbf{z})-\beta\right|>\delta\right)$
$\displaystyle\leq\mathbb{P}\left(\sup_{\tilde{y}\in\mathbb{R}}\left|\widehat{F}(\tilde{y}|\mathbf{z})-F(\tilde{y}|\mathbf{z})\right|>\delta\right).$

The uniform convergence of $\widehat{F}(\cdot|\tilde{x})$ in probability gives
the result. ∎

###### Proof of Corollary 4.8.

From the definition of $\hat{C}_{t-1}^{\alpha}$ in (3), we have

$\displaystyle\mathbb{P}\left(Y_{t}\in\hat{C}_{t-1}^{\alpha}(X_{t})\Big{|}X_{t}\right)$
$\displaystyle=\mathbb{P}\left(\tilde{Y}_{n+1}\in\left[\widehat{Q}_{\beta^{*}}(\tilde{X}_{n+1}),\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{X}_{n+1})\right]\Big{|}\,\tilde{X}_{n+1}\right)$
$\displaystyle=F\left(\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right)-F\left(\widehat{Q}_{\beta^{*}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right).$

By Theorem 4.7, we have the consistency of $\widehat{Q}_{\beta}$ for all
$\beta\in(0,1)$. On that, the continuous mapping theorem and Assumption 4.3
gives

$\displaystyle
F\left(\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right)-F\left(\widehat{Q}_{\beta^{*}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right)$
$\displaystyle\to
F\left(Q_{1-\alpha+\beta^{*}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right)-F\left(Q_{\beta^{*}}(\tilde{X}_{n+1})\Big{|}\tilde{X}_{n+1}\right)$
$\displaystyle=1-\alpha,$

where the convergence is in probability. ∎

###### Proof of Theorem 4.9.

From the definition of $\widehat{C}_{t-1}^{\alpha}$ in (3), we have

$\displaystyle\left|\mathbb{P}\left(Y_{t}\in\widehat{C}_{t-1}^{\alpha}(x)\Big{|}X_{t}=x\right)-(1-\alpha)\right|$
$\displaystyle=\left|\mathbb{P}\left(\tilde{Y}_{n+1}\in\left[\widehat{Q}_{\beta^{*}}(\tilde{x}),\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{x})\right]\Big{|}\
\tilde{X}_{n+1}=\tilde{x}\right)-(1-\alpha)\right|$
$\displaystyle=\left|F\left(\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)-F\left(\widehat{Q}_{\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)-(1-\alpha)\right|$
$\displaystyle\leq\left|\widehat{F}\left(\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)-(1-\alpha+\beta^{*})\right|+\left|\widehat{F}\left(\widehat{Q}_{\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)-\beta^{*}\right|$
$\displaystyle\quad+\left|F\left(\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)-\widehat{F}\left(\widehat{Q}_{1-\alpha+\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)\right|$
$\displaystyle\quad+\left|F\left(\widehat{Q}_{\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)-\widehat{F}\left(\widehat{Q}_{\beta^{*}}(\tilde{x})\Big{|}\tilde{X}_{n+1}=\tilde{x}\right)\right|$
$\displaystyle\leq 2\Delta(\tilde{x})+\mathcal{O}_{p}((nh^{w})^{-1/2}+h^{2}),$

where the last inequality comes from (A.7) and the definition of the discrete
gap $\Delta$. ∎

##  Appendix C Additional experiment results

We show additional rolling comparison on the solar (Figure A.1) and wind
dataset (Figure A.2), where the setup and conclusions remain consistent with
those on the electric dataset (Figure 2).

(a) Rolling coverage (boxplot)

(b) Width of prediction intervals (boxplot)

(c) Rolling coverage over prediction time indices

Figure A.1: Rolling coverage and width comparison on the solar dataset by
different methods.

(a) Rolling coverage (boxplot)

(b) Width of prediction intervals (boxplot)

(c) Rolling coverage over prediction time indices

Figure A.2: Rolling coverage and width comparison on the wind dataset by
different methods.